Capillary waves with surface viscosity

Experiments over the last 50 years have suggested a tentative correlation between the surface (shear) viscosity and the stability of a foam or emulsion. We examine this link theoretically using small-amplitude capillary waves in the presence of a surfactant solution of dilute concentrations where the associated Marangoni and surface viscosity effects are modelled via the Boussinesq-Scriven formulation. The resulting integro-differential initial value problem is solved analytically and surface viscosity is found to contribute an overall damping effect on the amplitude of the capillary wave with varying degrees depending on the lengthscale of the system. Numerically, we find the critical damping wavelength to increase for increasing surface concentration but the rate of increase remains different for both the surface viscosity and the Marangoni effect.


Introduction
Capillary waves on a viscous fluid interface have recently been observed (Aarts et al. 2004) to induce the spontaneous breakup of a thin liquid film and controls the inherent stochastic process of the sub-micron rupture event. Unlike gravity waves, these capillary waves have a short wavelength where the restoring force of surface tension dominates over the influence of gravity and can be found in the study of small-lengthscale interfacial phenomena; for instance, thin liquid films (Scheludko 1967) and droplet coalescence (Blanchette & Bigioni 2006). It is apparent that variations in surface tension can have dramatic knock-on effects on the dynamics of the capillary waves, with applications found both in surface chemistry (Edwards et al. 1991) and interfacial fluid dynamics (Levich & Krylov 1969).
Surface-active materials, or surfactants, often lead to the formation of foams and emulsions by lowering the surface tension of a liquid interface (Batchelor et al. 2003;Levich & Krylov 1969;Edwards et al. 1991). Gradients of surfactant concentration (and therefore the surface tension coefficient) caused by dilatational deformations induce the Marangoni stress, which acts to oppose the changes in surface area and slows down the drainage and rupture processes of a thin liquid film. Moreover, the two-dimensional surfactant monolayer displays the rheological response whereby shearing deformations can introduce an extra surface shear viscosity. In addition, a source of surface dilatational viscosity can result from the inherent compressibility of the two-dimensional surfactant monolayer (Zell et al. 2014); in direct contrast with the incompressible Newtonian bulk fluid which can be characterised entirely by a single viscosity parameter µ. Furthermore, the dissipative nature of the surfactant adsorption-desorption kinetic process can also contribute towards the effective surface dilatational viscosity (Lucassen & Hansen 1966). With multiple sources of surface viscosities, we henceforth denote the effective surface dilatational and shear viscosity by µ d and µ s , respectively. Finally, we note that the magnitude of µ d and µ s need not be comparable (Djabbarah & Wasan 1982), since they are each responsible for different physical processes. The physical manifestation of surface viscosity and its measurement remain controversial and subtle; for decades of literature cannot agree on measurements of µ s and µ d , the chief difficulty lies with the fact that not only are surface viscosity and Marangoni effects intimately intwined (Levich & Krylov 1969;Scheid et al. 2010;Langevin 2014), experiments give the total characteristics for both the surface and bulk phases simultaneously and it is not trivial to extract the surface information a-priori of the establishment of a particular surface model.
For insoluble surface-active (surfactants) solutions, the intrinsic surface shear viscosity is clearly defined. However, for soluble surfactant solutions, in particular sodium dodecyl sulfate (SDS), the presence of a 3-dimensional sublayer adjacent to the surface alters the rates of surface deformation (Stevenson 2005), which may explain the numerous inconsistencies in the reported literature of the magnitudes of surface shear viscosity of surfactants. Some progress has been made recently, namely by the experimental work of Zell et al. (2014), in which the use of microbutton surface rheometry appears to yield relatively unambiguous measurements of the surface shear viscosity µ s of SDS. They report an upper bound of µ s ∼ O(10 −8 Nsm −1 ), which suggests that surface shear viscosity need not be the dominant surface phenomena and that Marangoni effects and surface dilational viscosity may also be in effect.
In insoluble surfactant solutions, the surface shear viscosity is often much higher than O(10 −8 Nsm −1 ), in particular, in the case of 1-eicosanol, it is found (Zell et al. 2014;Gavranovic et al. 2006) to be at least 10 3 -10 4 times higher than that of soluble SDS solutions. Moreover, recent numerical (Gounley et al. 2016) and experimental studies conclude that surface viscosity effects in insoluble surfactants can give rise to noticeable behaviours on the resulting dynamics, which cannot otherwise be fully understood if we considered the Marangoni effect alone (Ponce-Torres et al. 2017). In this paper, we shall investigate both the effects of Marangoni and surface viscosity in insoluble surfactant solutions with a particular focus on the dynamics of very thin films with capillary waves close to critical damping. For such a thin film geometry of high wavenumber, we may consider a two-dimensional flow structure as well as a low Reynolds number under the Stokes' limit. Foreshadowed by the previous numerical work (Gounley et al. 2016;Ponce-Torres et al. 2017), we anticipate a similar importance of the Marangoni and surface viscosity effects on the capillary wave in the two-dimensional thin film case under the Stokes' limit.
In §2-4, we extend the previous work of Prosperetti (1976) and Shen et al. (2017) to incorporate both surface shear and dilatational viscosity to the leading-order, as described by the Boussinesq-Scriven model, into the dynamics of small-amplitude capillary waves. We delineate the effects of the convective-diffusive Marangoni stresses with surface viscosity effects in §5. In §6, we obtained an analytical form of the critical damping wavelength for the clean case considering only the bulk fluid viscosity. In §7, we outline a numerical method to calculate the damping ratio of a general higher-ordered system and construct a minimal pole matrix to encode the information of the poles of the system with significant residue. Under this approach, we identify the transition point of the wave from an underdamped to overdamped state of a general system and obtain numerically the correction to this critical wavelength by surface viscosity and Marangoni effects. The article is concluded in §8.

Boussinesq-Scriven surface viscosity
Under the Boussinesq-Scriven model of surface viscosity (Scriven 1960;Aris 1963;Slattery et al. 2007), the surface stress boundary conditions at the interface between two Newtonian fluids can be written as where T is the viscous stress tensor, ∇ s = P · ∇ is the surface gradient operator for the projection tensor P = I − nn T with normal vector n, [·] denotes the jump in magnitude across the interface and σ s is the surface viscous stress tensor defined by and is the surface rate of deformation tensor. The divergence of σ s may be written (Scriven 1960) in the form are the mean and Gaussian curvatures of a surface, respectively, and σ is the surface tension coefficient. Neglecting higher-order terms, the leading-order surface stress boundary condition in the context of small-amplitude capillary waves takes the reduced form [n · T] = ∇ sσ + 2Hσn, (2.7) whereσ is the surface tension augmented with the leading-order surface viscosity contribution given byσ Using the equation of motion derived in §3, the leading-order surface viscosity effect on the small-amplitude capillary wave can naturally be characterised (Lopez & Hirsa 1998) by the non-dimensional Boussinesq number where Bq d = µ d k/µ and Bq s = µ s k/µ are the Boussinesq dilatational and shear numbers, respectively, for dynamic viscosity µ and wavenumber k. More explicitly, surface viscosity can be modelled to be proportional to the surfactant concentration (Ponce-Torres et al. 2017). In the case of detergents, experimental work by Brown et al. (1953) suggests the bi-partisan action of the special solute pairs present in the detergent; where the primary constituent provides a large reservoir of surface-active material while the secondary constituent, lesser in amount, forms surface films of high viscosity. However, in the leading-order dynamics of the Boussinesq-Scriven formulation, surface viscosity is shown in §3 to not depend explicitly on the surfactant concentration and enters only implicitly via the surface tension coefficient.
The other non-dimensional numbers of the system which arise naturally in the equation of motion are the viscosity (ǫ), surfactant diffusivity (ς) and surfactant strength (β) parameters given by where D s denotes the coefficient of surface diffusivity, ν = µ/ρ is the kinematic viscosity for fluid density ρ, α = |dσ/dΓ | is the gradient of surface tension coefficient, ω is frequency of the capillary wave and Γ 0 is the initial surfactant concentration, which is assumed to be much less than the critical micelle concentration (cmc). In this system, these parameters act as the effective Reynolds, Schmidt and Marangoni numbers, respectively.

Equations of motion
The dynamics of an incompressible fluid of viscosity µ and density ρ in regions of Reynolds number Re = U λ/ν ≪ 1 satisfies the Stokes equation where u = (u, v) is the two-dimensional fluid velocity field, p is the pressure and F = −gj is the external (gravitational) force, with j denoting the upward unit vector in the ydirection, and g is the gravitational acceleration. The small-amplitude capillary wave is given at the free surface F by the standing wave where a(t) is the non-linear, time-dependent wave amplitude which satisfies the smallamplitude conditions that a ≪ λ = 2π/k and da/dt ≪ v c = ω/k, where λ is the wavelength and v c is the phase velocity. For vanishing Gaussian curvature in a two-dimensional space, the leading-order tangential and normal stress components, T and T ⊥ and the kinematic condition are given by respectively. The leading-order normal and tangent vectors are Similar to the small-amplitude condition, we consider a small departure from the equilibrium surface tension and let the coefficient of surface tension σ, to be defined via a linear equation of state where σ 0 is the initial surface tension coefficient, Γ (x, t) is the (dilute) concentration of a surfactant solution where adsorptive-desorptive processes are neglected. Using the wave-form Γ (x, t) − Γ 0 =Γ (t) cos kx, the governing equation for surfactant concentration along a two-dimensional deforming surface (Stone 1990) is given bỹ to the leading order, ω z (x, y, t) = Ω(y, t) sin kx the z-component of the vorticity, * is the convolution operator and F (t) is the auxiliary function The velocity and the pressure can be decomposed into inviscid and viscous parts, i.e.
with well-known solutions (Lamb 1932) dt 2 e ky cos kx , (3.14) where u ′ = ∇φ. The viscous component (u ′′ , p ′′ ) satisfies the Stokes problem which is solved by introducing the streamfunction ψ defined by u ′′ = (ψ y , −ψ x ). Taking the curl of Eq. (3.15) gives the bi-harmonic equation Writing ψ = Ψ (y, t) sin kx, the Stokes problem yields the solution where we have the viscous pressure correction p ′′ (x, y, t) = µΩ(0, t)e ky cos kx and the boundary vorticity Henceforth, using non-dimensional variables τ = ωt, ǫ = νk 2 /ω andΩ = Ω(0, t)/ω, the boundary vorticity becomes the integral equatioñ where we have for A = a/a 0 is the dimensionless amplitude, δ = a 0 k and˙= d/dτ denotes the nondimensional temporal derivative. Substituting the pressure and the velocity into Eq. (3.5) and Eq. (3.10) gives the simultaneous equation Equations (3.23) and (3.24) provides us with a dynamic equation system for the amplitude and the surfactant concentration, the solution of which we outline in the next section.
4. Solution of the simultaneous integro-differential equation where −z i are the roots of the polynomial Q(s 1/2 ). In the absence of the Marangoni effect, the surface viscosity case is given bŷ (4.9) By comparison with Lagrange polynomial interpolation, we have is the n-th order cyclic polynomial given by It follows that by comparing Eq. (4.7) and Eq. (4.11), the coefficients c i are it follows (see Appendix A) that the condition deg Q − deg P = 2 (4.14) implies Z(n, 0) = 0, where deg X is the degree of the polynomial X. Taking the inverse Laplace transform of Eq. (4.8) gives Finally, the non-dimensional amplitude is given by (4.17)

Surface viscosity effects on the wave amplitude
As shown in §4, the leading-order surface viscosity effects on the dynamics of smallamplitude capillary waves are characterised by the parameter B and similarly, the Marangoni effect can be reduced to the non-dimensional variables ς and β given in §2.
In what follows, we use water under room temperature and pressure (rtp), i.e. at 25 • C, with density ρ = 10 3 kgm −3 , surface tension σ = 7.2 × 10 −2 Nm −1 and viscosity µ = 8.9 × 10 −4 Pa s, as a test system with surfactant Schmidt number Sc = ν/D s = 10 4 . We define λ (0) c to be the wavelength for which the capillary wave undergoes critical damping, henceforth known as the critical wavelength. The superscript denotes the clean case which we understand as a system without the addition of surface-active material, i.e. β = B = 0. We will look at this critical wavelength in more detail in §7, but here we note a harmonic oscillator approximation of λ In figure 1, we compare the effect of surface viscosity with that of the equivalent bulk viscosity ν, and the Marangoni effect with that of simple reductions in the surface tension coefficient σ = σ 0 −αΓ 0 . Here we use the phrase equivalent to denote the quantity of bulk viscosity and simple reductions in surface tension which results in an identical effect on the overall wave amplitude due to surface viscosity and the Marangoni effect, respectively, under the limit of either a large or small wavelength. For wavelengths near critical damping, in figure 1(a), increasing B exhibits a relatively large difference to equivalent increases in ν, as compared to the case for λ = 950λ dynamic Marangoni effect vastly overshadows the static changes in the surface tension coefficient. Henceforth, we can approximate the Marangoni effect on the wave amplitude near the critical wavelength with the equivalent reduction in σ.
In figure 2, we see in more detail the mechanisms of surface viscosity and the Marangoni effect at altering wave amplitudes for different wavelengths. Both effects admit relatively self-similar solutions and that increasing either the surface viscosity or the surfactant concentration increases damping as expected. One of the explicit differences is that the Marangoni effect reduces surface tension and, thus, lowers the frequency ω, while surface viscosity leaves ω unchanged. This is due to surface viscosity being dependent on the surfactant concentration only to the linear order, and does not feature in the leading-order nonlinear amplitude equation in Eq.(3.23) explicitly, as noted previously. The consequence on the amplitude of this ω-lowering is a horizontal drift of the waveform for systems with the Marangoni effect, as evident from figures 2(d) and 2(f), as opposed to relatively centred waveforms in the cases of surface viscosity in figures 2(c) and 2(e).
We also observe that the surface viscosity effect weakens for large wavelengths λ ≫ λ c but is very potent for small wavelengths λ ∼ λ c . A particular consequence of this potency is its ability to alter the onset behaviour in interfacial phenomena, a number of which occurs near the region of the critical wavelength. The stochastic nature of many of the interfacial instabilities (Aarts et al. 2004) are often kickstarted by the small-amplitude local disturbances on the interface. Hence a small change in surface material, and thus the wave damping, can have a significant effect in starting or delaying the initialisation process of more complex phenomena. Furthermore, it would be of interest to obtain the modifications to the critical wavelength upon the addition of a small amount surface active material. However, we need to consider the definition of the critical wavelength for a higher-order system as it is not as readily defined as in a second-order system.

Capillary wave dispersion and the critical wavelength
To quantify the changes due to the presence of the Marangoni and surface viscosity effects to the critical wavelength, we must first obtain the critical wavelength in the clean case. Following Lamb (1932), the general dispersion relation for an interface with both the Marangoni effect and surface viscosity can be found (derivation in Appendix B) to take the form and W 0 (Z; ǫ) = Z 4 + 2Z 2 − 4Z + 1 + 1 ǫ 2 . (6.3) Specialising to the clean case, Eq. (6.1) reduces to as derived from linearised hydrodynamics (Levich 1962;Lamb 1932). The dispersion relation admits solution of the form where the polynomial f(ǫ) is given by For ǫ ≪ ǫ ⋆ , the wave frequency can be written as (6.8) and the damping coefficient can be extracted as Im(ω ′ ) ∼ 2ǫ, where ǫ ⋆ ∈ R + is the transition value defined by the (largest positive) root of the polynomial f, i.e.
This agrees with the asymptotic approximations by Levich (1962) which suggests that increasing viscosity would decrease damping for ǫ ≫ ǫ ⋆ . Using Eq. (6.18), the analytical critical wavelength (the value of λ with associated damping ratio ζ = 1) of the capillary wave is given by For water under rtp., we have λ ⋆ c . = 40.1894nm. In comparison, a harmonic oscillator approximation (Denner 2016) in Eq. (5.1) gives the result λ (0) c . = 40.9838nm. Consider the relative error between the harmonic oscillator and the analytical critical wavelengths we observe the system is largely second-order in the neighbourhood of the critical wavelength, as the harmonic oscillator value of λ c is within 2% of the analytical value from the wave dispersion.

Damping ratio for a generalised system
For systems of a higher order, the damping ratio ζ is not naturally defined and we usually inspect the root-locus diagram in order to decompose the system into a sum of first-and second-ordered systems to provide an estimate calculation. Here we consider a numerical method to obtain an equivalent damping ratio, whereby ζ 1 when the area of the amplitude below the settling value vanishes for all time almost everywhere. The critical wavelength is then the supremum of the set of wavelengths such that the above property holds. We express this as where H(x) is the Heaviside step function. For underdamped waves, even for second-order systems, the logarithmic decrement or fractional overshoot methods tend to break down or become less accurate near regions of critical damping. Hence, to determine the damping ratios in the neighbourhood of ζ ≈ 1, we adapt the area method in Eq. (7.1). Consider that the area under the t-axis is given by the function where Λ(ζ, t) satisfies the normalised harmonic oscillator equation The generalised (numerical) damping ratio for ζ 1 can then be obtained by the inverse operation where X is the area under the t-axis of a generalised system. This numerical method agrees well with logarithmic decrement and fractional overshoot schemes in the relevant underdamped regimes. In cases where the higher order system can readily be approximated locally by a second-order system, we note that its dominant poles have a larger residue and time constant t c = 1/(ζ n ω n ) relative to other poles, where ζ n and ω n are the damping ratio and frequency associated with each pole. In cases that are not clear cut, i.e. where all the poles are closer together with t c and residues of a similar magnitude, the numerical definition of the damping ratio above only provides an estimate of the true damping ratio and we need to examine the poles in more detail.
To encode such information into a convenient form, we construct the minimal pole matrix {ζ(λ)}. To decompose the system, we first consider the minimal realisation of the transfer function. For a general system, let Θ 1 = {q i ∈ C : Q(−q i ) = 0} and Θ 2 = {p j ∈ C : P(−p j ) = 0} be the set of poles and zeros of the transfer function tf(s) ≡ P(s)/Q(s), which can be written in the form, Applying the pole-zero cancellation procedure, we obtain the minimal realisation of the transfer function, henceforth known as the minimal transfer function mtf(s), defined by where the polynomials P m (s) and Q m (s) satisfy deg(Q m + P m ) deg(P + Q) (7.8) and Θ 1,m (ϑ) ⊆ Θ 1 is the set of poles of Q with significant residues (tolerance of order ϑ) Θ 1,m (ϑ) = q i,m ∈ C : Q(−q i,m ) = 0, |res(q j )| max qi∈Θ1 |res(q i )| = O(ϑ) . (7.9) Returning to the construction of the minimal pole matrix, we let the first column of the matrix illustrate the order of the poles of the minimal transfer function in dots form; the second column considers the relative magnitudes of their time constant t c ; the third column compares their relative residues; the fourth column gives the damping ratios ζ n associated with each pole. Furthermore, to the right of the line separator, we provide an estimated equivalent second-order damping ratio of the entire system using the area numerical method described previously. We say that a system is second-order dominant if one set of complex conjugate poles dominate the other poles (i.e. having the largest t c and residue). In diagrammatic form for ϑ = 1, we have from which we observe that the system is second-order dominant since the set of complex conjugate poles with associate damping ratio 0.45 dominates (in the sense of residue) the other and, thus, we can deduce that the true damping ratio of the system is close to the approximate second-order value.
We take this analysis to the region near critical damping, i.e. for ζ → 1 + and ζ → 1 − . In the clean case, we take the analytical result λ ⋆ c = 2πl vc /ǫ ⋆2 to be the definition of the critical (damping) wavelength. The relevant minimal pole matrices take the form (7.13) We can see that this transition from λ ⋆+ c → λ ⋆− c for the clean case boasts a transformation of the complex conjugate into two separate first-order poles. Or in diagrammatic form, we have (7.14) i.e. a 2 2 to 1 2 transition. Extending to contaminated systems, we summarise the results of the minimal pole matrices of the system at the critical transition in table 1, where the notation 1 a 2 b denotes a system with a first-order and b second-order poles with significant relative residue. We note that the effect of surface viscosity is to introduce an extra first-order pole (with unit damping ratio) to the system, while the Marangoni effect does not change the pole composition for the underdamped region before the critical damping transition. Moreover, we observe at this critical damping transition that the system enters overdamped regime if its minimal pole matrix is of the 1 2 -type irrespective of its type in the underdamped regime prior to the transition. Henceforth, we shall define a generalised higher-order (capillary) wave to be in overall overdamped motion if its minimal pole matrix is of the 1 2 -type.
Using the definition of the overdamped regime for a generalised higher-order system, we determine numerically the corrections to the critical wavelength for increasing surface viscosity (through B), surfactant concentration (through β) and bulk viscosity through θ, where ν = (1 + θ)ν 0 (7.15) in figure 3. We note that the curves denoting surface viscosity and the Marangoni effect intersect near B ≃ 0.65 and that while the Marangoni curve is roughly exponential, the surface viscosity curve is the sum of two exponential functions. Moreover, we observe that a small amount of surface viscosity present in the system has an amplified effect on the system and that a 7-fold increase in critical wavelength for B = 1 results in very different dynamics and mechanisms as the sub-100nm brings forward the possibility of long-range molecular interactions as well as the hydrodynamics. Also of consideration is the proximity of the critical wavelength to the wavelengths of the visible spectrum of light which allows thin film behaviours to be captured by light scattering methods.
Hence the presence of surfactants could determine whether or not we would be within such a range to allow interferometry techniques.
Comparisons of the range under with experimental and computational results on surface viscosity can be made using values reported previously in the literature. Experimentally, Kanner & Glass (1969) summarises in table 2 the surface viscosities for both surfactant and polymeric films, where, for a dilute amount of sodium lauryl sulfate and polydimethylsiloxane in particular, the surface viscosity corresponds to a lower bound of B = O(1) for k = O(10 6 ). A similar correspondence can be found with the upper bound surface shear viscosity of O(10 −8 Nsm −1 ) found in (Zell et al. 2014) for soluble surfactants. We note that this measurement does not include surface dilatational viscosity and so corresponds to the case (Lucassen & Hansen 1966) where diffusional transport between the surface and the bulk is neglected, and assumes that the bulk viscosity and density are constant right up to the interface. More recently, Gounley et al. (2016) characterises the influence of both shear and dilatational surface viscosity on droplets in shear flow in the range B = O(10 0 ) to O(10 1 ) for a range of capillary numbers.
Beyond the cmc value, the Marangoni effect should in principle have no overall contribution to the capillary wave; the dotted curve in figure 3 would end abruptly at the cmc value. The combination of surface viscosity together with the Marangoni effect is however not straightforward; as in previous experimental studies (Brown et al. 1953;Kanner & Glass 1969), surface viscosity also appears to alter the ability of the Marangoni effect to lower surface tension. It would therefore be fruitful in a future contribution to investigate this surfactant interference mechanism through a more systematic experimental and theoretical study. In particular, a numerical approach similar to that of Sinclair et al. (2018) could include the usage of a nonlinear equation of state for the surface tension coefficient σ. The effect of the deviation from the linear equation of state on the amplitude of the capillary wave would aid the analysis near the cmc value of the surfactant solution.

Conclusion
In this work the surface viscosity effect has been incorporated into the integrodifferential initial value problem describing the wave dynamics of small-amplitude capillary waves via the Boussinesq-Scriven surface model. We have shown that, particularly at lengthscales close to the critical damping wavelength, a very small amount of surface viscosity can dramatically increase the critical wavelength of the capillary waves, in contrast with the Marangoni effect which becomes prominent at larger wavelengths. In view of the important role that capillary waves play in inducing the rupture process of thin films (Aarts et al. 2004), we anticipate the various interfacial phenomena controlling the wave dynamics at the very minute lengthscale to contribute towards the understanding of the stability of foams with non-trivial surface viscosity. In particular, the correction of the critical wavelength due to surface viscosity and Marangoni effects, which we summarised in figure 3 using numerical methods, is bound to alter the onset of fluid instabilities for very thin liquid films. It is also useful towards the optimisation of additives to achieve the desired increases in the critical wavelength. Finally, we expect the concept of a critical damping wavelength and its correction by surface material to be useful to a further number of general interfacial phenomena, such as the onset of thin film quasi-elastic wrinkling and Faraday-like instabilities in the same lengthscale.
where P(s, m) is a polynomial of order m in s, is a polynomial of order n > m in s with distinct roots q i and 1 i1<i2<···<i k n q i1 q i2 · · · q i k = (−1) k ς n−k .
(A 4) Rewritingf(s) using a partial fraction decomposition, we havê and taking an inverse Laplace transform gives where q i are roots of the polynomial Q(s, n) and Z(n, j) =