# The Onset Acceleration for Surfactant Covered Faraday Waves

Stephen L. Strickland1  Karen E. Daniels2  and Michael Shearer3  1Dept. of
Physics, Samford University, Homewood, AL, USA 2Dept. of Physics, NC State
University, Raleigh, NC, USA 3Dept. of Mathematics, NC State University,
Raleigh, NC, USA

###### Abstract

Faraday waves are gravity-capillary waves that emerge on the surface of a
vertically vibrated fluid when the energy injected via vibration exceeds the
energy lost due to viscous dissipation. Because this dissipation primarily
occurs in the free surface boundary layer, their emergence is particularly
sensitive to free surface properties including the surface tension,
elasticity, and viscosity of surfactants present at the free surface. We study
this sensitivity by considering a Newtonian fluid bath covered by an insoluble
surfactant subject to vertical vibrations which produce sub-harmonic Faraday
waves. By assuming a finite-depth, infinite-breadth, low-viscosity bulk fluid
and accounting for surface tension, Marangoni, and Boussinesq effects, we
derive an expression for the onset acceleration up to second order in the
expansion parameter $\Upsilon=\sqrt{\tfrac{1}{\mathcal{R}e}}$. We recover the
results of previous numerical investigations, but only by modifying the
Marangoni and Boussinesq numbers to account for the low-viscosity limit. The
analytic expression allows us to consider a range of parameters not previously
studied, including a wide variety of fluid depths and driving frequencies. In
addition, we uncover regions of parameter space for which our model predicts
that the addition of surfactant would lower, rather than elevate, the onset
acceleration. We discuss the possible use of this model in developing a
surface viscometer for surfactant monolayers.

###### keywords:

Faraday waves, Instability, Interfacial Flows

##  1 Introduction

When faced with a roaring ocean, Roman sailors would break open casks of oil,
spilling the contents into the sea, and would ride in a patch of quiet oil-
covered water until the storm abated. This calming effect of an oil layer on
ocean waves (gravity-capillary waves) has been reported by Pliney the elder
(Pliny the Elder, ), publicized by Benjamin Franklin Franklin et al. (1774),
and utilized by Shields at Aberdeen Harbor Aitken (1882) as a way of keeping
sea-faring craft safe. More recently, the dynamic effects of surface materials
on gravity-capillary waves have become useful for the remote detection of
crude oil spills Cini et al. (1983); Brekke and Solberg (2005); Ghanmi et al.
(2015), detection of biological molecules Picard and Davoust (2007, 2009),
measurement of bulk and interfacial rheology Lucassen-Reynders and Lucassen
(1970); Douady (1990); Jiang et al. (1993); Raynal et al. (1999); Saylor et
al. (2000); Behroozi et al. (2007); Shao et al. (2018); Lau et al. (2020), and
the patterning of interfaces Henderson et al. (1991); Henderson (1998); Wright
and Saylor (2003).

In these applications, a layer of surfactant reduces the surface tension
($\sigma$) of the bulk fluid by an amount that depends upon the surface
density ($\Gamma$) of the surfactant, and typically the surface tension
decreases monotonically as surfactant density increases. An inhomogeneous
distribution of surfactant will result in surface tension gradients (Marangoni
stresses) which drive flow in the bulk fluid. This flow then transports the
surfactant molecules, modifying their spatial distribution. If left
unperturbed, the coupled surfactant-fluid system would reach an equilibrium
for which the surfactant becomes homogeneously distributed and the Marangoni
stresses vanish via diffusion.

A traveling gravity-capillary wave will compress and expand the surfactant-
covered interface, giving rise to Marangoni stresses (Lange and Hühnerfuss,
1984). These stresses in turn result in a viscous boundary layer at the fluid
surface where dissipation of the wave’s energy is enhanced. The energy
dissipation is often made apparent through the exponential decay of the wave’s
amplitude as it propagates (Reynolds, 1880; Levich, 1941; Dorrestein, 1951;
Case and Parkinson, 1957; Goodrich, 1961; Davies and Vose, 1965; Lucassen-
Reynders and Lucassen, 1970; Jiang et al., 1993; Saylor et al., 2000; Behroozi
et al., 2007).

For linear small-amplitude gravity-capillary waves, the energy dissipation
rate is characterized by the damping parameter $\delta$ and is related to the
surface elastic modulus (a.k.a. surface dilational modulus or Gibbs’
elasticity), $\varepsilon_{0}=-\Gamma_{0}\tfrac{d\sigma}{d\Gamma}$ where
$\Gamma_{0}$ is the equilibrium mean surfactant density. Contemporary
theoretical and experimental research (Reynolds, 1880; Levich, 1941;
Dorrestein, 1951; Case and Parkinson, 1957; Goodrich, 1961; Davies and Vose,
1965; Miles, 1967; Lucassen-Reynders and Lucassen, 1970; Jiang et al., 1993;
Saylor et al., 2000; Behroozi et al., 2007) has shown that the damping
$\delta$ increases non-monotonically as a function of $\varepsilon_{0}$. For
$\varepsilon_{0}=0$, $\delta$ measures the bulk damping effect in the absence
of surfactant. As $\varepsilon_{0}$ is increased (typically by adding more
surfactant), $\delta$ reaches a maximum that can be an order of magnitude
larger than the surfactant-free bulk damping. For larger $\varepsilon_{0}$,
$\delta$ decreases to a value that is roughly half of its peak value
(Reynolds, 1880; Levich, 1941; Dorrestein, 1951; Case and Parkinson, 1957;
Goodrich, 1961; Davies and Vose, 1965; Miles, 1967; Lucassen-Reynders and
Lucassen, 1970; Jiang et al., 1993; Saylor et al., 2000; Behroozi et al.,
2007).

Small-amplitude standing gravity-capillary waves, known as Faraday waves
(Faraday, 1831), emerge when a fluid bath is vertically vibrated above an
angular frequency $\omega_{c}$, provided the acceleration amplitude $a$ meets
or exceeds the onset acceleration $a_{c}$ for that frequency. The emergent
wave, with angular frequency $\omega_{0}$, can either be harmonic
($\omega_{0}=\omega$) or sub-harmonic ($\omega_{0}=\tfrac{1}{2}\omega$). In
this work, we focus exclusively on sub-harmonic Faraday waves.

We understand the emergence of the Faraday waves from an energy-balance
standpoint: the fluid bath dissipates energy in every wave mode due to the
viscosity of the fluid. On the other hand, vertical vibration injects energy
into all wave modes, but not uniformly. When the driving amplitude is less
than the onset acceleration, the energy dissipation exceeds energy injection
in all modes so that no wave emerges, but at the onset acceleration $a_{c}$, a
single wave mode has more energy injected than dissipated, while all others
remain dissipative. Therefore, a wave emerges with a selected pattern of
wavenumber $k_{c}$ (Edwards and Fauve, 1994; Gollub, 2006; Ibrahim, 2015).

Because the energy dissipation largely occurs in the boundary layer near the
free surface, the onset acceleration for Faraday waves is very sensitive to
the presence of surface stresses such as due to surface tension, surface
elasticity, and surface viscosity. Because surfactants modify these stresses,
the onset of Faraday waves can serve as an effective indicator of the presence
of a surfactant as well as a means of measuring the rheological properties of
the surfactant layer. The effects of soluble surfactant on the Faraday wave
onset have been observed (Ballesta and Manneville, 2005) as has the effect of
insoluble surfactant on the damping rates of Faraday waves (Henderson et al.,
1991; Henderson, 1998).

Benjamin and Ursell (1954); Kumar and Tuckerman (1994); Chen and Viñals
(1999), and Kumar and Matar (2002, 2004a, 2004b) have made theoretical
predictions for these non-linear waves that relate the onset acceleration to
the rate of energy dissipation in the system. These theoretical predictions
for the onset acceleration are typically formulated by using Floquet analysis,
first applied by Kumar and Tuckerman (1994), which results in a recursion
relation whose truncation is often solved with numerical techniques for a pre-
specified parameter regime. This combination of Floquet analysis and numerical
solvers has been expanded to consider surfactant effects by Ubal et al.
(2005a, b, c); Giavedoni and Ubal (2007); Kumar and Matar (2004b); Mikishev et
al. (2016), who found that the onset acceleration is sensitive to the
rheological properties of the surfactant layer in much the same way that the
viscous damping parameter $\delta$ for linear gravity-capillary waves depends
upon the surface elasticity $\varepsilon_{0}$. Using a different approach (a
purely analytic technique), Chen and Viñals (1999) considered Faraday waves on
an infinite-depth surfactant-free fluid and derived an exact expression for
the onset acceleration and wave number of the emergent Faraday waves. Chen and
Viñals also started with Floquet analysis yielding a recursion relation, but
instead of solving this relation numerically, they considered the weak
viscosity limit, expanding the driving acceleration and wave number in terms
of $\gamma=\tfrac{1}{\mathcal{R}e}$ and solving for the coefficients of the
expansion.

In this paper, we extend the techniques of Chen and Viñals (1999) into the
finite-depth low-viscosity regime with surfactant, and we show that our
analysis improves upon the numerical predictions of Kumar and Matar (2004a)
and Giavedoni and Ubal (2007). In section §2, we present the parameterization
of our system and the linearized governing equations for our model. The
techniques for solving these equations are in section §3 while the general
solution and special cases are in section §4. In section §5, we compare our
analytic solution to the results of previous numerical investigations. Novel
features of the onset acceleration are in section §6, and a possible
application of this model for developing a surface viscometer is detailed in
section §7.

##  2 Parameterization & Governing Equations

Figure 1:  Schematic showing a horizontally-infinite layer of incompressible fluid with mean finite depth $H$ covered with a surfactant layer (green) subject to vertical vibrations of angular frequency $\omega$ with an amplitude $a$ driven by the container floor (black). The flat equilibrium fluid surface (dashed line) is taken as $z=0$ and the perturbation of that surface is $\zeta$. The bulk fluid is assumed to have constant density $\rho$, dynamic viscosity $\mu$, and flow velocity $\vec{u}$. The free surface have spatiotemporally-varying surfactant mass density $\Gamma$, equilibrium surface tension $\sigma_{0}$, surface elasticity $\varepsilon_{0}$, surface dilational viscosity $\Lambda$, surface shear viscosity $M$, and surface diffusivity $D$.  Symbol | Description
---|---
$\left[x,y,z\right]$ | Cartesian coordinates
$\vec{r}_{H}$ | horizontal coordinates
$H$ | bulk fluid depth
$\zeta$ | surface elevation
$\vec{u}=\left[u,v,w\right]$ | flow velocity
$\rho$ | bulk fluid density
$\mu$ | bulk fluid dynamic viscosity
$\delta$ | damping parameter
$\Gamma$ , $\Gamma_{0}$ | surfactant surface density & mean equilibrium surface density
$\sigma$ , $\sigma_{0}$ | surface tension & mean equilibrium surface tension
$\varepsilon_{0}$ | surface elastic modulus
$\Lambda$ | surface dilational viscosity
$M$ | surface shear viscosity
$\Omega=\Lambda+2M$ | combined surface viscosity
$D$ | surface diffusivity
$g$ | gravitational constant
$a$ | driving acceleration amplitude
$\omega$ | driving frequency
$\omega_{0}=\tfrac{1}{2}\omega$ | sub-harmonic Faraday wave frequency
$a_{c}$ | critical (onset) acceleration
$k_{c}$ | critical (onset) wavenumber
$k_{0}$ | wavenumber for gravity-capillary waves on an unvibrated fluid
$l_{0}=1/k_{0}$ | lengthscale
Table 1:  Physical parameters  Definition | Description
---|---
$\mathcal{F}=\frac{\displaystyle a}{\displaystyle g}$ | dimensionless acceleration amplitude
$G=\frac{\displaystyle g}{\displaystyle l_{0}\omega_{0}^{2}}$ | dimensionless gravitational acceleration
$\Sigma=\frac{\displaystyle\sigma_{0}}{\displaystyle\rho\omega_{0}^{2}l_{0}^{3}}$ | dimensionless surface tension
$\frac{\displaystyle G}{\displaystyle\Sigma}=\frac{\displaystyle\rho gl_{0}^{2}}{\displaystyle\sigma_{0}}$ | Bond number
$\mathcal{R}e=\frac{\displaystyle\rho\omega_{0}l_{0}^{2}}{\displaystyle\mu}$ | bulk Reynolds number
$\Upsilon=\sqrt{\frac{\displaystyle\mu}{\displaystyle\rho\omega_{0}l_{0}^{2}}}=\sqrt{\frac{\displaystyle 1}{\displaystyle\mathcal{R}e}}$ | the expansion parameter
$\mathcal{C}a=\frac{\displaystyle 1}{\displaystyle\mathcal{R}e\Sigma}=\frac{\displaystyle\mu\omega_{0}l_{0}}{\displaystyle\sigma_{0}}$ | capillary number
$\mathcal{M}=\frac{\displaystyle\varepsilon_{0}/\sigma_{0}}{\displaystyle\mathcal{C}a}=\frac{\displaystyle\varepsilon_{0}}{\displaystyle\mu\omega_{0}l_{0}}$ | Marangoni number
$\mathcal{M}^{\dagger}=\Upsilon\mathcal{M}=\frac{\displaystyle\varepsilon_{0}}{\displaystyle\sqrt{\rho\mu\omega_{0}^{3}l_{0}^{4}}}$ | modified Marangoni number
$\mathcal{P}e=\frac{\displaystyle\omega_{0}l_{0}^{2}}{\displaystyle D}$ | surface Peclet number
$\mathcal{B}_{D}=\frac{\displaystyle\Lambda+M}{\displaystyle\mu l_{0}}$ | dilational Boussinesq number
$\mathcal{B}_{S}=\frac{\displaystyle M}{\displaystyle\mu l_{0}}$ | shear Boussinesq number
$\mathcal{B}=\mathcal{B}_{D}+\mathcal{B}_{S}=\frac{\displaystyle\Omega}{\displaystyle\mu l_{0}}$ | combined Boussinesq number
$\mathcal{B}^{\dagger}=\Upsilon\mathcal{B}=\frac{\displaystyle\Omega}{\displaystyle\sqrt{\rho\mu\omega_{0}l_{0}^{4}}}$ | modified Boussinesq number
Table 2:  Dimensionless numbers

We consider, as shown in figure 1 and parameterized in table 1, an
incompressible Newtonian bulk fluid of infinite horizontal extent and finite
depth $H$ with flow velocity $\vec{u}(\vec{r},t)=[u,v,w]$, density $\rho$, and
dynamic viscosity $\mu$. The elevation of the free surface is given by
$z=\zeta(\vec{r}_{H},t)$ where we denote the horizontal coordinates as
$\vec{r}_{H}$. We will use the subscript H to indicate a projection onto the
horizontal plane. The dynamics of the air above the free surface are taken to
be negligible.

The surfactant monolayer at the free surface is treated as a 2-dimensional
Newtonian fluid. At equilibrium the mean surfactant mass density is uniformly
$\Gamma_{0}$ which determines the mean surface tension
$\sigma_{0}(\Gamma_{0})$, mean surface dilational viscosity
$\Lambda(\Gamma_{0})$, mean surface shear viscosity $M(\Gamma_{0})$, and mean
surface diffusivity $D(\Gamma_{0})$. When the surface is dynamic, the
instantaneous local surfactant mass density is $\Gamma(\vec{r}_{H},t)$ thereby
inducing small negligible perturbations in $\Lambda$, $M$, and $D$ away from
their equilibrium values. Variations in the surface tension
$\sigma=\sigma(\Gamma)$ are significant however; in linearizing the governing
equations, effects of surface tension are divided into the equilibrium surface
tension $\sigma_{0}$ and the surface elasticity
$\varepsilon_{0}=-\Gamma\tfrac{d\sigma}{d\Gamma}$ which will both be
considered constants.

To ensure that the unperturbed fluid surface remains at $z=0$ during
vibration, we will analyze the system in the non-inertial reference frame that
is co-vibrated with the container floor. In effect, this gives a gravitational
body-force term $\vec{g}(t)=-\left(g+a\cos(\omega t)\right)\hat{z}$ in which
$g=9.8$ m/s2 is the gravitational constant and $a$ is the acceleration
amplitude of the vibration.

In developing the model, we non-dimensionalize the governing equations by
choosing scales for time, length, and mass. We scale time by
$\tfrac{1}{\omega_{0}}$, the Faraday wave period. The length scale is set by
$l_{0}=\tfrac{1}{k_{0}}$, where $k_{0}$ is the wavenumber for a gravity-
capillary wave on an unvibrated fluid as given by the finite-depth Kelvin
dispersion relation
$1=\left(\tfrac{g}{\omega_{0}^{2}}k_{0}+\tfrac{\sigma_{0}}{\rho\omega_{0}^{2}}k_{0}^{3}\right)\tanh(k_{0}H)$.
We scale mass by $\rho l_{0}^{3}$. Velocities are scaled by $l_{0}\omega$, and
$\zeta$ and $H$ are scaled by a factor of $l_{0}$. We also scale the
surfactant density by $\Gamma_{0}$ and pressure by $\rho
l_{0}^{2}\omega_{0}^{2}$. This process results in the standard dimensionless
numbers defined in table 2. For the remainder of this manuscript, $\vec{u}$,
$p$, $k$, $H$, $\zeta$, and $\Gamma$ will refer to the dimensionless versions
of these quantities.

To analyze the system at onset, we linearize the governing equations and
boundary conditions around the equilibrium $\vec{u}=0$, $\zeta=0$, $\Gamma=1,$
which refer to as the trivial solution. The bulk fluid satisfies the
linearized incompressible Navier-Stokes equations,

$\partial_{t}\vec{u}=-\vec{\nabla}p+\frac{1}{\mathcal{R}e}\nabla^{2}\vec{u}-G\left(1+\mathcal{F}\cos(2t)\right)\hat{z}$
(1) $\vec{\nabla}\cdot\vec{u}=0$ (2)

with Reynolds number $\mathcal{R}e$, gravitation number $G$, and dimensionless
driving acceleration $\mathcal{F}$. The bulk fluid also satisfies the no-slip
boundary condition at the floor $z=-H:$

$\vec{u}=0,\ \ \ z=-H.$ (3)

At the free surface $z=\zeta(\vec{r}_{H},t)$, we have further boundary
conditions: the kinematic boundary condition,

$\partial_{t}\zeta=w,$ (4)

and the surface continuity equation, expressing the advection and diffusion of
the surfactant,

$0=\partial_{t}\Gamma+\vec{\nabla}_{H}\cdot\vec{u}_{H}-\frac{1}{\mathcal{P}e}\nabla_{S}^{2}\Gamma$
(5)

with Peclet number $\mathcal{P}e$.

Since the mass of the surfactant monolayer is negligible, the surface
tangential and normal stress boundary conditions contain no inertial or
gravitational terms for the surfactant,

$0=-\left[\begin{array}[]{c}\partial_{x}w+\partial_{z}u\\\
\partial_{y}w+\partial_{z}v\end{array}\right]-\mathcal{M}\vec{\nabla}_{H}\Gamma+\mathcal{B}_{D}\vec{\nabla}_{H}\vec{\nabla}_{H}\cdot\vec{u}_{H}+\mathcal{B}_{S}\vec{\nabla}_{H}^{2}\vec{u}_{H},$
(6) $0=\frac{1}{\Sigma}p-2\mathcal{C}a\,\partial_{z}w-\nabla_{H}^{2}\zeta$ (7)

with Marangoni number $\mathcal{M}$ and dilational and shear Boussinesq
numbers $\mathcal{B}_{D}$ and $\mathcal{B}_{S}$ respectively. General
governing equations for a surfactant-covered surface are given by Scriven
(1960), with later corrections by Waxman (1984).

The equations can be reduced and expressed in terms of $w$, $\zeta$, and
$\Gamma$ exclusively. In the reduction, the dilational and shear Boussinesq
effects can be combined, so it is convenient to define an effective surface
viscosity $\Omega=\Lambda+2M$ and an effective Boussinesq number
$\mathcal{B}=\mathcal{B}_{D}+\mathcal{B}_{S}$. The equations and boundary
conditions are then

$0=\left[\partial_{t}(\nabla_{H}^{2}+\partial_{zz})-\frac{1}{\mathcal{R}e}(\nabla_{H}^{2}+\partial_{zz})^{2}\right]w$
(8a) $w=0;\ \ \ \partial_{z}w=0,\ \ {z=-H}$ (8b) $\partial_{t}\zeta=w,\ \ \
{z=\zeta}$ (8c)
$\partial_{t}\Gamma-\partial_{z}w-\frac{1}{\mathcal{P}e}\nabla_{H}^{2}\Gamma=0,\
\ \ {z=\zeta}$ (8d)
$-\left[\nabla_{H}^{2}-\partial_{zz}\right]w-\mathcal{M}\nabla_{H}^{2}\Gamma-\mathcal{B}\nabla_{H}^{2}\partial_{z}w=0,\
\ \ {z=\zeta}$ (8e)
$\mathcal{C}a~{}\partial_{z}\left(\mathcal{R}e\partial_{t}-\left(3\nabla_{H}^{2}+\partial_{zz}\right)\right)w-\tfrac{G}{\Sigma}\left(1+\mathcal{F}\cos(2t)\right)\nabla_{H}^{2}\zeta+\nabla_{H}^{2}\nabla_{H}^{2}\zeta=0,\
\ \ {z=\zeta}$ (8f)

with Capillary number $\mathcal{C}a$, dimensionless surface tension $\Sigma$,
and Bond number $\tfrac{G}{\Sigma}$.

##  3 Technique for Finding the Onset Acceleration

The onset acceleration $\mathcal{F}_{c}=\tfrac{a_{c}}{g}$ is the minimum value
of the driving parameter $\mathcal{F}$ for which $w$, $\zeta$, and $\Gamma$
are non-trivial. As detailed in appendix A, we apply an ansatz and solve for
non-trivial $w$, $\zeta$, and $\Gamma$ up to a family of constants $\zeta_{j}$
which are the amplitudes of each wave mode. We will see that these wave mode
amplitudes are coupled in a way that allows us to solve for the onset
acceleration.

For a given wavenumber $k$, the ansatz becomes:

$\displaystyle w$
$\displaystyle=\cos(\vec{k}\cdot\vec{r}_{H})\sum_{j\in\mathbb{Z}_{\text{odd}}}ije^{ijt}\left(\mathcal{A}_{j}\sinh(kz)+\mathcal{B}_{j}\cosh(kz)+\mathcal{C}_{j}\sinh(q_{j}z)+\mathcal{D}_{j}\cosh(q_{j}z)\right)$
(9) $\displaystyle\zeta$
$\displaystyle=\cos(\vec{k}\cdot\vec{r}_{H})\sum_{j\in\mathbb{Z}_{\text{odd}}}\zeta_{j}e^{ijt}$
$\displaystyle\Gamma$
$\displaystyle=1+\cos(\vec{k}\cdot\vec{r}_{H})\sum_{j\in\mathbb{Z}_{\text{odd}}}\Gamma_{j}e^{ijt}$

where $q^{2}_{j}=k^{2}+ij\mathcal{R}e$. The series for $w$ represents the
general solution of the fourth order partial differential equation (8a), which
is satisfied by any choice of coefficients $\mathcal{A}_{j}$,
$\mathcal{B}_{j}$, $\mathcal{C}_{j}$, and $\mathcal{D}_{j}$. Equation (8d)
then gives a homogeneous equation expressing the coefficients $\Gamma_{j}$ as
linear combinations of $\mathcal{A}_{j}$, $\mathcal{B}_{j}$,
$\mathcal{C}_{j}$, and $\mathcal{D}_{j}$. Together, the homogeneous equations
(8b), (8c), and (8e) form a linear system showing that the coefficients
$\mathcal{A}_{j}$, $\mathcal{B}_{j}$, $\mathcal{C}_{j}$, and $\mathcal{D}_{j}$
not only depend on $k,$ but also are all proportional to $\zeta_{j}$.

We are left with equation (8f), in which the temporal modes are coupled due to
the forcing term $\mathcal{F}\cos(2t)$. In this term,
$\zeta_{j}e^{ijt}\cos(2t)$ splits into $\zeta_{j}e^{i(j+2)t}$ and
$\zeta_{j}e^{i(j-2)t}$, so equation (8f) becomes a linear difference equation
for the sequence of coefficients $\\{\zeta_{j}\\}$:

$0=-\zeta_{j}H_{j}+\zeta_{j-2}\mathcal{F}+\zeta_{j+2}\mathcal{F}.$ (10)

Here, the term $\zeta_{j}H_{j}$ is given by the formula:

$\zeta_{j}H_{j}=-\frac{2}{G}\left[\zeta_{j}G+\zeta_{j}\Sigma
k^{2}+\frac{ij}{k\mathcal{R}e}\left(\mathcal{A}_{j}(k^{2}+q_{j}^{2})+2kq_{j}\mathcal{C}_{j}\right)\right].$
(11)

Since $\mathcal{A}_{j},\mathcal{C}_{j}$ are proportional to $\zeta_{j}$, the
right hand side of equation (11) is also proportional to $\zeta_{j}$.
Consequently, each $H_{j}$ depends only on the wavenumber $k.$

Following Chen and Viñals (1999), to solve for the onset acceleration
$\mathcal{F}_{c}$ (the minimum $\mathcal{F}$), we truncate this difference
equations as follows. We note that when a wave of wavenumber $k$ is driven
near its onset, it oscillates sub-harmonically indicating that $\zeta_{1}$ is
the most significant contributor with nearly all other modes being negligible.
But as the driving exceeds the onset, higher order frequencies emerge making
the higher order $\zeta_{j}$ more significant. Since we want to solve for the
onset, we can consider the $\zeta_{j}$ to approach $0$ as $j\to\infty$,
suggesting that truncation at large $j=n$ will provide increasingly accurate
estimates. Starting with $\zeta_{n}$, we use the recursion relation to solve
for $\zeta_{j},\ j\leq n:$

$\displaystyle\zeta_{n}=$ $\displaystyle\zeta_{n-2}\frac{\mathcal{F}}{H_{n}}$
$\displaystyle\zeta_{n-2}=$
$\displaystyle\zeta_{n-4}\frac{\mathcal{F}}{H_{n-2}-\frac{\mathcal{F}^{2}}{H_{n}}}$
$\displaystyle\ldots$ (12)

Extending to $\zeta_{1}$, we obtain:

$\zeta_{1}=\zeta_{-1}\frac{\mathcal{F}}{H_{1}-\frac{\mathcal{F}^{2}}{H_{3}-\frac{\mathcal{F}^{2}}{H_{5}-\ldots}}}$
(13)

Because $\zeta$ (which measures the displacement of the fluid surface) is real
valued, $\zeta_{-1}=\zeta_{1}^{*}$ and therefore
$\tfrac{\zeta_{1}}{\zeta_{-1}}=e^{i\phi}$ where $\phi=2\arg(\zeta_{1})$
measures the phase difference between the driving vibration (of frequency
$\omega$) and the wave oscillation (of frequency
$\omega_{0}=\tfrac{1}{2}\omega$). Eliminating the $\zeta_{j}$, we obtain an
expression relating $\mathcal{F}$ and $k:$

$\mathcal{F}=e^{i\phi}\left(H_{1}-\frac{\mathcal{F}^{2}}{H_{3}-\frac{\mathcal{F}^{2}}{H_{5}-\ldots}}\right).$
(14)

Since the coefficients $H_{j}$ depend only on the wavenumber $k,$ this
equation defines the driving acceleration $\mathcal{F}$ implicitly as a
function of $k.$ The onset acceleration $\mathcal{F}_{c}=\mathcal{F}(k_{c})$
is the minimum of $\mathcal{F}(k),$ with onset wavenumber $k=k_{c}.$

In finding $\mathcal{F}_{c}$, one could proceed with a first derivative test
as done by Chen and Viñals (1999). Here, we instead consider that, as with a
driven damped oscillator, resonance occurs when the driver optimally injects
energy into the oscillator. Often expressed in terms of a phase difference of
$\tfrac{\pi}{2}$ between the forcing and the position of a driven oscillator,
a similar phase difference happens with Faraday waves near onset. As shown in
Figure 1 of Douady et al. (1989), the vibration acceleration and the Faraday
wave have a phase difference $\phi=\tfrac{\pi}{2}$, and on close examination,
the solution reported by Chen and Viñals (1999) exhibits this same phase
difference. Consequently, at the onset acceleration
$\mathcal{F}=\mathcal{F}_{c}$, we have
$\tfrac{\zeta_{1}}{\zeta_{-1}}=e^{i\phi}=i$, and (14) becomes:

$0=i\mathcal{F}_{c}+H_{1}-\frac{\mathcal{F}_{c}^{2}}{H_{3}-\frac{\mathcal{F}_{c}^{2}}{H_{5}-\ldots}}.$
(15)

To solve equation (15) in the case of weak bulk viscosity, we express the
onset wavenumber $k_{c}$ and onset acceleration $\mathcal{F}_{c}$ as power
series in the parameter $\Upsilon=\sqrt{\tfrac{1}{\mathcal{R}e}}$. In the
limit $\Upsilon\to 0$, we observe $k_{c}\to 1$ and $\mathcal{F}_{c}\to 0$, so
that the power series become:

$\mathcal{F}_{c}=\sum_{n=1}^{\infty}\alpha_{n}\Upsilon^{n},\ \ \ \ \ \
k_{c}=1+\sum_{n=1}^{\infty}\beta_{n}\Upsilon^{n}.$ (16)

In considering the weak viscosity limit, the behaviors of the Marangoni and
Boussinesq numbers need careful attention. A naive consideration of Table 2
would suggest that
$\mathcal{M}\sim\mathcal{B}\sim\tfrac{1}{\mu}\sim\mathcal{O}(\Upsilon^{-2})$;
however, in using the method of dominant balance (see appendix B), we find
$\mathcal{M}\sim\mathcal{B}\sim\Upsilon^{-1}$ and therefore define modified
Marangoni and modified Boussinesq numbers $\mathcal{M}^{\dagger}$ and
$\mathcal{B}^{\dagger}$ as:

$\mathcal{M}^{\dagger}=\frac{\varepsilon_{0}}{\sqrt{\mu\rho\omega_{0}^{3}l_{0}^{4}}}\
\ \ ,\ \ \
\mathcal{B}^{\dagger}=\frac{\Omega}{\sqrt{\mu\rho\omega_{0}l_{0}^{4}}}$ (17)

so that $\mathcal{M}=\Upsilon^{-1}\mathcal{M}^{\dagger}$ and
$\mathcal{B}=\Upsilon^{-1}\mathcal{B}^{\dagger}$ with
$\mathcal{M}^{\dagger}\sim\mathcal{B}^{\dagger}\sim\mathcal{O}(1)$. The
particular choice $\mathcal{M}=\Upsilon^{-1}\mathcal{M}^{\dagger}$ provides a
way of putting the surfactant effects as $\mathcal{O}(1)$ in $H_{j}$ rather
than $\mathcal{O}(\Upsilon)$. If one were to take
$\mathcal{M}=\mathcal{M}^{\dagger}$, the onset acceleration would
monotonically increase without bound as the Marangoni number increases.
Alternatively, if one were to take
$\mathcal{M}=\Upsilon^{-2}\mathcal{M}^{\dagger}$, then all of the onset
acceleration terms second order and higher would be infinite. In choosing
$\mathcal{M}=\Upsilon^{-1}\mathcal{M}^{\dagger}$, we obtain quantitative
agreement with previous numerical work as we will show in section §5.

At this point, we solve for the $\alpha_{n}$ and $\beta_{n}$ by substitution
into eqn 15 and collect terms of like order in $\Upsilon$. We find that the
complex valued equations permit real valued solutions for $\alpha_{n}$ and
$\beta_{n}$. For reference, appendix B shows the $H_{j}$ expanded in terms of
$\Upsilon$.

##  4 The solution

Before presenting the general solution for the onset acceleration, it is
reassuring to consider the solution in a few specific cases that have already
been well-studied. The original analytical work by Chen and Viñals (1999)
considered the infinite-depth surfactant-free problem, so in section §4.1, we
will show that our analysis recovers their result. In section §4.2, we will
extend the surfactant-free problem to the finite-depth limit. In section §4.3,
we will consider the infinite-depth surfactant problem since the Marangoni,
Boussinesq, and Peclet effects are easier to identify. We will then present
the complete solution for a finite-depth surfactant-covered fluid in section
§4.4. For all cases, Mathematica was used to help calculate the coefficients
of the series expansions.

###  4.1 The surfactant-free infinite-depth case

We obtain the surfactant-free infinite-depth case by letting
$\mathcal{M}^{\dagger}\rightarrow 0$, $\mathcal{B}^{\dagger}\rightarrow 0$,
and $H\rightarrow\infty$. In this case, eqn (15) becomes:

$\displaystyle 0=$ $\displaystyle-\tfrac{2}{G}\left(-1+G+\Sigma\right)$ (18)
$\displaystyle+\Upsilon\left(i\alpha_{1}-\tfrac{2}{G}\beta_{1}\left(1+2\Sigma\right)\right)$
$\displaystyle+\Upsilon^{2}\left(i\left(\alpha_{2}-\tfrac{8}{G}\right)-\tfrac{2}{G}\beta_{2}(1+2\Sigma)\right)$
$\displaystyle+\Upsilon^{3}\left(i\left(\alpha_{3}+\tfrac{4\sqrt{2}}{G}\right)-\tfrac{2}{G}(\beta_{3}(1+2\Sigma)-2\sqrt{2})\right)$
$\displaystyle+\Upsilon^{4}\left(i\alpha_{4}-\tfrac{2}{G}(\beta_{4}(1+2\Sigma)+4+\tfrac{G^{2}}{32}\alpha_{2}^{2})\right)$
$\displaystyle+\Upsilon^{5}\left(i\left(\alpha_{5}-\tfrac{2\sqrt{2}+8\beta_{3}}{G}\right)-\tfrac{2}{G}(\beta_{5}(1+2\Sigma)-\sqrt{2}+\tfrac{G^{2}}{16}\alpha_{2}\alpha_{3})\right)$
$\displaystyle+\mathcal{O}(\Upsilon^{6})$

Because $\Upsilon$ is arbitrary, each term must independently vanish. The
zeroth order term is a reiteration of the infinite-depth Kelvin dispersion
relation, and the higher order terms permit solutions for the $\alpha_{j}$ ,
$\beta_{j}$. The onset acceleration and wavenumber thus become:

$\mathcal{F}_{c}=\frac{1}{G}\left[8\Upsilon^{2}-4\sqrt{2}\Upsilon^{3}+\frac{2\sqrt{2}(11-2G)}{(3-2G)}\Upsilon^{5}+\mathcal{O}(\Upsilon^{6})\right]$
(19a)
$k_{c}=\left[1+\frac{2\sqrt{2}}{3-2G}\Upsilon^{3}-\frac{6}{3-2G}\Upsilon^{4}+\frac{3\sqrt{2}}{3-2G}\Upsilon^{5}+\mathcal{O}(\Upsilon^{6})\right]$
(19b)

The expression for the onset acceleration $\mathcal{F}_{c}$ is identical to
Chen and Viñals (1999), where $\Upsilon=\sqrt{\gamma/2}$ in their notation.
However, their expression for the wavenumber $k_{c}$ is the same only up to
third order, so the formula in Chen and Viñals (1999) is not a valid solution
to equation (15).

###  4.2 The surfactant-free finite-depth case

The surfactant-free finite-depth case is achieved by letting
$\mathcal{M}^{\dagger}\rightarrow 0$ and $\mathcal{B}^{\dagger}\rightarrow 0$
while keeping $H$ arbitrary. The onset acceleration becomes:

$\mathcal{F}_{c}=\frac{1}{G}\left[\sqrt{2}\operatorname{csch}^{2}(H)\Upsilon+4\coth(H)\frac{4\Sigma\cosh(2H)+\cosh(3H)\operatorname{csch}(H)+4H-2\Sigma}{\mathcal{Q}_{H}}\Upsilon^{2}+\mathcal{O}(\Upsilon^{3})\right]$

(20)

where $\mathcal{Q}_{H}=2\Sigma\cosh(2H)+\sinh(2H)+2H-2\Sigma$.

In the infinite-depth limit ($H\to\infty$), the first-order term vanishes, and
the second order term collapses to $\frac{8}{G}\Upsilon^{2}$, in agreement
with (19a).

###  4.3 The surfactant-covered infinite-depth case

We obtain the surfactant-covered infinite-depth case by letting
$H\rightarrow\infty$ while keeping $\mathcal{M}^{\dagger}$ and
$\mathcal{B}^{\dagger}$ arbitrary. The onset acceleration becomes:

$\displaystyle\mathcal{F}_{c}=\frac{1}{G}\left[\sqrt{2}\left(\frac{\mathcal{Q}_{S}-1+\frac{\sqrt{2}\mathcal{M}^{\dagger}}{1+\tfrac{1}{\mathcal{P}e^{2}}}}{\mathcal{Q}_{S}}\right)\Upsilon+\left(\frac{2\Sigma\mathcal{N}_{1}+\mathcal{N}_{2}}{\left(2\Sigma+1\right)\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(\mathcal{Q}_{S}\right)^{3}}\right)\Upsilon^{2}+\mathcal{O}(\Upsilon^{3})\right]$
(21)

where

$\mathcal{Q}_{S}=1+\sqrt{2}\mathcal{B}^{\dagger}+\mathcal{B}^{\dagger
2}+\frac{\mathcal{M}^{\dagger}}{\mathcal{P}e\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)}\left(\sqrt{2}+2\mathcal{B}^{\dagger}-\sqrt{2}\mathcal{P}e+\mathcal{M}^{\dagger}\mathcal{P}e\right)$

and $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ (which are polynomials of
$\mathcal{M}^{\dagger}$, $\mathcal{B}^{\dagger}$, $\tfrac{1}{\mathcal{P}e}$,
and $1+\tfrac{1}{\mathcal{P}e^{2}}$) are printed in appendix C.3.

In the surfactant-free limit ($\mathcal{M}^{\dagger}\to 0$ and
$\mathcal{B}^{\dagger}\to 0$), we find

$\displaystyle\mathcal{Q}_{S}$ $\displaystyle\sim 1$
$\displaystyle\tfrac{\mathcal{N}_{1}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}}\sim\tfrac{\mathcal{N}_{2}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}}$
$\displaystyle\sim 8$

Thus the first-order term vanishes and the coefficient of the second-order
term approaches $\frac{16\Sigma+8}{G(2\Sigma+1)}\Upsilon^{2}$. Hence, equation
(21) agrees with equation (19a) in this limit.

###  4.4 The general solution

The general expression for the onset acceleration is:

$\begin{aligned}
\mathcal{F}_{c}=&\frac{1}{G}\left[\sqrt{2}\left(\frac{\mathcal{Q}_{S}\operatorname{csch}^{2}(H)+\left(\mathcal{Q}_{S}-1+\frac{\sqrt{2}\mathcal{M}^{\dagger}}{1+\tfrac{1}{\mathcal{P}e^{2}}}\right)\coth^{2}(H)}{\mathcal{Q}_{S}}\right)\Upsilon\right.\\\
&\hskip
20.00003pt+\left(\coth(H)\frac{\cosh(2H)\left(4\Sigma\mathcal{L}_{1}+\mathcal{L}_{2}\coth(H)\right)+\cosh(3H)\operatorname{csch}(H)\mathcal{L}_{3}+4H\mathcal{L}_{4}-2\Sigma\mathcal{L}_{5}+\coth(H)\mathcal{L}_{6}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(\mathcal{Q}_{S}\right)^{3}\mathcal{Q}_{H}}\right)\Upsilon^{2}\\\
&\left.\hskip 20.00003pt+\mathcal{O}(\Upsilon^{3})\right]\end{aligned}$

(22)

where the coefficients $\mathcal{L}_{1}$, $\mathcal{L}_{2}$,
$\mathcal{L}_{3}$, $\mathcal{L}_{4}$, $\mathcal{L}_{5}$, and $\mathcal{L}_{6}$
(which are polynomials of $\mathcal{M}^{\dagger}$, $\mathcal{B}^{\dagger}$,
$\tfrac{1}{\mathcal{P}e}$, and $1+\tfrac{1}{\mathcal{P}e^{2}}$) are printed in
appendix C.4.

In the surfactant-free limit, as $\mathcal{Q}_{S}\rightarrow 1$, the first
order term approaches $\tfrac{\sqrt{2}\operatorname{csch}^{2}(H)}{G}\Upsilon$.
Noting that the $\mathcal{L}$ coefficients approach:

$\displaystyle\tfrac{\mathcal{L}_{1}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}}\sim\tfrac{\mathcal{L}_{3}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}}\sim\tfrac{\mathcal{L}_{4}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}}\sim\tfrac{\mathcal{L}_{5}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}}$
$\displaystyle\sim 4$
$\displaystyle\tfrac{\mathcal{L}_{2}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}}\sim\tfrac{\mathcal{L}_{6}}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}}$
$\displaystyle\sim 0,$

the second order term approaches
$\tfrac{4\coth(H)}{G}\tfrac{4\Sigma\cosh(2H)+\cosh(3H)\operatorname{csch}(H)+4H-2\Sigma}{\mathcal{Q}_{H}}\Upsilon^{2}$.
In this limit, the general solution agrees with eqn 20.

In the infinite-depth limit, the first order term approaches
$\tfrac{\sqrt{2}}{G\mathcal{Q}_{S}}\left(\mathcal{Q}_{S}-1+\tfrac{\sqrt{2}\mathcal{M}^{\dagger}}{1+\tfrac{1}{\mathcal{P}e^{2}}}\right)\Upsilon$.
Noting that $\mathcal{L}_{1}=\tfrac{1}{2}\mathcal{N}_{1}$ and
$\mathcal{L}_{2}=\mathcal{N}_{2}-2\mathcal{L}_{3}$, the second-order term
asymptotically approaches

$\frac{\tfrac{1}{2}e^{2H}(2\Sigma\mathcal{N}_{1}+\mathcal{N}_{2}-2\mathcal{L}_{3})+e^{2H}\mathcal{L}_{3}}{(1+\tfrac{1}{\mathcal{P}e})^{2}\mathcal{Q}_{S}^{3}\tfrac{1}{2}(2\Sigma+1)e^{2H}},$

showing the general solution to agree with eqn 21.

##  5 Comparing to previous numerical results

Symbol | Kumar and Matar (2004a) | Giavedoni and Ubal (2007)
---|---|---
$\omega_{0}$ (rad/s) | $\pi\times 60$ | $\pi\times 120$
$\rho$ (kg/m3) | $1\times 10^{3}$ | $1\times 10^{3}$
$\mu$ (kg/m/s) | $1\times 10^{-3}$ | $1\times 10^{-3}$
$\sigma_{0}$ (kg/s2) | $30\times 10^{-2}$ | $70\times 10^{-2}$
$g$ (m/s2) | $9.81$ | $9.80$
$H$ (m) | $1\times 10^{-2}$ | $1.5\times 10^{-3}$
$\varepsilon_{0}$ (kg/s2) | $0-10^{-1}$ | $10^{-6}-10^{-1}$
$D$ (m/s2) | $10^{-9}-10^{-2}$ | $10^{-9}-10^{-3}$
$\Omega$ (kg/s) | N/A | $0-10^{-3}$
Table 3:  The parameter values and ranges used to compare our general solution
to the numerical results of Kumar and Matar (2004a) and Giavedoni and Ubal
(2007). The results of this comparison are shown in figures 2 and 3.  Figure
2:  Comparison of the analytic solution (22) to numerical solutions of Kumar
and Matar (2004a). The onset acceleration $\mathcal{F}_{c}$ is plotted against
the modified Marangoni number $\mathcal{M}^{\dagger}$ for
$\mathcal{B}^{\dagger}=0$ for a range of Peclet numbers $\mathcal{P}e$. The
numerical results (finely dotted lines) are nearly indistinguishable from the
corresponding analytic solution (solid lines).  Figure 3:  Comparison of the
analytic solution (22) to numerical solutions of Giavedoni and Ubal (2007).
(a) $\mathcal{F}_{c}$ vs $\mathcal{M}^{\dagger}$ with
$\mathcal{B}^{\dagger}=0$ and $\mathcal{P}e=7.991\times 10^{4}$. The numerical
result (finely dotted black line) is very close to the analytic solution (22)
(solid red line). The first-order term is shown as a dashed red curve. The
corresponding infinite-depth analytic approximation (21) is shown as a solid
yellow curve, and the first order contribution is a dashed yellow curve. (b)
$\mathcal{F}_{c}$ vs $\mathcal{B}^{\dagger}$ with $\mathcal{M}^{\dagger}=0$
and $\mathcal{P}e=7.991\times 10^{4}$. The numerical result (finely dotted
black line) is very close to the analytic solution (22) (solid blue line). The
first-order term is shown as a dashed blue curve. The corresponding infinite-
depth analytic approximation (21) is shown as a solid cyan curve, and the
first order contribution is a dashed cyan curve. (c) $\mathcal{F}_{c}$ vs
$\mathcal{M}^{\dagger}$ for $\mathcal{B}^{\dagger}=0$ and a range of
$\mathcal{P}e.$ Numerical results are finely dotted curves, barely
distiguishable from the corresponding analytic solutions (22) shown as solid
curves. (d) $\mathcal{F}_{c}$ vs $\mathcal{M}^{\dagger}$ for
$\mathcal{P}e=7.991\times 10^{4}$ and a range of $\mathcal{B}^{\dagger}.$
Numerical results are finely dotted curves, barely distiguishable from the
corresponding analytic solutions (22) shown as solid curves.  Figure 4:
Behavior of $\mathcal{F}_{c}$ in the
$\mathcal{M}^{\dagger}$-$\mathcal{B}^{\dagger}$ plane as the Peclet number is
decreased: (a) $\mathcal{P}e=7.991\times 10^{4}$ (b) $\mathcal{P}e=7.991\times
10^{0}$ (c) $\mathcal{P}e=7.991\times 10^{-1}$ (d) $\mathcal{P}e=7.991\times
10^{-2}$. These graphs show the role that surface diffusion plays in reducing,
broadening, and moving the maximum to higher $\mathcal{M}^{\dagger}$. All of
these graphs consider a finite-depth water-like bulk fluid using the same
physical parameters as Fig. 3.

Kumar and Matar (2004a, b); Ubal et al. (2005a, c), and Giavedoni and Ubal
(2007) studied the finite-depth surfactant-covered Faraday wave problem,
examining the effect of surfactants on the onset acceleration, wave number,
and the phase shift between the surfactant distribution and the surface
topography. Kumar and Matar (2004a, b) and Ubal et al. (2005a) accounted for
Marangoni and Peclet effects while Ubal et al. (2005c) considered Boussinesq
effects only. Giavedoni and Ubal (2007) generalize results of Ubal et al.
(2005a) and Ubal et al. (2005c), accounting for all three effects. In
establishing the efficacy of our analytic solution, we compare to Kumar and
Matar (2004a, b) and Giavedoni and Ubal (2007). The parameters for these
studies are shown in Table 3, and include a wide range of surface
elasticities, viscosities, and diffusivities but only a few values for the
fluid depth, viscosity, surface tension, density, and driving frequency.

Before comparing results, it is worth contrasting our analytic approach to the
recursion relation in eqn 10 with previous numerical approaches. In our
analysis, we truncated the recursion relation at an arbitrarily large $n$,
established a base case, and Taylor-expanded the problem in the weak-viscosity
limit to second order in $\Upsilon$; the numerical approaches are based on
truncating the relation at $j=10$, casting the relation into matrix form, and
numerically solving the remaining expression as an eigenvalue problem,
deducing the onset acceleration from the eigenvalue. Neither Kumar and Matar
(2004a) nor Giavedoni and Ubal (2007) explicitly considered the weak-viscosity
limit, but based on the parameters reported, they studied $\Upsilon=0.0696$
and $\Upsilon=0.0644$ respectively, well within the weak-viscosity limit.

In order to compare results, we converted the $\mathcal{M}$ and $\mathcal{B}$
from these numerical studies first to physical parameter values
$\varepsilon_{0}$, $D$, and $\Omega$ and then into the $\mathcal{M}^{\dagger}$
and $\mathcal{B}^{\dagger}$ used here. This process requires careful
accounting, not only in the additional use of $\Upsilon$ in the definitions of
$\mathcal{M}^{\dagger}$ and $\mathcal{B}^{\dagger}$, but also in the
definitions of the length and time scales used to non-dimensionalize the
problem. Although we use the wavelength from the finite-depth Kelvin
dispersion relation as the length scale, Kumar and Matar (2004a) used the
fluid depth $H$, and Giavedoni and Ubal (2007) used
$l_{0}=\tfrac{g}{\omega^{2}}+\sqrt[3]{\tfrac{\sigma_{0}}{\rho\omega^{2}}}$.

Figure 2 shows that our analysis is in quantitative agreement with the results
of Kumar and Matar (2004a) for a wide range of $\mathcal{M}^{\dagger}$ and
$\mathcal{P}e$. These curves clearly show that when diffusion is negligible,
the onset acceleration rapidly increases with surface elasticity up to a
maximum that is significantly larger than the surfactant-free case. The onset
acceleration then decreases to nearly half its peak value, very similar to the
energy damping rate in linear gravity-capillary waves as discussed in section
§1. Surface diffusivity acts to reduce, broaden, and shift the peak to higher
values of Marangoni number.

Figure 3 shows that our analysis also agrees with the results of Giavedoni and
Ubal (2007) for five decades of $\mathcal{M}^{\dagger}$, eight decades of
$\mathcal{B}^{\dagger}$, and six decades of $\mathcal{P}e$. Figure 3 (a) and
(c) show the dependence of $\mathcal{F}_{c}$ vs $\mathcal{M}^{\dagger}$ and
$\mathcal{P}e$, and although the fluid depth, driving frequency, and surface
tension are different than in Fig. 2, the trends in the plots are the same.
Figure 3 (b) shows $\mathcal{F}_{c}$ vs $\mathcal{B}^{\dagger}$ which exhibits
a steady rise to a plateau. Figure 3 (d) shows how surface viscosity also
reduces, broadens, and moves the peak in the $\mathcal{F}_{c}$ vs
$\mathcal{M}^{\dagger}$ plot. The surface viscosity also increases the onset
acceleration at low $\mathcal{M}^{\dagger}$.

To help visualize the behavior of $\mathcal{F}_{c}$ in the
$\mathcal{M}^{\dagger}$-$\mathcal{B}^{\dagger}$ plane, we have added figure 4.
Although no new comparison is made in this figure, it uses the same physical
parameters as Fig. 3. This visualization clearly shows that
$\lim_{\mathcal{M}^{\dagger}\rightarrow\infty}\mathcal{F}_{c}=\lim_{\mathcal{B}^{\dagger}\rightarrow\infty}\mathcal{F}_{c}$.
Further, as $\mathcal{P}e\rightarrow 0$, the peak lessens, broadens, shifts to
higher values of $\mathcal{M}^{\dagger}$, and ultimately vanishes.

Figure 3 (a,b) additionally show the behavior of the first-order and second-
order terms of the analytic solution (22) and infinite-depth analytic
approximation (21). The first-order terms captures the overall trends,
particularly the maximum when $\mathcal{M}^{\dagger}\approx 1$ and the rise
when $\mathcal{B}^{\dagger}\approx 0.5$, but generally under-predict the onset
acceleration. The second-order terms correct the under-prediction.

##  6 Additional features of the analytic solution

Figure 5:  Finite-depth Kelvin dispersion relation wave speed
$c=l_{0}\omega_{0}$ in the $H$-$\omega$ plane using water-like conditions
where the values for the physical parameters are: $\rho=1\times 10^{3}$
$\tfrac{\text{kg}}{\text{m}^{3}}$, $\mu=1\times 10^{-3}$
$\tfrac{\text{kg}}{\text{m s}}$, $\sigma_{0}=70\times 10^{-3}$
$\tfrac{\text{N}}{\text{m}}$, and $g=9.8$ $\tfrac{\text{m}}{\text{s}^{2}}$.
These plots show the same dispersion relation from two different perspectives.
Fig (a) is a top view of the $H$-$\omega$ plane while fig (b) uses the same
perspective of the plane as shown in figures 6,7,9. A few regions have been
indicated: $\mathcal{R}_{G}$ are gravity waves, $\mathcal{R}_{GH}$ are depth-
restricted gravity waves, $\mathcal{R}_{\Sigma}$ are capillary waves, and
$\mathcal{R}_{\Sigma H}$ are depth-restricted capillary waves.  Figure 6:
Behavior of $\mathcal{F}_{c}$ and the expansion parameter $\Upsilon$ in the
$H$-$\omega$ plane using the same surfactant-free water-like conditions as in
figure 5. The same regions from fig 5 are evident in the behavior of
$\mathcal{F}_{c}$ in that at the boundaries of each region, $\mathcal{F}_{c}$
exhibits significant curvature. The two blue dots indicate where our
comparison with Kumar and Matar (2004a) and Giavedoni and Ubal (2007)
occurred.  Figure 7:  Behavior of $\mathcal{F}_{c}$ in the $H$-$\omega$ plane
as surface elasticity is increased. These figures consider the same water-like
conditions as in figure 5 and 6 but now with a surfactant that only affects
the surface elasticity and does not diffuse. The left column of figures show
$\mathcal{F}_{c}$ while the right column shows the ratio of $\mathcal{F}_{c}$
to $\mathcal{F}_{c0}$ (the corresponding surfactant-free onset acceleration
from fig. 6). A ratio of $10^{0}=1$ means that the onset acceleration is
indistinguishable from the surfactant-free case. The eight cyan dots indicate
locations on the $H$-$\omega$ plane corresponding to the subplots in figure 8.
Figure 8:  Behavior of $\mathcal{F}_{c}$ in the
$\mathcal{M}^{\dagger}$-$\mathcal{B}^{\dagger}$ plane for several $H$ and
$\omega$. Each plot corresponds to a cyan dot in fig 7 (a). Each column
corresponds to a driving frequency and each row corresponds to a fluid depth.
The left column has a frequency of $\omega=2\pi\times 0.25$ rad/s while the
right has $\omega=2\pi\times 10$ rad/s. The depths of each row is: (a,e)
$H=1\times 10^{-5}$ m (b,f) $H=1\times 10^{-4}$ m (c,g) $H=1\times 10^{-2}$ m
(d,h) $H=1\times 10^{0}$ m. Although all of the plots use the same coloration,
the color scaling is unique for each plot as indicated by the colorbars.
Figure 9:  Behavior of $\mathcal{F}_{c}$ in the $H$-$\omega$ plane as surface
viscosity is increased. These figures consider the same water-like conditions
as in figures 5, 6, and 7 but now with a surfactant that only affects the
surface viscosity and does not diffuse. The left column of figures show
$\mathcal{F}_{c}$ while the right column shows the ratio of $\mathcal{F}_{c}$
to $\mathcal{F}_{c0}$ (the corresponding surfactant-free onset acceleration
from fig. 6). A ratio of $10^{0}=1$ means that the onset acceleration is
indistinguishable from the surfactant-free case.

The results presented in section §5 cover a wide range of
$\mathcal{M}^{\dagger}$, $\mathcal{B}^{\dagger}$, and $\mathcal{P}e$. Table 3
shows that this range was achieved by varying the surface parameters, surface
elasticity, surface viscosity, and surface diffusivity, while keeping constant
the parameters for the bulk fluid. Because of our analytic treatment, we are
able to efficiently explore the behavior of $\mathcal{F}_{c}$ for a wide range
of surface and bulk parameters. In this section, we will consider the behavior
of $\mathcal{F}_{c}$ in the $H$-$\omega$ plane, and observe that in some
regions of this plane, increasing $\mathcal{M}^{\dagger}$ can significantly
decrease the onset acceleration. We will also explore the dependence of
$\mathcal{F}_{c}$ on $\mu$, $\rho$, and $\sigma_{0}$.

The $H$ dependence of $\mathcal{F}_{c}$ (and similarly for the $\omega$
dependence) arises in two distinct ways. Although the general solution (eqn
22) explicitly references $H$, the fluid depth also affects the finite-depth
Kelvin dispersion relation which is used to determine the characteristic
length scale $l_{0}$ which is incorporated into nearly all of the
dimensionless numbers. Consequently, probing the depth-dependence cannot be
done by merely plotting eqn 22 while holding the dimensionless numbers
constant. Similarly, $\omega$ directly contributes to the dimensionless
numbers via $\omega_{0}$, and it also affects the dimensionless numbers
through $l_{0}$ via the same dispersion relation.

To illuminate the role of the finite-depth Kelvin dispersion relation, fig 5
shows the wave speed $c=l_{0}\omega_{0}$ in the $H$-$\omega$ plane for
terrestrial water ($\rho=1000$ kg/m3, $\mu=0.001$ kg/m/s, $\sigma_{0}=0.07$
N/m, $g=9.8$ m/s2). This dispersion relation strictly applies to linear
gravity-capillary waves, but in the limit $\Upsilon\rightarrow 0$, the
dispersion relation for non-linear Faraday waves approaches these plots. These
plots clearly show the gravity wave region $\mathcal{R}_{G}$ and the capillary
wave region $\mathcal{R}_{\Sigma}$ as well as two new regions which we refer
to as depth-restricted gravity waves $\mathcal{R}_{GH}$ and depth-restricted
capillary waves $\mathcal{R}_{\Sigma H}$ since the finite depth of the
container results in slower waves. This dispersion relation is so significant
that all of these regions are apparent in the behavior of $\mathcal{F}_{c}$ in
the $H$-$\omega$ plane.

Figure 6 (a) shows $\mathcal{F}_{c}$ in the $H$-$\omega$ plane for surfactant-
free water. In later figures, we use this surfactant-free onset acceleration
as a reference and denote it as $\mathcal{F}_{c0}$. In this figure, we have
marked the four regions from figure 5 on the plot. Although the most eye-
catching feature is the valley that traces along the
$\mathcal{R}_{G}$-$\mathcal{R}_{GH}$ border, each border either coincides with
or is next to significant curvature. Since the plot is logarithmic on all
axes, any planar surfaces indicate power-law behavior, and any curvature from
a plane corresponds to a change in the exponent. Figure 6 (b) shows the low-
viscosity expansion parameter $\Upsilon$ in the same $H$-$\omega$ plane. In
considering water, all points shown on the $H$-$\omega$ plane have
$\Upsilon<1$ with the maximum value of $\Upsilon=0.3458$ occurring at the
shallowest depth and highest driving frequency. The two blue dots in figure 6
show the locations in the $H$-$\omega$ plane where we compared our analysis to
Kumar and Matar (2004a) and Giavedoni and Ubal (2007) which lie in the
$\mathcal{R}_{\Sigma}$ and $\mathcal{R}_{\Sigma H}$ regions respectively.

Figure 7 shows the effect of surface elasticity on the onset acceleration by
showing a progression of graphs of $\mathcal{F}_{c}$ and
$\tfrac{\mathcal{F}_{c}}{\mathcal{F}_{c0}}$ in the $H$-$\omega$ plane. For
small elasticities, a new wedge-shaped region appears at shallow depths and
mid-range driving frequencies. We will refer to this region as the elasticity-
affected region $\mathcal{R}_{\varepsilon}$. As the surface elasticity
increases, $\mathcal{R}_{\varepsilon}$ descends the graph, pressing towards
deeper depths. The tip of the wedge follows the
$\mathcal{R}_{GH}$-$\mathcal{R}_{\Sigma H}$ boundary all the way to the
quadruple point where the four types of waves meet. Within
$\mathcal{R}_{\varepsilon}$, the onset acceleration are elevated above the
surfactant-free behavior, but at the boundary, the onset acceleration may
either increase or decrease. Figure 7 shows that the boundary tends to
decrease the onset acceleration for shallow systems and increase the onset
acceleration for deep systems. Figure 8 shows this unusual effect of surface
elasticity in more detail by plotting $\mathcal{F}_{c}$ in the
$\mathcal{M}^{\dagger}$-$\mathcal{B}^{\dagger}$ plane for eight different
locations of the $H$-$\omega$ plane, the $H$-$\omega$ location of each plot is
shown in figure 7 (a). Each column corresponds to a driving frequency and each
row corresponds to a fluid depth. Notably, the application of a surfactant can
decrease the onset acceleration by more than an order of magnitude. Although
the effects of diffusion are not shown in fig 7, increasing the diffusivity
only results in a lessening of the Marangoni-induced extremes (peaks and
valleys) and a broadening of the boundary around $\mathcal{R}_{\varepsilon}$.
Diffusivity does not affect the dependence of the onset acceleration on the
Boussinesq number.

Figure 9 shows the effect of surface viscosity on the onset acceleration by
showing a progression of graphs of $\mathcal{F}_{c}$ and
$\tfrac{\mathcal{F}_{c}}{\mathcal{F}_{c0}}$ in the $H$-$\omega$ plane. Unlike
surface elasticity, increasing surface viscosity will only increase the onset
acceleration. Weak surface viscosities will elevate the onset acceleration at
high frequency across all depths, and further increasing the surface viscosity
shifts these effects to lower frequencies.

Figure 10:  Dependence of $\mathcal{F}_{c}$ on the bulk viscosity. Similar to
the parameters used in Giavedoni and Ubal (2007), the system parameters are
$\rho=1\times 10^{3}$ kg/m3, $\sigma_{0}=70\times 10^{-3}$ N/m, and $g=9.8$
m/s2. (a,b) consider depth-restricted gravity waves where $H=1\times 10^{-3}$
m and $\omega=2\pi\times 0.25$ rad/s while (c,d) consider depth-restricted
capillary waves where $H=1.5\times 10^{-3}$ m and $\omega=2\pi\times 120$
rad/s. (a,c) consider a range of $\varepsilon_{0}$ with $\Omega=0$ while (b,d)
consider a range of $\Omega$ with $\varepsilon_{0}=0$. The limiting case of
surfactant-free is shown by a dashed black line while the case of a
surfactant-saturated surface is shown by a red dashed line.  Figure 11:
Dependence of $\mathcal{F}_{c}$ on the surface tension. Similar to the
parameters used in Giavedoni and Ubal (2007), the system parameters are
$\rho=1\times 10^{3}$ kg/m3, $\mu=1\times 10^{-3}$ kg/m/s, and $g=9.8$ m/s2.
(a,b) consider depth-restricted gravity waves where $H=1\times 10^{-3}$ m and
$\omega=2\pi\times 0.25$ rad/s while (c,d) consider depth-restricted capillary
waves where $H=1.5\times 10^{-3}$ m and $\omega=2\pi\times 120$ rad/s. (a,c)
consider a range of $\varepsilon_{0}$ with $\Omega=0$ while (b,d) consider a
range of $\Omega$ with $\varepsilon_{0}=0$. The limiting case of surfactant-
free is shown by a dashed black line while the case of a surfactant-saturated
surface is shown by a red dashed line.  Figure 12:  Dependence of
$\mathcal{F}_{c}$ on the bulk density. Similar to the parameters used in
Giavedoni and Ubal (2007), the system parameters are $\mu=1\times 10^{-3}$
kg/m/s, $\sigma_{0}=70\times 10^{-3}$ N/m, and $g=9.8$ m/s2. (a,b) consider
depth-restricted gravity waves where $H=1\times 10^{-3}$ m and
$\omega=2\pi\times 0.25$ rad/s while (c,d) consider depth-restricted capillary
waves where $H=1.5\times 10^{-3}$ m and $\omega=2\pi\times 120$ rad/s. (a,c)
consider a range of $\varepsilon_{0}$ with $\Omega=0$ while (b,d) consider a
range of $\Omega$ with $\varepsilon_{0}=0$. The limiting case of surfactant-
free is shown by a dashed black line while the case of a surfactant-saturated
surface is shown by a red dashed line.

Figures 10, 11, and 12 show the effects of bulk viscosity, surface tension,
and bulk density respectively. The parameters are based on water where
$\mu=10^{-3}$ kg/m/s, $\sigma_{0}=70\times 10^{-3}$ N/m, and $\rho=10^{3}$
kg/m3, and in each figure, we vary a single parameter. Frames (a,b) of these
figures consider depth-restricted gravity waves ($\omega=2\pi\times 0.25$
rad/s and $H=10^{-3}$ m), while frames (c,d) consider depth-restricted
capillary waves ($\omega=2\pi\times 120$ rad/s and $H=1.5\times 10^{-3}$ m).

Figure 10 graphs the onset acceleration $\mathcal{F}_{c}$ vs the bulk
viscosity $\mu$ showing that increasing the bulk viscosity generally increases
the onset acceleration. Although the Boussinesq effects shown in frames (b,d)
always increase the onset acceleration, the Marangoni effects shown in frames
(a,c) will sometimes decrease the onset acceleration as noted before. Further,
the presence of surface elasticity can result in cases where a more viscous
fluid could have a lower onset acceleration than a less viscous fluid. These
figures also show the extreme cases of a surfactant-free fluid (black dashed
line) and a surfactant-saturated fluid (red dashed line). For gravity waves, a
surfactant-saturated fluid always has a higher onset acceleration than a
surfactant-free fluid; however, for capillary waves, there is a crossover
where a saturated surface will onset Faraday waves more readily than a
surfactant-free surface. This crossover is for large viscosities, and since
this model is designed for the low-viscosity limit and only tested against
numerical results using a viscosity of $10^{-3}$ kg/m/s, the crossover may be
a limitation of our second-order analysis. However, Suman and Kumar (2008)
reported similar numerical results for high-viscosity systems, finding that
the onset acceleration for an inertial-less surfactant-free system would be
infinite, but Marangoni stresses allow Faraday waves to emerge, thereby
preferentially decreasing the onset acceleration.

Fig 11 graphs the onset acceleration $\mathcal{F}_{c}$ vs equilibrium surface
tension $\sigma_{0}$. The gravity-waves in frames (a,b) are unaffected by
surface tension while the capillary waves in frames (c,d) do respond to the
surface tension. In comparing our results with Giavedoni and Ubal (2007), we
used a surface tension of $70\times 10^{-3}$ N/m where increasing surface
elasticity and surface viscosity would preferentially increase the onset
acceleration; however upon decreasing the surface tension, we find another
crossover where increasing these surface parameters will decrease the onset
acceleration.

For completeness, we have included figure 12 which graphs the onset
acceleration $\mathcal{F}_{c}$ vs bulk fluid density $\rho$. The onset
acceleration generally decreases as bulk density increases.

##  7 Conclusion and Discussion

We have derived an analytic expression for the onset acceleration for Faraday
waves in a finite-depth infinite-breadth low-viscosity surfactant-covered
fluid. We have shown that this expression accurately reproduces the results of
previous numerical works. Our analysis required a novel definition of the
Marangoni and Boussinesq numbers to handle the low-viscosity limit as the
standard definition result in unbounded behavior for $\mathcal{F}_{c}$.

We have also shown that for shallow systems, the model model makes an
unexpected prediction: adding a surfactant to a shallow system can lower the
onset acceleration for Faraday waves. In context of the energy-balance
perspective of the emergence of Faraday waves, one would expect that by
increasing the surface elasticity, which introduces a new viscous boundary
layer at the free surface, one would increase the viscous dissipation and
thereby increase the onset acceleration. However, there are cases where
increasing the bulk viscosity in the presence of a surfactant reduces the
onset acceleration. These unexpected results may be related to the work by
Suman and Kumar (2008) where Marangoni effects in an inertial-less system act
to destabilize the system.

We conclude by noting the potential utility of our analysis in determining the
surface rheology of a surfactant monolayer. Despite the myriad of surface
rheometers that utilize macroscopic systems (Fuller and Vermant, 2012;
Jaensson and Vermant, 2018) and microscopic systems (Samaniuk and Vermant,
2014), transverse and longitudinal capillary waves have long been used to
probe the surface dilational viscosity (Lemaire and Langevin, 1992; Buzza et
al., 1998; Saylor et al., 2000), an historically difficult measurement.
Measuring the onset of Faraday waves is ideal for accessing the dilational
viscosity since (i) no mechanical probe is introduced to the system’s surface
and (ii) at onset there are minimal surfactant concentration gradients, two
key challenges that plague other measurement techniques (Fuller and Vermant,
2012). Further, the detection of Faraday waves requires a minimum of technical
equipment.

Our analysis works effectively in finite-depth systems where the bulk fluid
has a viscosity comparable to water (or less). With our general solution, one
could measure the onset acceleration for a range of frequencies and fit for
the surface rheological parameters. Marangoni and Boussinesq effects have
different frequency dependencies and can be readily distinguished. In fitting,
one would determine $\Omega=\Lambda+2M$ rather than the surface dilational
viscosity itself; however, with a surface shear viscometer, one could then
deduce the dilational viscosity using Faraday waves.

##  8 Acknowledgments

We would also like to thank Lake Bookman for the many helpful discussions in
the early attempts to formulate the theoretical framework. We would also like
to thank the NSF for grant # NSF DMS-0604047 and NSF DMS-0968258.

Declaration of Interests: The authors report no conflict of interest.

##  Appendix A Integrating the Governing Equations

In §2, we presented the governing equations for our model (eqns 8). Here, we
solve for $w$, $\zeta$, and $\Gamma$ at the moment that the Faraday waves
emerge, the moment when these functions become non-trivial. In this appendix,
we will present an ansatz and solve for all but the final constants of
integration, the wave mode amplitudes $\zeta_{j}$. The analysis of these final
constants of integration is addressed in §3 and will yield an expression for
the onset acceleration.

We use the following dimensionless ansatz:

$\begin{split}w&=\cos(\vec{k}\cdot\vec{r}_{H})\sum_{j\in\mathbb{Z}_{\text{odd}}}ijw_{j}(z)e^{ijt}\\\
\zeta&=\cos(\vec{k}\cdot\vec{r}_{H})\sum_{j\in\mathbb{Z}_{\text{odd}}}\zeta_{j}e^{ijt}\\\
\Gamma&=1+\cos(\vec{k}\cdot\vec{r}_{H})\sum_{j\in\mathbb{Z}_{\text{odd}}}\Gamma_{j}e^{ijt}\end{split}$

The dimensionless wave number $\vec{k}$ is not identically $1$ as the Faraday
wavenumber is not equal to the wavenumber from the Kelvin dispersion relation.

With this ansatz, eqn 8a yields a family of 4th order linear ODEs for the
$w_{j}(z)$ which can be readily solved.

$\left[ij\mathcal{R}e\left(-k^{2}+\partial_{zz}\right)-\left(-k^{2}+\partial_{zz}\right)^{2}\right]w_{j}(z)=0$
$w_{j}(z)=\mathcal{A}_{j}\sinh(kz)+\mathcal{B}_{j}\cosh(kz)+\mathcal{C}_{j}\sinh(q_{j}z)+\mathcal{D}_{j}\cosh(q_{j}z)$

where $q_{j}^{2}=k^{2}+ij\mathcal{R}e$. The coefficients $\mathcal{A}_{j}$,
$\mathcal{B}_{j}$, $\mathcal{C}_{j}$, and $\mathcal{D}_{j}$ are constants of
integration which will be proportional to $\zeta_{j}$.

The surfactant continuity equation (eqn 8d) yields the coefficients of the
surfactant distribution $\Gamma_{j}$.

$\Gamma_{j}=\frac{k\mathcal{A}_{j}+q_{j}\mathcal{C}_{j}}{ij+\frac{k^{2}}{\mathcal{P}e}}$

Eqns 8b, 8c, and 8e yield a system of equations for $\mathcal{A}_{j}$,
$\mathcal{B}_{j}$, $\mathcal{C}_{j}$, and $\mathcal{D}_{j}$.

$\displaystyle 0$
$\displaystyle=\mathcal{A}_{j}\sinh(-kH)+\mathcal{B}_{j}\cosh(-kH)+\mathcal{C}_{j}\sinh(-q_{j}H)+\mathcal{D}_{j}\cosh(-q_{j}H)$
$\displaystyle 0$
$\displaystyle=k\mathcal{A}_{j}\cosh(-kH)+k\mathcal{B}_{j}\sinh(-kH)+q_{j}\mathcal{C}_{j}\cosh(-q_{j}H)+q_{j}\mathcal{D}_{j}\sinh(-q_{j}H)$
$\displaystyle 0$ $\displaystyle=\zeta_{j}-\mathcal{B}_{j}-\mathcal{D}_{j}$
$\displaystyle 0$
$\displaystyle=\zeta_{j}k^{2}+\left(k\mathcal{A}_{j}+q_{j}\mathcal{C}_{j}\right)\mathcal{S}_{j}+(k^{2}\mathcal{B}_{j}+q_{j}^{2}\mathcal{D}_{j})$

The solutions are:

$\displaystyle\mathcal{A}_{j}$
$\displaystyle=\zeta_{j}\frac{\mathcal{S}_{j}q_{j}\mathcal{P}_{1j}-2k^{2}q_{j}\mathcal{P}_{3j}-(k^{2}+q_{j}^{2})\left(k\mathcal{P}_{4j}-q_{j}\right)}{\mathcal{Q}}$
$\displaystyle\mathcal{B}_{j}$
$\displaystyle=-\zeta_{j}\frac{\mathcal{S}_{j}q_{j}\left(k\left(1-\mathcal{P}_{3j}\right)-q_{j}\mathcal{P}_{4j}\right)-(k^{2}+q_{j}^{2})\mathcal{P}_{2j}}{\mathcal{Q}}$
$\displaystyle\mathcal{C}_{j}$
$\displaystyle=-\zeta_{j}\frac{\mathcal{S}_{j}k\mathcal{P}_{1j}-2k^{2}\left(k-q_{j}\mathcal{P}_{4j}\right)+k(k^{2}+q_{j}^{2})\mathcal{P}_{3j}}{\mathcal{Q}}$
$\displaystyle\mathcal{D}_{j}$
$\displaystyle=\zeta_{j}\frac{\mathcal{S}_{j}k\left(k\mathcal{P}_{4j}-q_{j}\left(1-\mathcal{P}_{3j}\right)\right)-2k^{2}\mathcal{P}_{2j}}{\mathcal{Q}}$

where

$\displaystyle\mathcal{S}_{j}$
$\displaystyle=k^{2}\left(\frac{\mathcal{M}}{ij+\frac{k^{2}}{\mathcal{P}e}}+\mathcal{B}\right)$
$\displaystyle\mathcal{P}_{1j}$ $\displaystyle=q_{j}\tanh(Hq_{j})-k\tanh(Hk)$
$\displaystyle\mathcal{P}_{2j}$ $\displaystyle=q_{j}\tanh(Hk)-k\tanh(Hq_{j})$
$\displaystyle\mathcal{P}_{3j}$
$\displaystyle=\operatorname{sech}(Hk)\operatorname{sech}(Hq_{j})$
$\displaystyle\mathcal{P}_{4j}$ $\displaystyle=\tanh(Hk)\tanh(Hq_{j})$
$\displaystyle\mathcal{Q}$
$\displaystyle=-(k^{2}-q_{j}^{2})\mathcal{P}_{2j}+\mathcal{S}_{j}\left(\left(k^{2}+q_{j}^{2}\right)\mathcal{P}_{4j}-2kq_{j}\left(1-\mathcal{P}_{3j}\right)\right)$

In compiling all of these steps, we obtain the ansatz listed in eqn 9 where
the wave amplitude $\zeta_{j}$ is the only remaining unsolved constant of
integration.

##  Appendix B Dominant Balance in the Weak-Viscosity Limit

In §3 we developed the central problem of this manuscript and mapped the route
to the solution by considering a weak-viscosity fluid. In expanding the
pertinent quantities in terms of the expansion parameter
$\Upsilon=\sqrt{\tfrac{1}{\mathcal{R}e}}$, we noted that that Marangoni and
Boussinesq numbers had to be rescaled per the method of dominant balance.
Here, we detail the expansions for the pertinent quantities and their
asymptotic behavior in the weak-viscosity limit.

The quantities $\mathcal{F}_{c}$, $k_{c}$, $q_{j}$, and $\mathcal{S}_{j}$ are
expanded as:

$\displaystyle\mathcal{F}_{c}$
$\displaystyle=\sum_{n=1}^{\infty}\alpha_{n}\Upsilon^{n}\sim\mathcal{O}(\Upsilon)$
$\displaystyle k_{c}$
$\displaystyle=1+\sum_{n=1}^{\infty}\beta_{n}\Upsilon^{n}\sim\mathcal{O}(1)$
$\displaystyle q_{j}^{2}$
$\displaystyle=k^{2}+ij\tfrac{1}{\Upsilon^{2}}\sim\mathcal{O}(\tfrac{1}{\Upsilon^{2}})$
$\displaystyle\mathcal{S}_{j}$
$\displaystyle=\frac{k^{2}}{\Upsilon}\left(\frac{\mathcal{M}^{\dagger}}{ij+\frac{k^{2}}{\mathcal{P}e}}+\mathcal{B}^{\dagger}\right)\sim\mathcal{O}(\tfrac{1}{\Upsilon})$

In the low-viscosity limit,

$\displaystyle\tanh(Hq_{j})$ $\displaystyle\sim
1+\mathcal{O}(e^{-\tfrac{H}{\Upsilon}})\rightarrow 1$
$\displaystyle\cosh(Hq_{j})$
$\displaystyle\sim\sinh(Hq_{j})\sim\mathcal{O}(e^{\tfrac{H}{\Upsilon}})\rightarrow\infty$
$\displaystyle\operatorname{sech}(Hq_{j})$
$\displaystyle\sim\operatorname{csch}(Hq_{j})\sim\mathcal{O}(e^{-\tfrac{H}{\Upsilon}})\rightarrow
0$

which simplifies the parameters $\mathcal{P}_{1j}$, $\mathcal{P}_{2j}$,
$\mathcal{P}_{3j}$, and $\mathcal{P}_{4j}$:

$\displaystyle\mathcal{P}_{1j}$ $\displaystyle\rightarrow
q_{j}-k\tanh(Hk)\sim\mathcal{O}(\tfrac{1}{\Upsilon})$
$\displaystyle\mathcal{P}_{2j}$ $\displaystyle\rightarrow
q_{j}\tanh(Hk)-k\sim\mathcal{O}(\tfrac{1}{\Upsilon})$
$\displaystyle\mathcal{P}_{3j}$ $\displaystyle\rightarrow 0\sim\mathcal{O}(0)$
$\displaystyle\mathcal{P}_{4j}$
$\displaystyle\rightarrow\tanh(Hk)\sim\mathcal{O}(1)$

The constants of integration can then be expressed as:

$\displaystyle\mathcal{Q}$
$\displaystyle=\mathcal{S}_{j}\left(\left(k^{2}+q_{j}^{2}\right)\mathcal{P}_{4j}-2kq_{j}\right)+\tfrac{ij}{\Upsilon^{2}}\mathcal{P}_{2j}\sim\mathcal{O}(\tfrac{1}{\Upsilon^{3}})$
$\displaystyle\frac{\mathcal{A}_{j}}{\zeta_{j}}$
$\displaystyle=\frac{\mathcal{S}_{j}q_{j}+(k^{2}+q_{j}^{2})}{\mathcal{Q}}\mathcal{P}_{1j}\sim\mathcal{O}(1)$
$\displaystyle\frac{\mathcal{B}_{j}}{\zeta_{j}}$
$\displaystyle=\frac{\mathcal{S}_{j}q_{j}+(k^{2}+q_{j}^{2})}{\mathcal{Q}}\mathcal{P}_{2j}\sim\mathcal{O}(1)$
$\displaystyle\frac{\mathcal{C}_{j}}{\zeta_{j}}$
$\displaystyle=\frac{\mathcal{D}_{j}}{\zeta_{j}}=-\frac{\mathcal{S}_{j}k\mathcal{P}_{1j}+2k^{2}\mathcal{P}_{2j}}{\mathcal{Q}}\sim\mathcal{O}(\Upsilon)$

The coupling coefficient $H_{j}$ becomes:

$H_{j}=-\frac{2}{G}\left[G+\Sigma
k^{2}-\frac{j^{2}\mathcal{S}_{j}q_{j}\mathcal{P}_{1j}}{k\mathcal{Q}}+\frac{ij\Upsilon^{2}}{k\mathcal{Q}}\left((k^{2}+q_{j}^{2})^{2}\mathcal{P}_{1j}-4k^{3}q_{j}\mathcal{P}_{2j}\right)\right]$

In this form, one can easily check that $H_{j}\sim\mathcal{O}(1)$. As
mentioned in §3, the definitions of $\mathcal{M}^{\dagger}$ and
$\mathcal{B}^{\dagger}$ ensure that the surfactant effects in the third term
$\tfrac{j^{2}\mathcal{S}_{j}q_{j}\mathcal{P}_{1j}}{k\mathcal{Q}}$ are
$\mathcal{O}(1)$.

##  Appendix C Onset Acceleration

Here we present the expressions for the onset acceleration of all the cases
given in §4. We also fully define all of the coefficients in the expressions.
These expressions were calculated with the aid of Mathematica.

###  C.1 The surfactant-free infinite-depth limit

$\mathcal{F}_{c}=\frac{1}{G}\left[8\Upsilon^{2}-4\sqrt{2}\Upsilon^{3}+\frac{2\sqrt{2}(11-2G)}{(3-2G)}\Upsilon^{5}+\mathcal{O}(\Upsilon^{6})\right]$
(23)

###  C.2 The surfactant-free finite-depth limit

$\mathcal{F}_{c}=\frac{1}{G}\left[\sqrt{2}\operatorname{csch}^{2}(H)\Upsilon+\frac{4\coth(H)\left(4\Sigma\cosh(2H)+\cosh(3H)\operatorname{csch}(H)+4H-2\Sigma\right)}{\mathcal{Q}_{H}}\Upsilon^{2}+\mathcal{O}(\Upsilon^{3})\right]$

(24)

where

$\displaystyle\mathcal{Q}_{H}$
$\displaystyle=2\Sigma\cosh(2H)+\sinh(2H)+2H-2\Sigma$

###  C.3 The infinite-depth surfactant limit

$\displaystyle\mathcal{F}_{c}=\frac{1}{G}\left[\sqrt{2}\left(\frac{\mathcal{Q}_{S}-1+\frac{\sqrt{2}\mathcal{M}^{\dagger}}{1+\tfrac{1}{\mathcal{P}e^{2}}}}{\mathcal{Q}_{S}}\right)\Upsilon+\left(\frac{2\Sigma\mathcal{N}_{1}+\mathcal{N}_{2}}{\left(2\Sigma+1\right)\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(\mathcal{Q}_{S}\right)^{3}}\right)\Upsilon^{2}+\mathcal{O}(\Upsilon^{3})\right]$
(25)

where the constants $\mathcal{Q}_{S}$, $\mathcal{N}_{1}$, and
$\mathcal{N}_{2}$ are below:

$\displaystyle\mathcal{Q}_{S}$
$\displaystyle=1+\sqrt{2}\mathcal{B}^{\dagger}+\mathcal{B}^{\dagger
2}+\frac{\mathcal{M}^{\dagger}}{\mathcal{P}e\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)}\left(\sqrt{2}+2\mathcal{B}^{\dagger}-\sqrt{2}\mathcal{P}e+\mathcal{M}^{\dagger}\mathcal{P}e\right)$
$\displaystyle\mathcal{N}_{1}$
$\displaystyle=\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(8+20\sqrt{2}\mathcal{B}^{\dagger}+48\mathcal{B}^{\dagger
2}+34\sqrt{2}\mathcal{B}^{\dagger 3}+30\mathcal{B}^{\dagger
4}+8\sqrt{2}\mathcal{B}^{\dagger 5}+2\mathcal{B}^{\dagger 6}\right)$
$\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{M}^{\dagger}\left(-20\sqrt{2}-68\mathcal{B}^{\dagger}-54\sqrt{2}\mathcal{B}^{\dagger
2}-44\mathcal{B}^{\dagger 3}-8\sqrt{2}\mathcal{B}^{\dagger 4}\right)$
$\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{M}^{\dagger
2}\left(48+102\sqrt{2}\mathcal{B}^{\dagger}+180\mathcal{B}^{\dagger
2}+80\sqrt{2}\mathcal{B}^{\dagger 3}+30\mathcal{B}^{\dagger 4}\right)$
$\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger
2}\left(-48\sqrt{2}\mathcal{B}^{\dagger}-120\mathcal{B}^{\dagger
2}-64\sqrt{2}\mathcal{B}^{\dagger 3}-24\mathcal{B}^{\dagger 4}\right)$
$\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger
3}\left(-54\sqrt{2}-132\mathcal{B}^{\dagger}-48\sqrt{2}\mathcal{B}^{\dagger
2}\right)+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger
4}\left(30+40\sqrt{2}\mathcal{B}^{\dagger}+30\mathcal{B}^{\dagger 2}\right)$
$\displaystyle\hskip 10.00002pt+\mathcal{M}^{\dagger
3}\left(20\sqrt{2}+88\mathcal{B}^{\dagger}+32\sqrt{2}\mathcal{B}^{\dagger
2}\right)+\mathcal{M}^{\dagger
4}\left(-32\sqrt{2}\mathcal{B}^{\dagger}-24\mathcal{B}^{\dagger
2}\right)+\mathcal{M}^{\dagger 5}\left(-8\sqrt{2}\right)+\mathcal{M}^{\dagger
6}\left(2\right)$ $\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{M}^{\dagger}\tfrac{1}{\mathcal{P}e}\left(+20\sqrt{2}+96\mathcal{B}^{\dagger}+102\sqrt{2}\mathcal{B}^{\dagger
2}+120\mathcal{B}^{\dagger 3}+40\sqrt{2}\mathcal{B}^{\dagger
4}+12\mathcal{B}^{\dagger 5}\right)$ $\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger
2}\tfrac{1}{\mathcal{P}e}\left(-68-108\sqrt{2}\mathcal{B}^{\dagger}-132\mathcal{B}^{\dagger
2}-32\sqrt{2}\mathcal{B}^{\dagger 3}\right)$ $\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger
3}\tfrac{1}{\mathcal{P}e}\left(34\sqrt{2}+120\mathcal{B}^{\dagger}+80\sqrt{2}\mathcal{B}^{\dagger
2}+40\mathcal{B}^{\dagger 3}\right)$ $\displaystyle\hskip
10.00002pt+\mathcal{M}^{\dagger
3}\tfrac{1}{\mathcal{P}e}\left(20\sqrt{2}-32\sqrt{2}\mathcal{B}^{\dagger
2}-16\mathcal{B}^{\dagger 3}\right)+\mathcal{M}^{\dagger
4}\tfrac{1}{\mathcal{P}e}\left(-44-32\sqrt{2}\mathcal{B}^{\dagger}\right)+\mathcal{M}^{\dagger
5}\tfrac{1}{\mathcal{P}e}\left(8\sqrt{2}+12\mathcal{B}^{\dagger}\right)$
$\displaystyle\mathcal{N}_{2}$
$\displaystyle=\mathcal{N}_{1}+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\mathcal{B}^{\dagger
3}\left(2\sqrt{2}+4\mathcal{B}^{\dagger}\right)$
$\displaystyle+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}\left(-4+2\sqrt{2}\mathcal{M}^{\dagger}-4\sqrt{2}\mathcal{B}^{\dagger}+12\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}-8\mathcal{B}^{\dagger
2}-6\sqrt{2}\mathcal{B}^{\dagger 3}\right)$
$\displaystyle+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}\left(16\sqrt{2}\mathcal{M}^{\dagger}-16\mathcal{M}^{\dagger
2}+24\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}-24\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger
2}+8\mathcal{B}^{\dagger 2}+4\sqrt{2}\mathcal{B}^{\dagger 3}\right)$
$\displaystyle+\mathcal{M}^{\dagger
2}\left(-4\sqrt{2}\mathcal{M}^{\dagger}+8\mathcal{M}^{\dagger
2}-2\sqrt{2}\mathcal{M}^{\dagger
3}-16\sqrt{2}\mathcal{B}^{\dagger}+16\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}-24\mathcal{B}^{\dagger
2}+20\sqrt{2}\mathcal{B}^{\dagger 2}\mathcal{M}^{\dagger}\right)$
$\displaystyle+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}\tfrac{1}{\mathcal{P}e}\left(4\sqrt{2}\mathcal{B}^{\dagger}+12\mathcal{B}^{\dagger
2}\right)$
$\displaystyle+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}\tfrac{1}{\mathcal{P}e}\left(-4\sqrt{2}\mathcal{M}^{\dagger}+4\mathcal{M}^{\dagger
2}+4\sqrt{2}\mathcal{B}^{\dagger}-20\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}+8\mathcal{B}^{\dagger
2}-20\sqrt{2}\mathcal{B}^{\dagger 2}\mathcal{M}^{\dagger}\right)$
$\displaystyle+\mathcal{M}^{\dagger
2}\tfrac{1}{\mathcal{P}e}\left(-8+8\sqrt{2}\mathcal{M}^{\dagger}-4\mathcal{M}^{\dagger
2}-8\sqrt{2}\mathcal{B}^{\dagger}+32\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}-12\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger
2}+16\mathcal{B}^{\dagger 2}+8\sqrt{2}\mathcal{B}^{\dagger 3}\right)$

###  C.4 The general solution

$\begin{aligned}
\mathcal{F}_{c}=&\frac{1}{G}\left[\sqrt{2}\left(\frac{\mathcal{Q}_{S}\operatorname{csch}^{2}(H)+\left(\mathcal{Q}_{S}-1+\frac{\sqrt{2}\mathcal{M}^{\dagger}}{1+\tfrac{1}{\mathcal{P}e^{2}}}\right)\coth^{2}(H)}{\mathcal{Q}_{S}}\right)\Upsilon\right.\\\
&\hskip
20.00003pt+\left(\frac{\coth(H)\left(\cosh(2H)\left(4\Sigma\mathcal{L}_{1}+\mathcal{L}_{2}\coth(H)\right)+\cosh(3H)\operatorname{csch}(H)\mathcal{L}_{3}+4H\mathcal{L}_{4}-2\Sigma\mathcal{L}_{5}+\coth(H)\mathcal{L}_{6}\right)}{\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(\mathcal{Q}_{S}\right)^{3}\mathcal{Q}_{H}}\right)\Upsilon^{2}\\\
&\left.\hskip 20.00003pt+\mathcal{O}(\Upsilon^{3})\right]\end{aligned}$

(27)

where

$\displaystyle\mathcal{L}_{1}$ $\displaystyle=\tfrac{1}{2}\mathcal{N}_{1}$
$\displaystyle\mathcal{L}_{2}$
$\displaystyle=\mathcal{N}_{2}-2\mathcal{L}_{3}$
$\displaystyle\mathcal{L}_{3}$
$\displaystyle=\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(4+17\mathcal{B}^{\dagger
4}+\mathcal{B}^{\dagger 6}\right)$ $\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{M}^{\dagger}\left(-10\sqrt{2}+24\mathcal{M}^{\dagger}+(10\sqrt{2})/\mathcal{P}e-7\sqrt{2}\mathcal{B}^{\dagger
4}+15\mathcal{M}^{\dagger}\mathcal{B}^{\dagger
4}+\tfrac{20\sqrt{2}}{\mathcal{P}e}\mathcal{B}^{\dagger 4}\right)$
$\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger}\left(2\sqrt{2}\mathcal{B}^{\dagger
4}-12\mathcal{M}^{\dagger}\mathcal{B}^{\dagger
4}-27\sqrt{2}\mathcal{M}^{\dagger 2}+15\mathcal{M}^{\dagger
3}+\tfrac{-34\mathcal{M}^{\dagger}+17\sqrt{2}\mathcal{M}^{\dagger
2}}{\mathcal{P}e}\right)$ $\displaystyle\hskip 10.00002pt+\mathcal{M}^{\dagger
2}\left(8\sqrt{2}\mathcal{M}^{\dagger}+4\mathcal{M}^{\dagger
2}-5\sqrt{2}\mathcal{M}^{\dagger 3}+\mathcal{M}^{\dagger
4}+\tfrac{-4+14\sqrt{2}\mathcal{M}^{\dagger}-24\mathcal{M}^{\dagger
2}+4\sqrt{2}\mathcal{M}^{\dagger 3}}{\mathcal{P}e}\right)$
$\displaystyle\mathcal{L}_{4}$
$\displaystyle=\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(4+10\sqrt{2}\mathcal{B}^{\dagger}+23\mathcal{B}^{\dagger
2}+14\sqrt{2}\mathcal{B}^{\dagger 3}+8\mathcal{B}^{\dagger
4}-\mathcal{B}^{\dagger 6}\right)$ $\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{M}^{\dagger}\left(-9\sqrt{2}+23\mathcal{M}^{\dagger}-28\mathcal{B}^{\dagger}+42\sqrt{2}\mathcal{M}^{\dagger}\mathcal{B}^{\dagger}-18\sqrt{2}\mathcal{B}^{\dagger
2}+48\mathcal{M}^{\dagger}\mathcal{B}^{\dagger 2}\right.$ $\displaystyle\hskip
100.00015pt\left.-8\mathcal{B}^{\dagger 3}+\sqrt{2}\mathcal{B}^{\dagger
4}-15\mathcal{M}^{\dagger}\mathcal{B}^{\dagger
4}+\tfrac{10\sqrt{2}+46\mathcal{B}^{\dagger}+42\sqrt{2}\mathcal{B}^{\dagger
2}+32\mathcal{B}^{\dagger 3}-6\mathcal{B}^{\dagger 5}}{\mathcal{P}e}\right)$
$\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger
2}\left(4+18\sqrt{2}\mathcal{M}^{\dagger}-8\mathcal{M}^{\dagger
2}+24\sqrt{2}\mathcal{B}^{\dagger}+24\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}+36\mathcal{B}^{\dagger
2}\right.$ $\displaystyle\hskip
100.00015pt\left.-6\sqrt{2}\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}+15\mathcal{B}^{\dagger 2}\mathcal{M}^{\dagger
2}-12\mathcal{B}^{\dagger 4}\right.$ $\displaystyle\hskip
100.00015pt\left.+\tfrac{28-14\sqrt{2}\mathcal{M}^{\dagger}+36\sqrt{2}\mathcal{B}^{\dagger}-32\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}+24\mathcal{B}^{\dagger
2}-4\sqrt{2}\mathcal{B}^{\dagger 3}+20\mathcal{B}^{\dagger
3}\mathcal{M}^{\dagger}}{\mathcal{P}e}\right)$ $\displaystyle\hskip
10.00002pt+\mathcal{M}^{\dagger
3}\left(8\sqrt{2}-4\mathcal{M}^{\dagger}+\sqrt{2}\mathcal{M}^{\dagger
2}-\mathcal{M}^{\dagger
3}+16\mathcal{B}^{\dagger}-4\sqrt{2}\mathcal{B}^{\dagger 2}\right.$
$\displaystyle\hskip
100.00015pt\left.+12\mathcal{M}^{\dagger}\mathcal{B}^{\dagger
2}+\tfrac{4\sqrt{2}-8\mathcal{M}^{\dagger}-8\mathcal{B}^{\dagger}+4\sqrt{2}\mathcal{M}^{\dagger}\mathcal{B}^{\dagger}-6\mathcal{M}^{\dagger
2}\mathcal{B}^{\dagger}+8\mathcal{B}^{\dagger 3}}{\mathcal{P}e}\right)$
$\displaystyle\mathcal{L}_{5}$
$\displaystyle=\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(4+4\sqrt{2}\mathcal{B}^{\dagger}-12\mathcal{B}^{\dagger
2}-34\sqrt{2}\mathcal{B}^{\dagger 3}-66\mathcal{B}^{\dagger
4}-32\sqrt{2}\mathcal{B}^{\dagger 5}-14\mathcal{B}^{\dagger 6}\right)$
$\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{M}^{\dagger}\left(-4\sqrt{2}+20\mathcal{B}^{\dagger}+54\sqrt{2}\mathcal{B}^{\dagger
2}+92\mathcal{B}^{\dagger 3}+32\sqrt{2}\mathcal{B}^{\dagger
4}-12\mathcal{M}^{\dagger}-102\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}\right.$
$\displaystyle\hskip 100.00015pt\left.-396\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}-320\sqrt{2}\mathcal{B}^{\dagger
3}\mathcal{M}^{\dagger}-210\mathcal{B}^{\dagger
4}\mathcal{M}^{\dagger}\right.$ $\displaystyle\hskip
100.00015pt\left.+\tfrac{4\sqrt{2}-24\mathcal{B}^{\dagger}-102\sqrt{2}\mathcal{B}^{\dagger
2}-264\mathcal{B}^{\dagger 3}-160\sqrt{2}\mathcal{B}^{\dagger
4}-84\mathcal{B}^{\dagger 5}}{\mathcal{P}e}\right)$ $\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger
2}\left(-48\sqrt{2}\mathcal{B}^{\dagger}-264\mathcal{B}^{\dagger
2}-256\sqrt{2}\mathcal{B}^{\dagger 3}-168\mathcal{B}^{\dagger
4}-54\sqrt{2}\mathcal{M}^{\dagger}-276\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}\right.$
$\displaystyle\hskip 100.00015pt\left.-192\sqrt{2}\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}+66\mathcal{M}^{\dagger
2}+160\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger
2}+210\mathcal{B}^{\dagger 2}\mathcal{M}^{\dagger 2}\right.$
$\displaystyle\hskip
100.00015pt\left.+\tfrac{-20-108\sqrt{2}\mathcal{B}^{\dagger}-276\mathcal{B}^{\dagger
2}-128\sqrt{2}\mathcal{B}^{\dagger
3}+34\sqrt{2}\mathcal{M}^{\dagger}+264\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}+320\sqrt{2}\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}+280\mathcal{B}^{\dagger
3}\mathcal{M}^{\dagger}}{\mathcal{P}e}\right)$ $\displaystyle\hskip
10.00002pt+\mathcal{M}^{\dagger
3}\left(-20\sqrt{2}-184\mathcal{B}^{\dagger}-128\sqrt{2}\mathcal{B}^{\dagger
2}+128\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}+168\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}+32\sqrt{2}\mathcal{M}^{\dagger
2}-14\mathcal{M}^{\dagger 3}\right.$ $\displaystyle\hskip
100.00015pt\left.+\tfrac{-20\sqrt{2}+128\sqrt{2}\mathcal{B}^{\dagger
2}+112\mathcal{B}^{\dagger
3}+92\mathcal{M}^{\dagger}+128\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}-32\sqrt{2}\mathcal{M}^{\dagger
2}-84\mathcal{B}^{\dagger}\mathcal{M}^{\dagger 2}}{\mathcal{P}e}\right)$
$\displaystyle\mathcal{L}_{6}$
$\displaystyle=\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{3}\left(28\mathcal{B}^{\dagger
2}+48\sqrt{2}\mathcal{B}^{\dagger 3}+95\mathcal{B}^{\dagger
4}+32\sqrt{2}\mathcal{B}^{\dagger 5}+15\mathcal{B}^{\dagger 6}\right)$
$\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)^{2}\mathcal{M}^{\dagger}\left(-6\sqrt{2}-48\mathcal{B}^{\dagger}-94\sqrt{2}\mathcal{B}^{\dagger
2}-140\mathcal{B}^{\dagger 3}-57\sqrt{2}\mathcal{B}^{\dagger
4}+36\mathcal{M}^{\dagger}+116\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}\right.$
$\displaystyle\hskip 100.00015pt\left.+432\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}+320\sqrt{2}\mathcal{B}^{\dagger
3}\mathcal{M}^{\dagger}+225\mathcal{B}^{\dagger
4}\mathcal{M}^{\dagger}\right.$ $\displaystyle\hskip
100.00015pt\left.+\tfrac{+6\sqrt{2}+40\mathcal{B}^{\dagger}+130\sqrt{2}\mathcal{B}^{\dagger
2}+300\mathcal{B}^{\dagger 3}+180\sqrt{2}\mathcal{B}^{\dagger
4}+84\mathcal{B}^{\dagger 5}}{\mathcal{P}e}\right)$ $\displaystyle\hskip
10.00002pt+\left(1+\tfrac{1}{\mathcal{P}e^{2}}\right)\mathcal{M}^{\dagger}\left(16\mathcal{B}^{\dagger}+24\sqrt{2}\mathcal{B}^{\dagger
2}+40\mathcal{B}^{\dagger 3}+14\sqrt{2}\mathcal{B}^{\dagger
4}+24\mathcal{M}^{\dagger}+16\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}-192\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}\right.$ $\displaystyle\hskip
100.00015pt\left.-256\sqrt{2}\mathcal{B}^{\dagger
3}\mathcal{M}^{\dagger}-180\mathcal{B}^{\dagger
4}\mathcal{M}^{\dagger}-93\sqrt{2}\mathcal{M}^{\dagger
2}-364\mathcal{B}^{\dagger}\mathcal{M}^{\dagger 2}\right.$
$\displaystyle\hskip 100.00015pt\left.-264\sqrt{2}\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger 2}+81\mathcal{M}^{\dagger
3}+160\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger
3}+210\mathcal{B}^{\dagger 2}\mathcal{M}^{\dagger 3}\right.$
$\displaystyle\hskip
100.00015pt\left.+\tfrac{+8\sqrt{2}+32\mathcal{B}^{\dagger}+28\sqrt{2}\mathcal{B}^{\dagger
2}+24\mathcal{B}^{\dagger
3}-62\mathcal{M}^{\dagger}-160\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}-392\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}-188\sqrt{2}\mathcal{B}^{\dagger
3}\mathcal{M}^{\dagger}}{\mathcal{P}e}\right.$ $\displaystyle\hskip
100.00015pt\left.+\tfrac{+51\sqrt{2}\mathcal{M}^{\dagger
2}+276\mathcal{B}^{\dagger}\mathcal{M}^{\dagger
2}+320\sqrt{2}\mathcal{B}^{\dagger 2}\mathcal{M}^{\dagger
2}+280\mathcal{B}^{\dagger 3}\mathcal{M}^{\dagger 2}}{\mathcal{P}e}\right)$
$\displaystyle\hskip 10.00002pt+\mathcal{M}^{\dagger
2}\left(-16-64\sqrt{2}\mathcal{B}^{\dagger}-72\mathcal{B}^{\dagger
2}+24\sqrt{2}\mathcal{M}^{\dagger}+264\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}+188\sqrt{2}\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}+28\mathcal{M}^{\dagger 2}\right.$ $\displaystyle\hskip
100.00015pt\left.-128\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger
2}-168\mathcal{B}^{\dagger 2}\mathcal{M}^{\dagger
2}-43\sqrt{2}\mathcal{M}^{\dagger 3}+15\mathcal{M}^{\dagger 4}\right.$
$\displaystyle\hskip
100.00015pt\left.+\tfrac{-28-8\sqrt{2}\mathcal{B}^{\dagger}+64\mathcal{B}^{\dagger
2}+24\sqrt{2}\mathcal{B}^{\dagger
3}+66\sqrt{2}\mathcal{M}^{\dagger}+96\mathcal{B}^{\dagger}\mathcal{M}^{\dagger}-128\sqrt{2}\mathcal{B}^{\dagger
2}\mathcal{M}^{\dagger}}{\mathcal{P}e}\right.$ $\displaystyle\hskip
100.00015pt\left.+\tfrac{-112\mathcal{B}^{\dagger
3}\mathcal{M}^{\dagger}-136\mathcal{M}^{\dagger
2}-164\sqrt{2}\mathcal{B}^{\dagger}\mathcal{M}^{\dagger
2}+36\sqrt{2}\mathcal{M}^{\dagger
3}+84\mathcal{B}^{\dagger}\mathcal{M}^{\dagger 3}}{\mathcal{P}e}\right)$

## References

  * [1] Pliny the Elder.  _Historia Naturalis_.
  * Franklin et al. [1774] B. Franklin, W. Brownrigg, and Farish.  Of the Stilling of Waves by means of Oil.  _Philos. Trans. Roy. Soc. Lond._ , 64:445–460, 1774.
  * Aitken [1882] J. Aitken.  On the effect of oil on a stormy sea.  _Proc. Roy. Soc. Edin._ , 12:56–75, 1882.
  * Cini et al. [1983] R. Cini, P. P. Lombardini, and H. Hühnerfuss.  Remote sensing of marine slicks utilizing their influence on wave spectra.  _Int. J. Remote Sens._ , 4(1):101–110, 1983.
  * Brekke and Solberg [2005] C. Brekke and A. H. S. Solberg.  Oil spill detection by satellite remote sensing.  _Remote Sens. Environ._ , 95(1):1–13, 2005.
  * Ghanmi et al. [2015] Helmi Ghanmi, Ali Khenchaf, and Fabrice Comblet.  A new method for reliable detection of polluted sea surfaces.  _SPIE Newsroom_ , pages 2–6, sep 2015.
  * Picard and Davoust [2007] C. Picard and L. Davoust.  Resonance frequencies of meniscus waves as a physical mechanism for a DNA biosensor.  _Langmuir_ , 23(3):1394–1402, 2007.
  * Picard and Davoust [2009] C. Picard and L. Davoust.  A “fluid cantilever” to detect amphiphilic biomolecules.  _Colloid. Surface. A_ , 343(1-3):12–19, 2009\.
  * Lucassen-Reynders and Lucassen [1970] E.H. Lucassen-Reynders and J. Lucassen.  Properties of capillary waves.  _Adv. Colloid Interfac._ , 2(4):347–395, 1970\.
  * Douady [1990] S. Douady.  Experimental study of the Faraday instability.  _J. Fluid Mech._ , 221:383–409, 1990.
  * Jiang et al. [1993] Q. Jiang, Y. C. Chiew, and J. E. Valentini.  The study of surface dilational properties of nonionic surfactant solutions by propagation of electrocapillary waves.  _J. Colloid. Interf. Sci._ , 155(1):8–15, 1993\.
  * Raynal et al. [1999] F. Raynal, S. Kumar, and S. Fauve.  Faraday instability with a polymer solution.  _Eur. Phys. J. B_ , (9):175–178, 1999.
  * Saylor et al. [2000] J. R. Saylor, A. J. Szeri, and G. P. Foulks.  Measurement of surfactant properties using a circular capillary wave field.  _Exp. Fluids_ , 29(6):509–518, 2000.
  * Behroozi et al. [2007] P. Behroozi, K. Cordray, W. Griffin, and F. Behroozi.  The calming effect of oil on water.  _Am. J. Phys._ , 75(5):407–414, 2007.
  * Shao et al. [2018] X. Shao, J. R. Saylor, and J. B. Bostwick.  Extracting the surface tension of soft gels from elastocapillary wave behavior.  _Soft Matter_ , 14(36):7347–7353, 2018.
  * Lau et al. [2020] Y. M. Lau, J. Westerweel, and W. Van De Water.  Using Faraday Waves to Measure Interfacial Tension.  _Langmuir_ , 36(21):5872–5879, 2020.
  * Henderson et al. [1991] D. M. Henderson, K. Larsson, and Y. K. Rao.  A study of wheat storage protein monolayers by Faraday wave damping.  _Langmuir_ , 7(11):2731–2736, 1991.
  * Henderson [1998] D. M. Henderson.  Effects of surfactants on Faraday-wave dynamics.  _J. Fluid Mech._ , 365:89–107, 1998.
  * Wright and Saylor [2003] P. H. Wright and J. R. Saylor.  Patterning of particulate films using Faraday waves.  _Rev. Sci. Instrum._ , 74(9):4063–4070, 2003\.
  * Lange and Hühnerfuss [1984] P. A. Lange and H. Hühnerfuss.  Horizontal Surface Tension Gradients Induced in Monolayers by Gravity Water Wave Action.  _J. Phys. Oceanogr._ , 14(10):1620–1628, 1984\.
  * Reynolds [1880] O. Reynolds.  On the effect of oil in destroying waves on the surface of water.  _Brit. Assoc. Rept. Papers_ , 50:489–490, 1880.
  * Levich [1941] V. Levich.  The damping of waves by surface-active substances I.  _Acta Physicochim. URS._ , 14(3):307–320, 1941\.
  * Dorrestein [1951] R. Dorrestein.  General linearized theory of the effect of surface films on water ripples.  _Proc. K. Ned. Akad. B-Ph._ , 54:250–272, 1951.
  * Case and Parkinson [1957] K. M. Case and W. C. Parkinson.  Damping of surface waves in an incompressible liquid.  _J. Fluid Mech._ , 2(02):172–184, 1957.
  * Goodrich [1961] F. C. Goodrich.  The Mathematical Theory of Capillarity. II.  _Proc. Roy. Soc. A_ , 260(1303):490–502, 1961\.
  * Davies and Vose [1965] J. T. Davies and R. W. Vose.  On the Damping of Capillary Waves by Surface Films.  _Proc. Roy. Soc. A_ , 286(1405):218–234, 1965\.
  * Miles [1967] J. W. Miles.  Surface-wave damping in closed basins.  _Proc. Roy. Soc. A_ , 297(1451):459–475, 1967\.
  * Faraday [1831] M. Faraday.  On the forms and states assumed by Fluids in contact with vibrating elastic surfaces.  _Philos. Trans. Roy. Soc. Lond._ , 121:319–340, 1831.
  * Edwards and Fauve [1994] W. S. Edwards and S. Fauve.  Patterns and quasi-patterns in the Faraday experiment.  _J. Fluid Mech._ , 278:123–148, 1994.
  * Gollub [2006] Jerry P Gollub.  Dynamics of Spatio-Temporal Cellular Structures.  In Innocent Mutabazi, José Eduardo Wesfreid, and Etienne Guyon, editors, _Dynamics of spatio-temporal cellular structures: Henri Benard centenary review_ , volume 207 of _Springer Tracts in Modern Physics_ , pages 213–224. Springer New York, New York, NY, 2006.
  * Ibrahim [2015] R. A. Ibrahim.  Recent advances in physics of fluid parametric sloshing and related problems.  _J. Fluid. Eng. - T. ASME_ , 137(9), 2015.
  * Ballesta and Manneville [2005] P. Ballesta and S. Manneville.  Signature of elasticity in the Faraday instability.  _Phys. Rev. E_ , 71(2):026308, 2005.
  * Benjamin and Ursell [1954] T. B. Benjamin and F. Ursell.  The Stability of the Plane Free Surface of a Liquid in Vertical Periodic Motion.  _Proc. Roy. Soc. A_ , 225(1163):505–515, 1954\.
  * Kumar and Tuckerman [1994] K. Kumar and L. S. Tuckerman.  Parametric instability of the interface between two fluids.  _J. Fluid Mech._ , 279(-1):49–68, 1994.
  * Chen and Viñals [1999] P. Chen and J. Viñals.  Amplitude equation and pattern selection in Faraday waves.  _Phys. Rev. E_ , 60(1):559–570, 1999.
  * Kumar and Matar [2002] S. Kumar and O. K. Matar.  Parametrically driven surface waves in surfactant-covered liquids.  _Proc. Roy. Soc. A_ , 458(2027):2815–2828, 2002\.
  * Kumar and Matar [2004a] S. Kumar and O. K. Matar.  On the Faraday instability in a surfactant-covered liquid.  _Phys. Fluids_ , 16(1):39–46, 2004a.
  * Kumar and Matar [2004b] S. Kumar and O. K. Matar.  Erratum: “On the Faraday instability in a surfactant-covered liquid” [Phys. Fluids 16, 39 (2004)].  _Phys. Fluids_ , 16(8):3239, 2004b.
  * Ubal et al. [2005a] S. Ubal, M. D. Giavedoni, and F. A. Saita.  Elastic effects of an insoluble surfactant on the onset of two-dimensional Faraday waves: a numerical experiment.  _J. Fluid Mech._ , 524:305–329, 2005a.
  * Ubal et al. [2005b] S. Ubal, M. D. Giavedoni, and F. A. Saita.  The formation of Faraday waves on a liquid covered with an insoluble surfactant: influence of the surface equation of state.  _Lat. Am. Appl. Res._ , 35(1):59–66, 2005b.
  * Ubal et al. [2005c] S. Ubal, M. D. Giavedoni, and F. A. Saita.  Influence of surface viscosity on two-dimensional Faraday waves.  _Ind. Eng. Chem. Res._ , 44(4):1090–1099, 2005c.
  * Giavedoni and Ubal [2007] M. D. Giavedoni and S. Ubal.  Onset of Faraday Waves in a Liquid Layer Covered with a Surfactant with Elastic and Viscous Properties.  _Ind. Eng. Chem. Res._ , 46(15):5228–5237, 2007\.
  * Mikishev et al. [2016] A. B. Mikishev, B. A. Friedman, and A. A. Nepomnyashchy.  Generation of transverse waves in a liquid layer with insoluble surfactant subjected to temperature gradient.  _Fluid Dyn. Res._ , 48(6), 2016.
  * Scriven [1960] L.E. Scriven.  Dynamics of a fluid interface Equation of motion for Newtonian surface fluids.  _Chem. Eng. Sci._ , 12(2):98–108, 1960.
  * Waxman [1984] A. M. Waxman.  Dynamics of a couple-stress fluid membrane.  _Stud. Appl. Math._ , 70(1):63–86, 1984.
  * Douady et al. [1989] S. Douady, S. Fauve, and O. Thual.  Oscillatory Phase Modulation of Parametrically Forced Surface Waves.  _Europhys. Lett._ , 10(4):309–315, 1989.
  * Suman and Kumar [2008] B. Suman and S. Kumar.  Surfactant- and elasticity-induced inertialess instabilities in vertically vibrated liquids.  _J. Fluid Mech._ , 610:407–423, 2008.
  * Fuller and Vermant [2012] G. G. Fuller and J. Vermant.  Complex Fluid-Fluid Interfaces: Rheology and Structure.  _Annu. Rev. Chem. Biomol._ , 3(1):519–543, 2012\.
  * Jaensson and Vermant [2018] N. Jaensson and J. Vermant.  Tensiometry and rheology of complex interfaces.  _Curr. Opin. Colloid Interf. Sci._ , 37:136–150, 2018.
  * Samaniuk and Vermant [2014] Joseph R. Samaniuk and Jan Vermant.  Micro and macrorheology at fluid-fluid interfaces.  _Soft matter_ , 10(36):7023–33, 2014.
  * Lemaire and Langevin [1992] C. Lemaire and D. Langevin.  Longitudinal surface waves at liquid interfaces: Measurement of monolayer viscoelasticity.  _Colloid. Surface._ , 65(2-3):101–112, 1992.
  * Buzza et al. [1998] D. M. A. Buzza, J. L. Jones, T. C. B. McLeish, and R. W. Richards.  Theory of surface light scattering from a fluid–fluid interface with adsorbed polymeric surfactants.  _J. Chem. Phys._ , 109(12):5008–5024, 1998.