2.1. First-order equations

To obtain the first-order equations with their boundary and initial conditions, the above series expansions are introduced into the respective equations, and all the terms with the parameter $\eta _0$ to the first power are collected. The first-order continuity and momentum equations read

(2.15)$$\begin{gather} \frac{1}{r^{2}} \frac{\partial}{\partial r} ( r^{2} u_{r1} ) + \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} ( u_{\theta 1} \sin \theta ) = 0 , \end{gather}$$

(2.16)$$\begin{gather}\nu_{0}\frac{\partial u_{r1}}{\partial \tau} ={-} \frac{\partial p_1}{\partial r} , \end{gather}$$

(2.17)$$\begin{gather}\nu_{0}\frac{\partial u_{\theta 1}}{\partial \tau} ={-} \frac{1}{r} \frac{\partial p_1}{\partial \theta}. \end{gather}$$

For the kinematic and normal-stress boundary conditions of first order at $r = 1$ we obtain

(2.18)$$\begin{gather} u_{r1} = \nu_{0}\frac{\partial \eta_1}{\partial \tau}, \end{gather}$$

(2.19)$$\begin{gather}-p_1 - \left( 2 \eta_1 + \frac{\partial \eta_1}{\partial \theta} \cot \theta + \frac{\partial^{2} \eta_1}{\partial \theta^{2}} \right) = 0, \end{gather}$$

respectively. Furthermore, the initial conditions of first order are

(2.20a,b)\begin{equation} \eta_1(\theta,0) = P_m (\cos\theta) \quad \mbox{and} \quad \nu_{0}\frac{\partial \eta_1}{\partial \tau} (\theta,0) = 0.\end{equation}

The first initial condition determines the initial shape of the deformed drop, governed by a Legendre polynomial of order $m$.

2.2. Second-order equations

The second-order equations with their boundary and initial conditions are obtained as all the terms from the expansion with the parameter $\eta _0$ to the second power. The second-order continuity and momentum equations read

(2.21)$$\begin{gather} \frac{1}{r^{2}} \frac{\partial}{\partial r} ( r^{2} u_{r2} ) + \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} ( u_{\theta 2} \sin \theta ) = 0, \end{gather}$$

(2.22)$$\begin{gather}\nu_{0} \frac{\partial u_{r2}}{\partial \tau} + \frac{\partial p_2}{\partial r} = \frac{u_{\theta 1}^{2}}{r} - \nu_{1} \frac{\partial u_{r1}}{\partial \tau} - u_{r1} \frac{\partial u_{r1}}{\partial r} - \frac{u_{\theta 1}}{r}\frac{\partial u_{r1}}{\partial \theta}, \end{gather}$$

(2.23)$$\begin{gather}\nu_{0}\frac{\partial u_{\theta 2}}{\partial \tau} + \frac{1}{r} \frac{\partial p_2}{\partial \theta} ={-} \nu_{1}\frac{\partial u_{\theta 1}}{\partial \tau} - u_{r1} \frac{\partial u_{\theta 1}}{\partial r} - \frac{u_{\theta 1}}{r} \frac{\partial u_{\theta 1}}{\partial \theta} - \frac{u_{r1} u_{\theta 1}}{r}. \end{gather}$$

The kinematic and zero normal-stress boundary conditions of second order to be satisfied at $r = 1$ are

(2.24)$$\begin{gather} u_{r2} - \nu_{0}\frac{\partial \eta_2}{\partial \tau} = \nu_{1}\frac{\partial \eta_1}{\partial \tau} + \frac{u_{\theta 1}}{r} \frac{\partial \eta_1}{\partial \theta} - \eta_1 \frac{\partial u_{r1}}{\partial r}, \end{gather}$$

(2.25)$$\begin{gather}-p_2 - \left( 2 \eta_2 + \frac{\partial \eta_2}{\partial \theta} \cot \theta + \frac{\partial^{2} \eta_2}{\partial \theta^{2}} \right) = \eta_1 \frac{\partial p_1}{\partial r} - 2 \eta_1^{2} - 2 \eta_1 \frac{\partial \eta_1}{\partial \theta} \cot \theta - 2 \eta_1 \frac{\partial ^{2} \eta_1}{\partial \theta^{2}}. \end{gather}$$

The initial conditions of second order are

(2.26a,b)\begin{equation} \eta_2(\theta,0) ={-} \frac{1}{2m + 1} \quad \mbox{and} \quad \nu_{0}\frac{\partial \eta_2}{\partial \tau} (\theta,0) + \nu_{1}\frac{\partial \eta_1}{\partial \tau} (\theta,0) = 0.\end{equation}

The first initial condition ensures volume conservation for the deformed drop to second order. Solving the sets of first- and second-order equations together with their boundary and initial conditions will reveal the nonlinear effect of unequal time spent in the prolate and oblate shapes as a function of the initial deformation amplitude. These solutions reveal the frequency decrease only for initial deformation modes $m > 2$. To find the frequency decrease for the two-lobed mode $m=2$, the theory must be taken to third order, as will be shown below.

2.3. Third-order equations

The third-order continuity and momentum equations read

(2.27) \begin{gather} \frac{1}{r^{2}} \frac{\partial}{\partial r} ( r^{2} u_{r3} ) + \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} ( u_{\theta 3} \sin \theta ) = 0,\\\nonumber\end{gather}

(2.28) \begin{align} &\nu_{0}\frac{\partial u_{r 3}}{\partial \tau} + \frac{\partial p_{3}}{\partial r} ={-} \nu_{1}\frac{\partial u_{r 2}}{\partial \tau} - \nu_{2}\frac{\partial u_{r 1}}{\partial \tau} + 2\frac{u_{\theta 1}u_{\theta 2}}{r} \nonumber\\ & \quad\qquad\qquad\qquad - u_{r1}\frac{\partial u_{r2}}{\partial r} - u_{r2}\frac{\partial u_{r1}}{\partial r} - \frac{u_{\theta 1}}{r}\frac{\partial u_{r2}}{\partial \theta} - \frac{u_{\theta 2}}{r}\frac{\partial u_{r1}}{\partial \theta}, \\ &\nu_{0}\frac{\partial u_{\theta 3}}{\partial \tau} + \frac{1}{r}\frac{\partial p_{3}}{\partial \theta} ={-} \nu_{1}\frac{\partial u_{\theta 2}}{\partial \tau} - \nu_{2}\frac{\partial u_{\theta 1}}{\partial \tau} -\frac{u_{r1}u_{\theta 2}}{r} - \frac{u_{r2}u_{\theta 1}}{r} \nonumber\\ &\qquad\qquad\qquad\qquad\!\! - u_{r1}\frac{\partial u_{\theta 2}}{\partial r} - u_{r2}\frac{\partial u_{\theta 1}}{\partial r} - \frac{u_{\theta 1}}{r}\frac{\partial u_{\theta 2}}{\partial \theta} - \frac{u_{\theta 2}}{r}\frac{\partial u_{\theta 1}}{\partial \theta}. \end{align}

The third-order kinematic and zero normal-stress boundary conditions, which must be satisfied at $r=1$, are

(2.29) \begin{align} & u_{r 3} - \nu_{0}\frac{\partial \eta_{3}}{\partial \tau} = \nu_{1}\frac{\partial \eta_{2}}{\partial \tau} + \nu_{2}\frac{\partial \eta_{1}}{\partial \tau} + \frac{u_{\theta 1}}{r} \frac{\partial \eta_2}{\partial \theta} + \frac{u_{\theta 2}}{r} \frac{\partial \eta_1}{\partial \theta} + \frac{\partial}{\partial r}\left(\frac{u_{\theta 1}}{r}\right)\frac{\partial \eta_{1}}{\partial \theta}\eta_{1} \nonumber\\ &\qquad\qquad\qquad - \frac{\partial u_{r2}}{\partial r}\eta_1 - \frac{\partial u_{r1}}{\partial r} \eta_2 - \frac{1}{2}\frac{\partial^{2} u_{r1}}{\partial r^{2}}\eta_{1}^{2},\end{align}

(2.30)\begin{align} &-p_{3} - \left( 2 \eta_3 + \frac{\partial \eta_3}{\partial \theta} \cot \theta + \frac{\partial^{2} \eta_3}{\partial \theta^{2}} \right) = \frac{\partial p_{1}}{\partial r}\eta_{2} + \frac{\partial p_{2}}{\partial r}\eta_{1} + \frac{1}{2}\frac{\partial^{2}p_{1}}{\partial r^{2}}\eta_{1}^{2} \nonumber\\ &\quad - \eta_{1}^{3} - 2\eta_{1}\left(\eta_{2} + \frac{\partial \eta_{2}}{\partial \theta} \cot \theta + \frac{\partial^{2} \eta_{2}}{\partial \theta^{2}} \right) - 2\eta_{2}\left(\eta_{1} + \frac{\partial \eta_{1}}{\partial \theta} \cot \theta + \frac{\partial^{2} \eta_{1}}{\partial \theta^{2}} \right) \nonumber\\ &\quad - \frac{1}{2}\left(\frac{\partial\eta_{1}}{\partial\theta}\right)^{2}\left(\frac{\partial \eta_{1}}{\partial\theta}\cot\theta + 3\frac{\partial^{2}\eta_{1}}{\partial\theta^{2}}\right) + 3\eta_{1}^{2}\left(\eta_{1} + \frac{\partial\eta_{1}}{\partial\theta}\cot\theta + \frac{\partial^{2}\eta_{1}}{\partial\theta^{2}}\right), \end{align}

respectively. Furthermore, the third-order initial conditions are

(2.32a,b)\begin{align} \eta_3(\theta,0) ={-}\frac{1}{3}\begin{pmatrix} m & m & m \\ 0 & 0 & 0 \end{pmatrix}^{2} \quad \mbox{and} \quad \left(\nu_{0}\frac{\partial \eta_3}{\partial \tau}+\nu_{1}\frac{\partial \eta_2}{\partial \tau}+\nu_{2}\frac{\partial \eta_1}{\partial \tau}\right) (\theta,0) = 0,\end{align}

where the first initial condition is described in terms of the 3-$j$ symbols (see the supplementary material available at https://doi.org/10.1017/jfm.2021.568). Solving the sets of equations up to third order together with their boundary and initial conditions will reveal the weakly nonlinear dependency of the angular oscillation frequency on the amplitude of the initial drop deformation. This effect becomes visible due to the strained time introduced for suppressing secular terms in solutions.

3.1. First-order solutions

The first-order equations describe the linear problem. The well known solutions in Rayleigh (Reference Rayleigh1879) and Lamb (Reference Lamb1881) for the inviscid drop in a vacuum are expected to be recovered. Since we investigate two-dimensional flow fields, the Stokesian stream function is applied for determining the first-order velocity field. We set the coefficient $\nu _{0}$ in the strained time to unity. The stream function $\psi (r,\theta ,\tau )$ is defined by its relations to the two velocity components $u_{r1}$ and $u_{\theta 1}$ as (Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot1962)

(3.1a,b)\begin{equation} u_{r 1} ={-} \frac{1}{r^{2} \sin \theta} \frac{\partial \psi}{\partial \theta} \quad \mbox{and} \quad u_{\theta 1} = \frac{1}{r \sin \theta} \frac{\partial \psi}{\partial r}. \end{equation}

These first-order velocity components satisfy the continuity equation identically.

The first-order drop surface deformation is governed by the Legendre polynomial of the initial deformation. The solution is therefore sought in the form

(3.2)\begin{equation} \eta_1 = {\hat \eta}_1 P_m (\cos\theta) \, \textrm{e}^{- \alpha_m \tau}, \end{equation}

with the first-order initial surface amplitude ${\hat \eta }_1$ and the first-order angular frequency $\alpha _m$ for the deformation mode $m$.

Taking the curl of the vectorial first-order momentum equation with the components (2.16) and (2.17), the partial differential equation

(3.3)\begin{equation} \frac{\partial^{2} \psi}{\partial r^{2}} + \frac{\sin\theta}{r^{2}}\frac{\partial}{\partial\theta} \left( \frac{1}{\sin\theta} \frac{\partial \psi}{\partial \theta} \right) = 0, \end{equation}

for the stream function is obtained. Its solution, which is regular at $r=0$, reads (Brenn Reference Brenn2017)

(3.4)\begin{equation} \psi = C_{1m} r^{m+1} \sin^{2} \theta ~P'_m (\cos\theta) \, \textrm{e}^{-\alpha_m \tau}, \end{equation}

where the prime denotes the derivative of the Legendre polynomial $P_m$ with respect to its argument. The radial and angular components of the first-order velocity vector follow as

(3.5)\begin{equation} u_{r 1} ={-} C_{1m} r^{m-1} m(m+1) P_m (\cos\theta) \, \textrm{e}^{-\alpha_{m} \tau} \end{equation}

and

(3.6)\begin{equation} u_{\theta 1} ={-} C_{1m} r^{m-1} (m+1) P_{m}^{1} (\cos\theta) \, \textrm{e}^{-\alpha_{m} \tau}, \end{equation}

respectively. The function $P_{m}^{1}$ is a Legendre function of degree $m$ and order one. The integration constant $C_{1m}$ is determined by the first-order kinematic boundary condition as

(3.7)\begin{equation} C_{1m} = \frac{\alpha_{m} {\hat \eta}_{1}}{m(m+1)}.\end{equation}

The first-order pressure field is obtained by integration of one component of the momentum equation as

(3.8)\begin{equation} p_1 ={-} C_{1m} \alpha_{m} (m+1) r^{m} P_m (\cos\theta) \, \textrm{e}^{-\alpha_{m} \tau}. \end{equation}

The characteristic equation of the drop, which determines the complex angular frequency $\alpha _m$, follows from the first-order zero normal stress boundary condition (2.19), as obtained in Rayleigh (Reference Rayleigh1879). Its complex conjugate solutions are

(3.9a,b)\begin{align} \alpha_{m}^{(p)} = i \sqrt{m(m-1)(m+2)} =: i\, \alpha_{m,0}, \qquad \alpha_{m}^{(n)} ={-}i \sqrt{m(m-1)(m+2)} =: -i\, \alpha_{m,0},\end{align}

where the superscripts $(p)$ and $(n)$ represent the positive and negative solutions, respectively. The existence of pairs of solutions of the characteristic equation requires representation in the first-order solutions. We therefore formulate the first-order deformation of the drop surface as

(3.10)\begin{equation} \eta_1 (\theta, \tau) = ( {\hat \eta}_1^{(p)} \, \textrm{e}^{- \alpha_m^{(p)} \tau} + {\hat \eta}_1^{(n)} \, \textrm{e}^{- \alpha_m^{(n)} \tau} ) P_m (\cos\theta). \end{equation}

The first-order initial conditions (2.20a,b) require that initially the drop surface shape is governed by the function $P_m (\cos\theta )$, and that the surface is at rest. These conditions reveal the deformation amplitudes

(3.11a,b)\begin{equation} {\hat \eta}_1^{(p)} ={-} \frac{\alpha_m^{(n)}}{\alpha_m^{(p)} - \alpha_m^{(n)}} = \frac{1}{2} \quad \mbox{ and} \quad {\hat \eta}_1^{(n)} = \frac{\alpha_m^{(p)}}{\alpha_m^{(p)} - \alpha_m^{(n)}} = \frac{1}{2}.\end{equation}

The first-order radial and polar velocity components are represented as

(3.12)\begin{equation} u_{r 1}(r, \theta, \tau)={-} (C_{1m}^{(p)}\, \textrm{e}^{-\alpha_{m}^{(p)} \tau} + C_{1m}^{(n)}\, \textrm{e}^{-\alpha_{m}^{(n)} \tau}) r^{m-1} m(m+1) P_m (\cos\theta) \end{equation}

and

(3.13)\begin{equation} u_{\theta 1}(r, \theta, \tau) ={-} (C_{1m}^{(p)}\, \textrm{e}^{-\alpha_{m}^{(p)} \tau} + C_{1m}^{(n)}\, \textrm{e}^{-\alpha_{m}^{(n)} \tau}) r^{m-1} (m+1) P_{m}^{1} (\cos\theta), \end{equation}

respectively. The first-order pressure field is

(3.14)\begin{equation} p_1 (r, \theta, \tau)={-} (C_{1m}^{(p)} \alpha_{m}^{(p)} \, \textrm{e}^{-\alpha_{m}^{(p)} \tau} + C_{1m}^{(n)} \alpha_{m}^{(n)} \, \textrm{e}^{-\alpha_{m}^{(n)} \tau}) (m+1) r^{m} P_m (\cos\theta) . \end{equation}

With account for the complex conjugate frequencies in (3.9a,b), the first-order drop surface shape, velocity components and pressure field simplify to the forms

(3.15)$$\begin{gather} u_{r 1}(r, \theta, \tau) ={-}\, \alpha_{m,0} r^{m-1} P_m\, (\cos\theta) \sin (\alpha_{m,0} \tau), \end{gather}$$

(3.16)$$\begin{gather}u_{\theta 1}(r, \theta, \tau) ={-} \frac{\alpha_{m,0}}{m} r^{m-1} P_{m}^{1} (\cos\theta) \sin (\alpha_{m,0} \tau), \end{gather}$$

(3.17)$$\begin{gather}p_{1} (r, \theta, \tau) = (m-1)(m+2) r^{m} P_{m} (\cos\theta) \cos (\alpha_{m,0} \tau), \end{gather}$$

(3.18)$$\begin{gather}\eta_1 (\theta, \tau) = P_m (\cos\theta) \cos (\alpha_{m,0} \tau). \end{gather}$$

3.2. Second-order solutions

We start developing the second-order solutions from pressure. The general solution for the second-order pressure is sought as a sum

(3.19)\begin{equation} p_2 (r,\theta,\tau) = {p}_{21}(r,\theta,\tau) + {p}_{22}(r,\theta,\tau), \end{equation}

of two contributions, where subscript ‘21’ indicates the solution of the system of second-order equations including the products of first-order terms, and subscript ‘22’ the solutions of the homogeneous system. In the second-order approximation, the coefficient $\nu _{1}$ in the strained time is unknown and needs to be determined.

The solutions with subscripts ‘21’ are first determined. It is convenient to form the divergence of the vectorial second-order momentum equation with the components (2.22) and (2.23) to eliminate the second-order velocity and the terms exhibiting $\nu _{1}$ with use of the continuity equation (2.21). This yields the differential equation for the second-order pressure $p_{21}$,

(3.20)\begin{align} &\frac{1}{r^{2}} \frac{\partial}{\partial r} \left( r^{2} \frac{\partial p_{21}}{ \partial r} \right) + \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} \left( \frac{\sin \theta}{r} \frac{\partial p_{21}}{ \partial \theta} \right) \nonumber\\ &\quad ={-} \frac{1}{r^{2}} \frac{\partial}{\partial r} \left[ r^{2} \left( u_{r1} \frac{\partial u_{r1}}{\partial r} + \frac{u_{\theta 1}}{r} \frac{\partial u_{r1}}{\partial \theta} - \frac{u_{\theta 1}^{2}}{r} \right) \right] \nonumber\\ &\qquad - \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} \left[ \left( u_{r1} \frac{\partial u_{\theta 1}}{\partial r} + \frac{u_{\theta 1}}{r} \frac{\partial u_{\theta 1}}{\partial \theta} + \frac{u_{r1} u_{\theta 1}}{r} \right) \sin \theta \right], \end{align}

which we rewrite, using the Lamé identity and accounting for the curl-free first-order velocity field, in the form

(3.21)\begin{equation} \varDelta [ p_{21} + {\boldsymbol v}_1^{2} / 2 ] = 0 ,\end{equation}

which states that the modified pressure ${\mathcal {P}}_{21} = p_{21} + {\boldsymbol v}_1^{2}/2$ is potential. The squared first-order velocity vector determines the time dependency of pressure $p_{21}$.

The solution of the Laplace equation for the modified pressure ${\mathcal {P}}_{21}$, formulated in the contribution $p_{21}$ to the second-order pressure, accounting for all the time dependencies arising from the square of the first-order velocity vector may be written as

(3.22)\begin{align} p_{21}(r,\theta,\tau) & ={-} \sum _{l=0}^{L} (C_{21l}^{(p)}\, \textrm{e}^{{-}2\alpha_{m}^{(p)}\tau} + C_{21l}^{(n)}\, \textrm{e}^{{-}2\alpha_{m}^{(n)}\tau} + C_{21l}^{(pn)}\, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{m}^{(n)}) \tau} ) r^{l}P_{l}(\cos\theta) \nonumber\\ &\quad - \frac{1}{2}(u_{r1}^{2} + u_{\theta 1}^{2}). \end{align}

The coefficients $C_{21l}^{(pn)}$ relating to time-independent terms in $p_{21}$ are set to zero for all $l$. This treatment of these coefficients is justified after the calculation of the radial velocity in the next step. All the summations arising in the solutions henceforth extend to an upper summation index symbolised by an uppercase letter. Rewriting the radial momentum equation (2.22) with the modified pressure for the contribution ‘21’, and knowing that the first-order velocity field is irrotational, one obtains

(3.23)\begin{equation} \frac{\partial u_{r21}}{\partial \tau} ={-} \frac{\partial {\mathcal{P}}_{21}}{\partial r} - \nu_{1} \frac{\partial u_{r1}}{\partial \tau}. \end{equation}

From this equation, the radial velocity component $u_{r21}$ is calculated by integration with respect to time, revealing

(3.24)\begin{equation} u_{r21}(r,\theta,\tau) ={-} \sum _{l=0}^{L} \left(\frac{C_{21l}^{(p)}}{2 \alpha_{m}^{(p)}}\, \textrm{e}^{{-}2\alpha_{m}^{(p)} \tau} + \frac{C_{21l}^{(n)}}{2 \alpha_{m}^{(n)}} \, \textrm{e}^{{-}2\alpha_{m}^{(n)} \tau}\right) l r^{l-1} P_{l}(\cos\theta) - \nu_{1} u_{r1}. \end{equation}

Non-zero coefficients $C_{21l}^{(pn)}$ in the pressure $p_{21}$ would have produced linearly time-dependent secular terms, which would grow in time. Setting the coefficients $C_{21l}^{(pn)}$ for all $l$ to zero in (3.22) prevents this. The velocity component $u_{\theta 21}$ is then obtained from the continuity equation and reads

(3.25)\begin{equation} u_{\theta 21}(r,\theta,\tau ) ={-}\sum _{l=0}^{L}\left(\frac{C_{21l}^{(p)}}{2 \alpha_{m}^{(p)}} \, \textrm{e}^{{-}2\alpha_{m}^{(p)} \tau} + \frac{C_{21l}^{(n)}}{2 \alpha_{m}^{(n)}} \, \textrm{e}^{{-}2\alpha_{m}^{(n)} \tau}\right) r^{l-1}P_{l}^{1}(\cos\theta) - \nu_{1} u_{\theta 1}. \end{equation}

Following the work by Tsamopoulos & Brown (Reference Tsamopoulos and Brown1983, Reference Tsamopoulos and Brown1984), the surface deformation $\eta _{21}$ is written in a general form as

(3.26)\begin{equation} \eta_{21}(\theta,\tau ) = \sum_{l=0}^{L} \delta_{2l}(\tau) P_{l}(\cos\theta) , \end{equation}

where the time dependency is governed by the unknown functions $\delta _{2l} (\tau )$, and the deformed shape by all possible Legendre polynomials. The reason for the latter is that, in the zero normal-stress boundary condition, pressure and surface deformation are coupled. The functions $\delta _{2l}(\tau )$ are determined for every $l$ using the boundary conditions (2.24) and (2.25). For doing this, $u_{r21}$ and $\eta _{21}$ from (3.24) and (3.26), respectively, are substituted into the kinematic boundary condition (2.24). Its derivative with respect to the strained time $\tau$ is taken. We substitute the pressure $p_{21}$ and $\eta _{21}$, from (3.22) and (3.26), respectively, into the zero normal-stress boundary condition (2.25). Since we deal with the boundary conditions, $r=1$ in both equations. We proceed by multiplying the two equations with a Legendre polynomial $P_{j}$ and integrating the results between the south and north poles of the drop. Multiplying the equation resulting from the zero-normal stress boundary condition by the natural number $j$ and subtracting the two resulting equations eliminates the coefficients $C_{21j}$ for all time dependencies, yielding the ordinary differential equation

(3.27)\begin{align} & {\ddot \delta}_{2j}(\tau) + j(\,j - 1)(\,j + 2)\delta_{2j}(\tau) =\left[ - \nu_{1} ({\alpha_{m}^{(p)}}^{2} \, \textrm{e}^{-\alpha_{m}^{(p)} \tau} + {\alpha_{m}^{(n)}}^{2} \, \textrm{e}^{-\alpha_{m}^{(n)} \tau}) \int_{{-}1}^{1} P_{m} P_{j}\, {\textrm{d}\kern0.7pt x} \right. \nonumber\\ &\quad - \left. \int_{{-}1}^{1} \frac{\partial }{\partial \tau}\left( \frac{u_{\theta 1}}{r} \frac{\partial \eta_1}{\partial \theta} - \eta_1 \frac{\partial u_{r1}}{\partial r}\right) \right|_{r=1} P_{j}\, {\textrm{d}\kern0.7pt x} \nonumber\\ &\quad + \int_{{-}1}^{1} \left.\left.\left(Rz_{21} - \frac{1}{2}(u_{r1}^{2} + u_{\theta 1}^{2})\right) \right|_{r=1} j P_{j}\, {\textrm{d}\kern0.7pt x} \right] \frac{2j + 1}{2} , \end{align}

for $\delta _{2j}(\tau )$ (Tsamopoulos & Brown Reference Tsamopoulos and Brown1983, Reference Tsamopoulos and Brown1984), where $Rz_{21}$ is the right-hand side of the zero normal-stress boundary condition (2.25). For the differential equation (3.27), we are interested in the particular solutions only, because the homogeneous solutions will be accounted for by the second-order contribution ‘22’. The particular solutions need to be found for every $j$ between zero and $L$. Introducing back all particular solutions $\delta _{2j}(\tau )$ of (3.27) into (3.26), we can write the final form as

(3.28)\begin{align} \eta_{21}(\theta, \tau ) &= \sum _{j=0}^{L} (H_{21j}^{(p)}\, \textrm{e}^{{-}2\alpha_{m}^{(p)} \tau} + H_{21j}^{(n)}\, \textrm{e}^{{-}2\alpha_{m}^{(n)} \tau} + H_{21j}^{(pn)}\, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{m}^{(n)}) \tau}) P_{j}(\cos\theta) \nonumber\\ &\quad + \frac{\nu_{1}}{4} [ ( 1 + 2 \alpha_{m}^{(p)} \tau ) \, \textrm{e}^{-\alpha_{m}^{(p)} \tau} + (1+2\alpha_{m}^{(n)}\tau )\, \textrm{e}^{-\alpha_{m}^{(n)} \tau}] P_{m}, \end{align}

where the coefficients $H_{21j}^{(p)}$, $H_{21j}^{(n)}$ and $H_{21j}^{(pn)}$ are known. The second part of the solution (3.28), with $\nu _1$ in front, exhibits secular terms growing in time. Due to this unphysical behaviour, these terms must be suppressed by setting $\nu _{1} = 0$. With the surface deformation $\eta _{21}$ known, the coefficients $C_{21l}$ in $u_{r21}$, $u_{\theta 21}$ and $p_{21}$ are readily determined, for every time dependency and for all $l$ between zero and $L$, using one of the boundary conditions (2.24) or (2.25). The details of the calculation of $C_{21l}$ are presented in the supplementary material. Now, we can write the contributions ‘$21$’ to the second-order solutions in the form

(3.29)$$\begin{gather} u_{r21}(r,\theta,\tau) = \sum _{l=0}^{L} \frac{C_{21l}}{\alpha_{m,0}} \sin(2 \alpha_{m,0} \tau ) l r^{l-1} P_{l}(\cos\theta), \end{gather}$$

(3.30)$$\begin{gather}u_{\theta 21}(r,\theta,\tau) = \sum _{l=0}^{L} \frac{C_{21l}}{\alpha_{m,0}} \sin(2 \alpha_{m,0} \tau ) r^{l-1}P_{l}^{1}(\cos\theta), \end{gather}$$

(3.31)$$\begin{gather}p_{21}(r,\theta,\tau) ={-} 2 \sum _{l=0}^{L} C_{21l} \cos(2 \alpha_{m,0} \tau ) r^{l}P_{l}(\cos\theta) - \frac{1}{2}(u_{r1}^{2} + u_{\theta 1}^{2}), \end{gather}$$

(3.32)$$\begin{gather}\eta_{21}(\theta,\tau) = \sum _{l=0}^{L} (2 H_{21l} \cos(2 \alpha_{m,0} \tau ) + H_{21l}^{(pn)}) P_{l}(\cos\theta). \end{gather}$$

The second contributions to the second-order solutions, with subscript ‘22’, are determined from the homogeneous forms of (2.21)–(2.23). This set of equations has the same structure as at first order, where the complex angular frequency $\alpha _{m}$ is replaced by $\alpha _{2k}$ (with the number ‘2’ in the subscript indicating the second-order solution, and $k$ a deformation mode number) and the deformation amplitude of the drop surface $\hat {\eta }_{1}$ by $\hat {\eta }_{22k}$. Both replaced variables appear under a sum with the summation index $k$, implying that all possible deformation modes $k$ affect the results.

The contributions ‘22’ to the second-order velocity, pressure and deformation are obtained as linear combinations of eigensolutions of the underlying differential equations with the summation index $k$. Their forms are similar to the first-order solutions and read

(3.33)\begin{align} u_{r22}(r,\theta,\tau ) & ={-}\sum _{k=0}^{K}(C_{22k}^{(p)}\, \textrm{e}^{-\alpha_{2k}^{(p)}\tau} + C_{22k}^{(n)}\, \textrm{e}^{-\alpha_{2k}^{(n)}\tau}) \nonumber\\ &\quad \times k(k+1)r^{k-1}P_{k}(\cos\theta) ,\\ u_{\theta 22}(r,\theta, \tau) & ={-}\sum _{k=0}^{K}(C_{22k}^{(p)}\, \textrm{e}^{-\alpha_{2k}^{(p)}\tau} + C_{22k}^{(n)}\, \textrm{e}^{-\alpha_{2k}^{(n)}\tau}) \nonumber\\ &\quad \times (k+1) r^{k-1}P_{k}^{1}(\cos\theta) ,\\ p_{22}(r,\theta,\tau) &={-}\sum _{k=0}^{K}(\alpha_{2k}^{(p)}C_{22k}^{(p)}\, \textrm{e}^{-\alpha_{2k}^{(p)}\tau} + \alpha_{2k}^{(n)} C_{22k}^{(n)}\, \textrm{e}^{-\alpha_{2k}^{(n)}\tau}) \nonumber\\ &\quad \times (k+1)r^{k}P_{k}(\cos\theta), \end{align}

(3.34)\begin{align} \eta_{22}(\theta,\tau) &= \sum _{k=0}^{K}(\hat{\eta}_{22k}^{(p)}\, \textrm{e}^{-\alpha_{2k}^{(p)}\tau} + \hat{\eta}_{22k}^{(n)} \, \textrm{e}^{-\alpha_{2k}^{(n)}\tau}) P_{k}(\cos\theta). \end{align}

The two new complex conjugate angular frequencies $\alpha _{2k}$, obtained as solutions of the homogeneous zero normal-stress boundary condition, for all summation indices $k$ take the forms obtained in Rayleigh (Reference Rayleigh1879):

(3.37a,b)\begin{align} \alpha_{2k}^{(p)} = i \sqrt{k(k-1)(k+2)} =: i\, \alpha_{2k,0}, \qquad \alpha_{2k}^{(n)} ={-}i \sqrt{k(k-1)(k+2)} =: -i\, \alpha_{2k,0}.\end{align}

The calculations of the coefficients $C_{22k}$ and of the amplitudes $\hat {\eta }_{22k}$ are detailed in the supplementary material.

The contributions ‘22’ to the second-order solutions are obtained in the simplified forms

(3.38)$$\begin{gather} u_{r22}(r,\theta,\tau ) ={-} 2 \sum _{k=0}^{K} {\hat \eta}_{22k} \alpha_{2k,0} r^{k-1}P_{k} (\cos\theta) \sin(\alpha_{2k,0} \tau ) , \end{gather}$$

(3.39)$$\begin{gather}u_{\theta 22}(r,\theta, \tau) ={-} 2 \sum _{k=1}^{K} \frac{{\hat \eta}_{22k} \alpha_{2k,0}}{k} r^{k-1}P_{k}^{1}(\cos\theta) \sin(\alpha_{2k,0} \tau ) , \end{gather}$$

(3.40)$$\begin{gather}p_{22}(r,\theta,\tau) = 2 \sum _{k=0}^{K}{\hat \eta}_{22k} (k - 1)(k + 2) r^{k}P_{k}(\cos\theta) \cos(\alpha_{2k,0} \tau), \end{gather}$$

(3.41)$$\begin{gather}\eta_{22}(\theta,\tau) = 2 \sum _{k=0}^{K}{\hat \eta}_{22k} P_{k}(\cos\theta) \cos(\alpha_{2k,0} \tau ). \end{gather}$$

3.3. Third-order solutions

We start developing the third-order solutions also from pressure. The general solution for the third-order pressure is sought as a sum

(3.42)\begin{equation} p_{3}(r,\theta,\tau) = p_{31}(r,\theta,\tau) + p_{32}(r,\theta,\tau), \end{equation}

of two contributions, where subscript ‘31’ indicates the solution of the third-order equations system including the products of first- and second-order terms, while the subscript ‘32’ denotes the solutions of the homogeneous system of the governing equations with its boundary and initial conditions. The coefficient $\nu _{2}$ in the strained time coordinate is unknown and needs to be determined.

The solutions with the subscript ‘31’ are first determined. Taking the divergence of the vectorial third-order momentum equation with the components (2.28) and (2.29), and using the continuity equation (2.27), the third-order velocities are eliminated. The result is the following partial differential equation for the third-order pressure $p_{31}$:

(3.43) \begin{align} & \frac{1}{r^{2}} \frac{\partial}{\partial r} \left( r^{2} \frac{\partial p_{31}}{ \partial r} \right) + \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} \left( \frac{\sin \theta}{r} \frac{\partial p_{31}}{ \partial \theta} \right) \nonumber\\ & \quad ={-} \frac{1}{r^{2}} \frac{\partial}{\partial r} \left[ r^{2} \left( u_{r1} \frac{\partial u_{r2}}{\partial r} + u_{r2} \frac{\partial u_{r1}}{\partial r} + \frac{u_{\theta 1}}{r} \frac{\partial u_{r2}}{\partial \theta} +\frac{u_{\theta 2}}{r} \frac{\partial u_{r1}}{\partial \theta} - 2\frac{u_{\theta 1} u_{\theta 2}}{r} \right) \right] \nonumber\\ &\qquad - \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} \left[ \left( u_{r1} \frac{\partial u_{\theta 2}}{\partial r} + u_{r2} \frac{\partial u_{\theta 1}}{\partial r} + \frac{u_{\theta 2}}{r} \frac{\partial u_{\theta 1}}{\partial \theta} + \frac{u_{\theta 1}}{r} \frac{\partial u_{\theta 2}}{\partial \theta} + \frac{u_{r2} u_{\theta 1}}{r} + \frac{u_{r1} u_{\theta 2}}{r} \right) \sin \theta \right], \end{align}

which we write in a compact form as

(3.44)\begin{equation} \varDelta p_{31} ={-} \textrm{div} [(\boldsymbol{u}_{1}\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u}_{2} + (\boldsymbol{u}_{2}\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u}_{1}]. \end{equation}

The equation is rewritten using the dot product rule to yield

(3.45)\begin{equation} \varDelta {\mathcal{P}}_{31} = \textrm{div} [\boldsymbol{u}_{1}\times(\boldsymbol{\nabla}\times\boldsymbol{u}_{2}) + \boldsymbol{u}_{2}\times(\boldsymbol{\nabla}\times\boldsymbol{u}_{1})], \end{equation}

where the modified pressure is defined as ${\mathcal {P}}_{31} = p_{31} + \boldsymbol {u}_{1}\boldsymbol {\cdot } \boldsymbol {u}_{2}$. Evidently the right-hand side is zero, since the velocity fields involved are irrotational. We determine the time dependencies of the third-order pressure $p_{31}$ following the same procedure as for the second-order modified pressure. Accounting for all the time dependencies arising from the inner product of the first- and second-order velocity fields we write

(3.46)\begin{align} p_{31}(r,\theta,\tau) &={-}\sum_{h=0}^{H} r^{h}P_{h}(x) \left( C_{31h}^{(p)} \, \textrm{e}^{{-}3\alpha_{m}^{(p)}\tau} + C_{31h}^{(n)} \, \textrm{e}^{{-}3\alpha_{m}^{(n)}\tau} \vphantom{\sum_{k=0}^{K}} \right.\nonumber\\ &\quad + C_{31h}^{(ppn)} \, \textrm{e}^{-(2\alpha_{m}^{(p)} + \alpha_{m}^{(n)})\tau} + C_{31h}^{(pnn)} \, \textrm{e}^{-(\alpha_{m}^{(p)} + 2\alpha_{m}^{(n)})\tau} \nonumber\\ &\quad + \sum_{k=0}^{K} C_{31hk}^{(pp)}\, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{2k}^{(p)})\tau} + \sum_{k=0}^{K} C_{31hk}^{(nn)}\, \textrm{e}^{-(\alpha_{m}^{(n)} + \alpha_{2k}^{(n)})\tau} \nonumber\\ &\quad + \left.\sum_{k=0}^{K} C_{31hk}^{(pn)}\, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{2k}^{(n)})\tau} + \sum_{k=0}^{K} C_{31hk}^{(np)}\, \textrm{e}^{-(\alpha_{m}^{(n)} + \alpha_{2k}^{(p)})\tau} \right) - \boldsymbol{u}_{1}\boldsymbol{\cdot} \boldsymbol{u}_{2} . \end{align}

The coefficients $C_{31hm}^{(pn)}$ and $C_{31hm}^{(np)}$ relating to time-independent terms in $p_{31}$ are set to zero for all $h$, following the same reasoning as in the second-order approximation. The contributions ‘31’ to the velocity field are determined by the same procedure as for second order. Rewriting the radial momentum equation (2.28), using that the first- and second-order velocity fields are irrotational, and with the definition of the third-order modified pressure ${\mathcal {P}}_{31}$, one obtains

(3.47)\begin{equation} \frac{\partial u_{r31}}{\partial \tau} ={-}\frac{\partial{\mathcal{P}}_{31}}{\partial r} - \nu_{2} \frac{\partial u_{r1}}{\partial \tau}. \end{equation}

Integration of this momentum equation with respect to time, with the modified pressure ${\mathcal {P}}_{31}$ known, yields the radial velocity component $u_{31}$ for every time dependency

(3.48)\begin{align} u_{r31}(r,\theta,\tau) &={-} \sum _{h=0}^{H} \left( \frac{C_{31h}^{(p)}}{3\alpha_{m}^{(p)}} \, \textrm{e}^{{-}3\alpha_{m}^{(p)}\tau} + \frac{C_{31h}^{(n)}}{3\alpha_{m}^{(n)}} \, \textrm{e}^{{-}3\alpha_{m}^{(n)}\tau} + \frac{C_{31h}^{(ppn)}}{2\alpha_{m}^{(p)} + \alpha_{m}^{(n)}} \, \textrm{e}^{-(2\alpha_{m}^{(p)}+\alpha_{m}^{(n)})\tau} \right. \nonumber\\ &\quad + \frac{C_{31h}^{(pnn)}}{\alpha_{m}^{(p)} + 2\alpha_{m}^{(n)}} \, \textrm{e}^{-(\alpha_{m}^{(p)} + 2\alpha_{m}^{(n)})\tau} + \sum_{k=0}^{K} \frac{C_{31hk}^{(pp)}}{\alpha_{m}^{(p)} + \alpha_{2k}^{(p)}} \, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{2k}^{(p)})\tau} \nonumber\\ &\quad + \sum_{k=0}^{K} \frac{C_{31hk}^{(nn)}}{\alpha_{m}^{(n)} + \alpha_{2k}^{(n)}} \, \textrm{e}^{-(\alpha_{m}^{(n)} + \alpha_{2k}^{(n)})\tau} + \sum_{k=0\backslash m}^{K} \frac{C_{31hk}^{(pn)}}{\alpha_{m}^{(p)} + \alpha_{2k}^{(n)}} \, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{2k}^{(n)})\tau} \nonumber\\ &\quad + \left.\sum_{k=0\backslash m}^{K} \frac{C_{31hk}^{(np)}}{\alpha_{m}^{(n)} + \alpha_{2k}^{(p)}} \, \textrm{e}^{-(\alpha_{m}^{(n)} + \alpha_{2k}^{(p)})\tau}\right)h r^{h-1} P_{h}(\cos\theta) - \nu_{2} u_{r1}. \end{align}

Non-zero coefficients $C_{31hm}^{(pn)}$ and $C_{31hm}^{(np)}$ in the pressure $p_{31}$ would have produced linearly time-dependent secular terms, which would grow in time, just as seen at second order. These secular terms are suppressed by this treatment. The polar velocity component is then determined using the continuity equation (2.27) and reads

(3.49)\begin{align} u_{\theta 31}(r,\theta,\tau) &={-} \sum _{h=0}^{H} \left( \frac{C_{31h}^{(p)}}{3\alpha_{m}^{(p)}} \, \textrm{e}^{{-}3\alpha_{m}^{(p)}\tau} + \frac{C_{31h}^{(n)}}{3\alpha_{m}^{(n)}} \, \textrm{e}^{{-}3\alpha_{m}^{(n)}\tau} + \frac{C_{31h}^{(ppn)}}{2\alpha_{m}^{(p)} + \alpha_{m}^{(n)}} \, \textrm{e}^{-(2\alpha_{m}^{(p)}+\alpha_{m}^{(n)})\tau} \right. \nonumber\\ &\quad + \frac{C_{31h}^{(pnn)}}{\alpha_{m}^{(p)} + 2\alpha_{m}^{(n)}} \, \textrm{e}^{-(\alpha_{m}^{(p)} + 2\alpha_{m}^{(n)})\tau} + \sum_{k=0}^{K} \frac{C_{31hk}^{(pp)}}{\alpha_{m}^{(p)} + \alpha_{2k}^{(p)}} \, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{2k}^{(p)})\tau} \nonumber\\ &\quad + \sum_{k=0}^{K} \frac{C_{31hk}^{(nn)}}{\alpha_{m}^{(n)} + \alpha_{2k}^{(n)}} \, \textrm{e}^{-(\alpha_{m}^{(n)} + \alpha_{2k}^{(n)})\tau} + \sum_{k=0\backslash m}^{K} \frac{C_{31hk}^{(pn)}}{\alpha_{m}^{(p)} + \alpha_{2k}^{(n)}} \, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{2k}^{(n)})\tau} \nonumber\\ &\quad + \left.\sum_{k=0\backslash m}^{K} \frac{C_{31hk}^{(np)}}{\alpha_{m}^{(n)} + \alpha_{2k}^{(p)}} \, \textrm{e}^{-(\alpha_{m}^{(n)} + \alpha_{2k}^{(p)})\tau}\right) r^{h-1} P_{h}^{1}(\cos\theta) - \nu_{2} u_{\theta 1}. \end{align}

For formulating the surface deformation $\eta _{31}$, we proceed as in the second-order approximation. Following the work by Tsamopoulos & Brown (Reference Tsamopoulos and Brown1983, Reference Tsamopoulos and Brown1984), the surface deformation is written in a general form as

(3.50)\begin{equation} \eta_{31}(\theta,\tau ) = \sum_{h=0}^{H} \delta_{3h}(\tau) P_{h}(\cos\theta), \end{equation}

where the time dependency is governed by the unknown functions $\delta _{3h} (\tau )$, and the deformed shape by all possible Legendre polynomials. We determine the functions $\delta _{3h} (\tau )$ by deriving the governing ordinary differential equation, following the same procedure as in the second-order approximation. The procedure eliminates the coefficients $C_{31j}$ and $C_{31jk}$ for all time dependencies from the kinematic and zero normal-stress boundary conditions at third order, (2.30) and (2.31), yielding the ordinary differential equation (Tsamopoulos & Brown Reference Tsamopoulos and Brown1983, Reference Tsamopoulos and Brown1984)

(3.51)\begin{align} & {\ddot \delta}_{3j}(\tau) + j(\,j - 1)(\,j + 2)\delta_{3j}(\tau) = \frac{2j + 1}{2} \left[\left. - \int_{{-}1}^{1} \frac{\partial Rk_{31}}{\partial \tau}\right|_{r=1} P_{j}\, {\textrm{d}\kern0.7pt x} \right.\nonumber\\ &\quad - \left.\left.\frac{\nu_{2}}{2} ({\alpha_{m}^{(p)}}^{2} \, \textrm{e}^{-\alpha_{m}^{(p)} \tau} + {\alpha_{m}^{(n)}}^{2} \, \textrm{e}^{-\alpha_{m}^{(n)} \tau}) \int_{{-}1}^{1} P_{m} P_{j}\, {\textrm{d}\kern0.7pt x} + \int_{{-}1}^{1} (Rz_{31} - \boldsymbol{u}_1 \boldsymbol{\cdot} \boldsymbol{u}_2)\right|_{r=1} j P_{j}\, {\textrm{d}\kern0.7pt x}\right], \end{align}

for $\delta _{3j} (\tau )$, where $Rk_{31}$ and $Rz_{31}$ are the right-hand sides of the two boundary conditions (2.30) and (2.31), respectively. We keep this compact formulation of the equation, since the right-hand sides of the boundary conditions consist of many terms. The calculations and symbolic derivations were performed with the software Mathematica. We are only interested in the particular solutions of (3.51). The general solution of the homogeneous equation will be accounted for by the contribution ‘32’, which solves the homogeneous system of equations. The solutions need to be found for every $j$ between zero and $H$. Introducing back all the particular solutions $\delta _{3j}(\tau )$ of the differential equation (3.51) into (3.50), we can write the final form as

(3.52)\begin{align} \eta_{31}(\theta,\tau) &= \sum _{j=0}^{H} \left( H_{31j}^{(p)}\, \textrm{e}^{{-}3\alpha_{m}^{(p)}\tau} + H_{31j}^{(n)}\, \textrm{e}^{{-}3\alpha_{m}^{(n)}\tau} \vphantom{\sum_{k=0}^{K}}\right. \nonumber\\ &\quad + H_{31j}^{(ppn)}\, \textrm{e}^{-(2\alpha_{m}^{(p)} + \alpha_{m}^{(n)})\tau} + H_{31j}^{(pnn)}\, \textrm{e}^{-(\alpha_{m}^{(p)} + 2\alpha_{m}^{(n)})\tau} \nonumber\\ &\quad + \sum_{k=0}^{K} H_{31jk}^{(pp)}\, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{2k}^{(p)})\tau} + \sum_{k=0}^{K} H_{31jk}^{(nn)}\, \textrm{e}^{-(\alpha_{m}^{(n)} + \alpha_{2k}^{(n)})\tau} \nonumber\\ &\quad + \left.\sum_{k=0}^{K} H_{31jk}^{(pn)}\, \textrm{e}^{-(\alpha_{m}^{(p)} + \alpha_{2k}^{(n)})\tau} + \sum_{k=0}^{K} H_{31jk}^{(np)}\, \textrm{e}^{-(\alpha_{m}^{(n)} + \alpha_{2k}^{(p)})\tau} \right) P_{j}(\cos\theta), \end{align}

where the coefficients $H_{31j}$ and $H_{31jk}$ for all time dependencies are known. Among them, the coefficients shown in table 1 are set to zero, since the corresponding solutions exhibit secular terms. These solutions must therefore be suppressed as unphysical, and their structure would not correspond to those represented by the above equation. The phenomenon arises because the initial deformation of the mode $m$ excites higher-order modes $i$, which can be interpreted as a resonance. The secular terms appear in the three cases where

(3.53) \begin{equation} \left. \begin{array}{c@{}} i(i-1)(i+2)/\alpha_{m}^{2} \\ i(i-1)(i+2)/4\alpha_{m}^{2} \\ i(i-1)(i+2)/9\alpha_{m}^{2} \end{array} \right\} ={-}1. \end{equation}

Equation (3.51) then assumes the general form

(3.54)\begin{equation} \ddot{\delta}_{3i}^{*}(\tau) - c^{2} \delta_{3i}^{*} (\tau) = A \, \textrm{e}^{{-}c \tau}, \end{equation}

with the particular solution

(3.55)\begin{equation} \delta_{3i}^{*} (\tau) ={-}\frac{A}{4c^{2}} (1+2 c \tau) \, \textrm{e}^{{-}c\tau}, \end{equation}

where the notation with the asterisk is introduced to indicate the exceptional cases of (3.51) and its solutions. The solutions exhibit in front of the exponential function a linear secular term which grows in time. Secular solutions relating to the first part of the statement (3.53), for which $i=m$, are due to functions of time on the right-hand sides of (3.51) leading to resonance with the frequency $\alpha _m^{(p)}$ or $\alpha _m^{(n)}$ appearing in the exponentials. They contain the coefficient $\nu _2$ in a linear combination with other, known, quantities represented by the factor $A$ in (3.55). For suppressing all these secular solutions, the values of the coefficient $\nu _{2}$ are determined for each mode $m$ of initial drop deformation such that the solutions vanish. This results in a straining of time, which may be seen equivalent to a change of the oscillation frequency against the first-order motion. The corresponding coefficients are all those listed in table 1 for the various modes $m$, except $H_{314}^{(p)}$ and $H_{314}^{(n)}$.

Table 1.

The coefficients $H_{31j}$ and $H_{31jk}$ in the surface deformation $\eta _{31}$ for $m = 2,3$ and $4$ set to zero since they relate to secular solutions.

The second part of the statement (3.53) appears only for $m=5$ when $i=8$, and for $m=10$ when $i=16$, corresponding to solutions with the coefficients $H_{3185}^{(pp)}$, $H_{3185}^{(nn)}$ and $H_{311610}^{(pp)}$, $H_{311610}^{(nn)}$, respectively. These cases are not treated here, since our analysis is restricted to modes less than $m=5$. The third part of the statement appears only for $m=2$ when $i=4$, with the corresponding two coefficients $H_{314}^{(p)}$ and $H_{314}^{(n)}$. The secular solutions $\delta _{34}^{*}(\tau )$ are due to functions of time on the right-hand sides of (3.51) leading to resonance with the frequency $3\alpha _m^{(p)}$ or $3\alpha _m^{(n)}$ appearing in the exponentials. In these secular solutions $\delta _{34}^{*}(\tau )$, the coefficient $\nu _2$ does not occur, so that it cannot be used for suppressing them. Therefore, the coefficients $H_{314}^{(p)}$ and $H_{314}^{(n)}$ corresponding to this case are set to zero to suppress this secular solution.

With the surface deformation $\eta _{31}$ known, the coefficients $C_{31j}$ and $C_{31jk}$ for all time dependencies are calculated using one of the boundary conditions (2.30) or (2.31). The procedure of this coefficient calculation is identical to the second-order approximation. We omit the details of this derivation because it involves large expressions. The final representations of the pressure (3.46), radial velocity (3.48) and polar velocity (3.49) ‘31’ are given in simplified forms as

(3.56) \begin{align} p_{31}(r,\theta,\tau) &={-} 2 \sum_{h=0}^{H} r^{h}P_{h}(x) \left[ C_{31h}^{(1)} \cos(3 \alpha_{m,0} \tau) + C_{31h}^{(2)}\cos(\alpha_{m,0} \tau) \vphantom{\sum_{k=0}^{K}}\right. \nonumber\\ &\quad + \sum_{k=0}^{K} C_{31hk}^{(3)} \cos((\alpha_{2k,0} + \alpha_{m,0})\tau) \nonumber\\ &\quad + \left.\sum_{k=0 \backslash m}^{K} C_{31hk}^{(4)} \cos((\alpha_{2k,0} - \alpha_{m,0}) \tau )\right] - \boldsymbol{u}_{1}\boldsymbol{\cdot} \boldsymbol{u}_{2} ,\end{align}

(3.57) \begin{align} u_{r31}(r,\theta,\tau) &= 2 \sum _{h=0}^{H} h r^{h-1} P_{h}(\cos\theta) \left[\frac{C_{31h}^{(1)}}{3 \alpha_{m,0}}\sin(3\alpha_{m,0} \tau) + \frac{C_{31h}^{(2)}}{\alpha_{m,0}}\sin (\alpha_{m,0}\tau)\right. \nonumber\\ &\quad + \sum_{k=0}^{K} \frac{C_{31hk}^{(3)}}{\alpha_{2k,0} + \alpha_{m,0}}\sin((\alpha_{2k,0} + \alpha_{m,0})\tau) \nonumber\\ &\quad + \left. \sum_{k=0\backslash m}^{K} \frac{C_{31hk}^{(4)}}{\alpha_{2k,0} - \alpha_{m,0}} \sin((\alpha_{2k,0} - \alpha_{m,0})\tau) \right] - \nu_{2}u_{r1}, \end{align}

(3.58) \begin{align} u_{\theta 31}(r,\theta,\tau) & = 2 \sum _{h=0}^{H} r^{h-1} P_{h}^{1}(\cos\theta) \left[\frac{C_{31h}^{(1)}}{3 \alpha_{m,0}}\sin(3\alpha_{m,0} \tau) + \frac{C_{31h}^{(2)}}{\alpha_{m,0}}\sin (\alpha_{m,0}\tau)\right. \nonumber\\ &\quad + \sum_{k=0}^{K} \frac{C_{31hk}^{(3)}}{\alpha_{2k,0} + \alpha_{m,0}}\sin((\alpha_{2k,0} + \alpha_{m,0})\tau)\nonumber\\ &\quad + \left.\sum_{k=0\backslash m}^{K} \frac{C_{31hk}^{(4)}}{\alpha_{2k,0} - \alpha_{m,0}} \sin((\alpha_{2k,0} - \alpha_{m,0})\tau) \right] - \nu_{2}u_{\theta 1}. \end{align}

The related drop surface deformation reads

(3.59)\begin{align} \eta_{31}(\theta,\tau) & = 2 \sum _{h=0}^{H} P_{h}(\cos\theta) \left[ H_{31h}^{(1)}\cos(3\alpha_{m,0} \tau ) + H_{31h}^{(2)}\cos(\alpha_{m,0}\tau) \vphantom{\sum_{k=0}^{K}} \right.\nonumber\\ &\quad + \left.\sum_{k=0}^{K} H_{31hk}^{(3)}\cos((\alpha_{2k,0} + \alpha_{m,0})\tau) + \sum_{k=0}^{K} H_{31hk}^{(4)}\cos((\alpha_{2k,0} - \alpha_{m,0})\tau) \right], \end{align}

where the definitions of the coefficients for the different time dependencies are given in the supplementary material.

Having determined the contributions ‘31’ to the third-order solutions, we move to the contributions ‘32.’ The homogeneous system of the third-order governing equations together with their boundary conditions determine these contributions. The homogeneous system consists of the third-order governing equations (2.27)–(2.29). Its solutions must satisfy the third-order homogeneous boundary conditions (2.30) and (2.31). The solutions are linear combinations of eigensolutions with index $h$ which are similar to the first order, and the contributions ‘$22$’ to the second-order solutions. The velocity components, the pressure field and the deformation of the drop surface read

(3.60) \begin{align} u_{r32}(r,\theta,\tau) & ={-}\sum _{h=0}^{H}(C_{32k}^{(p)}\, \textrm{e}^{-\alpha_{3h}^{(p)}\tau} + C_{32h}^{(n)}\, \textrm{e}^{-\alpha_{3h}^{(n)}\tau}) \nonumber\\ &\quad \times h(h+1)r^{h-1}P_{h}(\cos\theta), \end{align}

(3.61)\begin{align} u_{\theta 32}(r,\theta,\tau) &={-}\sum _{h=0}^{H}(C_{32h}^{(p)}\, \textrm{e}^{-\alpha_{3h}^{(p)}\tau} + C_{32h}^{(n)}\, \textrm{e}^{-\alpha_{3h}^{(n)}\tau}) \nonumber\\ &\quad \times (h+1)r^{h-1}P_{h}^{1}(\cos\theta), \end{align}

(3.62)\begin{align} p_{32}(r,\theta,\tau) & ={-}\sum _{h=0}^{H}(\alpha_{3h}^{(p)}C_{32h}^{(p)}\, \textrm{e}^{-\alpha_{3h}^{(p)}\tau} + \alpha_{3h}^{(n)} C_{32h}^{(n)}\, \textrm{e}^{-\alpha_{3h}^{(n)}\tau}) \nonumber\\ &\quad \times (h+1) r^{h}P_{h}(\cos\theta), \end{align}

(3.63)\begin{align} \eta_{32}(\theta,\tau) &= \sum _{h=0}^{H}(\hat{\eta}_{32h}^{(p)}\, \textrm{e}^{-\alpha_{3h}^{(p)}\tau} + \hat{\eta}_{32h}^{(n)}\, \textrm{e}^{-\alpha_{3h}^{(n)}\tau})P_{h}(\cos\theta). \end{align}

Two new complex conjugate angular frequencies $\alpha _{3h}$, which are determined from the homogeneous third-order zero normal-stress boundary condition, take the form obtained in Rayleigh (Reference Rayleigh1879). They are determined for all the summation indices $h$ as

(3.64a,b)\begin{align} \alpha_{3h}^{(p)} = i \sqrt{h(h-1)(h+2)} =: i\, \alpha_{3h,0}, \qquad \alpha_{3h}^{(n)} ={-}i \sqrt{h(h-1)(h+2)} =: -i\, \alpha_{3h,0}.\end{align}

The calculations of the coefficients $C_{32h}$ and of the amplitudes $\hat {\eta }_{32h}$ are detailed in the supplementary material.

The contributions ‘32’ to the third-order solutions are obtained in the simplified forms

(3.65)\begin{gather} u_{r32}(r,\theta,\tau ) ={-} 2 \sum _{h=0}^{H} {\hat \eta}_{32h} \alpha_{3h,0} r^{h-1}P_{h} (\cos\theta) \sin(\alpha_{3h,0} \tau ), \end{gather}

(3.66)\begin{gather}u_{\theta 32}(r,\theta, \tau) ={-} 2 \sum _{h=1}^{H} \frac{{\hat \eta}_{32h} \alpha_{3h,0}}{h} r^{h-1}P_{h}^{1}(\cos\theta) \sin(\alpha_{3h,0} \tau ), \end{gather}

(3.67)\begin{gather}p_{32}(r,\theta,\tau) = 2 \sum _{h=0}^{H}{\hat \eta}_{32h} (h - 1)(h + 2) r^{h}P_{h}(\cos\theta) \cos(\alpha_{3h,0} \tau), \end{gather}

(3.68)\begin{gather}\eta_{32}(\theta,\tau) = 2 \sum _{h=0}^{H}{\hat \eta}_{32h} P_{h}(\cos\theta) \cos(\alpha_{3h,0} \tau ). \end{gather}

The solutions developed in this section represent the response of an inviscid oscillating drop to an initial deformation of its surface from the spherical shape governed by one pure deformation mode.