Instabilities at a sheared interface over a liquid laden with soluble surfactant

The linear stability of a semi-infinite fluid undergoing a shearing motion over a fluid layer that is laden with soluble surfactant and that is bounded below by a plane wall is investigated under conditions of Stokes flow. While it is known that this configuration is unstable in the presence of an insoluble surfactant, it is shown via a linear stability analysis that surfactant solubility has a stabilising effect on the flow. As the solubility increases, large-wavelength perturbations are stabilised first, leaving open the possibility of mid-wave instability for moderate surfactant solubilities, and the flow is fully stabilised when the solubility exceeds a threshold value. The predictions of the linear stability analysis are supported by an energy budget analysis which is also used to determine the key physical effects responsible for the (de)stabilisation. Asymptotic expansions performed for long-wavelength perturbations turn out to be non-uniform in the insoluble surfactant limit. In keeping with the findings for insoluble surfactant obtained by Pozrikidis & Hill (IMA J Appl Math 76:859–875, 2011), the presence of the wall is found to be a crucial factor in the instability.


Introduction
It is well known that an insoluble surfactant can destabilise the interface between two fluids subjected to a shear flow, even if the fluids are stably stratified and have no inertia (which is essential for the development of interfacial instabilities in shear flows devoid of surfactant [1] [2][3][4]). The onset of this surfactant-induced instability depends crucially on the presence of a wall bounding one of the fluids [5]: with no wall and two semi-infinite fluids, the flow is stable. However, if a second wall is also present, as is the case for channel flow, then instability is not always possible. For example, [2] showed that when the interface lies in the middle of the channel, the flow is stable. Generally speaking inertia exacerbates any instability, but may have a stabilising influence under some conditions [6][7][8].
Recently the present authors showed that for the inertialess channel flow problem, the ability of the surfactant to dissolve in one of the fluids can either enhance or suppress interfacial instability for certain fluids and/or surfactant properties [9,10]. Studies of the related problem of film flow down an inclined substrate have found surfactant solubility to have a destabilising effect [11][12][13]. It has been suggested that the destabilisation due to solubility is connected to a phase shift between the interfacial disturbance and the surfactant mass flux onto/off the interface: this allows a redistribution of surfactant at the interface and it thereby mitigates the Marangoni force that is helping to stabilise the flow [13]. A similar explanation was provided for the stabilisation of the flow-induced Marangoni instability found in a vertically falling film loaded with insoluble surfactant in [14].
In this work, we consider a sheared interface separating a semi-infinite fluid and a liquid layer of different viscosity that is loaded with surfactant and that is bounded below by a solid wall. Our main aim is to investigate the influence of surfactant solubility on the linear stability of the flow. The effects of inertia and gravity are ignored. Instability at the interface is induced due to the underlying shear flow and the disparity in fluid viscosities (which results in a velocity jump at the perturbed interface), as well as the Marangoni forces generated because of the presence of surfactant at the interface that leads to local surface tension gradients [15]. Here, the surfactant can be transferred into the liquid layer via a kinetic flux driving the adsorption and desorption of molecules to and from the bulk fluid [16]. A linear stability analysis is carried out and an eigenvalue problem for the complex wave speed is formulated and solved numerically. An asymptotic approximation valid for perturbations of large wavelength is also examined to provide leading-order approximations for the growth rates. Moreover, an energy budget analysis in the manner of [17] is used to tease out the important physical effects driving the instability. A brief discussion of the case of two semi-infinite fluids, which is obtained when the wall is removed, is also presented and supported by analytical and numerical results.
The outline of the rest of the paper is as follows: in Sect. 2, the full non-linear system governing the dynamics of the flow is provided and is later linearised via an energy budget analysis. The disturbances are then written in a normal mode form to obtain a linearised eigenvalue system for the growth rate. A long-wave approximation of the growth rate is also presented. In Sect. 3, growth rate curves are provided as well as data based on energy budget analysis. The stability of the unbounded problem is also discussed. A summary of the main results is given in Sect. 4.

Linear stability analysis
We consider the evolution of an interface between an incompressible liquid of viscosity μ 1 and density ρ, which is laden with surfactant and which is confined below by a wall, and a semi-infinite fluid above of the same density and viscosity μ 2 , as shown in Fig. 1. The bulk surfactant concentration is assumed to not exceed the critical micelle concentration. At equilibrium, the film has a uniform thickness h 0 and a steady shear flow is imposed. The problem is cast in non-dimensional form by scaling lengths with h 0 , surface tension with the clean value γ 0 that prevails in the absence of surfactant, velocities with γ 0 /μ 1 , time with h 0 μ 1 /γ 0 , pressures with γ 0 / h 0 , the interface surfactant concentration by the maximum packing Γ ∞ , the bulk concentration with Γ ∞ / h 0 , and the mass flux by Γ ∞ γ 0 / h 0 μ 1 . The flow in each fluid region is described by equations for the conservation of mass and momentum. In this study the aim is to focus on the effect of surfactant solubility on the stability of the flow; since the two fluids are assumed to have the same density, the effect of gravity is negligible. The Reynolds number in each layer is defined by Re i = ργ 0 h 0 /μ 2 i , where i = 1, 2 is used to label the fluids (see Fig. 1), and is assumed to be small. Consequently the flow is governed by the Stokes momentum equation and the continuity equation, are the pressure and velocity fields in the two fluids, the gradient operator is defined in the cartesian coordinate system shown in Fig. 1 by ∇ = (∂ x , ∂ y ), and m i = 1 + (m − 1)(i − 1), with m the viscosity ratio (defined later in (7)). The evolution of the interface at y = h(x, t) is governed by the kinematic condition Subscripts will be used through this paper to either denote variables in the respective fluid (if variable i = 1, 2 is used) or partial derivatives with respect to the shown variable (where t is time). The no-slip and no-penetration conditions, respectively, require that u 1 = (0, 0) at the wall, and that velocity in the upper fluid u 2 is bounded in the far-field as y → ∞. Continuity of velocity, requiring that u 1 = u 2 , is imposed at the interface y = h(x, t) together with the continuity of normal and tangential stresses where H = 1 + h 2 x and the notation [F i ] 1 2 = F 1 − F 2 is used. The surface tension in conditions (3) is given by the Langmuir equation of state which is γ = 1 + Ma ln (1 − Γ ) [16,18] and is seen to depend on the interfacial surfactant concentration Γ , the evolution of which is described by the convection-diffusion equation [10,19] 1 with u I = u 1 (x, y = h(x, t), t). Here, the adsorption/desorption flux J b is defined by [9] where the bulk concentration C is evaluated immediately below the interface at y = h − . A positive/negative value of J b corresponds to adsorption/desorption of surfactant onto/from the interface. In the film the transport of surfactant molecules is governed by the advection-diffusion equation and associated boundary conditions, The two boundary conditions ensure that there is no flux of surfactant through the wall and that the flux of surfactant onto the interface matches the adsorption/desorption flux J b defined in (5).
The non-dimensional parameters that appear in the equations above are the viscosity ratio, the Marangoni, Biot, and Péclet numbers (surface and bulk), as well as a solubility parameter, respectively, defined by where k a , k d are the adsorption/desorption kinetic rates, respectively, D s,b are the surface and bulk diffusivities, R is the ideal gas constant, and T is the absolute temperature. We note that the case of an insoluble surfactant is recovered in one of the limits B → 0 or R b → ∞.

Energy budgets
The flow is disturbed by introducing perturbations to the steady state of the form where the steady state is denoted by the overbar notation and the perturbations are assumed to be small,f f .
Here, f stands for the various variables involved in the problem, i.e. u i , p i , h, Γ , C. The steady horizontal velocity is linear in each fluid and given bȳ where s is the shear rate at the interface, while the vertical velocities and pressures are constant in both fluids, v 1,2 (y) = 0,p 1,2 (y) = p 0 . The equilibrium state for the surfactant can be found by setting J b = 0 in (5) and is given by uniform concentrations 0 <Γ < 1 and The corresponding equilibrium surface tension isγ = 1 + Ma ln 1 −Γ . We note that normally the interfacial concentrationΓ is prescribed in computations.
To assess the stability properties of the flow and identify the dominant mechanism responsible for interfacial instability, we perform an energy budget analysis. Since the Reynolds number is zero in the present study, the equation for the kinetic energy of the flow does not offer insight into the instability mechanism, but we nevertheless provide it for completeness in Appendix A. Insight can be obtained, however, by examining an equation for the energy of concentration perturbations. To derive this, the linearised convection-diffusion equation for the concentration in the film is multiplied by the perturbationĈ and integrated over the film flow region y ∈ [0, 1]. The entire equation is then integrated in one horizontal wavelength x ∈ [0, ], where = 2π/k and k is the wavenumber, and spatial periodicity is assumed. Invoking the flux conditions at the wall and interface in (6), the energy budget equation takes the form Equation (10b) offers insight into the growth of concentration perturbations and it can be used to identify an instability in Stokes flow [20]. Instability is indicated by a positive value of the ENC term. The diffusion term DIFC is always negative, while the term FLXC is associated with surfactant flux between the bulk and the interface; if the latter effect is sufficiently strong then instability occurs. Following a similar procedure for the corresponding surfactant transport equation at the interface, but now integrating in the x direction only, the energy equation for the interfacial concentration is found to be Instability is present if the term ENG is positive. On the right-hand side in (11b), the term GSH plays a role only in the presence of the background shear, the GU1 term supplies energy due to the perturbed horizontal velocity gradient at the interface, while the term DIFG captures the dampening effect of interfacial diffusion and is always negative. A positive/negative value of FLXG corresponds to a net adsorption/desorption of surfactant onto/from the interface.
The relative sizes of the various terms in (10b) and (11b) can offer insight into the dominant physical mechanisms that provoke instability. Where terms are positive/negative their effect is destabilising/stabilising.

Normal modes
To perform a normal mode analysis it is convenient to introduce the streamfunctions ψ i , defined so that We represent all of the disturbance quantities in the normal mode formf (x, y, t) = f (y)e ik(x−ct) , where k is the real wavenumber to be specified, c is the complex wave speed to be found, and the tilde-decorated variables are eigenfunctions to be determined. The stability of the flow is controlled by the sign of the growth rate λ ≡ kIm(c); in particular, instability occurs if λ > 0.
Substituting the normal mode forms into the governing equations (1)-(6) and linearising, we obtain an eigenvalue problem for the complex wave speed c. The streamfunctions satisfy the ordinary differential equatioñ with the boundary conditions wherec = c −ū 1 (1) is the perturbed wave speed. The jump term in condition (12c) arises due to viscosity stratification. The linearised forms of the normal and tangential interfacial stress balances (3) are The surfactant concentration perturbation in the bulk satisfies and the linearised form of the interfacial surfactant concentration equation (4) is The general solution for the perturbation streamfunctions satisfying Eq. (12a) is for coefficients α j,i , j = 1, 2, 3, 4, i = 1, 2, to be determined. The far-field conditions in (12b) yield α 1,2 = α 2,2 = 0. The general solution for the bulk concentration equation (12f) is where Ai and Bi are the linearly independent solutions of the Airy equation, and b 1,2 are arbitrary constants. Given the exact solutions (13), (14), conditions (12) may be assembled to form the linear system

Long-wave approximation
Useful insight into the behaviour of the complex wave speed can be gained by considering perturbations of large wavelength. Assuming that k 1, we introduce the expansion for the complex wave speed substitute it into the transcendental equation det(M) = 0, and solve at successive orders in k to determine c 0 , c 1 , etc. In doing so we identify two modes which we term the primary modes, for each of which c 0 is real. The two corresponding values of c 1 are both purely imaginary and yield the leading-order approximation to the growth rate λ ≈ k 2 Im(c 1 ).
Besides the two normal modes that are captured by the expansion (16) there is an infinite number of other modes, which we term secondary modes, whose existence may be attributed to the freedom of the bulk surfactant concentration to disperse in the liquid film. These secondary modes satisfy long-wave expansions different to (16); and, notably, Im(c) < 0 at k = 0.
One of the two primary modes is always stable and satisfies Im c The first-order contribution to the second primary mode has imaginary part Im(c where and it is either stable or unstable depending on the parameters. In particular, noting that the denominator in (18b) is positive, the Descartes' rule of signs implies that there is exactly one positive value of R b at which the numerator, viewed as a polynomial in R b , changes sign and therefore the flow changes stability. For sufficiently large R b we find that Im c (2) 1 > 0 so that this mode is unstable in the case of weak solubility, which is qualitatively consistent with established results for insoluble surfactant [3,5]. Since Im c the second primary mode is evidently stable in the strong solubility limit R b → 0. We conclude that surfactant solubility has a stabilising influence on the flow. The fact that neither of the two leading-order expressions for the growth rate in (17) and (18b) depend on the viscosity ratio m is striking but it is in line with similar long-wave calculations for insoluble surfactant (e.g. for channel flow [3]). It is intriguing to note that Im c and hence, curiously, neither of the leading-order growth rates of the two primary modes are consistent with the growth rates presented by Pozrikidis & Hill [5] for insoluble surfactant in the limit B → 0 or R b → ∞. This apparent contradiction may be resolved by noting that the expansion (16) breaks down in either of the limits B → 0 or R b → ∞; instead the correct expansion should proceed in powers of k 1/2 . The situation is reminiscent of the non-uniformity in the long-wave expansion that has been identified in the limits s → 0, Ma → ∞, m → 1, and R b → R * b , where R * b is a certain finite value, for the similar flow in a channel (see [4,9]). We emphasise that the failure to achieve consistency with the results for insoluble surfactant in the limit B → 0 or R b → ∞ is a facet of the current long-wave approach; taking either of these limits in the solution to the full transcendental equation det(M) = 0 stemming from (15) recovers the growth of [5], as is discussed in Appendix B.

Results and discussion
We seek numerical solutions for the complex wave speed that satisfies system (15) by solving the linear eigenvalue problem (12) numerically using a Chebyshev collocation method [9,21]. In practice, it is convenient to solve the problem in a channel of large aspect ratio 0 ≤ y ≤ n for a suitably large value of n and impose no-slip and no-flux conditions at the upper wall. We map each fluid region to the canonical domain Y ∈ [−1, 1] by taking and The interface is therefore located at Y = −1, found by setting y = 1 in each fluid region. In our computations to be presented below we take n = 100, which we found to be sufficiently large for the presence of the upper wall to have a negligible effect on the results. Unless specified otherwise, the results presented in this section will take the following parameter values: s = 1, m = 0.5, Ma = 0.1,Γ = 0.5, Pe s = 10 8 , Pe b = 100 (the surface Péclet number is taken to be large due to the small surface diffusion rate D s which is normally found in practice). The two parameters that control the surfactant solubility will typically assume values in the ranges 0 ≤ B ≤ 1 and 0 < R b < ∞ and will be altered accordingly in order to investigate the impact of the solubility on interfacial stability. We begin by comparing the numerically computed growth rates to the leading-order long-wave approximations presented in Sect. 2.3 for the two primary modes. In Fig. 2 the first few modes calculated using the Chebyshev method are shown with solid lines for two B values and R b = 1, and the long-wave predictions for the largest value of B are shown with symbols (circles for mode 1 in Eq. (17) and squares for mode 2 in Eq. (18b)). The thin lines indicate growth rates for B = 0.01 so  We note that only two modes pass through the origin at k = 0 in Fig. 2, namely the primary modes identified in Sect. 2.3, and it is these that are responsible for any instability. The remaining modes are secondary modes for which Im(c) < 0 at k = 0 and they are stable over the entire wave number range. In fact, the secondary modes were found to be stable at all k for all parameter sets that we examined; in Appendix C, we prove the stability of these modes in the case of B = 0. ' Figure 3 demonstrates the effect of varying the Biot number B on the growth rates. The dominant mode for the insoluble case (B = 0) is unstable over the range 0 ≤ k ≤ k c for k c ≈ 1.4, but is stabilised for any B > 0 at sufficiently small wavenumbers. For weakly soluble surfactant, that is sufficiently small B, mid-wave instability occurs over a window of wavenumbers away from zero (this type of instability has also been found in similar systems, see for example [9,20,22]). For larger values of the Biot number such as B = 0.1, the system is stable over the entire wavenumber range. Figure 4 illustrates how a similar stabilisation occurs by instead decreasing the solubility parameter R b , in which case the maximum value of the growth rate is seen to monotonically decrease for decreasing values of R b . Here the system is unstable for large enough R b and is completely stabilised when R b drops below the threshold value R b ≈ 2.4 and the surfactant solubility is sufficiently strong. Figure 5 shows the interface concentration energy budget (11b) for the case examined in Fig. 3. The eigenfunctions are calculated by solving the linear system (12) using the Chebyshev method and applying the normal mode form of the disturbances into the energy equations (10b) and (11b). All of the terms in (11b) have been normalised with the magnitude of the diffusive term DIFC in (10b) (note that DIFG is very small owing to the large Péclet number Pe s and as such is not suitable for the normalisation). For the case B = 0.01 instability occurs when ENG switches sign at k ≈ 0.19, consistent with Fig. 3; the flow is stable for B = 0.1 (notice also that the terms presented for B = 0.1 are approximately five times smaller in amplitude than those for B = 0.01). The flux term FLXG is everywhere negative in both cases: surfactant is overall desorbed from the interface with a stabilising effect on the flow. Evidently the GSH term is both positive and larger in amplitude than the rest of the terms in the energy equation (11b), indicating that it is the underlying shear that is driving the instability. Note that as expected the term ENC in the bulk energy budget (10b) becomes positive when ENG becomes positive, but this is not shown.

Unbounded flow
The unbounded problem with two semi-infinite fluids separated by an interface was first investigated by [23] in the absence of surfactant. They found that instability only occurs in the presence of inertia and/or unstable density stratification. For Stokes flow with equal-density fluids, [5] showed that the flow is stable in the presence of an insoluble surfactant. Exactly two normal modes exist in this case-a hydrodynamic mode and a Marangoni modeboth of which can be found in a closed form.
The unbounded case of two semi-infinite fluids is effectively attained in the present formulation by considering a perturbation of sufficiently short wavelength. There are then only two length scales in the problem provided by the groupings D s μ 1 /γ 0 and D b μ 1 /γ 0 , the latter associated with surfactant solubility. Thus for insoluble surfactant with negligible surface diffusivity, D s = 0, there is no natural length scale and the growth rates are independent of the wave number. Even with negligible surface diffusivity (in our calculations we have assumed Pe s = 10 8 1), the growth rates for the soluble unbounded problem do depend on the wavenumber underscoring the importance of the length scale D b μ 1 /γ 0 . We note that the hydrodynamic mode is identical to that obtained in the case of a clean flow.
For two semi-infinite fluids, the coefficient b 2 in (14) must vanish to ensure a bounded concentration in the far-field. Looking at the form of ζ in (14) suggests the scaling c ∼ k −1/3 in order to retain solubility effects in the analysis. We can then expand the wave speed as c = c 0 k −1/3 + · · · in the transcendental equation, and we find that to satisfy the equation at leading order we need It is known that the derivative of the Airy function Ai (z) only has zeros on the negative real axis and there are an infinite number of them [24]. Consider this infinite set of roots z n = −ξ n < 0, n = 1, 2, . . ., with ξ n > 0. The leading-order wave speed c 0 can then be written as c 0 = −Aξ n e iπ/6 , where coefficient A = s 2/3 Pe is a positive constant. Hence any solution for c 0 lies on a line in the third quadrant, at an angle of 30 • with the negative real axis. The leading-order growth rate is therefore stable since Im(c 0 ) = −Aξ n /2 < 0 for all values of n.
We have provided above some analytical evidence that the growth rates of an unbounded flow are stable for small values of k. For arbitrary values of k, the stability of the growth rates is supported by the numerical results presented earlier that correspond to a semi-infinite flow. In particular, it has been demonstrated that all growth rates are stable for sufficiently short waves (this is also true for Fig. 2 at larger k not shown in the plot).

Conclusions
We have investigated the linear stability of an interface between two sheared viscous fluids above a flat wall, in the case when the thinner fluid layer is filled with a soluble surfactant. In particular, we have found that a soluble surfactant has a stabilising impact on the flow. We performed a linear stability analysis for Stokes flow and analysed  [5], in the case of insoluble surfactant the flow is always unstable to some perturbation of finite wavelength. We have found that the flow is completely stabilised when the surfactant is sufficiently soluble. Long-wave disturbances are stabilised by a relatively low level of solubility, but mid-wave instability is supported for weakly soluble surfactant; that is, perturbations of moderate wavelengths may be unstable.
The results discussed above were obtained numerically for disturbances of arbitrary wavelength by solving a linear eigenvalue problem for the wave speed of the perturbations. We have also examined the linear system analytically in the case of large-wavelength disturbances using an asymptotic approach. Somewhat unexpectedly, the leading-order forms of the dominant growth rates do not reduce to those found for an insoluble surfactant when the appropriate limits are applied, a discrepancy that we attribute to a non-uniformity in the relevant asymptotic expansion.
We have also carried out an energy budget analysis which provided insight into which of the terms in the surfactant transport equation are key for (de)stabilisation when this is supported. Specifically, we have demonstrated that a term proportional to the shear rate at the interface dominates over all other terms and hence, in keeping with previous studies with insoluble surfactant, it is shear that provides the key component to drive the instability.
It is interesting to contrast the results of this study with those in [12,13] for a film flow down an inclined substrate, where it was found that solubility has a destabilising effect. Therefore different flow fields disrupt the uniform surfactant distribution at the interface in different ways, leading to Marangoni stresses that may be stabilising (gravity-driven film flow with inertia) or destabilising (shear-driven film flow, as here). However, in both cases the effect of increasing solubility is to mitigate the respective trend, i.e. the stabilising/destabilising effect becomes weaker.
Finally, we have examined the importance of the wall in determining the stability of the flow. In common with the findings of Pozrikidis & Hill [5] for insoluble surfactant, we have found that the presence of the wall is necessary for interfacial instability when the surfactant is soluble. This outcome is suggested by our numerical and analytical results (based on asymptotic analysis for long-wave perturbations) for an unbounded flow. The physical role played by the wall in destabilising the flow remains an open question and is the subject of our current investigations.
cf. [25][26][27]. The first term on the left-hand side, KIN, is the spatially averaged rate of change of the kinetic energy of the disturbances (and is proportional to the growth rate, [17]), while the second term, REY, is the energy transfer from the basic flow to the disturbances via the Reynolds stresses [25]. Both of these terms are proportional to the Reynolds number in the film, which is minuscule due to the small fluid thickness; while these terms are left in the energy budget for completeness, in practice they are both zero for the problem considered here and so Eq. (24b) cannot be used to determine whether the flow is stable or unstable. The terms on the right-hand side are explained as follows: term DIS is viscous dissipation and is always negative; term TEN represents the rate of work done against the interfacial tension in deforming the interface [17,26] (a term involving the hydrostatic pressure gradient would also be added here had the effect of density stratification been also considered); term TAN represents the transfer of energy from the basic flow to the perturbed flow at the interface [23]-this term is non-zero since the viscosities are unequal, m = 1; and finally, term MAR is due to the Marangoni stresses at the interface.

Appendix B: Insoluble surfactant limits
Limit B → 0 In this limit the mass flux J b tends to zero and the adsorption/desorption process to/from the interface is suspended. The two normal modes obtained via linear stability analysis are the same as those coming from the insoluble surfactant case, as expected. Now there is an infinite number of additional modes that are always stable; the origin of these modes is from the bulk equation, which is now decoupled from the rest of the system and is solved with homogeneous boundary conditionsC (0) =C (1) = 0.
Limit R b → ∞ In this case, the linearised interfacial flux becomes The second boundary condition in (12g) then becomesC (1) + Pe b B R b (1 −Γ )C(1) = 0, hence in the limit R b → ∞ the condition simplifies toC(1) = 0. Consequently, none of the flux terms will remain in the linearised transport equation for the surfactant at the interface (since the only surviving term, seen in Eq. (25), isC(1) = 0 as found above), and hence the interfacial equation is reduced to that for insoluble surfactant. Finally, we note that the bulk equation is again decoupled from the rest of the system but is now solved with boundary conditions C (0) =C(1) = 0, in which case all modes found numerically are seen to not pass through the origin k = 0.

Appendix C: Normal modes of the bulk equation for B = 0
In the case of insoluble surfactant obtained by setting B = 0, the linearised bulk equation (12f) is solved with homogeneous boundary conditionsC (0) =C (1) = 0. We will show that all modes originating from the bulk concentration system are stable. Multiplying (12f) by the complex conjugateC * , integrating over the vertical fluid domain [0, 1], applying integration by parts, and using the boundary conditions, result in therefore it is clear that Re(c) > 0 and Im(c) < 0, and hence the growth rates have the property λ = k Im(c) ≤ 0 for all k ≥ 0.