2.1 Wave dynamics

We employ the governing equations derived by Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015), who considered an incompressible viscous fluid in an infinite domain, with a small vortical boundary layer at the surface. Since a bouncing droplet emits radially symmetric waves, it is natural to write the spatial system in cylindrical polar coordinates $(r,\unicode[STIX]{x1D703},z)$ . Assuming that the waves and fluid velocity are small with a shallow wave slope, the velocity potential $\unicode[STIX]{x1D719}(r,\unicode[STIX]{x1D703},z,t)$ and wave perturbation $\unicode[STIX]{x1D702}(r,\unicode[STIX]{x1D703},t)$ at time $t>0$ satisfy the linearised system

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D719}+\unicode[STIX]{x1D719}_{zz},\quad z\leqslant 0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{t}=-g_{\unicode[STIX]{x1D6E4}}(t)\unicode[STIX]{x1D702}+\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D702}+2\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D719}-\frac{1}{\unicode[STIX]{x1D70C}}P(r,\unicode[STIX]{x1D703},t),\quad z=0, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{t}=\unicode[STIX]{x1D719}_{z}+2\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D702},\quad z=0, & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\rightarrow \mathbf{0}$ and $\unicode[STIX]{x1D702}\rightarrow 0$ in the far field. In the above, we have constant surface tension $\unicode[STIX]{x1D70E}$ , density $\unicode[STIX]{x1D70C}$ and kinematic viscosity $\unicode[STIX]{x1D708}$ , where $\unicode[STIX]{x1D6FB}_{H}^{2}$ is the horizontal Laplacian.

In the vibrating frame, the effective gravity is $g_{\unicode[STIX]{x1D6E4}}(t)\equiv g(1-\unicode[STIX]{x1D6E4}\cos (\unicode[STIX]{x1D714}_{0}t))$ , where $\unicode[STIX]{x1D714}_{0}/(2\unicode[STIX]{x03C0})$ is the vibration frequency and $\unicode[STIX]{x1D6E4}$ is the ratio of maximum vibration acceleration relative to gravity $g$ . As the droplet impacts the surface, the externally applied pressure on the bath is given by $P(r,\unicode[STIX]{x1D703},t)$ , which we prescribe in § 3. Coupled with droplet dynamics, Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015) solved this system numerically, but we make several simplifying assumptions to analyse a wide range of dynamics observed in an unbounded domain.

2.2 Droplet dynamics

By adapting the work of Moláček & Bush (Reference Moláček and Bush2013b ), we model the droplet dynamics as a rigid sphere of mass $m$ under ballistic motion centred at the point $(\boldsymbol{x},z)=(\boldsymbol{X}(t),Z(t))$ in Cartesian coordinates. This neglects droplet deformation, which is a reasonable assumption for the small droplets considered (Moláček & Bush Reference Moláček and Bush2013b ). During flight, the droplet experiences an aerodynamical force described by Stokes drag, namely $-\unicode[STIX]{x1D708}_{p}\boldsymbol{X}^{\prime }(t)$ and $-\unicode[STIX]{x1D708}_{p}Z^{\prime }(t)$ in the horizontal and vertical directions respectively. Here, $\unicode[STIX]{x1D708}_{p}=6\unicode[STIX]{x03C0}R_{0}\unicode[STIX]{x1D707}_{air}$ for droplet radius $R_{0}$ and air viscosity $\unicode[STIX]{x1D707}_{air}$ (Moláček & Bush Reference Moláček and Bush2013b ).

The droplet experiences a force of magnitude $f(t)\geqslant 0$ due to the lubrication air layer during impact ( $f=0$ during flight). As the droplet radius $R_{0}$ is much smaller than a typical wavelength $\unicode[STIX]{x1D706}$ ( $R_{0}/\unicode[STIX]{x1D706}\ll 1$ ), the forces experienced from the wave field interaction are localised to the point $\boldsymbol{x}=\boldsymbol{X}(t)$ on the fluid surface. These are modelled as an impulsive force $f(t)\hat{\boldsymbol{n}}(\boldsymbol{X}(t),t)$ and shear forces described by the skidding friction $-c\sqrt{\unicode[STIX]{x1D70C}R_{0}/\unicode[STIX]{x1D70E}}f(t)\boldsymbol{X}^{\prime }(t)$ . Here, $\hat{\boldsymbol{n}}=(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D702},1)/(1+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2})^{-1/2}$ is the unit normal away from the fluid surface and $c$ is the dimensionless skidding friction coefficient introduced by Moláček & Bush (Reference Moláček and Bush2013b ) (this is discussed in § 3.3). As the fluid model assumes a shallow wave slope $|\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|\ll 1$ , we approximate $\hat{\boldsymbol{n}}\sim (-\unicode[STIX]{x1D735}\unicode[STIX]{x1D702},1)$ .

For analogies to quantum mechanics, experiments are construed to subject the horizontal motion of the droplet to an external potential ${\mathcal{V}}(\boldsymbol{X}(t),t)$ , such as the dynamics in a rotating bath (Fort et al. Reference Fort, Eddi, Boudaoud, Moukhtar and Couder2010) and a horizontal harmonic potential well (Perrard et al. Reference Perrard, Labousse, Miskin, Fort and Couder2014b ). On combining these forces, conservation of momentum supplies

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle m\boldsymbol{X}^{\prime \prime }(t)+c\sqrt{\unicode[STIX]{x1D70C}R_{0}/\unicode[STIX]{x1D70E}}f(t)\boldsymbol{X}^{\prime }(t)+\unicode[STIX]{x1D708}_{p}\boldsymbol{X}^{\prime }(t)=-\unicode[STIX]{x1D735}{\mathcal{V}}(\boldsymbol{X}(t),t)-f(t)\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}(\boldsymbol{X}(t),t), & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle mZ^{\prime \prime }(t)+\unicode[STIX]{x1D708}_{p}Z^{\prime }(t)=-mg_{\unicode[STIX]{x1D6E4}}(t)+f(t), & \displaystyle\end{eqnarray}$$

written in the vibrating frame. In what follows, we take ${\mathcal{V}}\equiv 0$ for free walking dynamics, or ${\mathcal{V}}=(1/2)\unicode[STIX]{x1D705}|\boldsymbol{X}(t)|^{2}$ for dynamics in a harmonic potential well with an adjustable spring constant $\unicode[STIX]{x1D705}\geqslant 0$ . The analysis for other potentials follows akin to this work.

For simplicity, we assume that the bouncing droplet lies in the prevalent $(2,1)$ mode, which implies that both $Z(t)$ and $f(t)$ are periodic with subharmonic period $T=4\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{0}$ . By integrating (2.5) over an impact period and exploiting periodicity, $f(t)$ must satisfy

(2.6) $$\begin{eqnarray}mgT=\int _{\unicode[STIX]{x1D70F}}^{\unicode[STIX]{x1D70F}+T}f(t)\,\text{d}t,\quad \forall \unicode[STIX]{x1D70F}\geqslant 0,\end{eqnarray}$$

as derived by Moláček & Bush (Reference Moláček and Bush2013b ). Although we could carefully model the force $f(t)$ as a response to the impact dynamics (e.g. as a logarithmic spring (Moláček & Bush Reference Moláček and Bush2013a )), we use condition (2.6) to prescribe a reasonable choice of $f(t)$ in § 3.2.

Table 1. Fixed variables used in this model.

3.1 Basis function expansion

The assumption of infinite depth allows us to adapt the ideas of Benjamin & Ursell (Reference Benjamin and Ursell1954), whereby we decompose the wave perturbation $\unicode[STIX]{x1D702}$ and velocity potential $\unicode[STIX]{x1D719}$ in terms of orthogonal eigenfunctions of Laplace’s equation. Specifically, we use Bessel functions $\text{J}_{m}(\cdot )$ to pose the expansions

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}(r,\unicode[STIX]{x1D703},t)=\mathop{\sum }_{m=0}^{\infty }\int _{0}^{\infty }k(a_{m}(t;k)\unicode[STIX]{x1D6F7}_{m}(r,\unicode[STIX]{x1D703};k)+b_{m}(t;k)\unicode[STIX]{x1D6F9}_{m}(r,\unicode[STIX]{x1D703};k))\,\text{d}k, & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}(r,\unicode[STIX]{x1D703},z,t)=\mathop{\sum }_{m=0}^{\infty }\int _{0}^{\infty }k\text{e}^{kz}(c_{m}(t;k)\unicode[STIX]{x1D6F7}_{m}(r,\unicode[STIX]{x1D703};k)+d_{m}(t;k)\unicode[STIX]{x1D6F9}_{m}(r,\unicode[STIX]{x1D703};k))\,\text{d}k, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6F7}_{m}(r,\unicode[STIX]{x1D703};k)\equiv \text{J}_{m}(kr)\cos (m\unicode[STIX]{x1D703})$ and $\unicode[STIX]{x1D6F9}_{m}(r,\unicode[STIX]{x1D703};k)\equiv \text{J}_{m}(kr)\sin (m\unicode[STIX]{x1D703})$ satisfy

(3.9a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D6F7}_{m}(r,\unicode[STIX]{x1D703};k)=-k^{2}\unicode[STIX]{x1D6F7}_{m}(r,\unicode[STIX]{x1D703};k)\quad \text{and}\quad \unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D6F9}_{m}(r,\unicode[STIX]{x1D703};k)=-k^{2}\unicode[STIX]{x1D6F9}_{m}(r,\unicode[STIX]{x1D703};k).\end{eqnarray}$$

The coefficients $a_{m}$ , $b_{m}$ , $c_{m}$ and $d_{m}$ may be determined on substitution into (3.1)–(3.3). For clarity, we set $b_{0}\equiv d_{0}\equiv 0$ since $\unicode[STIX]{x1D6F9}_{0}(r,\unicode[STIX]{x1D703};k)\equiv 0$ for all $k\geqslant 0$ .

By choice of the orthogonal basis functions in (3.8), the continuity equation (2.1) is automatically satisfied. For horizontal droplet position $(r,\unicode[STIX]{x1D703})=(r_{d}(t),\unicode[STIX]{x1D703}_{d}(t))$ at time $t>0$ and small droplets, we model the droplet impacts at a point. Thus, we prescribe the pressure

(3.10) $$\begin{eqnarray}P(r,\unicode[STIX]{x1D703},t)=f(t)\frac{1}{r}\unicode[STIX]{x1D6FF}(r-r_{d}(t))\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{d}(t)),\end{eqnarray}$$

where $f(t)$ is the force applied by the droplet on the surface, which is zero during droplet flight. By exploiting the closure relation $\int _{0}^{\infty }kr\text{J}_{m}(kr)\text{J}_{m}(\unicode[STIX]{x1D709}r)\,\text{d}r=\unicode[STIX]{x1D6FF}(k-\unicode[STIX]{x1D709})$ and trigonometric orthogonality, we expand

(3.11) $$\begin{eqnarray}P(r,\unicode[STIX]{x1D703},t)=\mathop{\sum }_{m=0}^{\infty }\int _{0}^{\infty }k(p_{m}(t;k)\unicode[STIX]{x1D6F7}_{m}(r,\unicode[STIX]{x1D703};k)+q_{m}(t;k)\unicode[STIX]{x1D6F9}_{m}(r,\unicode[STIX]{x1D703};k))\,\text{d}k,\end{eqnarray}$$

where

(3.12a,b ) $$\begin{eqnarray}p_{m}(t;k)=\frac{1}{W_{m}}f(t)\unicode[STIX]{x1D6F7}_{m}(r_{d}(t),\unicode[STIX]{x1D703}_{d}(t);k)\quad \text{and}\quad q_{m}(t;k)=\frac{1}{W_{m}}f(t)\unicode[STIX]{x1D6F9}_{m}(r_{d}(t),\unicode[STIX]{x1D703}_{d}(t);k),\end{eqnarray}$$

with eigenfunction norms $W_{m}=\unicode[STIX]{x03C0}$ if $m>0$ and $W_{0}=2\unicode[STIX]{x03C0}$ . It should be noted that $q_{0}\equiv 0$ .

We substitute (3.7)–(3.11) into (3.2)–(3.3) and eliminate $c_{m}$ and $d_{m}$ in favour of $a_{m}$ and $b_{m}$ . By orthogonality, we obtain a system of inhomogeneous damped Mathieu equations,

(3.13a,b ) $$\begin{eqnarray}{\mathcal{L}}_{k}a_{m}(t;k)=-kMGp_{m}(t;k)\quad \text{and}\quad {\mathcal{L}}_{k}b_{m}(t;k)=-kMGq_{m}(t;k),\end{eqnarray}$$

where

(3.14) $$\begin{eqnarray}{\mathcal{L}}_{k}\equiv \unicode[STIX]{x2202}_{tt}+2\unicode[STIX]{x1D6FE}(k)\unicode[STIX]{x2202}_{t}+(\unicode[STIX]{x1D6FE}^{2}(k)+\unicode[STIX]{x1D714}^{2}(k)-\unicode[STIX]{x1D6E4}\unicode[STIX]{x1D714}_{g}^{2}(k)\cos (4\unicode[STIX]{x03C0}t+\unicode[STIX]{x1D6FD}))\end{eqnarray}$$

is the wavenumber-dependent homogeneous damped Mathieu differential operator and

(3.15a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}(k)\equiv 2\unicode[STIX]{x1D716}k^{2},\quad \unicode[STIX]{x1D714}^{2}(k)\equiv Gk+Bk^{3},\quad \unicode[STIX]{x1D714}_{g}^{2}(k)\equiv Gk.\end{eqnarray}$$

The $\unicode[STIX]{x1D6FE}^{2}(k)$ term gives the additional damping from the vortical boundary layer; this is not present in the work of Eddi et al. (Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011), who instead replace $\unicode[STIX]{x1D716}$ with a phenomenological damping term to match the Faraday threshold.

3.2 Instantaneous impacts

To close the system, it remains to prescribe $f(t)$ so that (3.5) is satisfied. To alleviate the difficulty of analysing inhomogeneous Mathieu equations, we simply model the impacts as instantaneous. This gives homogeneous equations for the waves and droplet during flight, and jump conditions at impact. This is the opposite limit to that of Oza et al. (Reference Oza, Rosales and Bush2013), where the droplet glides across the surface of the bath with constant forcing $f(t)$ . Assuming that the first impact is at $t=0$ , we define

(3.16) $$\begin{eqnarray}f(t)=\mathop{\sum }_{n=0}^{\infty }\unicode[STIX]{x1D6FF}(t-t_{n}),\end{eqnarray}$$

where $t_{n}\equiv n$ for all non-negative integers $n$ and $\unicode[STIX]{x1D6FF}(\cdot )$ is the Dirac delta function. This satisfies periodicity and the integral condition (3.5). We now exploit properties of $\unicode[STIX]{x1D6FF}(\cdot )$ to replace the inhomogeneities in (3.13) with jump conditions.

For physical consistency, we assume that the droplet position $\boldsymbol{X}(t)$ and wave amplitudes ( $a_{m}(t;k)$ and $b_{m}(t;k)$ ) are continuous across all impacts. We denote jumps $[Q(t_{n})]_{-}^{+}\equiv Q(t_{n}^{+})-Q(t_{n}^{-})$ for any function of time $Q(t)$ , where $t_{n}^{+}$ and $t_{n}^{-}$ are the right and left limits of $t_{n}$ respectively. We now integrate the governing equations (3.13) over $t\in [t_{n}^{-},t_{n}^{+}]$ . For all $n\geqslant 0$ , the first equation gives

(3.17) $$\begin{eqnarray}\int _{t_{n}^{-}}^{t_{n}^{+}}{\mathcal{L}}_{k}a_{m}(t;k)\,\text{d}t=\frac{-k}{W_{m}}MG\int _{t_{n}^{-}}^{t_{n}^{+}}\unicode[STIX]{x1D6FF}(t-t_{n})\unicode[STIX]{x1D6F7}_{m}(r_{d}(t),\unicode[STIX]{x1D703}_{d}(t);k)\,\text{d}t.\end{eqnarray}$$

As the droplet position is assumed to be continuous, the sifting property of $\unicode[STIX]{x1D6FF}(\cdot )$ can be applied to the right-hand side. We use the continuity of $a_{m}(t;k)$ to simplify the left-hand side. Hence, (3.13) supply jump conditions

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle [a_{m}^{\prime }(t_{n};k)]_{-}^{+}=-P_{m}(k)\unicode[STIX]{x1D6F7}_{m}(r_{d}(t_{n}),\unicode[STIX]{x1D703}_{d}(t_{n});k), & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle [b_{m}^{\prime }(t_{n};k)]_{-}^{+}=-P_{m}(k)\unicode[STIX]{x1D6F9}_{m}(r_{d}(t_{n}),\unicode[STIX]{x1D703}_{d}(t_{n});k), & \displaystyle\end{eqnarray}$$

where a prime denotes the partial derivative with respect to $t$ and $P_{m}(k)=kMG/W_{m}$ .

For the droplet dynamics (3.4), we first consider a single impact at $t=t_{\star }$ . Hence,

(3.20) $$\begin{eqnarray}\boldsymbol{X}^{\prime \prime }(t)+c\sqrt{R/B}G\unicode[STIX]{x1D6FF}(t-t_{\star })\boldsymbol{X}^{\prime }(t)+\tilde{\unicode[STIX]{x1D708}}_{p}\boldsymbol{X}^{\prime }(t)+\tilde{\unicode[STIX]{x1D705}}\boldsymbol{X}(t)=-G\unicode[STIX]{x1D6FF}(t-t_{\star })\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}(\boldsymbol{X}(t),t).\end{eqnarray}$$

For $c=0$ , we may proceed as above to find the jump in $\boldsymbol{X}^{\prime }(t)$ at $t=t_{\star }$ . The case $c>0$ is more delicate; the sifting property cannot be applied as $\boldsymbol{X}^{\prime }(t)$ is discontinuous at $t=t_{\star }$ . Following the method of Catllá et al. (Reference Catllá, Schaeffer, Witelski, Monson and Lin2008), we replace $\unicode[STIX]{x1D6FF}(t-t_{\star })$ with

(3.21) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}}(t-t_{\star })\equiv \unicode[STIX]{x1D700}^{-1}\unicode[STIX]{x1D711}((t-t_{\star })\unicode[STIX]{x1D700}^{-1}),\quad \forall \unicode[STIX]{x1D700}>0,\end{eqnarray}$$

where $\unicode[STIX]{x1D711}(\unicode[STIX]{x1D70F})\geqslant 0$ for all $\unicode[STIX]{x1D70F}\in \mathbb{R}$ and $\int _{-\infty }^{\infty }\unicode[STIX]{x1D711}(\unicode[STIX]{x1D70F})\,\text{d}\unicode[STIX]{x1D70F}=1$ . Hence, the functions $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}}(\cdot )\rightarrow \unicode[STIX]{x1D6FF}(\cdot )$ pointwise as $\unicode[STIX]{x1D700}\rightarrow 0$ , except at the jump discontinuity in $\unicode[STIX]{x1D6FF}(\cdot )$ . This is an appealing formulation, as, physically, no impact is actually instantaneous – it just occurs over a much faster time scale than the dynamics of the rest of the system.

The idea is to find a solution to (3.20) when $t>t_{\star }$ with $\unicode[STIX]{x1D6FF}(\cdot )$ replaced by $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}}(\cdot )$ , and then consider the limit as $\unicode[STIX]{x1D700}\rightarrow 0$ . This determines $\boldsymbol{X}^{\prime }(t_{\star }^{+})$ . However, for $t<t_{\star }$ , (3.20) can be solved directly without approximation (as the $\unicode[STIX]{x1D6FF}(\cdot )$ terms vanish), which supplies $\boldsymbol{X}^{\prime }(t_{\star }^{-})$ . Hence, the jump $[\boldsymbol{X}^{\prime }(t_{\star })]_{-}^{+}$ is obtained, which is independent of $\unicode[STIX]{x1D711}(\cdot )$ . By generalising the proof of Catllá et al. (Reference Catllá, Schaeffer, Witelski, Monson and Lin2008) (as outlined in appendix A) and extending to periodic impacts, we obtain the jump condition shown in (3.27) below.

3.3 Model summary

We have the dimensionless system

(3.22) $$\begin{eqnarray}\displaystyle \displaystyle 0 & = & \displaystyle {\mathcal{L}}_{k}a_{m}(t;k),\quad \forall t\neq t_{n},\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle \displaystyle 0 & = & \displaystyle {\mathcal{L}}_{k}b_{m}(t;k),\quad \forall t\neq t_{n},\end{eqnarray}$$

(3.24) $$\begin{eqnarray}\displaystyle \displaystyle \mathbf{0} & = & \displaystyle \boldsymbol{X}^{\prime \prime }(t)+\tilde{\unicode[STIX]{x1D708}}_{p}\boldsymbol{X}^{\prime }(t)+\tilde{\unicode[STIX]{x1D705}}\boldsymbol{X}(t),\quad \forall t\neq t_{n},\end{eqnarray}$$

(3.25) $$\begin{eqnarray}\displaystyle \displaystyle [a_{m}^{\prime }(t_{n};k)]_{-}^{+} & = & \displaystyle -P_{m}(k)\unicode[STIX]{x1D6F7}_{m}(\boldsymbol{X}(t_{n});k),\end{eqnarray}$$

(3.26) $$\begin{eqnarray}\displaystyle \displaystyle [b_{m}^{\prime }(t_{n};k)]_{-}^{+} & = & \displaystyle -P_{m}(k)\unicode[STIX]{x1D6F9}_{m}(\boldsymbol{X}(t_{n});k),\end{eqnarray}$$

(3.27) $$\begin{eqnarray}\displaystyle \displaystyle [\boldsymbol{X}^{\prime }(t_{n})]_{-}^{+} & = & \displaystyle -F(c)\left(\frac{1}{c}\sqrt{\frac{B}{R}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}(\boldsymbol{X}(t_{n}),t_{n})+\boldsymbol{X}^{\prime }(t_{n}^{-})\right),\end{eqnarray}$$

where $P_{m}(k)=kMG/W_{m}$ and $F(c)=1-\exp (-cG\sqrt{R/B})$ . One should note the abuse of notation $\unicode[STIX]{x1D6F7}_{m}(\boldsymbol{X}(t_{n});k)\equiv \unicode[STIX]{x1D6F7}_{m}(r_{d}(t_{n}),\unicode[STIX]{x1D703}_{d}(t_{n});k)$ (and similarly for $\unicode[STIX]{x1D6F9}_{m}$ ). The system of jump conditions is self-consistent, with both $\unicode[STIX]{x1D702}$ and $\boldsymbol{X}$ continuous across impacts.

This model has two undefined parameters: the skidding friction $c$ and the phase shift $\unicode[STIX]{x1D6FD}$ . From the theoretical calculations of Moláček & Bush (Reference Moláček and Bush2013b ), $c\approx 0.3$ , but the authors consider $c\in [0.17,0.33]$ for simulations. Oza et al. (Reference Oza, Rosales and Bush2013) use the lower bound $c=0.17$ in the stroboscopic approximation model, which has the opposite impact duration limit to our model. Hence, it is natural to choose the upper bound $c=0.33$ , which is fixed throughout the paper. The phase shift $\unicode[STIX]{x1D6FD}$ between bath vibrations and droplet impacts for periodic states arises naturally in models where the droplet vertical dynamics are explicitly modelled (Moláček & Bush Reference Moláček and Bush2013b ; Milewski et al. Reference Milewski, Galeano-Rios, Nachbin and Bush2015). As we restrict the droplet to periodic impacts, we must choose $\unicode[STIX]{x1D6FD}$ . We focus on the prevalent $(2,1)^{2}$ walking mode, where $\unicode[STIX]{x1D6FD}$ has a typical value of $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x03C0}$ (Milewski et al. Reference Milewski, Galeano-Rios, Nachbin and Bush2015), which is fixed henceforth. For later works, $c$ and $\unicode[STIX]{x1D6FD}$ may be tuned when comparing with experimental data.

3.4 Faraday instability

Following Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015) to determine the subharmonic Faraday instability, we look for subharmonic solutions to ${\mathcal{L}}_{k}a(t;k)=0$ of the form $a(t;k)=A\cos (2\unicode[STIX]{x03C0}t)+B\sin (2\unicode[STIX]{x03C0}t)$ . After substitution, higher-order frequencies are neglected, which is equivalent to truncating the Hill matrix. For a non-trivial system, we require

(3.28) $$\begin{eqnarray}(w^{2}(k)+\unicode[STIX]{x1D6FE}^{2}(k)-4\unicode[STIX]{x03C0}^{2})^{2}+16\unicode[STIX]{x03C0}^{2}\unicode[STIX]{x1D6FE}^{2}(k)-{\textstyle \frac{1}{4}}w_{g}^{4}(k)\unicode[STIX]{x1D6E4}^{2}=0.\end{eqnarray}$$

The function $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}(k)$ is globally minimised at $k=k_{F}$ , where for $\unicode[STIX]{x1D6E4}<\unicode[STIX]{x1D6E4}(k_{F})\equiv \unicode[STIX]{x1D6E4}_{F}$ the subharmonic solutions are stable for any $k>0$ . By our dimensionless scaling, $\unicode[STIX]{x1D706}_{F}\approx 1$ , where $k_{F}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}_{F}$ . Although it is known that $\unicode[STIX]{x1D6E4}_{F}$ obtained in this model is not accurate when compared with experiments (Milewski et al. Reference Milewski, Galeano-Rios, Nachbin and Bush2015), it is usual to use $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}$ as a controlling parameter (Eddi et al. Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011). In fact, we show in § 5.3 that the wave field temporal decay rate may be written as a function of $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}$ , which justifies our approach for a theoretical investigation.

3.5 Discrete-time model

As (3.22)–(3.27) form a homogeneous system with jump conditions, it is natural to reformulate as a discrete-time system. We denote $\boldsymbol{a}_{n}(k)\equiv (a_{0}(t_{n};k),a_{1}(t_{n};k),\ldots )^{\text{T}}$ and $\boldsymbol{a}_{n}^{\prime }(k)\equiv (a_{0}^{\prime }(t_{n}^{+};k),a_{1}^{\prime }(t_{n}^{+};k),\ldots )^{\text{T}}$ , and similarly for $\boldsymbol{b}_{n}(k)$ and $\boldsymbol{b}_{n}^{\prime }(k)$ . We also write the eigenfunctions as vectors $\unicode[STIX]{x1D731}(\cdot ;k)=(\unicode[STIX]{x1D6F7}_{0}(\cdot ;k),\unicode[STIX]{x1D6F7}_{1}(\cdot ;k),\ldots )^{\text{T}}$ , and similarly for $\unicode[STIX]{x1D733}(\cdot ;k)$ , with $\unicode[STIX]{x1D64B}(k)$ a diagonal matrix with elements $P_{m}(k)$ . By periodicity of the Mathieu operator ${\mathcal{L}}_{k}$ , we numerically construct the principal fundamental matrix $\unicode[STIX]{x1D648}_{k}(\unicode[STIX]{x1D6E4})$ over the interval $(0,1)$ . As ${\mathcal{L}}_{k}$ is independent of the Bessel order, $\unicode[STIX]{x1D648}_{k}$ has block diagonal form

(3.29) $$\begin{eqnarray}\unicode[STIX]{x1D648}_{k}(\unicode[STIX]{x1D6E4})=\left[\begin{array}{@{}cc@{}}m_{11}(\unicode[STIX]{x1D6E4};k)\unicode[STIX]{x1D644} & m_{12}(\unicode[STIX]{x1D6E4};k)\unicode[STIX]{x1D644}\\ m_{21}(\unicode[STIX]{x1D6E4};k)\unicode[STIX]{x1D644} & m_{22}(\unicode[STIX]{x1D6E4};k)\unicode[STIX]{x1D644}\end{array}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the infinite identity matrix. The principal fundamental matrix $\unicode[STIX]{x1D641}(\tilde{\unicode[STIX]{x1D705}})$ corresponding to the droplet dynamics is constructed by analytically solving (3.24).

We reformulate (3.22)–(3.27) as an efficient one-step map with stages for droplet position $\boldsymbol{X}_{n}\equiv \boldsymbol{X}(t_{n})$ and velocity $\boldsymbol{X}_{n}^{\prime }\equiv \boldsymbol{X}^{\prime }(t_{n}^{+})$ .

1. (i) Use $\boldsymbol{X}_{n}$ and $\boldsymbol{X}_{n}^{\prime }$ to compute $\boldsymbol{X}_{n+1}$ and $\tilde{\boldsymbol{X}}^{\prime }\equiv \boldsymbol{X}^{\prime }(t_{n+1}^{-})$ :

   (3.30) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}\boldsymbol{X}_{n+1}\\ \tilde{\boldsymbol{X}}^{\prime }\end{array}\right]=\unicode[STIX]{x1D641}(\tilde{\unicode[STIX]{x1D705}})\left[\begin{array}{@{}c@{}}\boldsymbol{X}_{n}\\ \boldsymbol{X}_{n}^{\prime }\end{array}\right].\end{eqnarray}$$

2. (ii) Update wave amplitudes including jump conditions $\forall k>0$ (similarly for $\boldsymbol{b}_{n}(k)$ ):

   (3.31) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}\boldsymbol{a}_{n+1}(k)\\ \boldsymbol{a}_{n+1}^{\prime }(k)\end{array}\right]=\unicode[STIX]{x1D648}_{k}(\unicode[STIX]{x1D6E4})\left[\begin{array}{@{}c@{}}\boldsymbol{a}_{n}(k)\\ \boldsymbol{a}_{n}^{\prime }(k)\end{array}\right]-\left[\begin{array}{@{}c@{}}\mathbf{0}\\ \unicode[STIX]{x1D64B}(k)\unicode[STIX]{x1D731}(\boldsymbol{X}_{n+1};k)\end{array}\right].\end{eqnarray}$$

3. (iii) Apply droplet jump conditions:

   (3.32) $$\begin{eqnarray}\boldsymbol{X}_{n+1}^{\prime }=(1-F(c))\tilde{\boldsymbol{X}}^{\prime }-\frac{F(c)}{c}\sqrt{\frac{B}{R}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}(\boldsymbol{X}_{n+1},t_{n+1}).\end{eqnarray}$$

3.6 Numerical implementation

The remainder of this work involves simulating and finding fixed points of the system (3.30)–(3.32), both of which require suitable truncations in the wavenumber $k>0$ and Bessel mode $m\in \mathbb{N}$ . For any impact time $t_{n}$ , the wavenumbers are generally peaked around $k=k_{F}$ , with the peak becoming narrower as $\unicode[STIX]{x1D6E4}\rightarrow \unicode[STIX]{x1D6E4}_{F}^{-}$ , which can be analysed from Floquet analysis of the damped Mathieu equation. This determines the refinement and truncation in $k$ , which is successively reduced until the change in the numerical solution of the fixed points becomes negligible. Away from the Faraday threshold, a reasonably coarse mesh is sufficient, with $\unicode[STIX]{x1D6FF}k\sim 0.2$ and $k\in [k_{F}/2,3k_{F}/2]$ , where $k_{F}$ is a mesh point. Integrals over $k$ (e.g. for finding $\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}$ ) are evaluated using the trapezium rule, which is well suited to capturing the peaked behaviour in $k$ .

Truncation of the Bessel modes follows from the asymptotic behaviour of Bessel functions, namely $\text{J}_{m}(z)\sim (1/m!)(z/2)^{m}$ for $0<z\ll \sqrt{m+1}$ . Hence, for orbital solutions or simulations with a central force, the maximum radial extent of the droplet can be estimated, which provides a good guide for the truncation $m\in \{0,1,\ldots ,m^{\star }\}$ . For walking dynamics, the Floquet exponents provide an estimate of how the temporal decay affects the number of past impacts that contribute to the current wave field (this is the system ‘memory’, as discussed in § 5.3), where we typically have $m^{\star }=15$ for $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}\approx 0.81$ but $m^{\star }\approx 50$ for $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}\approx 0.95$ . The accuracy of this truncation can be easily checked a posteriori, with $m^{\star }$ increased until there is a negligible change in the numerical solution.

For simulation efficiency, Bessel functions are only evaluated once per impact period, with derivatives calculated using the identity $\text{J}_{m}^{\prime }(z)=(1/2)(\text{J}_{m-1}(z)-\text{J}_{m+1}(z))$ . We typically simulate 1000 impacts per second on a standard desktop machine using MATLAB.

4.1 Stability analysis

To determine the walking threshold $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{W}$ , we perform linear stability analysis of the periodic bouncing system. The aim is to construct a one-step matrix map $\unicode[STIX]{x1D64F}(\unicode[STIX]{x1D6E4})$ for the perturbed system, where stability is determined by the spectral radius $\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D64F})$ . We denote the steady state by $\boldsymbol{X}=\hat{\boldsymbol{X}}$ , $a_{m}=\hat{a}_{m}$ and $b_{m}=\hat{b}_{m}$ , where, for bouncing at the origin, $\hat{\boldsymbol{X}}\equiv \mathbf{0}$ , $\hat{a}_{m}\equiv \hat{b}_{m}\equiv 0$ for all $m\geqslant 1$ , and $\hat{b}_{0}\equiv 0$ . We then consider small perturbations

(4.5a-c ) $$\begin{eqnarray}\boldsymbol{X}(t)=\hat{\boldsymbol{X}}(t)+\unicode[STIX]{x1D750}(t),\quad a_{m}(t;k)=\hat{a}_{m}(t;k)+\tilde{a}_{m}(t;k),\quad b_{m}(t;k)=\hat{b}_{m}(t;k)+\tilde{b}_{m}(t;k),\end{eqnarray}$$

where we assume that $|\tilde{a}_{m}/\hat{a}_{0}|\sim |\tilde{b}_{m}/\hat{a}_{0}|\sim ||\unicode[STIX]{x1D750}||\ll 1$ . By including an explicit perturbation to the wave field, we consider a more general perturbation than Oza et al. (Reference Oza, Rosales and Bush2013).

The system is linearised via the jump conditions (3.25)–(3.27), giving

(4.6) $$\begin{eqnarray}\displaystyle [\tilde{a}_{m}^{\prime }(t_{n};k)]_{-}^{+} & = & \displaystyle -P_{m}(k)\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}_{m}(\mathbf{0};k)^{\text{T}}\unicode[STIX]{x1D750}(t_{n}),\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle [\tilde{b}_{m}^{\prime }(t_{n};k)]_{-}^{+} & = & \displaystyle -P_{m}(k)\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}_{m}(\mathbf{0};k)^{\text{T}}\unicode[STIX]{x1D750}(t_{n}),\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle [\unicode[STIX]{x1D750}^{\prime }(t_{n})]_{-}^{+} & = & \displaystyle -F(c)\left(\unicode[STIX]{x1D750}^{\prime }(t_{n}^{-})+\frac{1}{c}\sqrt{\frac{B}{R}}(\unicode[STIX]{x1D643}(\unicode[STIX]{x1D6E4})\unicode[STIX]{x1D750}(t_{n})+\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D702}}(\mathbf{0},t_{n}))\right).\end{eqnarray}$$

The Hessian matrix of the steady-state wave field at droplet impact is $\unicode[STIX]{x1D643}(\unicode[STIX]{x1D6E4})$ and

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D702}}(\mathbf{0},t_{n})=\mathop{\sum }_{m=0}^{\infty }\int _{0}^{\infty }k(\tilde{a}_{m}(t_{n};k)\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}_{m}(\mathbf{0};k)+\tilde{b}_{m}(t_{n};k)\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}_{m}(\mathbf{0};k))\,\text{d}k.\end{eqnarray}$$

Equations (4.6)–(4.8) are simplified considerably by noting that

(4.10a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}_{1}(\mathbf{0};k)=(k/2,0)^{\text{T}},\quad \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}_{1}(\mathbf{0};k)=(0,k/2)^{\text{T}},\quad \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}_{m}(\mathbf{0};k)\equiv \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}_{m}(\mathbf{0};k)\equiv \mathbf{0},\end{eqnarray}$$

for $m\neq 1$ . Hence, for $m\neq 1$ , the jump conditions for $\tilde{a}_{m}$ , $\tilde{b}_{m}$ and $\unicode[STIX]{x1D750}$ are all independent of each other (to $O(||\unicode[STIX]{x1D750}||)$ ). Therefore, the perturbations $\tilde{a}_{m}$ , $\tilde{b}_{m}$ ( $m\neq 1$ ) each decouple from the system. As they are unexcited solutions to the damped Mathieu equation with $\unicode[STIX]{x1D6E4}<\unicode[STIX]{x1D6E4}_{F}$ , they are stable and hence are neglected from the stability analysis.

Using the linearised jump conditions with matrices $\unicode[STIX]{x1D648}_{k}(\unicode[STIX]{x1D6E4})$ and $\unicode[STIX]{x1D641}(0)$ , we construct a discrete-time linear map for all variables, given by the matrix $\unicode[STIX]{x1D64F}(\unicode[STIX]{x1D6E4})$ . The system is neutrally stable if the spectral radius $\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D64F})=1$ , and unstable if $\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D64F})>1$ . The walking threshold $\unicode[STIX]{x1D6E4}_{W}$ is the largest $\unicode[STIX]{x1D6E4}$ such that $\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D64F})=1$ . By the translational invariance of the system, there always exists an eigenvalue $\unicode[STIX]{x1D707}$ of $\unicode[STIX]{x1D64F}$ such that $\unicode[STIX]{x1D707}=1$ , which prevents us from obtaining asymptotically stable solutions $(\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D64F})<1)$ in the absence of a central force.

5.1 Stability analysis

The stability analysis of a walker is similar to that of a bouncer and therefore we do not give the details. The main difference is that we use Graf’s addition theorem (Abramowitz & Stegun Reference Abramowitz and Stegun1964) to map the wave amplitude perturbation variables so that they are centred on the steady-state droplet position at each impact.

Following Oza et al. (Reference Oza, Rosales and Bush2013), we consider general perturbations, and perturbations to the droplet in the direction of motion. In the former case, the walking is unstable. Physically, this corresponds to a new walker forming after an initial transient, but in a new direction. The latter case is achieved by noting that $\unicode[STIX]{x1D6F9}_{m}\equiv 0$ along the $x$ -axis, so the linearised jumps for the perturbation coefficients $\tilde{b}_{m}$ will be zero for an in-line perturbation. Furthermore, $\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D6F9}_{m}\equiv 0$ (for all $m$ ) along the $x$ -axis, so the $\tilde{b}_{m}$ terms do not contribute to the linearised droplet perturbation jump condition. Hence, the $\tilde{b}_{m}$ terms decouple from the perturbed system and can be neglected as they are stable solutions to the Mathieu equation. In this case, the walker is neutrally stable, where an eigenvalue of 1 always exists due to the translational invariance of the problem. This agrees with experiments and also the model of Oza et al. (Reference Oza, Rosales and Bush2013), whose analysis is restricted to lateral and in-line perturbations.

5.2 Slow-walking analysis

As observed by Protière et al. (Reference Protière, Boudaoud and Couder2006), there is a supercritical bifurcation at the walking threshold, where the walking speed grows like $(\unicode[STIX]{x1D6E4}-\unicode[STIX]{x1D6E4}_{W})^{1/2}$ . This behaviour has been captured in the trajectory equation of Oza et al. (Reference Oza, Rosales and Bush2013). We explore the slow-walking speed in terms of the distance between impacts $\unicode[STIX]{x1D6FF}x$ , which is the average walking speed in dimensionless variables. For $0<\unicode[STIX]{x1D6FF}x\ll 1$ , we consider asymptotic expansions of the form

(5.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{W}+\unicode[STIX]{x1D6FF}x\unicode[STIX]{x1D6E4}_{1}+\unicode[STIX]{x1D6FF}x^{2}\unicode[STIX]{x1D6E4}_{2}+\cdots \,, & \displaystyle\end{eqnarray}$$

(5.12) $$\begin{eqnarray}\displaystyle & \displaystyle a_{m}(t;k)=a_{m}^{(0)}(t;k)+\unicode[STIX]{x1D6FF}xa_{m}^{(1)}(t;k)+\unicode[STIX]{x1D6FF}x^{2}a_{m}^{(2)}(t;k)+\cdots \,, & \displaystyle\end{eqnarray}$$

for all $m\geqslant 0$ . By the Frobenius series for Bessel functions, we have

(5.13a ) $$\begin{eqnarray}\text{J}_{0}(k\unicode[STIX]{x1D6FF}x)\sim 1-\left(\frac{k}{2}\right)^{2}\unicode[STIX]{x1D6FF}x^{2}+O(\unicode[STIX]{x1D6FF}x^{4}),\end{eqnarray}$$

(5.13b ) $$\begin{eqnarray}\text{J}_{1}(k\unicode[STIX]{x1D6FF}x)\sim \left(\frac{k}{2}\right)\unicode[STIX]{x1D6FF}x-\frac{1}{2}\left(\frac{k}{2}\right)^{3}\unicode[STIX]{x1D6FF}x^{3}+O(\unicode[STIX]{x1D6FF}x^{5}),\end{eqnarray}$$

(5.14a,b ) $$\begin{eqnarray}\text{J}_{2}(k\unicode[STIX]{x1D6FF}x)\sim \frac{1}{2}\left(\frac{k}{2}\right)^{2}\unicode[STIX]{x1D6FF}x^{2}+O(\unicode[STIX]{x1D6FF}x)^{4},\quad \text{J}_{3}(k\unicode[STIX]{x1D6FF}x)\sim \frac{1}{6}\left(\frac{k}{2}\right)^{3}\unicode[STIX]{x1D6FF}x^{3}+O(\unicode[STIX]{x1D6FF}x^{5}).\end{eqnarray}$$

These expansions are valid for $k/2=O(1)$ . As the least-damped wavenumber $k_{F}$ satisfies $k_{F}/2\approx \unicode[STIX]{x03C0}=O(1)$ , and wavenumbers far away from $k_{F}$ correspond to negligible wave amplitudes, this assumption is valid. We prove that $\unicode[STIX]{x1D6E4}_{1}=0$ and verify numerically that $\unicode[STIX]{x1D6E4}_{2}>0$ , which gives the desired bifurcation.

Substituting the asymptotic expansions into the map (5.1)–(5.2), the Mathieu equation (3.22) and jump condition (3.25), and equating powers of $\unicode[STIX]{x1D6FF}x$ , a system of differential equations with initial and jump conditions is obtained, all of the general form

(5.15) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{L}}_{k}^{\unicode[STIX]{x1D6E4}_{W}}a(t;k)={\mathcal{H}},\quad a(1;k)=a(0;k)+{\mathcal{B}},\quad a^{\prime }(1^{+};k)=a^{\prime }(0^{+};k)+{\mathcal{B}}^{\prime },\\ \displaystyle \quad [a^{\prime }(1;k)]_{-}^{+}={\mathcal{J}}.\end{array}\right\}\end{eqnarray}$$

Here, ${\mathcal{H}}$ , ${\mathcal{B}}$ , ${\mathcal{B}}^{\prime }$ and ${\mathcal{J}}$ are inhomogeneities, and ${\mathcal{L}}_{k}^{\unicode[STIX]{x1D6E4}_{W}}$ is the damped Mathieu differential operator with $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{W}$ . By (5.3)–(5.4), the map for the wave amplitude derivatives is the same as that for the wave amplitudes themselves. Hence, all terms in ${\mathcal{B}}$ of the form $a(0;k)$ are replaced by $a^{\prime }(0^{+};k)$ in ${\mathcal{B}}^{\prime }$ , so for notational efficiency, we do not state ${\mathcal{B}}^{\prime }$ explicitly.

First, we note that for $\unicode[STIX]{x1D6FF}x=0$ , we have a bouncer at the walking threshold, so the wave field is radially symmetric, giving $a_{m}^{(0)}\equiv 0$ for $m\geqslant 1$ . This is verified by noting that ${\mathcal{H}}_{m}^{(0)}={\mathcal{B}}_{m}^{(0)}={\mathcal{J}}_{m}^{(0)}=0$ for all $m\geqslant 1$ , as the only order-1 term in the Bessel expansions is in $\text{J}_{0}$ . This yields a periodic unexcited solution to the damped Mathieu equation, which must be the zero solution.

To find $\unicode[STIX]{x1D6E4}_{W}$ , we first note that $a_{0}^{(0)}$ and $a_{1}^{(1)}$ are governed by inhomogeneities

(5.16a-c ) $$\begin{eqnarray}\displaystyle {\mathcal{H}}_{0}^{(0)}=0,\quad {\mathcal{B}}_{0}^{(0)}=0,\quad {\mathcal{J}}_{0}^{(0)}=-P_{0}(k),\end{eqnarray}$$

(5.17a-c ) $$\begin{eqnarray}{\mathcal{H}}_{1}^{(1)}=0,\quad {\mathcal{B}}_{1}^{(1)}=ka_{0}^{(0)}(0;k),\quad {\mathcal{J}}_{1}^{(1)}=-P_{1}(k)\,(k/2).\end{eqnarray}$$

Equating $\unicode[STIX]{x1D6FF}x$ from the walking condition (5.9), we require $\unicode[STIX]{x1D6E4}_{W}$ such that

(5.18) $$\begin{eqnarray}1=\frac{-\tilde{\unicode[STIX]{x1D708}}_{p}^{-1}(1-\text{e}^{-\tilde{\unicode[STIX]{x1D708}}_{p}})}{1-\text{e}^{-\tilde{\unicode[STIX]{x1D708}}_{p}}(1-F(c))}\frac{F(c)}{c}\sqrt{\frac{B}{R}}\int _{0}^{\infty }\frac{k^{2}}{2}a_{1}^{(1)}(0;k)\,\text{d}k.\end{eqnarray}$$

Therefore, (5.16)–(5.18) is a closed system for $a_{0}^{(0)}$ , $a_{1}^{(1)}$ and $\unicode[STIX]{x1D6E4}_{W}$ . This nonlinear analysis gives a consistent walking threshold $\unicode[STIX]{x1D6E4}_{W}$ to the linear stability analysis in §§ 4.1 and 5.1.

Given the solution to the above system, we may obtain equations for $\unicode[STIX]{x1D6E4}_{1}$ , which is coupled with equations for $a_{0}^{(1)}$ and $a_{1}^{(2)}$ . By vanishing inhomogeneities, $a_{2}^{(1)}\equiv 0$ , so

(5.19a-c ) $$\begin{eqnarray}{\mathcal{H}}_{0}^{(1)}=\unicode[STIX]{x1D6E4}_{1}\unicode[STIX]{x1D714}_{g}^{2}(k)\cos (4\unicode[STIX]{x03C0}t+\unicode[STIX]{x1D6FD})a_{0}^{(0)}(t;k),\quad {\mathcal{B}}_{0}^{(1)}=0,\quad {\mathcal{J}}_{0}^{(1)}=0,\end{eqnarray}$$

(5.20a-c ) $$\begin{eqnarray}{\mathcal{H}}_{1}^{(2)}=\unicode[STIX]{x1D6E4}_{1}\unicode[STIX]{x1D714}_{g}^{2}(k)\cos (4\unicode[STIX]{x03C0}t+\unicode[STIX]{x1D6FD})a_{1}^{(1)}(t;k),\quad {\mathcal{B}}_{1}^{(2)}=ka_{0}^{(1)}(0;k),\quad {\mathcal{J}}_{1}^{(2)}=0.\end{eqnarray}$$

The $\unicode[STIX]{x1D6FF}x^{2}$ term from (5.9) completes the system with

(5.21) $$\begin{eqnarray}0=\int _{0}^{\infty }\frac{k^{2}}{2}a_{1}^{(2)}(0;k)\,\text{d}k.\end{eqnarray}$$

We now exploit the directional symmetry of the steady walker to show that $\unicode[STIX]{x1D6E4}_{1}=0$ . The map (5.1)–(5.2) also holds for $\unicode[STIX]{x1D6FF}x<0$ , which corresponds to a droplet walking in the negative $x$ -direction. Since $a_{0}$ is the coefficient of the radially symmetric $\text{J}_{0}(kr)$ basis function, $a_{0}$ must be even in $\unicode[STIX]{x1D6FF}x$ for the wave field to be the same for both walking directions. Hence, we require $a_{0}^{(1)}\equiv 0$ on $(0,1)$ . By (5.19), this is achieved if and only if $\unicode[STIX]{x1D6E4}_{1}=0$ . With these vanishing terms, equation (5.20) becomes

(5.22a-c ) $$\begin{eqnarray}{\mathcal{H}}_{1}^{(2)}=0,\quad {\mathcal{B}}_{1}^{(2)}=0,\quad {\mathcal{J}}_{1}^{(2)}=0,\end{eqnarray}$$

yielding $a_{1}^{(2)}\equiv 0$ on $(0,1)$ , which satisfies (5.21). Neglecting terms of $o(\unicode[STIX]{x1D6FF}x^{2})$ , we have

(5.23) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}x=\unicode[STIX]{x1D6E4}_{2}^{-1/2}\sqrt{\unicode[STIX]{x1D6E4}-\unicode[STIX]{x1D6E4}_{W}},\end{eqnarray}$$

which is valid for $0<\unicode[STIX]{x1D6FF}x\ll 1$ . Equations for determining $\unicode[STIX]{x1D6E4}_{2}>0$ are given in appendix C. The solution is shown in figure 1(a) (grey curve).

5.3 Wave field of a single impact

Previous models for the wave field have superposed given radially symmetric sources centred at each impact with exponential temporal decay (Fort et al. Reference Fort, Eddi, Boudaoud, Moukhtar and Couder2010; Eddi et al. Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011; Moláček & Bush Reference Moláček and Bush2013b ; Oza et al. Reference Oza, Rosales and Bush2013). In our model, the jump conditions (3.25)–(3.26) at time $t_{n}$ are in fact equivalent to $[\unicode[STIX]{x1D702}(\boldsymbol{x},t_{n})]_{-}^{+}=\bar{\unicode[STIX]{x1D702}}(|\boldsymbol{x}-\boldsymbol{X}(t_{n})|,0)$ , where

(5.24) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D702}}(r,t)=\int _{0}^{\infty }k\bar{a}(t;k)\text{J}_{0}(kr)\,\text{d}k.\end{eqnarray}$$

For all $k\geqslant 0$ , $\bar{a}(t;k)$ is governed by ${\mathcal{L}}_{k}\bar{a}=0$ for $t>0$ , with initial conditions

(5.25a,b ) $$\begin{eqnarray}\bar{a}(0;k)=0,\quad \bar{a}^{\prime }(0;k)=-P_{0}(k),\end{eqnarray}$$

with $P_{0}(k)=kMG/(2\unicode[STIX]{x03C0})$ . This is quickly verified using Graf’s addition theorem (Abramowitz & Stegun Reference Abramowitz and Stegun1964). As $\bar{\unicode[STIX]{x1D702}}(r,0)=0$ for all $r$ , the wave field is continuous across the droplet impact, which is not true in the aforementioned models. Although $\bar{\unicode[STIX]{x1D702}}_{t}(r,0)=\infty$ for all $r\geqslant 0$ (from the $\unicode[STIX]{x1D6FF}(\cdot )$ -function forcing), the rapid temporal decay of $\bar{a}(t;k)$ for large $k$ gives a finite solution for any $t>0$ . Furthermore, $\bar{\unicode[STIX]{x1D702}}$ depends only on the current position of the droplet and is independent of its velocity and path memory.

The development of this wave field over time is shown in figure 2(a). A decaying capillary wave propagates away from the droplet, which excites a field of standing Faraday waves. The wave field generated is similar to that of Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015) near the origin, but differs slightly at the wave front over time. This dissimilarity may be explained partially by the assumption of instantaneous point impacts and the truncation of the wavenumbers. As the wave speed is much larger than the walking speed, this difference has a negligible effect on the resulting dynamics.

Figure 2. (a) The wave field of a single impact at times $1,\ldots ,8$ , which are normalised to have the same amplitude at the origin. A damped travelling capillary wave propagates from the impact, exciting a standing field of Faraday waves in its wake ( $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.994$ ). (b) Spatial decay length $l_{d}$ of a walker as a function of $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}$ and memory $M_{e}$ (inset (c)).

By superposition of such wave fields, we may now express the wave field $\unicode[STIX]{x1D702}$ by

(5.26) $$\begin{eqnarray}\unicode[STIX]{x1D702}(\boldsymbol{x},t)=\mathop{\sum }_{n=0}^{\left\lfloor t\right\rfloor }\bar{\unicode[STIX]{x1D702}}(|\boldsymbol{x}-\boldsymbol{X}(t_{n})|,t-t_{n})H(t-t_{n}),\end{eqnarray}$$

where $H(\cdot )$ is the Heaviside function. Following Moláček & Bush (Reference Moláček and Bush2013b ), we define the dimensionless memory $M_{e}(\unicode[STIX]{x1D6E4})$ , which is the number of impacts until $\bar{\unicode[STIX]{x1D702}}$ becomes exponentially small. Using the calculations of Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015), for $0<\unicode[STIX]{x1D6E4}_{F}-\unicode[STIX]{x1D6E4}\ll 1$ ,

(5.27) $$\begin{eqnarray}M_{e}(\unicode[STIX]{x1D6E4})=\frac{T_{d}(\unicode[STIX]{x1D6E4})}{1-\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}},\quad \text{where }T_{d}(\unicode[STIX]{x1D6E4})=\frac{16\unicode[STIX]{x1D716}(\unicode[STIX]{x1D714}^{2}(k_{c}(\unicode[STIX]{x1D6E4}))+\unicode[STIX]{x1D6FE}^{2}(k_{c}(\unicode[STIX]{x1D6E4}))+4\unicode[STIX]{x03C0}^{2})}{G^{2}\unicode[STIX]{x1D6E4}_{F}^{2}},\end{eqnarray}$$

is the dimensionless decay rate of the waves ( $\unicode[STIX]{x1D6E4}_{F}=5.13$ for this fluid model). The critical wavenumber $k_{c}$ ( $\unicode[STIX]{x1D6E4}$ ) is the least-damped wavenumber given $\unicode[STIX]{x1D6E4}$ . This dependence is weak, with a 1 % variation in $k_{c}(\unicode[STIX]{x1D6E4})$ away from $k_{c}(\unicode[STIX]{x1D6E4}_{F})=k_{F}$ for a range of $\unicode[STIX]{x1D6E4}$ , resulting in a small variation to $T_{d}(\unicode[STIX]{x1D6E4})$ (typically $T_{d}(\unicode[STIX]{x1D6E4})\approx 0.6$ ). This form is slightly more sophisticated than that of Moláček & Bush (Reference Moláček and Bush2013b ), where $T_{d}$ is assumed to be constant. From numerical verification, this is also a reasonable approximation for the range of $\unicode[STIX]{x1D6E4}$ considered.

From figure 2(a), the wave field approaches a Bessel-like function as time increases. As $k_{F}$ is the least-damped wavenumber, $\bar{a}(t;k)$ is peaked about $k=k_{F}$ for large $t$ . Hence, the integral for $\bar{\unicode[STIX]{x1D702}}(r,t)$ may be approximated by $\bar{\unicode[STIX]{x1D702}}(r,t)=a_{0}(t;k_{F})\text{J}_{0}(k_{F}r)$ for large $t$ . The function $a_{0}(t;k_{F})$ has a Floquet exponent with real part $-1/M_{e}$ , where the underlying periodic behaviour is well approximated by a sinusoid. This long-time approximation is used by Moláček & Bush (Reference Moláček and Bush2013b ) and Oza et al. (Reference Oza, Rosales and Bush2013), and prevents the occurrence of a travelling capillary wave formed at each impact and a Doppler shift. As our model generates a similar wave field to that of Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015), we refer the reader to their paper for a comprehensive comparison between wave field models.

The wavenumbers near to $k_{F}$ also prove to be significant in leading to the experimentally observed spatial damping (Eddi et al. Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011), and are computed numerically Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015). We repeat their statistical calculation for our model, with the resulting exponential decay length $l_{d}$ shown in figure 2(b). Specifically, we compute the wave field $\unicode[STIX]{x1D702}$ for a walker, and find the extrema to $|\unicode[STIX]{x1D702}|$ in the direction perpendicular to travel. As the far-field spatial decay of a Bessel function is ${\sim}1/\sqrt{r}$ , we fit a nonlinear model of the form $\unicode[STIX]{x1D703}_{1}\exp (-r/\unicode[STIX]{x1D703}_{2})/\sqrt{r}$ for parameters $(\unicode[STIX]{x1D703}_{1},\unicode[STIX]{x1D703}_{2})$ , where $l_{d}=\unicode[STIX]{x1D703}_{2}$ , which is an approximate envelope for the wave field. Typically, we consider $r\in (0,10]$ , which gives enough points for a reasonable statistical fit. As this approach contains some numerical error, the resulting solution is not smooth.

The spatial decay length $l_{d}$ increases with $\unicode[STIX]{x1D6E4}$ , and grows linearly with memory as $M_{e}\rightarrow \infty$ . However, due to the fixed impact phase, we do not obtain jumps in the decay length as the impacts switch between $(m,n)$ modes (Milewski et al. Reference Milewski, Galeano-Rios, Nachbin and Bush2015). As Eddi et al. (Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011) and Borghesi et al. (Reference Borghesi, Moukhtar, Labousse, Eddi, Fort and Couder2014) used fixed dimensionless values of $l_{d}=1.6$ and $2.5$ respectively, they did not capture this intrinsic change in the wave field as $\unicode[STIX]{x1D6E4}$ varies.

5.4 Doppler effect

As shown by Eddi et al. (Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011), a Doppler effect is observed in the wave field of the walker; the waves ahead of the droplet are compressed, yet are elongated behind. This phenomenon was observed in the simulations of Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015). Moreover, the authors proved that this cannot occur for wave amplitudes of the form $a(t;k)=\unicode[STIX]{x1D6FC}(t)\unicode[STIX]{x1D6FF}(k-k_{F})$ , which was used in Oza et al. (Reference Oza, Rosales and Bush2013). However, as we maintain a significant range of $k$ in our numerical solution, we observe a Doppler effect far from the droplet (figure 3). Our model exhibits a smooth yet nonlinear change in the wavelength due to the Doppler effect, which differs from the linear prediction of Eddi et al. (Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011).

Figure 3. Evidence of a Doppler effect. (a) Wave in direction of travel (black) and transverse wave (grey) recentred to have the same peaks. The arrow determines the direction of travel, where $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.9$ . Waves are elongated behind the droplet and compressed ahead. (b) Measured wavelengths behind $(+)$ and ahead $(\times )$ of the droplet. The grey lines are the theoretical predictions made by Eddi et al. (Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011).

Figure 4. Wave fields corresponding to different steady-state dynamics at droplet impact. The droplet is denoted by a blue circle for in-phase impacts and a red circle when in flight for antiphase dynamics. Walking droplets for (a) $\unicode[STIX]{x1D6FF}x=0.065$ with $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.9$ and (b) $\unicode[STIX]{x1D6FF}x=0.08$ with $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.96$ . Two anticlockwise orbiting droplets with $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.91$ for (c) in-phase and (d) antiphase impacts. A single droplet orbiting anticlockwise under a central force for $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.975$ with orbit radius (e) $R_{d}=0.45$ and (f) $R_{d}=0.95$ .

5.5 Summary of steady walking dynamics

To conclude this section, a bouncing droplet destabilises via a pitchfork bifurcation to a steady walker, whose lower-energy wave field automatically captures the exponential spatial damping correction observed in experiments (Eddi et al. Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011). The jump conditions are equivalent to adding a radially symmetric propagating wave field $\bar{\unicode[STIX]{x1D702}}$ at each impact, whose temporal equidistant superposition yields a Doppler shift, with a weakly nonlinear dependence on $\unicode[STIX]{x1D6FF}x$ . Furthermore, the temporal decay of the system may be described by its memory $M_{e}(\unicode[STIX]{x1D6E4})$ , which allows for comparison between previous models (Fort et al. Reference Fort, Eddi, Boudaoud, Moukhtar and Couder2010; Moláček & Bush Reference Moláček and Bush2013b ). The wave fields corresponding to the dynamics of a walker are easily computed and are shown in figures 4(a) and 4(b) for two different walking speeds. For large $\unicode[STIX]{x1D6E4}$ , an interference pattern forms behind the droplet (Eddi et al. Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011).

6.1 In-phase orbiting

We consider two orbiting droplets about the origin; both droplets impact simultaneously on a radius $R_{d}$ , with the droplets and wave field rotating by an angle $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}>0$ between impacts. Moreover, we pose that the angular difference between $\boldsymbol{X}^{(1)}(t_{n})$ and $\boldsymbol{X}^{(2)}(t_{n})$ is $\unicode[STIX]{x03C0}$ for all $t_{n}$ . This allows us to model $\boldsymbol{X}^{(1)}(t)$ explicitly, with the contribution of $\boldsymbol{X}^{(2)}(t_{n})$ to the wave amplitude jump conditions treated implicitly. For notational convenience, we write $\boldsymbol{X}(t)\equiv \boldsymbol{X}^{(1)}(t)$ , where $\boldsymbol{X}(t_{n})$ has radial component $r_{d}(t_{n})=R_{d}$ and angular component $\unicode[STIX]{x1D703}_{d}(t_{n})=\unicode[STIX]{x1D703}_{n}\equiv n\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}$ (by rotational invariance). The droplet motion is piecewise linear, so $r_{d}(t)\neq R_{d}$ for $t\neq t_{n}$ .

As the implicit second droplet has angular component $(\unicode[STIX]{x1D703}_{n}+\unicode[STIX]{x03C0})$ at time $t_{n}$ , (6.1) gives

(6.2) $$\begin{eqnarray}\displaystyle [a_{m}^{\prime }(t_{n};k)]_{-}^{+} & = & \displaystyle -P_{m}(k)\unicode[STIX]{x1D6F7}_{m}(\boldsymbol{X}(t_{n});k)(1+(-1)^{m}),\end{eqnarray}$$

(6.3) $$\begin{eqnarray}\displaystyle [b_{m}^{\prime }(t_{n};k)]_{-}^{+} & = & \displaystyle -P_{m}(k)\unicode[STIX]{x1D6F9}_{m}(\boldsymbol{X}(t_{n});k)(1+(-1)^{m}).\end{eqnarray}$$

For $m$ odd, $a_{m}$ and $b_{m}$ are never excited during the orbit, giving $a_{m}\equiv b_{m}\equiv 0$ . This ensures that the wave field with respect to each droplet is the same.

Denoting $\boldsymbol{c}_{m}(t;k)=(a_{m}(t;k),b_{m}(t;k))^{\text{T}}$ , the wave field rotates with the droplets if

(6.4a,b ) $$\begin{eqnarray}\boldsymbol{c}_{m}(t_{n+1};k)=\unicode[STIX]{x1D63F}(m\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})\boldsymbol{c}_{m}(t_{n};k),\quad \boldsymbol{c}_{m}^{\prime }(t_{n+1}^{+};k)=\unicode[STIX]{x1D63F}(m\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})\boldsymbol{c}_{m}^{\prime }(t_{n}^{+};k),\end{eqnarray}$$

where $\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D711})\in \mathbb{R}^{2\times 2}$ is the rotation matrix for angle $\unicode[STIX]{x1D711}$ . The droplet rotation requires

(6.5a,b ) $$\begin{eqnarray}\boldsymbol{X}(t_{n+1})=\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})\boldsymbol{X}(t_{n}),\quad \boldsymbol{X}^{\prime }(t_{n+1}^{+})=\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})\boldsymbol{X}^{\prime }(t_{n}^{+}).\end{eqnarray}$$

As with steady walking, we solve the droplet dynamics (3.24)–(3.27) with $\tilde{\unicode[STIX]{x1D705}}=0$ for $t\in [t_{n},t_{n+1})$ subject to (6.5), and eliminate $\boldsymbol{X}^{\prime }(t_{n}^{+})$ in favour of $\boldsymbol{X}(t_{n})$ . This gives

(6.6) $$\begin{eqnarray}(\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})-(1-F(c))\text{e}^{-\tilde{\unicode[STIX]{x1D708}}_{p}}\unicode[STIX]{x1D644}_{2})(\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})-\unicode[STIX]{x1D644}_{2})\boldsymbol{X}(t_{n})=\frac{-F(c)}{c}\frac{1-\text{e}^{-\tilde{\unicode[STIX]{x1D708}}_{p}}}{\tilde{\unicode[STIX]{x1D708}}_{p}}\sqrt{\frac{B}{R}}\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})(\unicode[STIX]{x1D735}\unicode[STIX]{x1D702})_{n},\end{eqnarray}$$

where $(\unicode[STIX]{x1D735}\unicode[STIX]{x1D702})_{n}\equiv \unicode[STIX]{x1D735}\unicode[STIX]{x1D702}(\boldsymbol{X}(t_{n}),t_{n})$ , and $(\unicode[STIX]{x1D735}\unicode[STIX]{x1D702})_{n+1}=\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})(\unicode[STIX]{x1D735}\unicode[STIX]{x1D702})_{n}$ due to the rotation.

We need to find $R_{d}$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}$ such that the wave field amplitudes governed by (3.22)–(3.23) with jump conditions (6.2)–(6.3) satisfy the rotation conditions (6.4) and droplet condition (6.6). By periodicity, it is sufficient to solve for $n=0$ (so $t\in [0,1]$ ), where $\boldsymbol{X}(t_{0})=(R_{d},0)^{\text{T}}$ . This is formulated as many $4\times 4$ linear systems subject to (6.6), akin to finding the walking states. An example in-phase wave field is given in figure 4( $c$ ).

The stability analysis is similar to that of the walking states with two independent droplet perturbations, although we rotate the domain by $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}$ between iterations so that the steady-state droplet always lies at $(r,\unicode[STIX]{x1D703})=(R_{d},0)$ in the current domain. It should be noted that the wave amplitudes $a_{m}$ and $b_{m}$ are required for all $m$ (not just even indices), and the system is truncated for sufficiently large $m$ given $R_{d}$ , by the shape of the large-order Bessel functions near zero (it should be recalled that $\text{J}_{m}(z)\sim (1/m!)(z/2)^{m}$ for $0<z\ll \sqrt{m+1}$ ).

Figure 5. (a) Bifurcation diagram for the diameters $D$ of in-phase $I_{n}$ (light grey background) and antiphase $A_{n}$ (white background) orbiting pairs. The dark grey regions give the bounds for two droplets in a stable wobbly orbit. This has an upper bound given by the unstable upper branch. The curves correspond to stability types in table 2. (b) Example wobbly in-phase anticlockwise orbiters with $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.83$ , which corresponds to the black cross in (a).

The orbit diameters $D\equiv 2R_{d}$ are shown in figure 5(a), where curves with a light grey background correspond to in-phase orbiters. The diameters of the stable solutions are in good agreement with the experimental values recorded by Couder et al. (Reference Couder, Protière, Fort and Boudaoud2005b ), namely $D_{n}=(n-\bar{\unicode[STIX]{x1D716}})$ , where $\bar{\unicode[STIX]{x1D716}}\approx 0.2$ and $n$ is an integer. However, the co-existence of stable $n=1,2,3$ orbits occurs only for a limited range of $\unicode[STIX]{x1D6E4}$ . This system also contains wobbly orbits, which appear when the steady orbit destabilises via a complex conjugate pair of eigenvalues. Due to nonlinear effects, the droplets remain locked in a periodic wobbly orbit (see figure 5(b)). The minimum and maximum distances apart reached by the droplets are indicated by the dark grey regions in figure 5(a), and are obtained by simulating away from the steady states. We observe that the maximum distance apart is bounded above by the corresponding unstable upper branch to each solution, which sheds some light on why wobbly solutions only exist over a limited range of $\unicode[STIX]{x1D6E4}$ .

For further insight, we consider the solution at $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{F}$ with wavenumber $k=k_{F}$ . The Floquet exponents are $\unicode[STIX]{x1D707}_{1}=0$ and $\unicode[STIX]{x1D707}_{2}=-2\unicode[STIX]{x1D6FE}(k_{F})$ , giving the eigenvalues of the Mathieu fundamental matrix $\unicode[STIX]{x1D648}_{k_{F}}(\unicode[STIX]{x1D6E4}_{F})$ as $\unicode[STIX]{x1D70C}_{1}\equiv \exp (\unicode[STIX]{x1D707}_{1})=1$ and $\unicode[STIX]{x1D70C}_{2}\equiv \exp (\unicode[STIX]{x1D707}_{2})$ . Periodicity requires

(6.7) $$\begin{eqnarray}(\unicode[STIX]{x1D644}_{2}-\unicode[STIX]{x1D648}_{k_{F}}(\unicode[STIX]{x1D6E4}_{F}))\left(\begin{array}{@{}c@{}}a_{0}(0;k_{F})\\ a_{0}^{\prime }(0;k_{F})\end{array}\right)=-\left(\begin{array}{@{}c@{}}0\\ 2P_{0}(k_{F})\text{J}_{0}(k_{F}R_{d})\end{array}\right),\end{eqnarray}$$

where $(\unicode[STIX]{x1D644}_{2}-\unicode[STIX]{x1D648}_{k_{F}}(\unicode[STIX]{x1D6E4}_{F}))$ is singular. A candidate radius $R_{d}$ satisfies $\text{J}_{0}(k_{F}R_{d})=0$ , with a corresponding wave amplitude vector from the matrix null space. In fact, such radii satisfy the full problem and correspond to the lower branches in figure 5(a) as $\unicode[STIX]{x1D6E4}\rightarrow \unicode[STIX]{x1D6E4}_{F}$ .

6.2 Antiphase orbiting

For antiphase orbiting of two droplets, we impose that $\boldsymbol{X}^{(1)}(t_{n})$ and $\boldsymbol{X}^{(2)}(t_{n+1/2})$ both lie on the radius $R_{d}$ with angular components $\unicode[STIX]{x1D703}_{n}$ and $(\unicode[STIX]{x1D703}_{n+1/2}+\unicode[STIX]{x03C0})$ respectively, where $t_{n+1/2}\equiv t_{n}+1/2$ and $\unicode[STIX]{x1D703}_{n+1/2}\equiv \unicode[STIX]{x1D703}_{n}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}/2.$ We construct the wave field to rotate by $\unicode[STIX]{x03C0}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}/2$ in half of an impact period, which ensures that the wave field is the same with respect to each droplet at impact. The antiphase jump conditions are

(6.8) $$\begin{eqnarray}\displaystyle & \displaystyle [a_{m}^{\prime }(t_{n+1/2};k)]_{-}^{+}=-P_{m}(k)(-1)^{m}\unicode[STIX]{x1D6F7}_{m}(R_{d},\unicode[STIX]{x1D703}_{n+1/2};k), & \displaystyle\end{eqnarray}$$

(6.9) $$\begin{eqnarray}\displaystyle & \displaystyle [b_{m}^{\prime }(t_{n+1/2};k)]_{-}^{+}=-P_{m}(k)(-1)^{m}\unicode[STIX]{x1D6F9}_{m}(R_{d},\unicode[STIX]{x1D703}_{n+1/2};k). & \displaystyle\end{eqnarray}$$

Hence, both odd and even $m$ are required for antiphase orbiters.

It remains to find $R_{d}$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}$ such that (6.6) is satisfied. This is coupled with Mathieu equations (3.22)–(3.23), jump conditions (6.8)–(6.9) and wave amplitude maps (6.4), but with $t_{n+1}\mapsto t_{n+1/2}$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}\mapsto \unicode[STIX]{x03C0}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}/2$ . An example wave field is shown in figure 4( $d$ ). The stability analysis requires a two-stage transition map to account for the antiphase impacts.

The quantisation of antiphase orbiters is given by the regions of white background shown in figure 5(a). Although the reported orbit quantisation $D_{n}^{\prime }=(n+1/2-\bar{\unicode[STIX]{x1D716}})$ is obtained (Couder et al. Reference Couder, Protière, Fort and Boudaoud2005b ), the range of stable solutions with $\unicode[STIX]{x1D6E4}$ is limited, particularly for the smallest orbit. We note that the stability and existence of solutions are dependent on the skidding friction $c$ and impact phase $\unicode[STIX]{x1D6FD}$ , which are both fixed for this theoretical investigation. We expect that these results will be particularly sensitive to the impact phase.

6.3 Promenade pairs

For simplicity, we only analyse rectilinear promenade pairs. To exploit the symmetry of two in-phase droplets, we fix the droplets to walk in parallel along lines $y=\pm \unicode[STIX]{x1D6FF}y$ in the direction of increasing $x$ . To find $\unicode[STIX]{x1D6FF}y$ , we require $\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D702}=0$ at each drop, and the walking speed $\unicode[STIX]{x1D6FF}x$ is given by an equation analogous to (5.9). This couples with the wave amplitude maps (5.3)–(5.4) and multiple-droplet jump conditions (6.1), with $\boldsymbol{X}^{(1)}(t_{n})=(n\unicode[STIX]{x1D6FF}x,\unicode[STIX]{x1D6FF}y)$ and $\boldsymbol{X}^{(2)}(t_{n})=(n\unicode[STIX]{x1D6FF}x,-\unicode[STIX]{x1D6FF}y)$ . The stability analysis follows similarly to the walking case, yet the two droplets have independent perturbations.

Figure 6. Bifurcation diagram for (a) promenade pairs and (b) two-droplet trains as a function of distance apart $D$ and speed relative to a single walker $\unicode[STIX]{x1D6FF}x/\unicode[STIX]{x1D6FF}x_{1}$ . Grey and white backgrounds denote in-phase and antiphase dynamics respectively. The thin grey lines connect points of equal $\tilde{\unicode[STIX]{x1D6E4}}\equiv \unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}$ . The thick lines correspond to stability types in table 2.

The extension to antiphase promenaders is analogous to that for orbiters. In half of an impact period, we require the wave field to be translated $\unicode[STIX]{x1D6FF}x/2$ along the $x$ -axis and reflected about $y=0$ . This ensures that each droplet receives the same kick at impact.

Results are shown in figure 6(a), where the grey and white backgrounds correspond to in-phase and antiphase dynamics, respectively. We first observe that the speed of a promenading pair is less than that of a single walker, as reported by Borghesi et al. (Reference Borghesi, Moukhtar, Labousse, Eddi, Fort and Couder2014). The quantised distance between droplets $D\equiv 2\unicode[STIX]{x1D6FF}y$ depends weakly on $\unicode[STIX]{x1D6E4}$ , which was not reported in Borghesi et al. (Reference Borghesi, Moukhtar, Labousse, Eddi, Fort and Couder2014). Instead, the authors stated that $D$ has approximate in-phase values of $0.6,1.6,2.6,\ldots$ and antiphase values of $1.1,2.1,3.1,\ldots$ , which is consistent with our findings.

Crucially, we note that straight-line promenade pairs are unstable for all quantisations and all values of $\unicode[STIX]{x1D6E4}$ . When there is an oscillatory instability (with no real eigenvalue), we expect the droplets to become bound in parallel motion, but with a transverse oscillation, as reported in experiments (Protière et al. Reference Protière, Boudaoud and Couder2006; Borghesi et al. Reference Borghesi, Moukhtar, Labousse, Eddi, Fort and Couder2014). This is demonstrated in figure 7, where we simulate from the unstable steady state (with no real unstable eigenvalues), and nonlinear effects keep the droplets in a steady oscillatory bound state. Such a result was speculated on by Borghesi et al. (Reference Borghesi, Moukhtar, Labousse, Eddi, Fort and Couder2014), but is hard to observe experimentally due to the confines of a finite domain. The size of these oscillations depends on the skidding friction $c$ and phase shift $\unicode[STIX]{x1D6FD}$ free parameters, which would both be tuned for a more careful comparison with experimental data. Larger $c$ increases the solution space for oscillatory pairs, and the transverse extent of the oscillations is reduced.

To shed light on the instability of straight-line promenaders, we construct a new wave field based on the superposition of two parallel walkers (from § 5), and compare the average wave field energy over one impact period (appendix B). The distance between each walker is given by the continuous parameter $D$ . The antiphase case requires a mid-period impact for the second droplet; this presents a minor difficulty as we do not know the speed of the newly constructed walking pair, so the impact location is unknown (it should be recalled that a promenading pair travels slower than a single walker). For ease, we assume that the pair of walkers travel at the same speed as a single walker, which determines the impact location. We note that any error for this impact is small, and as the energy is a quadratic function of the wave amplitudes, this gives a negligible change to the energy.

Results of this calculation are shown in figure 8. Intuitively, we expect the energy of the walking pair to be minimised at the same distance as that of the promenading pair. However, the walking pair do not obey the condition $\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D702}=0$ at each droplet, which is required for a parallel procession. Hence, additional energy must be stored in the wave field in the case of the quantised promenaders, whose energy range (with $\unicode[STIX]{x1D6E4}$ ) is denoted by the thick black lines in figure 8. The difference in energy decreases as $D$ increases due to the spatial decay of the wave field of the walker. Physically, the wave field would tend to adopt the lowest-energy state of two walkers; however, as this motion cannot be rectilinear (since $\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D702}\neq 0$ at each droplet impact), a transverse oscillation ensues.

Figure 7. Example transition from an unstable straight-line promenade (starting at $x=0$ ) to a stable oscillating promenade, with the speed given by the grey scale bar.

Figure 8. Wave field energy $E/E_{W}$ (relative to a single walker) for two in-phase (black) and antiphase (grey) parallel walkers at a distance $D$ apart for different values of $\unicode[STIX]{x1D6E4}$ . The larger values of $\unicode[STIX]{x1D6E4}$ give the most extreme variations in energy ( $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.88,0.9,\ldots ,0.98$ ). The thick black lines give the energy of the corresponding quantised parallel promenade solutions. The required interaction energy for the promenade mode partially explains its instability.

6.4 Droplet trains

Trains of droplets may be studied in a similar fashion to promenades, except that all droplets lie on the $x$ -axis and are spaced $\unicode[STIX]{x1D6FF}s$ apart. This requires constant $\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D702}$ at each droplet, with the in-phase and antiphase results shown in figure 6(b), where we denote $D\equiv \unicode[STIX]{x1D6FF}s$ . For sufficiently large $\unicode[STIX]{x1D6E4}$ , two droplets may walk faster than a single droplet, which was also observed for experiments in a confined annulus (Filoux, Hubert & Vandewalle Reference Filoux, Hubert and Vandewalle2015). The stability types reported in figure 6(b) are for general perturbations, although the droplet trains are neutrally stable in the direction of travel (like a single walking droplet). Interestingly, the presence of a second droplet stabilises the system to general perturbations for some antiphase trains. A possible extension is to consider trains of multiple droplets, as considered in an annulus by Filoux et al. (Reference Filoux, Hubert and Vandewalle2015).

7.1 Circular orbits

For a single droplet in a circular orbit, the wave amplitudes must satisfy the rotational maps (6.4a,b ), and are governed by (3.22)–(3.23) with jump conditions (3.25)–(3.26). However, as $\tilde{\unicode[STIX]{x1D705}}>0$ , the droplet rotation condition (6.6) becomes

(7.1) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle F(c)\left(\frac{\unicode[STIX]{x1D707}_{+}\text{e}^{\unicode[STIX]{x1D707}_{+}}-\unicode[STIX]{x1D707}_{-}\text{e}^{\unicode[STIX]{x1D707}_{-}}}{\unicode[STIX]{x1D707}_{+}-\unicode[STIX]{x1D707}_{-}}\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})-\text{e}^{-\tilde{\unicode[STIX]{x1D708}}_{p}}\unicode[STIX]{x1D644}_{2}\right)\boldsymbol{X}(t_{n})+(\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})-\text{e}^{\unicode[STIX]{x1D707}_{+}}\unicode[STIX]{x1D644}_{2})(\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})-\text{e}^{\unicode[STIX]{x1D707}_{-}}\unicode[STIX]{x1D644}_{2})\boldsymbol{X}(t_{n})\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad =\frac{-F(c)}{c}\sqrt{\frac{B}{R}}\left(\frac{\text{e}^{\unicode[STIX]{x1D707}_{+}}-\text{e}^{\unicode[STIX]{x1D707}_{-}}}{\unicode[STIX]{x1D707}_{+}-\unicode[STIX]{x1D707}_{-}}\right)\unicode[STIX]{x1D63F}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})(\unicode[STIX]{x1D735}\unicode[STIX]{x1D702})_{n},\end{eqnarray}$$

where

(7.2) $$\begin{eqnarray}\unicode[STIX]{x1D707}_{\pm }\equiv {\textstyle \frac{1}{2}}(-\tilde{\unicode[STIX]{x1D708}}_{p}\pm \sqrt{\tilde{\unicode[STIX]{x1D708}}_{p}^{2}-4\tilde{\unicode[STIX]{x1D705}}}).\end{eqnarray}$$

When $\unicode[STIX]{x1D707}_{\pm }\in \mathbb{C}$ , we use Euler’s formula to ensure a real solution, whose stability is analysed analogously to § 6.1. Example wave fields are given in figure 4(e,f).

We characterise this system with $\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D6E4})\equiv \unicode[STIX]{x1D6FF}x(\unicode[STIX]{x1D6E4})/\sqrt{\tilde{\unicode[STIX]{x1D705}}}$ , where $\unicode[STIX]{x1D6FF}x(\unicode[STIX]{x1D6E4})$ is the steady walking speed for $\tilde{\unicode[STIX]{x1D705}}=0$ . For a classical orbit, balancing the spring constant $\tilde{\unicode[STIX]{x1D705}}$ and centripetal acceleration gives the orbit radius $R_{d}=\unicode[STIX]{x1D6EC}$ . However, the droplet interaction with the wave field leads to a more complicated radius structure as $\unicode[STIX]{x1D6E4}$ increases (figure 9), where the unstable orbits require a larger average wave field energy. As with the in-phase orbiters, radii $R_{d}$ satisfying $\text{J}_{0}(k_{F}R_{d})=0$ are solutions when $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{F}$ , which approximately correspond to the upper plateaus in figure 9(b). Noting that $\unicode[STIX]{x1D6EC}=\infty$ corresponds to $\tilde{\unicode[STIX]{x1D705}}=0$ , figure 9(b) suggests that orbital solutions in the absence of a central force are possible for sufficiently high $\unicode[STIX]{x1D6E4}$ . In the present parameter regime, self-orbiters are not stable; simulation from the self-orbiter initial conditions ( $R_{d}\approx 0.42$ ) yields only two complete orbits at $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.994$ before the trajectory diverges.

When the circular orbits destabilise, a range of new orbits appear, with common examples shown in figure 10. In § 7.2, we see how these orbits manifest themselves at even larger $\unicode[STIX]{x1D6E4}$ , where all orbits are unstable.

Figure 9. Bifurcation diagram for orbit radius $R_{d}$ of a droplet in a harmonic potential well of width $\unicode[STIX]{x1D6EC}$ . In (a), $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.916$ (almost linear) and $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.975$ (snaked curve). In (b), $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.994$ . The lines correspond to stability types in table 2.

Figure 10. Example trajectories (lemniscate, trefoil and butterfly) for $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}\approx 0.975$ and different values of $\unicode[STIX]{x1D6EC}$ . The colour bars indicate the speed of the droplet.

7.2 Towards a double quantisation

As shown in the experiments of Perrard et al. (Reference Perrard, Labousse, Miskin, Fort and Couder2014b ), these exotic orbits exhibit a double quantisation in their mean radius and angular momentum in their stable states. The quantised states follow the same selection rule as the energy and angular momentum of a quantum particle in a two-dimensional potential well. Specifically, the radius and angular momentum form a $(p,q)$ lattice, where $p$ is an integer and $q\in \{-p,-p+2,\ldots ,p-2,p\}$ . These points correspond to the grey crosses $(\times )$ in figures 11 and 12.

Although we obtain a similar quantisation for stable orbits, the methods of Perrard et al. (Reference Perrard, Labousse, Miskin, Fort and Couder2014b ) rely on a range of $\unicode[STIX]{x1D6E4}$ and $\unicode[STIX]{x1D6EC}$ in an unsystematic way. Instead, we explore the dynamics in the chaotic regime for fixed $\unicode[STIX]{x1D6E4}$ , whereby the observed orbits destabilise but continue to be recognisable over short time intervals. This allows for a novel and methodical analysis of the long-time dynamics, yielding a similar double quantisation.

Figure 11. (a) Example clustering for a single trajectory with $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.984$ and $\unicode[STIX]{x1D6EC}=1.23$ . The light blue circles are averages over sections and the grey crosses form the approximate lattice for the quantum double quantisation. The red dots are the cluster centroids for this example (initially, $R_{d}=1.75$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}>0$ ). (b) The corresponding trajectory (grey) with example coloured sub-trajectories shown over short time intervals.

Figure 12. Double quantisation for a droplet in a harmonic potential for $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.984$ . The grey crosses form the quantum double quantisation lattice and the red dots are the cluster centroids for simulations over all initial radii $R_{d}\in \{0.4,0.425,\ldots ,2\}$ and corresponding values of $\unicode[STIX]{x1D6EC}$ , with both $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}>0$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}<0$ initially.

For a given circular orbit radius $R_{d}>0$ , there exist unique $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}>0$ , $\unicode[STIX]{x1D6EC}>0$ and the corresponding wave field (see the example in figure 9(a) with $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.975$ ). This provides a well-defined systematic initial condition for numerical simulation. We simulate over $(N_{0}+N)$ impacts, where the first $N_{0}$ impacts are a ‘burn-in’ period (omitted from further analysis). This is the time allowed for the droplet to destabilise from its circular orbit and begin to explore the underlying probability distribution of the dynamics. Throughout this work, we fix $N_{0}=3000$ impacts and $N=3\times 10^{4}$ , where $N_{0}$ is estimated by considering the unstable eigenvalues of the system and $N$ is large enough so that the resulting probability distribution is ergodic. For stable circular orbits, the droplet remains on its initial trajectory for all impacts. We perform this process for $R_{d}\in \{0.4,0.425,\ldots ,2\}$ , which removes the stable orbits for small $\unicode[STIX]{x1D6EC}$ in which the central force dominates the trajectory. Simulation from initial conditions with $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}<0$ does not yield a corresponding symmetric trajectory (as the dynamics are chaotic), but the statistics are unchanged.

It remains to analyse the data from these trajectories. Following Perrard et al. (Reference Perrard, Labousse, Miskin, Fort and Couder2014b ), we define the dimensionless mean radius $\bar{R}$ and mean angular momentum $\bar{L}_{z}$ as

(7.3a,b ) $$\begin{eqnarray}\bar{R}=\frac{1}{\sqrt{N}}\left(\mathop{\sum }_{n=1}^{N}r_{d}(t_{n+N_{0}})^{2}\right)^{1/2},\quad \bar{L}_{z}=\frac{1}{N}\mathop{\sum }_{n=1}^{N}\left(\boldsymbol{X}(t_{n+N_{0}})\times \frac{1}{\unicode[STIX]{x1D6FF}x(\unicode[STIX]{x1D6E4})}\boldsymbol{X}^{\prime }(t_{n+N_{0}}^{+})\right)\boldsymbol{\cdot }\boldsymbol{e}_{z},\end{eqnarray}$$

where $r_{d}(t)$ is the droplet radius at time $t$ and $\times$ denotes the vector product.

When simulating from an oscillatory instability, we expect to obtain wobbling circular orbits. However, when a real unstable eigenvalue exists, the trajectory is pushed away from the circular orbit, and more exotic dynamics occur. From the work of Perrard et al. (Reference Perrard, Labousse, Miskin, Fort and Couder2014b ), we expect to obtain trajectories that resemble lemniscates and trefoils on a short time scale, but appear chaotic on a long time scale (see the examples in figure 10). Although taking averages over all $N$ impacts makes some progress towards a double quantisation, trajectories with switching states (typically between lemniscates and ovals) blur this picture. As these states are non-constant and unstable, it is difficult to identify switching times without hand-picking the data (as in Perrard et al. (Reference Perrard, Labousse, Miskin, Fort and Couder2014b )), so we seek an alternative approach.

The answer lies in a further understanding of stable lemniscate and trefoil orbits. Our systematic analysis starts off by segmenting the long trajectory into shorter sections between the radial maxima of impacts. This yields well-defined cutting times. For a fictitious stable lemniscate, the mean radius and angular momentum are the same between two radial maxima compared with the whole trajectory (we have ‘half’ of a lemniscate). Similarly, the radial maxima break a fictitious stable trefoil into three identical sections (up to a rotation of $2\unicode[STIX]{x03C0}/3$ ). For a stable circular orbit, each point is defined as a radial maximum. Hence, we proceed by splitting each trajectory into intervals between radial maxima, and compute $\bar{R}$ and $\bar{L}_{z}$ over each interval.

For each trajectory, we may visualise the resulting interval mean data by plotting a scatter graph in the $(\bar{L}_{z},\bar{R})$ -plane (see the example in figure 11). Each trajectory gives rise to a well-separated set of clusters of mean radius and angular momentum points, where multiple clusters only arise in the case of switching states. This motivates cluster analysis as a suitable statistical approach to completing the double quantisation.

To calculate the clusters for each simulation, we use $K$ -means clustering, where $K$ is a predetermined number of clusters (Hartigan Reference Hartigan1975; Hartigan & Wong Reference Hartigan and Wong1979). The aim is to globally minimise the cost function, which is the total Euclidean distance between points. Each sub-trajectory is weighted by its number of impacts; this prevents smaller trajectories from being favoured (as radial maxima are achieved more frequently). The algorithm takes random initial locations for the cluster means (centroids) and is guaranteed to converge to a local minimum of the cost function. A sufficiently large number of repetitions is likely to locate the global minimum, which determines the locations of the centroids. When $K=1$ , this is simply the global mean of the data.

The issue with all clustering techniques is determining the value of $K$ . As it happens, the cases considered for this double quantisation can almost be identified by eye; the example in figure 11(a) has $K=5$ well-separated clusters. However, statistical techniques also exist as a guide for the number of clusters, such as the elbow method, silhouettes and cross-validation (Kodinariya & Makwana Reference Kodinariya and Makwana2013). However, the great benefit of clustering in this application is that the technique automatically groups together similar trajectories and averages over them. The outlier trajectories (say the transitions between lemniscates and ovals) are infrequent so do not skew the data.

By combining the centroids from all of the simulations, we obtain a double quantisation in $(\bar{L}_{z},\bar{R})$ in the chaotic regime ( $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.984$ ), as shown in figure 12. We note that Perrard et al. (Reference Perrard, Labousse, Miskin, Fort and Couder2014b ) also observed a less sharp quantisation for $\bar{R}>2$ . There is a notable dearth of centroids around $(\pm 0.5,1.5)$ ; the loop tractories corresponding to this quantisation are unstable and only exist for much larger values of $\unicode[STIX]{x1D6E4}$ considered. The loops are the combination of straight-line solutions and spin states, whereby a single droplet may orbit itself in the absence of a central force. The existence of these unstable spin states can be seen from figure 9(b) for large $\unicode[STIX]{x1D6EC}$ (corresponding to small $\tilde{\unicode[STIX]{x1D705}}$ ). However, the centroids at approximately $(\pm 1,1.5)$ are stable trefoil-like features not reported in Perrard et al. (Reference Perrard, Labousse, Miskin, Fort and Couder2014b ). A more thorough investigation may lead to the discovery of new trajectories that alter the shape of the double quantisation.