Stability of liquid films covered by a carpet of self-propelled surfactant particles
aa r X i v : . [ c ond - m a t . s o f t ] A ug Stability of liquid films covered by a carpet of self-propelled surfactant particles
Andrey Pototsky, Uwe Thiele, ∗ and Holger Stark † Department of Mathematics, Faculty of Science Engineering and Technology,Swinburne University of Technology, Hawthorn, Victoria, 3122, Australia Institut f¨ur Theoretische Physik, Westf¨alische Wilhelms-Universit¨at M¨unster,Wilhelm Klemm Str. 9, D-48149 M¨unster, Germany Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin, Hardenbergstrasse 36, 10623, Berlin, Germany
We consider a carpet of self-propelled particles at the liquid-gas interface of a liquid film on a solidsubstrate. The particles excert an excess pressure on the interface and also move along the interfacewhile the swimming direction performs rotational diffusion. We study the intricate influence of theseself-propelled insoluble surfactants on the stability of the film surface and show that depending onthe strength of in-surface rotational diffusion and the absolute value of the in-surface swimmingvelocity several characteristic instability modes can occur. In particular, rotational diffusion caneither stabilize the film or induce instabilities of different character.
PACS numbers: 68.15.+e, 05.40.-a, 05.65.+b
The understanding of the physical principles of the mo-tion of self-propelled particles in viscous fluids [1–6], ei-ther in the bulk or at interfaces is of primary importancefor an increasing number of applications in microfluidicsand medicine [7–9]. A particularly interesting emergingapplication of such swimmers are biocoatings formed us-ing a suspension of living cells that is deposited ontoa solid substrate before the solvent is removed, e.g., byevaporation. This technique is used to fabricate bacterialcarpets consisting of living bacteria with rotating flagellathat are attached head down to a polymer layer [10]. Thecreated homogeneous monolayer of living cells is seen asa prototype of a novel biomaterial with remarkable ap-plications, e.g., as artificial skin, self-cleaning coating,or biosensor [11–15]. Beside applications in biotechnol-ogy, free liquid-gas interfaces loaded with motile bacteriaoccur naturally, for example, at the sea surface [16]. Mi-croswimmers at the interface of a thin liquid film alsoshow interesting collective phenomena since even in a di-lute suspension they interact with each other through thesurface flow field initiated by gradients in the surface ten-sion and curvatures in the height profile. In the following,we present an analysis of the stability of such thin filmsand demonstrate the subtle influence of in-plane swim-ming velocity and rotational diffusivity.The proximity of swimmers to the liquid-gas inter-face inevitably modifies the local surface tension, whichdepends on the swimmer concentration similar to pas-sive surfactant molecules and (nano-)particles [17–20].A gradient in the surface tension due to a non-uniformconcentration generates fluid flow at the surface (solu-tocapillary Marangoni effect), a phenomenon well stud-ied for passive surfactants. The impact of self-propelledsurfactants, i.e., surfactants that are capable to moveautonomously, on the dynamics of liquid-gas interfaces(free surfaces), has been studied only in one special case.Namely, Ref. [21] investigates a monolayer of insolubleswimmers that are adsorbed at the free surface of a liquid film and exclusively swim into the direction perpendic-ular to the surface. That is, they are “head up”, at alltimes and their motion along the liquid-gas interface isas for passive particles. As a result, the swimmers gen-erate an excess normal pressure and fluctuations in theirdensity may “bulge” the interface locally. The inducedfluid flow moves additional swimmers towards the bulgeand increases the local excess pressure further. The inter-face becomes unstable if the combined stabilizing effectof translational diffusion and Marangoni flow towards re-gions of smaller swimmer concentrations is too weak [21].In this letter we show that self-propelled motion ofthe swimmer parallel to the liquid-gas interface (in-planemotion) together with rotational diffusion has a profoundand non-trivial effect on the stability of the film. Depend-ing on the strength of rotational diffusion and swimmingvelocity, the in-plane motion can stabilize a flat film orinduce film instabilities of different character.Consider a liquid film on a smooth homogeneous solidsubstrate with free liquid-gas interface. The liquid-gas in-terface is loaded with self-propelled particles, each char-acterized by a unit vector p that gives the instantaneousdirection of swimming with swimming velocity v . Formean film thicknesses below several hundred microm-eters, the deformations of the liquid-gas interface arelong-wave and can be described in the long-wave or lu-brication approximation [22]. The position of a par-ticle moving along the liquid-gas interface is given bythe three-dimensional (3d) vector [ r ( t ) , h ( x, y, t )], where r ( t ) = [ x ( t ) , y ( t )] is a two-dimensional (2d) position vec-tor and h ( x, y, t ) is the local film height. In long-waveapproximation, the projections of the orientation vec-tor p orthogonal and parallel to the liquid-gas interface,are approximated by its z -component, p ⊥ ≈ p z , and by p k ≈ ( p x , p y ), respectively. Thus, the overdamped equa- perturbation wave vector -0.3-0.2-0.100.1 M a x i m a l g r o w t h r a t e V=00.51.01.21.5 5.0 (a) (b)
FIG. 1: (Color online) (a) Liquid-gas interface of a liquidfilm, loaded with self-propelled surfactant particles. Theirswimming orientations are indicated by the unit vectors p ,shown by red arrows. (b) Influence of the in-line motion ofself-propelled particles on film stability. Largest growth rateof a perturbation around a flat film plotted versus wave num-ber for different values of the in-plane self-propulsion velocity V as indicated. Other parameters are D = 1 and d = 0. Solid(dashed) lines correspond to oscillatory (steady) dynamics. tions of motion for the i -th particle become˙ r i = v ( p k ) i + U i + χ i ( t ) , ˙ p i = (cid:20) η i + 12 Ω i (cid:21) × p i , (1)where χ i ( t ) with h χ i ( t ) ⊗ χ k ( t ′ ) i = 2 M k B T δ ik δ ( t − t ′ ) and η i with h η i ( t ) ⊗ η k ( t ′ ) i = 2 D r δ ik δ ( t − t ′ ) repre-sent translational (2d) and rotational (3d) noise, respec-tively [23]. They both influence the dynamics of the self-propelled surfactants at the liquid-gas interface. M is thetranslational mobility and D r the rotational diffusivity ofthe swimmers (see note [24]), U i = [ U x ( x i , y i ) , U y ( x i , y i )]is the surface velocity of the fluid, and Ω i = ∇ × U i =Ω i e z with Ω = ∂ x U y − ∂ y U x denotes the vorticity at point( x i , y i ). For ξ i = η i = 0, Eqs. (1) reduce to the case stud-ied in Ref. [25], where a swimmer moves in a prescribedtime-independent Poiseuille flow field U . In our case,the flow field is initiated by variations in the liquid-filmheight and by the Marangoni effect.So far, Eqs. (1) describe the dynamics of the full ori-entation vector p at the interface. To reduce the dimen-sionality of the problem, we proceed to decouple the dy-namics of its in-plane component p k from the dynamicsof the vertical orientation p ⊥ . In spherical coordinates,we have p ⊥ ≈ p z = cos θ and p k ≈ ( p x , p y ) = sin θ q with q = (cos φ, sin φ ). Assuming that the characteris-tic rotational diffusion time of the swimmers τ ∼ D − r is much smaller than the characteristic relaxation timeof the film thickness fluctuations (a plausible assumptionfor the quasi-stationary lubrication approximation), onemay safely assume that the distribution of the verticalorientation of swimmers p z adjusts “instantaneously” tosome stationary distribution P ( p z ). This allows us toaverage Eqs. (1) over the angle θ to obtain˙ r i = v q i + U i + χ i ( t ) , ˙ φ i = 12 | Ω i | + ξ i ( t ) , (2) where v = v h p − p z i , with h p − p z i = R p − p z P ( p z ) dp z representing the average abso-lute value of p k , and ξ ( t ) is a Gaussian white noisewith h ξ i ( t ) ξ k ( t ′ ) i = 2 D r δ ik δ ( t − t ′ ). Note that Eqs. (2)describe the 2d dynamics of a self-propelled particlewith propulsion speed v moving in a fluid with velocityfield U = ( U x , U y ).The Smoluchowski equation for the surface particledensity ρ ( r , φ, t ) corresponding to Eqs. (2) is a continu-ity equation with translational and orientational currentsthat, as usual, contain drift and diffusional contributions: ∂ρ∂t + ∇ · ( vρ q + U ρ ) + Ω2 ∂ρ∂φ − D r ∂ ρ∂φ − M k B T ∆ ρ = 0 . (3)Here, ∇ and ∆ denote, respectively, the nabla andLaplace operator in positional coordinates. Eq. (3) iscoupled to the thin film equation via the solutocapillaryMarangoni effect. Following Refs. [21, 22], the equationfor the local film thickness h ( x, y, t ) is given by ∂h∂t + ∇ · (cid:18) h µ ∇ [Σ ∆ h + α h ρ i ] (cid:19) + ∇ · (cid:18) h µ ∇ Σ (cid:19) = 0 , (4)where µ is the dynamic viscosity, Σ( x, y, t ) is theconcentration-dependent surface tension, Σ is the con-stant reference surface tension without surfactants (seenote [26]), and α h ρ i describes the additional pressure ex-erted by the swimming particles onto the liquid-gas inter-face as in Ref. [21]. Here, h ρ i ( x, y, t ) = R π ρ ( x, y, φ, t ) dφ denotes the local particle density averaged over allswimming directions and α = h p z i v /M with h p z i = R P ( p z ) p z dp z is the negative of the force the interfaceexerts on the swimming particles in order to stop themin vertical direction [27]. Finally, the fluid velocity fieldat the free surface reads [22] U = hµ ∇ Σ + h µ ∇ (Σ ∆ h + α h ρ i ) . (5)To close our set of equations, we link surface tensionΣ to particle density ρ . As any passive surfactant, self-propelled particles at the interface modify the local sur-face tension Σ. Here, we assume that the concentrationis low, i.e., the surfactant particles are in a 2d gaseousstate implying that Σ depends linearly on the local par-ticle concentration h ρ i (cf. Ref. [19]):Σ = Σ − Γ h ρ i . (6)Typically, Γ >
0, i.e., the surface tension decreases withincreasing h ρ i and the Marangoni flow is directed towardslower particle concentration.Equations (3) to (6) form a closed system of nonlin-ear integro-differential equations for the two scalar fields h ( x, y, t ) and ρ ( x, y, φ, t ). It is important to note thatfor v = 0 and α = 0, the density is independent of theangle φ , i.e., ρ = ρ ( x, y, t ) and Eqs. (3) to (6) reduce tothe two usual coupled equations for film height h ( x, y, t )and surfactant concentration ρ ( x, y, t ) for passive insolu-ble surfactant [22, 28].To non-dimensionalize, we use h and h p Σ / Γ ρ as the vertical and horizontal length scale, respec-tively, µh Σ / (Γ ρ ) as the time scale and the direction-averaged density of swimmers in the homogeneous state h ρ i = R π ρ / (2 π ) dφ = ρ as the density scale.We introduce the dimensionless in-plane self-propulsionvelocity V = vµ Σ / / (Γ ρ ) / , rotational diffusivity D = D r h µ Σ / (Γ ρ ) [29], surface diffusivity d = k B T M µ/ ( h ρ Γ), and excess pressure β = αh / Γ. Thedimensionless equations of motion derived from Eqs. (3),(4), and (5) are summarized in the Supplement.The trivial homogeneous stationary state correspondsto a flat film covered uniformly by particles, i.e., ρ =1 / (2 π ) and h = 1. Any small amplitude perturbation ofthe trivial uniform state can be represented as δh ( r , t ) = Z ˆ h ( k ) e γ ( k ) t e I kr d k , (7) δρ ( r , φ, t ) = lim N →∞ π N X n = − N e Inφ Z W n ( k ) e γ ( k ) t e I kr d k , with small amplitudes ˆ h ( k ) and W n ( k ), the wave vector k = ( k x , k y ), and the growth rate γ ( k ). For any fixednumber N of the Fourier modes, the growth rate γ ( k ) ofthe fastest growing perturbation is found by solving theeigenvalue problem obtained by linearizing the dimen-sionless Eqs. (3), (4), and (5), as outlined in the Supple-ment. For the parameter values used here, the resultshave converged for N = 10.In what follows, we set d = 0 and focus on the moststriking effects of the in-plane motion and in-plane rota-tional diffusion of the swimmers on film stability. Thechoice of d = 0 is motivated by estimating the ratio V / ( Dd ) = v / ( D r M k B T ) ≫ R = 1 µ m and v = 1 µ m / s.As the eigenvalues γ ( k ) do only depend on the absolutevalue of k but not on its direction (spatial isotropy), wecompute the growth rate Re ( γ ) as a function of k .As a reference we set the excess pressure parameter β = 3 to be slightly above the critical value of β c = 2,where the long-wave steady instability sets in for V = 0.This case corresponds to particles that only swim perpen-dicularly to the free surface as studied in Ref. [21]. Thecorresponding dispersion relation Re ( γ ) vs. k is plottedas heavy solid line in Fig. 1(b) and indicates that the ex-cess pressure due to the upwards swimming surfactantsdrives the film unstable. Fig. 1(b) also shows how thefilm is stabilized by the in-plane self-propelled particlemotion, i.e., when increasing V from zero (here, at fixedrotational diffusion D = 1). We understand this stabi-lization qualitatively since active motion along the inter-face acts like the stabilizing translational diffusion with V /2D V k -0.2-0.100.10.2 R e ( γ ) k -0.2-0.100.1 R e ( γ ) k D h m i n , m a x (a) (b)(c) (d)
456 7 8 stable stable
FIG. 2: (Color online) (a) Stability diagram in the plane( d eff = V / D, V ) at fixed β = 3 and d = 0. The inset showsthe time evolution of the minimum and maximum of the localfilm thickness in the Smooth Particle Dynamics simulationsfor V = 0 (dashed lines) and V = 5 (solid lines) at D = 1.(b) and (c) The dispersion relations (maximal Re ( γ ( k ))) for V = 1 and different D as indicated by numbers near eachcurve: 1: D = 0 . D = 0 .
1, 3: D = 0 .
3, 4: D = 0 .
5, 5: D = 0 .
7, 6: D = 1, 7: D = 2, and 8: D = 10. Solid (dashed)lines correspond to oscillatory (steady) perturbations. Theheavy solid line represents D → ∞ and coincides with thecurve for V = 0 in Fig. 1(b). (d) Stability diagram in theplane ( D, k ). The solid line represents the level line of thedispersion relation at Re( γ ) = 0 .
03 to illustrate its bimodalcharacter. an effective diffusivity d eff = V / D on time scales largerthan the orientational correlation time D − [30–32]. Wefurther note that at low V the dispersion relation hastwo maxima that correspond to a steady (growing dropsand holes) and an oscillatory instability (traveling waves)mode that dominate at V = 0 and V = 0 .
5, respectively.The onset of the steady instability occurs at V ≈ V, d eff )-parameter plane [Fig. 2(a)] provides quantitative insightinto the stabilization of the liquid film. It reveals an inter-mediate range of d eff = V / D where the film is stable.A small effective diffusivity d eff (large D ) cannot stabi-lize the film. The film becomes stable when d eff exceedsa threshold value, which exactly agrees with the stabi-lizing translational diffusivity obtained in Ref. [21] forswimmers with purely head-up orientation (red dashedline). Indeed, by expanding the full density ρ ( x, y, φ, t )into angular moments and deriving dynamic equationsfor the moments from Eq. (3), one can formulate a Smolu-chowski equation for h ρ i ( x, y, t ) = R π ρ ( x, y, φ, t ) dφ ontime scales larger than D − and on large length scales,where the active motion only contributes to an effectivetranslational diffusion constant d + d eff [33, 34]. So, oneobtains the density equation employed in Ref. [28] forpassive surfactants and in Ref. [21] for swimmers withpurely head-up orientation.The selected dispersion relations for the growth ratein Fig. 2(b) and (c) at V = 1 reveal a small wave lengthinstability for D ≈ D ≈ .
6, which cor-responds to the unstable region in (a) for large d eff . Astability analysis of our dynamic equations for D = 0at a fixed delta-peaked distribution of the swimming di-rection reveals a growing density and height modulationwave that travels with the swimming speed V . This im-plies that dense swimmer regions do not disperse andrender the film unstable. With increasing D the swim-ming direction starts to diffuse what suppresses the trav-eling waves more and more (depending on k ) until at D ≈ . k, D )-plane in Fig. 2(d) illustrates the qualita-tively different character of the instabilities at low andlarge D . Above the stable range at intermediate D , thefilm becomes again unstable when the effective transla-tional diffusivity V / D becomes too small. The onsetoccurs at k = 0 as a steady instability, however, withfurther increasing D the dispersion relation develops abimodal character as indicated in Fig. 2(d) by the shownlevel line. Note, finally, that at large D [ V /D → L × L square boxwith periodic boundary conditions, coupled to the equa-tions of motion [Eqs. (2)] for n = 500 individual pointswimmers. Their discrete spatial distribution is trans-lated into a smooth particle density function ρ ( x, y, t )employing the method of the Smooth Particle Dynamics[35] (for details see the Supplement).Fixing the remaining parameters as in Fig.1(b), wevary the self-propulsion velocity V . From the linear sta-bility analysis we expect that the flat film is linearly un-stable for V = 0 and linearly stable for V = 5. Startingfrom randomly located and oriented particles on the flatfilm surface, we observe a time evolution that is noisy inboth cases due to the coupling of the continuum equationto the discrete particle dynamics. This can be appreci-ated in the inset of Fig. 2(a) showing the time evolutionof the minimum and maximum of the local film thick-ness for V = 0 (dashed lines) and V = 5 (solid lines).Corresponding movies are available in the Supplement.One discerns a clear difference between the two cases:For V = 0, the swimmer-induced instability of the flatfilm is apparent. The amplitude of the surface deflections first grows before at later times it varies (in a potentiallychaotic way) about a maximum that is by about a factor2 larger than the mean film thickness h = 1. This resem-bles a regime of interacting nonlinear ’traveling waves’ aspredicted in Ref. [21]. However, for swimmers with in-plane motility (here, V = 5), one finds that after a smallinitial growth the deformation amplitude soon saturatesand then fluctuates about a (small) finite amplitude thatis by about a factor 2 smaller than h = 1. This indicatesthat the instability is strongly suppressed. The remain-ing fluctuations of the film surface are a consequence ofthe hybrid calculations where the discrete stochastic dy-namics of the finite number of swimmers is coupled tothe continuum model for the evolution of the film height.In conclusion, we have shown that rotational diffusioncan play a distinctive role in the motion of self-propelledsurfactants on liquid films. Depending on its relativestrength, it can stabilize or destabilize the film and mayalso change the qualitative character of the instabilityand, in consequence, the nonlinear behaviour. Morespecifically, our analysis shows the flat film can be desta-bilized according to two different scenario: (i) throughthe steady long wave instability at V / D = 0 . D andthe in-plane velocity V , reported here. As an example,we estimate the dimensionless V and D for self-propelledJanus particles [36] in a water film at room temperature:thus we take µ = 10 − kg m − s − , Σ = 10 − N m − ,Γ = 10 kg m s − mol − (from [21]), v = 10 − ms − , R = 10 − m (from [36]). In addition, we estimate the(largest possible) average density ρ ≈ ( πR ) − ≈ − mol m − and take h = 10 − m. Taking into account thatfor a spherical Brownian particle D r = 3 k B T M/ R and M = (6 πµR ) − we obtain V ∼ and V /D ∼ − .These estimates show that self-propelled Janus particlesin a 100 µ m thick water film approximately fall ontothe vertical dashed line in Fig. 2(a). Consequently, byfine tuning the radius R of the particles ( V ∝ R and V /D ∝ R ) or the self-propulsion velocity v , one caninduce or suppress the long-wave steady instability of thefilm.A.P. thanks the research training group GRK 1558funded by DFG for financial support. ∗ Electronic address: [email protected];URL: † Electronic address: [email protected] [1] E. Lauga and T. R. Powers, Rep. Prog. Phys. , 096601(2009).[2] D. Saintillan and M. J. Shelley, C. R. Physique , 497(2013).[3] I. Aronson, ibid , 518 (2013).[4] R. Kapral, J. Chem. Phys. , 020901 (2013).[5] A. Z¨ottl and H. Stark, Phys. Rev. Lett. , 118101(2014).[6] M. Hennes, K. Wolff, and H. Stark, Phys. Rev. Lett. ,238104 (2014).[7] S. K. Chang, V. N. Paunov, D. N. Petsev, and O. D.Velev, Nature Materials , 235 (2007).[8] M. S. Baker, V. Yadav, S. Ayusman, and S. T. Phillips,Angew. Chem. Int. Ed. , 10295 (2013).[9] J. D. Yadav, V. Freedman, M. Grinstaff, and S. Ayus-man, Angew. Chem. Int. Ed. , 10997 (2013).[10] N. Darnton, L. Turner, K. Breuer, and H. C. Berg, Bio-physical J. , 1863 (2004).[11] K. Efimenko, J. Finlay, M. Callow, J. Callow, and G. J.,ACS Appl. Mater Interfaces , 1031 (2009).[12] Q. Xu, C. Barrios, T. Cutright, and B. Newby, Environ.Sci. Pollut. Res. , 278 (2005).[13] R. F. Bennett, Case Study of an Environmental Contam-inant (Cambridge University Press, United Kingdom,1996), p. 21.[14] K. D. Weiss, Prog. Polym. Sci. , 203 (1997).[15] M. Fidaleo, S. Charaniya, C. Solheid, U. Diel,M. Laudon, and M. C. Scriven, L. E. ans Flickinger,Biotechnol. Bioeng. , 446 (2006).[16] M. I. Gladushev, Biophysics of the Surface Microlayer ofAquadic Ecosystems (IWA Publishing, London, 2002).[17] M. J. Rosen and J. T. Kunjappu,
Surfactants and Inter-facial Phenomena (John Wiley & Sons., Hoboken, NewJersey, 2012).[18] B. Binks, Curr. Opin. Colloid Interface Sci. , 21 (2002).[19] U. Thiele, A. J. Archer, and M. Plapp, Phys. Fluids ,102107 (2012).[20] V. Garbin, J. C. Crocker, and K. J. Stebe, Langmuir ,1663 (2012).[21] S. Alonso and A. Mikhailov, Phys. Rev. E , 061906(2009).[22] A. Oron, S. Davis, and S. Bankoff, Rev. Mod. Phys. ,931 (1997).[23] Here ⊗ stands for the tensor product between two vec- tors, is the unit 2 × is the unit 3 × , 218104(2012).[26] The capillary pressure term Σ ∆ h is a direct consequenceof the long-wave approximation [22]. However, to recoverthe gradient dynamics structure discussed in Ref. [19]from Eqs. (3) and (4) in the ’passive’ limit v, α →
0, theterm ∇ · (Σ ∇ h ) would be used. The difference betweenthe two forms is negligible as it is only O ( ε ), where ε isthe smallness parameter of the long-wave approximation.[27] Note that α depends on the mobility of the particle M ,which in its turn is linked to the the particle radius viathe Stokes’ law, i.e. M = (6 πµR ) − for spherical particlesof radius R .[28] L. W. Schwartz, D. E. Weidner, and R. R. Eley, Langmuir , 3690 (1995).[29] Note that for the dimensionless thermal rotational dif-fusion D = k B T σ / (8 π Γ ρ ) × h /R , where R is theradius of the spherical swimmer.[30] J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough,R. Vafabakhsh, and R. Golestanian, Phys. Rev. Lett. ,048102 (2007).[31] J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet,Phys. Rev. Lett. , 088304 (2010).[32] M. Enculescu and H. Stark, Phys. Rev. Lett. , 058301(2011).[33] R. Golestanian, Phys. Rev. Lett. , 038303 (2012).[34] O. Pohl and H. Stark, Phys. Rev. Lett. , 238303(2014).[35] W. G. Hoover, Advanced Series in Nonlinear Dynam-ics Volume 25. Smooth Particle Applied Mechanics. TheState of the Art. (World Scientific Publishing Co. Pte.Ltd., Singapore, 2006).[36] I. Buttinoni, F. K¨ummel, G. Volpe, and C. Bechinger, J.Phys.: Cond. Mat. , 284129 (2012).[37] A. Pototsky and H. Stark, Europhys. Lett.98