2.1. Governing equations

To model the horizontal motion of the walker, we utilise the stroboscopic trajectory equation developed by Oza et al. (Reference Oza, Rosales and Bush2013, Reference Oza, Harris, Rosales and Bush2014a), whereby the pilot-wave system is time-averaged over one bouncing period, $T_F = 2/f$ (Moláček & Bush Reference Moláček and Bush2013b). The droplet's horizontal position $\boldsymbol {x}_p(t)$ thus evolves over time $t$ according to (Oza et al. Reference Oza, Rosales and Bush2013, Reference Oza, Harris, Rosales and Bush2014a,Reference Oza, Wind-Willassen, Harris, Rosales and Bushb):

(2.1a)\begin{equation} m\ddot{\boldsymbol{x}}_p + D\dot{\boldsymbol{x}}_p ={-}mg\,\boldsymbol{\nabla} h(\boldsymbol{x}_p(t),t) - 2m\boldsymbol{\varOmega}\times\dot{\boldsymbol{x}}_p. \end{equation}

The drop is propelled by the wave force $-mg\,\boldsymbol {\nabla } h({\boldsymbol {x}}_p(t), t)$, and also responds to the linear drag force $-D\dot {\boldsymbol {x}}_p$ and the Coriolis force $-2m\boldsymbol {\varOmega }\times \dot {\boldsymbol {x}}_p$. The accompanying pilot wave,

(2.1b)\begin{equation} h(\boldsymbol{x},t) = \frac{A}{T_F}\int_{-\infty}^t \mathrm{J}_0(k_F\,|\boldsymbol{x} - \boldsymbol{x}_p(s)|)\,\mathrm{e}^{-(t-s)/T_M}\,\mathrm{d}s, \end{equation}

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 $T_M$. 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 $k_F$, $D$ and $A$, are defined in terms of physical quantities in Appendix A.

We project the pilot wave onto the droplet's path to yield an integro-differential trajectory equation for the droplet (Oza et al. Reference Oza, Rosales and Bush2013) that may be expressed in dimensionless variables as (Oza et al. Reference Oza, Harris, Rosales and Bush2014a; Oza, Rosales & Bush Reference Oza, Rosales and Bush2018)

(2.2)\begin{equation} \kappa_0\ddot{\hat{\boldsymbol{x}}}_p + \dot{\hat{\boldsymbol{x}}}_p = 2\int_{-\infty}^{\hat t} \frac{\mathrm{J}_1(|\hat{\boldsymbol{x}}_p(\hat t) - \hat{\boldsymbol{x}}_p(s)|)}{|\hat{\boldsymbol{x}}_p(\hat t) - \hat{\boldsymbol{x}}_p(s)|}\,(\hat{\boldsymbol{x}}_p(\hat t) - \hat{\boldsymbol{x}}_p(s))\,\mathrm{e}^{-\mu(\hat t-s)}\,\mathrm{d}s - \hat{\boldsymbol{\varOmega}}\times\dot{\hat{\boldsymbol{x}}}_p, \end{equation}

where $\hat {\boldsymbol {x}}_p = k_F\boldsymbol {x}_p$, $\hat {t} = t/T_W$, and $T_W = \sqrt {2D T_F/mg A k_F^2}$ is the memory time at the onset of walking, $\gamma = \gamma _W$ (Oza et al. Reference Oza, Rosales and Bush2013; Durey et al. Reference Durey, Turton and Bush2020b). The dimensionless parameters $\mu = T_W/T_M > 0$ and $\kappa _0 = m/D T_W$ describe the wave decay rate and the relative importance of inertial and drag forces, respectively, and $\hat {\boldsymbol {\varOmega }} = 2m\boldsymbol {\varOmega }/D = \hat {\varOmega }\boldsymbol {e}_z$ is the dimensionless rotation vector (see table 1).

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

We characterise the pilot-wave dynamics in terms of the dimensionless vibration parameter $\varGamma = (\gamma - \gamma _W)/(\gamma _F - \gamma _W) = 1-\mu$ (Bush Reference Bush2015; Oza et al. Reference Oza, Rosales and Bush2018; Durey et al. Reference Durey, Turton and Bush2020b), which increases with increasing path memory. We note that $\varGamma = 0$ corresponds to the walking threshold in the absence of a Coriolis force ($\gamma = \gamma _W$), while $\varGamma = 1$ corresponds to the Faraday threshold ($\gamma = \gamma _F$), and thus infinite path memory (Bush Reference Bush2015). The experimental parameter regime of Harris & Bush (Reference Harris and Bush2014) corresponds to $\kappa _0 \approx 1.6$. We note that typically, $\kappa _0$ takes values in the range $0.8 \lesssim \kappa _0 \lesssim 1.6$ in the laboratory; likewise, the dimensionless rotation rate $\hat {\varOmega }$ is restricted to the interval $0 \leq |\hat {\varOmega }| \lesssim 1.3$ (Harris & Bush Reference Harris and Bush2014; Oza et al. Reference Oza, Harris, Rosales and Bush2014a, Reference Oza, Rosales and Bush2018). Henceforth, we thus treat $\kappa _0$ and $\hat {\varOmega }$ as $O(1)$ quantities, whose influence on the pilot-wave dynamics we characterise through systematic asymptotic analysis.

2.2. Orbital dynamics

We characterise orbits in terms of their radius $r_0$ and angular frequency $\omega > 0$. By omitting hats and substituting $\boldsymbol {x}_p(t) = r_0(\cos (\omega t),\sin (\omega t))$ into the trajectory equation (2.2), we express the radial and tangential force balances as (Oza et al. Reference Oza, Harris, Rosales and Bush2014a)

(2.3a)\begin{gather} -\kappa_0 r_0\omega^2 = 2\int_0^\infty \mathrm{J}_1\left(2r_0\sin\left(\frac{\omega t}{2}\right)\right)\sin\left(\frac{\omega t}{2}\right)\mathrm{e}^{-\mu t}\,\mathrm{d}t + \varOmega r_0\omega , \end{gather}

(2.3b)\begin{gather}r_0\omega = 2\int_0^\infty \mathrm{J}_1\left(2r_0\sin\left(\frac{\omega t}{2}\right)\right)\cos\left(\frac{\omega t}{2}\right)\mathrm{e}^{-\mu t}\,\mathrm{d}t, \end{gather}

which may be solved for $r_0$ and $\omega$ given $\kappa _0$, $\mu$ and $\varOmega$.

In our dimensionless notation, a suitable proxy for the orbital memory $M_e^O$ is $\omega _{{orb}} = \omega /\mu$, which is the ratio of the wave decay time scale $\mu ^{-1}$ to the orbital time scale $\omega ^{-1} \sim r_0/u_0$, where $u_0$ is the steady walking speed in the absence of bath rotation (Oza et al. Reference Oza, Rosales and Bush2013; Durey et al. Reference Durey, Turton and Bush2020b). The orbital speed $U = r_0\omega$ typically remains close to the free walking speed $u_0$ and satisfies $U < \sqrt {2}$ for all parameter values (see § 4.1). As $U$ depends only weakly on the orbital radius at fixed memory, we note that $\omega _{{orb}} = U/(r_0\mu )$ decreases with increasing $r_0$.

2.3. Orbital stability

In order to characterise the droplet's response to perturbations from a circular orbit, we apply linear stability analysis. Following the framework developed by Oza et al. (Reference Oza, Harris, Rosales and Bush2014a), we linearise the trajectory equation (2.2) about the orbital solution expressed by (2.3). Specifically, we write

(2.4a)\begin{equation} \boldsymbol{x}_p(t) = r_p(t)\left(\cos\theta_p(t),\sin\theta_p(t) \right), \end{equation}

where $r_p(t)$ and $\theta _p(t)$ are the time-varying radial and angular polar coordinates of the droplet's position, respectively. For a small perturbation from an orbital trajectory, we consider solutions of the form

(2.4b)\begin{equation} r_p(t) = r_0 + \epsilon\,r_1(t) \quad \mathrm{and}\quad \theta_p(t) = \omega t + \epsilon\,\theta_1(t), \end{equation}

where $r_0$ and $\omega$ satisfy the orbital equations (2.3), and $\epsilon \ll 1$ is a small parameter. We substitute (2.4) into (2.2), retain terms to $O(\epsilon )$, and then take the Laplace transform of the resultant linear equations. It follows that the perturbed trajectory's asymptotic complex growth rates $s$ satisfy $F(s) = 0$, where

(2.5)\begin{equation} F(s) = \mathscr{A}(s)\,\mathscr{D}(s) + \mathscr{B}(s)\,\mathscr{C}(s), \end{equation}

and the stability coefficients are defined (in a form equivalent to Oza et al. Reference Oza, Harris, Rosales and Bush2014a, Reference Oza, Rosales and Bush2018) as

(2.6a)\begin{gather} \mathscr{A}(s) = \kappa_0(s^2 - 2\omega^2) + \mu + s - 2\varOmega\omega + \mathcal{C}_0(s) + \mathcal{I}_1(s) - 2\,\mathcal{I}_0(0), \end{gather}

(2.6b)\begin{gather}\mathscr{B}(s) = 2\kappa_0\omega s + \varOmega s - \mu(\kappa_0\omega + \varOmega) - \mathcal{S}_0(s), \end{gather}

(2.6c)\begin{gather}\mathscr{C}(s) = 2\kappa_0\omega s + 2\omega + \varOmega s + \mu(\kappa_0\omega + \varOmega) - \mathcal{S}_0(s), \end{gather}

(2.6d)\begin{gather}\mathscr{D}(s) = \kappa_0 s^2 + s - \mu + \mathcal{C}_0(s) - \mathcal{I}_1(s). \end{gather}

Of particular interest in our investigation are the integrals (defined for $\text {Re}(s) > -\mu$ and any integer $m \geq 0$)

(2.7a)\begin{gather} \mathcal{I}_m(s) = \int_0^\infty \mathrm{J}_{2m}\left(2r_0\sin\left(\frac{\omega t}{2}\right)\right)\mathrm{e}^{-(\mu+s)t}\,\mathrm{d}t, \end{gather}

(2.7b)\begin{gather}\mathcal{C}_m(s) = \int_0^\infty \mathrm{J}_{2m}\left(2r_0\sin\left(\frac{\omega t}{2}\right)\right)\cos(\omega t)\,\mathrm{e}^{-(\mu+s)t}\,\mathrm{d}t, \end{gather}

(2.7c)\begin{gather}\mathcal{S}_m(s) = \int_0^\infty \mathrm{J}_{2m}\left(2r_0\sin\left(\frac{\omega t}{2}\right)\right)\sin(\omega t)\,\mathrm{e}^{-(\mu+s)t}\,\mathrm{d}t, \end{gather}

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 of Bessel functions of complex order. Specifically, we derive in Appendix B the closed-form expression

(2.8)\begin{equation} \mathcal{I}_m(s) = \frac{\rm \pi}{\omega}\,\mathrm{J}_{m+\mathrm{i}\eta}(r_0)\,\mathrm{J}_{m-\mathrm{i}\eta}(r_0)\operatorname{csch}({\rm \pi}\eta), \end{equation}

where $\eta = (\mu + s)/\omega$. Moreover, by representing $\cos (\omega t)$ and $\sin (\omega t)$ in terms of complex exponential functions, we deduce that

(2.9a,b)\begin{equation} \mathcal{C}_m(s) = \frac{1}{2}\left(\mathcal{I}_m(s+\mathrm{i}\omega) + \mathcal{I}_m(s-\mathrm{i}\omega)\right) \quad \mathrm{and}\quad \mathcal{S}_m(s) = \frac{1}{2\mathrm{i}}\left(\mathcal{I}_m(s-\mathrm{i} \omega) - \mathcal{I}_m(s+\mathrm{i}\omega)\right), \end{equation}

which we use to derive similar closed form formulae for $\mathcal {C}_m$ and $\mathcal {S}_m$. Using (2.8)–(2.9a,b), we derive in Appendix B simplified expressions for each of the stability integrals appearing in (2.6) in terms of products of Bessel functions of the first kind, $\mathrm {J}_\nu (r_0)$, and their derivatives $\mathrm {J}_\nu '(r_0)$, where the complex order $\nu$ takes values $\nu \in \{\pm \mathrm {i}(\mu + s)/\omega \}$. We then utilise 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, we parametrise the orbital dynamics entirely in terms of the radius $r_0$ (Oza Reference Oza2014), thereby effectively eliminating $\varOmega$ from the stability problem. This elimination process is achieved by first recasting the radial force balance (2.3a) as

(2.10)\begin{equation} \varOmega ={-}\kappa_0\omega - \frac{2}{r_0 \omega }\int_0^\infty \mathrm{J}_1\left(2r_0\sin\left(\frac{\omega t}{2}\right)\right)\sin\left(\frac{\omega t}{2}\right)\mathrm{e}^{-\mu t}\,\mathrm{d}t, \end{equation}

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

(2.11)\begin{equation} \int_0^\infty \mathrm{J}_1\left(2r_0\sin\left(\frac{\omega t}{2}\right)\right)\sin\left(\frac{\omega t}{2}\right)\mathrm{e}^{-\mu t}\,\mathrm{d}t ={-}\frac{1}{2}\, {\frac{\partial^{}{\mathcal{I}_0(0)}}{\partial{r_0}^{}}}. \end{equation}

By combining (2.10) and (2.11), we eliminate $\varOmega$ in the stability coefficients (2.6), yielding

(2.12a)\begin{gather} \mathscr{A}(s) = \kappa_0 s^2 - \frac{2}{r_0}\, {\frac{\partial^{}{\mathcal{I}_0(0)}}{\partial{r_0}^{}}} - 2I_0(0) + \mu + s + \mathcal{C}_0(s) + \mathcal{I}_1(s), \end{gather}

(2.12b)\begin{gather}\mathscr{B}(s) = \kappa_0\omega s - \frac{\mu-s}{r_0\omega}\, {\frac{\partial^{}{\mathcal{I}_0(0)}}{\partial{r_0}^{}}} - \mathcal{S}_0(s), \end{gather}

(2.12c)\begin{gather}\mathscr{C}(s) = \kappa_0\omega s + 2\omega + \frac{\mu+s}{r_0\omega}\, {\frac{\partial^{}{\mathcal{I}_0(0)}}{\partial{r_0}^{}}} - \mathcal{S}_0(s), \end{gather}

(2.12d)\begin{gather}\mathscr{D}(s) = \kappa_0 s^2 + s - \mu + \mathcal{C}_0(s) - \mathcal{I}_1(s). \end{gather}

Finally, we reduce the tangential force balance by integrating (2.3b) by parts, from which it follows that the orbital speed $U = r_0\omega$ satisfies (Oza et al. Reference Oza, Harris, Rosales and Bush2014a)

(2.13)\begin{equation} \mathcal{I}_0(0) = \frac{1}{\mu}\left(1 - \frac{r_0^2\omega^2}{2}\right). \end{equation}

For any given $r_0 > 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 $\text {Re}(s) > 0$. By denoting $s_*$ as the unstable root with largest real part, the instability is monotonic if $\text {Im}(s_*) = 0$, and oscillatory otherwise. The stability function $F$ has trivial eigenvalues at $0$ and $\pm \mathrm {i}\omega$, corresponding to rotational and translational invariance of the orbital motion, respectively (Oza et al. Reference Oza, Harris, Rosales and Bush2014a). It follows, therefore, that the non-trivial roots of the stability problem satisfy $G(s) = 0$, where

(2.14)\begin{equation} G(s) = \frac{\mathscr{A}(s)\,\mathscr{D}(s) + \mathscr{B}(s)\,\mathscr{C}(s)}{s(s^2+\omega^2)}. \end{equation}

We apply the method of Delves & Lyness (Reference Delves and Lyness1967) to find the roots of $G$ in the domain over which $G$ is analytic, i.e. $\text {Re}(s) > -\mu$. To ascertain whether a particular orbital state is stable or unstable, typically we utilise a rectangular integration contour spanning the domain $\text {Re}(s) \in [0,20]$ and $\text {Im}(s) \in [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\leq \varGamma \leq 0.99$). This approach differs from that of Oza et al. (Reference Oza, Rosales and Bush2018), who instead applied the methodology of Delves & Lyness (Reference Delves and Lyness1967) to $F$, integrating $F'/F$ over a deformed contour specifically chosen to avoid the trivial zeros at $s=0$ and $s=\pm \mathrm {i}\omega$. 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. (Reference Oza, Harris, Rosales and Bush2014a, Reference Oza, Rosales and Bush2018) in presenting the dependence of the orbital radius on the bath rotation rate for $\varGamma = 0.8$. In figure 2(b), we summarise the stability behaviour for all $\varGamma$. 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 as a wobbling orbit, in which the orbital centre remains approximately constant, but the radius of curvature exhibits small-amplitude oscillations with a frequency approximately twice that of the orbital frequency, as reported in the experiments of Harris & Bush (Reference Harris and Bush2014) and the numerical simulations of Oza et al. (Reference Oza, Wind-Willassen, Harris, Rosales and Bush2014b) (see figure 1d). Furthermore, the instability 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 the neighbouring wobbling instability tongues.

Figure 2. Orbital stability and the onset of quantisation for $\kappa _0 = 1.6$, where $\kappa _0 = m/(DT_W)$ is the inertial coefficient appearing in (2.2). (a) The dependence of orbital radius on bath rotation rate for $\varGamma = 0.8$, where $\varGamma = (\gamma - \gamma _W)/(\gamma _F - \gamma _W)$ is the dimensionless vibration parameter. (b) The delineation of orbital stability for any given radius and memory. The stability boundary is highlighted in white, and the yellow line at $\varGamma = 0.8$ corresponds to the curve in (a). In both plots, stable orbital states are indicated in blue, while oscillatory and monotonic 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 the order of the Faraday wavelength (Durey et al. Reference Durey, Turton and Bush2020b), which is subdominant for $r_0/\lambda _F \lesssim 4$, and truncates the instability tongues for $r_0/\lambda _F \gtrsim 4$.

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 arising near the tip of each instability tongue. In particular, we determine the main contributions to each stability integral arising along a stability boundary ($s = \mathrm {i} S$ with $S$ real) for large orbital radius ($r_0 \gg 1$), for which each integrand is highly oscillatory. Using (2.9a,b) to express $\mathcal {C}_m$ and $\mathcal {S}_m$ in terms of $\mathcal {I}_m$, we henceforth focus our attention on the study of $\mathcal {I}_m$. Furthermore, by recognising that $\mathcal {I}_m$ is a Laplace transform of a periodic function, we reduce the integral (2.7a) to

(2.15)\begin{equation} \mathcal{I}_m(\mathrm{i} S) = \frac{\mathcal{L}_m(\xi)}{\omega\left(1 - \mathrm{e}^{{-}2{\rm \pi}(\beta + \mathrm{i} \xi)}\right)}, \quad \mathrm{where}\ \mathcal{L}_m(\xi) = \int_0^{2{\rm \pi}}\mathrm{J}_{2m}\left(2r_0\sin\left(\frac{\theta}{2}\right)\right)\mathrm{e}^{-(\beta+\mathrm{i}\xi)\theta}\,\mathrm{d}\theta. \end{equation}

Here, $\beta = \mu /\omega$ is the inverse orbital memory, and $\xi = S/\omega$ is the scaled destabilisation frequency (see table 1). Based on experimental and numerical observations of monotonic and wobbling instabilities (Harris & Bush Reference Harris and Bush2014; Oza et al. Reference Oza, Wind-Willassen, Harris, Rosales and Bush2014b), we assume henceforth that the destabilisation frequency is comparable to the orbital frequency (i.e. $\xi = O(1)$). Notably, our analysis does not account for the instability associated with in-line speed oscillations arising at larger orbital radii (denoted by the pink curve in figure 2b), for which $S\sim U$, or $\xi \sim r_0$ (Durey et al. Reference Durey, Turton and Bush2020b). Finally, by writing $\beta = \mu r_0/U$ and noting that $\mu$ 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 $\beta$ is at most of $O(r_0)$ 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 $\beta$ and $r_0$, namely $\beta = O(\ln r_0)$, where $\mathrm {ln}$ denotes the natural logarithm.

Before proceeding with the asymptotic expansions, we provide a physical interpretation for the integral $\mathcal {L}_m(\xi )$. The argument $2r_0\sin ({\theta }/{2}) \geq 0$ is the length of the chord spanning two points lying an angle $\theta$ apart on a circle of radius $r_0$. This distance reflects the influence of 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 $\beta > 0$. Notably, $\mathrm {e}^{-{\rm \pi} \beta }$ 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, $\mathrm {e}^{-2{\rm \pi} \beta }$ determines the wave damping factor over a complete orbital period. Finally, the factor $\mathrm {e}^{-\mathrm {i} \xi \theta }$ accounts for oscillations in the perturbed droplet trajectory. For large $r_0$, the integrand of $\mathcal {L}_m(\xi )$ is generally highly oscillatory, with dominant contributions arising over 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 $\mathrm {J}_{2m}$ vanishes, i.e. at $\theta =0$ and $\theta =2{\rm \pi}$ (Bleistein & Handelsman Reference Bleistein and Handelsman1975). The internal points of stationary phase arise when the argument of the Bessel function is stationary, i.e. at $\theta = {\rm \pi}$. We now proceed to determine the magnitude of the contribution made by each critical point.

We first examine the contributions to $\mathcal {L}_m$ arising about $\theta = 0$ and $\theta = 2{\rm \pi}$, which we denote by $\mathcal {L}_{m,0}$ and $\mathcal {L}_{m, 2{\rm \pi} }$, respectively. We derive in Appendix C the leading-order contribution $\mathcal {L}_{m,0} = O(r_0^{-1})$, which is valid when $\xi = O(1)$ and $\beta$ is of maximum size $O(r_0)$; both of these conditions are met near the tip of each instability tongue. Using the structure of the integrand of $\mathcal {L}_m(\xi )$, we determine similarly that the leading-order contribution about $\theta =2{\rm \pi}$ satisfies $\mathcal {L}_{m,2{\rm \pi} }(\xi ) = \mathrm {e}^{-2{\rm \pi} (\beta + \mathrm {i} \xi )}\mathcal {L}_{m,0}(\xi )$; thus the relative size of $\mathcal {L}_{m,0}(\xi )$ and $\mathcal {L}_{m,2{\rm \pi} }(\xi )$ is controlled by the orbital damping factor $\mathrm {e}^{-2{\rm \pi} \beta }$. Finally, we use the method of stationary phase (see Appendix C) to determine that the interior point contribution to $\mathcal {L}_m$ at ${\rm \pi}$ has magnitude

(2.16)\begin{equation} \mathcal{L}_{m,{\rm \pi}} = O\left(\frac{\mathrm{e}^{-{\rm \pi}\beta}}{r_0}\right). \end{equation}

The relative weight of the integral contributions about $\theta = 0$, ${\rm \pi}$ and $2{\rm \pi}$ seemingly decreases consecutively by a factor $\mathrm {e}^{-{\rm \pi} \beta }$. However, the contribution about $\theta = {\rm \pi}$ becomes significant when $\mathrm {e}^{{\rm \pi} \beta }$ scales algebraically with $r_0$, corresponding to a strong influence of the waves generated diametrically opposite the droplet's current position, as is characteristic of high orbital memory (Fort et al. Reference Fort, Eddi, Moukhtar, Boudaoud and Couder2010). 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 $\mathrm {e}^{{\rm \pi} \beta } = O(r_0^2)$, which motivates the asymptotic scaling relationships utilised in the forthcoming analysis (§ 3). In fact, we observe that each instability tongue is bounded above in memory by either the oscillatory instability threshold satisfying $\mathrm {e}^{{\rm \pi} \beta } = O(r_0)$ (gold dashed line), or the in-line speed oscillation instability threshold (purple dashed line) arising for free walkers (Bacot et al. Reference Bacot, Perrard, Labousse, Couder and Fort2019; Hubert et al. Reference Hubert, Labousse, Perrard, Labousse, Vandewalle and Couder2019; Durey et al. Reference Durey, Turton and Bush2020b). Detailing the latter instability (denoted by the pink curve in figure 2b), which manifests as in-line speed oscillations along the circular orbit with amplitudes of the order of the Faraday wavelength, and arises for orbits so large as to be inaccessible within the laboratory, will be the subject of a future investigation.

Figure 3. Envelopes of the stability boundary, the scalings of which are deduced from our stationary phase analysis (see § 2.4). The blue curve corresponds to the white stability boundary in figure 2(b). The green and red lines correspond to the wobbling and monotonic instability envelopes. The gold and purple dashed lines denote the upper bounds on the existence of stable circular orbits, corresponding to the transition from wobbling to monotonic instabilities, and the onset of in-line speed oscillations (pink curve in figure 2b), respectively. The half-orbit wave damping factor $\mathrm {e}^{-{\rm \pi} \beta }$ scales as $(r_0 k_F)^{-2}$ for the green and red lines, and as $(r_0 k_F)^{-1}$ for the gold line.

3.1. Asymptotic expansion

We proceed by using the asymptotic expansions of the integrals in (2.12) to identify the imaginary roots of the stability function $G$ defined in (2.14) when $r_0 \gg 1$. Specifically, we seek the roots that minimise the critical memory of instability, which arise along the stability boundary in figure 2(b). We utilise asymptotic expansions for the closed-form expressions of the stability integrals (see Appendix B), which involve Bessel functions $\mathrm {J}_\nu (r_0)$ of the first kind with complex order $\nu$, and their derivatives with respect to argument $\mathrm {J}_\nu '(r_0)$ (see (B7)). As the complex order takes values $\nu \in \{\pm \mathrm {i} \beta, \pm \mathrm {i} (\beta \pm \mathrm {i} \xi )\}$ along the stability boundary, where $\beta = O(\ln (r_0))$ and $\xi = O(1)$ for $r_0 \gg 1$, the argument of each Bessel function is asymptotically large relative to its order. We may thus expand each of the stability coefficients in (2.12) when $r_0 \gg 1$, utilising the dominant balance $\mathrm {e}^{{\rm \pi} \beta } = O(r_0^2)$, with details presented in Appendix D. Likewise, we use the large-argument expansions of the Bessel functions, valid when $\beta = O(\ln (r_0)) \ll \sqrt {r_0}$, to deduce that

(3.1)\begin{equation} U^2 = 2\left(1 - \frac{\beta}{r_0}\right) + O\left(\frac{1}{r_0^3}\right), \end{equation}

which will be utilised throughout the following analysis.

As our aim is to identify the imaginary roots of the equation $G(s) = 0$, we proceed by substituting the asymptotic approximations expressed in (D2) into (2.14). Moreover, we determine the orbital speed $U$, scaled destabilisation frequency $\xi$, and inverse orbital memory $\beta$, by means of an asymptotic expansion in terms of the small parameter $r_0^{-1}$, namely

(3.2a–c)\begin{equation} U = U_0 + \frac{U_1}{r_0} + O\left(\frac{1}{r_0^2}\right), \quad \xi = \xi_0 + \frac{\xi_1}{r_0} + O\left(\frac{1}{r_0^2}\right), \quad \beta = \beta_0 + \frac{\beta_1}{r_0} + O\left(\frac{1}{r_0^2}\right). \end{equation}

In § 3.2, we present the solution to the leading-order problem, for which we demonstrate systematically that $\xi = 2 + O(r_0^{-2})$ at the tip of each wobbling instability tongue, corresponding to a wobbling instability with wobbling angular frequency $2\omega$, the so-called 2-wobble (Harris & Bush Reference Harris and Bush2014; Oza et al. Reference Oza, Wind-Willassen, Harris, Rosales and Bush2014b; 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 the asymptotic 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 $U_0$, $\xi _0$ and $\beta _0$. By using the asymptotic relationship $\mathscr {D} = \mathrm {i}\xi \mathscr {C} + O(r_0^{-3})$ deduced in (D2), we find that the stability condition (2.14) satisfies

(3.3)\begin{equation} \frac{\mathscr{C}(\mathrm{i}\xi \mathscr{A} + \mathscr{B})}{\mathrm{i}\xi(\xi^2 - 1)} = O\left(\frac{1}{r_0^4}\right), \end{equation}

whereupon substituting the leading-order expressions for $\mathscr {A}$, $\mathscr {B}$, $U$, $\xi$ and $\beta$ from (D2) and (3.2a–c) results in

(3.4)\begin{equation} \frac{\mathscr{C}}{U_0}\left[\frac{2\sin(2r_0)}{ \xi_0^2 - 1}\left(\operatorname{csch}({\rm \pi}(\beta_0 + \mathrm{i}\xi_0)) + \operatorname{csch}({\rm \pi}\beta_0)\right) + \frac{\kappa_0 U_0^3 + 1}{r_0^2}\right] = O\left(\frac{1}{r_0^3}\right). \end{equation}

As $\mathscr {C} = 2U_0/r_0 + O(r_0^{-2})$ is non-zero to leading order (see (D2)), the leading-order solution 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 $\xi _0$ is an integer, whose possible values will be the focus of the remainder of this subsection. We proceed to eliminate possible integer values for $\xi _0$ by looking for the solutions to the stability problem that occur at the highest possible value of $\beta _0$, as these solutions correspond to the instabilities arising at lowest memory for a given orbital radius. Our analysis will show that only two solutions are possible: (i) $\xi _0 = 0$, corresponding to a monotonic instability; and (ii) $\xi _0 = 2$, corresponding to a $2\omega$ instability.

To explore the possibility of $\xi _0$ being odd, we first consider the limit $\xi _0\to 1$ in (3.4). By applying L'Hôpital's rule, we find that the leading-order stability condition reduces to

(3.5)\begin{equation} 1+\kappa_0 U_0^3 + \mathrm{i}{\rm \pi} r_0^2\sin(2r_0)\coth({\rm \pi}\beta_0)\operatorname{csch}({\rm \pi}\beta_0) = 0. \end{equation}

As the real parts cannot be balanced (since $U_0 > 0$), there are no solutions to this equation. Similarly, if $\xi _0$ were odd and not equal to 1, the leading-order stability condition (3.4) would become

(3.6)\begin{equation} 1 + \kappa_0 U_0^3 = 0, \end{equation}

which is also impossible to satisfy. We thus conclude that $\xi _0$ cannot be odd, meaning that the destabilisation frequency must be an even multiple of the orbital frequency.

To explore the possible even values of $\xi _0$, we denote $\xi _0 = 2n$ (where $n$ is an integer) and use the approximation $\sinh (x) \approx \cosh (x) \approx \tfrac {1}{2}\mathrm {e}^x$ for $x \gg 1$; as such, the leading-order stability condition (3.4) reduces to

(3.7)\begin{equation} \mathrm{e}^{{\rm \pi}\beta_0} = \frac{8r_0^2\sin(2r_0)}{(1-4n^2)(1 + \kappa_0 U_0^3)}. \end{equation}

To be consistent with the assumed scaling of $\mathrm {e}^{{\rm \pi} \beta _0} = O(r_0^2)$, we require $\sin (2r_0) = O(1)$. As $\sin (2r_0)$ can be either positive or negative, the lowest memory (or largest $\beta _0$) condition 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 $U_0 = \sqrt {2}$ from (3.1) and (3.2a–c), we thus arrive at the monotonic stability boundary

(3.8)\begin{equation} \mu_{{mon}} = \frac{\sqrt{2}}{{\rm \pi} r_0}\ln\left(\frac{8r_0^2\sin(2r_0)}{1 + 2\sqrt{2}\kappa_0}\right), \end{equation}

which is valid when $\sin (2r_0) > 0$ and $\sin (2r_0) = O(1)$. For the case $n \neq 0$, we observe that $1-4n^2 < 0$; we thus deduce the requirement $\sin (2r_0) < 0$. The magnitude of the right-hand side of (3.7) is then minimised at $n=1$, which corresponds to $\xi = 2$, or $s = 2\mathrm {i}\omega$. 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 a wobbling instability is

(3.9)\begin{equation} \mu_{{wob}} = \frac{\sqrt{2}}{{\rm \pi} r_0}\ln\left(-\frac{8r_0^2\sin(2r_0)}{3(1 + 2\sqrt{2}\kappa_0)}\right), \end{equation}

which is valid when $\sin (2r_0) < 0$ and $\sin (2r_0) = O(1)$.

The asymptotic expressions (3.8)–(3.9) for the instability memory give rise to alternation between wobbling ($\sin (2r_0) < 0$) and monotonic ($\sin (2r_0) > 0$) instabilities with increasing orbital radius. In addition, $\mu _{{mon}}$ and $\mu _{{wob}}$ are maximised (corresponding to the tip of each instability tongue) when, to leading order for large orbital radius, $\sin (2r_0) = 1$ and $\sin (2r_0) = -1$, respectively. These extrema thus determine the wobbling and monotonic envelopes

(3.10a,b)\begin{equation} \mu_{{wob}}^{{env}} = \frac{\sqrt{2}}{{\rm \pi} r_0}\ln\left(\frac{8r_0^2}{3(2\sqrt{2}\kappa_0 + 1)}\right) \quad \mathrm{and}\quad \mu_{{mon}}^{{env}} = \frac{\sqrt{2}}{{\rm \pi} r_0}\ln\left(\frac{8r_0^2}{2\sqrt{2}\kappa_0 + 1}\right), \end{equation}

which are represented in figure 5(a) by the green and red dashed curves, respectively. Notably, increasing $\kappa _0$ increases the critical memory of wobbling and monotonic instabilities. Furthermore, since $\mu _{mon}^{env} - \mu _{wob}^{env} = \sqrt {2}\ln (3)/({\rm \pi} r_0) > 0$, we conclude that 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.

Figure 5. Comparison of the numerical solution and asymptotic approximations to the wobbling and monotonic stability boundaries. (a,b) The critical dimensionless vibration parameter at onset, $\varGamma = (\gamma -\gamma _W)/(\gamma _F - \gamma _W)$, and (c,d) the instability frequency relative to the orbital frequency, $\xi = S/\omega$, for $\kappa _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 dashed curves. (c,d) The critical wobbling frequency and its asymptotic counterpart (orange). All the asymptotic results presented in (a–d) are valid when $\sin (2k_Fr_0) = O(1)$, $k_Fr_0$ is large, and $\xi = 0$ (monotonic instabilities) or $\xi \approx 2$ (2-wobbles). Discontinuities in the grey curves (c,d) reflect transitions between monotonic ($\xi = 0$) and wobbling ($\xi \approx 2$) instabilities.

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 the most unstable orbital radii. However, the leading-order analysis did not provide insight into the behaviour of the instability frequency along each wobbling tongue, specifically the extent to which wobbling instabilities are approximated by a 2-wobble. So as to investigate this behaviour and so determine $U_1$, $\xi _1$ and $\beta _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 expansions in (3.2a–c). The expanded tangential force balance (D2) shows that $U_1 = -\beta _0/U_0$. By substituting into the stability condition $\mathrm {i} \xi \mathscr {A} + \mathscr {B} = 0$, using (D2) and (3.2a–c), retaining next-order terms and solving for the real and imaginary parts, we deduce that wobbling instabilities ($\xi _0 \neq 0$) have the first-order correction

(3.11a,b)\begin{equation} \xi_1 ={-}\frac{4\beta_0}{\rm \pi}\cot(2r_0)\quad\mathrm{and}\quad \beta_1 = \frac{3\sqrt{2}\kappa_0 \beta_0}{{\rm \pi}(2\sqrt{2}\kappa_0 + 1)} + \left(\frac{64\beta_0+27{\rm \pi}}{12{\rm \pi}^2} - \frac{\beta_0^2}{\rm \pi}\right)\cot(2r_0). \end{equation}

Similarly, for the monotonic instability ($\xi _0 = 0$), one may deduce that $\xi _1 = 0$ and

(3.12)\begin{equation} \beta_1 = \frac{3\sqrt{2}\kappa_0 \beta_0}{{\rm \pi}(2\sqrt{2}\kappa_0 + 1)} + \left(\frac{2{\rm \pi}\beta_0 + 1}{4{\rm \pi}} - \frac{\beta_0^2}{\rm \pi}\right)\cot(2r_0). \end{equation}

We recall that the asymptotic solution is valid when $\sin (2r_0) = O(1)$. Notably, $\xi _1$ vanishes for wobbling instabilities when $\cos (2r_0) = 0$, corresponding to $s = \mathrm {i} \omega (2 + O(r_0^{-2}))$ at the most unstable orbital radii (see § 3.2). Wobbling instabilities are thus driven by a resonance between the destabilisation frequency and the orbital frequency, with a larger critical memory necessary to destabilise orbits whose instability frequency deviates from twice the orbital frequency.

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 in figures 5(a,c) and first-order corrections in figures 5(b,d). Notably, our asymptotic results are valid when $r_0 \gg 1$, $\sin (2 r_0) = O(1)$ and either $\xi = 0$ (corresponding to monotonic instabilities) or $\xi \approx 2$ (corresponding to wobbling instabilities). Despite these restrictions, the leading-order memory captures the main features of the stability tongues, and the asymptotic instability frequency closely matches the numerical behaviour for $\xi = 0$ or $\xi \approx 2$. The success of our asymptotic results even for orbits of moderate radius is rooted in the choice of expansion parameter, namely $k_F r_0$ in dimensional variables, which is assumed to be large. We note that all orbits presented in figure 5 satisfy $k_F r_0 \geq 2{\rm \pi}$ (or $r_0/\lambda _F \geq 1$), which is evidently sufficiently large for our asymptotic results to yield reasonable agreement. Finally, we note that a similarly favourable agreement between our asymptotic and numerical results was obtained across a wide range of $\kappa _0$ values, including for those inaccessible in the laboratory (see § 4.1).

In summary, our asymptotic results explain the preponderance of 2-wobbles, with an exact 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 the monotonic instabilities arising at a memory lower than that of wobbling instabilities (figure 3).

4.1. Energy-like heuristics

The hydrodynamic pilot-wave system is a driven-dissipative system: energy is supplied by the system vibration and ultimately lost through viscous dissipation. Nevertheless, quasi-steady and periodic dynamical states arise in which the energy input precisely balances that lost through dissipation, and energy is exchanged primarily between the bouncing drop and the Faraday wave field excited by its impact (Moláček & Bush Reference Moláček and Bush2013a,Reference Moláček and Bushb). For inviscid gravity-capillary waves, kinetic energy is exchanged with gravitational potential and surface energies. For our investigation, it thus suffices to characterise the system 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 of wave intensity $E$, as coined by Hubert et al. (Reference Hubert, Perrard, Vandewalle and Labousse2022), which acts as a suitable proxy for the sum of the wave field gravitational potential and surface energies. In dimensional units, we thus define

(4.1)\begin{equation} E = \lim_{R\to\infty}\frac{1}{R}\left[\int_{|\boldsymbol{x}| \leq R} \frac{1}{2} \rho g h^2 \,\mathrm{d}\kern0.06em \boldsymbol{x} + \int_{|\boldsymbol{x}| \leq R} \sigma \left(\sqrt{1 + |\boldsymbol{\nabla} h|^2} - 1\right) \,\mathrm{d}\kern0.06em \boldsymbol{x}\right], \end{equation}

where $\sigma$ and $\rho$ 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 $\tfrac {1}{2}|\boldsymbol {\nabla } h|^2$ (see Appendix E). By applying the divergence theorem to the surface energy contribution and exploiting the fact that the wave field is monochromatic with wavenumber $k_F$, we demonstrate in § E.1 that the contributions from gravitational and surface energies are proportional to each other, and (4.1) reduces to

(4.2)\begin{equation} E = (\rho g + \sigma k_F^2) \lim_{R\to\infty}\frac{1}{2R}\int_{|\boldsymbol{x}| \leq R} h^2(\boldsymbol{x}, t) \,\mathrm{d}\kern0.06em \boldsymbol{x}, \end{equation}

an expression proportional to that obtained by Labousse et al. (Reference Labousse, Perrard, Couder and Fort2016b) and Hubert et al. (Reference Hubert, Perrard, Vandewalle and Labousse2022), who neglected the contribution of surface tension.

To analyse the wave intensity, we transform (4.2) to dimensionless variables by defining $\hat {h} = h/h_0$ and $\hat {E} = E/E_0$, where $h_0 = A T_W/T_F$ and $E_0 = h_0^2 k_F^{-1}(\rho g + \sigma k_F^2)$ are the characteristic wave height and intensity. By once again taking $k_F^{-1}$ and $T_W$ as the units of length and time (see § 2.1), we define the dimensionless wave intensity by

(4.3)\begin{equation} \hat{E} = \lim_{R\to\infty}\frac{1}{2R}\int_{|\boldsymbol{x}| \leq R} \hat{h}^2(\boldsymbol{x},t)\,\mathrm{d}\kern0.06em \boldsymbol{x}, \end{equation}

where

(4.4)\begin{equation} \hat{h}(\boldsymbol{x},t) = \int_{-\infty}^t \mathrm{J}_0(|\boldsymbol{x} - \boldsymbol{x}_p(s)|)\,\mathrm{e}^{-\mu(t-s)}\,\mathrm{d}s \end{equation}

is the dimensionless form of the pilot wave.

For circular orbital motion, the dimensionless wave intensity takes the remarkably simple form (see § E.2)

(4.5)\begin{equation} \hat{E} = \frac{1}{\mu}\,\mathcal{I}_0(0), \end{equation}

where $\mathcal {I}_0(0)$ (defined in (2.7a)) is the amplitude of the wave field beneath the droplet (Durey & Bush Reference Durey and Bush2021). Since the wave intensity and droplet gravitational potential energy are proportional to each other for orbital motion, either quantity serves equally well as a diagnostic 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

(4.6)\begin{equation} \hat{E} = \frac{1}{\mu^2}\left(1 - \frac{U^2}{2}\right), \end{equation}

where $U = r_0 \omega$ is the orbital speed. For fixed wave decay rate $\mu$, 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 from bouncing to rectilinear walking (Durey & Milewski Reference Durey and Milewski2017). Finally, combining the bound $\hat {E} > 0$ with (4.6) supplies the upper bound on the orbital speed noted in § 2.2, namely $U < \sqrt {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 the orbital radius in figure 6. As the wave decay rate $\mu$ 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 the wave intensity. (We note that the other wobbling instabilities appear as side bands to the first monotonic instability tongue; see figure 2(b).) This correlation is thus indicative of an underlying orbital energy principle, according to which the magnitude of the wave intensity prescribes the stability of the corresponding circular orbit (Durey Reference Durey2018). 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 radii that locally minimise $\varGamma$ along the stability boundary (the white curve in figure 2b). For each critical value of $\varGamma$ (or $\mu = 1- \varGamma$), we compare the corresponding critical orbital radius to the extrema of the wave field intensity (blue curves), confirming our observation that maxima and minima of $E$ correspond approximately to the tips of the monotonic and wobbling instability tongues, respectively, with the agreement improving for larger orbital radius.

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

Figure 7. The dependence of the critical orbital radius on the inertia-to-drag ratio, $\kappa _0 = m/DT_W$, 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 and wobbling instabilities to those that maximise and minimise the wave intensity ((4.5), blue) and the derivative of the mean wave force (orange, § 4.2), respectively. We also compare the critical radii with those predicted by 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 $\mathrm {J}_0(k_F r_0)$ (dotted), $\mathrm {J}_1(k_F r_0)$ (dashed) and $\mathrm {J}_2(k_F r_0)$ (dot-dashed), whose relation to orbital stability was noted by Labousse et al. (Reference Labousse, Oza, Perrard and Bush2016a). With increasing orbital radius, the range of each vertical axis narrows, reflecting the improved prediction of each heuristic.

To further assess the utility of different energy-based heuristics, we also compare the critical radii to the zeros of $\mathrm {J}_0(r_0)$, $\mathrm {J}_1(r_0)$ and $\mathrm {J}_2(r_0)$, as suggested by Labousse et al. (Reference Labousse, Oza, Perrard and Bush2016a). As presented in figure 7, the monotonic instabilities appear at radii slightly smaller than the zeros of $\mathrm {J}_1(r_0)$, while the wobbling instabilities typically align closely with the zeros of $\mathrm {J}_2(r_0)$. We note, however, that the critical radii exhibit a weak dependence on $\kappa _0$ that is not captured by the zeros of Bessel functions. On the other hand, the zeros of $\mathrm {J}_1(r_0)$ approach the radii satisfying $\sin (2r_0) = 1$ for larger orbital radius, consistent with our asymptotic analysis of the monotonic instability (§ 3.2). Likewise, the zeros of $\mathrm {J}_0(r_0)$ and $\mathrm {J}_2(r_0)$ both approach the radii satisfying $\sin (2r_0) = -1$, in agreement with our analysis of wobbling instabilities (§ 3.2). We thus rationalise the Labousse et al. (Reference Labousse, Oza, Perrard and Bush2016a) 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 $\kappa _0$, and satisfying the anticipated $O(r_0^{-2})$ convergence as the orbital radius is increased. One limitation of this approach, however, is the absence of a local minimum for the first wobbling tongue when using the first-order correction (3.10a,b) (owing to the parasitic influence of the $\cot (2r_0)$ term). Instead, we compare the critical radius in this case to that computed from the leading-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 are moderately successful, with the agreement improving for larger orbital radius. The zeros of Bessel functions generally give better agreement with the numerical results than the extrema of the wave intensity, or, equivalently, the droplet's gravitational potential energy. Notably, incorporating the droplet's kinetic energy within this latter heuristic does not affect the critical orbital radii significantly, with the resultant curves generally being indistinguishable from the blue curves in figure 7. Finally, these heuristic arguments are limited by their inability to predict the critical memory at the onset of instability and to capture the dependence of the critical radii on $\kappa _0$. Both quantities may be computed accurately using our asymptotic framework.

4.2. The mean wave field

Although the zeros of $\mathrm {J}_1(r_0)$ and $\mathrm {J}_2(r_0)$ provide satisfactory agreement with the numerical results for the tip of each instability tongue (Labousse et al. Reference Labousse, Oza, Perrard and Bush2016a), 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. (Reference Labousse, Oza, Perrard and Bush2016a), 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-axisymmetric component that serves to propel the droplet at a constant horizontal speed (Bush, Oza & Moláček Reference Bush, Oza and Moláček2014; Labousse & Perrard Reference Labousse and Perrard2014), and a static axisymmetric component (the mean wave field) that imparts either an inward or outward radial force to the droplet (Labousse et al. Reference Labousse, Perrard, Couder and Fort2014; Durey et al. Reference Durey, Milewski and Bush2018). 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. (Reference Perrard, Labousse, Miskin, Fort and Couder2014b) 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 $\bar {h}(r)$ accompanying a droplet executing a circular orbit of radius $r_0$ about the origin is given by (Perrard et al. Reference Perrard, Labousse, Miskin, Fort and Couder2014b; Tambasco & Bush Reference Tambasco and Bush2018)

(4.7)\begin{equation} \bar{h}(r) = \frac{1}{\mu}\,\mathrm{J}_0(r_0)\,\mathrm{J}_0(r). \end{equation}

We thus deduce that the dimensionless radial wave force

(4.8)\begin{equation} \mathcal{F}(r_0) =\left. -2\,\frac{\mathrm{d} \bar{h}}{\mathrm{d} r}\right|_{r = r_0} \end{equation}

applied by the mean wave field along the droplet's trajectory is

(4.9)\begin{equation} \mathcal{F}(r_0) = \frac{2}{\mu}\,\mathrm{J}_0(r_0)\,\mathrm{J}_1(r_0) ={-}\frac{2\cos(2r_0)}{\mu{\rm \pi} r_0} + O\left(\frac{1}{r_0^2}\right), \end{equation}

with derivative

(4.10)\begin{equation} \mathcal{F}'(r_0) = \frac{4\sin(2r_0)}{\mu{\rm \pi} r_0} + O\left(\frac{1}{r_0^2}\right). \end{equation}

As shown in § 3.1, wobbling instabilities occur when $\sin (2r_0) < 0$, whereas monotonic instabilities occur when $\sin (2r_0) > 0$. Equation (4.10) suggests that for sufficiently large $r_0$, the type of instability exhibited by increasing the memory at constant $r_0$ is thus related to the derivative of the mean wave force. Specifically, if $\mathcal {F}'(r_0) < 0$, then the corresponding circular orbit destabilises via a wobbling instability. Conversely, if $\mathcal {F}'(r_0) > 0$, then the circular orbit destabilises via a monotonic instability.

This correlation between the sign of $\mathcal {F}'(r_0)$ 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 $r_p(t)$ has the form $r_p(t) = r_0 + \epsilon \,r_1(t)$, where $0 < \epsilon \ll 1$ is a small parameter, and $r_1(t)$ denotes the perturbation to the orbital radius. It follows that the force exerted by the mean wave field may be approximated by $\mathcal {F}(r_p) \approx \mathcal {F}(r_0) + \epsilon \,r_1 \mathcal {F}'(r_0)$; thus the direction and magnitude of the perturbed wave force is prescribed by $\mathcal {F}'(r_0)$. When $\mathcal {F}'(r_0) > 0$, an outward radial perturbation ($r_1 > 0$) results in an increase in the outward force that drives the droplet away from equilibrium. Similarly, an inward perturbation ($r_1 < 0$) decreases the outward force. When $\mathcal {F}'(r_0) > 0$, we may thus regard the mean wave force as a repulsive spring force that induces monotonic changes in the orbital radius. Conversely, when $\mathcal {F}'(r_0) < 0$, the mean wave force behaves like an attractive spring that opposes any perturbations in the orbital radius from $r_0$ and so induces oscillations in the orbital radius. This physical picture is consistent with the observation that the type of instability changes when $\mathcal {F}'(r_0)$ changes sign (see figure 8). Moreover, the mean wave field's opposition to perturbations from $r_0$ hinders the onset of the corresponding wobbling instabilities, which is consistent with wobbling instabilities occurring at higher memory than their monotonic counterparts. Finally, the magnitude of the radial wave force increases as $\mu$ decreases, thereby increasing the sensitivity to perturbation at higher memory.

Figure 8. Schematic representation of orbital instability in terms of the mean wave force. We present the dependence of $\tfrac {1}{2}\mu \,\mathcal {F}(r_0) = \mathrm {J}_0(k_F r_0)\,\mathrm {J}_1(k_F r_0)$, which is proportional to the mean wave force, on the orbital radius for $\kappa _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 like an attractive and repulsive spring, respectively, giving rise to wobbling and monotonic instabilities. The sign of the mean wave force is denoted by arrows, and its slope is maximised and minimised at the yellow and purple squares, respectively.

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 and posit that these will correspond to the radii marked by the largest relative change in the perturbation force, as characterised by the coefficient $\mathcal {F}'(r_0)$. The coefficient $\mathcal {F}'(r_0)$ is maximised in magnitude at critical radii $r_c$ satisfying $\mathcal {F}''(r_c) = 0$. Specifically, maxima in $\mathcal {F}'(r_0)$ correspond to monotonic instabilities, and minima in $\mathcal {F}'(r_0)$ correspond to wobbling instabilities. In figure 7, we observe that the most unstable radii predicted by this heuristic ($\mathcal {F}''(r_c) = 0$) exhibit excellent quantitative agreement with the numerical solution 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 and monotonic instabilities in terms of the force exerted on the droplet by the mean wave field. Specifically, our investigation indicates that the critical radii $r_c$ at the onset of instability approximately satisfy $\mathcal {F}''(r_c) = 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 and quasi-periodic nature of the stability boundary (see figure 2b) is correlated with the quasi-monochromatic mean wave field. As a caveat, our rationale is valid only for circular 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 the parameter $\kappa _0$. Nevertheless, it performs remarkably well over the range of values of the inertia-to-drag parameter $\kappa _0$ accessible in the laboratory.