2 Boussinesq–Scriven surface viscosity

Under the Boussinesq–Scriven model of surface viscosity (Scriven Reference Scriven1960; Aris Reference Aris1963; Slattery, Sagis & Oh Reference Slattery, Sagis and Oh2007), the surface stress boundary conditions at the interface between two Newtonian fluids can be written as

(2.1) $$\begin{eqnarray}[\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}]=\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{s},\end{eqnarray}$$

where $\unicode[STIX]{x1D64F}$ is the viscous stress tensor, $\unicode[STIX]{x1D735}_{\!s}=\unicode[STIX]{x1D64B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the surface gradient operator for the projection tensor $\unicode[STIX]{x1D64B}=\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n}^{T}$ with normal vector $\boldsymbol{n}$ , $[\cdot ]$ denotes the jump in magnitude across the interface and $\unicode[STIX]{x1D748}_{s}$ is the surface viscous stress tensor, defined by

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D748}_{s}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D64B}+(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{s})(\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }\boldsymbol{u}_{s})\unicode[STIX]{x1D64B}+2\unicode[STIX]{x1D707}_{s}\unicode[STIX]{x1D63F}_{s},\end{eqnarray}$$

where

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D63F}_{s}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D64B}\boldsymbol{ : }\unicode[STIX]{x1D735}_{\!s}\boldsymbol{u}_{s}+(\unicode[STIX]{x1D735}\boldsymbol{u}_{s})^{T}\boldsymbol{ : }\unicode[STIX]{x1D64B})\end{eqnarray}$$

is the surface rate of deformation tensor. The divergence of $\unicode[STIX]{x1D748}_{s}$ may be written (Scriven Reference Scriven1960) in the form

(2.4) $$\begin{eqnarray}\displaystyle [\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}] & = & \displaystyle \unicode[STIX]{x1D735}_{\!s}\unicode[STIX]{x1D70E}+(\unicode[STIX]{x1D707}_{d}+\unicode[STIX]{x1D707}_{s})\unicode[STIX]{x1D735}_{\!s}(\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }\boldsymbol{u}_{s})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D707}_{s}[2K\boldsymbol{u}_{s}+\boldsymbol{n}\times \unicode[STIX]{x1D735}_{\!s}(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!s}\times \boldsymbol{u}_{s})+2(\boldsymbol{n}\times \unicode[STIX]{x1D735}_{\!s}\boldsymbol{n}\times \boldsymbol{n})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!s}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})]\nonumber\\ \displaystyle & & \displaystyle +\,\boldsymbol{n}[2H\unicode[STIX]{x1D70E}+2H(\unicode[STIX]{x1D707}_{d}+\unicode[STIX]{x1D707}_{s})\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }\boldsymbol{u}_{s}-2\unicode[STIX]{x1D707}_{s}(\boldsymbol{n}\times \unicode[STIX]{x1D735}_{\!s}\boldsymbol{n}\times \boldsymbol{n})\boldsymbol{ : }\unicode[STIX]{x1D735}_{\!s}\boldsymbol{u}_{s}],\end{eqnarray}$$

where

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle 2H=-\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle 2K=-(\boldsymbol{n}\times \unicode[STIX]{x1D735}_{\!s}\boldsymbol{n}\times \boldsymbol{n})\boldsymbol{ : }\unicode[STIX]{x1D735}_{\!s}\boldsymbol{n} & \displaystyle\end{eqnarray}$$

are the mean and Gaussian curvatures of a surface respectively and $\unicode[STIX]{x1D70E}$ is the surface tension coefficient. Neglecting higher-order terms, the leading-order surface stress boundary condition in the context of small-amplitude capillary waves takes the reduced form

(2.7) $$\begin{eqnarray}[\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}]=\unicode[STIX]{x1D735}_{\!s}\tilde{\unicode[STIX]{x1D70E}}+2H\tilde{\unicode[STIX]{x1D70E}}\boldsymbol{n},\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D70E}}$ is the surface tension augmented with the leading-order surface viscosity contribution given by

(2.8) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D70E}}=\unicode[STIX]{x1D70E}+(\unicode[STIX]{x1D707}_{d}+\unicode[STIX]{x1D707}_{s})\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }\boldsymbol{u}_{s}.\end{eqnarray}$$

Using the equation of motion derived in § 3, the leading-order surface viscosity effect on the small-amplitude capillary wave can naturally be characterised (Lopez & Hirsa Reference Lopez and Hirsa1998) by the non-dimensional Boussinesq number

(2.9) $$\begin{eqnarray}B\equiv Bq_{d}+Bq_{s}=\left(\frac{\unicode[STIX]{x1D707}_{d}+\unicode[STIX]{x1D707}_{s}}{\unicode[STIX]{x1D707}}\right)k,\end{eqnarray}$$

where $Bq_{d}=\unicode[STIX]{x1D707}_{d}k/\unicode[STIX]{x1D707}$ and $Bq_{s}=\unicode[STIX]{x1D707}_{s}k/\unicode[STIX]{x1D707}$ are the Boussinesq dilatational and shear numbers respectively for dynamic viscosity $\unicode[STIX]{x1D707}$ and wavenumber $k$ . More explicitly, surface viscosity can be modelled to be proportional to the surfactant concentration (Ponce-Torres et al. Reference Ponce-Torres, Montanero, Herrada, Vega and Vega2017). In the case of detergents, experimental work by Brown, Thuman & McBain (Reference Brown, Thuman and McBain1953) suggests a bi-partisan action of the special solute pairs present in the detergent, where the primary constituent provides a large reservoir of surface-active material while the secondary constituent, lesser in amount, forms surface films of high viscosity. However, in the leading-order dynamics of the Boussinesq–Scriven formulation, surface viscosity is shown in § 3 to not depend explicitly on the surfactant concentration and enters only implicitly via the surface tension coefficient.

The other non-dimensional numbers of the system that arise naturally in the equation of motion are the viscosity ( $\unicode[STIX]{x1D716}$ ), surfactant diffusivity ( $\unicode[STIX]{x1D70D}$ ) and surfactant strength ( $\unicode[STIX]{x1D6FD}$ ) parameters, given by

(2.10) $$\begin{eqnarray}(\unicode[STIX]{x1D716},\unicode[STIX]{x1D70D},\unicode[STIX]{x1D6FD})=\frac{k}{\unicode[STIX]{x1D714}}\left(\unicode[STIX]{x1D708}k,D_{s}k,\frac{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6E4}_{0}}{\unicode[STIX]{x1D707}}\right),\end{eqnarray}$$

where $D_{s}$ denotes the coefficient of surface diffusivity, $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}$ is the kinematic viscosity for fluid density $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D6FC}=|\text{d}\unicode[STIX]{x1D70E}/\text{d}\unicode[STIX]{x1D6E4}|$ is the gradient of the surface tension coefficient, $\unicode[STIX]{x1D714}$ is the frequency of the capillary wave and $\unicode[STIX]{x1D6E4}_{0}$ is the initial surfactant concentration, which is assumed to be much less than the critical micelle concentration (cmc). In this system, these parameters act as the effective Reynolds, Schmidt and Marangoni numbers respectively.

3 Equations of motion

The dynamics of an incompressible fluid of viscosity $\unicode[STIX]{x1D707}$ and density $\unicode[STIX]{x1D70C}$ in a region of Reynolds number $\mathit{Re}=U\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}\ll 1$ satisfies the Stokes equation

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}(\boldsymbol{u}_{t}-\boldsymbol{F})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}=(u,v)$ is the two-dimensional fluid velocity field, $p$ is the pressure and $\boldsymbol{F}=-g\,\boldsymbol{j}$ is the external (gravitational) force, with $\boldsymbol{j}$ denoting the upward unit vector in the $y$ direction, and $g$ is the gravitational acceleration. The small-amplitude capillary wave is given at the free surface $F$ by the standing wave

(3.3) $$\begin{eqnarray}F(x,y,t)=y-a(t)\cos kx,\end{eqnarray}$$

where $a(t)$ is the nonlinear time-dependent wave amplitude which satisfies the small-amplitude conditions that $a\ll \unicode[STIX]{x1D706}\,=2\unicode[STIX]{x03C0}/k$ and $\text{d}a/\text{d}t\ll v_{c}=\unicode[STIX]{x1D714}/k$ , where $\unicode[STIX]{x1D706}$ is the wavelength and $v_{c}$ is the phase velocity.

For vanishing Gaussian curvature in a two-dimensional space, the leading-order tangential and normal stress components $T_{\Vert }$ and $T_{\bot }$ and the kinematic condition are given by

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle T_{\Vert }\equiv {\textstyle \frac{1}{2}}\unicode[STIX]{x1D707}(v_{x}+u_{y})=\unicode[STIX]{x1D735}_{\!s}\tilde{\unicode[STIX]{x1D70E}}, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle T_{\bot }\equiv -p+2\unicode[STIX]{x1D707}v_{y}=\tilde{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle F_{t}+vF_{y}=0 & \displaystyle\end{eqnarray}$$

respectively. The leading-order normal and tangent vectors are

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{n}\simeq (ak\sin kx,1), & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{t}\simeq (1,-ak\sin kx). & \displaystyle\end{eqnarray}$$

Similarly to the small-amplitude condition, we consider a small departure from the equilibrium surface tension and let the coefficient of surface tension $\unicode[STIX]{x1D70E}$ be defined via a linear equation of state,

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(x,t)=\unicode[STIX]{x1D70E}_{0}-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6E4}(x,t),\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{0}$ is the initial surface tension coefficient and $\unicode[STIX]{x1D6E4}(x,t)$ is the (dilute) concentration of a surfactant solution where adsorptive–desorptive processes are neglected.

Using the waveform $\unicode[STIX]{x1D6E4}(x,t)-\unicode[STIX]{x1D6E4}_{0}=\tilde{\unicode[STIX]{x1D6E4}}(t)\cos kx,$ the governing equation for surfactant concentration along a two-dimensional deforming surface (Stone Reference Stone1990) is given by

(3.10) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D6E4}}_{t}+k^{2}D_{s}\tilde{\unicode[STIX]{x1D6E4}}=k(a_{t}+\unicode[STIX]{x1D708}k\unicode[STIX]{x1D6FA}(0,t)\ast {\mathcal{F}}(t))\end{eqnarray}$$

to the leading order, $\unicode[STIX]{x1D714}_{z}(x,y,t)=\unicode[STIX]{x1D6FA}(y,t)\sin kx$ is the $z$ component of the vorticity, $\ast$ is the convolution operator and ${\mathcal{F}}(t)$ is the auxiliary function

(3.11) $$\begin{eqnarray}{\mathcal{F}}(t)=\frac{1}{\sqrt{\unicode[STIX]{x03C0}\unicode[STIX]{x1D708}k^{2}t}}\text{e}^{-\unicode[STIX]{x1D708}k^{2}t}-\text{erfc}\sqrt{\unicode[STIX]{x1D708}k^{2}t}.\end{eqnarray}$$

The fluid velocity and the pressure can be decomposed into the sum of an inviscid and a viscous part, i.e. $(\boldsymbol{u},p)=(\boldsymbol{u}^{\prime }+\boldsymbol{u}^{\prime \prime },\,p^{\prime }+p^{\prime \prime })$ . The inviscid part $(\boldsymbol{u}^{\prime },p^{\prime })$ satisfies the Euler problem

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}(\boldsymbol{u}_{t}^{\prime }-\boldsymbol{F})=-\unicode[STIX]{x1D735}p^{\prime }, & \displaystyle\end{eqnarray}$$

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle (F_{t}+v^{\prime }F_{y})|_{y=0}=0, & \displaystyle\end{eqnarray}$$

with well-known solutions (Lamb Reference Lamb1932)

(3.14) $$\begin{eqnarray}(\unicode[STIX]{x1D719},\,p^{\prime })=\left(\frac{1}{k}\frac{\text{d}a}{\text{d}t}\text{e}^{ky}\cos kx,\,-\unicode[STIX]{x1D70C}gy+\frac{\unicode[STIX]{x1D70C}}{k}\frac{\text{d}^{2}a}{\text{d}t^{2}}\text{e}^{ky}\cos kx\right),\end{eqnarray}$$

where $\boldsymbol{u}^{\prime }=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ . The viscous component $(\boldsymbol{u}^{\prime \prime },p^{\prime \prime })$ satisfies the Stokes problem

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\boldsymbol{u}_{t}^{\prime \prime }=\unicode[STIX]{x1D735}p^{\prime \prime }+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{\prime \prime }, & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle 0=v^{\prime }F_{y}|_{y=0}, & \displaystyle\end{eqnarray}$$

which is solved by introducing the streamfunction $\unicode[STIX]{x1D713}$ , defined by $\boldsymbol{u}^{\prime \prime }=(\unicode[STIX]{x1D713}_{y},-\unicode[STIX]{x1D713}_{x})$ . By taking the curl of (3.15), we obtain the bi-harmonic equation

(3.17) $$\begin{eqnarray}(\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6FB}^{2}-\unicode[STIX]{x1D6FB}^{4})\unicode[STIX]{x1D713}=0.\end{eqnarray}$$

By writing $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D6F9}(y,t)\sin kx$ , the Stokes problem yields the solution

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle 2k\unicode[STIX]{x1D6F9}=-\text{e}^{-ky}\int _{-\infty }^{y}\unicode[STIX]{x1D6FA}\text{e}^{ky^{\prime }}\,\text{d}y^{\prime }+\text{e}^{ky}\left(\int _{-\infty }^{0}\unicode[STIX]{x1D6FA}\text{e}^{ky^{\prime }}\text{d}y^{\prime }+\int _{0}^{y}\unicode[STIX]{x1D6FA}\text{e}^{-ky^{\prime }}\text{d}y^{\prime }\right), & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D6FA}(0,t)\ast \frac{y}{\sqrt{\unicode[STIX]{x03C0}\unicode[STIX]{x1D708}t^{3}}}\text{exp}\left(-\unicode[STIX]{x1D708}k^{2}t-\frac{y^{2}}{4\unicode[STIX]{x1D708}t}\right), & \displaystyle\end{eqnarray}$$

where we have the viscous pressure correction $p^{\prime \prime }(x,y,t)=\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FA}(0,t)\text{e}^{ky}\cos kx$ and the boundary vorticity

(3.20) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}(0,t)=2\left(-\frac{T_{\Vert }}{\unicode[STIX]{x1D707}}+v_{x}\right).\end{eqnarray}$$

Henceforth, using non-dimensional variables $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D714}t,$ $\unicode[STIX]{x1D716}=\unicode[STIX]{x1D708}k^{2}/\unicode[STIX]{x1D714}$ and $\tilde{\unicode[STIX]{x1D6FA}}=\unicode[STIX]{x1D6FA}(0,t)/\unicode[STIX]{x1D714}$ , the boundary vorticity becomes the integral equation

(3.21) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D6FA}}(0,\unicode[STIX]{x1D70F})=\text{f}(\unicode[STIX]{x1D70F})+2\unicode[STIX]{x1D716}\,B\,\tilde{\unicode[STIX]{x1D6FA}}(0,\unicode[STIX]{x1D70F}^{\prime })\ast {\mathcal{F}}(\unicode[STIX]{x1D70F}),\end{eqnarray}$$

where we have

(3.22) $$\begin{eqnarray}\text{f}(\unicode[STIX]{x1D70F})=-2[\unicode[STIX]{x1D6FD}\tilde{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D70F})+\unicode[STIX]{x1D6FF}(1+B){\dot{A}}].\end{eqnarray}$$

Here, $A=a/a_{0}$ is the dimensionless amplitude, $\unicode[STIX]{x1D6FF}=a_{0}k$ and $\dot{\,}=\text{d}/\text{d}\unicode[STIX]{x1D70F}$ denotes the non-dimensional temporal derivative. By substituting the pressure and the velocity into (3.5) and (3.10), we have the simultaneous equation

(3.23) $$\begin{eqnarray}\displaystyle & \displaystyle \ddot{A}+2\unicode[STIX]{x1D716}{\dot{A}}+A=\unicode[STIX]{x1D716}\tilde{\unicode[STIX]{x1D6FA}}(0,\unicode[STIX]{x1D70F})-2\unicode[STIX]{x1D716}^{2}\tilde{\unicode[STIX]{x1D6FA}}(0,\unicode[STIX]{x1D70F})\ast {\mathcal{F}}(\unicode[STIX]{x1D70F}), & \displaystyle\end{eqnarray}$$

(3.24) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\tilde{\unicode[STIX]{x1D6E4}}}+\unicode[STIX]{x1D70D}\tilde{\unicode[STIX]{x1D6E4}}=\unicode[STIX]{x1D6FF}{\dot{A}}+\unicode[STIX]{x1D716}\tilde{\unicode[STIX]{x1D6FA}}(0,\unicode[STIX]{x1D70F})\ast {\mathcal{F}}(\unicode[STIX]{x1D70F}). & \displaystyle\end{eqnarray}$$

Equations (3.23) and (3.24) provide us with a dynamic equation system for the amplitude and the surfactant concentration, the solution of which we outline in the next section.

4 Solution of the simultaneous integro-differential equation

Let $F(s)={\mathcal{L}}[A](s)$ , $G(s)={\mathcal{L}}[\tilde{\unicode[STIX]{x1D6E4}}](s)$ and $\hat{\unicode[STIX]{x1D6F1}}(s)=sF(s)-A_{0}$ be the Laplace transforms of $A(\unicode[STIX]{x1D70F}),\,\tilde{\unicode[STIX]{x1D6E4}}(\unicode[STIX]{x1D70F})$ and ${\dot{A}}(\unicode[STIX]{x1D70F})$ , and define the polynomial expressions $\unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D716}}^{(i)}\equiv \unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D716}}^{(i)}(s+\unicode[STIX]{x1D716})$ for $1\leqslant i\leqslant 6$ as

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D716}}^{(1)}=2[(s+\unicode[STIX]{x1D716})^{1/2}-\unicode[STIX]{x1D716}^{1/2}], & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D716}}^{(2)}=s\unicode[STIX]{x1D716}^{1/2}+B\unicode[STIX]{x1D716}\unicode[STIX]{x1D6E9}^{(1)}, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D716}}^{(3)}=(s+\unicode[STIX]{x1D70D})\unicode[STIX]{x1D6E9}^{(2)}+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D716}\unicode[STIX]{x1D6E9}^{(1)}, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D716}}^{(4)}=(s^{2}+2\unicode[STIX]{x1D716}s+1+2B^{\prime }\unicode[STIX]{x1D716}s)\unicode[STIX]{x1D6E9}^{(2)}-2\unicode[STIX]{x1D716}^{2}sB^{\prime 2}\unicode[STIX]{x1D6E9}^{(1)}, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D716}}^{(5)}=2\unicode[STIX]{x1D716}s\unicode[STIX]{x1D6FD}(B^{\prime }\unicode[STIX]{x1D716}\unicode[STIX]{x1D6E9}^{(1)}-\unicode[STIX]{x1D6E9}^{(2)}), & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D716}}^{(6)}=\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6E9}^{(2)}-B^{\prime }\unicode[STIX]{x1D716}\unicode[STIX]{x1D6E9}^{(1)}, & \displaystyle\end{eqnarray}$$

where $B^{\prime }=1+B$ . The rational function $\hat{\unicode[STIX]{x1D6F1}}_{\unicode[STIX]{x1D716}}=\hat{\unicode[STIX]{x1D6F1}}_{\unicode[STIX]{x1D716}}(s+\unicode[STIX]{x1D716})=\text{P}(s^{1/2})/\text{Q}(s^{1/2})$ is decomposed into its partial fraction

(4.7) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D6F1}}_{\unicode[STIX]{x1D716}}(s+\unicode[STIX]{x1D716}) & \equiv & \displaystyle \mathop{\sum }_{i=1}^{10}\frac{c_{i}}{s^{1/2}+z_{i}}\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{(U_{0}s-A_{0})\unicode[STIX]{x1D6E9}^{(2)}\unicode[STIX]{x1D6E9}^{(3)}+\tilde{\unicode[STIX]{x1D6E4}}_{0}\unicode[STIX]{x1D6E9}^{(2)}\unicode[STIX]{x1D6E9}^{(5)}}{\unicode[STIX]{x1D6E9}^{(4)}\unicode[STIX]{x1D6E9}^{(3)}-\tilde{\unicode[STIX]{x1D6E4}}_{0}\unicode[STIX]{x1D6E9}^{(5)}\unicode[STIX]{x1D6E9}^{(6)}},\end{eqnarray}$$

where $-z_{i}$ are the roots of the polynomial $\text{Q}(s^{1/2})$ . In the absence of the Marangoni effect, the surface viscosity case is given by

(4.9) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D6F1}}_{\unicode[STIX]{x1D716}}(s+\unicode[STIX]{x1D716})=\frac{(U_{0}s-A_{0})\unicode[STIX]{x1D6E9}^{(2)}}{(s^{2}+2\unicode[STIX]{x1D716}s+1+2B^{\prime }\unicode[STIX]{x1D716}s)\unicode[STIX]{x1D6E9}^{(2)}-2\unicode[STIX]{x1D716}^{2}sB^{\prime 2}\unicode[STIX]{x1D6E9}^{(1)}}.\end{eqnarray}$$

By comparison with Lagrange polynomial interpolation, we have

(4.10) $$\begin{eqnarray}\displaystyle \frac{\text{P}(s^{1/2})}{\text{Q}(s^{1/2})} & \equiv & \displaystyle \text{P}(s^{1/2})\mathop{\prod }_{i=1}^{k}\frac{1}{s^{1/2}+z_{i}}\end{eqnarray}$$

(4.11) $$\begin{eqnarray}\displaystyle & = & \displaystyle \mathop{\sum }_{i=1}^{k}\frac{\text{P}(-z_{i})}{\unicode[STIX]{x1D70E}_{i}^{(k)}(-z_{i})}\frac{1}{s^{1/2}+z_{i}},\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{i}^{(n)}$ is the $n$ th-order cyclic polynomial given by

(4.12) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{j}^{(n)}=\mathop{\prod }_{i=1}^{n-1}(z_{j+i\,\mathit{mod}(n)}-z_{j}).\end{eqnarray}$$

It follows that by comparing (4.7) and (4.11), the coefficients $c_{i}$ are $c_{i}=\text{P}(-z_{i})/\unicode[STIX]{x1D70E}_{i}^{(10)}(-z_{i})$ . Let

(4.13) $$\begin{eqnarray}\text{Z}(n,j)=\mathop{\sum }_{i=1}^{n}\frac{\text{P}(-z_{i})}{\unicode[STIX]{x1D70E}_{i}^{(n)}}(-z_{i})^{j};\end{eqnarray}$$

it follows (see appendix A) that the condition

(4.14) $$\begin{eqnarray}\deg \text{Q}-\deg \text{P}=2\end{eqnarray}$$

implies $\text{Z}(n,0)=0$ , where $\deg X$ is the degree of the polynomial $X$ . By taking the inverse Laplace transform of (4.8), we obtain

(4.15) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D70F})=\frac{\text{Z}(10,0)}{\sqrt{\unicode[STIX]{x03C0}t}}-\mathop{\sum }_{i=1}^{10}\frac{\text{P}(-z_{i})}{\unicode[STIX]{x1D70E}_{i}^{(10)}}z_{i}\text{e}^{z_{i}^{2}\unicode[STIX]{x1D70F}}\text{erfc}(z_{i}\unicode[STIX]{x1D70F}^{1/2}).\end{eqnarray}$$

Finally, the non-dimensional amplitude is given by

(4.16) $$\begin{eqnarray}A(\unicode[STIX]{x1D70F})=1+\mathop{\sum }_{i=1}^{10}\frac{z_{i}}{\unicode[STIX]{x1D70E}_{i}^{(10)}}\text{P}(-z_{i})\unicode[STIX]{x1D711}(z_{i},\unicode[STIX]{x1D70F};\unicode[STIX]{x1D716}),\end{eqnarray}$$

where $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D711}(z_{i},\unicode[STIX]{x1D70F};\unicode[STIX]{x1D716})$ satisfies

(4.17) $$\begin{eqnarray}\unicode[STIX]{x1D711}(z_{i},\unicode[STIX]{x1D70F};\unicode[STIX]{x1D716})=\frac{1}{z_{i}^{2}-\unicode[STIX]{x1D716}}\left(\text{e}^{(z_{i}^{2}-\unicode[STIX]{x1D716})\unicode[STIX]{x1D70F}}\text{erfc}(z_{i}\unicode[STIX]{x1D70F}^{1/2})+\frac{z_{i}}{\unicode[STIX]{x1D716}^{1/2}}\text{erf}[(\unicode[STIX]{x1D716}\unicode[STIX]{x1D70F})^{1/2}]-1\right).\end{eqnarray}$$

5 Surface viscosity effects on the wave amplitude

As shown in § 4, the leading-order surface viscosity effects on the dynamics of small-amplitude capillary waves are characterised by the parameter $B$ , and, similarly, the Marangoni effect can be reduced to the non-dimensional variables $\unicode[STIX]{x1D70D}$ and $\unicode[STIX]{x1D6FD}$ given in § 2. In what follows, we use water under room temperature and pressure (rtp), i.e. at $25\,^{\circ }\text{C}$ , with density $\unicode[STIX]{x1D70C}=10^{3}~\text{kg}~\text{m}^{-3}$ , surface tension $\unicode[STIX]{x1D70E}=7.2\times 10^{-2}~\text{N}~\text{m}^{-1}$ and viscosity $\unicode[STIX]{x1D707}=8.9\times 10^{-4}~\text{Pa}~\text{s}$ , as a test system with surfactant Schmidt number $\mathit{Sc}=\unicode[STIX]{x1D708}/D_{s}=10^{4}$ . We define $\unicode[STIX]{x1D706}_{c}^{(0)}$ to be the wavelength for which the capillary wave undergoes critical damping, henceforth known as the critical wavelength. The superscript denotes the clean case, which we understand as a system without the addition of surface-active material, i.e. $\unicode[STIX]{x1D6FD}=B=0$ . We will look at this critical wavelength in more detail in § 7, but here we note a harmonic oscillator approximation of $\unicode[STIX]{x1D706}_{c}^{(0)}$ (Denner Reference Denner2016) whereby

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{c}^{(0)}=2^{1/3}\unicode[STIX]{x03C0}l_{vc}/\unicode[STIX]{x1D6E9}\end{eqnarray}$$

for $\unicode[STIX]{x1D6E9}=1.0625$ and the viscocapillary length scale

(5.2) $$\begin{eqnarray}l_{vc}=\frac{\unicode[STIX]{x1D707}^{2}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}.\end{eqnarray}$$

In figure 1, we compare the effect of surface viscosity with that of the equivalent bulk viscosity $\unicode[STIX]{x1D708}$ , and the Marangoni effect with that of simple reductions in the surface tension coefficient $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{0}-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6E4}_{0}$ . Here, we use the phrase equivalent to denote the quantities of bulk viscosity and simple reductions in surface tension that result in an identical effect on the overall wave amplitude due to surface viscosity and the Marangoni effect respectively under the limit of either a large or a small wavelength. For wavelengths near critical damping, in figure 1(a), an increase in $B$ exhibits a relatively large difference from an equivalent increase in $\unicode[STIX]{x1D708}$ , as compared with the case for $\unicode[STIX]{x1D706}=950\unicode[STIX]{x1D706}_{c}^{(0)}$ in figure 1(d). This difference decreases as we increase the wavelength into less damped regions, as shown in figure 1(b,c). In contrast, the picture is reversed for the Marangoni effect, where equivalent changes in the surface tension coefficient are almost identical to increases in surfactant concentration (through $\unicode[STIX]{x1D6FD}$ ) near critical damping in figure 1(e), and their difference increases with larger wavelength from figure 1(f) to 1(h), where the dynamic Marangoni effect vastly overshadows the static changes in the surface tension coefficient. Henceforth, we can approximate the Marangoni effect on the wave amplitude near the critical wavelength with the equivalent reduction in $\unicode[STIX]{x1D70E}$ .

Figure 1. The wave amplitude as a function of the dimensionless time $\unicode[STIX]{x1D70F}$ for water at approximate room temperature and pressure for $\mathit{Sc}=10^{4}$ , comparing the influence of the surface viscosity and Marangoni effects for various values of $B$ (with $\unicode[STIX]{x1D6FD}=0$ ) and $\unicode[STIX]{x1D6FD}$ (with $B=0$ ). Here, $\unicode[STIX]{x1D708}_{0}=8.9\times 10^{-7}~\text{m}^{2}~\text{s}^{-1}$ and $\unicode[STIX]{x1D70E}_{0}=0.072~\text{N}~\text{m}^{-1}$ denote the baseline kinematic viscosity and surface tension respectively. The modified surface tension in (e,f) represents a clean system (i.e. $\unicode[STIX]{x1D6FD}=B=0$ ) with $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{0}-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6E4}_{0}$ .

Figure 2. The wave amplitude as a function of the dimensionless time $\unicode[STIX]{x1D70F}$ for water at approximate room temperature and pressure for $\mathit{Sc}=10^{4}$ , showing the influence of the surface viscosity and Marangoni effects for various values of $B$ (with $\unicode[STIX]{x1D6FD}=0$ ) and $\unicode[STIX]{x1D6FD}$ (with $B=0$ ) for different wavelengths $\unicode[STIX]{x1D706}$ .

In figure 2, we see in more detail the mechanisms of surface viscosity and the Marangoni effect at different wave amplitudes for different wavelengths. Both effects admit relatively self-similar solutions, and on increasing either the surface viscosity or the surfactant concentration, the damping is increased, as expected. One of the explicit differences is that the Marangoni effect reduces the surface tension and, thus, lowers the frequency $\unicode[STIX]{x1D714}$ , while the surface viscosity leaves $\unicode[STIX]{x1D714}$ unchanged. This is due to the surface viscosity being dependent on the surfactant concentration only to the linear order, and it does not feature in the leading-order nonlinear amplitude equation (3.23) explicitly, as noted previously. The consequence on the amplitude of this $\unicode[STIX]{x1D714}$ lowering is a horizontal drift of the waveform for systems with the Marangoni effect, as is evident from figure 2(d,f), as opposed to relatively centred waveforms in the cases of surface viscosity in figure 2(c,e).

We also observe that the surface viscosity effect weakens for large wavelengths $\unicode[STIX]{x1D706}\gg \unicode[STIX]{x1D706}_{c}$ but is very potent for small wavelengths $\unicode[STIX]{x1D706}\sim \unicode[STIX]{x1D706}_{c}$ . A particular consequence of this potency is its ability to alter the onset behaviour in interfacial phenomena, a number of which occur near the region of the critical wavelength. The stochastic nature of many of the interfacial instabilities (Aarts et al. Reference Aarts, Schmidt and Lekkerkerker2004) is often kickstarted by the small-amplitude local disturbances on the interface. Hence, a small change in surface material, and thus the wave damping, can have a significant effect in starting or delaying the initialisation process of more complex phenomena. Furthermore, it would be of interest to obtain the modification to the critical wavelength upon the addition of a small amount surface-active material. However, we need to consider the definition of the critical wavelength for a higher-order system as it is not as readily defined as in a second-order system.

6 Capillary wave dispersion and the critical wavelength

To quantify the changes to the critical wavelength due to the presence of the Marangoni and surface viscosity effects, we must first obtain the critical wavelength in the clean case. Following Lamb (Reference Lamb1932), the general dispersion relation for an interface with both the Marangoni effect and surface viscosity can be found (derivation in appendix B) to take the form

(6.1) $$\begin{eqnarray}W_{0}(Z;\unicode[STIX]{x1D716})\left(1+Z+B+\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D716}^{2}(Z^{2}-1)}\right)+\left(B+\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D716}^{2}(Z^{2}-1)}\right)(Z-1)^{3}=0,\end{eqnarray}$$

where

(6.2) $$\begin{eqnarray}\frac{\text{i}\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D716}}=Z^{2}-1\end{eqnarray}$$

and

(6.3) $$\begin{eqnarray}W_{0}(Z;\unicode[STIX]{x1D716})=Z^{4}+2Z^{2}-4Z+1+\frac{1}{\unicode[STIX]{x1D716}^{2}}.\end{eqnarray}$$

On specialising to the clean case, (6.1) reduces to

(6.4) $$\begin{eqnarray}(\text{i}\unicode[STIX]{x1D714}^{\prime }+2\unicode[STIX]{x1D716})^{2}-4\unicode[STIX]{x1D716}^{2}\left(1+\frac{\text{i}\unicode[STIX]{x1D714}^{\prime }}{\unicode[STIX]{x1D716}}\right)^{1/2}+1=0,\end{eqnarray}$$

as derived from linearised hydrodynamics (Lamb Reference Lamb1932; Levich Reference Levich1962). The dispersion relation admits a solution of the form

(6.5) $$\begin{eqnarray}\unicode[STIX]{x1D714}^{\prime }=2\text{i}\unicode[STIX]{x1D716}-\frac{1}{2}h-\left(1-\frac{1}{4}h^{2}-\frac{8\text{i}\unicode[STIX]{x1D716}^{3}}{h}\right)^{1/2},\end{eqnarray}$$

where $h^{2}=(1/3)-J_{+}^{1/3}-J_{-}^{1/3}$ for

(6.6) $$\begin{eqnarray}J_{\pm }=\frac{1}{27}-\frac{2}{3}\unicode[STIX]{x1D716}^{4}+2\unicode[STIX]{x1D716}^{6}\pm \frac{2}{3\sqrt{3}}\unicode[STIX]{x1D716}^{3}\sqrt{\text{f}(\unicode[STIX]{x1D716})},\end{eqnarray}$$

where the polynomial $\text{f}(\unicode[STIX]{x1D716})$ is given by

(6.7) $$\begin{eqnarray}\text{f}(\unicode[STIX]{x1D716})=11\unicode[STIX]{x1D716}^{6}-18\unicode[STIX]{x1D716}^{4}-\unicode[STIX]{x1D716}^{2}-1.\end{eqnarray}$$

For $\unicode[STIX]{x1D716}\ll \unicode[STIX]{x1D716}^{\star },$ the wave frequency can be written as

(6.8) $$\begin{eqnarray}\unicode[STIX]{x1D714}^{\prime }\sim \frac{1}{2}h^{2}-h+2\text{i}\unicode[STIX]{x1D716}\left(1+\frac{\unicode[STIX]{x1D716}^{2}}{h}\right)\end{eqnarray}$$

and the damping coefficient can be extracted as $\text{Im}(\unicode[STIX]{x1D714}^{\prime })\sim 2\unicode[STIX]{x1D716},$ where $\unicode[STIX]{x1D716}^{\star }\in \mathbb{R}^{+}$ is the transition value defined by the (largest positive) root of the polynomial $\text{f},$ i.e.

(6.9) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{\star }=\sup _{\mathbb{ R}^{+}}\{\unicode[STIX]{x1D716}:\text{f}(\unicode[STIX]{x1D716})=0\}.\end{eqnarray}$$

By solving (6.7) exactly, we obtain

(6.10) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{\star }=\left({\textstyle \frac{6}{11}}+{\textstyle \frac{1}{33}}\left({\textstyle \frac{3}{2}}(5571-341\sqrt{93})\right)^{1/3}+{\textstyle \frac{1}{33}}\left({\textstyle \frac{3}{2}}(5571+341\sqrt{93})\right)^{1/3}\right)^{1/2},\end{eqnarray}$$

with the numerical value $\unicode[STIX]{x1D716}^{\star }\simeq 1.3115.$ Reintroducing dimensional variables, we define $P,Q$ and the variable $K$ by

(6.11) $$\begin{eqnarray}\displaystyle & \displaystyle K=k^{\star }-\frac{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D716}^{\star 2}}{3\unicode[STIX]{x1D708}^{2}\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$

(6.12) $$\begin{eqnarray}\displaystyle & \displaystyle P=-\frac{\unicode[STIX]{x1D70E}^{2}\unicode[STIX]{x1D716}^{\star 4}}{3\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x1D708}^{4}}, & \displaystyle\end{eqnarray}$$

(6.13) $$\begin{eqnarray}\displaystyle & \displaystyle Q=-\left(\frac{\text{g}\unicode[STIX]{x1D716}^{\star 2}}{\unicode[STIX]{x1D708}^{2}}+\frac{2\unicode[STIX]{x1D70E}^{3}\unicode[STIX]{x1D70E}^{\star 6}}{27\unicode[STIX]{x1D70C}^{3}\unicode[STIX]{x1D708}^{6}}\right). & \displaystyle\end{eqnarray}$$

It follows from the definition of $\unicode[STIX]{x1D716}$ that the critical wavenumber $k^{\star }$ satisfies the cubic equation

(6.14) $$\begin{eqnarray}K^{3}+PK+Q=0.\end{eqnarray}$$

We require the real solution given by

(6.15) $$\begin{eqnarray}K=\left\{-{\textstyle \frac{1}{2}}Q+\unicode[STIX]{x1D6E5}^{1/2}\right\}^{1/3}+\left\{-{\textstyle \frac{1}{2}}Q-\unicode[STIX]{x1D6E5}^{1/2}\right\}^{1/3},\end{eqnarray}$$

where

(6.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}^{1/2} & \equiv & \displaystyle \left(\frac{Q^{2}}{4}+\frac{P^{3}}{27}\right)^{1/2}\end{eqnarray}$$

(6.17) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{\text{g}\unicode[STIX]{x1D716}^{\star 2}}{2\unicode[STIX]{x1D708}^{2}}\left(1+\frac{4\unicode[STIX]{x1D70E}^{3}\unicode[STIX]{x1D716}^{\star 4}}{27\unicode[STIX]{x1D70C}^{3}\unicode[STIX]{x1D708}^{4}\text{g}}\right)^{1/2}.\end{eqnarray}$$

By inspecting (6.15), the critical wavenumber under the limit $k\ll (\unicode[STIX]{x1D70C}\text{g}/\unicode[STIX]{x1D70E})^{1/2}$ is $k^{\star }\sim (\text{g}\unicode[STIX]{x1D716}^{\star 2}/\unicode[STIX]{x1D708}^{2})^{1/3},$ corresponding to the gravity-dominated regime with $\unicode[STIX]{x1D714}_{0}\sim (\text{g}k)^{1/2}$ . For $k\geqslant (\unicode[STIX]{x1D70C}\text{g}/\unicode[STIX]{x1D70E})^{1/2}$ , the critical wavenumber reduces to

(6.18) $$\begin{eqnarray}k^{\star }\sim \unicode[STIX]{x1D716}^{\star 2}\frac{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}^{2}}=\frac{\unicode[STIX]{x1D716}^{\star 2}}{l_{vc}},\end{eqnarray}$$

corresponding to the capillary-dominated regime with $\unicode[STIX]{x1D714}_{0}\sim (\unicode[STIX]{x1D70E}k^{3}/\unicode[STIX]{x1D70C})^{1/2}$ .

For $\unicode[STIX]{x1D716}\gg \unicode[STIX]{x1D716}^{\star }$ , we note that $\mathit{Re}(\unicode[STIX]{x1D714}^{\prime })=0$ and the system is in an overall overdamped regime. The damping ratio is given by expanding $\unicode[STIX]{x1D714}_{-}^{\prime }$ (since $\unicode[STIX]{x1D714}_{+}^{\prime }\gg \unicode[STIX]{x1D714}_{-}^{\prime }$ and would thus rapidly damp the motion) in (6.5) in ascending powers of $\unicode[STIX]{x1D716}^{-1}$ , i.e.

(6.19) $$\begin{eqnarray}\text{Im}(\unicode[STIX]{x1D714}^{\prime })\sim \frac{1}{2\unicode[STIX]{x1D716}}+O\left(\frac{1}{\unicode[STIX]{x1D716}^{2}}\right).\end{eqnarray}$$

This agrees with the asymptotic approximations by Levich (Reference Levich1962), which suggests that an increase of viscosity would decrease the damping for $\unicode[STIX]{x1D716}\gg \unicode[STIX]{x1D716}^{\star }$ .

Using (6.18), the analytical critical wavelength (the value of $\unicode[STIX]{x1D706}$ with associated damping ratio $\unicode[STIX]{x1D701}=1$ ) of the capillary wave is given by

(6.20) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{c}^{\star }=\frac{2\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D716}^{\star 2}}l_{vc}.\end{eqnarray}$$

For water under rtp, we have $\unicode[STIX]{x1D706}_{c}^{\star }\doteq 40.1894~\text{nm}$ . In comparison, a harmonic oscillator approximation (Denner Reference Denner2016) in (5.1) gives the result $\unicode[STIX]{x1D706}_{c}^{(0)}\doteq 40.9838~\text{nm}$ . We consider the relative error between the harmonic oscillator and the analytical critical wavelengths,

(6.21) $$\begin{eqnarray}\left|\frac{2^{1/3}}{\unicode[STIX]{x1D6E9}}-\frac{2}{\unicode[STIX]{x1D716}^{\star 2}}\right|\frac{\unicode[STIX]{x1D716}^{\star 2}}{2}\simeq 0.01977.\end{eqnarray}$$

We observe that the system is largely second-order in the neighbourhood of the critical wavelength, as the harmonic oscillator value of $\unicode[STIX]{x1D706}_{c}$ is within 2 % of the analytical value from the wave dispersion.

7 Damping ratio for a generalised system

For systems of a higher order, the damping ratio $\unicode[STIX]{x1D701}$ is not naturally defined and we usually inspect the root-locus diagram in order to decompose the system into a sum of first- and second-ordered systems to provide an estimate calculation. Here, we consider a numerical method to obtain an equivalent damping ratio, whereby $\unicode[STIX]{x1D701}\geqslant 1$ when the area of the amplitude below the settling value vanishes for all time almost everywhere. The critical wavelength is then the supremum of the set of wavelengths such that the above property holds. We express this as

(7.1) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{c}=\sup _{\unicode[STIX]{x1D706}\in \mathbb{R}^{+}}\left\{\unicode[STIX]{x1D706}:\lim _{T\rightarrow \infty }\int _{0}^{T}\text{A}(\unicode[STIX]{x1D706},\unicode[STIX]{x1D70F})[1-H(\text{A}(\unicode[STIX]{x1D706},\unicode[STIX]{x1D70F}))]\,\text{d}\unicode[STIX]{x1D70F}\rightarrow 0\,\text{a.e.}\right\},\end{eqnarray}$$

where $H(x)$ is the Heaviside step function.

For underdamped waves, even for second-order systems, logarithmic decrement or fractional overshoot methods tend to break down or become less accurate near regions of critical damping. Hence, to determine the damping ratios in the neighbourhood of $\unicode[STIX]{x1D701}\approx 1$ , we adapt the area method in (7.1). We consider that the area under the $t$ -axis is given by the function

(7.2) $$\begin{eqnarray}\unicode[STIX]{x1D6EF}(\unicode[STIX]{x1D701})=\int _{0}^{\infty }\,\text{d}t\,\{\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D701},t)\,[1-H(\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D701},t))]\},\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D701},t)$ satisfies the normalised harmonic oscillator equation

(7.3) $$\begin{eqnarray}\frac{\text{d}^{2}\unicode[STIX]{x1D6EC}}{\text{d}t^{2}}+2\unicode[STIX]{x1D701}\frac{\text{d}\unicode[STIX]{x1D6EC}}{\text{d}t}+\unicode[STIX]{x1D6EC}=0\quad \text{subject to }\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D701},0)=1,\quad \frac{\text{d}\unicode[STIX]{x1D6EC}}{\text{d}t}(\unicode[STIX]{x1D701},0)=0.\end{eqnarray}$$

The generalised (numerical) damping ratio for $\unicode[STIX]{x1D701}\leqslant 1$ can then be obtained by the inverse operation

(7.4) $$\begin{eqnarray}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D6EF}^{-1}(X),\end{eqnarray}$$

where $X$ is the area under the $t$ -axis of a generalised system. This numerical method agrees well with logarithmic decrement and fractional overshoot schemes in the relevant underdamped regimes.

In cases where the higher-order system can readily be approximated locally by a second-order system, we note that its dominant poles have a larger residue and time constant $t_{c}=1/(\unicode[STIX]{x1D701}_{n}\unicode[STIX]{x1D714}_{n})$ relative to other poles, where $\unicode[STIX]{x1D701}_{n}$ and $\unicode[STIX]{x1D714}_{n}$ are the damping ratio and frequency associated with each pole. In cases that are not clear cut, i.e. where all the poles are closer together with $t_{c}$ and residues of a similar magnitude, the numerical definition of the damping ratio above only provides an estimate of the true damping ratio and we need to examine the poles in more detail.

To encode such information into a convenient form, we construct the minimal pole matrix $\{\unicode[STIX]{x1D701}(\unicode[STIX]{x1D706})\}$ . To decompose the system, we first consider the minimal realisation of the transfer function. For a general system, let $\unicode[STIX]{x1D6E9}_{1}=\{q_{i}\in \mathbb{C}:~\text{Q}(-q_{i})=0\}$ and $\unicode[STIX]{x1D6E9}_{2}=\{p_{j}\in \mathbb{C}:\text{P}(-p_{j})=0\}$ be the sets of poles and zeros of the transfer function $\text{tf}(s)\equiv \text{P}(s)/\text{Q}(s)$ , which can be written in the form

(7.5) $$\begin{eqnarray}\text{tf}(s)=\mathop{\prod }_{q_{i}\in \unicode[STIX]{x1D6E9}_{1},p_{j}\in \unicode[STIX]{x1D6E9}_{2}}\left(\frac{s+p_{j}}{s+q_{i}}\right).\end{eqnarray}$$

By applying the pole-zero cancellation procedure, we obtain the minimal realisation of the transfer function, henceforth known as the minimal transfer function $\text{mtf}(s)$ , defined by

(7.6) $$\begin{eqnarray}\displaystyle \text{mtf}(s,\unicode[STIX]{x1D717}) & = & \displaystyle \frac{\text{P}_{m}(s)}{\text{Q}_{m}(s)}\end{eqnarray}$$

(7.7) $$\begin{eqnarray}\displaystyle & = & \displaystyle \mathop{\sum }_{q_{i,m}\in \unicode[STIX]{x1D6E9}_{1,m}(\unicode[STIX]{x1D717})}\frac{\text{res}(-q_{i,m})}{s+q_{i,m}},\end{eqnarray}$$

where the polynomials $\text{P}_{m}(s)$ and $\text{Q}_{m}(s)$ satisfy

(7.8) $$\begin{eqnarray}\deg (\text{Q}_{m}+\text{P}_{m})\leqslant \deg (\text{P}+\text{Q})\end{eqnarray}$$

and $\unicode[STIX]{x1D6E9}_{1,m}(\unicode[STIX]{x1D717})\subseteq \unicode[STIX]{x1D6E9}_{1}$ is the set of poles of $\text{Q}$ with significant residues (tolerance of order  $\unicode[STIX]{x1D717}$ ),

(7.9) $$\begin{eqnarray}\unicode[STIX]{x1D6E9}_{1,m}(\unicode[STIX]{x1D717})=\left\{q_{i,m}\in \mathbb{C}:\text{Q}(-q_{i,m})=0,~\frac{|\text{res}(q_{j})|}{\max _{q_{i}\in \unicode[STIX]{x1D6E9}_{1}}|\text{res}(q_{i})|}=O(\unicode[STIX]{x1D717})\right\}.\end{eqnarray}$$

Returning to the construction of the minimal pole matrix, we let the first column of the matrix illustrate the order of the poles of the minimal transfer function in dot form; the second column considers the relative magnitudes of their time constant $t_{c}$ ; the third column compares their relative residues; the fourth column gives the damping ratios $\unicode[STIX]{x1D701}_{n}$ associated with each pole. Furthermore, to the right of the line separator, we provide an estimated equivalent second-order damping ratio of the entire system using the area numerical method described previously. We say that a system is second-order dominant if one set of complex conjugate poles dominates the other poles (i.e. having the largest $t_{c}$ and residue). In diagrammatic form for $\unicode[STIX]{x1D717}=1$ , we have

(7.10)

For example, for $\unicode[STIX]{x1D706}=6.22\unicode[STIX]{x1D706}_{c}^{\star }$ at rtp, the clean case for water exhibits the following minimal pole matrix:

(7.11) $$\begin{eqnarray}\displaystyle \left\{\left.\begin{array}{@{}cccc@{}}\bullet \bullet & 1 & 1 & 0.75\\ \bullet \bullet & 1 & 2.30 & 0.45\end{array}\right|0.45\right\}, & & \displaystyle\end{eqnarray}$$

from which we observe that the system is second-order dominant since the set of complex conjugate poles with associated damping ratio 0.45 dominates (in the sense of residue) the other, and, thus, we can deduce that the true damping ratio of the system is close to the approximate second-order value.

We take this analysis to the region near critical damping, i.e. for $\unicode[STIX]{x1D701}\rightarrow 1^{+}$ and $\unicode[STIX]{x1D701}\rightarrow 1^{-}$ . In the clean case, we take the analytical result $\unicode[STIX]{x1D706}_{c}^{\star }=2\unicode[STIX]{x03C0}l_{vc}/\unicode[STIX]{x1D716}^{\star 2}$ to be the definition of the critical (damping) wavelength. The relevant minimal pole matrices take the form

(7.12) $$\begin{eqnarray}\displaystyle & \displaystyle \{\unicode[STIX]{x1D701}(\unicode[STIX]{x1D706}\rightarrow \unicode[STIX]{x1D706}_{c}^{\star +})\}=\left\{\begin{array}{@{}cccc@{}}\bullet \bullet & 1 & 1 & 1\\ \bullet \bullet & 1^{-} & 1^{+} & {<}1\end{array}\right\}, & \displaystyle\end{eqnarray}$$

(7.13) $$\begin{eqnarray}\displaystyle & \displaystyle \{\unicode[STIX]{x1D701}(\unicode[STIX]{x1D706}\rightarrow \unicode[STIX]{x1D706}_{c}^{\star -})\}=\left\{\begin{array}{@{}cccc@{}}\bullet & 1^{-} & 1 & 1\\ \bullet & 1 & 1^{-} & 1\end{array}\right\}. & \displaystyle\end{eqnarray}$$

We can see that this transition from $\unicode[STIX]{x1D706}_{c}^{\star +}\rightarrow \unicode[STIX]{x1D706}_{c}^{\star -}$ for the clean case boasts a transformation of the complex conjugate into two separate first-order poles, or, in diagrammatic form, we have

(7.14) $$\begin{eqnarray}\begin{array}{@{}c@{}}\bullet \bullet \\ \bullet \bullet \end{array}\rightarrow \begin{array}{@{}c@{}}\bullet \\ \bullet \end{array},\end{eqnarray}$$

i.e. a $2^{2}$ – $1^{2}$ transition.

Extending to contaminated systems, we summarise the results of the minimal pole matrices of the system at the critical transition in table 1, where the notation $1^{a}2^{b}$ denotes a system with $a$ first-order and $b$ second-order poles with significant relative residue. We note that the effect of surface viscosity is to introduce an extra first-order pole (with unit damping ratio) to the system, while the Marangoni effect does not change the pole composition for the underdamped region before the critical damping transition. Moreover, we observe at this critical damping transition that the system enters the overdamped regime if its minimal pole matrix is of the $1^{2}$ type irrespective of its type in the underdamped regime prior to the transition. Henceforth, we shall define a generalised higher-order (capillary) wave to be in overall overdamped motion if its minimal pole matrix is of the $1^{2}$ type.

Table 1. Schematic of poles at the transition from $\unicode[STIX]{x1D706}_{c}^{\star +}\rightarrow \unicode[STIX]{x1D706}_{c}^{\star -}$ for water at room temperature and pressure.

Using the definition of the overdamped regime for a generalised higher-order system, we determine numerically the corrections to the critical wavelength for increasing surface viscosity (through $B$ ), surfactant concentration (through $\unicode[STIX]{x1D6FD}$ ) and bulk viscosity (through $\unicode[STIX]{x1D703}$ , where

(7.15) $$\begin{eqnarray}\unicode[STIX]{x1D708}=(1+\unicode[STIX]{x1D703})\unicode[STIX]{x1D708}_{0}\end{eqnarray}$$

in figure 3). We note that the curves denoting surface viscosity and the Marangoni effect intersect near $B\simeq 0.65$ and that while the Marangoni curve is roughly exponential, the surface viscosity curve is the sum of two exponential functions. Moreover, we observe that a small amount of surface viscosity present in the system has an amplified effect on the system and that a sevenfold increase in critical wavelength for $B=1$ results in very different dynamics and mechanisms as the sub-100 nm brings forward the possibility of long-range molecular interactions as well as the hydrodynamics. Also of consideration is the proximity of the critical wavelength to the wavelengths of the visible spectrum of light, which allows thin-film behaviours to be captured by light scattering methods. Hence, the presence of surfactants could determine whether or not we would be within such a range as to allow interferometry techniques.

Figure 3. Independent corrections to the critical wavelength for increasing surface viscosity (via $B$ , with $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FD}=0$ ), bulk viscosity (via $\unicode[STIX]{x1D703}$ , where $\unicode[STIX]{x1D708}/\unicode[STIX]{x1D708}_{0}=1+\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D6FD}=B=0$ ) and Marangoni effect (through $\unicode[STIX]{x1D6FD}$ , with $\unicode[STIX]{x1D703}=B=0$ ).

Comparisons of the range with experimental and computational results on surface viscosity can be made using values reported previously in the literature. Experimentally, Kanner & Glass (Reference Kanner and Glass1969) summarised in table 2 the surface viscosities for both surfactant and polymeric films, where, for a dilute amount of sodium lauryl sulphate and polydimethylsiloxane in particular, the surface viscosity corresponds to a lower bound of $B=O(1)$ for $k=O(10^{6})$ . A similar correspondence can be found with the upper bound surface shear viscosity of $O(10^{-8}~\text{N}~\text{s}~\text{m}^{-1})$ found in Zell et al. (Reference Zell, Nowbahar, Mansard, Leal, Deshmukh, Mecca, Tucker and Squires2014) for soluble surfactants. We note that this measurement does not include surface dilatational viscosity and so corresponds to the case (Lucassen & Hansen Reference Lucassen and Hansen1966) where diffusional transport between the surface and the bulk is neglected, and assumes that the bulk viscosity and density are constant right up to the interface. More recently, Gounley et al. (Reference Gounley, Boedec, Jaeger and Leonetti2016) characterised the influence of both shear and dilatational surface viscosity on droplets in shear flow in the range $B=O(10^{0})$ – $O(10^{1})$ for a range of capillary numbers.

Beyond the cmc value, the Marangoni effect should in principle have no overall contribution to the capillary wave; the dotted curve in figure 3 would end abruptly at the cmc value. The combination of surface viscosity together with the Marangoni effect is, however, not straightforward; as in previous experimental studies (Brown et al. Reference Brown, Thuman and McBain1953; Kanner & Glass Reference Kanner and Glass1969), surface viscosity also appears to alter the ability of the Marangoni effect to lower surface tension. It would therefore be fruitful in a future contribution to investigate this surfactant interference mechanism through a more systematic experimental and theoretical study. In particular, a numerical approach similar to that of Sinclair, Levy & Daniels (Reference Sinclair, Levy and Daniels2018) could include the use of a nonlinear equation of state for the surface tension coefficient $\unicode[STIX]{x1D70E}$ . The effect of the deviation from the linear equation of state on the amplitude of the capillary wave would aid the analysis near the cmc value of the surfactant solution.