2.1 The Dirichlet-to-Neumann map

To numerically compute the DtN map in the corral geometry, we employ a domain decomposition method. In simple geometries, there are known analytical solutions to Laplace’s equation. In certain complex geometries that are composed of combinations of such simple geometries, the domain decomposition method is particularly effective since the known solutions form the building blocks of the DtN map. In our case, the corral is the union of two cylinders. By using a Hankel transform to solve Laplace’s equation analytically in each cylindrical sub-domain, and using appropriate matching, one obtains the required DtN map for the full domain. This method also applies to the annular cavity (Filoux et al. Reference Filoux, Hubert and Vandewalle2015) and submerged pillar (Harris et al. Reference Harris, Brun, Damiano, Faria and Bush2018) configurations.

Other problems can also be approached this way: when each sub-domain instead has a rectangular cross-section, such as for planar steps (Pucci et al. Reference Pucci, Sáenz, Faria and Bush2016), square corrals (Protière et al. Reference Protière, Boudaoud and Couder2006; Eddi et al. Reference Eddi, Fort, Moisy and Couder2009) or linear channels (Harris Reference Harris2015; Filoux et al. Reference Filoux, Hubert, Schlagheck and Vandewalle2017), a domain decomposition based on Fourier modes is appropriate. Furthermore, elliptical corrals (Sáenz et al. Reference Sáenz, Cristea-Platon and Bush2018) can be studied using elliptical coordinates (corresponding to Mathieu functions). The efficiency of the method can be further increased when each sub-geometry exhibits either rotational or translational invariance, as discussed below.

We first decompose the velocity potential, $\unicode[STIX]{x1D719}$, for each sub-domain

$$\begin{eqnarray}\unicode[STIX]{x1D719}=\left\{\begin{array}{@{}l@{}}\unicode[STIX]{x1D719}^{S},\quad (\boldsymbol{x},z)\in \unicode[STIX]{x1D6FA}_{S},\quad \\ \unicode[STIX]{x1D719}^{D},\quad (\boldsymbol{x},z)\in \unicode[STIX]{x1D6FA}_{D}.\quad \end{array}\right.\end{eqnarray}$$

Both $\unicode[STIX]{x1D719}^{S}$ and $\unicode[STIX]{x1D719}^{D}$ satisfy Laplace’s equation in their respective domains, equipped with boundary conditions consistent with (2.2). At each instant in time, the aim is to solve for $\unicode[STIX]{x1D719}$ given Dirichlet data $\unicode[STIX]{x1D719}|_{z=0}=a(\boldsymbol{x})$, from which we determine the corresponding Neumann data, $\unicode[STIX]{x1D719}_{z}|_{z=0}$. To solve for $\unicode[STIX]{x1D719}^{S}$ and $\unicode[STIX]{x1D719}^{D}$ concurrently, we introduce the unknown functions $b^{S}(\boldsymbol{x})$ and $b^{D}(\boldsymbol{x})$, satisfying

$$\begin{eqnarray}\unicode[STIX]{x1D719}^{S}=b^{S}\quad \text{and}\quad \unicode[STIX]{x1D719}^{D}=b^{D}\quad \text{on }z=-h_{s}.\end{eqnarray}$$

The functions $b^{S}(\boldsymbol{x})$ and $b^{D}(\boldsymbol{x})$ are chosen to be consistent with the fluid flow, specifically

$$\begin{eqnarray}\displaystyle & b^{S}=b^{D}\quad \text{for }|\boldsymbol{x}|<1, & \displaystyle \nonumber\\ \displaystyle & \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}^{S}=\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}^{D}\quad \text{on }z=-h_{s},|\boldsymbol{x}|<1, & \displaystyle \nonumber\\ \displaystyle & \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}^{S}=0\quad \text{on }z=-h_{s},|\boldsymbol{x}|\in [1,r_{\infty }). & \displaystyle \nonumber\end{eqnarray}$$

Upon solving for $b^{S}$, one can determine $\unicode[STIX]{x1D719}^{S}$ and hence $\unicode[STIX]{x1D719}_{z}^{S}|_{z=0}$, completing the DtN map.

The inversion for $b^{S}$ and $b^{D}$ can be simplified when the domain exhibits rotational or translational invariance, as is the case for circular corrals, infinitely long rectilinear channels and annular cavities. Upon appropriate discretisation, this invariance allows one to divide the domain decomposition method into a series of problems of much smaller size (for which computations may be parallelised), significantly reducing the computational cost. We proceed to describe this method for a circular corral, the domain of interest in this study. (The implementations of the DtN map for an annular channel and an infinite-length, rectilinear channel are presented in the supplementary material available at https://doi.org/10.1017/jfm.2020.140.)

To analytically solve for $\unicode[STIX]{x1D719}^{S}$ and $\unicode[STIX]{x1D719}^{D}$ in the circular corral geometry, we introduce cylindrical polar coordinates $(r,\unicode[STIX]{x1D703},z)$. We proceed by expanding $\unicode[STIX]{x1D719}^{S}$, $\unicode[STIX]{x1D719}^{D}$, $a$, $b^{S}$ and $b^{D}$ into Fourier angular modes, for example $\unicode[STIX]{x1D719}^{S}(r,\unicode[STIX]{x1D703},z)=\sum _{m=-\infty }^{\infty }\unicode[STIX]{x1D719}_{m}^{S}(r,z)\text{e}^{\text{i}m\unicode[STIX]{x1D703}}$, where $\text{i}^{2}=-1$. As $\unicode[STIX]{x1D719}^{S}$ and $\unicode[STIX]{x1D719}^{D}$ are both harmonic functions, it is necessary to solve

$$\begin{eqnarray}\displaystyle & {\mathcal{L}}_{m}\unicode[STIX]{x1D719}_{m}^{S}=0,\quad \text{for }(r,z)\in [0,r_{\infty })\times (-h_{s},0), & \displaystyle \nonumber\\ \displaystyle & {\mathcal{L}}_{m}\unicode[STIX]{x1D719}_{m}^{D}=0,\quad \text{for }(r,z)\in [0,1)\times (-h_{d},-h_{s}), & \displaystyle \nonumber\end{eqnarray}$$

where ${\mathcal{L}}_{m}=\unicode[STIX]{x2202}_{rr}+r^{-1}\unicode[STIX]{x2202}_{r}+\unicode[STIX]{x2202}_{zz}-(m/r)^{2}$ is Laplace’s operator for angular mode $m$. As depicted in figure 2, the boundary conditions for the shallow region are $\unicode[STIX]{x1D719}_{m}^{S}=0$ on $r=r_{\infty }$, $\unicode[STIX]{x1D719}_{m}^{S}=a_{m}$ on $z=0$, $\unicode[STIX]{x1D719}_{m}^{S}=b_{m}^{S}$ on $z=-h_{s}$, with $\unicode[STIX]{x1D719}_{m}^{S}$ bounded as $r\rightarrow 0$. Similarly, we require $\unicode[STIX]{x2202}_{r}\unicode[STIX]{x1D719}_{m}^{D}=0$ on $r=1$, $\unicode[STIX]{x1D719}_{m}^{D}=b_{m}^{D}$ on $z=-h_{s}$, $\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{m}^{D}=0$ on $z=-h_{d}$ and boundedness of $\unicode[STIX]{x1D719}_{m}^{D}$ as $r\rightarrow 0$.

Figure 2.

Cross-sectional schematic diagram of the domain decomposition method used to compute the DtN map. The circles indicate the $N_{S}$ radial collocation points, $\unicode[STIX]{x1D70C}_{j}$, with black circles for $j\leqslant N_{D}$ and white circles for $j>N_{D}$. The differential operator is ${\mathcal{L}}_{m}=\unicode[STIX]{x2202}_{rr}+r^{-1}\unicode[STIX]{x2202}_{r}+\unicode[STIX]{x2202}_{zz}-(m/r)^{2}$. At $r=0$, we seek a bounded solution.

We first solve for $\unicode[STIX]{x1D719}_{m}^{S}$ and $\unicode[STIX]{x1D719}_{m}^{D}$ before determining $b_{m}^{S}$ and $b_{m}^{D}$ according to the matching conditions. To solve ${\mathcal{L}}_{m}\unicode[STIX]{x1D719}_{m}^{S}=0$ with respect to the radial boundary conditions, we employ the Hankel transform $\unicode[STIX]{x1D719}_{m}^{S}(r,z)=\sum _{p=1}^{\infty }\unicode[STIX]{x1D719}_{mp}^{S}(z)\text{J}_{m}(k_{mp}^{S}r)$ (similarly for $b_{m}^{S}$ and $a_{m}$), where $\text{J}_{m}(x)$ is the Bessel function of the first kind of order $m$. The distinct wavenumbers $k_{mp}^{S}>0$ satisfy the Dirichlet boundary condition (2.2b) at $r=r_{\infty }$, specifically $\text{J}_{m}(k_{mp}^{S}r_{\infty })=0$. The analytic solution for $\unicode[STIX]{x1D719}_{m}^{S}(r,z)$ is then

(2.4)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{m}^{S}(r,z)=\mathop{\sum }_{p=1}^{\infty }\left[a_{mp}\frac{\sinh (k_{mp}^{S}(z+h_{s}))}{\sinh (k_{mp}^{S}h_{s})}-b_{mp}^{S}\frac{\sinh (k_{mp}^{S}z)}{\sinh (k_{mp}^{S}h_{s})}\right]\text{J}_{m}(k_{mp}^{S}r).\end{eqnarray}$$

In the deep region, the radial boundary conditions determine the Hankel transform $\unicode[STIX]{x1D719}_{m}^{D}(r,z)=\sum _{p=1}^{\infty }\unicode[STIX]{x1D719}_{mp}^{D}(z)\text{J}_{m}(k_{mp}^{D}r)$, where the distinct wavenumbers $k_{mp}^{D}\geqslant 0$ satisfy $\text{J}_{m}^{\prime }(k_{mp}^{D})=0$. This problem yields the well-known DtN map for finite depth, namely $\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{mp}^{D}(-h_{s})=b_{mp}^{D}k_{mp}^{D}\tanh (k_{mp}^{D}(h_{d}-h_{s}))$. Hence,

(2.5)$$\begin{eqnarray}\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{m}^{D}(r,-h_{s})=\mathop{\sum }_{p=1}^{\infty }b_{mp}^{D}k_{mp}^{D}\tanh (k_{mp}^{D}(h_{d}-h_{s}))\text{J}_{m}(k_{mp}^{D}r).\end{eqnarray}$$

To match the velocity potential along $z=-h_{s}$ for each angular mode $m$, we utilise $N_{S}$ radial collocation points $\unicode[STIX]{x1D70C}_{j}$, where $0\leqslant \unicode[STIX]{x1D70C}_{1}<\unicode[STIX]{x1D70C}_{2}<\cdots <\unicode[STIX]{x1D70C}_{N_{S}}<r_{\infty }$, which we define in appendix A. To satisfy the velocity continuity between the domains $\unicode[STIX]{x1D6FA}_{S}$ and $\unicode[STIX]{x1D6FA}_{D}$, as well as the no-flux condition through the base of the corral (see (2.2e)), at each radial collocation point, $\unicode[STIX]{x1D70C}_{j}$, we require

(2.6a)$$\begin{eqnarray}\displaystyle & b_{m}^{S}=b_{m}^{D},\quad \unicode[STIX]{x1D70C}_{j}\in [0,1), & \displaystyle\end{eqnarray}$$

(2.6b)$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{m}^{S}=\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{m}^{D},\quad \unicode[STIX]{x1D70C}_{j}\in [0,1),z=-h_{s}, & \displaystyle\end{eqnarray}$$

(2.6c)$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{m}^{S}=0,\quad \unicode[STIX]{x1D70C}_{j}\in [1,r_{\infty }),z=-h_{s}. & \displaystyle\end{eqnarray}$$

$N_{D}<N_{S}$

$\unicode[STIX]{x1D70C}_{j}<1$

$k_{mp}^{S}$

$k_{mp}^{D}$

$p\leqslant N_{S}$

$p\leqslant N_{D}$

To satisfy conditions (2.6b)–(2.6c), we evaluate (2.4)–(2.5) at the collocation points, yielding a system of $N_{S}$ equations for $b_{m1}^{S},\ldots ,b_{mN_{S}}^{S}$ and $b_{m1}^{D},\ldots ,b_{mN_{D}}^{D}$

(2.7)$$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{p=1}^{N_{S}}(a_{mp}-b_{mp}^{S}\cosh (k_{mp}^{S}h_{s}))\frac{k_{mp}^{S}\text{J}_{m}(k_{mp}^{S}\unicode[STIX]{x1D70C}_{j})}{\sinh (k_{mp}^{S}h_{s})}\nonumber\\ \displaystyle & & \displaystyle \quad =\left\{\begin{array}{@{}ll@{}}\displaystyle \mathop{\sum }_{p=1}^{N_{D}}b_{mp}^{D}k_{mp}^{D}\tanh (k_{mp}^{D}(h_{d}-h_{s}))\text{J}_{m}(k_{mp}^{D}\unicode[STIX]{x1D70C}_{j}),\quad & j\leqslant N_{D},\\ 0,\quad & j>N_{D}.\end{array}\right.\end{eqnarray}$$

To close the system, the continuity condition (2.6a) gives

(2.8)$$\begin{eqnarray}\mathop{\sum }_{p=1}^{N_{S}}b_{mp}^{S}\text{J}_{m}(k_{mp}^{S}\unicode[STIX]{x1D70C}_{j})=\mathop{\sum }_{p=1}^{N_{D}}b_{mp}^{D}\text{J}_{m}(k_{mp}^{D}\unicode[STIX]{x1D70C}_{j}),\quad \forall j\leqslant N_{D}.\end{eqnarray}$$

For each $m\in \mathbb{Z}$, equations (2.7)–(2.8) may be recast in the form

(2.9)$$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D63D} & \unicode[STIX]{x1D63E}\\ \unicode[STIX]{x1D63F} & \unicode[STIX]{x1D640}\end{array}\right)\left(\begin{array}{@{}c@{}}\boldsymbol{b}_{S}\\ \boldsymbol{b}_{D}\end{array}\right)=\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D63C}\boldsymbol{a}\\ \boldsymbol{ 0}\end{array}\right),\end{eqnarray}$$

where it remains to find vectors $\boldsymbol{b}_{S}$ and $\boldsymbol{b}_{D}$ for given $\boldsymbol{a}$. Following the numerical computation of $b_{mp}^{S}$ and $b_{mp}^{D}$, equation (2.4) determines the DtN map

$$\begin{eqnarray}\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{mp}^{S}(0)=k_{mp}^{S}(a_{mp}\coth (k_{mp}^{S}h_{s})-b_{mp}^{S}\,\text{cosech}(k_{mp}^{S}h_{s})).\end{eqnarray}$$

When repeatedly computing the DtN, we pre-compute the LU-decomposition (with partial pivoting) of the map’s mass matrix (see (2.9)).

The justification for truncating large wavenumbers relies on the shallow penetration depth of short waves, in which surface perturbations do not propagate to the sufficiently deep submerged topography. Although energy may be transferred from long to short waves across the step, the short waves continue to evolve as if the domain were of infinite depth and so do not significantly influence the wave evolution near the submerged step.

Figure 3.

The DtN map for $m=1$ with a step at $r=1$ (black dashed line), outer radius $r_{\infty }=5$ and depths $h_{d}=1$ and $h_{s}=0.1$. The submerged topography is not shown in the figure. (a) Example Dirichlet data $\unicode[STIX]{x1D719}_{1}^{S}(r,0)=\text{J}_{m}(k_{1,8}^{S}r)$. (b) The Neumann data, $\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{1}^{S}(r,0)$, for the corral geometry (black curve) computed using $N_{S}=512$ collocation points, $\unicode[STIX]{x1D70C}_{j}$. The corresponding data $\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{1}^{S}(r,0)$ for a cylindrical domain of radius $r_{\infty }$ with constant depth $h_{s}$ (light grey curve) and constant depth $h_{d}$ (dash-dotted curve). (c) Zoom-in of the rectangular box in panel $(b)$, showing that, for the corral, $\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}_{1}^{S}(r,0)$ coincides with the constant-depth cases away from the submerged step, with a smooth transition between the two solutions across $r=1$.

We present an example of the DtN map for the corral geometry in figure 3 for $m=1$ and given Dirichlet data. Sufficiently far away from the submerged step, the Neumann data for the corral coincides with Neumann data for a cylindrical domain of constant depth $h$, with $h=h_{d}$ for $r<1$, and $h=h_{s}$ for $r>1$. In the vicinity of the step, the corral solution smoothly transitions between the two constant-depth solutions, over a length scale determined by the spatial non-locality of the operator. The self-convergence of this method as $N_{S}\rightarrow \infty$ is discussed in § B.1, from which we deduce that $N_{S}=512$ is sufficient to accurately resolve the DtN map in the corral geometries considered.

Figure 4.

Comparison of the wave field generated by a single impact, evolved using our model and Faria’s model, for corral radius $R_{c}=14~\text{mm}$, frequency $f_{0}=80~\text{Hz}$ and $\unicode[STIX]{x1D6E4}=4.1$ (corresponding to $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.9$ in our model). The free-surface elevation is normalised by the droplet radius, $R_{0}=0.35~\text{mm}$. The impact was at position $\boldsymbol{x}/R_{c}=(0.5,0)^{\text{T}}$, indicated by the cross in (a,b), with phase $\unicode[STIX]{x1D703}_{I}=1.43\unicode[STIX]{x03C0}$. (a,b) The computations of our model are presented for $y>0$, and Faria’s model for $y<0$, at times (a) $t/T_{F}=4$ and (b) $t/T_{F}=8$. The white circle denotes the submerged step. (c,d) The $y=0$ cross-section of our model (black) and Faria’s model (grey). The data in (c,d) are identical to (a,b).

Unlike the work of Faria (Reference Faria2017), our computation of the DtN map does not invoke jumps in the velocity potential across the step, making this approach applicable to inviscid flows, without relying on viscosity to dampen spurious short waves. Moreover, Faria’s method is specifically designed to approximate the evolution of Faraday waves, whereas our formulation may be applied to a broad class of water wave problems. We explore the effect of these approximations in § 2.2.

2.2 Comparison to a Faraday phase-speed-matched model

We compare the predictions of our model, with the computation of the DtN map described in § 2.1, to that of Faria (Reference Faria2017). In Faria’s model, the DtN map is replaced by the expression $\unicode[STIX]{x1D719}_{z}|_{z=0}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(b\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|_{z=0})$, where $b(\boldsymbol{x})$ is chosen to match the phase speed of the Faraday wavenumber in regions of constant depth. In dimensional variables, $b(\boldsymbol{x})$ takes the value $\tanh (k_{D}H_{D})/k_{D}$ in the deep region and $\tanh (k_{S}H_{S})/k_{S}$ in the shallow layer, where, for given $\unicode[STIX]{x1D6E4}$, $k_{D}$ and $k_{S}$ are the critical Faraday wavenumbers for flat-bottom domains of depth $H_{D}$ and $H_{S}$, respectively (Milewski et al. Reference Milewski, Galeano-Rios, Nachbin and Bush2015; Faria Reference Faria2017). Other than the DtN map, the two models are identical, but their implementation differs significantly.

We consider the wave evolution arising from a single instantaneous point impact, which is the Green’s function for the corral geometry. The numerical implementation of our model is described in § 2.3, and Faria’s model is implemented using the Fourier pseudo-spectral method detailed in Faria (Reference Faria2017). As highlighted in figure 4, there are marked differences in the wave evolution, particularly at short times. To understand the origins of this discrepancy, we consider the same wave field generation but for a flat-bottom domain (with depth $H_{D}$ or $H_{S}$) and infinite horizontal extent, the implementation of which is as described by Durey & Milewski (Reference Durey and Milewski2017). As presented in figure 5, similar discrepancies persist in the absence of submerged topography at depth $H_{D}$, but less so at depth $H_{S}$ (since Faria’s DtN map is inspired by shallow water theory), highlighting a fundamental difference in the wave evolution of these two models in the deep region. This difference is significant in the resulting pilot-wave dynamics, as the droplet is confined to the deep region of the corral, for which the wave field discrepancy is greatest.

Figure 5.

Comparison of the axisymmetric wave field generated by a single impact over a flat-bottom domain with infinite horizontal extent, evolved using our model (black curves) and Faria’s model (grey curves) for $f_{0}=80~\text{Hz}$, $\unicode[STIX]{x1D6E4}=3.8$ and $\unicode[STIX]{x1D703}_{I}=1.43\unicode[STIX]{x03C0}$. The free surface is normalised by the droplet radius, $R_{0}=0.35~\text{mm}$. (a) Depth $H_{D}=6.6~\text{mm}$. (b) Depth $H_{S}=0.6~\text{mm}$. (c) Phase speed of the unforced waves ($\unicode[STIX]{x1D6E4}=0$) for our model ($C_{k}$, black curves) and Faria’s model ($C_{k}^{F}$, grey curves). Solid and dash-dotted curves correspond to depths $H_{D}$ and $H_{S}$, respectively, with respective intersections at $k=k_{D}$ and $k=k_{S}$.

The differences at short time may be explained through consideration of the phase speed, $C_{k}$, arising from the dispersion relation of (2.1) for wavenumber $k$ and $\unicode[STIX]{x1D6E4}=0$. Similarly, $C_{k}^{F}$ is the phase speed for Faria’s model. For constant depth $H$, we obtain

$$\begin{eqnarray}C_{k}=\sqrt{\frac{\tanh (kH)}{k}\left(g+\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}k^{2}\right)}\quad \text{and}\quad C_{k}^{F}=\sqrt{\frac{\tanh (k_{F}H)}{k_{F}}\left(g+\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}k^{2}\right)},\end{eqnarray}$$

where $k_{F}$ is the corresponding critical Faraday wavenumber ($k_{D}$ or $k_{S}$). As presented in figure 5(c), $C_{k}$ and $C_{k}^{F}$ agree only at $k=k_{F}$, where the difference is most apparent for deep fluid. Waves of intermediate length are only weakly damped at short times, playing a pivotal role in the wave field evolution. While the discrepancy between the two models may reduce when one considers subharmonic periodic impacts on the bath resonating with the Faraday frequency (Faria Reference Faria2017), short-time transient effects can strongly influence the pilot-wave dynamics of droplet–droplet interactions (Galeano-Rios, Couchman & Bush Reference Galeano-Rios, Couchman and Bush2018) or for dynamics very close to boundaries, where reflected transients are significant. Our method for computing the DtN map is thus particularly beneficial for modelling droplet–droplet dynamics confined by submerged topography, such as strings of propelling droplets (Filoux et al. Reference Filoux, Hubert and Vandewalle2015) and droplet lattices (Thomson, Couchman & Bush Reference Thomson, Couchman and Bush2020).

2.3 The discrete-time iterative map

We employ a similar decomposition to that of the previous section, expressing $\boldsymbol{x}$ as $(r,\unicode[STIX]{x1D703})$ in plane-polar coordinates. Hence, $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D719}|_{z=0}$ are recast as

$$\begin{eqnarray}\unicode[STIX]{x1D702}(\boldsymbol{x},t)=\mathop{\sum }_{m=-\infty }^{\infty }\mathop{\sum }_{p=1}^{\infty }\unicode[STIX]{x1D702}_{mp}(t)\unicode[STIX]{x1D6F7}_{mp}(\boldsymbol{x})\quad \text{and}\quad \unicode[STIX]{x1D719}(\boldsymbol{x},0,t)=\mathop{\sum }_{m=-\infty }^{\infty }\mathop{\sum }_{p=1}^{\infty }\unicode[STIX]{x1D719}_{mp}(t)\unicode[STIX]{x1D6F7}_{mp}(\boldsymbol{x}),\end{eqnarray}$$

with orthogonal basis functions $\unicode[STIX]{x1D6F7}_{mp}(\boldsymbol{x})=\text{J}_{m}(k_{mp}^{S}r)\text{e}^{\text{i}m\unicode[STIX]{x1D703}}$. For tractable notation, we henceforth remove the superscript from $k_{mp}^{S}$. The reality condition for $\unicode[STIX]{x1D702}$ is $\unicode[STIX]{x1D702}_{-m,p}=(-1)^{m}\unicode[STIX]{x1D702}_{m,p}^{\ast }$, where $\ast$ denotes complex conjugation, and similarly for $\unicode[STIX]{x1D719}$.

By similarly expressing the pressure, $P_{0}(\boldsymbol{x},t)$, in terms of the transformed variables $P_{mp}(t)$, the dimensionless governing fluid equations (2.1) reduce to the system

(2.10a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D719}_{mp}=-G(1-\unicode[STIX]{x1D6E4}\cos (4\unicode[STIX]{x03C0}t))\unicode[STIX]{x1D702}_{mp}-Bo\,k_{mp}^{2}\unicode[STIX]{x1D702}_{mp}-2Re^{-1}k_{mp}^{2}\unicode[STIX]{x1D719}_{mp}-P_{mp}, & \displaystyle\end{eqnarray}$$

(2.10b)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D702}_{mp}={\mathcal{N}}_{mp}[\unicode[STIX]{x1D753}_{m}]-2Re^{-1}k_{mp}^{2}\unicode[STIX]{x1D702}_{mp}, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D753}_{m}=(\unicode[STIX]{x1D719}_{m1},\unicode[STIX]{x1D719}_{m2},\ldots )^{\text{T}}$

${\mathcal{N}}_{mp}[\unicode[STIX]{x1D753}_{m}]$

$m$

$k_{mp}$

The Hankel transform of a point-force pressure distribution is

(2.11)$$\begin{eqnarray}P_{mp}(t)=\frac{f(t)}{2\unicode[STIX]{x03C0}C_{mp}}\unicode[STIX]{x1D6F7}_{mp}^{\ast }(\boldsymbol{X}(t)),\end{eqnarray}$$

where $C_{mp}=\int _{0}^{r_{\infty }}r\text{J}_{m}^{2}(k_{mp}r)\,\text{d}r=\frac{1}{2}(r_{\infty }\text{J}_{m}^{\prime }(k_{mp}r_{\infty }))^{2}$. Integrating the dynamic boundary condition (2.10a) across an impact at time $t_{n}$ yields the jump

(2.12)$$\begin{eqnarray}[\unicode[STIX]{x1D719}_{mp}(t_{n})]_{-}^{+}=-\frac{1}{2\unicode[STIX]{x03C0}C_{mp}}\unicode[STIX]{x1D6F7}_{mp}^{\ast }(\boldsymbol{X}(t_{n})),\end{eqnarray}$$

while during flight ($t\neq t_{n}$), system (2.10) simplifies to

(2.13a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D719}_{mp}=-G(1-\unicode[STIX]{x1D6E4}\cos (4\unicode[STIX]{x03C0}t))\unicode[STIX]{x1D702}_{mp}-Bo\,k_{mp}^{2}\unicode[STIX]{x1D702}_{mp}-2Re^{-1}k_{mp}^{2}\unicode[STIX]{x1D719}_{mp}, & \displaystyle\end{eqnarray}$$

(2.13b)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D702}_{mp}={\mathcal{N}}_{mp}[\unicode[STIX]{x1D753}_{m}]-2Re^{-1}k_{mp}^{2}\unicode[STIX]{x1D702}_{mp}. & \displaystyle\end{eqnarray}$$

As the governing fluid equations (2.13) are of Floquet form and the droplet impacts are modelled as instantaneous, we may evolve $\unicode[STIX]{x1D702}_{mp}$ and $\unicode[STIX]{x1D719}_{mp}$ between impacts ($t_{n}^{+}\mapsto t_{n+1}^{-}$) using fundamental matrices, $\unicode[STIX]{x1D648}_{m}(\unicode[STIX]{x1D6E4})$, for each angular mode, $m\in \mathbb{Z}$. We exploit periodicity and initialise (2.13) with $\unicode[STIX]{x1D702}_{mp}(t_{0})=(1-\unicode[STIX]{x1D70D})\unicode[STIX]{x1D6FF}_{pq}$ and $\unicode[STIX]{x1D719}_{mp}(t_{0})=\unicode[STIX]{x1D70D}\unicode[STIX]{x1D6FF}_{pq}$ for each $q\in \{1,\ldots ,P\}$ and $\unicode[STIX]{x1D70D}\in \{0,1\}$, where $P\leqslant N_{S}$ accounts for a further truncation of the highly damped short waves, which have a negligible effect on the pilot-wave dynamics (but do play a role in the computation of the DtN map). Suitable values of $P\leqslant N_{S}$ are discussed in § B.2. We then solve the governing equations (2.13) numerically for a given vibrational forcing $\unicode[STIX]{x1D6E4}<\unicode[STIX]{x1D6E4}_{F}$ to find $\unicode[STIX]{x1D702}_{mp}(t_{1})$ and $\unicode[STIX]{x1D719}_{mp}(t_{1})$ (for all $p=1,\ldots ,P$), corresponding to the entries of $\unicode[STIX]{x1D648}_{m}\in \mathbb{R}^{2P\times 2P}$ for column number $(q+P\unicode[STIX]{x1D70D})$. The numerical integration is performed using a fourth-order Runge–Kutta method with the integrating factor $\exp (2Re^{-1}k_{mp}^{2}(t-\unicode[STIX]{x1D70F}_{j}))$ for integrating from time $t=\unicode[STIX]{x1D70F}_{j}$ to $t=\unicode[STIX]{x1D70F}_{j+1}$.

When truncating to $P\leqslant N_{S}$ wavenumbers, we define vectors $\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{n},\unicode[STIX]{x1D753}_{m}^{n}\in \mathbb{C}^{P\times 1}$ as

$$\begin{eqnarray}\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{n}=(\unicode[STIX]{x1D702}_{m1}(t_{n}),\ldots ,\unicode[STIX]{x1D702}_{mP}(t_{n}))^{\text{T}}\quad \text{and}\quad \unicode[STIX]{x1D753}_{m}^{n}=(\unicode[STIX]{x1D719}_{m1}(t_{n}^{+}),\ldots ,\unicode[STIX]{x1D719}_{mP}(t_{n}^{+}))^{\text{T}},\end{eqnarray}$$

where $\unicode[STIX]{x1D702}$ is continuous across each impact, whereas we consider the velocity potential, $\unicode[STIX]{x1D719}$, in the instant just after impact. We also define $\unicode[STIX]{x1D64B}_{m}\in \mathbb{R}^{P\times P}$ as a diagonal matrix with elements $(2\unicode[STIX]{x03C0}C_{mp})^{-1}$ (the jump condition coefficients), and the row vector of basis functions $\unicode[STIX]{x1D731}_{m}(\boldsymbol{x})=(\unicode[STIX]{x1D6F7}_{m1}(\boldsymbol{x}),\ldots ,\unicode[STIX]{x1D6F7}_{mP}(\boldsymbol{x}))$. Combining the wave field evolution during flight (2.13) and at impact (2.12) yields (2.16b) in the iterative map below.

Under the imposed choice of $f(t)=\sum _{n=0}^{\infty }\unicode[STIX]{x1D6FF}(t-t_{n})$, equation (2.3) may be decomposed into the droplet’s horizontal motion during flight ($t\neq t_{n}$)

(2.14)$$\begin{eqnarray}\boldsymbol{X}^{\prime \prime }(t)+C_{a}\boldsymbol{X}^{\prime }(t)=\mathbf{0},\end{eqnarray}$$

with jumps in the droplet’s horizontal velocity at impact (Durey & Milewski Reference Durey and Milewski2017)

(2.15)$$\begin{eqnarray}[\boldsymbol{X}^{\prime }(t_{n})]_{-}^{+}=-(1-\text{e}^{-C_{i}})\left(\frac{\unicode[STIX]{x1D6FF}_{0}G}{C_{i}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}(\boldsymbol{X}(t_{n}),t_{n})+\boldsymbol{X}^{\prime }(t_{n}^{-})\right).\end{eqnarray}$$

To express the droplet dynamics (2.14)–(2.15) as a discrete-time map, we integrate (2.14) analytically and eliminate $\boldsymbol{X}^{\prime }(t_{n}^{\pm })$ in (2.15) in favour of $\boldsymbol{X}(t_{n})$ and $\boldsymbol{X}(t_{n+1})$. Then, by defining $\boldsymbol{X}_{n}=\boldsymbol{X}(t_{n})$ as the droplet position at impact and $\boldsymbol{U}_{n}=\boldsymbol{X}_{n+1}-\boldsymbol{X}_{n}$ as the mean velocity during flight, the iterative map for the corral pilot-wave dynamics is

(2.16a)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{X}_{n}-\boldsymbol{X}_{n-1}=\boldsymbol{U}_{n-1}, & \displaystyle\end{eqnarray}$$

(2.16b)$$\begin{eqnarray}\displaystyle & \displaystyle \left(\begin{array}{@{}c@{}}\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{n}\\ \unicode[STIX]{x1D753}_{m}^{n}\end{array}\right)=\unicode[STIX]{x1D648}_{m}(\unicode[STIX]{x1D6E4})\left(\begin{array}{@{}c@{}}\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{n-1}\\ \unicode[STIX]{x1D753}_{m}^{n-1}\end{array}\right)-\left(\begin{array}{@{}c@{}}\mathbf{0}\\ \unicode[STIX]{x1D64B}_{m}\unicode[STIX]{x1D731}_{m}^{\ast }(\boldsymbol{X}_{n})\end{array}\right), & \displaystyle\end{eqnarray}$$

(2.16c)$$\begin{eqnarray}\displaystyle & \displaystyle (\boldsymbol{U}_{n}-\boldsymbol{U}_{n-1})+{\mathcal{D}}\boldsymbol{U}_{n-1}=-{\mathcal{F}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}(\boldsymbol{X}_{n},t_{n}), & \displaystyle\end{eqnarray}$$

$\ast$

$\unicode[STIX]{x1D702}(\boldsymbol{x},t_{n})=\sum _{m=-\infty }^{\infty }\unicode[STIX]{x1D731}_{m}(\boldsymbol{x})\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{n}$

${\mathcal{D}}$

${\mathcal{F}}$

$$\begin{eqnarray}{\mathcal{D}}=1-\text{e}^{-(C_{i}+C_{a})}\quad \text{and}\quad {\mathcal{F}}=\frac{1-\text{e}^{-C_{i}}}{C_{i}}\frac{1-\text{e}^{-C_{a}}}{C_{a}}\unicode[STIX]{x1D6FF}_{0}G.\end{eqnarray}$$

MATLAB code for evolving (2.16) is provided in the supplementary material.

2.4 Summary

In § 2, we derived a first-principles fluid model for the corral geometry, which is coupled to a reduced model for the droplet’s evolution. To recast this pilot-wave system into a tractable form, we formulated a method for numerically computing the DtN map for the corral geometry, reducing the evolution of the system to flow on the surface. As discussed in §§ 2.1 and 2.2, this methodology supersedes the ad hoc formulation of Faria (Reference Faria2017) in many respects. We then utilised fundamental matrices to express the evolution of the system as a discrete-time iterative map, where all the information of the corral geometry, vibrational forcing and fluid properties is encoded within the fundamental matrices $\unicode[STIX]{x1D648}_{m}(\unicode[STIX]{x1D6E4})$. This discrete-time formulation is not only computationally more efficient (compared to repeatedly evolving the waves during flight (Faria Reference Faria2017)), but also provides a framework amenable for detailed analysis of the pilot-wave dynamics.

3.1 The Faraday threshold

In addition to computational efficiency, the fundamental matrices, $\unicode[STIX]{x1D648}_{m}(\unicode[STIX]{x1D6E4})$, also determine the subharmonic Faraday eigenmodes for given vibrational forcing $\unicode[STIX]{x1D6E4}<\unicode[STIX]{x1D6E4}_{F}$. Furthermore, the eigenvalues of $\unicode[STIX]{x1D648}_{m}$ give the decay rates of each fluid eigenmode. For each $m\in \mathbb{Z}$, we denote the dominant real eigenvalue as $\unicode[STIX]{x1D706}_{m}^{c}(\unicode[STIX]{x1D6E4})$, which satisfies

(3.1)$$\begin{eqnarray}\unicode[STIX]{x1D648}_{m}\left(\begin{array}{@{}c@{}}\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{c}\\ \unicode[STIX]{x1D753}_{m}^{c}\end{array}\right)=\unicode[STIX]{x1D706}_{m}^{c}\left(\begin{array}{@{}c@{}}\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{c}\\ \unicode[STIX]{x1D753}_{m}^{c}\end{array}\right).\end{eqnarray}$$

Each angular mode destabilises at the critical vibrational forcing $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{m}^{c}$, where $\unicode[STIX]{x1D706}_{m}^{c}(\unicode[STIX]{x1D6E4}_{m}^{c})=1$. This feature contrasts with earlier works (Durey & Milewski Reference Durey and Milewski2017), where the Faraday threshold was independent of the angular mode, $m$, for a fluid of infinite depth and horizontal extent. The corresponding Faraday threshold, $\unicode[STIX]{x1D6E4}_{F}$, for the corral geometry is thus the minimum for each angular mode, with $\unicode[STIX]{x1D6E4}_{F}=\min _{m}\unicode[STIX]{x1D6E4}_{m}^{c}.$ For $\unicode[STIX]{x1D6E4}<\unicode[STIX]{x1D6E4}_{m}^{c}$ and each angular mode $m$, the critical eigenvalue, $\unicode[STIX]{x1D706}_{m}^{c}(\unicode[STIX]{x1D6E4})$, determines the memory ${\mathcal{M}}_{m}(\unicode[STIX]{x1D6E4})=|\text{log}\,\unicode[STIX]{x1D706}_{m}^{c}(\unicode[STIX]{x1D6E4})|^{-1}$, a proxy for the number of prior impacts that significantly influence the system. The relative magnitude of each angular memory, ${\mathcal{M}}_{m}$, controls the form of the pilot-wave dynamics as $\unicode[STIX]{x1D6E4}\rightarrow \unicode[STIX]{x1D6E4}_{F}$, a concept revisited throughout this work. Moreover, the system memory, $Me(\unicode[STIX]{x1D6E4})=\max _{m}{\mathcal{M}}_{m}(\unicode[STIX]{x1D6E4})$, diverges in the limit $\unicode[STIX]{x1D6E4}\rightarrow \unicode[STIX]{x1D6E4}_{F}$.

3.2 Cavity modes

At the critical vibrational forcing $\unicode[STIX]{x1D6E4}_{m}^{c}$ for each angular mode $m$, the eigenvalue problem (3.1) yields the Faraday cavity modes $\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{c}$ and $\unicode[STIX]{x1D753}_{m}^{c}$. Inverting the Hankel transform back to physical space gives the corresponding radial part of each free-surface cavity mode, namely $\unicode[STIX]{x1D702}_{m}^{c}(r,t_{0})$ and $\unicode[STIX]{x1D719}_{m}^{c}(r,t_{0})$. It should be noted that, unlike the free-space dynamics explored in earlier works (Durey & Milewski Reference Durey and Milewski2017), the cavity modes are a superposition of Bessel functions with different wavenumbers. In particular, beyond the submerged step ($r\geqslant 1$), we observe exponential spatial decay in $\unicode[STIX]{x1D702}_{m}^{c}$ and $\unicode[STIX]{x1D719}_{m}^{c}$.

Figure 6.

The critical vibrational forcing, $\unicode[STIX]{x1D6E4}_{m}^{c}$, for $m=0$ (solid curve) and $m=1$ (dash-dotted curve), where the Faraday threshold is $\unicode[STIX]{x1D6E4}_{F}=\min _{m}\unicode[STIX]{x1D6E4}_{m}^{c}$. When $\unicode[STIX]{x1D6E4}_{0}^{c}<\unicode[STIX]{x1D6E4}_{1}^{c}$ (white background), the critical Faraday mode is axisymmetric, but has angular dependence for $\unicode[STIX]{x1D6E4}_{0}^{c}>\unicode[STIX]{x1D6E4}_{1}^{c}$ (grey background). (a) The dependence of $\unicode[STIX]{x1D6E4}_{m}^{c}$ on the vibrational frequency, $f_{0}$, for corral radius $R_{c}=10~\text{mm}$. The insets portray the real part of the critical cavity mode $\unicode[STIX]{x1D702}_{m}^{c}(r,0)\text{e}^{\text{i}m\unicode[STIX]{x1D703}}$ at $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{m}^{c}=\unicode[STIX]{x1D6E4}_{F}$ for frequencies $f_{0}=68~\text{Hz}$, $f_{0}=80~\text{Hz}$ and $f_{0}=92~\text{Hz}$ (arbitrary vertical units). The white circle indicates the submerged step. (b) The dependence of $\unicode[STIX]{x1D6E4}_{m}^{c}$ on the corral radius, $R_{c}$, for frequency $f_{0}=80~\text{Hz}$. The marker indicates half the Faraday wavelength, $\unicode[STIX]{x1D706}_{F}=5.02~\text{mm}$, for a circular cylinder of constant depth $H_{D}=6.6~\text{mm}$ (Milewski et al. Reference Milewski, Galeano-Rios, Nachbin and Bush2015). (c) The increase in critical forcing, $\unicode[STIX]{x0394}\unicode[STIX]{x1D6E4}_{m}^{c}=\unicode[STIX]{x1D6E4}_{m}^{c}-\unicode[STIX]{x1D6E4}_{F}$, for corral radius $R_{c}=11~\text{mm}$ with frequency $f_{0}=80~\text{Hz}$, where $\unicode[STIX]{x1D6E4}_{F}=\unicode[STIX]{x1D6E4}_{0}^{c}$. The dashed lines are guides for the eyes.

Figure 7.

Phase-space time evolution of the dimensionless cavity modes $\unicode[STIX]{x1D702}_{m}^{c}$ (black) and $\unicode[STIX]{x1D719}_{m}^{c}$ (grey) at $r=1$ (in arbitrary units). (a) $m=0$ and $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{0}^{c}=4.806$, (b) $m=1$ and $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{1}^{c}=4.754$. If a Dirichlet or Neumann approximation were adequate to describe the dynamics at $r=1$ then the trajectories would be approximately vertical or horizontal, respectively. The arrow denotes the time evolution. The reference dots denote the functions evaluated at time $t=0$. The frequency is $f_{0}=80~\text{Hz}$ with corral radius $R_{c}=10~\text{mm}$.

In figure 6, we present the critical vibrational forcing, $\unicode[STIX]{x1D6E4}_{m}^{c}$, as a function of vibrational frequency, $f_{0}$, and corral radius, $R_{c}$. Not only does changing the forcing frequency affect the Faraday wavelength, $\unicode[STIX]{x1D706}_{F}$, and the subharmonic Faraday period, $T_{F}=2/f_{0}$, but it also determines which angular mode, $m$, is the most unstable, i.e. the value of $m$ that minimises the critical forcing acceleration, $\unicode[STIX]{x1D6E4}_{m}^{c}$. For $73.5\lesssim f_{0}\lesssim 86.5~\text{Hz}$, the $m=1$ mode is excited beyond the critical forcing acceleration, whereas there exist surrounding regions where the critical mode is radially symmetric ($m=0$). In the parameter regimes considered, we observed that $\unicode[STIX]{x1D6E4}_{2m}^{c}\geqslant \unicode[STIX]{x1D6E4}_{2m-2}^{c}$ and $\unicode[STIX]{x1D6E4}_{2m+1}^{c}\geqslant \unicode[STIX]{x1D6E4}_{2m-1}^{c}$ for all $m\geqslant 1$ (see figure 6c), indicating that the Faraday instability, $\unicode[STIX]{x1D6E4}_{F}$, corresponds only to cavity modes with radial symmetry ($m=0$) or the lowest angular frequency $(m=1)$. This result differs from the Faraday instability in a circular cylinder, where cavity modes with larger angular frequencies may be the first to destabilise (Gilet Reference Gilet2016). We note, however, that there is no guarantee that the supercritical Faraday pattern predicted by the linear theory will agree with experiments, for which weakly nonlinear (Skeldon & Rucklidge Reference Skeldon and Rucklidge2015) and fully nonlinear (Périnet, Juric & Tuckerman Reference Périnet, Juric and Tuckerman2009; Kahouadji et al. Reference Kahouadji, Périnet, Tuckerman, Shin, Chergui and Juric2015) effects will be significant. Nevertheless, when only one angular mode is linearly unstable, we may expect that the nonlinear Faraday pattern has a similar form for a supercritical bifurcation at the Faraday threshold. Finally, it is important to note that, for subcritical vibrational forcing ($\unicode[STIX]{x1D6E4}<\unicode[STIX]{x1D6E4}_{F}$), we do not observe any consistent ordering of the memories ${\mathcal{M}}_{m}$.

To explore the feasibility of a simple effective boundary condition above the step, we consider the time evolution of the amplitude and gradient of the cavity modes, $\unicode[STIX]{x1D702}_{m}^{c}(r,t)$ and $\unicode[STIX]{x1D719}_{m}^{c}(r,t)$, evaluated at $r=1$, with the results presented in figure 7. At the critical vibrational forcing, $\unicode[STIX]{x1D6E4}_{m}^{c}$, the time evolution of the cavity modes is necessarily periodic, whereas the phase-space trajectory spirals towards the origin on the memory time scale, ${\mathcal{M}}_{m}$, for subcritical vibrations. Since the evolution of the free surface and velocity potential does not yield straight lines in phase space, considering effective Dirichlet (Blanchette Reference Blanchette2016; Filoux et al. Reference Filoux, Hubert, Schlagheck and Vandewalle2017), Neumann (Dubertrand et al. Reference Dubertrand, Hubert, Schlagheck, Vandewalle, Bastin and Martin2016; Gilet Reference Gilet2016) or Robin boundary conditions for $\unicode[STIX]{x1D702}_{m}^{c}$ and $\unicode[STIX]{x1D719}_{m}^{c}$ above the submerged step would be inadequate for capturing the evolution of the cavity modes. This result is an important feature of our study and highlights the benefits of our approach. A similar result was observed for the reflection of a walker off a submerged planar wall (Damiano et al. Reference Damiano, Brun, Harris, Galeano-Rios and Bush2016; Faria Reference Faria2017).

3.3 Relating the droplet statistics and the mean pilot wave

In certain circumstances, the mean pilot wave, averaged over the droplet’s trajectory, has a controlling influence on the droplet motion, determining the quantisation length of periodic trajectories (Labousse et al. Reference Labousse, Perrard, Couder and Fort2014; Perrard et al. Reference Perrard, Labousse, Miskin, Fort and Couder2014b) or acting as an effective potential in the long-path-memory limit (Cristea-Platon et al. Reference Cristea-Platon, Sáenz and Bush2018; Durey, Milewski & Bush Reference Durey, Milewski and Bush2018). Indeed, in an experimental study of elliptical corrals, Sáenz et al. (Reference Sáenz, Cristea-Platon and Bush2018) observed that the droplet’s statistical distribution and the accompanying mean pilot wave have a similar form, a result that we precisely define for the pilot-wave model considered herein.

As discussed by Durey et al. (Reference Durey, Milewski and Bush2018), several experimental and numerical studies provide evidence to support the validity of the ergodic hypothesis for confined, chaotic pilot-wave dynamics, typically leading to the emergence of a stationary probability density function, $\unicode[STIX]{x1D707}(\boldsymbol{x})$, for the droplet’s position. Under the assumptions of ergodicity and stationarity, we prove in § C.1 that the mean wave field

$$\begin{eqnarray}\bar{\unicode[STIX]{x1D702}}(\boldsymbol{x})=\lim _{N\rightarrow \infty }\frac{1}{N}\mathop{\sum }_{n=1}^{N}\unicode[STIX]{x1D702}(\boldsymbol{x},t_{n})\end{eqnarray}$$

satisfies

(3.2)$$\begin{eqnarray}\bar{\unicode[STIX]{x1D702}}(\boldsymbol{x})=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};\boldsymbol{y})\unicode[STIX]{x1D707}(\,\boldsymbol{y})\,\text{d}\boldsymbol{y},\end{eqnarray}$$

where the wave field generated by infinitely many prior impacts at position $\boldsymbol{x}=\boldsymbol{x}_{0}$ is $\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};\boldsymbol{x}_{0})$, the time-periodic Green’s function for the domain (see figure 8 and appendix C for details). In contrast, the pilot-wave dynamics might instead be periodic, in which case $\boldsymbol{X}_{n+{\mathcal{P}}}=\boldsymbol{X}_{n}$ for all $n$ for some ${\mathcal{P}}\in \mathbb{N}$. As shown in § C.2, equation (3.2) also applies in this setting when we instead define

$$\begin{eqnarray}\unicode[STIX]{x1D707}(\boldsymbol{x})=\frac{1}{{\mathcal{P}}}\mathop{\sum }_{n=1}^{{\mathcal{P}}}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{X}_{n})\quad \text{and}\quad \bar{\unicode[STIX]{x1D702}}(\boldsymbol{x})=\frac{1}{{\mathcal{P}}}\mathop{\sum }_{n=1}^{{\mathcal{P}}}\unicode[STIX]{x1D702}(\boldsymbol{x},t_{n}).\end{eqnarray}$$

While we utilise these two relationships to elucidate the pilot-wave dynamics explored herein, there might also exist regimes in which the dynamics is neither ergodic nor periodic. One may then also consider the mean pilot wave generated over a finite time and its convergence to the asymptotic result given in (3.2). This statistical relaxation time was explored numerically by Tambasco & Bush (Reference Tambasco and Bush2018), where it was observed that the relaxation time scale was comparable to the memory time scale.

Figure 8.

The wave field $\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};\boldsymbol{x}_{0})$ (upper panel) and its $y=0$ cross-section (lower panel). The white circles indicate the corral wall and the black dots denote the droplet position. We set $\boldsymbol{x}_{0}=(r_{0},0)^{\text{T}}$ with (a) $r_{0}=0$, (b) $r_{0}=0.4R_{c}$ and (c) $r_{0}=0.8R_{c}$. The corral radius is $R_{c}=14~\text{mm}$, with frequency $f_{0}=80~\text{Hz}$ and $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.965$, corresponding to memory $Me=19.3$. The droplet radius is $R_{0}=0.35~\text{mm}$ with impact phase $\unicode[STIX]{x1D703}_{I}=1.43\unicode[STIX]{x03C0}$. In the asymmetric cases, we observe three features of $\unicode[STIX]{x1D702}_{G}$: the magnitude of the wave decreases with $r_{0}$, the wave field gradient $\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}_{G}$ is generally non-zero at $\boldsymbol{x}=\boldsymbol{x}_{0}$ (due to the wave reflection off the submerged step) and the wave damping increases in the shallow region, as is most apparent in (c).

To rationalise this observation for dynamics in a circular corral, we consider the pilot-wave dynamics arising on an initially quiescent bath. Moreover, by conditioning over the first $N$ impacts of the droplet’s trajectory $\boldsymbol{X}_{0},\ldots ,\boldsymbol{X}_{N-1}$, we define the statistical distribution, $\unicode[STIX]{x1D707}_{N}(\boldsymbol{x})$, and accompanying mean wave field, $\bar{\unicode[STIX]{x1D702}}_{N}(\boldsymbol{x})$, as

$$\begin{eqnarray}\unicode[STIX]{x1D707}_{N}(\boldsymbol{x})=\frac{1}{N}\mathop{\sum }_{n=0}^{N-1}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{X}_{n})\quad \text{and}\quad \bar{\unicode[STIX]{x1D702}}_{N}(\boldsymbol{x})=\frac{1}{N}\mathop{\sum }_{n=0}^{N-1}\unicode[STIX]{x1D702}(\boldsymbol{x},t_{n}),\end{eqnarray}$$

respectively. As derived in § C.3, the trajectory-dependent mean pilot wave is $\bar{\unicode[STIX]{x1D702}}_{N}(\boldsymbol{x})=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};\boldsymbol{y})\unicode[STIX]{x1D707}_{N}(\,\boldsymbol{y})\,\text{d}\boldsymbol{y}+R_{N}(\boldsymbol{x})$, where the remainder term, $R_{N}(\boldsymbol{x})$, is

$$\begin{eqnarray}R_{N}(\boldsymbol{x})=\frac{1}{N}\unicode[STIX]{x1D702}(\boldsymbol{x},t_{N-1})-\frac{1}{N}\mathop{\sum }_{m=-\infty }^{\infty }\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D731}_{m}(\boldsymbol{x}) & \mathbf{0}\end{array}\right)(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D648}_{m})^{-1}\left(\begin{array}{@{}c@{}}\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{N-1}\\ \unicode[STIX]{x1D753}_{m}^{N-1}\end{array}\right).\end{eqnarray}$$

Strikingly, the remainder term, $R_{N}(\boldsymbol{x})$, depends only on the wave field at the most recent impact. Assuming that the system remains bounded for all time, the first term decays to zero as $N\rightarrow \infty$ with size $O(1/N)$. As the dominant real eigenvalue of $(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D648}_{m})^{-1}$ is ${\mathcal{M}}_{m}(1+O({\mathcal{M}}_{m}^{-1}))$ for ${\mathcal{M}}_{m}\gg 1$, and as $Me=\max _{m}{\mathcal{M}}_{m}$, we deduce that the second term in $R_{N}(\boldsymbol{x})$ has size $O(Me/N)$ as $N\rightarrow \infty$, rationalising the simulation results of Tambasco & Bush (Reference Tambasco and Bush2018). Moreover, the temporal decay of $R_{N}(\boldsymbol{x})$ highlights the robustness of the asymptotic relationship between $\bar{\unicode[STIX]{x1D702}}(\boldsymbol{x})$ and $\unicode[STIX]{x1D707}(\boldsymbol{x})$ given in (3.2).

3.4 The slow-moving limit: an effective potential

Cristea-Platon et al. (Reference Cristea-Platon, Sáenz and Bush2018) reported that when the droplet is confined to a corral, the pilot-wave system may exhibit a similar dynamics to the case where the droplet is instead subjected to a harmonic potential with no submerged topography (Perrard et al. Reference Perrard, Labousse, Fort and Couder2014a,Reference Perrard, Labousse, Miskin, Fort and Couderb). These two systems exhibit qualitatively different trapping mechanisms: in the case of an applied harmonic potential, the central force is applied directly to the droplet, whereas in a corral, the droplet is driven back towards the centre of the corral by the wave reflection off the submerged topography. To rationalise the origins of this similarity, we consider the form of the pilot wave for a droplet walking in the slow-moving limit. In this limit, the pilot wave at any given droplet position $\boldsymbol{x}_{0}$ is well approximated by the wave field $\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};\boldsymbol{x}_{0})$ of a stationary droplet bouncing at that same position (see figure 8). As the stationary droplet is positioned closer towards the submerged wall, the slope underneath the droplet deviates from zero, providing an inward radial force.

For a stationary bouncing at $\boldsymbol{x}_{0}=r\boldsymbol{e}_{r}$, where $\boldsymbol{e}_{r}$ is a radial unit vector, the dimensionless inward force acting on the droplet is $F_{G}(r)=\boldsymbol{e}_{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};r\boldsymbol{e}_{r})|_{\boldsymbol{x}=r\boldsymbol{e}_{r}}$. An effective dimensionless potential, ${\mathcal{V}}(r)$, satisfies ${\mathcal{V}}^{\prime }(r)=F_{G}(r)$. As presented in figure 9, ${\mathcal{V}}(r)$ is approximately quadratic in $r$ away from the origin, strikingly similar to a harmonic potential, suggesting a dynamical correspondence between these two systems in this particular limit. However, unlike an applied potential, the steepness of ${\mathcal{V}}(r)$ increases with the vibrational forcing and the origin is not always a minimum of the potential, as rationalised in § 4.3.

Figure 9.

Dimensionless effective potential, ${\mathcal{V}}(r)$, for a stationary droplet. The vibrational forcing is $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.80$ (black curve), $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.85$ (dash-dotted curve) and $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.90$ (light grey curve). (a) $R_{c}=12.5~\text{mm}$, (b) $R_{c}=14.0~\text{mm}$. Inset (i) ${\mathcal{V}}(r)$ for log–log axes, indicating quadratic-like growth away from the origin. Inset (ii) ${\mathcal{V}}(r)$ in the vicinity of the origin. The system parameters are $f_{0}=80~\text{Hz}$, $R_{0}=0.35~\text{mm}$ and $\unicode[STIX]{x1D703}_{I}=1.43\unicode[STIX]{x03C0}$.

3.5 Summary

In § 3, we explored the fundamental properties of this Faraday problem and pilot-wave system. We demonstrated that the form of the most unstable cavity mode depends on the corral geometry and the vibrational frequency. Furthermore, an effective boundary condition (in the form of a Dirichlet or Neumann condition above the step) is inadequate to describe the interaction with submerged topography, highlighting the limitations of prior works (Blanchette Reference Blanchette2016; Dubertrand et al. Reference Dubertrand, Hubert, Schlagheck, Vandewalle, Bastin and Martin2016; Gilet Reference Gilet2016; Filoux et al. Reference Filoux, Hubert, Schlagheck and Vandewalle2017). We then determined how the droplet statistics relate to the mean pilot wave, which may control the form of the pilot-wave dynamics in certain circumstances (Labousse et al. Reference Labousse, Perrard, Couder and Fort2014; Perrard et al. Reference Perrard, Labousse, Miskin, Fort and Couder2014b). Finally, we provided a rationale for the experimentally observed similitude between this corral system (Cristea-Platon et al. Reference Cristea-Platon, Sáenz and Bush2018) and the dynamics in a harmonic potential (Perrard et al. Reference Perrard, Labousse, Miskin, Fort and Couder2014b), demonstrating that the topography acts as a quadratic-like potential in the slow-moving limit. We proceed to explore the full range of pilot-wave dynamics predicted by this model, with the aim of rationalising further similarities and differences between the dynamics in a corral and a harmonic potential.

4.1 Linear stability analysis

To study the asymptotic linear stability of the orbital and bouncing states, we consider perturbations from the steady state of the form

$$\begin{eqnarray}\boldsymbol{X}_{n}=r_{0}\boldsymbol{e}_{r}+\unicode[STIX]{x1D716}\hat{\boldsymbol{X}}_{n},\quad \boldsymbol{\unicode[STIX]{x1D702}}_{m}^{n}=\boldsymbol{\unicode[STIX]{x1D702}}_{m}^{O}+\unicode[STIX]{x1D716}\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{m}^{n},\quad \unicode[STIX]{x1D753}_{m}^{n}=\unicode[STIX]{x1D753}_{m}^{O}+\unicode[STIX]{x1D716}\hat{\unicode[STIX]{x1D753}}_{m}^{n},\end{eqnarray}$$

where $0<\unicode[STIX]{x1D716}\ll 1$. Consequentially, the perturbed free surface at each impact is $\hat{\unicode[STIX]{x1D702}}(\boldsymbol{x},t_{n})=\sum _{m=-\infty }^{\infty }\unicode[STIX]{x1D731}_{m}(\boldsymbol{x})\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{m}^{n}$. By substituting this perturbation ansatz in the orbital trajectory equation (4.1), the perturbed variables evolve according to the linearised system

(4.5a)$$\begin{eqnarray}\displaystyle & \displaystyle \left(\begin{array}{@{}c@{}}\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{m}^{n}\\ \hat{\unicode[STIX]{x1D753}}_{m}^{n}\end{array}\right)=\text{e}^{\text{i}m\unicode[STIX]{x1D703}_{0}}\unicode[STIX]{x1D648}_{m}\left(\begin{array}{@{}c@{}}\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{m}^{n-1}\\ \hat{\unicode[STIX]{x1D753}}_{m}^{n-1}\end{array}\right)-\left(\begin{array}{@{}c@{}}\mathbf{0}\\ \unicode[STIX]{x1D64B}_{m}\unicode[STIX]{x1D645}_{m}^{\ast }(r_{0}\boldsymbol{ e}_{r})\hat{\boldsymbol{X}}_{n}\end{array}\right), & \displaystyle\end{eqnarray}$$

(4.5b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64D}_{+}\hat{\boldsymbol{X}}_{n+1}-2\hat{\boldsymbol{X}}_{n}+\unicode[STIX]{x1D64D}_{-}\hat{\boldsymbol{X}}_{n-1}=-{\mathcal{D}}(\hat{\boldsymbol{X}}_{n}-\unicode[STIX]{x1D64D}_{-}\hat{\boldsymbol{X}}_{n-1})-{\mathcal{F}}[\unicode[STIX]{x1D643}\hat{\boldsymbol{X}}_{n}+\unicode[STIX]{x1D735}\hat{\unicode[STIX]{x1D702}}(r_{0}\boldsymbol{e}_{r},t_{n})], & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D64D}_{\pm }=\unicode[STIX]{x1D64D}(\pm \unicode[STIX]{x1D703}_{0})$

$\unicode[STIX]{x1D643}$

$\unicode[STIX]{x1D702}_{O}(\boldsymbol{x})$

$\boldsymbol{x}=r_{0}\boldsymbol{e}_{r}$

$\unicode[STIX]{x1D643}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}^{\text{T}}\unicode[STIX]{x1D702}_{O}(r_{0}\boldsymbol{e}_{r})$

We now seek solutions to (4.5) where the perturbed variables evolve according to $\hat{\boldsymbol{X}}_{n}=\unicode[STIX]{x1D706}^{n}\hat{\boldsymbol{X}}_{0}$, $\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{m}^{n}=\unicode[STIX]{x1D706}^{n}\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{m}^{0}$ and $\hat{\unicode[STIX]{x1D753}}_{m}^{n}=\unicode[STIX]{x1D706}^{n}\hat{\unicode[STIX]{x1D753}}_{m}^{0}$ for all $n$, where $\unicode[STIX]{x1D706}\in \mathbb{C}$ is an eigenvalue of the perturbed system (to be determined). By substituting this ansatz into (4.5), we obtain

(4.6a)$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D706}\unicode[STIX]{x1D644}-\text{e}^{\text{i}m\unicode[STIX]{x1D703}_{0}}\unicode[STIX]{x1D648}_{m})\left(\begin{array}{@{}c@{}}\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{m}^{0}\\ \hat{\unicode[STIX]{x1D753}}_{m}^{0}\end{array}\right)=-\unicode[STIX]{x1D706}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D7EC}\\ \unicode[STIX]{x1D64B}_{m}\unicode[STIX]{x1D645}_{m}^{\ast }(r_{0}\boldsymbol{e}_{r})\end{array}\right)\hat{\boldsymbol{X}}_{0}, & \displaystyle\end{eqnarray}$$

(4.6b)$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D706}^{2}\unicode[STIX]{x1D64D}_{+}-2\unicode[STIX]{x1D706}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D64D}_{-})\hat{\boldsymbol{X}}_{0}=-{\mathcal{D}}(\unicode[STIX]{x1D706}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D64D}_{-})\hat{\boldsymbol{X}}_{0}-{\mathcal{F}}\unicode[STIX]{x1D706}[\unicode[STIX]{x1D643}\hat{\boldsymbol{X}}_{0}+\unicode[STIX]{x1D735}\hat{\unicode[STIX]{x1D702}}(r_{0}\boldsymbol{e}_{r},t_{0})]. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D735}\hat{\unicode[STIX]{x1D702}}(r_{0}\boldsymbol{e}_{r},t_{0})=-\unicode[STIX]{x1D706}\unicode[STIX]{x1D642}(\unicode[STIX]{x1D706})\hat{\boldsymbol{X}}_{0}$

$\unicode[STIX]{x1D642}\in \mathbb{C}^{2\times 2}$

$$\begin{eqnarray}\unicode[STIX]{x1D642}(\unicode[STIX]{x1D706})=\mathop{\sum }_{m=-\infty }^{\infty }\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D645}_{m}(r_{0}\boldsymbol{e}_{r}) & \unicode[STIX]{x1D7EC}\end{array}\right)(\unicode[STIX]{x1D706}\unicode[STIX]{x1D644}-\text{e}^{\text{i}m\unicode[STIX]{x1D703}_{0}}\unicode[STIX]{x1D648}_{m})^{-1}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D7EC}\\ \unicode[STIX]{x1D64B}_{m}\unicode[STIX]{x1D645}_{m}^{\ast }(r_{0}\boldsymbol{e}_{r})\end{array}\right).\end{eqnarray}$$

Substituting $\unicode[STIX]{x1D735}\hat{\unicode[STIX]{x1D702}}(r_{0}\boldsymbol{e}_{r},t_{0})$ into (4.6b) yields $\unicode[STIX]{x1D63C}(\unicode[STIX]{x1D706})\unicode[STIX]{x1D64D}_{-}\hat{\boldsymbol{X}}_{0}=\mathbf{0}$, where

(4.7)$$\begin{eqnarray}\unicode[STIX]{x1D63C}(\unicode[STIX]{x1D706})=(\unicode[STIX]{x1D706}\unicode[STIX]{x1D64D}_{+}-\unicode[STIX]{x1D644})^{2}+{\mathcal{D}}(\unicode[STIX]{x1D706}\unicode[STIX]{x1D64D}_{+}-\unicode[STIX]{x1D644})+\unicode[STIX]{x1D706}{\mathcal{F}}(\unicode[STIX]{x1D643}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D642}(\unicode[STIX]{x1D706}))\unicode[STIX]{x1D64D}_{+}.\end{eqnarray}$$

For non-trivial $\hat{\boldsymbol{X}}_{0}$, it remains to find an eigenvalue $\unicode[STIX]{x1D706}\in \mathbb{C}$ such that $\det \unicode[STIX]{x1D63C}(\unicode[STIX]{x1D706})=0$. The system is unstable if there exists $\unicode[STIX]{x1D706}$ such that $|\unicode[STIX]{x1D706}|>1$, and stable otherwise, where we note that rotational invariance (for $r_{0}>0$) supplies that $\unicode[STIX]{x1D706}=1$ is always an eigenvalue.

4.2 The high-memory limit: rationalising the orbital quantisation

The mean wave field, $\bar{\unicode[STIX]{x1D702}}(\boldsymbol{x})$, associated with a circular orbit with radius $r_{0}$ is

$$\begin{eqnarray}\bar{\unicode[STIX]{x1D702}}(\boldsymbol{x})=-\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D731}_{0}(\boldsymbol{x}) & \mathbf{0}\end{array}\right)(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D648}_{0})^{-1}\left(\begin{array}{@{}c@{}}\mathbf{0}\\ \unicode[STIX]{x1D64B}_{0}\unicode[STIX]{x1D731}_{0}^{\ast }(r_{0}\boldsymbol{e}_{r})\end{array}\right),\end{eqnarray}$$

which is exactly equal to the radially symmetric contribution to the orbital wave field, $\unicode[STIX]{x1D702}_{O}(\boldsymbol{x})$, at each impact (see (4.2a)). Hence, the contribution of the mean wave field to the inward radial kick imparted on the droplet at each impact is defined as $F_{O}(r_{0})=\boldsymbol{e}_{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\unicode[STIX]{x1D702}}(r_{0}\boldsymbol{e}_{r})$, where the function $F_{O}(r)$ is given by

(4.8)$$\begin{eqnarray}F_{O}(r)=-\boldsymbol{e}_{r}\boldsymbol{\cdot }\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D645}_{0}(r\boldsymbol{e}_{r}) & \unicode[STIX]{x1D7EC}\end{array}\right)(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D648}_{0})^{-1}\left(\begin{array}{@{}c@{}}\mathbf{0}\\ \unicode[STIX]{x1D64B}_{0}\unicode[STIX]{x1D731}_{0}^{\ast }(r\boldsymbol{e}_{r})\end{array}\right).\end{eqnarray}$$

We now use (4.8) to rationalise the emergence of orbital quantisation in the high-memory limit. We first consider the case where the $m=0$ mode destabilises first (so $\unicode[STIX]{x1D6E4}_{F}=\unicode[STIX]{x1D6E4}_{0}^{c}$ and $Me={\mathcal{M}}_{0}$ close to the Faraday threshold). In the limit $\unicode[STIX]{x1D6E4}\rightarrow \unicode[STIX]{x1D6E4}_{F}$, the dominant real eigenvalue of $(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D648}_{0})^{-1}$ is $Me(1+O(Me^{-1}))$ for $Me\gg 1$. Hence, the magnitude of the radial kick, $|F_{O}(r_{0})|$, grows without bound as $Me\rightarrow \infty$ unless the orbital radius, $r_{0}$, suppresses the singularity. This behaviour may be achieved via two different mechanisms: either the orbital radius, $r_{0}$, minimises the excitation of the radially symmetric mode, or the radially symmetric mode is excited at each impact yet $r_{0}$ lies at a radial extremum. Hence, in the case of an axisymmetric Faraday instability, the origin of orbital quantisation in a corral is analogous to dynamics in an unbounded domain: the corresponding expression for (4.8) in the limit $Me\rightarrow \infty$ may be approximated as $F_{O}(r)/Me\propto \text{J}_{0}(k_{F}r)\text{J}_{0}^{\prime }(k_{F}r)$, where $k_{F}$ is the critical Faraday wavenumber (Oza et al. Reference Oza, Harris, Rosales and Bush2014; Perrard et al. Reference Perrard, Labousse, Miskin, Fort and Couder2014b). Moreover, this quantisation mechanism depends on the droplet properties only through the impact phase, $\unicode[STIX]{x1D703}_{I}$, and is independent of ${\mathcal{D}}$ and ${\mathcal{F}}$.

When the decay time of the axisymmetric mode is instead subdominant to another mode (so $\unicode[STIX]{x1D6E4}_{F}<\unicode[STIX]{x1D6E4}_{0}^{c}$), one may expect a different form of the asymptotic orbital quantisation. Even though ${\mathcal{M}}_{0}$ achieves only a finite value at $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{F}$, the dominant eigenvalue of $(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D648}_{0})^{-1}$ remains sufficiently large to necessitate the suppression of the radial kick, resulting in orbital radii approximately coinciding with the zeros, $r_{\ast }$, of $F_{O}(r)$ as $\unicode[STIX]{x1D6E4}\rightarrow \unicode[STIX]{x1D6E4}_{F}$.

The zeros, $r_{\ast }$, of $F_{O}(r)$ in the limit $\unicode[STIX]{x1D6E4}\rightarrow \unicode[STIX]{x1D6E4}_{F}$ are denoted as black circles in figure 10. When the $m=0$ mode is most unstable ($R_{c}=14~\text{mm}$), the agreement between the asymptotic theory, $r_{\ast }$, and orbital radius, $r_{0}$, is excellent, while a qualitative agreement is obtained for the $m=1$ case ($R_{c}=12.5~\text{mm}$), as anticipated in the above argument. For circular orbits with an oscillatory eigenvalue (green curves), $r_{0}$ minimises the excitation of the mean wave field, $\bar{\unicode[STIX]{x1D702}}$, for which the exponential growth rate from the steady state is typically comparable to the orbital frequency. In contrast, orbits along the extrema of $\bar{\unicode[STIX]{x1D702}}$ have a real unstable eigenvalue (red curves), with a growth rate much larger than the orbital frequency. In the former case, $\bar{\unicode[STIX]{x1D702}}$ appears to control the pilot-wave dynamics through an underlying minimisation principle (Perrard et al. Reference Perrard, Labousse, Miskin, Fort and Couder2014b).

4.3 The bifurcation from bouncing to orbiting

In figure 10, we observe that circular orbits with infinitesimally small radii are permitted for $R_{c}=12.5~\text{mm}$, yet the orbital radius instead approaches a non-zero value for $R_{c}=14~\text{mm}$. We now rationalise these two bifurcations from bouncing to orbiting.

We first consider the limit of vanishing orbital radius, $r_{0}\rightarrow 0$, corresponding to a droplet repeatedly bouncing at the centre of the bath. The accompanying Faraday wave field, $\unicode[STIX]{x1D702}_{B}(\boldsymbol{x})$, for this stationary bouncing state is axisymmetric and related to the time-periodic Green’s function via $\unicode[STIX]{x1D702}_{B}(\boldsymbol{x})=\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};\mathbf{0})$. For weak vibrational forcing, the droplet is trapped by the submerged topography. However, as the vibrational forcing is increased beyond a critical threshold $\unicode[STIX]{x1D6E4}>\unicode[STIX]{x1D6E4}_{0}$ (to be determined), the axisymmetric bouncing state destabilises due to the increased influence of the waves generated by all prior impacts.

The critical threshold, $\unicode[STIX]{x1D6E4}_{0}$, is determined by utilising the stability analysis for circular orbits in the case $\unicode[STIX]{x1D703}_{0}=0$ and $r_{0}=0$. In this case, the rotational symmetry implies that the steady state Hessian matrix, $\unicode[STIX]{x1D643}$, is diagonal, with $\unicode[STIX]{x1D643}={\mathcal{H}}\unicode[STIX]{x1D644}$ and ${\mathcal{H}}=-{\textstyle \frac{1}{2}}\sum _{p=1}^{\infty }k_{0p}^{2}\unicode[STIX]{x1D702}_{0p}^{B}$. Furthermore, by noting that $\unicode[STIX]{x1D645}_{m}(\mathbf{0})=\unicode[STIX]{x1D7EC}$ for all $m$ such that $|m|\neq 1$, and

$$\begin{eqnarray}\unicode[STIX]{x1D645}_{\pm 1}(\mathbf{0})=\frac{1}{2}\left(\begin{array}{@{}c@{}}\pm 1\\ \text{i}\end{array}\right)\boldsymbol{k}_{1}^{\text{T}},\end{eqnarray}$$

where $\boldsymbol{k}_{1}=(k_{11},k_{12},\ldots )^{\text{T}}$ is a column vector of wavenumbers, we obtain $\unicode[STIX]{x1D642}={\mathcal{G}}\unicode[STIX]{x1D644}$, where

$$\begin{eqnarray}{\mathcal{G}}(\unicode[STIX]{x1D706})=\frac{1}{2}\left(\begin{array}{@{}cc@{}}\boldsymbol{k}_{1}^{\text{T}} & \mathbf{0}\end{array}\right)(\unicode[STIX]{x1D706}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D648}_{1})^{-1}\left(\begin{array}{@{}c@{}}\boldsymbol{ 0}\\ \unicode[STIX]{x1D64B}_{1}\boldsymbol{ k}_{1}\end{array}\right).\end{eqnarray}$$

As $\unicode[STIX]{x1D643}$ and $\unicode[STIX]{x1D642}$ both reduce to diagonal matrices, equation (4.7) yields $\unicode[STIX]{x1D63C}={\mathcal{A}}\unicode[STIX]{x1D644}$, where

(4.9)$$\begin{eqnarray}{\mathcal{A}}(\unicode[STIX]{x1D706})=(\unicode[STIX]{x1D706}-1)^{2}+{\mathcal{D}}(\unicode[STIX]{x1D706}-1)+\unicode[STIX]{x1D706}{\mathcal{F}}({\mathcal{H}}-\unicode[STIX]{x1D706}{\mathcal{G}}(\unicode[STIX]{x1D706})).\end{eqnarray}$$

The stability of the axisymmetric bouncing state is determined by the eigenvalues, $\unicode[STIX]{x1D706}$, satisfying ${\mathcal{A}}(\unicode[STIX]{x1D706})=0$. The instability threshold, $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{0}$, corresponds to when the dominant eigenvalue (which is either real or part of a complex-conjugate pair) satisfies $|\unicode[STIX]{x1D706}|=1$. As the topography prohibits translational invariance, $\unicode[STIX]{x1D706}=1$ is not typically an eigenvalue.

When the droplet instead exhibits stationary bouncing at a fixed position $r_{b}\boldsymbol{e}_{r}$, where $r_{b}>0$ and $\boldsymbol{e}_{r}$ is an arbitrary unit vector, the accompanying wave field is $\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};r_{b}\boldsymbol{e}_{r})$. The distance $r_{b}$ satisfies $F_{G}(r_{b})=0$, where $F_{G}(r)=\boldsymbol{e}_{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}_{G}(\boldsymbol{x};r\boldsymbol{e}_{r})|_{\boldsymbol{x}=r\boldsymbol{e}_{r}}$ is the radial slope. We note that $r_{b}$ depends only weakly on the vibrational forcing (see figure 10b). The function $F_{G}(r)$ has at most one non-zero root in the current parameter regime. A necessary condition for this root to exist appears to be $F_{G}^{\prime }(0)<0$ (see figure 12a), where $F_{G}^{\prime }(0)={\mathcal{H}}-{\mathcal{G}}(1)$. Hence, a second bouncing state emerges for geometries where the axisymmetric bouncing state satisfies ${\mathcal{H}}={\mathcal{G}}(1)$, which, from (4.9), corresponds to the instability threshold with $\unicode[STIX]{x1D706}=1$. This second bouncing state is present for the $R_{c}=14~\text{mm}$ corral (figure 10b) but absent for $R_{c}=12.5~\text{mm}$ (figure 10a). The existence of this off-centre bouncing state rationalises why a non-zero minimum orbital radius was reported in experiments (Cristea-Platon et al. Reference Cristea-Platon, Sáenz and Bush2018). Furthermore, from the linear stability analysis of circular orbits (and the simulations in § 5), we see that whenever an off-centre bouncing state exists, the corresponding orbital solution is unstable for small radii (with a complex-conjugate pair of unstable eigenvalues).

Figure 12.

Off-centre bouncing states for $R_{c}=14~\text{mm}$. (a) The function $F_{G}(r)$ for $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.85$ (black) and $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.95$ (grey). Inset: $F_{G}(r)$ has one positive root and $F_{G}^{\prime }(0)<0$. (b) The three unstable eigenvalues, $\unicode[STIX]{x1D706}$, that leave the unit disc (solid curve) nearly simultaneously as $\unicode[STIX]{x1D6E4}$ is increased (denoted by the arrows) from $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.853$ to $\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D6E4}_{F}=0.891$, with increments of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6E4}=0.002\unicode[STIX]{x1D6E4}_{F}$. The black and white circles are to aid visualisation.

The walking threshold, $\unicode[STIX]{x1D6E4}_{W}$, is defined as follows. In the regime where only an axisymmetric bouncing state exists, then $\unicode[STIX]{x1D6E4}_{W}=\unicode[STIX]{x1D6E4}_{0}$. In contrast, when off-centre bouncing states arise, the linear stability may be rationalised using the analysis for circular orbits with $\unicode[STIX]{x1D703}_{0}=0$ and $r_{0}=r_{b}$. Remarkably, the instability of the off-centre bouncing state occurs when three eigenvalues (each satisfying $\det \unicode[STIX]{x1D63C}=0$) almost simultaneously leave the unit disc at $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{W}$ (see figure 12b): two eigenvalues are complex-conjugate pairs that give rise to radial oscillations (in the direction of $\boldsymbol{e}_{r}$), and the remaining real eigenvalue yields transverse growth (in the direction $\boldsymbol{e}_{\unicode[STIX]{x1D703}}$). This behaviour provides strong indication of a codimension-two bifurcation at or near this point in $(\unicode[STIX]{x1D6E4},R_{c})$ parameter space. In fact, when an off-centre bouncing state exists, the corresponding axisymmetric bouncing state is unstable for all vibrational forcing, so $\unicode[STIX]{x1D6E4}_{0}=0$. Hence, at most one stable bouncing state may exist for all of the small ($R_{c}\leqslant 15~\text{mm}$) corral geometries considered.

We present the results of the walking threshold, $\unicode[STIX]{x1D6E4}_{W}/\unicode[STIX]{x1D6E4}_{F}$, in figure 13 (red curve). Off-centre bouncing states typically arise when $\unicode[STIX]{x1D6E4}_{W}/\unicode[STIX]{x1D6E4}_{F}$ is decreasing, although a more precise statement remains the subject of future work. Moreover, the walking threshold approaches a constant as the corral radius increases. To explore the pilot-wave dynamics that arises when the bouncing states and circular orbits are both unstable, we proceed to simulate the iterative map (2.16) for a range of corral radii and vibrational forcings.