J. Fluid Mech. (2023), vol. 973, A4, doi:10.1017/jfm.2023.742

Pilot-wave dynamics in a rotating frame: theonset of orbital instability
Nicholas Liu1 , Matthew Durey2 and John W.M. Bush1, †1 Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
2 School of Mathematics and Statistics, University of Glasgow, University Place, Glasgow G12 8QQ, UK

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press(Received 1 June 2023; revised 31 July 2023; accepted 1 September 2023)

We report the results of a theoretical investigation of the stability of a hydrodynamic
analogue of Landau levels, specifically circular orbits arising when a millimetric droplet
self-propels along the surface of a vibrating, rotating liquid bath. Our study elucidates
the form of the stability diagram characterising the critical memory at which circularorbits destabilise, and the form of instability. Particular attention is given to rationalisingobservations reported in prior experimental works, including the prevalence of resonant
wobbling instabilities, in which the instability frequency is approximately twice theorbital frequency. We also explore the physical mechanism responsible for the onset of
instability. Specifically, we compare the efficacy of different heuristic arguments proposed
in prior studies, including propositions that the most unstable orbits arise when their radii
correspond to the zeros of Bessel functions or when their associated wave intensity isextremised. We establish a new relation between orbital stability and the mean wave field,
which supersedes existing heuristic arguments and suggests a rationale for the alternate
wobbling and monotonic instabilities arising at onset as the orbital radius is increased
progressively.
Key words: Faraday waves, drops

1. Introduction

Millimetric droplets may bounce and self-propel along the surface of a vertically vibrating
liquid bath (Couder et al. 2005a,b). These walking droplets, or ‘walkers’, provide atangible macroscopic example of wave–particle duality and represent a classical realisationof a pilot-wave system of the form envisaged by de Broglie (1926, 1930) in the 1920s.
Since its discovery in 2005, the hydrodynamic pilot-wave system has provided the basis
for an unexpectedly long list of hydrodynamic quantum analogues (Bush 2015; Bush &Oza 2020). Notable examples include single-particle diffraction and interference (Couder† Email address for correspondence: bush@math.mit.edu© The Author(s), 2023. Published by Cambridge University Press

973 A4-1

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

N. Liu, M. Durey and J.W.M. Bush
& Fort 2006; Pucci et al. 2018; Ellegaard & Levinsen 2020), unpredictable tunnelling(Eddi et al. 2009; Tadrist et al. 2020), Friedel oscillations (Sáenz, Cristea-Platon & Bush2020), spin lattices (Sáenz et al. 2021), and quantum-like statistics in corrals (Harris et al.
2013; Sáenz, Cristea-Platon & Bush 2018). Of particular interest here are the quantisedorbits emerging when a droplet walks in the presence of an imposed force (Fort et al.
2010; Perrard et al. 2014b).
Faraday waves are generated at the free surface of a vibrating liquid bath when thebath’s vibrational acceleration exceeds a critical value known as the Faraday threshold(Benjamin & Ursell 1954; Miles & Henderson 1990). Walking droplets arise just below
this threshold, and so respond only to the waves generated by their own bouncing. A key
feature of the walker system is that the droplets bounce at the frequency of the bath’s most
unstable Faraday waves, specifically half the frequency of the bath’s vibrational forcing.
The resulting resonance between the bouncing droplet and the bath ensures that thepilot-wave field is quasi-monochromatic, with a dominant wavelength prescribed by the
Faraday wavelength. The dynamics is local in that the droplet responds to the local slopeof its guiding wave; however, as the slope is determined by the droplet’s past trajectory,
the droplet dynamics are non-Markovian. The path memory of the system determines themean number of prior impacts that contribute to the wave force acting on the droplet,and is prescribed by the proximity of the vibrational acceleration to the Faraday threshold(Eddi et al. 2011). The quantum features of the pilot-wave hydrodynamic system emerge inthe high-memory limit arising when the vibrational acceleration approaches the Faraday
threshold, and the pilot wave is most persistent.
Orbital quantisation is a canonical feature of the hydrodynamic pilot-wave system,and one of its most compelling as a quantum analogue. At sufficiently high memory,the quasi-monochromatic form of the guiding or ‘pilot’ wave constrains the droplet toquantised dynamical states (Fort et al. 2010; Oza et al. 2014a; Perrard et al. 2014b;
Labousse et al. 2016a). Orbital quantisation has been shown to arise for walkers in arotating frame (Fort et al. 2010; Harris & Bush 2014; Oza et al. 2014a), a simple harmonic
potential (Perrard et al. 2014a,b; Labousse et al. 2016a; Durey & Milewski 2017), andconfinement to a submerged well (Harris et al. 2013; Cristea-Platon, Sáenz & Bush 2018;
Durey, Milewski & Wang 2020a). In the first of these systems, quantised states consistof circular orbits (Fort et al. 2010). In the latter two systems, more complex orbits may
arise, including lemniscates and trefoils (Cristea-Platon et al. 2018; Durey et al. 2020a),and orbits are quantised in both energy and angular momentum (Perrard et al. 2014b). In
all three systems, in the high-memory limit the quantised orbits destabilise, and the droplet
switches intermittently between the accessible orbits, giving rise to statistics reminiscentof their quantum counterpart (Harris et al. 2013; Oza et al. 2014a,b; Labousse et al. 2016a;
Cristea-Platon et al. 2018).
The first investigation of droplets walking in a rotating frame was that of Fort et al.
(2010), who observed that the droplets move in circular orbits (see figure 1). At lower path
memory, the orbital radius depends continuously on the bath rotation rate in a manner
expected to arise for inertial orbits, from a balance between centripetal and Coriolis forces.
At higher memory, the geometric constraint imposed by the monochromatic Faraday wave
field restricts the permissible stable orbital radii, giving rise to an effective quantisation inorbital radii. The authors noted that the Coriolis force takes the same form as the Lorentz
force on a moving charged particle, and so drew the analogy between the quantised inertial
orbits arising in their system and the Landau levels arising when a charged quantum
particle moves in a uniform magnetic field. The observed orbital quantisation was captured
in their accompanying simulations (Fort et al. 2010).
973 A4-2

The onset of orbital instability
(a)

(c)

Ω

0.50
0.25
γ cos(2πft)

0

y/λF

–0.25

(b)

–0.50
–0.50 –0.25 00.25 0.50

x/λF
(d )
0.75
0.50
0.25

y/λF

0

–0.25
–0.50
–0.75
–0.6 –0.3 0 0.3 0.6

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Pressx/λF
Figure 1. (a) A schematic illustration of the physical system of interest, millimetric droplets walking on thesurface of a liquid bath vibrating vertically with acceleration γ cos(2πft) and rotating about the vertical axis
with angular frequency Ω. (b) In the rotating frame, the droplets may execute anticyclonic inertial circularorbits, in which the dominant force balance is between the outward inertial force and the inward Coriolis
force. This force balance is augmented by the wave force, which at high memory may give rise to orbits thatare quantised in radius (Fort et al. 2010; Harris & Bush 2014; Oza et al. 2014a). The wave field generated
by the walking droplet has a characteristic wavelength corresponding to the Faraday wavelength λF . As thememory is increased progressively, these orbits may go unstable via one of two mechanisms. (c) A circularorbit (black dashed circle) of radius r0 /λF = 0.60 destabilises by a monotonic instability (red) characterised
by initially exponential monotonic divergence from the circular path, the result being an orbit of smaller radius.
(d) A circular orbit (black dashed circle) of radius r0 /λF = 0.83 destabilises into a wobbling orbit (green),characterised by a wobbling frequency of approximately twice the orbital frequency, known as a 2-wobble.
In (c,d), the dimensionless vibration parameter Γ = (γ − γW )/(γF − γW ) takes values (c) Γ = 0.7 and(d) Γ = 0.8, where γW and γF are the walking and Faraday thresholds, respectively. Images (a,b) are adaptedfrom Harris & Bush (2014).

Harris & Bush (2014) revisited droplets walking in a rotating frame experimentally,and focused on the destabilisation of the quantised orbits and the onset of chaos as
the memory is increased progressively. They demonstrated the emergence of wobbling(see figure 1d), drifting, wobble-and-leap motions, and erratic trajectories at high path
memory. In the long-path-memory limit, all circular orbits are unstable, and the droplet
transitions intermittently between different quantised circular orbits; thus the histogramof the droplet’s radius of curvature is multimodal, with peaks at the orbital radii ofthe quantised circular orbits. The emerging physical picture in this and other pilot-wavehydrodynamic systems is one of a droplet exciting then navigating its own potential (Bush& Oza 2020).
The theoretical modelling of the walking droplet system developed by Moláček &Bush (2013a,b) formed the basis of the stroboscopic model (Oza, Rosales & Bush 2013;
973 A4-3

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

N. Liu, M. Durey and J.W.M. Bush
Oza et al. 2014a), in which the walking droplet is treated as a continuous source ofmonochromatic waves. The stroboscopic model has been successful in rationalising a
number of features of the hydrodynamic pilot-wave system, including the destabilisationof a bouncing droplet into a walker, and the stability of various static and dynamic boundstates (Bush & Oza 2020). Oza et al. (2014a) used the stroboscopic model to analyse
droplets walking in a rotating frame, and rationalised the onset of orbital quantisationin terms of orbital instability at certain radii. Moreover, simulation of the stroboscopic
model in a rotating frame revealed a variety of wobbling, drifting and quasi-periodic
trajectories (Oza et al. 2014b) consistent with experimental observations reported byHarris & Bush (2014). Both experimental and theoretical studies report that at theonset of wobbling, the wobbling frequency is approximately twice the orbital frequency(Harris & Bush 2014; Oza et al. 2014b). We rationalise such resonant instabilities
herein.
Several connections between stability and system energy have been proposed in thecontext of pilot-wave hydrodynamics. In particular, several investigations have suggested
that the wave field energy is decreased at the onset of several instabilities, including
the transition from bouncing to walking states (Durey & Milewski 2017), the transition
from parallel walkers to promenading pairs (Borghesi et al. 2014; Durey & Milewski2017), and the destabilisation of periodic orbits in a harmonic potential (Durey 2018;
Durey, Milewski & Bush 2018). Moreover, the stability of droplet lattices appears to berelated to the height of the local wave field averaged across all droplets in the lattice,a proxy for the mean droplet gravitational potential energy (Couchman & Bush 2020;
Thomson, Durey & Rosales 2020). Finally, it has been suggested that the onset of orbital
instability in a harmonic potential is controlled by the relative energy contribution of asmall number of wave modes (Labousse et al. 2016a). Here, we assess the value of theseenergy-based arguments for predicting the instability of orbital walker motion on a rotatingbath.
The mean wave field in pilot-wave hydrodynamics has proven to be a useful diagnostic
for interpreting both the dynamical and statistical behaviour of walking droplets; moreover,
it is of particular interest given its proposed relation to the quantum potential in Bohmianmechanics (Bush & Oza 2020). The influence of the mean pilot wave on the droplet
dynamics has been explored in a number of settings, including the quantisation of circularand exotic orbits arising in a harmonic potential (Labousse 2014; Perrard et al. 2014b).
The connection between the mean pilot wave and the emergent droplet statistics was firstreported in an experimental study of walker motion in an elliptical corral (Sáenz et al.
2018), and a mathematical link between the two developed by Durey et al. (2018, 2020a).
Here, we demonstrate the utility of the mean pilot-wave field as a diagnostic for orbital
stability in a rotating frame.
We present herein a theoretical investigation of the orbital instability of a droplet
walking in a rotating frame. We introduce the stroboscopic pilot-wave model in § 2,
and determine an alternative formulation for the orbital stability problem in terms ofvarious integrals that specify the influence of the system’s path memory. In § 3, wededuce analytical expressions for the critical memory at the onset of orbital instability,
and also rationalise the form of instability. For the case of wobbling instabilities,
the associated wobbling frequency is deduced. In § 4, we compare the efficacy ofvarious heuristic arguments for the onset of orbital instability and introduce a moreinsightful heuristic based on consideration of the mean wave field. Finally, in § 5,
we discuss the implications of our findings for a broader class of orbital pilot-wave
systems.
973 A4-4

The onset of orbital instability

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

2. Pilot-wave hydrodynamics

We consider the motion of a millimetric drop of mass m self-propelling along thesurface of a liquid bath vibrating vertically with frequency f and acceleration γ cos(2πft)(see figure 1a). The system is subject to a vertical gravitational acceleration −gez androtates about a vertical axis with angular velocity Ω = Ωez , where ez is the verticalunit vector. When the vibrational acceleration exceeds the Faraday threshold, γ > γF , thebath spontaneously destabilises to subharmonic Faraday waves with period TF = 2/f andwavelength λF = 2π/kF prescribed by the water-wave dispersion relation (Benjamin &Ursell 1954). The parameter regime of interest is γ < γF , where the fluid bath would
remain quiescent in the absence of the droplet. For vibrational acceleration beneath thebouncing threshold, γ < γB , the droplet coalesces into the bath; however, for γ > γB ,it bounces in place on the free surface. As γ is increased further, the drop achieves
resonance with the bath, bouncing at the Faraday frequency. Beyond the walking threshold,γ > γW , this resonant bouncing state destabilises into a dynamic walking state, and thewalking droplet generates a quasi-monochromatic wave field with the Faraday wavelength.
The wave decay time TM depends on the proximity to the Faraday threshold, specifically
TM = Td /(1 − γ /γF ), where Td is the viscous decay in the absence of vibrational forcing(Moláček & Bush 2013b). The memory parameter Me = TM /TF prescribes the numberof prior impacts relevant to the dynamics. In the absence of system rotation, the droplet
follows a rectilinear path at constant speed, u0 . In the rotating frame, the droplet instead
executes anticyclonic inertial orbits with radius r0 and angular frequency ω.
A key notion in our study is that of orbital memory, MeO , the ratio of the memory time
TM to the orbital period 2π/ω, which prescribes the extent to which the orbiting dropinteracts with its own wake (Oza et al. 2014a). When MeO  1, the pilot wave decays
over several orbital periods, the droplet interacts significantly with its own wake, and its
trajectory is strongly influenced by its history. It is in the regime of intermediate orbital
memory, MeO ∼ O(1), that orbital instability first arises, and so will be a focus of our study.
When MeO  1, the droplet is unperturbed by its wake, and instability emerges in the formof in-line speed oscillations that may arise for the rectilinear trajectory of a free walker(Bacot et al. 2019; Hubert et al. 2019; Durey, Turton & Bush 2020b).
We root our discussion in the parameter regime explored by Harris & Bush (2014),
who used a fluid of density 949 kg m−3 , kinematic viscosity 20 cSt, surface tension
0.0206 N m−1 and depth 4 mm, with vibrational frequency f = 80 Hz, and a droplet ofradius 0.4 mm whose free walking speed was approximately u0 = 11 mm s−1 .
2.1. Governing equations
To model the horizontal motion of the walker, we utilise the stroboscopic trajectory
equation developed by Oza et al. (2013, 2014a), whereby the pilot-wave system istime-averaged over one bouncing period, TF = 2/f (Moláček & Bush 2013b). The
droplet’s horizontal position xp (t) thus evolves over time t according to (Oza et al.
2013, 2014a,b):mẍp + Dẋp = −mg ∇h(xp (t), t) − 2mΩ × ẋp .

(2.1a)The drop is propelled by the wave force −mg ∇h(xp (t), t), and also responds to the linear
drag force −Dẋp and the Coriolis force −2mΩ × ẋp . The accompanying pilot wave,
 tAh(x, t) =J0 (kF |x − xp (s)|) e−(t−s)/TM ds,(2.1b)TF −∞973 A4-5

N. Liu, M. Durey and J.W.M. Bush
Dimensionless parameters

Definition

κ0 = m/(DTW )Ω̂ = 2mΩ/D
μ = TW /TMΓ = (γ − γW )/(γF − γW )r0ωωorb = ω/μ−1
β = μ/ω = ωorb
U = r0 ωs
S = Im(s)ξ = S/ω

Inertia-to-drag ratio
Dimensionless rotation vector
Wave decay rate
Vibration parameter
Orbital radius, normalised by kF−1
Orbital angular frequency
Orbital memory parameter
Inverse orbital memory parameter
Orbital speed
Asymptotic complex growth rate of perturbations
Destabilisation frequency
Destabilisation frequency relative to orbital frequency

Table 1. The dimensionless parameters appearing in the pilot-wave system (2.2) and subsequent analysis.

is modelled as a continuous superposition of axisymmetric waves of amplitude A centred
along the droplet’s path, decaying exponentially in time over the memory time scale TM .
The quasi-monochromatic form of the pilot-wave field imposes a geometric constraint
on the droplet’s motion whose effects are most pronounced at high memory, where the
Faraday waves are most persistent. The Faraday wavenumber, drag and wave amplitude
parameters, respectively kF , D and A, are defined in terms of physical quantities inAppendix A.
We project the pilot wave onto the droplet’s path to yield an integro-differential
trajectory equation for the droplet (Oza et al. 2013) that may be expressed in dimensionlessvariables as (Oza et al. 2014a; Oza, Rosales & Bush 2018)https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Pressκ0 x̂¨ p + x̂˙ p = 2 t̂

J1 (|x̂p (t̂) − x̂p (s)|)(x̂p (t̂) − x̂p (s)) e−μ(t̂−s) ds − Ω̂ × x̂˙ p ,−∞ |x̂p (t̂) − x̂p (s)|

(2.2)


where x̂p = kF xp , t̂ = t/TW , and TW = 2DTF /mgAkF2 is the memory time at the onset ofwalking, γ = γW (Oza et al. 2013; Durey et al. 2020b). The dimensionless parameters μ =TW /TM > 0 and κ0 = m/DTW describe the wave decay rate and the relative importanceof inertial and drag forces, respectively, and Ω̂ = 2mΩ/D = Ω̂ez is the dimensionlessrotation vector (see table 1).
We characterise the pilot-wave dynamics in terms of the dimensionless vibration
parameter Γ = (γ − γW )/(γF − γW ) = 1 − μ (Bush 2015; Oza et al. 2018; Durey et al.
2020b), which increases with increasing path memory. We note that Γ = 0 corresponds
to the walking threshold in the absence of a Coriolis force (γ = γW ), while Γ = 1
corresponds to the Faraday threshold (γ = γF ), and thus infinite path memory (Bush2015). The experimental parameter regime of Harris & Bush (2014) corresponds to κ0 ≈1.6. We note that typically, κ0 takes values in the range 0.8  κ0  1.6 in the laboratory;
likewise, the dimensionless rotation rate Ω̂ is restricted to the interval 0 ≤ |Ω̂|  1.3(Harris & Bush 2014; Oza et al. 2014a, 2018). Henceforth, we thus treat κ0 and Ω̂ as O(1)quantities, whose influence on the pilot-wave dynamics we characterise through systematic
asymptotic analysis.
973 A4-6

The onset of orbital instability
2.2. Orbital dynamics
We characterise orbits in terms of their radius r0 and angular frequency ω > 0. By omitting
hats and substituting xp (t) = r0 (cos(ωt), sin(ωt)) into the trajectory equation (2.2), we
express the radial and tangential force balances as (Oza et al. 2014a)
 ∞ 
 ωt 
 ωt 
−κ0 r0 ω2 = 2
J1 2r0 sin(2.3a)
sin
e−μt dt + Ωr0 ω,
2
2
0
 ∞ 
 ωt 
 ωt 
J1 2r0 sin(2.3b)
cos
e−μt dt,r0 ω = 2220which may be solved for r0 and ω given κ0 , μ and Ω.
In our dimensionless notation, a suitable proxy for the orbital memory MeO is ωorb =ω/μ, which is the ratio of the wave decay time scale μ−1 to the orbital time scale ω−1 ∼r0 /u0 , where u0 is the steady walking speed in the absence of bath rotation (Oza et al.
2013; Durey et al. 2020b). The orbital speed U = r0 ω typically remains close to the free√walking speed u0 and satisfies U < 2 for all parameter values (see § 4.1). As U depends
only weakly on the orbital radius at fixed memory, we note that ωorb = U/(r0 μ) decreases
with increasing r0 .
2.3. Orbital stability
In order to characterise the droplet’s response to perturbations from a circular orbit, weapply linear stability analysis. Following the framework developed by Oza et al. (2014a),
we linearise the trajectory equation (2.2) about the orbital solution expressed by (2.3).
Specifically, we write
xp (t) = rp (t) cos θp (t), sin θp (t) ,(2.4a)

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

where rp (t) and θp (t) are the time-varying radial and angular polar coordinates of thedroplet’s position, respectively. For a small perturbation from an orbital trajectory, we
consider solutions of the formrp (t) = r0 + r1 (t)andθp (t) = ωt + θ1 (t),(2.4b)

where r0 and ω satisfy the orbital equations (2.3), and  1 is a small parameter. We
substitute (2.4) into (2.2), retain terms to O( ), and then take the Laplace transform of theresultant linear equations. It follows that the perturbed trajectory’s asymptotic complex
growth rates s satisfy F(s) = 0, whereF(s) = A (s) D (s) + B (s) C (s),(2.5)

and the stability coefficients are defined (in a form equivalent to Oza et al. 2014a, 2018) as

A (s) = κ0 (s2 − 2ω2 ) + μ + s − 2Ωω + C0 (s) + I1 (s) − 2 I0 (0),(2.6a)

B (s) = 2κ0 ωs + Ωs − μ(κ0 ω + Ω) − S0 (s),(2.6b)

C (s) = 2κ0 ωs + 2ω + Ωs + μ(κ0 ω + Ω) − S0 (s),(2.6c)

D (s) = κ0 s2 + s − μ + C0 (s) − I1 (s).

(2.6d)
973 A4-7

N. Liu, M. Durey and J.W.M. Bush
Of particular interest in our investigation are the integrals (defined for Re(s) > −μ andany integer m ≥ 0) ∞
 ωt Im (s) =J2m 2r0 sin
e−(μ+s)t dt,
20 ∞ ωt 
Cm (s) =J2m 2r0 sincos(ωt) e−(μ+s)t dt,
2
0
 ∞

 ωt 
sin(ωt) e−(μ+s)t dt,
Sm (s) =J2m 2r0 sin
20(2.7a)(2.7b)(2.7c)

which encode the effects of memory on the stability problem, and present most of the
difficulty in solving the stability problem analytically. One important contribution of
our study is the exact analytical evaluation of the stability integrals (2.7) in terms ofBessel functions of complex order. Specifically, we derive in Appendix B the closed-form
expression

Im (s) =π
Jm+iη (r0 ) Jm−iη (r0 ) csch(πη),
ω

(2.8)

where η = (μ + s)/ω. Moreover, by representing cos(ωt) and sin(ωt) in terms of complex
exponential functions, we deduce that1
(Im (s − iω) − Im (s + iω)) ,2i(2.9a,b)which we use to derive similar closed form formulae for Cm and Sm . Using (2.8)–(2.9a,b),
we derive in Appendix B simplified expressions for each of the stability integrals appearingin (2.6) in terms of products of Bessel functions of the first kind, Jν (r0 ), and theirderivatives Jν (r0 ), where the complex order ν takes values ν ∈ {±i(μ + s)/ω}. We thenutilise asymptotic expansions of each integral evaluation to characterise orbital instability(§ 3).
Motivated by our exact analytical evaluation of the stability integrals, we seek to recast
the force balance equations (2.3) in a similar manner. To simplify our investigation, weparametrise the orbital dynamics entirely in terms of the radius r0 (Oza 2014), thereby
effectively eliminating Ω from the stability problem. This elimination process is achieved
by first recasting the radial force balance (2.3a) as

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

Cm (s) =1(Im (s + iω) + Im (s − iω))22
Ω = −κ0 ω −
r0 ω

 ∞
0

and Sm (s) = ωt 
 ωt 

J1 2r0 sinsin
e−μt dt,22(2.10)

where we observe that the integral in (2.10) may be expressed as ∞0

973 A4-8 ωt 
 ωt 
1 ∂ I0 (0)J1 2r0 sin.
sin
e−μt dt = −
2
2
2 ∂r0

(2.11)

The onset of orbital instability
By combining (2.10) and (2.11), we eliminate Ω in the stability coefficients (2.6), yielding
2 ∂ I0 (0)− 2I0 (0) + μ + s + C0 (s) + I1 (s),r0 ∂r0μ − s ∂ I0 (0)B (s) = κ0 ωs −− S0 (s),r0 ω∂r0
μ + s ∂ I0 (0)C (s) = κ0 ωs + 2ω +− S0 (s),r0 ω∂r0A (s) = κ0 s2 −D (s) = κ0 s2 + s − μ + C0 (s) − I1 (s).

(2.12a)(2.12b)(2.12c)(2.12d)Finally, we reduce the tangential force balance by integrating (2.3b) by parts, from which
it follows that the orbital speed U = r0 ω satisfies (Oza et al. 2014a)

r2 ω2
1
I0 (0) =1− 0.
(2.13)μ2
For any given r0 > 0, the orbital stability problem may be expressed solely in terms of the
reduced tangential force balance (2.13) and the stability condition F(s) = 0, both of which
are defined in terms of the stability integrals (2.7).
The orbital solution is unstable if there are any roots s of F satisfying Re(s) > 0. By
denoting s∗ as the unstable root with largest real part, the instability is monotonic if
Im(s∗ ) = 0, and oscillatory otherwise. The stability function F has trivial eigenvalues at 0
and ±iω, corresponding to rotational and translational invariance of the orbital motion,
respectively (Oza et al. 2014a). It follows, therefore, that the non-trivial roots of thestability problem satisfy G(s) = 0, where

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

G(s) =

A (s) D (s) + B (s) C (s).
s(s2 + ω2 )(2.14)

We apply the method of Delves & Lyness (1967) to find the roots of G in the domain
over which G is analytic, i.e. Re(s) > −μ. To ascertain whether a particular orbital stateis stable or unstable, typically we utilise a rectangular integration contour spanning thedomain Re(s) ∈ [0, 20] and Im(s) ∈ [0, 5], which we find to be sufficient for identifying
all roots with a positive real part across the experimentally based parameter regime
considered in this study (0 ≤ Γ ≤ 0.99). This approach differs from that of Oza et al.
(2018), who instead applied the methodology of Delves & Lyness (1967) to F, integrating
F /F over a deformed contour specifically chosen to avoid the trivial zeros at s = 0 ands = ±iω. The method presented here instead removes the singularities analytically, thereby
avoiding contour deformations near the trivial zeros of F; however, local Taylor series
approximations are necessary to avoid numerical difficulties arising sufficiently close to
the removable singularities of G.
In figure 2(a), we follow Oza et al. (2014a, 2018) in presenting the dependence ofthe orbital radius on the bath rotation rate for Γ = 0.8. In figure 2(b), we summarise
the stability behaviour for all Γ . As path memory is increased progressively, stable
circular orbits (blue) destabilise via either a monotonic (red, see figure 1c) or oscillatory(green, see figure 1d) instability mechanism. Associated monotonic (red) and oscillatory(green) instability ‘tongues’ emerge in the stability diagram, with the tip of each tongue
corresponding to the onset of a new instability. The blue regions between the instability
tongues correspond to regions of orbital quantisation. Notably, as memory is increased
beyond the tip of an oscillatory instability tongue, the orbital instability typically manifests
973 A4-9

N. Liu, M. Durey and J.W.M. Bush
(a)

(b)

6
5

0.9

r0/λF

4

Γ

3
2

0.7
0.6

1
0

0.8

0.5

1.0

1.5

2.0

2|Ω|λF /u0

2.5

3.0

0.5
0

1

2

3

4

5

6

r0/λF

Figure 2. Orbital stability and the onset of quantisation for κ0 = 1.6, where κ0 = m/(DTW ) is the inertial
coefficient appearing in (2.2). (a) The dependence of orbital radius on bath rotation rate for Γ = 0.8, whereΓ = (γ − γW )/(γF − γW ) is the dimensionless vibration parameter. (b) The delineation of orbital stability forany given radius and memory. The stability boundary is highlighted in white, and the yellow line at Γ = 0.8
corresponds to the curve in (a). In both plots, stable orbital states are indicated in blue, while oscillatory andmonotonic instabilities are highlighted in green and red, respectively. We note that in (b), the instability tongues
alternate between monotonic and wobbling as the orbital radius is increased progressively. The cyan dots in
the red and green regions denote the parameters corresponding to the monotonic and wobbling trajectories
presented in figures 1(c,d). The pink curve denotes an instability related to in-line speed oscillations of theorder of the Faraday wavelength (Durey et al. 2020b), which is subdominant for r0 /λF  4, and truncates theinstability tongues for r0 /λF  4.

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

as a wobbling orbit, in which the orbital centre remains approximately constant, but theradius of curvature exhibits small-amplitude oscillations with a frequency approximately
twice that of the orbital frequency, as reported in the experiments of Harris & Bush (2014)
and the numerical simulations of Oza et al. (2014b) (see figure 1d). Furthermore, theinstability tongues appear to have a periodic structure, with the critical memory increasing
with increasing orbital radius. We observe that the monotonic and wobbling instability
tongues are nested, with monotonic instability tongues forming at lower memory than theneighbouring wobbling instability tongues.
2.4. The onset of instability: asymptotic scaling relationships
Although the stability integrals (2.7) may be evaluated analytically (see (2.8)–(2.9a,b)),the purpose of this subsection is to motivate the asymptotic scaling relationships arisingnear the tip of each instability tongue. In particular, we determine the main contributions
to each stability integral arising along a stability boundary (s = iS with S real) for large
orbital radius (r0  1), for which each integrand is highly oscillatory. Using (2.9a,b) toexpress Cm and Sm in terms of Im , we henceforth focus our attention on the study of Im .
Furthermore, by recognising that Im is a Laplace transform of a periodic function, we
reduce the integral (2.7a) to 2π
Lm (ξ )
θ , where Lm (ξ ) =e−(β+iξ )θ dθ.
Im (iS) = 
J2m 2r0 sin−2π(β+iξ
)
2
ω 1−e
0
(2.15)
Here, β = μ/ω is the inverse orbital memory, and ξ = S/ω is the scaled destabilisation
frequency (see table 1). Based on experimental and numerical observations of monotonic
and wobbling instabilities (Harris & Bush 2014; Oza et al. 2014b), we assume henceforththat the destabilisation frequency is comparable to the orbital frequency (i.e. ξ = O(1)).
Notably, our analysis does not account for the instability associated with in-line speed
973 A4-10

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

The onset of orbital instability
oscillations arising at larger orbital radii (denoted by the pink curve in figure 2b), for which
S ∼ U, or ξ ∼ r0 (Durey et al. 2020b). Finally, by writing β = μr0 /U and noting that μ
at the tip of successive wobbling (or successive monotonic) instability tongues decreases
with increasing orbital radius (see figure 2b), we deduce that the magnitude of β is at
most of O(r0 ) when the orbital radius is large (since U = O(1) for all orbital radii). Our
analysis in this subsection determines the precise scaling relationship between β and r0 ,namely β = O(ln r0 ), where ln denotes the natural logarithm.
Before proceeding with the asymptotic expansions, we provide a physical interpretation
for the integral Lm (ξ ). The argument 2r0 sin(θ/2) ≥ 0 is the length of the chord spanningtwo points lying an angle θ apart on a circle of radius r0 . This distance reflects the influenceof the droplet’s path memory on the evolution of the perturbed trajectory, where the extent
of the path memory is controlled by the damping rate β > 0. Notably, e−πβ is the wave
damping factor over half an orbital period, which accounts for the contribution of waves
generated when the droplet was last diametrically opposite its current position; likewise,
e−2πβ determines the wave damping factor over a complete orbital period. Finally, thefactor e−iξ θ accounts for oscillations in the perturbed droplet trajectory. For large r0 , theintegrand of Lm (ξ ) is generally highly oscillatory, with dominant contributions arisingover non-oscillatory intervals centred about critical points; these critical points are either
internal points of stationary phase, or boundary points arising when the argument of J2mvanishes, i.e. at θ = 0 and θ = 2π (Bleistein & Handelsman 1975). The internal pointsof stationary phase arise when the argument of the Bessel function is stationary, i.e. atθ = π. We now proceed to determine the magnitude of the contribution made by eachcritical point.
We first examine the contributions to Lm arising about θ = 0 and θ = 2π, which wedenote by Lm,0 and Lm,2π , respectively. We derive in Appendix C the leading-ordercontribution Lm,0 = O(r0−1 ), which is valid when ξ = O(1) and β is of maximum sizeO(r0 ); both of these conditions are met near the tip of each instability tongue. Using
the structure of the integrand of Lm (ξ ), we determine similarly that the leading-ordercontribution about θ = 2π satisfies Lm,2π (ξ ) = e−2π(β+iξ ) Lm,0 (ξ ); thus the relative sizeof Lm,0 (ξ ) and Lm,2π (ξ ) is controlled by the orbital damping factor e−2πβ . Finally, we
use the method of stationary phase (see Appendix C) to determine that the interior point
contribution to Lm at π has magnitude

Lm,π = O

e−πβ
r0

.

(2.16)The relative weight of the integral contributions about θ = 0, π and 2π seemingly
decreases consecutively by a factor e−πβ . However, the contribution about θ = π becomessignificant when eπβ scales algebraically with r0 , corresponding to a strong influenceof the waves generated diametrically opposite the droplet’s current position, as ischaracteristic of high orbital memory (Fort et al. 2010). Indeed, as is evident in figure 3,
the tips of each instability tongue, both monotonic (red line) and oscillatory (green line),satisfy the asymptotic scaling relationship eπβ = O(r02 ), which motivates the asymptotic
scaling relationships utilised in the forthcoming analysis (§ 3). In fact, we observe thateach instability tongue is bounded above in memory by either the oscillatory instability
threshold satisfying eπβ = O(r0 ) (gold dashed line), or the in-line speed oscillation
instability threshold (purple dashed line) arising for free walkers (Bacot et al. 2019; Hubert
et al. 2019; Durey et al. 2020b). Detailing the latter instability (denoted by the pink curvein figure 2b), which manifests as in-line speed oscillations along the circular orbit with
973 A4-11

N. Liu, M. Durey and J.W.M. Bush
100

e–πβ

10–2
O(r0–1)O(r0–2)O(r0–2)
100.5

101.0

101.5

r0kF
Figure 3. Envelopes of the stability boundary, the scalings of which are deduced from our stationary phaseanalysis (see § 2.4). The blue curve corresponds to the white stability boundary in figure 2(b). The green and redlines correspond to the wobbling and monotonic instability envelopes. The gold and purple dashed lines denotethe upper bounds on the existence of stable circular orbits, corresponding to the transition from wobbling tomonotonic instabilities, and the onset of in-line speed oscillations (pink curve in figure 2b), respectively. The
half-orbit wave damping factor e−πβ scales as (r0 kF )−2 for the green and red lines, and as (r0 kF )−1 for thegold line.

amplitudes of the order of the Faraday wavelength, and arises for orbits so large as to beinaccessible within the laboratory, will be the subject of a future investigation.

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

3. The onset of instability

In figure 4, we compare the scaled destabilisation frequency ξ = S/ω, computed along
the stability boundary for the first wobbling instability tongue (see figure 2b), to that ofnonlinear wobbling states arising just beyond the instability threshold (see figure 1d), as
reported in the experimental study of Harris & Bush (2014) and the numerical simulationsof Oza et al. (2014b). The nonlinear wobbling frequency remains close to that predictedby the linear stability analysis, with the two coinciding in the small-wobbling-amplitude
limit. Furthermore, both frequencies remain close to twice the orbital frequency, which
serves to further motivate our analysis of resonant wobbling instabilities. We proceed
to elucidate these observations by means of systematic asymptotic analysis performed
along the stability boundary. Our analysis will also rationalise the quasi-periodicity andenvelopes of the instability tongues, and the influence of the inertial coefficient κ0 on
orbital stability.
We characterise the onset and form of each instability tongue using asymptotic
analysis valid when r0  1. For wobbling instabilities, the destabilisation frequency
S is comparable to the orbital frequency ω; thus we assume that the dimensionlesswobbling frequency is ξ = S/ω = O(1). Moreover, the dimensionless orbital speed U =r0 ω typically remains close to the free walking speed u0 = O(1) (Durey et al. 2020b).
We thus replace ω with U/r0 in the stability equations (2.12) and tangential force balance(2.13), and henceforth assume U = O(1). The dependence of the orbital memory ωorb
on the orbital radius is more subtle. As motivated in § 2.4 and evidenced in figure 3,
the tip of each instability tongue is characterised by the scaling relationship eπβ =−1), which represents the key dominant balance underpinning ourO(r02 ) (where β = ωorb
asymptotic analysis. Notably, this dominant balance implies that the orbital memory
973 A4-12

The onset of orbital instability
2.4
2.3

ξ

2.2
2.1
2.0
1.9

1.0

1.5

2.0

2.5

2|Ω|λF /u0
Figure 4. The dependence of the ratio of the wobbling frequency to the orbital frequency ξ on the bathrotation rate Ω, just beyond the onset of instability for the first wobbling tongue in figure 2(b). Data indicating
the wobbling frequencies are reported in the experimental study of Harris & Bush (2014, figure 7b) (redsquares, 0.961 ≤ γ /γF ≤ 0.978) and the numerical simulations of Oza et al. (2014b, figures 2c,d) (black
triangles, 0.952 ≤ γ /γF ≤ 0.967), corresponding to κ0 = 1.6, where κ0 = m/DTW . The blue curve indicates
our prediction for the wobbling frequency along the stability boundary (white curve in figure 2b), corresponding
to the small-wobbling-amplitude limit.

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

generally decreases with increasing orbital radius, such that ωorb ∼ (ln(r0 ))−1 whenr0  1.
3.1. Asymptotic expansion
We proceed by using the asymptotic expansions of the integrals in (2.12) to identify theimaginary roots of the stability function G defined in (2.14) when r0  1. Specifically,
we seek the roots that minimise the critical memory of instability, which arise along thestability boundary in figure 2(b). We utilise asymptotic expansions for the closed-form
expressions of the stability integrals (see Appendix B), which involve Bessel functionsJν (r0 ) of the first kind with complex order ν, and their derivatives with respect to argumentJν (r0 ) (see (B7)). As the complex order takes values ν ∈ {±iβ, ±i(β ± iξ )} along thestability boundary, where β = O(ln(r0 )) and ξ = O(1) for r0  1, the argument of each
Bessel function is asymptotically large relative to its order. We may thus expand each of thestability coefficients in (2.12) when r0  1, utilising the dominant balance eπβ = O(r02 ),with details presented in Appendix D. Likewise, we√use the large-argument expansions of
the Bessel functions, valid when β = O(ln(r0 ))  r0 , to deduce thatβ
1U2 = 2 1 −+O 3 ,(3.1)
r0
r0
which will be utilised throughout the following analysis.
As our aim is to identify the imaginary roots of the equation G(s) = 0, we proceedby substituting the asymptotic approximations expressed in (D2) into (2.14). Moreover,
we determine the orbital speed U, scaled destabilisation frequency ξ , and inverse orbital
memory β, by means of an asymptotic expansion in terms of the small parameter r0−1 ,namely



U1
1
ξ1
1
β1
1
+ O 2 , ξ = ξ0 ++ O 2 , β = β0 ++O 2 .
U = U0 +
r0
r0
r0
r0
r0
r0
(3.2a–c)973 A4-13

N. Liu, M. Durey and J.W.M. BushIn § 3.2, we present the solution to the leading-order problem, for which we demonstrate
systematically that ξ = 2 + O(r0−2 ) at the tip of each wobbling instability tongue,
corresponding to a wobbling instability with wobbling angular frequency 2ω, the so-called
2-wobble (Harris & Bush 2014; Oza et al. 2014b; see figure 2) and find the critical memory
for wobbling and monotonic instabilities. To determine how the scaled destabilisation
frequency varies away from the tip of each wobbling instability tongue, we extend theasymptotic procedure to incorporate higher-order corrections in § 3.3.
3.2. Leading-order solution
We proceed to determine the leading-order solution, corresponding to the values of U0 ,
ξ0 and β0 . By using the asymptotic relationship D = iξ C + O(r0−3 ) deduced in (D2), we
find that the stability condition (2.14) satisfies1
C (iξ A + B )=O 4 ,(3.3)
iξ(ξ 2 − 1)r0

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

whereupon substituting the leading-order expressions for A , B , U, ξ and β from (D2)and (3.2a–c) results inκ0 U03 + 1
C 2 sin(2r0 )1= O 3 . (3.4)(csch(π(β0 + iξ0 )) + csch(πβ0 )) +22U0ξ0 − 1
r0
r0
As C = 2U0 /r0 + O(r0−2 ) is non-zero to leading order (see (D2)), the leading-ordersolution to (3.4) may be found by setting the term in square brackets equal to zero.
Furthermore, the imaginary part of (3.4) can be satisfied only when the destabilisation
frequency ξ0 is an integer, whose possible values will be the focus of the remainder ofthis subsection. We proceed to eliminate possible integer values for ξ0 by looking for thesolutions to the stability problem that occur at the highest possible value of β0 , as thesesolutions correspond to the instabilities arising at lowest memory for a given orbital radius.
Our analysis will show that only two solutions are possible: (i) ξ0 = 0, corresponding to a
monotonic instability; and (ii) ξ0 = 2, corresponding to a 2ω instability.
To explore the possibility of ξ0 being odd, we first consider the limit ξ0 → 1 in (3.4).
By applying L’Hôpital’s rule, we find that the leading-order stability condition reduces to
1 + κ0 U03 + iπr02 sin(2r0 ) coth(πβ0 ) csch(πβ0 ) = 0.

(3.5)

As the real parts cannot be balanced (since U0 > 0), there are no solutions to this equation.
Similarly, if ξ0 were odd and not equal to 1, the leading-order stability condition (3.4)
would become(3.6)
1 + κ0 U03 = 0,
which is also impossible to satisfy. We thus conclude that ξ0 cannot be odd, meaning thatthe destabilisation frequency must be an even multiple of the orbital frequency.
To explore the possible even values of ξ0 , we denote ξ0 = 2n (where n is an integer)
and use the approximation sinh(x) ≈ cosh(x) ≈ 12 ex for x  1; as such, the leading-orderstability condition (3.4) reduces toeπβ0 =

8r02 sin(2r0 )(1 − 4n2 )(1 + κ0 U03 )

.

(3.7)

To be consistent with the assumed scaling of eπβ0 = O(r02 ), we require sin(2r0 ) = O(1).
As sin(2r0 ) can be either positive or negative, the lowest memory (or largest β0 ) condition
973 A4-14

The onset of orbital instability
requires maximising the magnitude of the right-hand side of (3.7). In the case n = 0, we√
have a monotonic instability: by noting that U0 = 2 from (3.1) and (3.2a–c), we thus
arrive at the monotonic stability boundary
8r02 sin(2r0 )2μmon =,ln√
πr0
1 + 2 2κ0√

(3.8)

which is valid when sin(2r0 ) > 0 and sin(2r0 ) = O(1). For the case n =/ 0, we observethat 1 − 4n2 < 0; we thus deduce the requirement sin(2r0 ) < 0. The magnitude of theright-hand side of (3.7) is then minimised at n = 1, which corresponds to ξ = 2, or s =2iω. We have thus demonstrated that the destabilisation frequency along wobbling stability
boundaries is approximately twice the orbital angular frequency. By substituting n = 1
into (3.7) and rearranging, we deduce that the corresponding critical wave decay rate for awobbling instability is8r02 sin(2r0 )2μwob =ln −,√πr03(1 + 2 2κ0 )
√

(3.9)

which is valid when sin(2r0 ) < 0 and sin(2r0 ) = O(1).
The asymptotic expressions (3.8)–(3.9) for the instability memory give rise toalternation between wobbling (sin(2r0 ) < 0) and monotonic (sin(2r0 ) > 0) instabilities
with increasing orbital radius. In addition, μmon and μwob are maximised (corresponding
to the tip of each instability tongue) when, to leading order for large orbital radius,sin(2r0 ) = 1 and sin(2r0 ) = −1, respectively. These extrema thus determine the wobbling
and monotonic envelopes

8r2
2
ln
√ 0
πr0
3(2 2κ0 + 1)√

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

μenv
wob =


8r2
2
ln √ 0,πr0
2 2κ0 + 1
√

and

μenv
mon =

(3.10a,b)

which are represented in figure 5(a) by the green and red dashed curves, respectively.
of wobbling and monotonic
Notably, increasing κ0 increases the critical memory√
env = 2 ln(3)/(πr ) > 0, we conclude that−μinstabilities. Furthermore, since μenv
0
mon
wob
the envelope of the monotonic instabilities arises at a lower memory than that of wobbling
instabilities, as evident from the red and green lines in figure 3.

3.3. First-order wobbling and monotonic solutions
The leading-order analysis presented in § 3.2 determined approximations for the wobbling
frequency and the critical memory at onset of instability, from which we deduced themost unstable orbital radii. However, the leading-order analysis did not provide insight
into the behaviour of the instability frequency along each wobbling tongue, specifically theextent to which wobbling instabilities are approximated by a 2-wobble. So as to investigate
this behaviour and so determine U1 , ξ1 and β1 , we proceed to solve the tangential force
balance and the stability problem in (3.1) and (3.3) to next order by using the expansionsin (3.2a–c). The expanded tangential force balance (D2) shows that U1 = −β0 /U0 . By
substituting into the stability condition iξ A + B = 0, using (D2) and (3.2a–c), retaining
next-order terms and solving for the real and imaginary parts, we deduce that wobbling
973 A4-15

N. Liu, M. Durey and J.W.M. Bush(a) 1.0(b) 1.0

0.9

0.9

0.8

0.8

Γ

0.7

(c)

1

2

3

0.7

4

(d )

3

1

2

1

2

3

4

3

4

3

2

2

1

1

ξ

0

1

2

r0 /λF3

40

r0 /λF

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

Figure 5. Comparison of the numerical solution and asymptotic approximations to the wobbling andmonotonic stability boundaries. (a,b) The critical dimensionless vibration parameter at onset, Γ = (γ −γW )/(γF − γW ), and (c,d) the instability frequency relative to the orbital frequency, ξ = S/ω, for κ0 = 1.6.
We compare the numerical solution (grey) with (a,c) the leading-order ((3.8)–(3.9)) and (b,d) the first-order((3.10a,b)–(3.12)) asymptotic results. (a,b) The critical memory compared to the asymptotic wobbling (green)and monotonic (red) instability boundaries. In (a), the leading-order envelopes (3.10a,b) are denoted by dashedcurves. (c,d) The critical wobbling frequency and its asymptotic counterpart (orange). All the asymptotic resultspresented in (a–d) are valid when sin(2kF r0 ) = O(1), kF r0 is large, and ξ = 0 (monotonic instabilities) orξ ≈ 2 (2-wobbles). Discontinuities in the grey curves (c,d) reflect transitions between monotonic (ξ = 0) andwobbling (ξ ≈ 2) instabilities.

instabilities (ξ0 =/ 0) have the first-order correction
4β0
cot(2r0 ) andξ1 = −π√
64β0 + 27π β02
3 2κ0 β0
+
β1 =
−
√
π
12π2
π(2 2κ0 + 1)

cot(2r0 ).

(3.11a,b)
Similarly, for the monotonic instability (ξ0 = 0), one may deduce that ξ1 = 0 and√
2πβ0 + 1 β02
3 2κ0 β0
+
−
β1 =
√
4π
π
π(2 2κ0 + 1)

cot(2r0 ).

(3.12)We recall that the asymptotic solution is valid when sin(2r0 ) = O(1). Notably, ξ1 vanishes
for wobbling instabilities when cos(2r0 ) = 0, corresponding to s = iω(2 + O(r0−2 )) at themost unstable orbital radii (see § 3.2). Wobbling instabilities are thus driven by a resonancebetween the destabilisation frequency and the orbital frequency, with a larger criticalmemory necessary to destabilise orbits whose instability frequency deviates from twice
the orbital frequency.
973 A4-16

The onset of orbital instability
3.4. Comparison to the numerical instability tongues
To buttress our asymptotic developments, we compare the numerically deduced instability
tongues to our asymptotic formulae in figure 5, with leading-order results presented infigures 5(a,c) and first-order corrections in figures 5(b,d). Notably, our asymptotic results
are valid when r0  1, sin(2r0 ) = O(1) and either ξ = 0 (corresponding to monotonic
instabilities) or ξ ≈ 2 (corresponding to wobbling instabilities). Despite these restrictions,the leading-order memory captures the main features of the stability tongues, and theasymptotic instability frequency closely matches the numerical behaviour for ξ = 0 or
ξ ≈ 2. The success of our asymptotic results even for orbits of moderate radius is rooted
in the choice of expansion parameter, namely kF r0 in dimensional variables, which isassumed to be large. We note that all orbits presented in figure 5 satisfy kF r0 ≥ 2π
(or r0 /λF ≥ 1), which is evidently sufficiently large for our asymptotic results to yieldreasonable agreement. Finally, we note that a similarly favourable agreement between our
asymptotic and numerical results was obtained across a wide range of κ0 values, including
for those inaccessible in the laboratory (see § 4.1).
In summary, our asymptotic results explain the preponderance of 2-wobbles, with anexact resonance arising at the most unstable radius of each wobbling instability tongue.
Moreover, we quantify the detuning from an exact resonance and the corresponding
increase in the critical memory for nearby orbital radii. Finally, our asymptotic results
demonstrate that the instability tongues alternate between wobbling and monotonic
instabilities as the orbital radius is increased progressively, with the envelope of themonotonic instabilities arising at a memory lower than that of wobbling instabilities(figure 3).

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

4. Physical interpretations of the wobbling and monotonic instability tongues

Having established concise asymptotic formulae describing the onset of wobbling andmonotonic instabilities, we now deepen our physical understanding by comparing theefficacy of different heuristic arguments for the critical radii at the onset of instability.
Specifically, we first compare the predictions of several energy-based heuristics suggested
in prior investigations to those of our numerical and asymptotic results (§ 4.1). We thenpropose a new heuristic based on the form of the mean pilot wave (§ 4.2), which we show
to be more fruitful than prior heuristics.

4.1. Energy-like heuristics
The hydrodynamic pilot-wave system is a driven-dissipative system: energy is suppliedby the system vibration and ultimately lost through viscous dissipation. Nevertheless,quasi-steady and periodic dynamical states arise in which the energy input preciselybalances that lost through dissipation, and energy is exchanged primarily between the
bouncing drop and the Faraday wave field excited by its impact (Moláček & Bush 2013a,b).
For inviscid gravity-capillary waves, kinetic energy is exchanged with gravitational
potential and surface energies. For our investigation, it thus suffices to characterise thesystem energetics in terms of the wave energy, which is computed readily from our model.
Owing to the slow spatial decay of the walker wave field (2.1b), the standard wave
energy integral diverges when integrating over the plane. Instead, we adopt the notion ofwave intensity E, as coined by Hubert et al. (2022), which acts as a suitable proxy for thesum of the wave field gravitational potential and surface energies. In dimensional units,
973 A4-17

N. Liu, M. Durey and J.W.M. Bush
we thus define
1
E = lim
R→∞ R



1
ρgh2 dx +|x|≤R 2


|x|≤R


σ
1 + |∇h|2 − 1


dx ,

(4.1)

where σ and ρ are the fluid surface tension coefficient and density, respectively, and E
has units of energy per unit length. For small wave slope, we henceforth approximate
the bracketed term in the second integral by 12 |∇h|2 (see Appendix E). By applying
the divergence theorem to the surface energy contribution and exploiting the fact thatthe wave field is monochromatic with wavenumber kF , we demonstrate in § E.1 that thecontributions from gravitational and surface energies are proportional to each other, and(4.1) reduces to12
h2 (x, t) dx,(4.2)
E = (ρg + σ kF ) lim
R→∞ 2R |x|≤Ran expression proportional to that obtained by Labousse et al. (2016b) and Hubert et al.
(2022), who neglected the contribution of surface tension.
To analyse the wave intensity, we transform (4.2) to dimensionless variables bydefining ĥ = h/h0 and Ê = E/E0 , where h0 = ATW /TF and E0 = h20 kF−1 (ρg + σ kF2 ) are
the characteristic wave height and intensity. By once again taking kF−1 and TW as the unitsof length and time (see § 2.1), we define the dimensionless wave intensity by1Ê = limĥ2 (x, t) dx,(4.3)
R→∞ 2R |x|≤R
where

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Pressĥ(x, t) = t−∞

J0 (|x − xp (s)|) e−μ(t−s) ds(4.4)

is the dimensionless form of the pilot wave.
For circular orbital motion, the dimensionless wave intensity takes the remarkably
simple form (see § E.2)
1
(4.5)
Ê = I0 (0),μwhere I0 (0) (defined in (2.7a)) is the amplitude of the wave field beneath the droplet(Durey & Bush 2021). Since the wave intensity and droplet gravitational potential energy
are proportional to each other for orbital motion, either quantity serves equally well as adiagnostic measure of orbital stability. We may now use the tangential force balance (2.13)
to deduce that the orbital wave intensity (4.5) reduces to
Ê =

U2
1
1
−
2
μ2

,

(4.6)

where U = r0 ω is the orbital speed. For fixed wave decay rate μ, we conclude that the wave
intensity is smaller for faster orbiters. Moreover, the transition from stationary bouncing(U = 0) to orbiting (U > 0) serves to decrease the wave energy, as does the transition frombouncing to rectilinear walking (Durey & Milewski 2017). Finally, combining the boundÊ > 0√with (4.6) supplies the upper bound on the orbital speed noted in § 2.2, namely
U < 2.
To explore the connection between wave intensity (or, equivalently, droplet gravitational
potential energy) and orbital stability, we present the dependence of the wave intensity on
973 A4-18

The onset of orbital instability

E/EW

(a)

E/EW

(c)

(b)

2.5

2.5

2.0

2.0

1.5

1.5

1.0

1.0

0.5

0.5

0

1

2

3

2.5

(d )

0

2.0

1.5

1.5

1.0

1.0

0.5

0.5
1

2

r0 /λF

3

2

3

2

3

2.5

2.0

0

1

0

1

r0 /λF

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

Figure 6. The dependence of the wave intensity E (see (4.6)), relative to that of a walker, EW , on the orbital
radius for κ0 = 1.6 and Γ values (a) 0.7, (b) 0.75, (c) 0.8, and (d) 0.85. The colour scheme is the same as infigure 2. The onsets of monotonic and wobbling instabilities are correlated with, respectively, the local maxima
and minima of the wave intensity E.

the orbital radius in figure 6. As the wave decay rate μ is decreased, we observe that
all monotonic instabilities arise in the vicinity of radii that maximise the wave intensity,
while wobbling instabilities generally arise close to the orbital radii that minimise thewave intensity. (We note that the other wobbling instabilities appear as side bands to thefirst monotonic instability tongue; see figure 2(b).) This correlation is thus indicative ofan underlying orbital energy principle, according to which the magnitude of the wave
intensity prescribes the stability of the corresponding circular orbit (Durey 2018). To test
this hypothesis, we present in figure 7 the orbital radius arising at the tip of each of the first
four wobbling (green curves) and monotonic (red curves) instability tongues, i.e. the radiithat locally minimise Γ along the stability boundary (the white curve in figure 2b). For
each critical value of Γ (or μ = 1 − Γ ), we compare the corresponding critical orbital
radius to the extrema of the wave field intensity (blue curves), confirming our observationthat maxima and minima of E correspond approximately to the tips of the monotonic andwobbling instability tongues, respectively, with the agreement improving for larger orbital
radius.
To further assess the utility of different energy-based heuristics, we also compare thecritical radii to the zeros of J0 (r0 ), J1 (r0 ) and J2 (r0 ), as suggested by Labousse et al.
(2016a). As presented in figure 7, the monotonic instabilities appear at radii slightly
smaller than the zeros of J1 (r0 ), while the wobbling instabilities typically align closelywith the zeros of J2 (r0 ). We note, however, that the critical radii exhibit a weak dependenceon κ0 that is not captured by the zeros of Bessel functions. On the other hand, the zerosof J1 (r0 ) approach the radii satisfying sin(2r0 ) = 1 for larger orbital radius, consistent
with our asymptotic analysis of the monotonic instability (§ 3.2). Likewise, the zeros ofJ0 (r0 ) and J2 (r0 ) both approach the radii satisfying sin(2r0 ) = −1, in agreement with our
973 A4-19

N. Liu, M. Durey and J.W.M. Bush
(a)

(b)
1.13

r0/λF

0.62

1.12

0.60

(d )

1.63

2.13

1.62

2.12

1.11

0.58
0.5 1.0 1.5 2.0 2.5(e) 0.88
r0/λF

(c)

(f)

0.5 1.0 1.5 2.0 2.5

(g)

1.37
1.35

0.82
0.5 1.0 1.5 2.0 2.5

κ0

(h)

1.87

0.86
0.84

1.61
0.5 1.0 1.5 2.0 2.5

0.5 1.0 1.5 2.0 2.5

2.37

1.86
2.36

1.85

1.33
0.5 1.0 1.5 2.0 2.5

κ0

0.5 1.0 1.5 2.0 2.5

κ0

0.5 1.0 1.5 2.0 2.5

κ0https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

Figure 7. The dependence of the critical orbital radius on the inertia-to-drag ratio, κ0 = m/DTW , for (a–d)
monotonic (red) and (e–h) wobbling (green) instabilities. Each plot corresponds to a different instability
tongue in figure 2(b), with radius increasing from left to right. We compare the critical radii of monotonic andwobbling instabilities to those that maximise and minimise the wave intensity ((4.5), blue) and the derivativeof the mean wave force (orange, § 4.2), respectively. We also compare the critical radii with those predictedby the asymptotic formulae for monotonic ((3.8) and (3.12)) and wobbling ((3.9) and (3.10a,b)) instabilities(purple). The black horizontal lines correspond to zeros of J0 (kF r0 ) (dotted), J1 (kF r0 ) (dashed) and J2 (kF r0 )(dot-dashed), whose relation to orbital stability was noted by Labousse et al. (2016a). With increasing orbital
radius, the range of each vertical axis narrows, reflecting the improved prediction of each heuristic.

analysis of wobbling instabilities (§ 3.2). We thus rationalise the Labousse et al. (2016a)
conjecture that the zeros of Bessel functions are likely to identify the loci of the most
unstable orbital radii.
Finally, we compare the critical orbital radii of the numerically computed stability
boundary with our asymptotic formulae derived in §§ 3.2 and 3.3. Specifically, we
minimise numerically the first-order solutions, corresponding to the combination of (3.8)
and (3.12) for the monotonic instability tongues (figures 7a–d), and (3.9) and (3.10a,b) for
the wobbling instability tongues (figures 7e–h). Our asymptotic results appear as purple
curves in figure 7, and agree favourably with the critical radii of the numerical stability
boundary, generally capturing the correct trend with increasing κ0 , and satisfying the
anticipated O(r0−2 ) convergence as the orbital radius is increased. One limitation of thisapproach, however, is the absence of a local minimum for the first wobbling tongue whenusing the first-order correction (3.10a,b) (owing to the parasitic influence of the cot(2r0 )term). Instead, we compare the critical radius in this case to that computed from theleading-order solution given in (3.9), which explains the larger discrepancy in figure 7(e).
Nevertheless, our large-radius asymptotic results work surprisingly well in this case, where
the orbital radius is relatively small.
In summary, we find that the simple heuristic criteria for the critical orbital radii aremoderately successful, with the agreement improving for larger orbital radius. The zerosof Bessel functions generally give better agreement with the numerical results than theextrema of the wave intensity, or, equivalently, the droplet’s gravitational potential energy.
Notably, incorporating the droplet’s kinetic energy within this latter heuristic does notaffect the critical orbital radii significantly, with the resultant curves generally being
973 A4-20

The onset of orbital instability
indistinguishable from the blue curves in figure 7. Finally, these heuristic arguments arelimited by their inability to predict the critical memory at the onset of instability andto capture the dependence of the critical radii on κ0 . Both quantities may be computed
accurately using our asymptotic framework.

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

4.2. The mean wave field
Although the zeros of J1 (r0 ) and J2 (r0 ) provide satisfactory agreement with the numericalresults for the tip of each instability tongue (Labousse et al. 2016a), this heuristic does
not provide a rationale for the type of instability. We proceed to develop a new dynamic
rationale that is asymptotically equivalent to the heuristic of Labousse et al. (2016a), yet
explains the alternation between monotonic and wobbling instabilities with increasing
orbital radius evident in figure 2(b).
We proceed by developing a dynamical interpretation of orbital instability in terms of
the force applied by the mean wave field, specifically that averaged over one orbital period.
One may decompose the orbital wave field into a continually evolving, non-axisymmetriccomponent that serves to propel the droplet at a constant horizontal speed (Bush, Oza &Moláček 2014; Labousse & Perrard 2014), and a static axisymmetric component (the mean
wave field) that imparts either an inward or outward radial force to the droplet (Labousse
et al. 2014; Durey et al. 2018). We note that a similar decomposition of the wave field
applies when considering small perturbations from orbital motion, for which the mean
wave field may now be regarded as a quasi-static potential. Notably, Perrard et al. (2014b)
used this potential to determine the radii of quantised circular orbits; in contrast, we use
this potential to explain the onset of orbital instability.
The dimensionless axisymmetric mean wave field h̄(r) accompanying a droplet
executing a circular orbit of radius r0 about the origin is given by (Perrard et al. 2014b;
Tambasco & Bush 2018)1h̄(r) = J0 (r0 ) J0 (r).
(4.7)
μ
We thus deduce that the dimensionless radial wave forcedh̄ 
F (r0 ) = −2 
dr r=r0(4.8)

applied by the mean wave field along the droplet’s trajectory is2
2 cos(2r0 )1
F (r0 ) = J0 (r0 ) J1 (r0 ) = −
+O 2
μ
μπr0
r0

with derivative
4 sin(2r0 )1
F (r0 ) =
+O 2
μπr0
r0

.

,

(4.9)

(4.10)As shown in § 3.1, wobbling instabilities occur when sin(2r0 ) < 0, whereas monotonic
instabilities occur when sin(2r0 ) > 0. Equation (4.10) suggests that for sufficiently large
r0 , the type of instability exhibited by increasing the memory at constant r0 is thus
related to the derivative of the mean wave force. Specifically, if F (r0 ) < 0, then thecorresponding circular orbit destabilises via a wobbling instability. Conversely, if F (r0 ) >0, then the circular orbit destabilises via a monotonic instability.
973 A4-21

N. Liu, M. Durey and J.W.M. Bush0.04

J0(kF r0) J1(kF r0)

0.02

0

–0.02

–0.04

1.8

2.1

2.4

r0 /λF

2.7

3.0

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

Figure 8. Schematic representation of orbital instability in terms of the mean wave force. We present thedependence of 12 μ F (r0 ) = J0 (kF r0 ) J1 (kF r0 ), which is proportional to the mean wave force, on the orbital
radius for κ0 = 1.6. The curve is colour-coded according to the type of instability arising along the stability
boundary (white curve in figure 2b). Along green and red portions of the curve, the mean wave field acts likean attractive and repulsive spring, respectively, giving rise to wobbling and monotonic instabilities. The sign ofthe mean wave force is denoted by arrows, and its slope is maximised and minimised at the yellow and purple
squares, respectively.

This correlation between the sign of F (r0 ) and the form of instability can be interpreted
physically through consideration of figure 8. As in our linear stability analysis (§ 2.3),
we posit that the instantaneous orbital radius rp (t) has the form rp (t) = r0 + r1 (t),where 0 <  1 is a small parameter, and r1 (t) denotes the perturbation to the orbital
radius. It follows that the force exerted by the mean wave field may be approximatedby F (rp ) ≈ F (r0 ) + r1 F (r0 ); thus the direction and magnitude of the perturbed wave
force is prescribed by F (r0 ). When F (r0 ) > 0, an outward radial perturbation (r1 > 0)
results in an increase in the outward force that drives the droplet away from equilibrium.
Similarly, an inward perturbation (r1 < 0) decreases the outward force. When F (r0 ) > 0,
we may thus regard the mean wave force as a repulsive spring force that induces monotonic
changes in the orbital radius. Conversely, when F (r0 ) < 0, the mean wave force behaves
like an attractive spring that opposes any perturbations in the orbital radius from r0 andso induces oscillations in the orbital radius. This physical picture is consistent with theobservation that the type of instability changes when F (r0 ) changes sign (see figure 8).
Moreover, the mean wave field’s opposition to perturbations from r0 hinders the onset ofthe corresponding wobbling instabilities, which is consistent with wobbling instabilities
occurring at higher memory than their monotonic counterparts. Finally, the magnitudeof the radial wave force increases as μ decreases, thereby increasing the sensitivity toperturbation at higher memory.
In an attempt to rationalise the radii of the most unstable circular orbits (corresponding
to the tips of the stability tongues in figure 2b), we follow the above argument andposit that these will correspond to the radii marked by the largest relative change in theperturbation force, as characterised by the coefficient F (r0 ). The coefficient F (r0 ) is
maximised in magnitude at critical radii rc satisfying F (rc ) = 0. Specifically, maxima
in F (r0 ) correspond to monotonic instabilities, and minima in F (r0 ) correspond to
973 A4-22

The onset of orbital instability
wobbling instabilities. In figure 7, we observe that the most unstable radii predicted bythis heuristic (F (rc ) = 0) exhibit excellent quantitative agreement with the numericalsolution for monotonic instabilities, and a favourable agreement for wobbling instabilities.
The limitation of this heuristic for wobbling instabilities is presumably rooted in small
deviations of the mean wave field from a quasi-static potential for oscillatory droplet
motion.
In summary, we have developed a new rationale for the onset of wobbling andmonotonic instabilities in terms of the force exerted on the droplet by the mean wave field.
Specifically, our investigation indicates that the critical radii rc at the onset of instability
approximately satisfy F (rc ) = 0. In contrast to the heuristic arguments presented in§ 4.1, our rationale explains the alternation of wobbling and monotonic instabilities
with increasing orbital radius. Moreover, our study suggests that the oscillatory andquasi-periodic nature of the stability boundary (see figure 2b) is correlated with thequasi-monochromatic mean wave field. As a caveat, our rationale is valid only forcircular orbits near the stability boundary, and so cannot be used to differentiate between
stable and unstable orbits. Like the heuristic arguments presented in § 4.1, our rationale
neither predicts the critical memory of instability nor accounts for the dependence on theparameter κ0 . Nevertheless, it performs remarkably well over the range of values of theinertia-to-drag parameter κ0 accessible in the laboratory.

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

5. Discussion

We have examined the stability of circular inertial orbits executed by droplets walking
in a rotating frame. We have developed an asymptotic framework for studying orbital
stability, specifically for characterising the critical memory at the onset of the wobbling
and monotonic instabilities. Our asymptotic model rationalises the repeating structureof the stability diagram in figure 2, which we demonstrate to be rooted in the periodic
nature of the mean wave field. Our theory not only rationalises the preponderance of
2-wobble instabilities, but also predicts that the destabilisation frequency is exactly twice
the angular frequency at the most unstable orbital radii, corresponding to the tips of thestability tongues. In rationalising the most unstable wobbling frequencies reported in theexperiments of Harris & Bush (2014) and simulations of Oza et al. (2014b), our study has
demonstrated the importance of resonant instability in orbital pilot-wave dynamics.
Our asymptotic results show that orbital instability is enhanced through resonance.
Resonant instabilities have been reported in other pilot-wave systems, including the in-line
speed oscillations of the free walker (Bacot et al. 2019; Durey et al. 2020b), 4ω orbital
instabilities in a linear central force (Kurianski, Oza & Bush 2017), and 2ω, 3ω and 4ω
instabilities for pairs of orbiting droplets (Oza et al. 2017; Durey 2018). This naturally
raises the question of the prevalence of resonant instabilities in other pilot-wave systems. In
subsequent work, we will consider the more complex orbits arising for a droplet in a linear
central force, where richer resonances are expected to arise. We expect the asymptotic
framework developed here to be well-suited to addressing this class of problems.
It is worth enumerating the limitations of our study. In our study, we have assumed
resonant bouncing at constant bouncing phase, and neglected the far-field exponential
decay of the wave field (Moláček & Bush 2013b; Turton, Couchman & Bush 2018; Tadrist
et al. 2018). Both of these approximations are known to break down for multi-droplet
systems (Oza et al. 2017; Arbelaiz, Oza & Bush 2018; Galeano-Rios et al. 2018;
Couchman, Turton & Bush 2019), and so might become significant at sufficiently highorbital memory. In addition, we have restricted our attention to the linear stabilityof circular orbits, and so not considered nonlinear effects such as the jump up/down
973 A4-23

N. Liu, M. Durey and J.W.M. Bushinstabilities reported by Harris & Bush (2014) and Oza et al. (2014b). Finally, we havenot considered instability in the large-radius limit inaccessible in the laboratory, which
one expects to be related to the instability of the free, rectilinear walking state (Bacot et al.
2019; Hubert et al. 2019; Durey et al. 2020b), a matter to be considered elsewhere.
We have also compared the success of various heuristics for predicting the onset oforbital instability in our system. Our study has lead us to introduce two new heuristics,
specifically the wave intensity and the mean wave force. We show that the orbital radii
corresponding to the onset of monotonic and wobbling instabilities generally arise near
the radii that maximise and minimise the wave intensity, respectively. We have also
demonstrated the equivalence of the wave intensity and wave height beneath the dropas proxies for the droplet energy; specifically, the wave energy is proportional to thegravitational potential energy of the drop. This result provides new insight into theobservations made by Couchman & Bush (2020) and Thomson et al. (2020) that a ring ofdroplets rearranges itself so as to minimise the mean gravitational potential of the droplets,suggesting that they are doing so in order to minimise the global wave energy.
Labousse et al. (2016a) proposed that wobbling and monotonic instabilities arise fromorbits that receive little wave energy from the J1 and J2 modes, and postulated that the most
unstable orbital radii occur at zeros of these Bessel functions. We have demonstrated that
this heuristic leads to impressive agreement with our numerical calculations over a wide
range of parameter values, and outperforms our wave intensity extremisation principle,
defined in § 4.1. However, the Labousse et al. (2016a) heuristic does not provide a physicalmechanism that distinguishes between the wobbling and monotonic instabilities. Our new
heuristic based on the form of the mean wave field (§ 4.2) indicates the correlation between
the type of instability and the derivative of the mean wave force, and so sheds light on thealternation between wobbling and monotonic instabilities with increasing orbital radius.
The general utility of mean-pilot-wave-based heuristics will be considered elsewhere.
Acknowledgements. The authors thank A. Oza and D. Harris for providing the images presented in figure 1
and the data presented in figure 4.
Funding. We gratefully acknowledge the generous financial support of the NSF through grantCMMI-2154151.

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

Declaration of interests. The authors report no conflict of interest.
Data availability statement. The data that support the findings of this study are available from thecorresponding author, J.W.M.B., upon reasonable request.
Author ORCIDs.
Nicholas Liu https://orcid.org/0000-0002-8962-8536;
Matthew Durey https://orcid.org/0000-0002-4232-1705;
John W.M. Bush https://orcid.org/0000-0002-7936-7256.
Author contributions. J.W.M.B. and M.D. conceived the project. N.L. was responsible for the mathematicaland numerical developments, which were overseen by M.D. All authors contributed to the writing.

Appendix A. Physical parameters

We here define the Faraday wavenumber kF , drag factor D, and wave amplitude parameter
A, in terms of physical quantities (Oza et al. 2013, 2014a). Specifically, the Faraday
wavenumber satisfies the gravity-capillary dispersion relation (Benjamin & Ursell 1954)

σ k3
(A1)
(πf )2 = gkF + F tanh(kF H),
ρ
973 A4-24

The onset of orbital instability
where σ is the coefficient of surface tension, ρ is the fluid density, and H is the fluid depth.
Furthermore, we define the drag factorρair gRD
ρRD
D = 0.17 mg+ 6πμair RD 1 +,(A2)σ
12μair f
and wave amplitude
2 1/2

RD kF νeff
1
kF RD
A= √
mgTF sin Φ,√22
2π 3kF RD + Bo σ TF

(A3)

where RD is the droplet radius, μair and ρair are the dynamic viscosity and density of
air, respectively, νeff is the effective fluid kinematic viscosity (Moláček & Bush 2013b),
Bo = ρgR2D /σ is the Bond number, and sin Φ = 0.2 is the droplet impact phase (Oza et al.
2013, 2014b).
Appendix B. Evaluation of stability integrals

We proceed to evaluate analytically the stability integrals defined in (2.7), where weleverage the relationship between the stability integrals (2.9a,b) to evaluate Sm and Cm in
terms of Im only. The evaluation of each stability integral hinges on the observation that
J2m (2r0 sin(x)) is even and π-periodic, and so may be expressed in terms of its bi-infinite
cosine series (derived from Gradshteyn & Ryzhik (2014, § 6.681)) as
J2m (2r0 sin(x)) =∞(−1)n Jm+n (r0 ) Jm−n (r0 ) cos(2nx).

(B1)

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Pressn=−∞

By performing the change of integration variable t → ωt in (2.7a) and using the cosine
series expansion (B1), the integral Im (s) may be recast as the infinite sum
 ∞
∞
1 
n
Im (s) =(−1) Jm+n (r0 ) Jm−n (r0 )cos(nt) e−ηt dt,(B2)ω n=−∞
0where η = (μ + s)/ω and Re(η) > 0. By integrating each of the resultant integrals
analytically, we deduce that the stability integral is equivalent to

Im (s) =∞η  (−1)n Jm+n (r0 ) Jm−n (r0 ).
ω n=−∞
η2 + n2(B3)This formulation is convenient for numerical computation of the stability integrals, andis used in the construction of orbital stability diagrams and the tracking of orbital
stability boundaries (e.g. figure 2). However, further analytical progress can be made byusing a partial fraction expansion in (B3) and applying the Lerche–Newberger sum rule
(Newberger 1982; Bakker & Temme 1984), giving rise to the closed-form expression givenin (2.8), which forms the basis of our analytical developments.
Notably, the sum (or difference) appearing in the Bessel function orders in (2.8) whenm = 1 complicates the application of large-argument asymptotic expansions in our study.
As such, we derive two Bessel function identities that assist with the elimination of thesum in the order, thereby allowing us to derive a set of convenient equations for further
973 A4-25

N. Liu, M. Durey and J.W.M. Bushanalysis. To proceed, we first recall the Bessel function Wronskian, recurrence relation andderivative expression given by (Abramowitz & Stegun 1948)
2 sin(νπ)
Jν (x) J1−ν (x) + J−ν (x) Jν−1 (x) =
,
(B4a)
πx
2ν
Jν−1 (x) + Jν+1 (x) =(B4b)Jν (x),
xJν−1 (x) − Jν+1 (x) = 2 Jν (x),(B4c)
respectively, which form the basis of our forthcoming algebraic manipulations. To
derive the first identity, we apply the recurrence relation (B4b), add and subtract(2ν/x) J−ν (x) Jν−1 (x) to and from the resultant expression, and simplify using theWronskian (B4a) and derivative relation (B4c) to obtain4ν sin(νπ) 4νJ1−ν (x) J−1+ν (x) − J1+ν (x) J−1−ν (x) =−(B5)J−ν (x) Jν (x).
xπx2
To derive the second identity, we differentiate the Wronskian identity (B4a) and apply therecurrence relation (B4b) to find that2ν sin(πν)
J1−ν (x) Jν (x) + J1+ν (x) J−ν (x) = −2 J−ν (x) Jν (x) −
.
(B6)
πx2
We now use the relations (B5) and (B6) in conjunction with (2.8) and (2.9a,b) to evaluate
the integrals and integral combinations arising in the stability problem (2.6), yielding
π csch(πβ)
I0 (0) =(B7a)Jiβ (r0 ) J−iβ (r0 ),ω−2π csch(πβ) i sinh(πβ)∂ I0 (0)=− J−iβ (r0 ) Jiβ (r0 ) ,
(B7b)
∂r0
ω
πr0

2πη csch(πη) i sinh(πη) J−iη (r0 ) Jiη (r0 ),(B7c)
S0 (s) =
−
ω
r0
πr02
2iηπ csch(πη) i sinh(πη)− iη Jiη (r0 ) J−iη (r0 ) ,
π
r02 ω

2π csch(πη) η sinh(πη)
C0 (s) + I1 (s) =+ Jiη (r0 ) J−iη (r0 ) ,ωπr02

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

C0 (s) − I1 (s) = −

(B7d)

(B7e)

where we recall that η = (μ + s)/ω, and β = μ/ω is the reciprocal of the orbital memory.
Appendix C. Stationary phase point contributions to stability integrals

We seek the contributions at each critical point to the integral Lm (ξ ), defined in (2.15),
for which we consider r0  1 and ξ = O(1). We first consider the contribution to Lm (ξ )near the edges of the integration region, namely θ = 0 and θ = 2π, which we denote by
Lm,0 (ξ ) and Lm,2π (ξ ), respectively. To determine Lm,0 (ξ ), we approximate sin(θ/2) ≈θ/2 in the argument of the Bessel function, and recast the integration region as ∞ ∞Lm,0 (ξ ) =J2m (r0 θ) e−(β+iξ )θ dθ −J2m (r0 θ) e−(β+iξ )θ dθ.
(C1)
0

2π

The first integral may be evaluated analytically (Abramowitz & Stegun 1948), and it
remains now to estimate the size of the second integral. We approximate the integrand973 A4-26

The onset of orbital instability
by utilising the large-argument Bessel function expansion (Abramowitz & Stegun 1948)
and noting that θ ≥ 2π across the integration domain; then, by evaluating the resultant
integral analytically, we obtain the approximate bound
 ∞



 ∞

e−2πβ
,
√
r0 r0
2π
2π
(C2)
where ϕm = 14 π(1 + 4m) is a constant phase shift. Provided that β is of size at
most O(r0 ), the analytical evaluation of the first integral in (C1) determines that
this term is the dominant contribution to Lm,0 . By neglecting the contribution fromthe second integral in (C1), we deduce that Lm,0 = O(r0−1 ). The periodicity of the
non-exponential portion of the integrand indicates that the contribution near θ = 2πsatisfies Lm,2π (ξ ) = Lm,0 (ξ ) e−2π(β+iξ ) . The factor e−2πβ ensures that the contribution
Lm,2π (ξ ) is negligible relative to Lm,0 (ξ ) and may henceforth be neglected.
We now determine the contributions arising near points of stationary phase. The internal
contribution is localised about θ = π, and we denote this contribution by Lm,π (ξ ).
To proceed, we deform the integration region to a small region about θ = π, namely
π − δ < θ < π + δ with 0 < δ  1, and define
J2m (r0 θ ) e−(β+iξ )θ

dθ ∼Lm,π (ξ ) =

1
π2 r0

cos(r0 θ − ϕm ) e−(β+iξ )θ dθ = O

 π+δ
π−δ

J2m 2r0 sinθ2e−(β+iξ )θ dθ.

(C3)

By applying the large-argument expansion for the Bessel function (Abramowitz & Stegun1948) and utilising the Taylor expansion sin(θ/2) ≈ 1 − 18 (θ − π)2 , we determine theapproximation

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Presse−(β+iξ )π  ±iϕm
Lm,π (ξ ) ∼ √e
2 πr0 ±

 π+δ
π−δ

e



∓i2r0 1− 81 (θ −π)2

dθ.

(C4)

The leading-order form of the resultant integral may be determined using method ofstationary phase for r0  1, giving rise to the contribution

Lm,π (ξ ) ∼ −2e−π(β+iξ )cos (2r0 + ϕm ) .
r0

(C5)

Appendix D. Asymptotic expansions for the wobbling regime

Along stability boundaries (with s = iξ ω), the stability integrals (2.7) are all evaluatedin terms of Bessel functions of complex order ±iζ and real argument r0 , where ζ ∈{β, β + iξ } (see (B7)). For wobbling and monotonic instabilities, corresponding to β =O(ln(r0 )) and ξ = O(1), the order of the Bessel function is much smaller than theargument. Thus we may expand each of the Bessel function products in (B7) by applying
the large-argument asymptotic expansion to each Bessel function (Abramowitz & Stegun
973 A4-27

N. Liu, M. Durey and J.W.M. Bush1948). Direct evaluation by Mathematica yields
cos(2r0 ) + i sinh(ζ π) cosh(ζ π) (4ζ 2 − 1) sin(2r0 )
−
+
πr0
2πr02
4πr02

3(1 + 4ζ 2 ) cosh(ζ π)
1
(D1a)
+
+O 3 ,4
16πr0r0
sin(2r0 ) + cosh(ζ π) 1 + 4ζ 2 cosh(ζ π)1Jiζ (r0 ) J−iζ (r0 ) =−+ cos(2r0 ) + O 3 ,
2
πr0
2r0
4πr0
r0
(D1b)
J−iζ (r0 ) Jiζ (r0 ) =Jiζ (r0 ) J−iζ (r0 ) =− sin(2r0 ) + cosh(ζ π) (3 + 4ζ 2 ) cosh(ζ π)
+
πr0
8πr03

(−3 + 4ζ 2 ) cos(2r0 )1++O 3 ,
2
4πr0
r0

(D1c)

where all terms necessary for the dominant balance eπβ = O(r02 ) have been retained,
and the expansions (D1) are valid when |ζ 2 + 14 |  r0 (Abramowitz & Stegun 1948).
By substituting these expansions into the analytical expressions for each of the stability
integrals (see (B7)), we deduce that the stability coefficients (see (2.12)) have theasymptotic formA =−2 sin(2r0 )(4β 2 − 7) cos(2r0 ) csch(πβ)(csch(π(β + iξ )) + csch(πβ)) +U
2Ur03 − ξ 2 κ0 ξ 2 U 2 cos(2r0 ) csch(π(β + iξ ))(−3 + 4β 2 + 8iβξ − 4ξ 2 )
−
+
2Ur0
Ur02
r02

1
+O 4 ,
r0

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

+

(D2a)
iκ0 ξ U 22iξ
2 cos(2r0 ) csch(πβ) (iξ − β)
−
+
Ur0
r02
Ur02


1
+ csch(π(β + iξ )) (β + iξ ) + O 4 ,r02U 2(β + iξ ) iκ0 ξ U 2
1
C =
−
+
+O 3 ,
2
2
r0
Ur0
r0
r0

2iξ U 2iξ(β + iξ ) κ0 ξ 2 U 2
1
D=−−+O 3 ,
2
2
r0
Ur0
r0
r0

B=

where we have again utilised the dominant balance eπβ = O(r02 ).
973 A4-28

(D2b)

(D2c)

(D2d)

The onset of orbital instability
Appendix E. Wave intensity

E.1. Surface energy contribution
We proceed to show that the contribution of surface energy to the wave intensity is
proportional to that of gravitational potential energy. For small wave slopes, we define
the contribution of surface energy as (in dimensional units)1
σ|∇h|2 dx.
(E1)
ES = lim
R→∞ R |x|≤R 2
By writing |∇h|2 = ∇ · (h ∇h) − h Δh in (E1), where Δ denotes the Laplacian operator,
and applying the divergence theorem to the first term, we obtain
 


1 σ
σ
ES = limh ∇h · n ds −h Δh dx ,(E2)R→∞ R 2 |x|=R|x|≤R 2where n is the outward-pointing unit vector to the circle |x| = R, and ds denotes line
integration. Owing to the far-field decay of the wave field, we see that the boundary integralin (E2) vanishes in the limit R → ∞, leaving
σ1
ES = lim− h Δh dx.
(E3)
R→∞ R |x|≤R
2
Finally, we recall that the wave field is composed of monochromatic Bessel functions withwavenumber kF . It follows that Δh = −kF2 h for our wave model, which allows us to reduce(E3) to the simplified form
1
ES = limR→∞ Rσ kF2 2h (x, t) dx,|x|≤R 2(E4)https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

which is proportional to the contribution from the gravitational potential energy.
E.2. Wave intensity for orbital motion
To compute the wave intensity for a given droplet trajectory, we first introduce the polar
coordinate system x = r(cos θ, sin θ) and decompose the dimensionless wave field (see(4.4)) asĥ(x, t) =∞an (t) Jn (r) einθ .

(E5)

n=−∞

By substituting (E5) into (4.4) and applying Graf’s addition theorem (Abramowitz &Stegun 1948), we find that each coefficient an (t) is defined by the integral (Labousse et al.
2014)
 t
an (t) =Jn (rp (s)) e−inθp (s) e−μ(t−s) ds,(E6)−∞

where xp (t) = rp (t) (cos θp (t), sin θp (t)) is the droplet position. Finally, by substituting(E5) into (4.3), applying Parseval’s theorem in the azimuthal direction, evaluating the973 A4-29

N. Liu, M. Durey and J.W.M. Bushintegral in the radial direction, and taking the limit R → ∞, we obtain (Labousse et al.
2016b; Hubert et al. 2022)∞
Ê(t) =|an (t)|2 .
(E7)
n=−∞

We proceed to evaluate the wave intensity for the wave field accompanying a droplet
orbiting the origin. By substituting rp (t) = r0 and θp (t) = ωt into (E6) and evaluating theresultant integral, we deduce thatan (t) =Jn (r0 )e−inωt.
μ − inω(E8)Upon substituting an (t) into (E7), we thus deduce
Ê =

∞


J2n (r0 ),μ2 + n2 ω2
n=−∞

(E9)

whereupon comparing this equation with (B3) determines the expression given in (4.5).

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

R EFERENCES
A BRAMOWITZ , M. & S TEGUN , I.A. 1948 Handbook of Mathematical Functions with Formulas, Graphs, andMathematical Tables, vol. 55. US Government Printing Office.
A RBELAIZ , J., O ZA , A.U. & B USH , J.W.M. 2018 Promenading pairs of walking droplets: dynamics andstability. Phys. Rev. Fluids 3, 013604.
BACOT, V., P ERRARD , S., L ABOUSSE , M., C OUDER , Y. & F ORT, E. 2019 Multistable free states of an active
particle from a coherent memory dynamics. Phys. Rev. Lett. 122, 104303.
BAKKER , M. & T EMME , N.M. 1984 Sum rule for products of Bessel functions: comments on a paper byNewberger. J. Math. Phys. 25 (5), 1266–1267.
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.
B LEISTEIN , N. & H ANDELSMAN , R.A. 1975 Asymptotic Expansions of Integrals. Ardent Media.
BORGHESI , C., M OUKHTAR , J., L ABOUSSE , M., E DDI , A., F ORT, E. & C OUDER , Y. 2014 Interaction of
two walkers: wave-mediated energy and force. Phys. Rev. E 90, 063017.
DE B ROGLIE , L. 1926 Ondes et Mouvements. Gauthier-Villars.
DE B ROGLIE , L. 1930 An Introduction to the Study of Wave Mechanics. Methuen & Co.
B USH , J.W.M. 2015 Pilot-wave hydrodynamics. Annu. Rev. Fluid Mech. 47, 269–292.
B USH , J.W.M. & O ZA , A.U. 2020 Hydrodynamic quantum analogs. Rep. Prog. Phys. 84, 017001.
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 OUCHMAN , M.M.P. & B USH , J.W.M. 2020 Free rings of bouncing droplets: stability and dynamics. J. Fluid
Mech. 903, A49.
C OUCHMAN , M.M.P., T URTON , S.E. & B USH , J.W.M. 2019 Bouncing phase variations in pilot-wavehydrodynamics and the stability of droplet pairs. J. Fluid Mech. 871, 212–243.
C OUDER , Y. & F ORT, E. 2006 Single particle diffraction and interference at a macroscopic scale. Phys. Rev.
Lett. 97, 154101.
C OUDER , Y., F ORT, E., GAUTIER , C.-H. & BOUDAOUD , A. 2005a From bouncing to floating:noncoalescence of drops on a fluid bath. Phys. Rev. Lett. 94, 177801.
C OUDER , Y., P ROTIÈRE , S., F ORT, E. & BOUDAOUD , A. 2005b Walking and orbiting droplets. Nature 437,
208.
C RISTEA-P LATON , T., S ÁENZ , P.J. & B USH , J.W.M. 2018 Walking droplets in a circular corral: quantisationand chaos. Chaos 28, 096116.
D ELVES , L.M. & LYNESS , J.N. 1967 A numerical method for locating the zeros of an analytic function. Math.
Comput. 21 (100), 543–560.
D UREY, M. 2018 Faraday wave-droplet dynamics: a hydrodynamic quantum analogue. PhD thesis, Universityof Bath.

973 A4-30

https://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

The onset of orbital instability
D UREY, M. & B USH , J.W.M. 2021 Classical pilot-wave dynamics: the free particle. Chaos 31, 033136.
D UREY, M. & M ILEWSKI , P.A. 2017 Faraday wave-droplet dynamics: discrete-time analysis. J. Fluid Mech.
821, 296–329.
D UREY, M., M ILEWSKI , P.A. & B USH , J.W.M. 2018 Dynamics, emergent statistics and the mean-pilot-wave
potential of walking droplets. Chaos 28, 096108.
D UREY, M., M ILEWSKI , P.A. & WANG , Z. 2020a Faraday pilot-wave dynamics in a circular corral. J. Fluid
Mech. 891, A3.
D UREY, M., T URTON , S.E. & B USH , J.W.M. 2020b Speed oscillations in classical pilot-wave dynamics.
Proc. R. Soc. A 476, 20190884.
E DDI , A., F ORT, E., M OISY, F. & C OUDER , Y. 2009 Unpredictable tunneling of a classical wave–particle
association. Phys. Rev. Lett. 102, 240401.
E DDI , A., S ULTAN , E., M OUKHTAR , J., F ORT, E., ROSSI , M. & C OUDER , Y. 2011 Information stored in
Faraday waves: the origin of a path memory. J. Fluid Mech. 674, 433–463.
E LLEGAARD , C. & L EVINSEN , M.T. 2020 Interaction of wave-driven particles with slit structures. Phys. Rev.
E 102, 023115.
F ORT, E., E DDI , A., M OUKHTAR , J., BOUDAOUD , A. & C OUDER , Y. 2010 Path-memory induced
quantization of classical orbits. Proc. Natl Acad. Sci. 107 (41), 17515–17520.
GALEANO -R IOS , C.A., C OUCHMAN , M.M.P., CALDAIROU, P. & B USH , J.W.M. 2018 Ratcheting droplet
pairs. Chaos 28, 096112.
G RADSHTEYN , I.S. & RYZHIK , I.M. 2014 Table of Integrals, Series, and Products. Academic Press.
H ARRIS , D.M. & B USH , J.W.M. 2014 Drops walking in a rotating frame: from quantized orbits to multimodal
statistics. J. Fluid Mech. 739, 444–464.
H ARRIS , D.M., M OUKHTAR , J., F ORT, E., C OUDER , Y. & B USH , J.W.M. 2013 Wavelike statistics frompilot-wave dynamics in a circular corral. Phys. Rev. E 88, 011001.
H UBERT, M., L ABOUSSE , M., P ERRARD , S., L ABOUSSE , M., VANDEWALLE , N. & C OUDER , Y. 2019Tunable bimodal explorations of space from memory-driven deterministic dynamics. Phys. Rev. E 100,
032201.
H UBERT, M., P ERRARD , S., VANDEWALLE , N. & L ABOUSSE , M. 2022 Overload wave-memory induces
amnesia of a self-propelled particle. Nat. Commun. 13 (1), 4357.
K URIANSKI , K.M., O ZA , A.U. & B USH , J.W.M. 2017 Simulations of pilot-wave dynamics in a simple
harmonic potential. Phys. Rev. Fluids 2, 113602.
L ABOUSSE , M. 2014 Étude d’une dynamique à mémoire de chemin: une expérimentation théorique. PhD
thesis, Université Pierre et Marie Curie UPMC Paris VI.
L ABOUSSE , M., O ZA , A.U., P ERRARD , S. & B USH , J.W.M. 2016a Pilot-wave dynamics in a harmonic
potential: quantization and stability of circular orbits. Phys. Rev. E 93, 033122.
L ABOUSSE , M. & P ERRARD , S. 2014 Non-Hamiltonian features of a classical pilot-wave dynamics. Phys.
Rev. E 90, 022913.
L ABOUSSE , M., P ERRARD , S., C OUDER , Y. & F ORT, E. 2014 Build-up of macroscopic eigenstates in amemory-based constrained system. New J. Phys. 16, 113027.
L ABOUSSE , M., P ERRARD , S., C OUDER , Y. & F ORT, E. 2016b Self-attraction into spinning eigenstates of amobile wave source by its emission back-reaction. Phys. Rev. E 94, 042224.
M ILES , J. & H ENDERSON , D. 1990 Parametrically forced surface waves. Annu. Rev. Fluid Mech. 22, 143–165.
M OLÁ ČEK , J. & B USH , J.W.M. 2013a Droplets bouncing on a vibrating fluid 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.
N EWBERGER , B.S. 1982 New sum rule for products of Bessel functions with application to plasma physics.
J. Math. Phys. 23 (7), 1278–1281.
O ZA , A.U. 2014 A trajectory equation for walking droplets: hydrodynamic pilot-wave theory. PhD thesis,
Massachusetts Institute of Technology.
O ZA , A.U., H ARRIS , D.M., ROSALES , R.R. & B USH , J.W.M. 2014a Pilot-wave dynamics in a 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., ROSALES , R.R. & B USH , J.W.M. 2018 Hydrodynamic spin states. Chaos 28, 096106.
O ZA , A.U., S IÉFERT, E., H ARRIS , D.M., M OLÁ ČEK , J. & B USH , J.W.M. 2017 Orbiting pairs of 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-wavehydrodynamics in a rotating frame: exotic orbits. Phys. Fluids 26, 082101.

973 A4-31

N. Liu, M. Durey and J.W.M. Bushhttps://doi.org/10.1017/jfm.2023.742 Published online by Cambridge University Press

P ERRARD , S., L ABOUSSE , M., F ORT, E. & C OUDER , Y. 2014a Chaos driven by interfering memory. Phys.
Rev. Lett. 113, 104101.
P ERRARD , S., L ABOUSSE , M., M ISKIN , M., F ORT, E. & C OUDER , Y. 2014b Self-organization into quantizedeigenstates of a classical wave-driven particle. Nat. Commun. 5, 3219.
P UCCI , G., H ARRIS , D.M., FARIA , L.M. & B USH , J.W.M. 2018 Walking droplets interacting with single
and double slits. J. Fluid Mech. 835, 1136–1156.
S ÁENZ , P.J., C RISTEA-P LATON , T. & B USH , J.W.M. 2018 Statistical projection effects in a hydrodynamic
pilot-wave system. Nat. Phys. 14, 315–319.
S ÁENZ , P.J., C RISTEA-P LATON , T. & B USH , J.W.M. 2020 A hydrodynamic analog of Friedel oscillations.
Sci. Adv. 6 (20), eaay9234.
S ÁENZ , P.J., P UCCI , G., T URTON , S.E., G OUJON , A., ROSALES , R.R., D UNKEL , J. & B USH , J.W.M.
2021 Emergent order in hydrodynamic spin lattices. Nature 596, 58–62.
TADRIST, L., G ILET, T., S CHLAGHECK , P. & B USH , J.W.M. 2020 Predictability in a hydrodynamic
pilot-wave system: resolution of walker tunneling. Phys. Rev. E 102, 013104.
TADRIST, L., S HIM , J.-B., G ILET, T. & S CHLAGHECK , P. 2018 Faraday instability and subthreshold Faraday
waves: surface waves emitted by walkers. J. Fluid Mech. 848, 906–945.
TAMBASCO , L.D. & B USH , J.W.M. 2018 Exploring orbital dynamics and trapping with a generalized
pilot-wave framework. Chaos 28, 096115.
T HOMSON , S.J., D UREY, M. & ROSALES , R.R. 2020 Collective vibrations of a hydrodynamic active lattice.
Proc. R. Soc. A 476, 20200155.
T URTON , S.E., C OUCHMAN , M.M.P. & B USH , J.W.M. 2018 A review of the theoretical modeling of walking
droplets: towards a generalized pilot-wave framework. Chaos 28, 096111.

973 A4-32