1. Introduction
Droplet manipulation on substrates is relevant to technologies such as inkjet printing (Calvert Reference Calvert2001), spray cooling (James, Smith & Glezer Reference James, Smith and Glezer2003; Vukasinovic, Smith & Glezer Reference Vukasinovic, Smith and Glezer2007b; Tsai et al. Reference Tsai, Lin, Mao and Tsai2012) and phase-change heat transfer (Daniel, Chaudhury & Chen Reference Daniel, Chaudhury and Chen2001; Macner, Daniel & Steen Reference Macner, Daniel and Steen2014), and is often facilitated by external driving forces associated with acoustic fields (Marston & Apfel Reference Marston and Apfel1980; Chastrette et al. Reference Chastrette, Baudoin, Brunet, Royon and Wunenburger2022), aerodynamic effects (Bouwhuis et al. Reference Bouwhuis, Winkels, Peters, Brunet, van Der Meer and Snoeijer2013), magnetic fields (Hill & Eaves Reference Hill and Eaves2012), electric fields (Feng & Beard Reference Feng and Beard1991; Oh, Ko & Kang Reference Oh, Ko and Kang2008) and substrate vibrations (Brunet, Eggers & Deegan Reference Brunet, Eggers and Deegan2007; Chang et al. Reference Chang, Bostwick, Daniel and Steen2015). The temporal response of the droplet is often determined by the particular type of applied forcing. In this paper, we study the oscillations of a sessile drop on a mechanically vibrated substrate whose axis of vibration is oriented at an angle 
$\chi$ from the normal. This gives rise to parametric oscillations typically associated with Faraday waves, which contribute to the unidirectional drop motion (Costalonga & Brunet Reference Costalonga and Brunet2020; Guo, Shui & Chen Reference Guo, Shui and Chen2020) and climbing drops (Brunet et al. Reference Brunet, Eggers and Deegan2007; Sartori et al. Reference Sartori, Quagliati, Varagnolo, Pierno, Mistura, Magaletti and Casciola2015).
Faraday (Reference Faraday1831) observed the formation of surface waves on a vertically vibrated liquid bath that oscillated with frequency 
$f$ one-half that of the driving frequency 
$f_{d}$, i.e. 
$f=f_d/2$, above a critical threshold driving acceleration 
$a$. This is termed a subharmonic response. Benjamin & Ursell (Reference Benjamin and Ursell1954) showed that Faraday waves were solutions of a generalized Mathieu equation,
\begin{equation} \ddot{\boldsymbol{y}}+\boldsymbol{c} \dot{\boldsymbol{y}}+ (\boldsymbol{P}-2A\boldsymbol{Q}\cos f_{d}t)\boldsymbol{y} =0, \end{equation}
where 
$\boldsymbol {y}$ is the surface disturbance. Floquet theory is often invoked to solve (1.1) and gives rise to a set of instability tongues for a subharmonic, harmonic and superharmonic wave response in which the wave is observed. There is a large volume of literature on Faraday waves on planar interfaces (Miles & Henderson Reference Miles and Henderson1990). However, there is limited research on Faraday waves on drops, particularly theoretical modelling. Ebo-Adou & Tuckerman (Reference Ebo-Adou and Tuckerman2016) use Floquet theory to predict the instability tongues for a free drop subject to radial forcing. Maksymov & Pototsky (Reference Maksymov and Pototsky2019) develop a simple hydrodynamic model of a one-dimensional pancake-shaped drop, and analytically determine the response amplitudes for the first harmonic and first superharmonic. To our knowledge, there has been no theoretical investigation of Faraday waves on a sessile drop with arbitrary contact angle 
$\alpha$ and slanting angle 
$\chi$. The aim of this study is to fill this gap in the literature. Unlike the case of radial forcing, slanted vibrations lead to a set of coupled Mathieu equations, which we solve using the Floquet theory method proposed by Kidambi (Reference Kidambi2013) for a brimful circular cylinder.
The frequency spectrum for the sessile drop exhibits a rich structure that has been predicted theoretically by Bostwick & Steen (Reference Bostwick and Steen2014) and verified experimentally under both terrestrial (Chang et al. Reference Chang, Bostwick, Daniel and Steen2015) and microgravity (McCraney et al. Reference McCraney, Kern, Bostwick, Daniel and Steen2022) conditions, and numerically by Sakakeeny & Ling (Reference Sakakeeny and Ling2020). Here, the drop will oscillate with frequency 
$\lambda _{\ell,m}$ and mode shape defined by the mode number pair 
$[\ell,m]$, where 
$\ell$ is the polar mode number, and 
$m$ is the azimuthal mode number, restricted such that 
$\ell +m=\textrm {even}$ due to symmetry considerations associated with enforcement of the no-penetration condition at the solid support. For reference, the 
$\ell +m=\textrm {odd}$ modes have been predicted recently by Ding & Bostwick (Reference Ding and Bostwick2022), and are associated with a flux condition on the solid support. Figure 1 shows typical mode shapes, which can be classified by their mode number pair 
$[\ell, m]$ into zonal 
$[\ell,0]$, lateral 
$[\ell,1]$, sectoral 
$[\ell,\ell ]$ and tesseral 
$[\ell,m\neq \ell ]$ shapes. The spectrum of an inviscid sessile drop depends upon the static contact angle 
$\alpha$ and contact-line boundary condition. The special case of a hemispherical 
$\alpha =90^\circ$ base state with fully mobile contact line exhibits a high degree of symmetry and inherits the Rayleigh–Lamb spectrum for the free drop, including the well-known spectral degeneracy with respect to azimuthal mode number 
$m$ (Rayleigh Reference Rayleigh1879; Lamb Reference Lamb1924). Deviations from this special base state lead to spectral splitting and spectral reordering, which have been organized using a symmetry-breaking perspective into a periodic table of droplet motions by Steen, Chang & Bostwick (Reference Steen, Chang and Bostwick2019). Recent work by Zhang, Zhou & Ding (Reference Zhang, Zhou and Ding2023) has revealed symmetry-breaking due to gravitational effects, demonstrating that the lowest-frequency mode is transformed from a zonal mode to a sectoral mode.
Sessile drop mode shape 
$[\ell,m]$ classification into (a) zonal 
$[\ell,0]$, (b) lateral 
$[\ell,1]$, and (c) sectoral 
$[\ell,\ell ]$ modes.
Sessile drop mobility is achieved once the pinning forces at the three-phase contact line are overcome. Daniel & Chaudhury (Reference Daniel and Chaudhury2002) and Daniel et al. (Reference Daniel, Sircar, Gliem and Chaudhury2004) showed that a drop placed on a chemically anisotropic substrate will undergo translational motion once the contact-angle hysteresis has been overcome so as to mobilize the contact line. Here, substrate vibration is shown to help mobilize the contact line and induce such motions through resonant shape oscillations. This effect can become so strong that a drop may be driven up an inclined substrate against the weight of gravity (Brunet et al. Reference Brunet, Eggers and Deegan2007; Sartori et al. Reference Sartori, Quagliati, Varagnolo, Pierno, Mistura, Magaletti and Casciola2015). Noblin, Kofman & Celestini (Reference Noblin, Kofman and Celestini2009) have shown that such motion can be decomposed into components associated with normal and tangential substrate vibrations. To explain and quantitatively predict directional drop motion, direct numerical simulations (Ding et al. Reference Ding, Zhu, Gao and Lu2018) and simplified flow models, e.g. lubrication theory (Savva & Kalliadasis Reference Savva and Kalliadasis2014) and shallow-water models (Benilov & Billingham Reference Benilov and Billingham2011; Bradshaw & Billingham Reference Bradshaw and Billingham2016), have been put forth to simulate a drop climbing up a vibrating inclined plate. Costalonga & Brunet (Reference Costalonga and Brunet2020) conducted quantitative experiments of directional motion of sessile drops induced by slanted vibrations to compare with existing models. It is generally accepted that to induce directional drop motion requires the simultaneous excitation of the asymmetric 
$[1,1]$ rocking mode and the symmetric 
$[2,0]$ pumping mode as the minimal ingredients (Noblin et al. Reference Noblin, Kofman and Celestini2009; John & Thiele Reference John and Thiele2010). It is possible to do so in an experiment with a single driving frequency. This is illustrated in figure 2, which plots the subharmonic and harmonic tongues for the rocking and pumping modes in the acceleration–frequency (
$\varLambda \unicode{x2013}\omega$) parameter space for a typical set of experimental conditions. Note the overlap of the instability tongues for the subharmonic rocking mode and harmonic pumping mode, which gives rise to drop motion that we will refer to hereafter as a mixed mode. Thus the mixed mode is defined as the superposition of a subharmonic rocking mode and harmonic pumping mode. For reference, other mixed mode pairs have been observed experimentally that consist of a harmonic zonal mode and subharmonic tesseral or sectoral mode (Chang et al. Reference Chang, Bostwick, Daniel and Steen2015, figure 12). What is also interesting from figure 2 is that the harmonic 
$[1,1]$ rocking mode is isolated, which makes it easy to observe in experiment. Given the interest in such motions, we focus on presenting theoretical predictions for the 
$[1,1]$ and 
$[2,0]$ modes, and conduct preliminary experiments to verify these theoretical predictions.
Typical instability tongues for the 
$[1,1]$ and 
$[2,0]$ modes plotted in the acceleration–frequency (
$\varLambda \unicode{x2013}\omega$) space. Here, the contact line is pinned, the slanting angle is 
$\chi =30^\circ$, and the contact angle is 
$\alpha =90^\circ$.
We begin this paper by describing our mathematical model of surface waves on a sessile drop subjected to slanted substrate vibrations. This gives rise to a set of coupled Mathieu equations, which are solved using Floquet theory to give the shape of the instability tongues, as they depend upon static contact angle 
$\alpha$, slanting angle 
$\chi$, and contact-line conditions. Experiments are performed to test theoretical predictions and are focused on the 
$[1,1]$ rocking and 
$[2,0]$ pumping modes, as relevant to drop transport on substrates. Theoretical predictions show good agreement with experimental observations for most cases, but there are exceptions, and we provide a rationalization for these differences. We finish by offering some concluding remarks.
2. Mathematical formulation
Consider a sessile drop defined by a base radius 
$r$ and contact angle 
$\alpha$ on a horizontal substrate, as shown in figure 3, and defined as
\begin{align} X(s,\varphi;\alpha)=\frac{\sin (s)}{\sin (\alpha)} \cos (\varphi),\quad Y(s,\varphi;\alpha)=\frac{\sin (s)}{\sin (\alpha)} \sin (\varphi),\quad Z(s;\alpha)= \frac{\cos (s) - \cos (\alpha)}{\sin (\alpha)},\end{align}
using arc-length-like 
$s\in [0,\alpha ]$ and azimuthal angle 
$\varphi \in [0,2{\rm \pi} ]$ as generalized surface coordinates (Bostwick & Steen Reference Bostwick and Steen2014). The liquid is immersed in a passive gas, has density 
$\varrho$, and is assumed to be incompressible and inviscid. Gravitational effects are neglected. The plate is driven by an applied acceleration 
$a_{tt}$ that is oriented at an angle 
$\chi$ away from the vertical. This induces an interface disturbance 
$\eta$ and corresponding free surface flow that is described by a velocity field 
$\boldsymbol {U}$ and pressure field 
$P$ that satisfy the continuity and Euler equations, i.e.
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{U}=0, \end{gather}
\begin{gather}\varrho\left(\frac{\partial}{\partial t}+\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{U}={-}\boldsymbol{\nabla}P+\varrho a_{tt}\boldsymbol{k}, \end{gather}respectively.
Definition sketch of a sessile drop on a substrate that is mechanically vibrated at angle 
$\chi$ defined from the normal in (a) two-dimensional planar and (b) three-dimensional perspective views.
At the free interface, the interface shape 
$\eta$ is coupled to the flow field through the kinematic condition
\begin{equation} \frac{\partial \eta}{\partial t}+\boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla}\eta= \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{n},\end{equation}and the Young–Laplace equation
\begin{equation} P/\sigma=2 H,\end{equation}
which relates the pressure jump across the interface to the mean curvature 
$H$ there, with 
$\sigma$ the surface tension. Finally, the no-penetration condition is applied on the solid boundary:
\begin{equation} \boldsymbol{U}\boldsymbol{\cdot}\hat{z}=0\lvert_{z=0}.\end{equation}
The orientation of the applied acceleration gives rise to both horizontal and vertical forcing, and can be expressed in component form as 
$a_{tt}\boldsymbol {k}=\cos \chi \,a_{vtt}\boldsymbol {z}+\sin \chi \, a_{htt}\boldsymbol {x}$. Here, we assume 
$a_{vtt}=a_{htt}=({\sigma }/{\varrho r^2})\varLambda \cos \omega t$, with 
$\varLambda \equiv A/r$, where 
$A$ is the driving amplitude, and 
$\omega$ is the driving frequency.
2.1. Linearized equations
Equation (2.2) admits the base state
\begin{equation} \hat{\boldsymbol{U}}=0,\quad \boldsymbol{\nabla}\hat{P}= \frac{\sigma}{r^2}\,\varLambda(\cos\chi\,\boldsymbol{z}+\sin \chi\,\boldsymbol{x})\cos \omega t . \end{equation}We perform a linear stability analysis about this base state by decomposing the velocity and pressure fields as
\begin{gather} \boldsymbol{U}=\hat{\boldsymbol{U}}+ \epsilon\,\boldsymbol{u}= \epsilon \boldsymbol{\nabla}\phi, \end{gather}
\begin{gather}P=\hat{P}+ \epsilon p, \end{gather}
where 
$\epsilon$ is a small parameter. Here, we have assumed an irrotational velocity field that is described by the velocity potential 
$\phi$, such that 
$\boldsymbol {u}=\boldsymbol {\nabla } \phi$. Substituting (2.7) into (2.2) leads to the linearized equations
\begin{gather} \nabla^2\phi=0, \end{gather}
\begin{gather}p={-}\varrho\,\frac{\partial \phi}{\partial t}, \end{gather}on the drop domain. Similarly, the kinematic condition (2.3) reduces to
\begin{equation} \frac{\partial \phi}{\partial n}=\frac{\partial \eta}{\partial t},\end{equation}and the no-penetration condition (2.5) reduces to
\begin{equation} \frac{\partial \phi}{\partial z}=0.\end{equation}The Young–Laplace equation (2.4) can be written as
\begin{equation} P|_{\rho=R+ \eta}\approx\hat{P}|_{\rho=R}+ p|_{\rho=R}+\eta \left.\left(\frac{\partial \hat{P}}{\partial \rho}\right)\right|_{\rho=R}= \sigma(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n})_{\rho=R+ \eta},\end{equation}
using a spherical coordinate system 
$(\rho,\theta, \varphi )$ with the static drop surface defined at 
$\rho =R$. These linearized equations can be simplified by defining 
$\hat {P}|_{\rho =R}=\sigma (\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {n})_{\rho =R}$ and 
$p|_{\rho =R}=-\varrho ({\partial \phi }/{\partial t})$, and performing some algebraic manipulations to get
\begin{align} & \varrho\,\frac{\partial \phi}{\partial t}+\sigma(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{n})_{\rho=R+ \eta}-\sigma(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n})_{\rho=R} \nonumber\\ &\quad - \frac{\sigma}{ r^2}\,\varLambda(\cos\chi \cos\theta+\sin\chi \sin\theta \cos\varphi)\cos\omega t\eta=0. \end{align}2.2. Dimensionless equations
Dimensionless variables are introduced,
\begin{equation} \rho^{*}=\rho/r,\quad \eta^{*}=\eta/r,\quad t^{*}=t\,\sqrt{\frac{\sigma}{\varrho r^3}},\quad \phi^{*}=\phi\,\sqrt{\frac{\varrho}{\sigma r}},\quad p^{*}=p\,\frac{r}{\sigma}, \end{equation}
and applied to the governing equations. Herein, we drop the 
$*$ and refer to dimensionless quantities. The resulting equations to be solved are Laplace's equation on the domain
\begin{equation} \nabla^2\phi=0, \end{equation}
with no-penetration condition on the solid support 
$z=0$,
\begin{equation} \frac{\partial \phi}{\partial z}=0, \end{equation}and kinematic condition
\begin{equation} \frac{\partial \phi}{\partial n}=\frac{\partial \eta}{\partial t} \end{equation}and Young–Laplace equation
\begin{align} & \frac{\partial \phi}{\partial t}-\sin^2 (\alpha)\left(\frac{\partial^2\eta}{\partial s^2}+\cot (s)\, \frac{\partial\eta}{\partial s}+2\eta+\frac{1}{\sin^2(s)}\,\frac{\partial^2\eta}{\partial \psi^2}\right) \nonumber\\ &\quad -\varLambda \eta (\cos\chi \cos\theta+\sin\chi \sin\theta \cos\varphi) \cos \omega t=0 \end{align}
on the free surface 
$\rho =R$. The governing equations are augmented with a contact-line boundary condition applied at 
$s=\alpha$, which we assume to be either pinned or free to move, i.e.
\begin{gather} \eta=0,\quad \text{pinned condition}, \end{gather}
\begin{gather}\frac{\partial \eta}{\partial s}-(\cos \alpha) \eta=0,\quad \text{free condition}, \end{gather}and an integral condition
\begin{equation} \int_{\varGamma}\frac{\partial \phi}{\partial n}\,{\rm d}\varGamma=0,\end{equation}necessary to ensure volume conservation.
2.3. Derivation of the coupled Mathieu equations
We seek a solution 
$(\eta, \phi )$ to the governing equations (2.14), (2.15), (2.16), (2.17), (2.18), (2.19) in the form
\begin{align} \eta(x,t)=\cos m\varphi\sum_{\ell =1}^{\infty}a_{\ell}(t)\,V^{m}_{\ell}(x),\quad \phi(\rho,\theta,t)=\cos m\varphi\sum_{\ell =1}^{\infty}b_{\ell}(t)\, \rho^{\ell}\,P^{m}_{\ell}(\cos\theta),\end{align}
where 
$x=\cos s$. Here, 
$V^{m}_{\ell }$ are a set of orthonormal basis functions chosen to satisfy the contact line (2.18) and integral conditions (2.19). The idea is to choose basis functions from the natural oscillations problem, as described by Bostwick & Steen (Reference Bostwick and Steen2014, (4.10)), as
\begin{equation} V^{m}_{\ell}(x)=\int_{b}^{1}G(x,y;m)\,(\rho^{\ell}\,P^{m}_{\ell}(y)+ \delta_{m,0}C)\,{{\rm d}y},\quad \ell=1,2,\ldots, N, \end{equation}where the Green's function is defined in
\begin{equation} G(x,y;m) =
\left\{\begin{array}{@{}ll}
\displaystyle\dfrac{1}{1-y^2}\,\dfrac{U(x;m)\,X(y;m)}{W(y;m)},
& b< y< x<1, \\
\displaystyle\dfrac{1}{1-y^2}\,\dfrac{U(y;m)\,X(x;m)}{W(y;m)},
& b< x< y<1 . \end{array}\right.
\end{equation}
Here, 
$b \equiv \cos \alpha$, and 
$U$ and 
$X$ are the homogeneous solutions of
\begin{equation} \left(2-\frac{m^2}{\sin^2(s)}\right)\frac{\partial \phi}{\partial n}+ \cot(s) \left(\frac{\partial \phi}{\partial n}\right)' + \left(\frac{\partial \phi}{\partial n} \right)''=0 \end{equation}that satisfy the associated boundary conditions
\begin{equation} U=y_{1}(x;m),\quad X=y_{2}(x;m)-\frac{\tau_{2}(m)}{\tau_{1}(m)}\,y_{1}(x;m), \end{equation}
and 
$W$ is the Wronskian of the solutions 
$U$ and 
$X$.
The Green's function is parametrized by the azimuthal mode number 
$m$ and the transformed contact angle 
$b$. The functions 
$y_{1}$ and 
$y_{2}$ in the solutions 
$U$ and 
$X$ are given by
\begin{align} \left.\begin{gathered}
y_{1}(x;0)=P_{1}(x),\quad y_{2}(x;0)=Q_{1}(x),\quad
y_{1}(x;1)=P_{1}^1(x),\quad y_{2}(x;1)=Q_{1}^1(x),\\
y_{1}(x;m\geq2)=(x+m)\left(\frac{1-x}{1+x}\right)^{m/2},\quad
y_{2}(x;m\geq2)=\frac{x-m}{2m(m^2-1)}\left(\frac{1+x}{1-x}\right)^{m/2},
\end{gathered}\right\}
\end{align}
where 
$P_{1}$, 
$Q_{1}$ and 
$P_{1}^1$, 
$Q_{1}^1$ are the Legendre functions of orders 0 and 1 and index 1, respectively. The parameters 
$\tau _{1}$ and 
$\tau _{2}$ are related to the contact-line boundary conditions, with
\begin{equation} \tau_{1}^{p}=y_{1}(b;m),\quad \tau_{2}^{p}=y_{2}(b;m),\end{equation}
for the pinned contact-line disturbance (superscript 
$p$), and
\begin{equation} \tau_{1}^{n}=y_{1}'(b;m)+\frac{b}{\sqrt{1-b^2}}\,y_{1}(b;m),\quad \tau_{2}^{n}=y_{2}'(b;m)+\frac{b}{\sqrt{1-b^2}}\,y_{2}(b;m),\end{equation}
for the natural contact-line disturbance (superscript 
$n$). The constant 
$C$ is defined as
\begin{equation} C={-}\frac{\displaystyle\int_{b}^{1}\displaystyle\int_{b}^{1}G(x,y;0)\, \phi(y)\,{\rm d} y\,{\rm d}{\kern0.8pt}x}{\displaystyle\int_{b}^{1}\int_{b}^{1}G(x,y;0)\,{\rm d} y\,{\rm d}{\kern0.8pt}x}. \end{equation}Substituting (2.20a,b) into (2.16) and (2.17) gives
\begin{equation} \sum_{\ell =1}^{\infty}b_{\ell}(t)\,\frac{\partial (\rho^{\ell}\,P^{m}_{\ell}(\cos\theta))}{\partial n}= \frac{{\rm d} a_{\ell}(t)}{{\rm d} t}\,V^{m}_{\ell}(x) \end{equation}and
\begin{align}
&\cos m \varphi \sum_{\ell =1}^{\infty}\frac{{\rm d}
b_{\ell}(t)}{{\rm d} t}\,
(\rho^{\ell}\,P^{m}_{\ell}(\cos\theta)) \nonumber\\
&\quad {}\times\left(-\sin^2\alpha\left((1-x^2)\,\frac{\partial^2
}{\partial x^2}-2x\, \frac{\partial}{\partial
x}+2-\frac{m^2}{1-x^2}\right)V^{m}_{\ell}(x)\right)a_{\ell}(t) \nonumber\\
&\quad {}- V^{m}_{\ell}(x)\,a_{\ell}(t)(\varLambda\cos\chi
\cos\theta+\varLambda\sin\chi \sin\theta
\cos\varphi)\cos\omega t=0,
\end{align}
respectively. We can project these equations onto the function space 
$\cos m\varphi \,V^{m}_{k}(x), k=1,$ 
$2,\ldots, N,$ via the inner product 
$\langle{\kern1.8pt}f,g\rangle =\int _{b}^{1}f(x)\,g(x)\,\textrm {d}{\kern0.8pt}x$, noting the identity
\begin{equation} \int_{0}^{2{\rm \pi}}\int_{b}^{1}\cos^2 m \varphi\,\varLambda_{0}\,V^{m}_{\ell}(x)\, V^{m}_{k}(x)\sin\chi \sin\theta \cos\varphi \cos\omega t=0. \end{equation}We project (2.28a) onto the function space to give
\begin{equation} b_{\ell}(t)=\sum^{\infty}_{k=1}\beta_{\ell k}\,\frac{{\rm d} a_{k}}{{\rm d} t},\end{equation}with
\begin{equation} \beta_{\ell k}=\left(\int_{b}^{1}\frac{\partial (\rho^{\ell}\,P^{m}_{\ell} (\cos\theta))}{\partial n}\,V^{m}_{k}(x)\,{{\rm d}{\kern0.8pt}x}\right)^{{-}1} \int_{b}^{1} V_{\ell}^{m}(x)\,V_{k}^{m}(x)\,{{\rm d}{\kern0.8pt}x}. \end{equation}
Similarly, projecting (2.28b) onto the function space 
$V^{m}_{\ell }$ yields
\begin{align} & \sum_{\ell =1}^{\infty}\frac{{\rm d} b_{\ell}(t)}{{\rm d} t}\,(\rho^{\ell}\,P^{m}_{\ell}(\cos\theta)) \nonumber\\ &\quad \times\left(-\sin^2\alpha\left((1-x^2)\,\frac{\partial^2 }{\partial x^2}-2x\, \frac{\partial}{\partial x}+2-\frac{m^2}{1-x^2}\right)V^{m}_{\ell}(x)\right)a_{\ell}(t) \nonumber\\ &\quad -\varLambda\cos\chi \cos\theta \cos\omega t\,V^{m}_{\ell}(x)\,a_{\ell}(t)=0. \end{align}
Finally, we can project (2.32) onto the function space 
$V^{m}_{p}$ and use (2.30) to obtain a system of 
$N$ coupled second-order ordinary differential equations:
\begin{equation} \sum_{\ell =1}^{N}\gamma_{p\ell}\sum^{N}_{k=1}\beta_{\ell k}\, \frac{{\rm d}^2 a_{k}}{{\rm d} t^2}+(\tau_{p\ell}-\varLambda\cos \chi\,H_{p\ell}\cos\omega t)a_{\ell}=0,\quad p=1,2,\ldots, N, \end{equation}with
\begin{gather} \gamma_{p\ell}=\int_{b}^{1}\rho^{\ell}\,P^{m}_{\ell}(\cos\theta)\,V^{m}_{p}(x)\,{{\rm d}{\kern0.8pt}x}, \end{gather}
\begin{gather}\tau_{p\ell}=\sin^2\alpha\int_{b}^{1}-\left((1-x^2)\,\frac{\partial^2 }{\partial x^2}-2x\, \frac{\partial}{\partial x}+2-\frac{m^2}{1-x^2}\right)V^{m}_{\ell}(x)\,V^{m}_{p}(x)\,{{\rm d}{\kern0.8pt}x}, \end{gather}
\begin{gather}H_{p\ell}=\int_{b}^{1} \cos\theta\,V^{m}_{\ell}(x)\,V^{m}_{p}(x)\,{{\rm d}{\kern0.8pt}x}. \end{gather} As is customary, we have kept only a finite number of modes for 
$N$ for numerical computation. Equation (2.33) can be written in matrix form as
\begin{equation}
\boldsymbol{\mathsf{A}}\,\frac{{\rm d}^2\boldsymbol{y}}{{\rm d}
t^2}+ (\hat{\boldsymbol{\mathsf{P}}}-2\varLambda\cos\chi\,\hat{\boldsymbol{\mathsf{Q}}}\cos\omega
t)\boldsymbol{y}=0, \end{equation}
where 
$\boldsymbol {y}=(a_{1},a_{2},\ldots,a_{N})$ is a coefficient vector, and the components of the matrices 
$\boldsymbol{\mathsf{A}}$ and 
$\hat {\boldsymbol{\mathsf{P}}}$ are given by
\begin{equation}
{\mathsf{A}}_{ij}=\sum_{\ell=1}^{N}\gamma_{i\ell}\beta_{\ell j},\quad
\hat{{\mathsf{P}}}_{ij}=\tau_{ij},\quad \hat{{\mathsf{Q}}}_{ij}={\mathsf{H}}_{ij}/2.
\end{equation}
To facilitate a solution, we rescale (2.35) by defining 
$T=\omega t/2$, 
$\boldsymbol{\mathsf{P}}=4\boldsymbol{\mathsf{A}}^{-1}\hat {\boldsymbol {\boldsymbol{\mathsf{P}}}}$, 
$\boldsymbol{\mathsf{Q}}=4\boldsymbol{\mathsf{A}}^{-1}\hat {\boldsymbol{\mathsf{Q}}}$, which gives
\begin{equation}
\omega^2\,\frac{{\rm d}^2\boldsymbol{y}}{{\rm d} T^2}+
(\boldsymbol{\mathsf{P}}-2\varLambda\cos\chi\,\boldsymbol{\mathsf{Q}}\cos
2T)\boldsymbol{y}=0.
\end{equation}
2.4. Solution method
Floquet theory is used to determine the stability of (2.37). We follow the method suggested by Kumar & Tuckerman (Reference Kumar and Tuckerman1994) and Kidambi (Reference Kidambi2013) by fixing the values of the Floquet exponent 
$\mu$ for both harmonic and subharmonic tongues to trace out the shape of the instability tongues. In general, this is a computationally efficient method.
We begin by seeking a solution of the form
\begin{equation} \boldsymbol{y}(T)={\rm e}^{\mu T}\left(\sum_{l={-}L}^{L}\boldsymbol{\xi}_{l}\exp({{\rm i}2lT})\right),\end{equation}
with 
$\boldsymbol {\xi }_{l}$ a vector with components 
$\xi _{l}^{n}$, 
$n=1,\ldots, N$. Substituting (2.38) into (2.37) gives
\begin{equation}
\sum_{l={-}L}^{L}\left(\omega^2(\mu+i2l)^2
\xi_{l}^{n}+ \sum_{j=1}^{N}({\mathsf{P}}_{nj}-2\varLambda
\cos\chi\,{\mathsf{Q}}_{nj}\cos2T)\xi_{l}^{j}\right)=0. \end{equation}
For a given 
$\varLambda$, we can compute the stability boundary by setting 
$\mu =0$ for the harmonic case and 
$\mu =i$ for the subharmonic case, and solving for the unknown 
$\omega$ from the generalized eigenvalue problem 
$(\boldsymbol {M}\omega ^2+\boldsymbol {W})\boldsymbol {v}=0$, where 
$\boldsymbol {v}=(\xi _{-L}^{1},\xi _{-L+1}^{1},\ldots,\xi _{L-1}^{1},\xi _{L}^{1}, \xi _{-L}^{2},\ldots,\xi _{L-1}^{N},\xi _{L}^{N})^{T}$.
3. Theoretical predictions
For the results presented here, we use a truncation 
$N=10$ and 
$L=20$ to produce sufficient convergence for the Floquet theory calculations. We begin by verifying our model against prior results in the literature. We then examine the limiting cases of horizontal and vertical vibrations, before discussing the more general case of slanted vibrations. Our focus is on the low-frequency 
$[1,1]$ and 
$[2,0]$ modes as the main contributors to directional drop transport, as discussed in the Introduction.
3.1. Radial vibration verification of solution method
To our knowledge, only Ebo-Adou & Tuckerman (Reference Ebo-Adou and Tuckerman2016) have considered the parametric oscillations of a free spherical drop, and they focus on the case where motion is driven by a time-harmonic radial acceleration. A sessile drop with hemispherical base state 
$\alpha =90^\circ$ and free contact line has a subset of motions that may be extended smoothly by symmetry operation to the full drop (Steen et al. Reference Steen, Chang and Bostwick2019). This allows us to verify our Floquet theory. It is straightforward to adapt our analysis to this limiting case, where (2.17) becomes
\begin{equation} \frac{\partial \phi}{\partial t}-\sin^2 (\alpha)\left(\frac{\partial^2\eta}{\partial s^2}+\cot (s)\, \frac{\partial\eta}{\partial s}+2\eta+\frac{1}{\sin^2(s)}\, \frac{\partial^2\eta}{\partial \psi^2}\right) -\varLambda \eta \cos \omega t=0. \end{equation}
The steps from § 2.3 may then be applied to generate a set of uncoupled Mathieu equations from which the subharmonic and harmonic tongues are computed by Floquet theory. Figure 4 compares the shape of the subharmonic instability tongues for the 
$[2,0]$, 
$[4,0]$, 
$[6,0]$ hemispherical sessile drop modes to those for the free drop, which reveals excellent agreement, verifying our Floquet theory for sessile drops.
Verification of Floquet theory for sessile drops with radial forcing. Subharmonic instability tongues for the 
$[2,0]$, 
$[4,0]$, 
$[6,0]$ hemispherical sessile drop modes with free contact line compared with those predicted by Ebo-Adou & Tuckerman (Reference Ebo-Adou and Tuckerman2016) for the full spherical drop.
3.2. Horizontal vibration
Horizontal vibration 
$\chi =90^\circ$ is a limiting case where (2.37) becomes
\begin{equation} \omega^2\,\frac{{\rm d}^2\boldsymbol{y}}{{\rm d} T^2}+\boldsymbol{\mathsf{P}}\boldsymbol{y}=0. \end{equation}
Here, there is no parametric term implying that the temporal response is harmonic, which has been observed experimentally by Khan & Eslamian (Reference Khan and Eslamian2019) for surface waves on a water layer with pinned contact line, and argued mathematically by Or (Reference Or1997) for a liquid film with infinite lateral extent. We similarly report a harmonic drop response in experiments for slanting angle close to 
$\chi =90^\circ$, which we will discuss in the next subsection. For reference, the frequency spectrum for the hemispherical sessile drop subject to horizontal forcing has been reported by Lyubimov, Lyubimova & Shklyaev (Reference Lyubimov, Lyubimova and Shklyaev2004).
3.3. Vertical vibration
Vertical vibration 
$\chi =0^\circ$ is another limiting case that is relevant to most experimental studies. Here, (2.37) becomes
\begin{equation} \omega^2\,\frac{d^2\boldsymbol{y}}{dT^2}+(\boldsymbol{\mathsf{P}}-2\varLambda\boldsymbol{\mathsf{Q}}\cos2 T)\boldsymbol{y}=0. \end{equation}
Figure 5 plots the instability tongues for the three modes 
$[2,0]$, 
$[1,1]$, 
$[2,2]$ with smallest natural frequencies for the case of vertical vibration. For both harmonic and subharmonic tongues, note that the pinned contact-line modes are shifted to higher frequency relative to the free contact-line modes. This is due to the fact that the pinned contact line is more constrained than the free contact line, which is known to lead to higher frequencies according to variational principles (Bostwick & Steen Reference Bostwick and Steen2015). Also, the instability tongues for the free contact line are notably wider than those for the pinned contact line, a feature that is more prominent for the harmonic tongues (cf. figure 5b). A comparison between modes shows that the 
$[2,0]$ mode has the largest instability tongue, which we attribute to the fact that zonal modes are dominated by vertical motions, making them easier to be excited by vertical vibration. Figure 6 shows how the shape of the instability tongues depends upon the contact angle 
$\alpha$. Here, the bandwidth increases with 
$\alpha$, suggesting that larger drop inertia makes it easier to excite waves.
(a) Subharmonic and (b) harmonic instability tongues for the 
$[2,0]$, 
$[1,1]$ and 
$[2,2]$ modes plotted in the acceleration–frequency (
$\varLambda - \omega$) space, contrasting pinned and free contact-line conditions for 
$\alpha =60^\circ$.
(a,c) Subharmonic and (b,d) harmonic instability tongues for the (a,b) 
$[2,0]$ and (c,d) 
$[1,1]$ modes plotted in the acceleration–frequency (
$\varLambda - \omega$) space, depending upon contact angle (
$\alpha =60^\circ, 90^\circ, 120^\circ$).
3.4. Slanted vibration
Equation (2.37) describes the most general case for slanted vibration, as defined by the slanting angle 
$\chi$. Here, the 
$\cos \chi$ coefficient on the parametric forcing term strongly affects the shape of the instability tongues. Our focus is on the 
$[1,1]$ harmonic mode, 
$[2,0]$ harmonic mode and 
$[1,1]$ subharmonic mode, as these are the easiest to excite in experiment for small driving frequency and related to translational drop motion. Figure 7 plots the instability tongues for a hemispherical 
$\alpha =90^\circ$ drop showing an overlap of the 
$[1,1]$ subharmonic mode with the 
$[2,0]$ harmonic mode. This occurs because the natural frequencies for these spatial modes are related by 
$2\omega _{1,1}\approx \omega _{2,0}$ for 
$\alpha =90^\circ$. This gives rise to the mixing mode, which we will discuss in detail in the next section when we describe our experimental results. We also note that the bandwidth for the instability tongues decreases with increasing slanting angle 
$\chi$, which makes sense given that the coefficient 
$\cos \chi$ on the parametric term also decreases with 
$\chi$.
Instability tongues for 
$[2,0]$ harmonic and 
$[1,1]$ subharmonic modes plotted in the acceleration–frequency (
$\varLambda - \omega$) space, depending upon the slanting angle 
$\chi =0^\circ$ and 
$60^\circ$, with pinned contact line and 
$\alpha =90^\circ$.
4. Experiment
Experiments were performed and compared with the theoretical predictions discussed previously. Given the recent interest in directional drop motion on substrates, we focus on characterizing the motions for the 
$[2,0]$ harmonic pumping mode and the 
$[1,1]$ subharmonic rocking mode for pinned contact-line conditions, as they depend upon the slanting angle 
$\chi$ and drop contact angle 
$\alpha$.
4.1. Experimental procedure
Experiments are performed using the set-up shown in figure 8. A droplet of distilled water is deposited onto a substrate through a needle using an AL-1000 syringe pump at low flow rate to induce pinch-off of a droplet of known volume 
$V$. The liquid properties of distilled water include the density 
$\rho =986.2\pm 2.4\,\textrm {kg}\,\textrm {m}^{-3}$, surface tension 
$\sigma = 71.06\pm 0.01\,\textrm {mN}\,\textrm {m}^{-1}$ and dynamic viscosity 
$\mu = 1\,\textrm {mPa}\,\textrm {s}$. The droplet contact line is pinned by a circular depression with depth 
$0.1\,\textrm {mm}$ and diameter 
$D=3.3\,\textrm {mm}$, which is machined into a solid piece of either aluminium or copper. The substrate is secured to a tilting angle table (TAT) (Sherline 3750-LAZ), which consists of two plates separated by angle 
$\chi$, the bottom of which is mounted on a Labworks ET-139 mechanical shaker. The shaker is driven by a Labworks power amplifier (PA-138) with feedback from an accelerometer (PCB J352C33) that is controlled by VibeLab Pro software (VL144x-7.09) to provide sinusoidal vibrations with acceleration 
$a$ and frequency 
$f$ at an angle 
$\chi \in [0^\circ, 90^\circ ]$ with limiting cases of vertical 
$\chi =0^\circ$ and horizontal 
$\chi =90^\circ$ vibration.
Experimental schematic shown in (a) front and (b) side views.
Droplet dynamics are captured by a Phantom VEO 410L high-speed camera at over 1000 fps. The camera is positioned on a three-axis stage and equipped with a LAOWA Ultra Macro 
$25\,\textrm {mm}$ lens. An optical diffuser (Edmund optics, 
$127\times 178\,\textrm {mm}$ sandblasted glass) is used to diffuse the backlight (MultiLED LT-V9-15). Both the shaker and the camera are placed on separate optic tables to minimize interference. The shaker rests on a 
$2.5\,\textrm {cm}$ thick vibration isolation pad for bench isolation. Image calibration is carried out through ImageJ using the outer diameter of the needle as the known length scale. An in-house MATLAB code is developed in which the drop image undergoes Gaussian filtering, a morphological closing operation using disk-shaped structuring elements, and finally binarization. The contact angle 
$\alpha$ is measured by fitting the interface shape to a second-order polynomial for the 50 closest points adjacent to the contact line, and computing the slope of the tangent line to give 
$\alpha$. Given 
$\alpha$ and the fixed wetted radius 
$r=D/2$, the drop volume can be determined by assuming a spherical cap shape to give 
$V$. The shape change dynamics are tracked using the binarized image from which the interface shape and centre of mass are readily computed. The resulting signal is subjected to a fast Fourier transform to extract the resonant frequency. Following each experiment, the droplet is removed using compressed air, and periodic cleaning of the substrate is performed using an ultrasonic cleaner (MGUC500).
4.2. Results
Two types of motion are observed in experiment: (1) low driving amplitude harmonic pure modes, and (2) high driving amplitude mixed modes. The latter involves the simultaneous excitation of a subharmonic mode and a harmonic mode at a single driving frequency. Here, we quantify the drop response for these motions, as it depends upon the experimental parameters.
4.2.1. Pumping and rocking modes
Pure 
$[2,0]$ pumping and 
$[1,1]$ rocking modes are excited at comparatively low driving acceleration (
$a=0.25\,g$ in this case), as shown in figures 9 and 10, respectively. The pumping mode is dominated by vertical motion, and the rocking mode by horizontal motion.
Harmonic pumping mode 
$[2,0]$ over a complete cycle of oscillation with period 
$T$, which has been excited at acceleration 
$a=0.25\,g$ and frequency 
$f=128\,\textrm {Hz}$ using vertical 
$\chi =0^\circ$ vibration, with contact angle 
$\alpha =60.7^\circ$ and volume 
$V=4.6$ ml.
Harmonic rocking mode 
$[1,1]$ over a complete cycle of oscillation with period 
$T$, which has been excited at acceleration 
$a=0.25\,g$ and frequency 
$f=63.5\,\textrm {Hz}$ using horizontal 
$\chi =90^\circ$ vibration, with contact angle 
$\alpha =60.7^\circ$ and volume 
$V=4.6$ ml.
4.2.2. Harmonic instability tongues
The harmonic rocking mode is isolated, as there is no adjacent mode to interact with. Here, the wave emerges gradually, making it experimentally challenging to identify the boundary of these tongues. This is in stark contrast to the threshold instability typically associated with the emergence of subharmonic Faraday waves. As such, we employ an artificial criterion to detect the harmonic drop response, which we use to present an equivalent harmonic instability tongue. The procedure that we use is to construct a frequency response curve by tracking the dynamic contact angles 
$\alpha _L,\alpha _R$ on the left and right sides of the drop during a frequency sweep with fixed driving acceleration, as shown in figure 11. This metric is more suitable for our experiments than tracking the centre of mass (COM) for two main reasons: (1) any rocking motion corresponds to a change in dynamic contact angles, and (2) this method focuses on accurately capturing the region near the contact points, whereas the calculation of the COM relies on the entire drop shape, making it susceptible to the influence of any defects in the interface shape due to light reflection. Figure 11(a) plots the time trace of the dynamic contact angles with the contact-angle difference 
$|\alpha _R-\alpha _L|$, an effective metric to quantify the drop response. Figure 11(b) plots the time trace of the contact-angle difference and identifies the local maxima (red dots) with associated average value (red solid line), as well as the average value of the overall signal (blue dashed line). This gives the average value representing the maximum difference in dynamic contact angles for that specific frequency, which we repeat across multiple frequencies (with maximum interval 
$2\,\textrm {Hz}$) to generate the frequency response curves shown in figure 11(c) for different slanting angles 
$\chi$. Each response curve exhibits a peak near the resonant frequency and an associated bandwidth, consistent with the resonance plots of Celestini & Kofman (Reference Celestini and Kofman2006). The bandwidth can be defined by choosing some fraction of the maximum drop response, for example 65 % shown as the green dashed line. We note that this choice is arbitrary and will change the corresponding value of the bandwidth, but it will not change the interpretation of the drop response if the same metric is used for all response curves. This procedure is repeated for various driving accelerations.
Procedure for obtaining the harmonic frequency response curve by measuring the time trace of (a) the dynamic contact angles at two sides of the drop, 
$\alpha _L,\alpha _R$, for 
$a= 0.25\,g$, 
$f=63\,\textrm {Hz}$ and 
$\chi = 60^\circ$. (b) The difference in contact angles 
$|\alpha _R-\alpha _L|$ is used as a metric for the frequency response, where the red dots denote the peak points of the signal, and the red line is the average value of these points, while the dashed blue line represents the overall average. (c) Frequency response curve plotting 
$|\alpha _R-\alpha _L|$ against driving frequency 
$f$ for two different slanting angles, 
$\chi = 30^\circ$ (squares) and 
$\chi = 60^\circ$ (circles). The green dashed line corresponds to 65 % of the peak value of the response curve at 
$\chi = 60^\circ$.
Figure 12 plots the bandwidth, computed using the procedure described in the previous paragraph, in the acceleration–frequency domain for two different slanting angles (
$\chi = 30^\circ$ and 
$60^\circ$) to produce an equivalent harmonic instability tongue. It is important to note that we have not observed a pronounced rocking mode below 
$a=0.05\,g$ and 
$a=0.1\,g$ for slanting angles 
$\chi = 60$ and 
$30^\circ$, respectively, for 
$\alpha = 60.7^\circ$ and 
$V=4.6$ ml. Below these critical accelerations, the size of the dynamic contact angle difference is comparable to the contact-angle hysteresis, and the periodic signature of the rocking mode is not observed in the corresponding time trace. Here, decreasing the slanting angle 
$\chi$ increases the threshold acceleration where the rocking mode is prominent, while also decreasing the bandwidth. This is because the horizontal component of the vibration force, which contributes to the rocking mode, decreases as the slanting angle decreases, shifting the response curves downwards (cf. figure 11c). Consequently, we observed no harmonic rocking mode at 
$\chi = 0^\circ$. Figure 6(d) shows that there exists a theoretical rocking mode at vertical forcing (
$\chi = 0^\circ$), but its bandwidth is so small that is difficult to observe in our experiments. Increasing the static contact angle 
$\alpha$ decreases the threshold acceleration and increases the bandwidth of the instability tongue, consistent with our theoretical results (cf. Figure 6d).
Harmonic instability tongue for the 
$[1,1]$ rocking mode with 
$\alpha = 60.7^\circ$ and 
$V=4.6$ ml for 
$\chi = 60^\circ$ (red) and 
$\chi =30^\circ$ (blue), and 
$\alpha = 74.2^\circ$ and 
$V=6.5$ ml for 
$\chi =30^\circ$ (green).
It is important to note that hysteresis is a factor in our measurements. To be precise, the measurements obtained during the sweeping process (over either frequency or acceleration) in the increasing direction differ from those obtained in the reverse direction due to the influence of wave size and nonlinearities (Benjamin & Ursell Reference Benjamin and Ursell1954). To mitigate this, we conduct each measurement using a new drop, with the sweep performed in the increasing direction.
4.2.3. Mixed mode
The close proximity of the instability tongues for the subharmonic 
$[1,1]$ rocking mode and the harmonic 
$[2,0]$ pumping mode leads to the mixed mode shown in figure 13. The observed subharmonic response implies that it takes time 
$2T$ for the drop to complete one cycle of oscillation, whereas the substrate oscillation completes one cycle in time 
$T$. A characteristic feature of the mixing mode is that the COM has both horizontal and vertical components, as shown in the time trace of figure 14(a), with the subharmonic behaviour attributed to the horizontal motion (i.e. rocking mode) of the drop COM, and the harmonic behaviour to the vertical motion (i.e. pumping mode) of the drop COM. This is more clearly illustrated by plotting the COM motion in the phase space, as shown in figure 14(b), which demonstrates that two cycles of vertical motion occur for every one cycle of horizontal motion. Note that the trajectory is not symmetric because the drop motion is biased in the forcing direction. The shape of this trajectory is expected to vary with the slanting angle, among other parameters, but such an exploration is beyond the scope of this paper.
Mixed mode excited with 
$a=2.2\,g$, 
$f=128\,\textrm {Hz}$, 
$\chi =60^\circ$, 
$\alpha =60.7^\circ$ and 
$V=4.6$ ml. The period of the substrate oscillation is 
$T$, whereas the period of the drop oscillation is 
$2T$.
The COM motion for the mixing mode plotted (a) as a time trace and (b) in the phase space, showing the trajectory corresponding to the experiment shown in figure 13.
The subharmonic instability tongues are obtained by conducting an acceleration sweep at a constant driving frequency. As the driving acceleration surpasses the threshold value, the COM of the drop experiences a sudden and pronounced horizontal oscillation corresponding to the emergence of the subharmonic 
$[1,1]$ rocking mode. This process is repeated over a range of driving frequencies. Figure 15(a) plots the instability tongues comparing experiment (dashed line type) to theory (solid line type), depending upon the slanting angle 
$\chi$. These curves shift upwards with increasing 
$\chi$, indicating that higher driving accelerations are required to excite the mixing mode. This trend is further highlighted in figure 15(b), which plots the normalized threshold acceleration 
$a_{th}/g$ for the mixed mode against 
$\chi$, showing an increasing trend. For the special case of horizontal forcing (
$\chi =90^\circ$), we did not observe a subharmonic 
$[1,1]$ rocking mode, which agrees well with our theoretical predictions (§ 3.2).
(a) Instability tongue for the mixed mode plotted in the acceleration–frequency (
$\varLambda - \omega$) space. Experimental data are given by the symbols. Theoretical predictions are shown in solid lines for the subharmonic 
$[1,1]$ rocking mode and harmonic 
$[2,0]$ pumping mode for 
$\alpha =60.7^\circ$ and 
$V=4.6$ ml with pinned contact line. (b) Threshold acceleration 
$a/g$ against slanting angle 
$\chi$ at the resonance frequency 
$f=128$ Hz.
Comparing experimental results with theoretical predictions, we observe that the right-hand side of the tongue shows good agreement with the subharmonic 
$[1,1]$ mode, whereas the left-hand side has deviations. In addition, theory does not predict an overlap of the instability tongues for the subharmonic 
$[1,1]$ rocking and harmonic 
$[2,0]$ pumping modes for these experimental conditions. This suggests that the harmonic 
$[2,0]$ pumping mode is not a Faraday wave, but rather a meniscus wave excited by low-amplitude vibrations. The interactions between this harmonic meniscus wave with the subharmonic rocking mode may manifest themselves through the shape of the instability tongue. Our model does not account for such interactions, which could explain the discrepancy with experiment. The presence of the meniscus wave is known to affect the shape of the Faraday wave instability tongue, as has been reported for waves in a cylindrical geometry by Batson, Zoueshtiagh & Narayanan (Reference Batson, Zoueshtiagh and Narayanan2013, figure 14). Low-amplitude meniscus waves have also been reported in the theoretical model of Maksymov & Pototsky (Reference Maksymov and Pototsky2019) for pancake-shaped sessile drops. These prior observations support our conclusion that the harmonic 
$[2,0]$ mode is indeed a meniscus wave.
5. Concluding remarks
We have derived a mathematical model of the parametric oscillations of the sessile drop subject to slanted vibrations, which gives rise to an infinite system of coupled Mathieu equations that is solved by Floquet theory. The limiting case of horizontal (
$\chi =90^\circ$) vibrations gives rise to a purely harmonic drop response, consistent with experimental observations (Sharp, Farmer & Kelly Reference Sharp, Farmer and Kelly2011; Khan & Eslamian Reference Khan and Eslamian2019), whereas the more general case 
$\chi \neq 90^\circ$ produces parametric oscillations defined by subharmonic and harmonic instability tongues in the driving frequency–amplitude space. We show how the shape of the instability tongues depends upon the static contact angle 
$\alpha$, slanting angle 
$\chi$, and contact-line boundary conditions.
Experiments have been performed to validate theoretical predictions focusing on the 
$[1,1]$ rocking mode, the 
$[2,0]$ pumping mode, and the mixed mode that is the superposition of a subharmonic rocking mode and harmonic pumping mode, the latter of which is highly relevant to directional drop motion (Brunet et al. Reference Brunet, Eggers and Deegan2007; Noblin et al. Reference Noblin, Kofman and Celestini2009; Guo et al. Reference Guo, Shui and Chen2020). The agreement between theoretical predictions and experimental observations is generally quite good, with the exception of the shape of the instability tongue for the mixed mode (cf. figure 15a). Here, the theoretical instability tongue for the harmonic 
$[2,0]$ pumping mode does not overlap with the instability tongue for the subharmonic 
$[1,1]$ rocking mode, as we both expect and observe experimentally for the mixing mode. We conjecture that this is because the harmonic pumping mode is not a true Faraday wave associated with a threshold instability, but rather an edge wave that can be excited at comparatively low driving acceleration. This has been observed experimentally for zonal modes in general (Vukasinovic, Smith & Glezer Reference Vukasinovic, Smith and Glezer2007a; Chang et al. Reference Chang, Bostwick, Daniel and Steen2015). Similar wave dynamics has been observed in Faraday waves in a cylindrical container, including axisymmetric harmonic edge waves excited at low driving amplitude (Shao et al. Reference Shao, Wilson, Bostwick and Saylor2021b), and mode mixing of a harmonic axisymmetric wave with a subharmonic asymmetric waves at a single driving frequency (Shao et al. Reference Shao, Wilson, Saylor and Bostwick2021c). Shao et al. (Reference Shao, Gabbard, Bostwick and Saylor2021a) conduct a set of experiments that clearly demonstrate that the presence of a meniscus at the container sidewall is responsible for generating harmonic edge waves, which are suppressed whenever the interface is perfectly flat. This observation can help to explain the reason why sessile drop zonal modes exhibit a harmonic response, and is perhaps not surprising given that a sessile drop necessarily has a meniscus.
Dissipation gives rise to a non-trivial threshold acceleration, which we do not predict given that our analysis is focused on inviscid fluids with either pinned or free contact-line boundary conditions. In the future, extensions to our model could be made to account for sources of dissipation, including bulk viscous dissipation (Prosperetti Reference Prosperetti1980; Strani & Sabetta Reference Strani and Sabetta1988; Bostwick & Steen Reference Bostwick and Steen2013) and contact-line dissipation (Davis Reference Davis1980; Bostwick & Steen Reference Bostwick and Steen2015). To incorporate the former into a model would require one to resolve the rotational component to the velocity field, while the latter can be accommodated through the Hocking (Reference Hocking1987) contact-line condition. Recent interest in complex fluids with a non-trivial rheology suggests model extensions to include viscoelastic effects (Tamim & Bostwick Reference Tamim and Bostwick2021), which are known to produce novel dynamics in classical Faraday wave systems (Müller & Zimmermann Reference Müller and Zimmermann1999).
Funding
J.B.B. acknowledges support from NSF grant CMMI-1935590.
Declaration of interests
The authors report no conflict of interest.
Open access
