On a toroidal method to solve the sessile-drop oscillation problem

Abstract We present a fully analytical solution for the natural oscillation of an inviscid sessile drop with small Bond number (surface tension dominates gravity) and a fixed contact line on a flat horizontal plate. The governing equations are expressed in terms of the toroidal coordinate system which yields solutions involving hypergeometric functions. Resonant frequencies are identified for zonal, sectoral and tesseral vibration modes. The predictions show excellent agreement with experimental data reported in the literature, particularly for flatter drops (lower $\theta _{c}$, but not so low as to incur significant viscous dissipation) and higher modes of vibration.

inviscid drop. While the free drop is generally assumed to be spherical, a sessile drop takes the form of a spherical cap when surface tension dominates gravity (i.e. √ γ /ρg c, where γ is the surface tension, ρ is the density and c is the contact radius of the drop). To find natural frequencies of the latter, analytical models in the literature either converted the geometry to a simplified form (replacing the planar substrate by a spherical one, Strani & Sabetta 1984) or developed a solution using spherical coordinates (Bostwick & Steen 2014). Although the former approach leads to a highly simplified physical model, the latter requires hybrid analytical-numerical schemes: neither is then suitably accurate and accessible for use by an experimentalist. , Steen, Chang & Bostwick (2016) have presented a detailed account of the underlying physics and mechanics of this problem. The contact angle θ c (and shape) of the drop is established at static equilibrium by balancing the liquid-gas, liquid-solid and solid-gas interfacial tensions. The drop stability is determined by the behaviour of the contact line (CL), via its speed u CL . Stick-slip behaviour of the CL (Shaikeea et al. 2017) gives rise to hysteresis, which is captured using a CL model. In the 'Hocking condition' presented by Davis (1980), contact-angle deviations are expressed in the form Δθ c ∝ u CL , with a constant of proportionality Λ which quantifies the CL resistance. This phenomenological parameter characterises the CL mobility; Λ = 0 corresponds to a fully mobile CL and Λ = ∞ to a pinned CL. In the current work, the toroidal framework imposes the pinned CL condition on the problem (see § 2.4).
We present here an analytical solution to this long-standing problem by using a toroidal coordinate system. The liquid-vapour and liquid-solid boundaries of a spherical cap, δD f and δD s (cf. figure 1a), correspond to a pair of β-coordinate curves in this system, where the boundary conditions can be directly expressed without any geometric conversions or complex computations. Solving the hydrodynamic equations in this framework requires the use of hypergeometric functions, which ultimately yield a fully analytical solution in the form of (2.18). The importance of choosing this framework to solve the sessile-drop evaporation problem was first presented by Popov (2005) and we believe this is the first time it has been extended to the oscillating sessile drop. Bostwick & Steen (2014) (hereafter referred to as Bo-St) presented a hybrid analytical-numerical model which solves the same problem and employs inverse operators to find the solution. Theirs is the most comprehensive investigation of the sessile-drop oscillation problem to date. They considered different types of vibration mode shapes, namely zonal, sectoral and tesseral, which were subsequently validated experimentally by Chang et al. (2015). In the current work, the resonant frequencies for the mode shapes discussed by Bo-St are calculated and compared with the experimental data reported in the literature.
The purpose of this work is to show that our model, based on toroidal coordinates, yields fully analytical solutions for the case of an inviscid drop with fixed CL in the shape of a spherical cap. Stating the hydrodynamic equations with boundary conditions, we perform an eigenmode analysis to find the solution ( § 2). This model is then used to identify resonant frequencies for zonal, sectoral and tesseral vibration modes ( § 3). Its predictions are compared with experimental data reported in the literature. Future work and possible extensions of this model are discussed in ( § 4).

Sessile-drop geometry
The liquid-vapour interface of a sessile drop with contact angle θ c ∈ (0, π) can be expressed in toroidal coordinates as r Figure 1. (a) Three-dimensional schematic of toroidal coordinate system r = r (α, β, ϕ) overlaid on a sessile drop. Based on Li, Kar & Kumar (2019). (b) Diametral section of the drop (blue shaded region) with toroidal gridlines embedded into it. On a red circle, β is constant, on a purple circle, α is constant. Defining expressions for α and β are also displayed.
varies along the surface ∂D f , where β ∈ [0, π] is the angle subtended by foci F 1 , F 2 on ∂D f and ϕ ∈ [0, 2π] varies in the azimuthal direction. The equilibrium (base) state Γ of the drop can be defined as (2.1) A small perturbation η (α, ϕ, t) on Γ (with the CL being fixed) leads to a competition between the drop's inertia and capillarity, and the resulting motion is oscillatory in nature (cf. figure 1b). These disturbances are often expressed in terms of where c is the drop contact radius and h α , h β and h ϕ are the scale factors of the toroidal system. Here the prime notation indicates that the variables are dimensional. In the following text, prime notation will be dropped from dimensionless variables, except for density ρ, surface tension γ and contact radius c of the drop. The scale factor quantifies the change in position of a point on changing one of its coordinates, so a Δα change in α (keeping other coordinates constant) corresponds to a change in distance alongα of h α Δα (cf. figure 1a).

Equations and boundary conditions
The flow is assumed to be incompressible and irrotational. The velocity field v is described as v = ∇ψ , where the velocity potential ψ satisfies Laplace's equation in the drop domain D. The equation becomes closed form when subject to the no-penetration condition at the substrate ∂D s , and a free-surface kinematic boundary condition at the interface ∂D f , where the normal velocity is set equal to the time derivative of perturbation. For an inviscid fluid, applying linear wave theory, the pressure field is described by the momentum equation (2.6) A small disturbance η to the equilibrium surface Γ causes a deviation from the initially spherical shape which is described by the modified Laplace equation where k 1 , k 2 are the principal curvatures, Δ T is the Laplace-Beltrami operator; definitions are given in Appendix A. In subsequent sections, we have replaced the term cosh α − cos β with b(α, β) while simplifying the terms involving scale factors h α , h β and h ϕ .

Curvatures and Laplace-Beltrami operator for toroidal coordinates
The first and second fundamental forms of a surface allow the calculation of curvature and Laplace-Beltrami operators, respectively, for a parametric surface x(u 1 , u 2 ). The coefficients for first fundamental form are given by the metric tensor (Kreyszig 1959). The derivation of the principal curvatures and the Laplace-Beltrami operator from the coefficients E, F and G is given in Appendix A.
An important consequence of (2.14) is that y = 0 at the CL, which arises because P → 0 as α → ∞, and so use of the toroidal coordinate system imposes the fixed CL condition on the problem. The mobility of the CL is defined as 1/Λ by Bo-St, which is zero for an immobile CL and infinite for a fully mobile CL. Only the immobile CL case, Λ = 0, is considered in the current work.
Substituting the above equation into (2.10d) gives, at β = β 0 , This can be rearranged as where I and II are, respectively, The term I is equivalent to (v 2 − 1 4 )P (see Lebedev 1965, p. 224). In fact, an analogous simplification is performed while deriving an expression for the eigenfrequencies of a free spherical drop in Rayleigh's derivation (see Landau & Lifshitz 1987, p. 246). Further simplification of the right-hand side of (2.16) gives −λ 2 = 2 sin 2 β 0 −b 2 (α, β 0 ) τ 2 + 1 4 where a single or double subscript α on a function denotes single or double derivative of the function with respect to α. The expressions for T α , T αα , P and P α (which fall under the class of hypergeometric functions) are given in Appendix B.

Results
The variation of dimensionless frequency λ with contact angle θ c = π − β 0 is determined by solving (2.18). Previous studies such as Bo-St classified the vibrational modes as zonal (l = 0), sectoral (τ = l) and tesseral (l / = 0, τ / = l). Results are presented for each type of mode in turn.
3.1. Zonal (l = 0) modes When the disturbance of the interface is axisymmetric, the mode shapes are termed zonal. For a sessile drop of fixed contact radius c, increasing the contact angle θ c increases the volume of the drop (inertia) and thus decreases the frequency λ (cf. figure 2a). There is excellent agreement between the model and the data of Chang et al. (2015) in this figure, particularly at higher mode numbers. For instance, for τ = 10 and θ c = 40 • , our model Further comparisons of predicted zonal mode frequencies with experimental measurements are shown for the data sets reported by Chang et al. (2013) and Mettu & Chaudhury (2012) in figures 2(b) and 2(c), respectively. In the former, the experimental values lie within the range of theoretical frequencies calculated for the range of contact angles θ c involved. For the higher modes, τ = 8 and 10, the frequencies lie at the upper end of the theoretical span, where there were a limited number of data points as these require larger droplets ( 5μL (see Mettu & Chaudhury 2012, figure 4a), whereas lower modes (τ = 2, 4, 6) were experimentally accessible for droplets with smaller volume, ≤5μL. In addition, there is a small increase in slope for τ = 10 (compared to slope at τ = 6); this feature is also present in figure 2(a) for θ c ≈ 65 • . In figure 2(c), the width of the predicted frequency band is small and lies at the lower end of the range of observed frequencies.
A possible explanation for this mismatch is that the model neglects contributions from viscous effects. Chang et al. (2013) reported that the bandwidth of predicted frequencies increased when viscous contributions were added (noting that the dimensional frequency is plotted here). Chang et al. (2015) subsequently showed that the viscous contribution is characterised by the Ohnesorge number, Oh = μ/ √ ρcγ , and even at a small value of Oh = 0.003 for water (instead of Oh = 0 for the inviscid case) the resonant peak changed from an infinite to a finite value and, thus, increased the bandwidth of predicted frequency (see Chang et al. 2015, p. 446). The effect of viscosity on a drop undergoing oscillations of arbitrary amplitude has been discussed both for free drops and sessile/pendant drops (see Basaran 1992 andBasaran 1997). For the latter case, it has been reported that as the viscosity increases, the resonant frequency also increases, so that excluding viscous effects can lead to predicted frequencies lying at the lower end of the observed spectrum.  Figure 3. Effect of contact angle on dimensionless frequency for sectoral modes, τ = l, for (a) [5,5], (b) [7,7] and (c) [9,9]. Solid loci show the solutions to (2.18), dashed loci are the results presented by Bostwick & Steen (2014). Symbols indicate experimental data reported by Chang et al. (2015).  Bostwick & Steen (2014). Symbols show experimental data reported by Chang et al. (2015). The shaded region in (d) represents the range of frequencies calculated using VPF theory by Chang et al. (2015) for water, with substrate forcing and viscosity included.
3.2. Sectoral (τ = l, l / = 0) modes A non-axisymmetric mode with wavenumber pair [τ, l] has l longitudinal intersections and (τ − l)/2 latitudinal intersections (or τ − l nodes on the interface) with the undisturbed interface Γ (see Bostwick & Steen 2014, p. 19). A sectoral mode, with τ = l, is a special case where there are only longitudinal intersections. Figure 3 compares the experimental frequencies reported by Chang et al. (2015) with our model and the Bo-St model. There is good agreement with our model for τ = 9. For τ = 5 and 7, the two models bracket the data.
3.3. Tesseral (τ / = l, l / = 0) modes A tesseral mode shape with wavenumber pair [τ, l] has non-zero longitudinal and latitudinal intersections because τ / = l. Figure 4 compares the results for our model and the Bo-St model in a similar fashion to the sectoral mode. For the τ = 9 cases, our model agrees with the experimental data quite well for all θ c values investigated. For θ c ≤ 65 • , the Bo-St model does not capture the observed trend, for 50 • and l = 7, it overpredicts by a factor of 1.17, whereas our model underpredicts slightly by a factor of 1.05. For τ = 7, there is good agreement between the experimental data and both models as θ c decreases from 140 • to 70 • , below which our model continues to perform well and Bo-St starts to diverge. For τ = 5, our model captures the frequencies at low θ c whereas the Bo-St model works well at higher values.
It should be noted that our model cannot predict the frequencies for small contact angles, because in this case a larger fraction of the sessile-drop volume lies within the solid-liquid boundary layer and drop-substrate interactions then cause damping of oscillations. The range of contact angles suggested to avoid boundary layer viscous dissipation effects, discussed in Sharp (2012), is 30 • -150 • .
For the lowest mode, [5,3], there is a slight overprediction for larger contact angles. This can be attributed to the assumption of a pinned CL in the current work. If the CL is, instead, assumed to be mobile and not pinned, the slope of the frequency-versus-contact-angle curve will decrease at larger contact angles (see Bostwick & Steen 2014, figures 10 and 11). This represents a limitation in the current model in that mobile CL behaviour is not readily incorporated in the toroidal coordinate approach.

Discussion and conclusions
The superior performance of our model for lower θ c and higher modes is probably the result of using toroidal coordinates, which fit the sessile drop naturally. An interesting physical insight from this work is that the slope of the λ versus θ c curve decreases as θ c decreases; the curve almost reaching a plateau. This is also suggested by experiments. The physical models present in literature incorporate bulk dissipation and CL (Davis) dissipation, e.g. , to account for this plateau. On the one hand, the current toroidal model, while established on zero viscosity and fixed CL assumptions, can still predict this plateau with good success. On the other hand, incorporating more dissipation terms will improve this model and bring more understanding of observations such as mode mixing and mode competition . Note that the strength of this inviscid theory coupled with an appropriate coordinate system points to the importance of choosing a framework which maps the complicated geometries of physics problems perfectly, as previously done in Fokas & Nachbin (2012) and Richardson (1992).
There are exceptions, e.g. figures 3(a) and 4(a), and we here consider whether the mismatch between the predictions of the model and the experimental data could arise from the assumptions made in our model. The larger error incurred by the Bo-St model can be attributed to the approach used to enforce the no-penetration condition. Earlier works on constrained drops (e.g. Ramalingam, Ramkrishna &Basaran 2012 andProsperetti 2012) essentially used the Lagrange multiplier method to enforce the no-penetration condition at the pinning circle and these methods permitted a discontinuity in the interface shape at the contact point. The Bo-St method does not allow a discontinuity at the pinning sites, which leads to overprediction of the frequency (see Bostwick & Steen 2015, p. 558).
The model considers the sessile drop on a substrate as a mass-spring system, where viscous effects and substrate-drop interactions are neglected. These assumptions were also made in the Bo-St model and were subsequently relaxed in their subsequent work, for example, , where they studied damping for viscous drops (with fixed CL) undergoing substrate-forced oscillation. To extend our model along the lines discussed in Chang et al. (2015), viscous contributions could be incorporated by adding a damping term, iλ C [y], to the right-hand side of (2.15), where C is the dissipation operator and = μ/(ρcγ ) 1/2 is the Ohnesorge number. Substrate-drop interactions can be modelled using two main assumptions: (i) constant contact radius; and (ii) modelling the substrate forcing via the bulk pressure in the drop in the form of Faraday oscillations. With regard to (i), it should be noted that contact-angle hysteresis on the modes cannot be incorporated using the toroidal framework presented here because it requires the incorporation of a dynamic CL condition (Bostwick & Steen 2014). With (ii), the substrate contribution is incorporated by adding a term F 0 e iλt to the right-hand side of (2.15), where λ is the frequency of substrate forcing (not the natural frequency) and F 0 its amplitude. Chang et al. (2015) used these assumptions and incorporated the aforementioned effects in their viscous potential flow (VPF) theory. The envelope of solutions which they obtained is shown as a shaded region in figure 4(d) and it spans Bo-St and our model (both inviscid). It is, thus, expected that the addition of viscous and substrate contributions to our model will modify (2.18) and increase the bandwidth of predicted frequencies. This is the subject of ongoing work, where the aim is to identify the contributions of viscous damping and substrate forcing, and thereby establish when significant differences will arise from the inviscid model.
The description of an oscillating drop presented here is not suitable for cases where the drop shape is influenced by gravity, which is quantified by Bo (ratio of gravity to surface tension). As the volume of the drop increases, the drop shape changes from that of a truncated sphere (Bo = 0), towards being ellipsoidal (0 < Bo < 5) until it forms a flat puddle (Bo > 5), with uniform depth except near the edges (Lubarda & Talke 2011). We believe that it should be possible to model a flattened drop using confocal ellipsoidal coordinates system in the 0.5 < Bo < 5 regime. Finding resonant frequencies for a flattened drop will allow us to extend our Bo = 0 theory to 0 < Bo < 5, and compare the results with the theory for flattened drops presented by Noblin, Buguin & Brochard-Wyart (2004) where the drop is modelled as a liquid bath and the resonant frequency is that associated with a standing wave on its interface. This work introduces, for the first time, an analytical solution to the sessile-drop oscillation problem. The superiority of this model lies in the fact that its predictions work well for lower contact angles (<75 • ) compared with the existing models. It also predicts a decrease in slope as θ c decreases, which is consistent with experiments. The behaviour at lower contact angles (<30 • ) remains to be experimentally validated (and physically understood) for all types of modes.
To summarise, our model provides a concise solution to the sessile-drop oscillation problem which opens a new window to researchers interested in this and related problems. A clear next step could be to test the θ c < 30 • regime experimentally, model a drop being vibrated on an inclined plane (see Brunet, Eggers & Deegan 2007) and extend the model to larger drops by including the effects of gravity.