The Marangoni effect on small-amplitude capillary waves in viscous fluids

We derive a general integrodifferential equation for the transient behaviour of small-amplitude capillary waves on the planar surface of a viscous fluid in the presence of the Marangoni effect. The equation is solved for an insoluble surfactant solution in concentration below the critical micelle concentration (cmc) undergoing convective-diffusive surface transport. The special case of a diffusion-driven surfactant is considered near the the critical damping wavelength. The Marangoni effect is shown to contribute to the overall damping mechanism and a first-order term correction to the critical wavelength with respect to the surfactant concentration difference and the Schmidt number is proposed.


I. INTRODUCTION
On the free surface between any two immiscible phases, external perturbations such as thermally induced motion or vibrations can cause statistical fluctuations [1], or interfacial waves, leading to rough surfaces [1,2]. Long waves with wavelength λ √ σ/(ρg) ≡ l * (where ρ is fluid density, g is gravitational acceleration, and σ is surface tension coefficient) are gravitydriven and known as gravity waves. Short waves with wavelength λ l * are capillary waves where the restoring effects of surface tension dominate gravity. Furthermore, in regions of large viscosity, the energy dissipation is also affected by the distribution of fluid vorticity [3], leading to an increasing viscous attenuation of the capillary wave [4] for decreasing wavelength. The evolution of such waves is critical to the study of film stability, e.g., thin film rupture [5], droplet coalescence [6], polymer films [7], and foam dynamics [8], as well as capillarydriven phenomena, such as the breakup of liquid sheets, jets [9,10], bridges [11], ligaments [10], and curtains [12].
In addition, interfaces found in nature and engineering applications are often contaminated with surface-active substances. These surfactants can significantly alter interfacial flow structures and energy dissipation through the Marangoni effect as a result of the generated gradients in surface tension [13,14]. Hence, an understanding of how the Marangoni effect changes the damping of the interface through the attenuation of capillary waves could extend the classification of surface roughness [2] to a greater number of industrial and biological films and membranes.
The classical treatment of the capillary wave considers the interfacial displacement as an infinite superposition of modes, then assuming each component contributes according to an equipartition theorem [1,15] or approximating the solution as the least-damped mode [4,16]. In either case, while the long-time behavior is sufficiently captured as the waves are damped out and converge towards a single mode, the transient behavior remains difficult to extract because of its dependence on multiple of such modes. To circumvent this difficulty, Prosperetti [3] recast the equation of motion for viscous fluids as an integro-differential equation solved exactly as an initial value problem. The long-time connection of the solution with the normal mode approach was also demonstrated through the complex characteristic equation for the discrete wave spectrum [4] along with the equivalence of irrotational flows with the short-time approximation of the solution. Motivated by this approach, we consider in this paper the capillary waves in viscous fluids with a significant Marangoni effect in the presence of a surfactant solution.
As a prelude to our discussion, we consider the implicit dispersion relation of the capillary wave, derived by linear response theory or linearized hydrodynamics [17,18] under the assumption of weak viscous damping, where ν is the kinematic viscosity, k is wave number, ω is angular frequency of the capillary wave, and is the angular frequency on an interface of an inviscid ideal fluid. For large wave number k 1, ω 0 (σ k 3 /ρ) 1/2 , at which the effects of gravity are negligible. By inspection, Eq. (1) suggests a critical wave number k c , where there is a transition from the propagative and underdamped mode to an overdamped mode [17] in which the effect of viscosity dominates the oscillatory motion of the interface. More recently, Denner [19] proposed a rational parametrization for the dispersion in the underdamped regime, including critical damping, that leads to a self-similar solution of the frequency dispersion of capillary waves in viscous fluids.
The region near the critical wavelength λ c = 2π/k c is particularly interesting due to its proximity to the threshold whereby complex rupture behaviors are exhibited by biological and industrial interfacial systems [5,14,20,21]. Therefore, the evolution of λ c as a result of the Marangoni effect could potentially provide clues to the onset of these phenomena.
For constant surface tension, the wave number k c can be related on dimensional grounds to the quantity δ ν = (ν/ω) 1/2 , which is the characteristic thickness of the penetration of vorticity that is generated at the interface boundary due to the inability of the initially irrotational flow to sustain the nonzero tangential boundary stress condition [22]. More precisely, it has been shown [19,23] that k c ∼ 1/l vc , where is a viscocapillary length scale under the regime Oh = 1, where is the Ohnesorge number. The quantity t σ = ρ l 3 /σ is the capillary time scale proportional to the period of an undamped capillary wave and t μ = ρ l 2 /μ is the viscous time scale associated with the length of time taken for momentum to diffuse through a characteristic length scale l.
In the presence of a variable surface tension, λ c must increase since the motion driven by surface tension gradients is always dissipative in nature due to the induced viscous shear stresses. Consider the stress boundary condition where T is the viscous stress tensor, n the normal vector, and 2H = −∇ s · n the mean curvature, where is the surface gradient operator, M is the Marangoni term, and [·] denotes the jump across the interface of the quantity within.
For systems with surface-active substances, Levich and Krylov [24] suggested M = −μ s ∇ 2 u s , where u s = (I − nn) · u is the surface velocity and μ s is a phenomenological surface viscosity coefficient. This term explicitly accounts for the energy dissipation at the interface due to the irreversible processes involved, and hence a natural consequence is that the Marangoni effect directly contributes to the overall damping of the system. However, the phenomenological and model-specific nature of μ s does not easily yield to general analysis, and so in this paper we shall adopt the classical approach commonly used in the literature [13] of letting M = −∇ s σ , coupled with an equation of state linking surface tension σ and surfactant concentrations (or another field such as temperature). Specializing to the case of a predominantly diffusive surfactant solution, the critical wave number k c is now dimensionally related to the ratio δ ν /δ D , where δ D = (D/ω) 1/2 is the thickness of the mass transfer boundary layer with surfactant diffusivity D, more commonly known as its square the Schmidt number,

II. EQUATION OF MOTION
The dynamics of the fluid of viscosity μ and density ρ in regions of Reynolds number Re 1 are described by the time-dependent Stokes' equation, where u is the two-dimensional fluid velocity field (u,v), p is the pressure, and gj is the gravity term, with j denoting the upward unit vector. The free surface is given by the standing wave where the nonlinear, time-dependent wave amplitude a(t) satisfies the small-amplitude conditions that for wavelength λ, wave number k, angular frequency ω, and phase velocity v p . The linearized tangential and normal stress components, T and T ⊥ , respectively, as well as the kinematic conditions are given on the surface y = 0 by where n (ak sin kx,1) and t (1, − ak sin kx) are the leading-order normal and tangent vectors. The surface tension coefficient σ is modeled by a linear equation of state, where α = |dσ/d (x,y,t)| is the magnitude of surface tension gradient and (x,y,t) is a function which depends on the mode of the Marangoni effect (e.g., thermal, soluto, or electric etc., [14]). In this paper, we shall consider the soluto-Marangoni effect of an interface in the presence of a surfactant solution. Henceforth, = (x,t) is the surface concentration of an insoluble surfactant solution and D = D s denotes the surface diffusivity coefficient. Following [3], we decompose the velocity and pressure into harmonic and viscous corrections, where the harmonic component (u ,p ) assumes the velocity field for a harmonic function φ satisfying the hydrostatic potential problem 053110-2 with known solution [3,4] On the other hand, the incompressible viscous correction component (u ,p ) satisfies the Stokes' problem and the incompressibility condition ∇ · u = 0. Taking the curl of Eq. (24) yields the vorticity equation where ω z is given on the free surface f(x,t) = a cos kx by using the kinematic condition and the insoluble surfactant condition which neglects any surface adsorption processes. Satisfying Eq. (24) with the stream function ψ defined by we can write such that Eq. (16) takes on the biharmonic form Furthermore, writing and assuming that is bounded as y → −∞ gives the Green's function solution of the form where satisfies the second-order equation with the initial condition (y,0) = 0. Writing the surfactant concentration in the wave form then yields the boundary condition The Laplace transform of Eq. (36) then giveŝ wheref is the Laplace transform of f (t).
Integrating Eq. (35), eliminating ∂ /∂t using Eq. (36) and integration by parts yield and the viscous pressure is given by Substituting Eq. (8) into the normal stress condition gives the integro-differential equation of motion in the nondimensional form where is the complementary error function, τ = ωt is the nondimensional time, τ is the integration variable (a nondimensional time discussed in Sec. III), = νk 2 /ω is the nondimensional kinematic viscosity, and the nondimensional surfactant diffusivity and surfactant strength parameters are given by In the absence of the Marangoni effect, i.e., for β = 0, Eq. (43) reduces tö as derived in [3], where˙= d/dτ , * is the convolution operator, and F(τ ) is the auxiliary function 053110-3

III. NONLINEAR VISCOUS DISSIPATION VIA FRACTIONAL INTEGRAL
The viscous equation of motion in Eq. (46) can be rewritten in the form where J(τ ) is the forcing function given by and G(τ ) is a local viscous dissipation density function defined by Recall the definition for the left-handed Riemann-Liouville fractional integral of order [25,26] for t > 0, ∈ C via where f (τ ) is a function in time, is the γ function and is a nonlinear time which satisfies the scaling property such that for the local linear time variables τ = ωt and τ = ωt, we have We can now interpret the effect of viscosity on the system as introducing a forcing term and both a linear and a nonlinear viscous dissipation term. Recasting the irrotationally inviscid surface wave system into the well-known Hamiltonian form [27,28], the governing equation in the fluid region ϒ is written where H is the Hamiltonian defined by and h = h(x,t), = φ(x,h,t) are the canonical coordinate, momentum variables denoting fluid height and interfacial fluid velocity potential, respectively, for x = (x,y). The term 4 ȧ can then be viewed as the linear part the dissipation function Q(h, ) of the system as we equate the rate of loss of energy [28,29], with the rate of dissipation in incompressible fluids under irrotational flow given [22] by where q 2 = (∂φ/∂x i ) 2 for velocity potential φ.
The fractional integral term in Eq. (48) suggests that the nonlinear part of the dissipation function Q(h, ) is the exponential decay of G with decay rate (which is half the rate obtained [4] for the inviscid flow with velocity potential φ = k −1ȧ e ky cos kx). Let then d(τ ) can be interpreted as the actual nonlinear viscous dissipation for which the local values of the dissipation density G(τ )e − (τ −τ ) and the local linear time τ are encoded with respect to the global inhomogeneous time scale given by the function τ ν (τ ; τ, 1 2 ). From Fig. 1, τ ν is approximately linear for small time τ and the exponentially-decaying fractional integral term can be shown [3] to be negligible. As τ increases, the time scale τ ν visibly slows down and deviates from the linear time. This slowdown for large τ could explain why the integral term in Eq. (43) cannot be neglected for τ −1 , since the exponential decay of the local viscous dissipation is progressively retarded for increasing τ .
Incorporating the Marangoni effect, the equation of motion Eq. (43) takes the same form as Eq. (48), where the forcing function J and the local viscous dissipation density function  G are transformed via the map defined by where

IV. SURFACTANT MODELING: MARANGONI CONVECTION
To model the surfactant transport, we consider a mass balance [30] at the interface for an insoluble surfactant of concentration (x,t). The convective-diffusive equation takes the leading-order form Assuming the wave form and the initial condition (x,t = 0) = 0 cmc , where cmc is the critical micelle concentration for surfactant, Eq. (62) reduces to where The resulting equation of motion is given by the simultaneous integro-differential equation where δ = a 0 k, ζ = Dk 2 /ω for surface diffusivity coefficient D. Let Let σ (n) i be the nth-order cyclic polynomial given by Decomposingˆ into partial fractions yieldŝ where −z i are the roots of the polynomial Q(s 1/2 ) and by comparison with Lagrange interpolation, the coefficients c i are given by the expression Taking the inverse Laplace transform of Eq. (70) gives where we have used the fact (shown in Appendix A) that which implies Z(n,0) = 0, where Hence, the amplitude is where ϕ = ϕ(z i ,τ ; ) is given by In the next section, we derive the explicit form of the amplitude solution in the special case near the critical wavelength where convective effects are neglected.

V. THE EXACT DIFFUSIVE-MARANGONI AMPLITUDE SOLUTION
Consider the solution This yields the interfacial vorticity whereχ and F ν (s) is the Laplace transform of the solution of the viscous amplitude equation in Eq. (46), which satisfies the equation Evaluating termwise, the inverse Laplace transform of where z i are the roots of the polynomial P (z) = z 4 + 2 z 2 + 4 3/2 z + 2 + 1.
To simplify Eq. (86), consider the sum then [3] gives the relations Similarly, we have the relations (derivation in the Appendix B) where Therefore, the solution to the full amplitude equation in Eq. (43) is given by where a ν (τ ) = 4a 0 is the solution to the viscous amplitude equation [3] in Eq. (46). In the region near the critical wavelength λ (0) c , we note that Marangoni diffusion is significant compared to Marangoni convection for Sc 10 4 in the regions near the critical wavelength. To examine this in more detail, we define the Marangoni criticality coefficient ς by where ζ (0) denotes the numerically calculated damping ratio of the capillary wave in the absence of Marangoni effect and ζ (1) ,ζ (2) are the numerical damping ratios (1) with diffusional Marangoni effect and (2) with full convectivediffusional Marangoni effect, respectively. At ς = 1 2 , the effects of Marangoni convection contribute equally towards the global damping ratio as the Marangoni diffusion effects. Consequently, in the region where ς > 1 2 the diffusional Marangoni effect contributes more towards the overall damping ratio than the convectional Marangoni effect, and vice versa for ς < 1 2 . Furthermore, the regions where ς > 0.95 are defined as the Marangoni diffusion regime, where Marangoni diffusion dominates compared with the convection effects and is responsible for more than 95% of the total Marangoni effect. Similarly, the region for which ς < 0.05 denotes the Marangoni convection regime, which is not specifically investigated in this study. Under room temperature and pressure (rtp) with parameters density ρ = 10 3 kg m −3 , surface tension σ = 7.2 × 10 −2 Nm −1 , and viscosity μ = 10 −3 Pa s, we plot ς for water with respect to the nondimensional wavelength λ/λ (0) c normalized using the critical damping wavelength λ (0) c in the absence of Marangoni effect, for Schmidt number Sc = 100 in Fig. 1(a) and Sc = 500 in Fig. 1(b) with various Marangoni numbers Ma > 500. We notice that in both cases the diffusional effects as compared to the convection effects remain significant for roughly 2.5 critical wavelengths and dominate near the critical wavelength for Sc ∼ O(10 2 ) and  Ma ∼ O(10 2 ) as shown in Fig. 1(a); it is only when we increase Sc that the diffusional Marangoni effect noticeably lessens near λ (0) c . This dominance of Marangoni diffusion over convection breaks down for Sc O(10 2 ) and for wavelengths λ 2.5λ (0) c , where Marangoni convection rapidly becomes more significant than diffusional effects. In the next section, we shall focus on the diffusion solution in Eq. (95) to analyze how the critical wavelength evolves for increasing Ma within the Marangoni diffusion regime.

VI. RESULTS AND DISCUSSION
Interactions of the viscosity and the Marangoni effect within the fluid yield different consequences on the interface. While both effects act as damping mechanisms to the surface waves [24]; increasing the viscosity retards the rate at which vorticity generated at the boundary enters the bulk [4,22], whereas increasing the concentration of surface-active substance endows the interface with an effective surface shear and dilatational elasticity that works to suppress surface motion [21,24].
We explore this viscosity-Marangoni interaction in the free-oscillation (i.e., u 0 = 0) solution in Eq. (95) by plotting the amplitude function a(τ ) in Fig. 2 for various Marangoni numbers and two viscosity values. Starting from the critically damped region with nondimensional wavelength λ/λ (0) c = 1.5543 in Fig. 2 where d = b/(2 √ mc) is the damping ratio for viscous damping coefficient b, ω 0 = √ c/m is the undamped frequency, and f (τ ) = 0. Equating the spring constant c with the surface tension σ , let m = ρ/k 3 and consider the viscous damping coefficient be defined as b = μL μ , where L μ = √ 2μ/ρω 0 is the viscous damping length, Denner [19] proposed the critical for a single fluid with a free surface, where = 1.0625 is a constant and l vc = μ 2 /ρσ is the viscocapillary length scale. It has been shown [23] recently that this definition of the critical wavelength λ (0) c also holds for capillary waves (under constant surface tension) with a finite amplitude. Switching on the Marangoni effect, we have the linearized external forcing term f (τ ) = −2 βe −ζ τ in Eq. (98), which solves to give where B = (1 − d 2 ) 1/2 , C = 2 βω 0 /(ζ 2 − 2ζ dω 0 + ω 2 0 ), and A 0 = a 0 / cos φ 0 for phase angle φ 0 satisfying −Bω 0 tan φ 0 = (u 0 /a 0 ) + 2 . Note that the ratio 2 /ζ = 2Sc of the exponentials in the linearized solution Eq. (100) determines whether the motion favors the sinusoidal viscous term or the second term which arises from the Marangoni correction. Moreover, the viscous damping factor d deviates from the viscous case due to the presence of the surfactant solution. We investigate this damping factor below. First, consider a range of Sc from orders 10 1 to 10 4 , which is typical for chemical compounds in liquids under room temperature [30]; we plot in Fig. 3 the ratio of the surface tension difference and the initial surface  To derive an analytical scaling for the critical wavelength λ (Sc) c for different values of Sc, we consider that the viscous damping length L μ , which is the depth of penetration of the vorticity generated by a capillary wave [19,29], must be augmented by a suitable diffusion damping length L D ∼ √ D/ω 0 . Dimensional analysis then gives where θ = θ (Sc). Consider the leading-order expansion which provides a first-order Marangoni correction term to the critical wavelength λ (0) c = 2 1/3 πl vc / obtained by Denner [19].

VII. CONCLUSIONS
To conclude, we derived a generalized integro-differential initial value problem for the wave amplitude of surface capillary waves in the presence of the Marangoni effect which is solved exactly for a surfactant solution with concentration much less than the cmc under convective-diffusive surface transport. In particular, we investigated the diffusively dominated region near the critical wavelength λ (0) c and identified a first-order correction of λ (Sc) c from λ (0) c as a function of the Schmidt number Sc and = σ/σ 0 the ratio of surface tension difference and the initial surface tension σ 0 . This first-order correction provides an initial glimpse into the change of fundamental properties due to the Marangoni effect near the critical wavelength λ (0) c .