1. Introduction
On vertically vibrating a bath of silicone oil at frequency 
$f$, a droplet of the same oil can be made to bounce indefinitely on the free surface of the liquid (Walker Reference Walker1978; Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005a). As the amplitude of the forcing is increased, the bouncing droplet destabilises and transitions to a steady walking state (Couder et al. Reference Couder, Protiere, Fort and Boudaoud2005b). The walking droplet, also called a ‘walker’, emerges just below the Faraday instability threshold (Faraday Reference Faraday1831), above which the whole surface becomes unstable to standing Faraday waves oscillating at the subharmonic frequency 
$f/2$. On each bounce, the walker generates a localised damped Faraday wave on the fluid surface. It then interacts with these waves on subsequent bounces, giving rise to a self-propelled wave–droplet entity. Intriguingly, such walkers have been shown to mimic several peculiar behaviours that were previously thought to be exclusive to the quantum world. These include orbital quantisation in rotating frames (Fort et al. Reference Fort, Eddi, Boudaoud, Moukhtar and Couder2010; Harris & Bush Reference Harris and Bush2014; Oza et al. Reference Oza, Harris, Rosales and Bush2014) and harmonic potentials (Perrard et al. Reference Perrard, Labousse, Fort and Couder2014a,Reference Perrard, Labousse, Miskin, Fort and Couderb; Labousse et al. Reference Labousse, Oza, Perrard and Bush2016), Zeeman splitting in rotating frames (Eddi et al. Reference Eddi, Moukhtar, Perrard, Fort and Couder2012; Oza, Rosales & Bush Reference Oza, Rosales and Bush2018), wave-like statistical behaviour in confined geometries (Harris et al. Reference Harris, Moukhtar, Fort, Couder and Bush2013; Gilet Reference Gilet2016; Cristea-Platon, Sáenz & Bush Reference Cristea-Platon, Sáenz and Bush2018; Sáenz, Cristea-Platon & Bush Reference Sáenz, Cristea-Platon and Bush2018; Durey, Milewski & Wang Reference Durey, Milewski and Wang2020) as well as in an open system (Sáenz, Cristea-Platon & Bush Reference Sáenz, Cristea-Platon and Bush2020) and tunnelling across submerged barriers (Eddi et al. Reference Eddi, Fort, Moisy and Couder2009; Nachbin, Milewski & Bush Reference Nachbin, Milewski and Bush2017; Tadrist et al. Reference Tadrist, Gilet, Schlagheck and Bush2020). They have also been predicted to show anomalous two-droplet correlations (Nachbin Reference Nachbin2018; Valani, Slim & Simula Reference Valani, Slim and Simula2018). A detailed review of hydrodynamic quantum analogues of walking droplets is provided by Bush (Reference Bush2015) and Bush et al. (Reference Bush, Couder, Gilet, Milewski and Nachbin2018).
Recently, a new class of walking droplets, coined superwalkers, has been observed (Valani, Slim & Simula Reference Valani, Slim and Simula2019). These emerge when the bath is driven simultaneously at two frequencies, 
$f$ and 
$f/2$, with a relative phase difference 
${\rm \Delta} \phi$. For a commonly studied system with silicone oil of 
$20\ \text {cSt}$ viscosity, single-frequency driving at 
$f=80\ \text {Hz}$ produces walkers with radii between 0.3 and 
$0.5\ \text {mm}$, and walking speeds up to 
$15\ \text {mm}\ \text {s}^{-1}$, with speed typically increasing with size (Moláček & Bush Reference Moláček and Bush2013b; Wind-Willassen et al. Reference Wind-Willassen, Moláček, Harris and Bush2013). In the same system with two-frequency driving at 
$f=80$ and 
$f/2=40$ Hz, superwalkers can be significantly larger than walkers with radii up to 
$1.4\ \text {mm}$ and walking speeds up to 
$50\ \text {mm}\ \text {s}^{-1}$ (Valani et al. Reference Valani, Slim and Simula2019). Intriguingly, the walking speed and the vertical dynamics of superwalkers are strongly dependent on the phase difference 
${\rm \Delta} \phi$, with peak superwalking speed occurring near 
${\rm \Delta} \phi =140^{\circ }$, while near 
${\rm \Delta} \phi =45^{\circ }$ they only bounce or may even coalesce. For a fixed phase difference, smaller superwalkers typically behave very similarly to walkers, with their speed increasing with their size and impacting the surface once every two up-and-down cycles of the bath. Conversely, the speed of larger superwalkers decreases with their size. These large superwalkers appear to impact the bath twice every two up-and-down cycles of the bath and have prolonged contact with the bath, with the largest ones hardly lifting from the surface. Using sophisticated numerical simulations, Galeano-Rios, Milewski & Vanden-Broeck (Reference Galeano-Rios, Milewski and Vanden-Broeck2019) were able to replicate superwalking behaviour for a single droplet of moderate radius 
$R=0.68$ mm, and reported a good match in the superwalking speed between their simulation and the experiments of Valani et al. (Reference Valani, Slim and Simula2019). Although these two studies describe the characteristics of superwalkers, an understanding of the mechanism that enables superwalking is still lacking. In this paper, our aim is to understand this underlying mechanism by adapting the theoretical models used for walkers driven with a single frequency, to superwalkers driven with two frequencies.
Over the years, many theoretical models have been developed to describe the dynamics of a walker. These range from phenomenological stroboscopic models that only capture the horizontal dynamics to sophisticated models that resolve the vertical and horizontal dynamics and the detailed evolution of the surface waves created by the walker. Intermediate complexity models that resolve the vertical and horizontal dynamics but assume a predetermined form for the standing wave generated by the droplet on each impact have been widely used. In this latter category, Moláček & Bush (Reference Moláček and Bush2013a) modelled the vertical bouncing dynamics of the droplet using a linear spring model, inspired by the investigations of Gilet & Bush (Reference Gilet and Bush2009a) and Gilet & Bush (Reference Gilet and Bush2009b) on droplets bouncing on a soap film. They also developed a nonlinear logarithmic spring model, although it is not clear whether this model is more accurate and hence the linear spring model is often used for simplicity (Couchman, Turton & Bush Reference Couchman, Turton and Bush2019). Moláček & Bush (Reference Moláček and Bush2013b) and Couchman et al. (Reference Couchman, Turton and Bush2019) coupled these vertical spring models with a horizontal dynamics model comprising a propulsive force from the impact and a lumped drag force consisting of aerodynamic drag and momentum drag during contact. The propulsive force is the horizontal component of the normal force, which arises because the small-amplitude waves incline the free surface. The waves are modelled by a linear superposition of predetermined standing waves generated by the droplet on each impact. Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015) further refined the modelling by coupling the vertical and horizontal dynamics models of Moláček & Bush (Reference Moláček and Bush2013a,Reference Moláček and Bushb) to a quasi-potential, weakly viscous wave model for wave generation and evolution. Galeano-Rios, Milewski & Vanden-Broeck (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) developed a more complete model for the vertical dynamics of the droplet by imposing a kinematic match between the motion of the free surface and that of the impacting droplet, modelled as a solid sphere. Combining this vertical dynamics model with the free surface evolution of Milewski et al. (Reference Milewski, Galeano-Rios, Nachbin and Bush2015), Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2019) were able to obtain an accurate model for walking droplets free of any impact parametrisation. Such models give an accurate description of the system at the time scale of a single bounce but they become inefficient in calculating a droplet's horizontal trajectory over long times. Hence, stroboscopic models that average over the walker's periodic vertical motion but capture the horizontal motion have been developed to investigate the horizontal dynamics of walkers over long time scales. Many of these stroboscopic models use a predetermined form of the standing wave field (Protière, Boudaoud & Couder Reference Protière, Boudaoud and Couder2006; Oza, Rosales & Bush Reference Oza, Rosales and Bush2013; Oza et al. Reference Oza, Siéfert, Harris, Moláček and Bush2017; Arbelaiz, Oza & Bush Reference Arbelaiz, Oza and Bush2018; Couchman et al. Reference Couchman, Turton and Bush2019) while the model of Durey & Milewski (Reference Durey and Milewski2017) uses a more sophisticated wave model to accurately capture the droplet's wave field. A comprehensive review of the different walker models is given by Turton, Couchman & Bush (Reference Turton, Couchman and Bush2018).
In this study, we couple the vertical and horizontal dynamics models of Moláček & Bush (Reference Moláček and Bush2013a,Reference Moláček and Bushb) along with a new model for the wave field of a superwalker to understand and rationalise superwalking. In § 2, we provide a summary of the theoretical model, explore the wave field for two-frequency driving and describe the nomenclature we use to describe the bouncing modes. In § 3, we show that adding the second driving frequency with an appropriate phase difference raises every second peak and lowers the intermediate peaks of the bath's motion, and that this allows larger droplets to leap over the intermediate peaks, thereby enabling superwalking. We also show the importance of the phase difference 
${\rm \Delta} \phi$ in the emergence of superwalking and compare the results from simulations with the experiments of Valani et al. (Reference Valani, Slim and Simula2019). We discuss and conclude the study in § 4.
2. Theoretical model
Consider a droplet of mass 
$m$ and radius 
$R$ bouncing on a bath of liquid of density 
$\rho$, kinematic viscosity 
$\nu$ and surface tension 
$\sigma$. The bath is vibrating vertically with acceleration 
$\gamma (t)=\varGamma _f{g}\sin (\varOmega t)+\varGamma _{f/2}{g}\sin (\varOmega t/2+{\rm \Delta} \phi )$. Here, 
$\varOmega =2{\rm \pi} f$ is the angular frequency, 
$\varGamma _f$ and 
$\varGamma _{f/2}$ are the acceleration amplitudes of the two frequencies relative to gravity 
$g$ and 
${\rm \Delta} \phi$ is the relative phase difference between them. This configuration is shown schematically in figure 1. We describe our system in the oscillating frame of the bath by horizontal coordinates 
$\boldsymbol {x} = (x,y)$ and vertical coordinate 
$z$, with 
$z=0$ chosen to coincide with the undeformed surface of the bath. In this frame, the centre of mass of the droplet is located at a horizontal position 
$\boldsymbol {x}_d$ and the south pole of the droplet at a vertical position 
$z_d$ such that 
$z_d=0$ would represent initiation of contact with the undeformed surface of the bath. The free surface elevation of the liquid filling the bath is at 
$z=h(\boldsymbol {x},t)$.
(a) Schematic of the system consisting of a bath of liquid vibrated vertically with acceleration 
$\gamma (t)$ and a droplet of the same liquid walking horizontally at velocity 
$\dot {\boldsymbol {x}}_d$ and located at vertical position 
$z_d$ relative to the free surface of the liquid at rest. Panel (b) shows the vertical and horizontal forces acting on the droplet in the comoving frame of the bath. In the vertical direction, the droplet experiences an effective gravity, 
$-m[{g}+\gamma (t)]$, and an upward normal force, 
$F_N(t)$, during contact. In the horizontal direction the droplet experiences a propulsive force, 
$-F_N(t) \nabla h(\boldsymbol {x}_d,t)$, during contact due to the slope of the wave field and a lumped drag force composed of momentum loss during contact, 
$-D_{mom}\dot {\boldsymbol {x}}_d$ and air drag, 
$-D_{air}\dot {\boldsymbol {x}}_d$.
2.1. Vertical dynamics
We simulate the vertical droplet dynamics using the linear spring model of Moláček & Bush (Reference Moláček and Bush2013a) adapted for two-frequency driving. Using this model, the vertical equation of motion of the droplet in the comoving frame of the bath is given by
\begin{equation} m\ddot{z}_d=-m[{g}+\gamma(t)]+{F_N}(t), \end{equation}where the first term on the right-hand side is the effective gravitational force on the droplet in the bath's frame of reference. The second term on the right-hand side is the normal force imparted to the droplet during contact with the liquid surface. This normal force is calculated by modelling interaction with the bath as a linear spring and damper according to,
\begin{equation} {F_N}(t)=H(-\bar{z}_d)\ max\left(-k_s\bar{z}_d-b\dot{\bar{z}}_d,0\right). \end{equation}
Here, 
$H$ is the Heaviside step function and 
$\bar {z}_d=z_d-h(\boldsymbol {x}_d,t)$ is the height of the droplet's south pole above the free surface of the bath (from here on simply referred to as the height of the droplet). The maximum condition in (2.2) ensures a non-negative reaction force on the droplet during contact. The constants 
$k_s$ and 
$b$ are the spring constant and damping force coefficient, respectively. We discuss appropriate values for these constants in § 2.4.
2.2. Horizontal dynamics
To describe the horizontal dynamics of the walking droplet, we use the model of Moláček & Bush (Reference Moláček and Bush2013b) for which the horizontal equation of motion is
\begin{equation} m\ddot{\boldsymbol{x}}_d=-\left[D_{mom}(t)+D_{air}\right] \dot{\boldsymbol{x}}_d-{F_N}(t) \boldsymbol{\nabla} h(\boldsymbol{x}_d,t). \end{equation}
The term in square brackets on the right-hand side is the total instantaneous drag force, composed of momentum loss during contact, 
$D_{mom}(t)=C\sqrt {\rho R/\sigma }{F_N}(t)$, and an air drag of the form 
$D_{air}=6{\rm \pi} R \mu _a$. Here, 
$\mu _a$ is the dynamic viscosity of air and 
$C$ is the contact drag coefficient, an adjustable parameter. We discuss appropriate values for this coefficient in § 2.4. The final term on the right-hand side is the horizontal component of the contact force arising from the small slope, 
$|\boldsymbol {\nabla } h(\boldsymbol {x}_d,t)|\ll 1$, of the underlying wave field during contact.
2.3. Wave field
The free surface 
$z=h(\boldsymbol {x},t)$ is calculated as the linear superposition of all the individual waves generated by the droplet on its previous bounces,
\begin{equation} h(\boldsymbol{x},t)=\sum_{n} h_n(\boldsymbol{x},t), \end{equation}
where 
$h_n(\boldsymbol {x},t)$ is the wave field generated by the 
$n$th bounce at location 
$\boldsymbol {x}_n$ and time 
$t_n$.
The individual waves generated by the droplet on each bounce are localised, damped standing Faraday waves. Various different models of the waveform have been developed to describe a single impact of a walker (Eddi et al. Reference Eddi, Sultan, Moukhtar, Fort, Rossi and Couder2011; Moláček & Bush Reference Moláček and Bush2013b; Milewski et al. Reference Milewski, Galeano-Rios, Nachbin and Bush2015; Tadrist et al. Reference Tadrist, Shim, Gilet and Schlagheck2018). One of the most commonly used wave models is that of Moláček & Bush (Reference Moláček and Bush2013b), given by
\begin{equation} h^{(M)}_n(\boldsymbol{x},t)=\frac{A^{(M)}}{\sqrt{t-t_n}}\cos \left(\frac{\varOmega t}{2}\right)\textrm{J}_0(k_F|\boldsymbol{x}-\boldsymbol{x}_n|) \exp{\left[-\frac{t-t_n}{T_F {Me}^{(M)}}\right]}, \end{equation}
where 
$k_F$ is the Faraday wavenumber, 
$T_F=4{\rm \pi} /\varOmega$ is the Faraday period, and 
${Me}^{(M)}=T_d/T_F(1-\varGamma _f/\varGamma _F)$ is the memory parameter that determines the proximity to the Faraday threshold. In this latter expression, 
$T_d=1/\nu _e k_F^2$ is the time constant for wave decay, 
$\nu _e$ is the effective kinematic viscosity and 
$\varGamma _F$ is the dimensionless acceleration amplitude at the Faraday threshold for single-frequency driving at frequency 
$f$. This model describes a wave with the shape of a Bessel function of the first kind, 
$\textrm {J}_0$, that oscillates at the subharmonic frequency 
$f/2$ and decays exponentially in time with a decay constant inversely proportional to the memory. The location and instant of the droplet's impact are approximated respectively by,
\begin{equation} \boldsymbol{x}_n=\int_{t_n^i}^{t_n^c}\boldsymbol{x}_d(t')F_N(t')\,\text{d}t' \Bigg{/}\int_{t_n^i}^{t_n^c}F_N(t')\,\text{d}t'\quad \text{and}\quad t_n=\int_{t_n^i}^{t_n^c}t'F_N(t')\,\text{d}t'\Bigg{/}\int_{t_n^i}^{t_n^c}F_N(t')\,\text{d}t', \end{equation}
where 
$t_n^i$ and 
$t_n^c$ are the times of initiation and completion of the 
$n$th impact. The equation for the wave amplitude coefficient 
$A^{(M)}$ is
\begin{equation} A^{(M)}=\sqrt{\frac{2\nu_e}{\rm \pi}}\frac{k_F^3}{3\sigma k_F^ 2 + \rho {g}} \int_{t_n^i}^{t_n^c}\sin\left(\frac{\varOmega t'}{2}\right)F_N(t') \,\text{d}t'. \end{equation} A detailed theoretical study of the wave field generated by a single bounce of a walker was undertaken by Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018). They derived the following improved waveform for the wave generated by an instantaneous impact of a walker of force strength 
$F_0$ at location 
$\boldsymbol {x}_n$ and time 
$t_n$,
\begin{align} h^{(T)}_n(\boldsymbol{x},t)=\frac{A_0^{(T)}}{\sqrt{t-t_n}}\cos \left(\frac{\varOmega t}{2}+\theta^{+}_F\right)\textrm{J}_0(k_F|\boldsymbol{x}-\boldsymbol{x}_n|) \exp{\left[-\frac{(t-t_n)}{T_F {Me}^{(T)}}-\frac{T_F|\boldsymbol{x}-\boldsymbol{x}_n|^2}{8{\rm \pi} D(t-t_n)}\right]}. \end{align} In this expression, the memory parameter is given by 
${Me}^{(T)}=-1/2{\rm \pi} \delta ^{+}_F$ with 
$\delta ^{+}_F$ the dimensionless decay rate of the longest-lived Faraday wave. This improved form of the wave field has two new additions: (i) the phase shift 
$\theta ^{+}_F$ of the Faraday waves relative to the driving signal and (ii) an exponential spatial decay with diffusive spreading (with a diffusion coefficient 
$D$). Note that similar additions can also be obtained following the derivation of Moláček & Bush (Reference Moláček and Bush2013b) by including higher-order terms in their decay rate expansions. The amplitude coefficient 
$A_0^{(T)}$ takes the form,
\begin{equation} A_0^{(T)}=\sqrt{\frac{2{\rm \pi}}{\varOmega^5 D}}\frac{2k_F^2}{{\rm \pi}\rho}B^{+}_F(t_n)F_0, \end{equation}
where 
$B^{+}_F(t_n)$ is a function that prescribes the amplitude based on the instant of impact 
$t_n$ and is given by
\begin{equation} B^{+}_F(t_n) =\frac{-2\cos(\varOmega t_n/2+\theta_{F}^{-})}{(\delta^{+}_F-\delta^{-}_F) [\cos(\varOmega t_n+\theta^{+}_{F}+\theta^{-}_{F})+\cos(\theta^{+}_{F}-\theta^{-}_{F})]-2\sin(\theta^{+}_{F}-\theta^{-}_{F})}, \end{equation}
where 
$\theta ^{-}_{F}$ and 
$\delta ^{-}_{F}$ are the phase shift and decay rate of a companion, short-lived Faraday wave. The reader is referred to Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) for further details on these parameters. To extend this model to a finite contact time, we follow the suggestion in Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) of using Duhamel's principle and the approach used in appendix A.4 of Moláček & Bush (Reference Moláček and Bush2013b), and integrate the impulse response with a time varying contact force 
$F_N(t)$. This results in replacing the amplitude coefficient 
$A_0^{(T)}$ in (2.8) by
\begin{equation} A^{(T)}=\sqrt{\frac{2{\rm \pi}}{\varOmega^3 D}}\frac{k_F^2}{{\rm \pi} \rho}\int_{t_n^i}^{t_n^c}B^{+}_F(t') F_N(t') \,\text{d}t', \end{equation}
and replacing the initial impact location 
$\boldsymbol {x}_n$ and time 
$t_n$ with their respective weighted averages as defined in (2.6a,b). Moreover, in (2.11), the amplitude prescribing function 
$B_F^{+}(t')$ is identical to the one presented in (2.10).
To model the wave field generated by a superwalker with two-frequency driving, we use a similar approach to that of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018). A detailed derivation is presented in appendix A. The final form of the wave field for the case of the bath being driven at 
$f=80$ and 
$f/2=40$ Hz is
\begin{align} & h_n^{({SW})}(\boldsymbol{x},t) = A_{40}\frac{\cos(\varOmega t/2+\theta_{F40}^{+})}{\sqrt{t-t_n}} \text{J}_0(k_{F40}|\boldsymbol{x}-{\boldsymbol{x}}_n|) \exp{\left[ -\frac{t-{t}_n}{T_F {Me}_{40}} - \frac{T_F |\boldsymbol{x}-{\boldsymbol{x}}_n|^2}{8 {\rm \pi}D_{40} (t-{t}_n)} \right]} \nonumber\\ &\quad + A_{20}\frac{\cos(\varOmega t /4 + \theta_{F20}^{+})}{\sqrt{t-t_n}} \text{J}_0(k_{F20}|\boldsymbol{x}-{\boldsymbol{x}}_n|)\exp{\left[ -\frac{t-{t}_n}{T_F {Me}_{20}} - \frac{T_F |\boldsymbol{x}-{\boldsymbol{x}}_n|^2}{8 {\rm \pi}D_{20} (t-t_n)} \right]}, \end{align}
where the impact location 
$\boldsymbol {x}_n$ and the time of impact 
$t_n$ are calculated using (2.6a,b), and the amplitude coefficients are given by
\begin{gather} A_{40}=\sqrt{\frac{2{\rm \pi}}{\varOmega^3 D_{40}}} \frac{k_{F40}^2}{{\rm \pi} \rho}\int_{t_n^i}^{t_n^c}B^{+}_{F40}(t') F_N(t')\,\text{d}t', \end{gather}
\begin{gather}A_{20}=\sqrt{\frac{2{\rm \pi}}{\varOmega^3 D_{20}}}\frac{k_{F20}^2}{{\rm \pi} \rho} \int_{t_n^i}^{t_n^c}B^{+}_{F20}(t') F_N(t')\,\text{d}t'. \end{gather}
The interpretation of (2.12) is that a droplet bouncing under the prescribed two-frequency driving excites two dominant waves at wavenumbers 
$k_{F40}$ and 
$k_{F20}$, corresponding to frequencies of 
$40$ and 
$20$ Hz. These waves decay in time at rates 
$\text {Re}(\delta ^{+}_{F40})$ and 
$\text {Re}(\delta ^{+}_{F20})$, which determine the corresponding memory parameters 
${Me}_{40}=-1/2{\rm \pi} \text {Re}(\delta ^{+}_{F40})$ and 
${Me}_{20}=-1/2{\rm \pi} \text {Re}(\delta ^{+}_{F20})$. Here, 
$\text {Re}(\cdot )$ denotes the real part of the complex argument. The waves also spread diffusively with diffusion coefficients 
$D_{40}$ and 
$D_{20}$, and have phase shifts 
$\theta _{F40}^{+}$ and 
$\theta _{F20}^{+}$ with respect to the driving signal. We refer the reader to appendix A for explicit equations for these parameters.
Comparing the superwalker wave field in (2.12) to that of a walker in (2.8) leads to two key observations: (i) Both models have a wave at frequency 
$f/2=40$ Hz. We note that Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) derived (2.8) by considering a cosine form of driving while we have considered a sine form of driving to be consistent with experiments of Valani et al. (Reference Valani, Slim and Simula2019). This results in a constant shift of 
${\rm \pi} /4$ in the phase shift 
$\theta ^{+}_F$ in (2.8) which has been taken into account when comparing results. (ii) An additional wave of frequency 
$f/4=20$ Hz appears in the wave field of a superwalker. However, in the region of 
$(\varGamma _{80},\varGamma _{40})$ parameter space where superwalking is realised, typically the amplitude of the 40 Hz wave, 
$A_{40}$, is 
$4$ to 
$8$ times that of the 20 Hz wave, 
$A_{20}$. Thus in general, our new two-frequency wave model is not appreciably different from the single-frequency model of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018). This is illustrated further in figure 2 where the wave fields predicted using the models of Moláček & Bush (Reference Moláček and Bush2013b) in (2.5), Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) in (2.8) and the superwalker wave field in (2.12) are shown for an instantaneous impact at time 
$0.22T_F$, corresponding to the typical impact phase for superwalkers, with an appreciable 
$\varGamma _{f/2}$ component. The waves from our new two-frequency model (2.12) and the single-frequency (Tadrist et al. Reference Tadrist, Shim, Gilet and Schlagheck2018) model (2.8) are quantitatively similar (figures 2a and 2c–e), except near times when the overall wave field is quite flat and is changing rapidly (figure 2b). The comparison with the single-frequency (Moláček & Bush Reference Moláček and Bush2013b) model appears poorer, however, the difference is primarily in the far field and arises from the absence of diffusive spatial decay in this model. In the near-field region of primary interest for walking, all three models are quantitatively similar with a maximum relative error of approximately 
$20\,\%$, as shown in figures 2(e) and 2(f). Moreover, as shown in figure 2(f), the relative height difference at the impact location between the waves of the Moláček & Bush (Reference Moláček and Bush2013b) model and the Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) model is sinusoidal due to the added phase shift of 
$\theta _F^{+}\approx -4^{\circ }$ in the Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) model for the chosen parameters in figure 2, and the relative height difference at the impact location between the superwalker wave and that of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) reveals the added 
$20$ Hz wave in the superwalker wave field. Overall these results suggest that the wave fields observed for two-frequency and single-frequency driving remain very similar, an observation that was also made qualitatively by Valani et al. (Reference Valani, Slim and Simula2019). In § 3, we present results using our new two-frequency model, but note that results using either the Moláček & Bush (Reference Moláček and Bush2013b) model or the Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) model are comparable; we provide details in appendix B.
Comparison of the wave fields generated by an instantaneous impact at 
$x=0$ and at time 
$t_i=0.22T_F$ for typical superwalker parameter values. The wave fields from the Moláček & Bush (Reference Moláček and Bush2013b) model (green dashed-dotted curve), Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) model (maroon dotted curve) and the superwalker model of this work (blue solid curve) are shown at times (a) 
$t_1=t_i+0.23T_F$, (b) 
$t_2=t_i+0.57T_F$, (c) 
$t_3=t_i+0.76T_F$ and (d) 
$t_4=t_i+1.00T_F$. The evolution of the absolute wave height 
$h$ at 
$x=0$ from an impact at 
$t_i$ (vertical red dashed line) is shown in (e) and the relative wave height 
${\rm \Delta} h/h_*^{(T)}$ with respect to the Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) model is shown in (f). Here, 
$h_*^{(T)}$ is the wave field from Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) model in (2.8) excluding the cosine term to avoid singularities in 
${\rm \Delta} h/h_*^{(T)}$, and 
${\rm \Delta} h = h^{(SW)} - h^{(T)}$ or 
$h^{(M)} - h^{(T)}$. The parameters are 
$\varGamma _{80}=3.8$, 
$\varGamma _{40}=0.6$ and 
${\rm \Delta} \phi =130^{\circ }$.
2.4. Numerical method and parameter values
We solved (2.1) and (2.3) using the leap-frog method (Sprott Reference Sprott2003), a modified version of the Euler method where the new horizontal and vertical positions of the droplet are calculated using the old velocities and then the new velocities are calculated using the new positions. To increase computational speed, we only stored the waves generated by the 
$100$ most recent bounces of the droplet and discard the earlier ones, which have typically decayed to below 
$10^{-5}$ of their initial amplitude. The simulations were initialised with 
$\boldsymbol {x}_d=(0,0)\ \text {mm}$, 
$\dot {\boldsymbol {x}}_d=(1,0)$, 
$\dot {z}_d=0\ \textrm {mm}\ \textrm {s}^{-1}$ and six different equally spaced vertical positions 
$z_d= (0,2,4,6,8,10)R$. Multiple initial conditions were used so that different modes existing at the same parameter values are likely to be captured.
The physical parameters were fixed to match the experiments of Valani et al. (Reference Valani, Slim and Simula2019): 
$\rho =950\ \text {kg}\ \text {m}^{-3}$, 
$\nu =20$ cSt, 
$\sigma =20.6\ \text {mN}\ \text {m}^{-1}$ and 
$f=80$ Hz. There are three adjustable parameters in the model: the spring constant of the bath 
$k_s$, the damping coefficient of the bath 
$b$ and the dimensionless contact drag coefficient 
$C$. The dimensionless parameters corresponding to 
$k_s$ and 
$b$ are 
$K=k_s/m\omega _d^2$ and 
$B=b/m\omega _d$, where 
$\omega _d=\sqrt {\sigma /\rho R^3}$ is the droplet's characteristic internal vibration frequency (Moláček & Bush Reference Moláček and Bush2013a). For walkers, appropriate values were determined by fitting to experimental data (Moláček & Bush Reference Moláček and Bush2013a,Reference Moláček and Bushb) and typical values are 
$K=0.59$ and 
$B=0.48$ (Couchman et al. Reference Couchman, Turton and Bush2019), and 
$C=0.17$ (Moláček & Bush Reference Moláček and Bush2013b). For superwalkers, we also set 
$C=0.17$, but adjust 
$K$ and 
$B$ to best fit the available experimental data. The details of this fit are described in appendix C. We use both constant values of 
$K=0.70$ and 
$B=0.60$, as well as allowing the parameter 
$K$ to vary with droplet radius 
$R$ according to
\begin{equation} K=1.06\sqrt{{Bo}}+0.37 \end{equation}
with a fixed 
$B=0.60$, where 
${Bo}=\rho g R^2/\sigma$ is the Bond number of the droplet. We refer the reader to appendix C for more details on how this relationship was obtained. We note that these values give a good match with the experiments of Valani et al. (Reference Valani, Slim and Simula2019); however, the qualitative behaviour of the results remains unchanged for a range of 
$K$ and 
$B$ values.
2.5. Description of vertical bouncing modes
The vertical bouncing dynamics is crucially important for the existence and characteristics of superwalkers. Here, we provide a description of the nomenclature we use to distinguish the different vertical bouncing modes.
We follow Valani et al. (Reference Valani, Slim and Simula2019) and use the notation 
$(l,m,n)$ to indicate that the droplet impacts the surface 
$n$ times during 
$m$ oscillation periods of the bath at frequency 
$f$, which equals 
$l$ oscillation periods of the bath at frequency 
$f/2$. For single-frequency driving, the index 
$l$ is dropped. For normal walking droplets, one of the most commonly observed mode is 
$(2,1)$, with the droplets leaping over every second peak in the bath's motion. After Moláček & Bush (Reference Moláček and Bush2013a), we distinguish two different styles of 
$(2,1)$ walking and corresponding 
$(1,2,1)$ superwalking, with a high-bouncing, short-contact mode denoted by 
$(2,1)^{H}$ and 
$(1,2,1)^{H}$, and a low-bouncing, long-contact mode denoted by 
$(2,1)^{L}$ and 
$(1,2,1)^{L}$. Droplets that have two peaks in the normal force during contact are classified as 
$(1,2,1)^{{L}}$ while those that have only one peak are classified as 
$(1,2,1)^{{H}}$ (Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2019). Another commonly observed mode is 
$(2,2)$ and corresponding 
$(1,2,2)$, in which the droplets no longer are able to leap over intermediate peaks and contact the bath twice, typically a high bounce and a low bounce, every two up-and-down cycles of the bath. Note that experimentally it is difficult to distinguish between a 
$(2,1)^{L}$ and 
$(2,2)$ mode (see figures 7 and 8 of Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2019). A less commonly observed mode is 
$(4,2)$ and corresponding 
$(2,4,2)$, in which the droplets leap over every second peak, but each bounce in a pair has different amplitude. Finally, bouncing modes with no discernible periodicity or those with periodic contact but aperiodic modulation of peak bouncing heights are referred to as chaotic modes.
3. Emergence of superwalking
To illustrate the emergence of superwalking and its relationship with normal walking, we begin by describing the dynamics of a relatively small normal walker with the bath driven at a single frequency of 
$f=80$ Hz and acceleration amplitude 
$\varGamma _{80}=3.8$ (compared to a Faraday threshold 
$\varGamma _{F80}=4.15$). This results in a normal walker that is bouncing in a 
$(2,1)^{{H}}$ mode (see figure 3a). The 
$(2,1)$ bouncing mode is crucial for walking as the droplet is bouncing at the same frequency as the frequency of the subharmonic Faraday waves that emerge beyond the Faraday instability threshold. Thus the droplet's bouncing is in resonance with the damped Faraday waves it generates and with which it interacts. For slightly larger droplets (see figure 3b), the same 
$(2,1)^{{H}}$ bouncing mode is maintained but the height of the bounces are reduced, while for larger droplets still, the bounces reduce in height to such an extent that the droplet can no longer leap over the second peak in the bath's motion. For the chosen parameters, this results in the droplet transitioning to a chaotic bouncing mode and no longer walking (figure 3c).
Emergence of a superwalker. (a–c) Vertical dynamics of a walker of radius (a) 
$R=0.36$ mm and (b) 
$R=0.40$ mm bouncing in a 
$(2,1)^{H}$ mode, and a bigger non-walking droplet of radius (c) 
$R=0.54$ mm bouncing in a chaotic mode. Here, the bath is driven at a single frequency of 
$f=80$ Hz with acceleration amplitude 
$\varGamma _{80}=3.8$. (d) Vertical dynamics of a superwalker of radius 
$R=0.54$ mm bouncing in a 
$(1,2,1)^{{H}}$ mode. Here the bath is driven at 
$f=80$ and 
$f/2=40$ Hz with phase difference 
${\rm \Delta} \phi =130^{\circ }$ and acceleration amplitudes 
$\varGamma _{80}=3.8$ and 
$\varGamma _{40}=0.6$. In (a–d), the solid black curves indicate the bath motion, 
$\mathcal {B}(t)=-(\varGamma _f{g}/\varOmega ^2)\sin (\varOmega t)-(4\varGamma _{f/2}{g}/\varOmega ^2)\sin (\varOmega t/2+{\rm \Delta} \phi )$, the coloured curves represent the motion of the south pole of the droplet, 
$z_d(t)+\mathcal {B}(t)$, and the filled blue regions illustrate the motion of the liquid surface, 
$h(\boldsymbol {x}_d,t)+\mathcal {B}(t)$, all in the laboratory frame. The grey regions indicate times at which the droplet is in contact with the bath. Panel (e) shows a schematic of the speed-size characteristics for the droplets in (a–d). Here, the values of the parameters 
$K$ and 
$B$ are both fixed to 
$0.70$ and 
$0.60$ respectively.
In contrast, figure 3(d) shows the vertical dynamics of the same-sized droplet as in figure 3(c) with the addition of the subharmonic frequency 
$f/2=40$ Hz and amplitude 
$\varGamma _{40}=0.6$ (compared to a Faraday threshold 
$\varGamma _{F40}=1.22$) at a phase difference of 
${\rm \Delta} \phi =130^\circ$. This additional subharmonic driving raises every second peak and lowers the intermediate peaks in both the bath's and the waves’ motion. This allows the bigger droplet to clear every second peak in the bath's motion and settle in a 
$(1,2,1)^{{H}}$ bouncing mode, effectively identical to the 
$(2,1)^{{H}}$ mode for a walker, and results in the emergence of a superwalker. This jump from a walker to a superwalker is shown schematically on the speed-size curve in figure 3(e).
3.1. Importance of the phase difference between the two driving frequencies
The phase difference between the two driving frequencies controls the relative height of the two peaks in one full cycle of the periodic bath motion, equivalently two up-and-down cycles, and it is therefore a crucial parameter for the emergence of superwalking. Figure 4(a) shows the walking speed 
$u$ as a function of the phase difference 
${\rm \Delta} \phi$ for a fixed-sized droplet that is too large to walk with single-frequency driving (the largest droplet shown in figure 3). The different vertical modes at different 
${\rm \Delta} \phi$ are shown in figure 4(b). Depending on the phase difference, the droplet either bounces without walking or it superwalks. In the bouncing regime, 
$20^{\circ }\lesssim {\rm \Delta} \phi \lesssim 90^{\circ }$, the droplet's vertical dynamics appears chaotic. This can be attributed to the height difference 
${\rm \Delta} \mathcal {B}$ between successive peaks in the bath's motion being small (see dashed curve in figure 4a) and hence the droplet behaves similarly to the single-frequency case (see figure 3c). Conversely, regions of high superwalking speed are associated with a large height difference 
${\rm \Delta} \mathcal {B}$ between the two peaks in the bath's motion and a droplet can easily settle in a 
$(1,2,1)$ bouncing mode.
Effect of phase difference on superwalking behaviour. (a) Walking speed 
$u$ as a function of the phase difference 
${\rm \Delta} \phi$ for a superwalker of radius 
$R=0.54$ mm with 
$\varGamma _{80}=3.8$ and 
$\varGamma _{40}=0.6$. The solid curve represents results from numerical simulations with colours indicating different bouncing modes: 
$(1,2,1)^{{L}}$ in green, 
$(1,2,1)^{{H}}$ in blue and chaotic in purple. The experimental results of Valani et al. (Reference Valani, Slim and Simula2019) are shown by points, with the style of marker indicating the bouncing modes: 
$(1,2,2)$ are red circles 
$\bullet$, 
$(1,2,1)^{{H}}$ are blue triangles 
$\blacktriangle$, transition between a 
$(1,2,1)^{{H}}$ and a 
$(1,2,2)$ mode are grey squares 
$\blacksquare$, and chaotic are purple asterisks 
$*$. The dashed curve indicates the height difference 
${\rm \Delta} \mathcal {B}$ between consecutive peaks in one period of the bath motion. The data to the right of the vertical dotted line are repeated. Panel (b) shows bouncing modes obtained for different values of 
${\rm \Delta} \phi$ from (a). In this panel, the grey regions indicate times at which the droplet is in contact with the bath. The parameters 
$K$ and 
$B$ are fixed to 
$0.70$ and 
$0.60$ respectively.
The predicted speeds from the numerical simulations agree well with experiments. The chaotic mode in the bouncing regime and the 
$(1,2,1)^{{H}}$ bouncing mode in the superwalking regime are also observed at parameter values comparable to those in experiments. The numerically observed 
$(1,2,1)^{{L}}$ superwalkers were not reported in experiments, instead 
$(1,2,2)$ modes were observed at the corresponding parameter values. However, as noted in § 2.5, it is difficult to distinguish between a 
$(1,2,1)^{L}$ and a 
$(1,2,2)$ mode experimentally. Hence, it is not clear whether all the 
$(1,2,2)$ superwalkers reported by Valani et al. (Reference Valani, Slim and Simula2019) are truly 
$(1,2,2)$ superwalkers or if some may in fact be 
$(1,2,1)^{{L}}$ superwalkers.
3.2. Speed-size characteristics of superwalking droplets
In the size range for which walkers exist, their walking speed typically increases with their size (Moláček & Bush Reference Moláček and Bush2013b). For superwalkers, two trends were observed in experiments: an ascending branch for smaller superwalkers where the speed increases with size, followed by a descending branch for larger superwalkers where the speed decreases with size (Valani et al. Reference Valani, Slim and Simula2019). Figure 5 shows the speed-size characteristics of simulated superwalkers at 
$\varGamma _{80}=3.8$ and 
${\rm \Delta} \phi =130^{\circ }$ for a range of 
$\varGamma _{40}$ values.
Speed of a superwalker as a function of its size at fixed 
$\varGamma _{80}=3.8$ and 
${\rm \Delta} \phi =130^{\circ }$. (a) Comparison of the speed-size characteristics of a droplet from numerical simulations (solid curves) with experimental results of Valani et al. (Reference Valani, Slim and Simula2019) (black circles with empty circles indicating coalescence) and the stroboscopic model of Oza et al. (Reference Oza, Rosales and Bush2013) (dashed curve) at 
$\varGamma _{40}=0.6$. For the stroboscopic model, we set the adjustable parameters of the impact phase and the non-dimensional drag coefficient to 
$\sin (\varPhi ) = 0.2$ and 
$C = 0.17$ respectively, while the other parameters were chosen to match the experiments of Valani et al. (Reference Valani, Slim and Simula2019). The black horizontal bars indicate where different bouncing modes in experiments were observed. Panels (c), (d) and (e) show the speed-size characteristics at 
$\varGamma _{40}=0$, 
$\varGamma _{40}=0.3$ and 
$\varGamma _{40}=1$ respectively. In each panel the grey curve is for fixed parameter values of 
$K=0.70$ and 
$B=0.60$, and multicoloured curve represents when 
$K$ is varied as a linear function of the droplet radius 
$R$ according to (2.14) with a fixed 
$B=0.60$. The colours on this curve represent a chaotic bouncing mode in purple, 
$(2,4,2)$ mode in yellow, 
$(1,2,1)^{{H}}$ mode in blue, 
$(1,2,1)^{{L}}$ mode in green and 
$(1,2,2)$ mode in red. Termination of the solid curves indicate coalescence. The typical bouncing modes from (a) at different droplet radii are shown in (b). In this panel, the grey regions indicate times at which the droplet is in contact with the bath.
We begin by focusing on the comparison for the ascending branch. Simulated superwalking speeds for different droplet radii for constant 
$K=0.70$ and 
$B=0.60$ (grey curves) as used in the simulations shown in figures 3 and 4, and 
$K$ linearly increasing function of droplet radius as in (2.14) with a fixed 
$B=0.60$ (coloured curves) are shown. We refer the reader to appendix C for details on this linear relationship. Both the superwalking speed and the bouncing mode are captured well for both combinations for small- to moderate-sized superwalkers, and this is generally true for a broad range of 
$K$ and 
$B$ values (see appendix C). By allowing 
$K$ to vary linearly with the droplet radius 
$R$ (coloured curve), we obtain a better fit for droplets on the ascending branch at relatively high 
$\varGamma _{40}$ values (see figure 5e). Focusing on the vertical dynamics for this fit when 
$\varGamma _{40}=0.6$ (see figure 5a), we find that superwalkers on this branch universally impact the bath once every two up and down cycles of the bath's motion. For the smallest superwalkers, the amplitude of the bounces is chaotic. As the droplet size increases, there is a transition to a 
$(2,4,2)$ mode in a narrow region near 
$R=0.51$ mm. Beyond this, the droplets bounce in a 
$(1,2,1)^{{H}}$ mode (blue) for the remainder of the ascending branch. This agrees well with the experimental results of Valani et al. (Reference Valani, Slim and Simula2019) where chaotic and 
$(1,2,1)^{{H}}$ bouncing modes were also reported on the ascending branch.
Simulations of larger droplets that lie on the descending branch in experiments reveal that the model is unable to capture the larger superwalkers. We have explored different constant values of 
$K$ and 
$B$ as well as varying 
$K$ and 
$B$ as a function of 
$R$ but were unable to obtain a better fit to the experimental superwalking speeds on this branch than the relatively poor fits shown in figure 5. However, we note that the bouncing modes predicted from simulations on the descending branch are comparable with experimental observations. For the curves shown, the superwalkers on the descending branch bounce typically bounces in a 
$(1,2,1)^{{L}}$ mode. Although only the 
$(1,2,2)$ mode was reported in experiments, as previously mentioned, 
$(1,2,1)^{{L}}$ and 
$(1,2,2)$ are similar and have been difficult to distinguish in experiments.
3.3. Superwalking behaviour in the 
$(\varGamma _{80},\varGamma _{40})$ parameter space
By fixing the phase difference 
${\rm \Delta} \phi$ and the droplet radius 
$R$, we can explore the vertical and horizontal dynamics of a droplet in the parameter space formed by the two acceleration amplitudes 
$\varGamma _{80}$ and 
$\varGamma _{40}$. We choose a droplet radius of 
$R=0.60$ mm and a phase difference of 
${\rm \Delta} \phi =130^{\circ }$ to compare our results with the experimental results of Valani et al. (Reference Valani, Slim and Simula2019). Figure 6(a) shows the region of parameter space where bouncing (lighter shades) and superwalking (darker shades) are observed as well as the bouncing modes (different colours) observed in those regions. Regions of bouncing (empty circles) and superwalking (filled circles) that were identified in the experiments of Valani et al. (Reference Valani, Slim and Simula2019) are also shown. We find an excellent agreement in the transition boundary from bouncing to superwalking. Moreover, we identify that the superwalking region is dominated by the 
$(1,2,1)$ bouncing mode with a 
$(1,2,1)^{{L}}$ mode when 
$\varGamma _{40}$ is small and a 
$(1,2,1)^{{H}}$ mode when 
$\varGamma _{40}$ is large. In contrast, the bouncing mode is nearly independent of 
$\varGamma _{80}$ at a fixed 
$\varGamma _{40}$ except at relatively high 
$\varGamma _{80}$ values.
Superwalking behaviour in the 
$(\varGamma _{80},\varGamma _{40})$ parameter space. (a) Bouncing modes shown as different colours for a droplet of radius 
$R=0.60$ mm in the 
$(\varGamma _{80},\varGamma _{40})$ parameter space with multiple colours at the same 
$(\varGamma _{80},\varGamma _{40})$ values indicating multiple bouncing modes that were observed at the same 
$(\varGamma _{80},\varGamma _{40})$ values. The lighter shade of each colour indicates bouncing and the darker shade is where superwalking is observed with walking speed indicated by dotted constant speed contours in mm/s. The markers indicate bouncing (empty circles) and superwalking (filled circles) for a droplet of radius 
$R=(0.60\pm 0.05)$ mm from the experiments of Valani et al. (Reference Valani, Slim and Simula2019). (b) A vertical slice of the parameter space in (a) (solid line) showing walking speed 
$u$ as a function of 
$\varGamma _{40}$ at a fixed 
$\varGamma _{80}=3.8$. The solid curve is the result from simulations with colours indicating bouncing modes and the filled black markers are the experimental walking speeds from Valani et al. (Reference Valani, Slim and Simula2019) for a droplet of radius 
$R=(0.63\pm 0.03)$ mm. The grey shaded region indicates the jump in walking speed for this droplet when 
$\varGamma _{40}$ is appreciable. The different bouncing modes at different 
$\varGamma _{40}$ values are shown in (c) with the grey regions in this panel indicating contact with the bath. The phase difference is fixed to 
${\rm \Delta} \phi =130^{\circ }$. The parameters 
$K$ and 
$B$ are fixed to 
$0.70$ and 
$0.60$ respectively.
To understand how the superwalking speed 
$u$ changes as a function of 
$\varGamma _{40}$, we show a slice of the 
$(\varGamma _{80},\varGamma _{40})$ parameter space at 
$\varGamma _{80}=3.8$ in figure 6(b) with corresponding bouncing modes in figure 6(c). We find that the walking speed is initially zero for all 
$\varGamma _{40}\lesssim 0.3$ before increasing rapidly with 
$\varGamma _{40}$ to a peak value near 
$\varGamma _{40}=0.7$ and then marginally decreasing. This illustrates the rather abrupt rise in walking speed that occurs once the asymmetry between the heights of succeeding peaks in the bath's and waves’ motion is sufficient. Comparison of this numerically simulated walking speed 
$u$ versus acceleration amplitude 
$\varGamma _{40}$ curve with that obtained from the experiments of Valani et al. (Reference Valani, Slim and Simula2019) for a droplet radius of 
$R=(0.63\pm 0.03)$ mm, shows good agreement highlighting the success of the present model for small- to moderate-sized superwalkers.
4. Discussion and conclusions
We have studied the dynamics of bouncing droplets on a vibrating liquid bath under two-frequency driving using the theoretical model of Moláček & Bush (Reference Moláček and Bush2013a,Reference Moláček and Bushb) and a new model for the wave field to understand the emergence of superwalkers. We have shown that two-frequency driving at 
$f$ and 
$f/2$ with an appropriately chosen phase difference 
${\rm \Delta} \phi$ lifts every second peak and lowers the intermediate peaks in the bath's motion. This allows larger droplets to bounce in a resonant 
$(1,2,1)$ mode where they can efficiently excite damped subharmonic Faraday waves that enable them to superwalk. We note that superwalking would not be expected for two arbitrary frequency combinations, as the lowering of every second peak is crucial for droplets to remain in a 
$(1,2,1)$ mode. For example, for two-frequency driving at 
$f=80$ and 
$4f/5=64$ Hz, Sampara & Gilet (Reference Sampara and Gilet2016) reported chaotic bouncing modes with irregular walking at typical speeds of only 
$5\ \text {mm}\ \textrm {s}^{-1}$.
We have shown that the phase difference 
${\rm \Delta} \phi$ plays a crucial role in the dynamics of superwalking droplets because it controls the relative amplitudes of two succeeding peaks in one full cycle of the bath's motion. Fast superwalking occurs for phase differences between 
$130^\circ$ and 
$180^\circ$ where there is a larger difference between these amplitudes, while phase differences around 
$45^\circ$, where the amplitude difference between succeeding peaks is small, correspond to stationary bouncing or coalescence. Extending these observations, it will be interesting to use our model to explore a system with slightly detuned driving frequencies of 
$f$ and 
$f/2+\epsilon$ with 
$\epsilon \ll f/2$. This is equivalent to driving at frequencies of 
$f$ and 
$f/2$ with a slowly varying phase difference 
${\rm \Delta} \phi (t)={\rm \Delta} \phi + 2{\rm \pi} \epsilon t$. Such traversal of the phase difference gives rise to new types of locomotion in the walking droplet system where, e.g. the droplet alternates periodically between superwalking and stationary bouncing. Such behaviour was observed experimentally by Valani et al. (Reference Valani, Slim and Simula2019) and was coined stop-and-go motion. We aim to discuss the details of such motion in a future work.
On comparing the speed-size characteristics of simulated superwalkers with the experimental results of Valani et al. (Reference Valani, Slim and Simula2019), we find excellent agreement on the ascending branch, with 
$(1,2,1)^{{H}}$ superwalkers primarily observed. These observations also explain the good agreement noted by Valani et al. (Reference Valani, Slim and Simula2019) between superwalking speeds obtained in experiments and those predicted using the stroboscopic model of Oza et al. (Reference Oza, Rosales and Bush2013) (dashed curve in figure 5a). The latter is a reduced form of the full Moláček & Bush (Reference Moláček and Bush2013a,Reference Moláček and Bushb) model predicated on a 
$(2,1)^{H}$ bouncing mode and our two-frequency model would reduce to essentially the same model for such modes.
The superwalking speed of larger superwalkers is not captured well by the current model. This suggests that the model does not include the fundamental mechanism that allows the largest superwalkers to walk, and even exist. Indeed, Valani et al. (Reference Valani, Slim and Simula2019) noted that the largest superwalkers on the descending branch undergo significant internal deformations (Valani et al. Reference Valani, Slim and Simula2019). In appendix B, we present results incorporating deformation of the droplets by modelling them as a vertical spring following Blanchette (Reference Blanchette2016) and Gilet et al. (Reference Gilet, Terwagne, Vandewalle and Dorbolo2008), and find this to have a limited effect on the speed-size curve. We also present results using the nonlinear logarithmic spring model of Moláček & Bush (Reference Moláček and Bush2013b) for the vertical dynamics with no better success. Another observation made by Valani et al. (Reference Valani, Slim and Simula2019) was that larger superwalkers have a prolonged contact time with the bath. This prolonged contact time, potentially in combination with internal deformation, may change the wave field in the vicinity of the droplet and the long-time approximation of the standing wave field in (2.12) may break down. Perhaps a more refined modelling of the system that incorporates the detailed contact interaction between the droplet and the bath, the wave evolution and droplet deformations might be required to capture the behaviour of these larger superwalkers.
Acknowledgements
We acknowledge financial support from an Australian Government Research Training Program (RTP) Scholarship (R.N.V.) and the Australian Research Council via the Future Fellowship Project No. FT180100020 (T.P.S.).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of waves generated by a superwalker
To derive the form of the surface waves generated by a superwalker, we closely follow the approach of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) who considered walkers driven at a single frequency. We consider an incompressible, Newtonian liquid in a bath that is infinitely large in horizontal extent and infinitely deep. The bath is subjected to periodic vertical vibrations that result in a modulation of the effective gravity in the frame of the bath 
${g}^*(t)={g}[1+\varGamma _{f} \sin (\varOmega t) + \varGamma _{f/2}\sin (\varOmega t/2 + {\rm \Delta} \phi )]$, where 
${g}$ is the gravitational acceleration constant and 
$\varGamma _f$ and 
$\varGamma _{f/2}$ are dimensionless acceleration amplitudes of the two frequencies with a relative phase difference of 
${\rm \Delta} \phi$. For the sake of notational clarity, we will refer to specific frequencies 
$f=\varOmega /2{\rm \pi} =80$ and 
$f/2=40$ Hz, but the derivation is general. The horizontal coordinates are 
$(x,y)$ and the vertical coordinate is 
$z$ with the origin located on the free surface of the undeformed liquid. The evolution of the liquid is governed by the incompressible Navier–Stokes equations,
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0 \quad \text{and} \quad \frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{v}=-\frac{1}{\rho}\boldsymbol{\nabla} P - {g}^*(t)\boldsymbol{\hat{k}}+\nu\nabla^2\boldsymbol{v}, \end{equation}
where 
$\boldsymbol {v}(\boldsymbol {r},t)$ is the velocity field, 
$P(\boldsymbol {r},t)$ is the pressure field relative to atmospheric pressure, 
$\boldsymbol {r}=(x,y,z)$ is the position vector, 
$\boldsymbol {\hat {k}}$ is a unit vector in the z direction, 
$\rho$ is the density and 
$\nu$ is the kinematic viscosity. At the free surface of the liquid, 
$z=h(x,y,t)$, the kinematic boundary condition implies
\begin{equation} \frac{\partial h}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} (h-z)=0, \end{equation}while balancing stresses requires
\begin{equation} -P\boldsymbol{\hat{n}}+\rho\nu\left(\boldsymbol{\nabla}\boldsymbol{v} + \boldsymbol{\nabla} \boldsymbol{v}^\textrm{T}\right)\boldsymbol{\cdot} \boldsymbol{\hat{n}}=-P^{{ext}}\boldsymbol{\hat{n}}-\sigma (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\hat{n}})\boldsymbol{\hat{n}}, \end{equation}
where 
$\sigma$ is the coefficient of surface tension, 
$P^{{ext}}(x,y,t)$ is the pressure exerted by the droplet on the free surface through the intervening air layer and 
$\boldsymbol {\hat {n}}$ is the unit normal out of the liquid. Assuming that the pressure distribution imparted by the droplet during contact is uniform in the contact region we get
\begin{equation} P^{{ext}}(x,y,t)=P^{{ext}}(t)=F_N(t)/{\rm \pi} w^2, \end{equation}
where 
$w$ is the effective radius of the contact area and 
$F_N(t)$ is the normal force as described in § 2.1.
We consider small perturbations from the stationary equilibrium state 
$\boldsymbol {{v}}=\boldsymbol {0}$, 
${{P}}=-\rho {g}^*(t) z$ and 
$h=0$. We linearise the above equations about this equilibrium state with pressure perturbation 
$p(\boldsymbol {r},t)$, velocity perturbation 
$\boldsymbol {v}(\boldsymbol {r},t)$ and free surface perturbation 
$h(\boldsymbol {r},t)$ and follow the steps outlined in § 2.1 of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018). For the remainder of this Appendix, we will use the dimensionless time 
$\tau =\varOmega t/2$ for ease of comparison with their equations and will revert back to the dimensional time 
$t$ in the final expressions. Taking a Fourier transform in 
$x$ and 
$y$ of the linearised equations, followed by a Laplace transform and some algebra, we obtain the equation
\begin{equation} f_k(s)h_{\boldsymbol{k},s}+\frac{2\varGamma_{80} {g} k}{\text{i}\varOmega^{2}}(h_{\boldsymbol{k},s-2\text{i}}- h_{\boldsymbol{k},s+2\text{i}})+\frac{2\varGamma_{40} {g} k}{\text{i}\varOmega^{2}}(h_{\boldsymbol{k},s-\text{i}} \,\text{e}^{\textrm{i}{\rm \Delta}\phi}-h_{\boldsymbol{k},s+\text{i}}\,\text{e}^{-\textrm{i}{\rm \Delta}\phi})+\varPi_{\boldsymbol{k},s}=0 \end{equation}
for the transformed free surface 
$h_{\boldsymbol {k},s}$. Here, the Fourier transform for an arbitrary variable 
$X(x,y,\tau )$ is defined as
\begin{equation} X_{\boldsymbol{k}}(\tau) = \int_0^\infty \int_0^\infty X(x,y,\tau) \exp[-\text{i}(k_x x + k_y y)]\,\text{d}x\,\text{d}y \end{equation}
with 
$k=|\boldsymbol {k}|$ and the Laplace transform is defined as
\begin{equation} X_{\boldsymbol{k},s} = \int_0^\infty X_{\boldsymbol{k}}(\tau) \,\text{e}^{-s\tau}\,\text{d}\tau. \end{equation}Furthermore,
\begin{equation} f_k(s) = (s+\gamma_k)^2-\gamma_k^{3/2}\sqrt{\gamma_k+2s}+\omega^2_{k}, \end{equation}
with 
$\gamma _k=4\nu k^2/\varOmega$ and 
$\omega ^2_k=4[{g} k+(\sigma /\rho )k^3]/\varOmega ^2$. This function and all its derivatives obey 
$f_k(\bar {z})=\overline {f_k(z)}$ where the overline denotes complex conjugation. The last term of (A 5) described the transformed pressure distribution from droplet's impact 
$\varPi _{\boldsymbol {k},s}=(4k/\varOmega ^2\rho )P^{{ext}}_{\boldsymbol {k},s}$. Using the definition of Fourier transform and (A 4) we get assuming 
$wk \ll 1$, 
$P^{{ext}}_{\boldsymbol {k},s}=P^{{ext}}(s)\int _0^w \text {J}_0(kr)r\,\text {d}r\approx F_N(s)/2{\rm \pi}$ with 
$r=\sqrt {x^2+y^2}$. Hence we obtain
\begin{equation} \varPi_{\boldsymbol{k},s}=\frac{2k}{{\rm \pi}\varOmega^2\rho}F_N(s), \end{equation}
where 
$P^{{ext}}(s)$ and 
$F_N(s)$ are the Laplace transforms of 
$P^{{ext}}(\tau )$ and 
$F_N(\tau )$ respectively. We note that (A 5) reduces to (2.2) of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) on setting 
$\varGamma _{40}=0$, with the caveat that our driving is a sine function while Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) use a cosine. We have chosen a sine function for driving for the sake of consistency with the experiments of Valani et al. (Reference Valani, Slim and Simula2019).
We first consider Faraday waves in the absence of external pressure perturbations, which reduces (A 5) to
\begin{equation} f_k(s)h_{\boldsymbol{k},s}+\frac{2\varGamma_{80} {g} k}{\text{i}\varOmega^{2}} (h_{\boldsymbol{k},s-2\text{i}}-h_{\boldsymbol{k},s+2\text{i}})+ \frac{2\varGamma_{40} {g} k}{\text{i}\varOmega^{2}} (h_{\boldsymbol{k},s-\text{i}}\,\text{e}^{\textrm{i}{\rm \Delta}\phi}-h_{\boldsymbol{k},s+\text{i}}\,\text{e}^{-\textrm{i}{\rm \Delta}\phi})=0. \end{equation}Due to the periodic driving of the system, a Floquet ansatz is appropriate in the time domain. The form we assume and its corresponding Laplace transform are given by (Besson, Edwards & Tuckerman Reference Besson, Edwards and Tuckerman1996)
\begin{equation} h_{\boldsymbol{k}}(\tau)=\sum_{l=-\infty}^{\infty} h^{(l)}_{\boldsymbol{k}}\,\text{e}^{\text{i}l\tau}\,\text{e}^{\delta_{k}\tau}\quad \text{and}\quad h_{\boldsymbol{k},s}=\sum_{l=-\infty}^{\infty} \frac{h^{(l)}_{\boldsymbol{k}}}{s-\text{i}l-\delta_k}. \end{equation}
Here, 
$\delta _k$ is a complex perturbation whose real part vanishes when the Faraday instability threshold is reached. Substituting this form into (A 10), we obtain
\begin{align} & \sum_{l=-\infty}^{\infty} h^{(l)}_{\boldsymbol{k}} \left[\frac{f_k(s)}{s-\text{i}l-\delta_k} - \text{i}\varGamma_{80}\beta_k \left( \frac{1}{s-\text{i}(l+2)-\delta_k} - \frac{1}{s-\text{i}(l-2)-\delta_k}\right) \right. \nonumber\\ &\quad \left.-\text{i}\varGamma_{40} \beta_k \left( \frac{1}{s-\text{i}(l+1)-\delta_k}\,\text{e}^{\text{i}{\rm \Delta}\phi} - \frac{1}{s-\text{i}(l-1)-\delta_k}\,\text{e}^{-\text{i}{\rm \Delta}\phi}\right) \right]=0, \end{align}
with 
$\beta _k=2{g}k/\varOmega ^2$. Using the Heaviside cover-up method (Thomas & Finney Reference Thomas and Finney1996) yields an infinite-dimensional linear system 
${\boldsymbol{\mathsf{A}}} \boldsymbol {h}=\boldsymbol {0}$ coupling the Floquet components together (Kumar & Tuckerman Reference Kumar and Tuckerman1994; Besson et al. Reference Besson, Edwards and Tuckerman1996; Kumar Reference Kumar1996; Tadrist et al. Reference Tadrist, Shim, Gilet and Schlagheck2018). Here 
${\boldsymbol{\mathsf{A}}}$ is the pentadiagonal matrix
\begin{align} {\boldsymbol{\mathsf{A}}} = \begin{bmatrix} \ddots & \ddots & \ddots & \vdots & \vdots & \vdots & \dots \\ \ddots & f_k(-2\text{i}+\delta_k) & \Upsilon_k & \alpha_k & 0 & 0 & \dots \\ \ddots & \bar{\Upsilon}_k & f_k(-1\text{i}+\delta_k) & \Upsilon_k & \alpha_k & 0 & \dots\\ \dots & \bar{\alpha}_k & \bar{\Upsilon}_k & f_k(\delta_k) & \Upsilon_k & \alpha_k & \dots\\ \dots & 0 & \bar{\alpha}_k & \bar{\Upsilon}_k & f_k(1\text{i}+\delta_k) & \Upsilon_k & \ddots\\ \dots & 0 & 0 & \bar{\alpha}_k & \bar{\Upsilon}_k & f_k(2\text{i}+\delta_k) & \ddots\\ \dots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots \end{bmatrix}, \end{align}
with 
$\alpha _k=\text {i}\varGamma _{80} \beta _k$ and 
$\Upsilon _k=\text {i}\varGamma _{40}\beta _k\,\text {e}^{-\text {i}{\rm \Delta} \phi }$, and 
$\boldsymbol {h}$ is a vector of the Floquet components 
$h^{(l)}_{\boldsymbol {k}}$. To obtain non-trivial solutions of this linear system requires
\begin{equation} \text{det}({\boldsymbol{\mathsf{A}}})=0. \end{equation}A.1. Decay rate of damped Faraday waves and the Faraday instability threshold
Solving (A 14) for fixed 
$\varGamma _{80}$, 
$\varGamma _{40}$ and 
${\rm \Delta} \phi$, we obtain 
$\delta _k$ as a function of the wavenumber 
$k$. Below the Faraday instability threshold, this corresponds to a decay rate for the waves 
$\text {Re}(\delta _k)$ and a dispersion relation 
$\text {Im}(\delta _k)$. Results for typical parameter values of superwalkers are shown in figure 7. Figure 7(a) shows the numerically converged 
$\text {Re}(\delta _k)$ as a function of 
$k$ (solid curves). We see two different Faraday windows, one in which the waves are locked at 
$\text {Im}(\delta _k)=1/2$ (the blue-shaded region in figures 7a,b) and one in which waves are locked at 
$\text {Im}(\delta _k)=0$ (the red-shaded region in figures 7a,c). In each of these windows, we see two different branches for the decay rate 
$\text {Re}(\delta _k)$, an upper branch corresponding to a slowly decaying wave and a lower branch corresponding to a more rapidly decaying wave.
Properties of two-frequency, damped Faraday waves. (a) Decay rate 
$\text {Re}(\delta _k)$ as a function of wavenumber 
$k$ for 
$\varGamma _{80}=3.8$, 
$\varGamma _{40}=0.6$ and 
${\rm \Delta} \phi =130^{\circ }$ using a 
$21$-mode truncation corresponding to 
$|l|\le 10$ (solid black curve). The blue and red dotted curves show the decay rate of the slowly decaying wave using the two-mode approximation 
$\text {Re}(\delta ^{+}_{k20})$ in the blue Faraday window and the two-mode approximation 
$\text {Re}(\delta ^{+}_{k40})$ in the red Faraday window, respectively. The grey dashed curves are second-order approximations to these decay rates at the peak values in each Faraday window. Panels (b,c) show the dispersion relation 
$\text {Im}(\delta _k)$ in the two Faraday windows. In (d,e), the magnitude of the amplitudes 
$h_{kF40}^{(l)}$ and 
$h_{kF20}^{(l)}$ of the different modes 
$l$ are shown at the most unstable wavenumbers in each Faraday window using the 
$21$-mode truncation, with the dominant modes coloured. These correspond to the eigenvectors of 
${\boldsymbol{\mathsf{A}}} \boldsymbol {h}=\boldsymbol {0}$ with eigenvalue 
$0$. Note that these amplitudes only yield information about the relative values of each mode.
To obtain analytical forms of these results, we truncate the infinite dimensional matrix equation to a few dominant modes. For the 
$\text {Im}(\delta _k)=0$ (red) Faraday window in figure 7(a), we find that the dominant contribution to the amplitude is from the two modes with 
$l=\pm 1$ (see figure 7e) corresponding to a frequency of 
$\pm 40$ Hz. Denoting the decay rate in this Faraday window by 
$\text {Re}(\delta _{k40})$ and using this two-mode approximation, (A 14) reduces to
\begin{equation} f_k(-\text{i}+\text{Re}(\delta_{k40}))f_k(\text{i}+\text{Re}(\delta_{k40}))-|\alpha_k|^2=0. \end{equation} We can obtain a good approximation to this decay rate by following an approach similar to § 2.2.2 of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) and expanding the function 
$f_k(\pm \text {i}+\text {Re}(\delta _{k40}))$ to second order in the small decay rate 
$\text {Re}(\delta _{k40})$ to get
\begin{equation} \text{Re}(\delta_{k40}^{\pm})=-\frac{b_{k}(\text{i})}{2a_{k}(\text{i})} \left(1 \mp \sqrt{1-\frac{4a_{k}(\text{i})c_{k}(\text{i},\alpha_k)}{b^2_{k}(\text{i})}} \right), \end{equation}where the functions
\begin{gather} a_k(u)=\dot{f}_k(u)\dot{f}_k(-u)+\tfrac{1}{2}\ddot{f}_k(u)f_k(-u)+\tfrac{1}{2}\ddot{f}_k(-u)f_k(u), \end{gather}
\begin{gather}b_k(u)=\dot{f}_k(u)f_k(-u)+\dot{f}_k(-u)f_k(u), \end{gather}
\begin{gather}c_k(u,Z)=f_k(u)f_k(-u)-|Z|^2. \end{gather} Here, 
$\text {Re}(\delta _{k40}^{+})$ and 
$\text {Re}(\delta _{k40}^{-})$ correspond to the decay rates of the slowly and quickly decaying wave respectively. This approximation for the slowly decaying wave 
$\text {Re}(\delta _{k40}^{+})$ is shown as a red, dotted curve in figure 7(a). We can further approximate this decay rate near the most unstable wavenumber 
$k_{F40}$ by
\begin{equation} \text{Re}(\delta_{k40}^{+}) \approx \text{Re}(\delta_{F40}^{+}) - D_{40} (k-k_{F40})^2, \end{equation}
where 
$\text {Re}(\delta _{F40}^{+})=\lim _{k\to k_{F40}}\text {Re}(\delta _{k40}^{+})$ and 
$D_{40}=-\frac {1}{2}\text {d}^2\text {Re}(\delta _{k40}^{+})/\text {d}k^2|_{k=k_{F40}}$ is the diffusion coefficient, both of which can be calculated from (A 16). This approximation of 
$\text {Re}(\delta _{k40}^{+})$ is shown as a grey, dashed curve in figure 7(a).
We follow a similar approach to obtain an analytical expression for the decay rate in the 
$\text {Im}(\delta _k)=1/2$ (blue) Faraday window in figure 7(a). In this window, the dominant contribution is from the 
$l=-1$ and 
$0$ modes (see figure 7d), corresponding to frequencies 
${\pm } 20$ Hz. Using this two-mode approximation and denoting the decay rate by 
$\text {Re}(\delta _{k20})$, (A 14) reduces to
\begin{equation} f_k(-\text{i}/2+\text{Re}(\delta_{k20}))f_k(\text{i}/2+\text{Re}(\delta_{k20}))-|\Upsilon_k|^2=0. \end{equation}
A good approximation for this decay rate is obtained by expanding the function 
$f_k(\pm \text {i}/2+\text {Re}(\delta _{k20}))$ to second order, giving us
\begin{equation} \text{Re}(\delta_{k20}^{\pm})=-\frac{b_{k}(\text{i}/2)}{2a_{k}(\text{i}/2)} \left(1 \mp \sqrt{1-\frac{4a_{k}(\text{i}/2)c_{k}(\text{i}/2,\Upsilon_k)}{b^2_{k}(\text{i}/2)}} \right), \end{equation}
where 
$\text {Re}(\delta _{k20}^{+})$ and 
$\text {Re}(\delta _{k20}^{-})$ correspond to the decay rates of the slowly and quickly decaying wave respectively. We can further approximate 
$\text {Re}(\delta _{k20}^{+})$ near the most unstable wavenumber 
$k_{F20}$ by
\begin{equation} \text{Re}(\delta^{+}_{k20}) \approx \text{Re}(\delta_{F20}^{+}) - D_{20} (k-k_{F20})^2, \end{equation}
where 
$\text {Re}(\delta _{F20}^{+})=\lim _{k\to k_{F20}}\text {Re}(\delta _{k20}^{+})$ and 
$D_{20}=-\frac {1}{2}\text {d}^2\text {Re}(\delta ^{+}_{k20})/\text {d}k^2|_{k=k_{F20}}$ is the diffusion coefficient corresponding to this Faraday window. These approximations of 
$\text {Re}(\delta ^{+}_{k20})$ from (A 22) and (A 23) are shown in figure 7(a) as a blue dotted and a grey dashed curve respectively.
When 
$\text {Re}(\delta _k)> 0$ for any wavenumber 
$k$, growing Faraday waves are predicted. For two-frequency driving at 
$f$ and 
$f/2$, either 
$f/2$ Faraday waves or 
$f/4$ Faraday waves can emerge depending on the relative strength of the acceleration amplitudes and the phase difference (Müller Reference Müller1993; Valani et al. Reference Valani, Slim and Simula2019). The marginal stability curves representing the acceleration amplitudes at onset of Faraday waves, 
$\varGamma _{F80}$ and 
$\varGamma _{F40}$, can be found by setting 
$\text {Re}(\delta _k)=0$ when solving (A 14). From figure 7(a), we see two Faraday windows where 
$\text {Re}(\delta _k)$ can potentially cross zero corresponding to either the 
$f/2$ instability of frequency 
$40$ Hz or the 
$f/4$ instability of frequency 
$20$ Hz. Figure 8(a) shows the comparison of the numerically converged marginal stability curve obtained using a 
$21$-mode truncation (yellow solid curve) and the two-mode approximation for the 
$20$ Hz (blue dashed curve) and 
$40$ Hz (red dashed curve) Faraday waves, with the experiments of Valani et al. (Reference Valani, Slim and Simula2019) (circles). Figure 8(b) shows the numerically converged marginal stability curves at different phase differences 
${\rm \Delta} \phi$. We note that changes in 
${\rm \Delta} \phi$ cause appreciable changes in the 
$20$ Hz Faraday threshold.
Faraday threshold curves for two-frequency driving. (a) Comparison of the Faraday threshold curves for 
${\rm \Delta} \phi =130^{\circ }$ obtained using 
$21$ modes (solid yellow curves) and using the two-mode approximations (dotted curves) together with the experimental results (circles) of Valani et al. (Reference Valani, Slim and Simula2019). For the latter, empty circles indicate that flat liquid surfaces were observed while filled circles indicate that Faraday waves were observed. (b) Faraday thresholds for different values of the phase difference 
${\rm \Delta} \phi$ using a 
$21$-mode calculation.
A.2. Amplitude and phase shift of damped Faraday waves
In figure 7(d,e), the relative amplitudes of the Floquet modes are shown for the slowest decaying modes in the 
$20$ Hz and 
$40$ Hz Faraday windows respectively. We now turn to calculating these amplitudes for our reduced-mode approximations and use these to obtain expressions for the wave profile generated by a single bounce of a droplet.
Using the two-mode approximation for the 
$40$ Hz window, we can write the Floquet ansatz in (A 11a,b) as
\begin{equation} h_{\boldsymbol{k}}(\tau) \approx \left(h^{(-1)}_k\,\text{e}^{-\text{i}\tau} + h^{(1)}_k\,\text{e}^{\text{i}\tau}\right) \text{e}^{\delta_{k}\tau}, \end{equation}
and the infinite-dimensional linear system 
${\boldsymbol{\mathsf{A}}} \boldsymbol {h}=\boldsymbol {0}$ reduces to a 
$2\times 2$ matrix system 
${\boldsymbol{\mathsf{A}}}_2 \boldsymbol {h}_2=\boldsymbol {0}$. Since the determinant of the matrix 
${\boldsymbol{\mathsf{A}}}_2$ is zero, we obtain the amplitudes 
$\boldsymbol {h}_2$ from the null space vector, which gives 
$h^{(-1)}_k=c\xi ^{\pm }_{40}$ with 
$\xi ^{\pm }_{40}=\text {i}\sqrt {\alpha _k/f_k(-\text {i}+\delta ^{\pm }_{k40})}$ and 
$h^{(1)}_k=\overline {h^{(-1)}_k}$, where 
$c$ is a free parameter. Substituting this solution into (A 24) and using 
$\text {Im}(\delta ^{\pm }_{k40})=0$ in this window, we obtain the wave forms
\begin{equation} h^{\pm}_{\boldsymbol{k}40}(\tau) = c \left[\xi^{\pm}_{40}\,\text{e}^{-\text{i}\tau} + \bar{\xi}^{\pm}_{40}\,\text{e}^{\text{i}\tau}\right] \exp({\text{Re}(\delta^{\pm}_{k40})\tau}). \end{equation}Thus, the total wave field in this Faraday window can be represented as
\begin{equation} h_{\boldsymbol{k}40}(\tau) = \zeta_{40}^{+} \exp({\text{Re}(\delta_{k40}^{+})\tau}) \cos(\tau+\theta_{k40}^{+}) + \zeta_{40}^{-} \exp({\text{Re}(\delta_{k40}^{-})\tau})\cos(\tau+\theta_{k40}^{-}), \end{equation}where
\begin{equation} \theta_{k40}^{\pm}=\arctan\left(\frac{-\text{Im}(\xi_{40}^{\pm})}{\text{Re}(\xi_{40}^{\pm})} \right). \end{equation}
and 
$\zeta ^{\pm }_{40}=2c|\xi ^{\pm }_{40}|$. These (A 27) and (A 26) are equivalent to (2.47) and (2.48) of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018). Similar to § 2.3 of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018), we now continue by modelling the temporal profile of a droplet's impact by a delta function. The corresponding pressure and force exerted by the droplet on the liquid is then
\begin{equation} \varPi_{\boldsymbol{k}}(\tau)=(2k/{\rm \pi}\varOmega^2\rho)F_N(\tau)=\delta(\tau - \tau_i)v_{\boldsymbol{k}}. \end{equation}
By integrating the time domain version of (A 5) across the delta kick, we find that 
$v_{\boldsymbol {k}}$ corresponds to negative change of velocity of 
$h_{\boldsymbol {k}}$ following impact. If the surface is perfectly flat and at rest before the impact, the wave profile is axisymmetric i.e. 
$h_{\boldsymbol {k}}(\tau )=h_{k}(\tau )$. Using the initial conditions as 
$\tau \xrightarrow {} \tau _i$ that 
$h_{\boldsymbol {k}}=0$ and 
${\partial h_{\boldsymbol {k}}}/{\partial \tau }=-v_k$ we get 
$\zeta _{40}^{\pm }=v_k \alpha _{40}^{\pm }$, where
\begin{align} \alpha_{40}^{\pm}= \frac{-2\cos(\tau_i+\theta_{k40}^{\mp})\exp({-\text{Re}(\delta_{k40}^{\pm})\tau_i})}{(\text{Re}(\delta_{k40}^{\pm})-\text{Re}(\delta_{k40}^{\mp}))(\cos(2\tau_i+\theta_{k40}^{\pm}+\theta_{k40}^{\mp})+\cos(\theta_{k40}^{\pm}-\theta_{k40}^{\mp}))-2\sin(\theta_{k40}^{\pm}-\theta_{k40}^{\mp})}. \end{align} Taking a similar approach, we can obtain an equation for the wave field by using the two-mode approximation in the 
$20$ Hz Faraday window. The two-mode approximation of (A 11a,b) gives
\begin{equation} h^{\pm}_{\boldsymbol{k}20}(\tau) = (h^{(-1)}_k\,\text{e}^{-\text{i}\tau/2} + h^{(0)}_k \,\text{e}^{\text{i}\tau/2}) \exp({\text{Re}(\delta^{\pm}_{k20})\tau}). \end{equation} Solving for the null space of the matrix equation we get 
$h^{(-1)}_k=c\xi ^{\pm }_{20}$ with 
$\xi ^{\pm }_{20}=\text {i}\sqrt {\Upsilon _k/f_k(-\text {i}/2+\text {Re}(\delta ^{\pm }_{k20}))}$ and 
$h^{(0)}_k=\overline {h^{(-1)}_k}$, where 
$c$ is a free parameter.
For this we obtain the amplitudes
\begin{equation} h^{\pm}_{\boldsymbol{k}20}(\tau) = c(\xi^{\pm}_{20}\,\text{e}^{-\text{i}\tau/2} + \bar{\xi}^{\pm}_{20}\,\text{e}^{\text{i}\tau/2}) \exp({\text{Re}(\delta^{\pm}_{k20})\tau}). \end{equation}Hence, we can express the total wave field for this Faraday window as
\begin{align} h_{\boldsymbol{k}20}(\tau) = \zeta_{20}^{+}\exp({\text{Re}(\delta_{k20}^{+})\tau})\cos(\tau/2 + \theta_{k20}^{+})+\zeta_{20}^{-}\exp({\text{Re}(\delta_{k20}^{-})\tau})\cos(\tau/2 + \theta_{k20}^{-}), \end{align}where
\begin{equation} \theta_{k20}^{\pm} = \arctan\left( \frac{-\text{Im}(\xi_{20}^{\pm})}{\text{Re}(\xi_{20}^{\pm})} \right). \end{equation}
Using the same initial conditions as for 
$40$ Hz waves we get 
$\zeta _{20}^{\pm }=v_k\alpha _{20}^{\pm }$, where
\begin{align} \alpha_{20}^{\pm}=\frac{-2\cos(\tau_i/2+\theta_{k20}^{\mp})\exp({-\text{Re}(\delta_{k20}^{\pm})\tau_i})}{(\text{Re}(\delta_{k20}^{\pm})-\text{Re}(\delta_{k20}^{\mp}))(\cos(\tau_i+\theta_{k20}^{\pm}+\theta_{k20}^{\mp}) +\cos(\theta_{k20}^{\pm}-\theta_{k20}^{\mp}))-2\sin(\theta_{k20}^{\pm}-\theta_{k20}^{\mp})}.\end{align}A.3. Wave field of a superwalker
For late times after the impact, 
$\tau \gg \tau _i$, and when the acceleration amplitudes are close to their respective Faraday thresholds, 
$\varGamma _{80}\lesssim \varGamma _{F80}$ and 
$\varGamma _{40}\lesssim \varGamma _{F40}$, the wave field is dominated by the slowly decaying Faraday waves from both the 
$40$ and 
$20$ Hz modes. Hence, we can approximate the final wave field generated by the impact of the droplet as
\begin{align} h_{\boldsymbol{k}}(\tau)= \alpha_{40}^{+}v_k \exp({\text{Re}(\delta_{k40}^{+})\tau}) \cos(\tau+\theta_{k40}^{+})+\alpha_{20}^{+}v_k\exp({\text{Re}(\delta_{k20}^{+})\tau})\cos(\tau/2 + \theta_{k20}^{+}). \end{align}Transforming back to the spatial domain with an inverse Fourier transform yields
\begin{equation} h(x,y,\tau)=\frac{1}{2{\rm \pi}}\int_0^\infty \int_0^\infty h_{\boldsymbol{k}}(\tau) \exp[{\text{i}(k_x x +k_y y)}] \,\text{d}k_x \,\text{d}k_y. \end{equation}Since the wave profile is radially symmetric, the above inverse Fourier transform reduces to an inverse Hankel transform,
\begin{equation} h(x,y,\tau)=\int_0^\infty h_{\boldsymbol{k}}(\tau) \text{J}_0(kr)k\,\text{d}k. \end{equation}Hence, the wave profile in the real space is given by
\begin{align} h(x,y,\tau) &= \int_0^\infty B_{k40}^{+}v_k \exp({\text{Re}(\delta_{k40}^{+})(\tau-\tau_i)}) \cos(\tau+\theta_{k40}^{+})\text{J}_0(kr)k\,\text{d}k \nonumber\\ &\quad +\int_0^\infty B_{k20}^{+}v_k\exp({\text{Re}(\delta_{k20}^{+})(\tau-\tau_i)})\cos(\tau/2 + \theta_{k20}^{+}) \text{J}_0(kr)k\,\text{d}k, \end{align}
where 
$B_{k40}^{+}=\alpha _{40}^{+}\exp ({\text {Re}(\delta _{k40}^{+})\tau _i})$ and 
$B_{k20}^{+}=\alpha _{20}^{+}\exp ({\text {Re}(\delta _{k20}^{+})\tau _i})$. Using the second-order expansion for 
$\text {Re}(\delta ^{+}_{k40})$ and 
$\text {Re}(\delta ^{+}_{k20})$ in (A 20) and (A 23), we get the following approximation to the above integral in the limit 
$\tau \xrightarrow {}\infty$ (for details see appendix C of Tadrist et al. Reference Tadrist, Shim, Gilet and Schlagheck2018)
\begin{align} & h(x,y,\tau)= \tilde{A}^{(0)}_{40}(\tau_i)\frac{\cos(\tau+\theta_{F40}^{+})}{\sqrt{\tau-\tau_i}} \text{J}_0(k_{F40}|\boldsymbol{x}-\boldsymbol{x}_i|)\exp{\left[ -\frac{\tau-\tau_i}{2{\rm \pi} {Me}_{40}} - \frac{|\boldsymbol{x}-\boldsymbol{x}_i|^2}{4 D_{40} (\tau-\tau_i)} \right]} \nonumber\\ &\quad + \tilde{A}^{(0)}_{20}(\tau_i)\frac{\cos(\tau/2 + \theta_{F20}^{+})}{\sqrt{\tau-\tau_i}} \text{J}_0(k_{F20}|\boldsymbol{x}-\boldsymbol{x}_i|)\exp{\left[ -\frac{\tau-\tau_i}{2{\rm \pi} {Me}_{20}} - \frac{|\boldsymbol{x}-\boldsymbol{x}_i|^2}{4 D_{20} (\tau-\tau_i)} \right]}, \end{align}
where 
$\boldsymbol {x}_i$ is the location of the impact and the memory parameters 
${Me}_{40}$ and 
${Me}_{20}$ are given by 
${Me}_{40}=-1/2{\rm \pi} \text {Re}(\delta _{F40}^{+})$ and 
${Me}_{20}=-1/2{\rm \pi} \text {Re}(\delta _{F20}^{+})$. Furthermore,
\begin{equation} \tilde{A}^{(0)}_{40}=k_{F40}\sqrt{\frac{\rm \pi}{D_{40}}}v_{k}B_{F40}^{+}(\tau_i)\quad \text{and}\quad \tilde{A}^{(0)}_{20}=k_{F20}\sqrt{\frac{\rm \pi}{D_{20}}}v_{k}B_{F20}^{+}(\tau_i). \end{equation}
To include a finite contact time, we follow the suggestion in Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) of using Duhamel's principle and the approach used in appendix A.4 of Moláček & Bush (Reference Moláček and Bush2013b), and integrate the impulse response with a time varying impact signal 
$\varPi _{\boldsymbol {k}}(\tau )$. This results in replacing the amplitude coefficients 
$\tilde {A}^{(0)}_{40}$ and 
$\tilde {A}^{(0)}_{20}$ by
\begin{gather} \tilde{A}_{40}=k_{F40}\sqrt{\frac{\rm \pi}{D_{40}}}\int_{\tau_n^i}^{\tau_n^c}B^{+}_{F40}(\tau') \varPi_{\boldsymbol{k}}(\tau') \,\text{d}\tau', \end{gather}
\begin{gather}\tilde{A}_{20}=k_{F20}\sqrt{\frac{\rm \pi}{D_{20}}}\int_{\tau_n^i}^{\tau_n^c}B^{+}_{F20}(\tau') \varPi_{\boldsymbol{k}}(\tau')\,\text{d}\tau'. \end{gather}
We change the dimensionless time 
$\tau$ back to dimensional time 
$t$ and replace 
$\varPi _{\boldsymbol {k}}(t)$ by 
$(2k/{\rm \pi} \varOmega ^2 \rho )F_N(t)$ using (A 9). We also replace the initial contact time 
$t_i$ and location of contact 
$\boldsymbol {x}_i$ by their weighted average values 
$t_n$ and 
$\boldsymbol {x}_n$ as given in (2.6a,b), and replace the dimensionless amplitudes 
$\tilde {A}_{40}$ and 
$\tilde {A}_{20}$ by 
$A_{40}=\sqrt {2/\varOmega }\tilde {A}_{40}$ and 
$A_{20}=\sqrt {2/\varOmega }\tilde {A}_{20}$ which gives (2.13b) and results in the wave field equation (2.12).
Appendix B. Comparison of different droplet models
To test the robustness of the superwalking behaviour, we explored superwalkers using alternative models for the vertical dynamics, the wave field generated, and adding droplet deformations to the model presented in § 2. Comparison of these models with the model presented in this paper and the experimental results of Valani et al. (Reference Valani, Slim and Simula2019) for a typical speed-size curve of superwalkers is presented in figure 9.
Comparison of the speed-size characteristics of different models at 
$\varGamma _{80}=3.8$, 
$\varGamma _{40}=0.6$ and 
${\rm \Delta} \phi =130^{\circ }$. In each panel, the black circles are experimental results of Valani et al. (Reference Valani, Slim and Simula2019), grey curves are the results using the model presented in this paper and the coloured curves are results from different models stated below with the colour indicating bouncing mode using the same conventions as in figure 5. Termination of the solid curves indicate coalescence. Results of using two alternative vertical spring models, a simple linear spring model and the logarithmic spring model of Moláček & Bush (Reference Moláček and Bush2013a), are shown in (a) and (b), respectively. Results obtained using a wave field from the model of Moláček & Bush (Reference Moláček and Bush2013b) and Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) are shown in (c) and (d), respectively. Results obtained by adding droplet deformation based on Blanchette (Reference Blanchette2016) and Gilet et al. (Reference Gilet, Terwagne, Vandewalle and Dorbolo2008) are shown in (e) and (f), respectively. For the grey curves and the coloured curves in all panels except (b), the linear spring model was used for the vertical dynamics with the parameters 
$K$ defined according to (2.14) and a fixed 
$B=0.60$.
Apart from the linear spring model used in this work, two alternative spring models for the vertical dynamics of a bouncing droplet were presented by Moláček & Bush (Reference Moláček and Bush2013a): (i) a simple linear spring model that does not restrict the normal force to be positive i.e. without the maximum condition in (2.2), and (ii) a logarithmic spring model, which can be implemented by replacing (2.1) with
\begin{equation} \left(1+\frac{C_3}{\ln^2{\left|\dfrac{C_1R}{\bar{z}_d}\right|}} \right)m\ddot{z}_d+ \frac{4}{3}\frac{{\rm \pi} \nu \rho R C_2}{\ln^2{\left|\dfrac{C_1R}{\bar{z}_d}\right|}} \dot{\bar{z}}_d+\dfrac{2{\rm \pi}\sigma}{\ln{\left|\dfrac{C_1R}{\bar{z}_d}\right|}}\bar{z}_d=-m[{g}+\gamma(t)], \end{equation}
when the droplet is in contact with the bath, and using 
$m\ddot {z}_d=-m[{g}+\gamma (t)]$ when the droplet is in the air. We fixed the parameter values to 
$C_1=2$, 
$C_2=12.5$ and 
$C_3=1.4$ which are typical values used for walkers (Moláček & Bush Reference Moláček and Bush2013a). Coupling these vertical dynamics models with the wave field and the horizontal dynamics described in § 2, we obtain the speed-size curves presented in figures 9(a) and 9(b).
Using the wave field of a walker from the Moláček & Bush (Reference Moláček and Bush2013b) model presented in (2.5) and the Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) model presented in (2.8) in place of the superwalker wave field that was used in this work, we obtain the speed-size curves shown in figures 9(c) and 9(d). These curves also show good match with the experiments on the ascending branch. We note that for a droplet in a 
$(1,2,1)^{{H}}$ bouncing mode, the subsequent bounce would occur one Faraday period after the initial impact. At this time, there is approximately a 10 % difference in the amplitudes between the three models, and a slightly greater difference in the gradients (see figure 2). This would suggest a comparable difference in the walking speeds. However, although in figure 9(c), the peak of the speed-size curve from the wave model of Moláček & Bush (Reference Moláček and Bush2013b) only goes up to approximately 
$17$ mm s
$^{-1}$ for the present choice of 
$K$ and 
$B$ values, we obtain a better fit to the experimental results by alternate choices of parameters 
$K$ and 
$B$. Hence, by tuning the 
$K$ and 
$B$ values and using the wave model of Moláček & Bush (Reference Moláček and Bush2013b), we can obtain good fit to the experimental data which is comparable to the fit obtained from the superwalker wave model. The speed-size curve from the wave model of Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) is identical to the curve from the superwalker wave model on the ascending branch. On the descending branch, we see that lower speeds are obtained from the Tadrist et al. (Reference Tadrist, Shim, Gilet and Schlagheck2018) model compared to the superwalker wave field. This shows that the added 
$20$ Hz waves seems to slightly speed up larger droplets on the descending branch in 
$(1,2,1)^{{L}}$ bouncing mode.
Finally, to account for droplet deformations, we couple the droplet deformation models of Blanchette (Reference Blanchette2016) and Gilet et al. (Reference Gilet, Terwagne, Vandewalle and Dorbolo2008) to the theoretical model presented in § 2. The additional droplet deformation equation for the model of Blanchette (Reference Blanchette2016) is
\begin{equation} m\ddot{R}_v+c_d\dot{R}_v+m\omega^2(R_v-R)=-F_N(t), \end{equation}
where 
$R_v$ is the vertical radius of the droplet, 
$c_d=3.8m \nu /R^2$ is the effective damping coefficient of the droplet deformation, and 
$\omega =\sqrt {N_{\omega }\sigma /\rho R^3}$ is the droplet's natural frequency with 
$N_{\omega }=5.84$. The model of Gilet et al. (Reference Gilet, Terwagne, Vandewalle and Dorbolo2008) after some algebra also reduces to an effectively similar equation for droplet deformations and is given by
\begin{equation} c_3m\ddot{R}_v+\frac{c_5\,m \nu}{R^2}\dot{R}_v+c_4\sigma(R_v-R)=-c_6 F_N(t), \end{equation}
where the parameters 
$c_3=0.1$, 
$c_4=10$, 
$c_5=3.3$ and 
$c_6=1$. While implementing both of these models, the criteria for contact changes from 
$\bar {z}_d\leq 0$ to 
$\bar {z}_d+R-R_v\leq 0$. Coupling these droplet deformation models to the theoretical model in § 2 results in the speed-size curves shown in figures 9(e) and 9(f). We see that the model of Gilet et al. (Reference Gilet, Terwagne, Vandewalle and Dorbolo2008) seems to have an insignificant effect on the speed-size characteristics with the curves completely overlapping each other. The model of Blanchette (Reference Blanchette2016) increases the walking speed of droplets in a small neighbourhood around 
$R=0.7$ mm but the model is still unable to capture the large superwalkers.
Appendix C. Determination of parameters $K$ and $B$
The theoretical model for simulating superwalkers presented in § 2 has three free parameters that are not known for superwalkers: (i) the dimensionless spring constant 
$K$ (ii) the dimensionless damping coefficient 
$B$ and (iii) the contact drag coefficient 
$C$. In our study we fixed 
$C=0.17$, a typical value that is used for walkers (Moláček & Bush Reference Moláček and Bush2013b). To determine values of 
$K$ and 
$B$, we simulated superwalkers in the 
$(K,B)$ parameter space and selected values that provide a good fit to the experimental results of Valani et al. (Reference Valani, Slim and Simula2019). We found that using constant values of 
$K=0.70$ and 
$B=0.60$ provided a reasonably good fit for small- to moderate-sized superwalkers on the ascending branch of the speed-size curves presented in figure 5, but failed for the largest superwalkers on the ascending branch for 
$\varGamma _{40}=1$. By allowing 
$K$ to vary linearly with the droplet radius 
$R$ while keeping 
$B$ fixed to 
$0.60$, we were able to obtain a better fit on the ascending branch for the results presented in figure 5. To arrive at this linear relationship, we simulated superwalkers for a fixed 
$\varGamma _{80}=3.8$, 
${\rm \Delta} \phi =130^{\circ }$ and four different values of 
$\varGamma _{40}=0$, 
$0.3$, 
$0.6$ and 
$1$. Droplet size that corresponds to the ascending branch in figure 5 were simulated. Typical graphs for the droplet speed and bouncing modes are shown in figure 10. For each droplet size and 
$\varGamma _{40}=0.6$ and 
$1$, the region of the 
$(K,B)$ parameter space where the relative difference between the speed of simulated superwalker and the corresponding experimental value of 
${\rm \Delta} u/u_{{exp}}=(u-u_{{exp}})/u_{{exp}}$ is within 
$20\,\%$ was determined and then a value of 
$K$ was selected from that region that matched with the experimentally observed bouncing mode. A linear best fit through all such 
$K$ values for different sized droplets results in one generic linear relationship given in (2.14).
Determination of the parameter values 
$K$ and 
$B$. Bouncing modes (markers) and relative difference between the numerical and the experimental values of the walking speed 
${\rm \Delta} u/u_{{exp}}=(u-u_{{exp}})/u_{{exp}}$ (contours) in the 
$(K,B)$ parameter space for three different droplet radii (a) 
$R=0.5$ mm, (b) 
$R=0.6$ mm and (c) 
$R=0.7$ mm at 
$\varGamma _{80}=3.8$, 
$\varGamma _{40}=1$ and 
${\rm \Delta} \phi =130^{\circ }$. In all the three panels, blue circles 
$\bullet$ are 
$(1,2,1)^{{H}}$, green triangles 
$\blacktriangle$ are 
$(1,2,1)^{{L}}$, yellow asterisks 
$*$ are (2,4,2) and purple squares 
$\blacksquare$ represent chaotic or other higher periodicity bouncing modes. The cyan solid lines represent the boundaries of the region inside which 
$|{\rm \Delta} u/u_{{exp}}|<20\,\%$ and the red circle corresponds to our chosen 
$K$ according to (2.14) and a fixed 
$B=0.60$.
