Rotational motion of a droplet induced by interfacial tension ∗ Ken H. Nagai Department of Physics, Graduate School of Science, The University of Tokyo, Tokyo, 113-0033, Japan 3 1 Fumi Takabatake and Masatoshi Ichikawa 0 2 Department of Physics, Kyoto University, Kyoto, 606-8502, Japan n a J Yutaka Sumino 0 1 Department of Physics, Faculty of Education, ] Aichi University of Education, Aichi 448-8542, Japan t f o s . Hiroyuki Kitahata t a m Department of Physics, Graduate School of Science, - d Chiba University, Chiba 263-8522, Japan and n o PRESTO, Japan Science and Technology Agency, Saitama, 332-0012, Japan c [ 1 Natsuhiko Yoshinaga v 7 WPI-AIMR, Tohoku University, Sendai 980-8577, Japan 4 0 (Dated: January 11, 2013) 2 . 1 Abstract 0 3 SpontaneousrotationofadropletinducedbytheMarangoniflowisanalyzedinatwo-dimensional 1 : v system. The droplet with the small particle which supplies a surfactant at the interface is consid- i X ered. We calculated flow field around the droplet using Stokes equation and found that advective r a nonlinearity breaks symmetryfor rotation. Theoretical calculation indicates that thedropletspon- taneously rotates whentheradiusofthedropletisan appropriatesize. Thetheoretical resultswere validated through comparison with the experiments. PACS numbers: 47.55.D-, 47.63.mf, 68.03.Cd ∗ [email protected] 1 I. INTRODUCTION Recently, self-propelled motion of both biological objects and nonbiological objects are extensively studied to find general aspects in self-propelled systems [1]. Collective behaviors such as nonequilibrium order-disorder transition of direction of motion are candidates for such general characteristics [2]. One of the best experimental systems to investigate gen- eral aspects is self-propelled objects driven by interfacial tension [3, 4]. Actually, several results on collective behavior of the objects driven by interfacial tension have been reported recently [5]. This is because this kind of systems usually consists of a small number of ele- ments and allows us to perform experiments under broad parameter space. Despite of this advantages, however, mechanical analysis for above droplet motion is still lacking. The above-described droplet motion driven by interfacial tension is observed even under isotropic conditions. Although a droplet itself has no asymmetry, symmetry of interfacial tension is broken through the advective nonlinearity. Actually, Toyota et al. reported spontaneous translational motion of a circular droplet [4]. Recently, theoretical analysis for such translational motion has been reported and the nonlinear effect which breaks the symmetry has been made clear [6]. Not only spontaneous translational motion but also spontaneous rotation has also been demonstrated by Takabatake et al. experimentally [7]. In this system, a droplet rotates even though the system is mirror symmetric. Although this system was analyzed using Langevin equation in the previous study, mechanical analysis is lacking and the physical mechanism of symmetry breaking is not clear. In this article, we analyzed the spontaneous rotation of a droplet using an advection-diffusion equation coupled with Stokes equation. We found that the droplet can rotate when the droplet has an appropriate radius. To validate our model, theoretical results were compared with the experimental ones. II. MODEL EQUATION We consider a circular undeformable oil droplet with a radius, R, in an aqueous phase in a two-dimensional space, as shown in Fig. 1(a). A small solid particle is attached at the interface and a surfactant that reduces interfacial tension is supplied to both the inside and the interface of the droplet from the particle. In the droplet, the surfactant is decomposed. 2 (a) (b-1) (b-2) aqueous phase R q oil FIG. 1. (Color Online) (a) Schematic diagram of the model in this article. (b) The flow field when the interface concentration field is c(θ). The images are stream line around the droplet in the frame where the particle and the center of mass of the droplet are fixed. The red line is the oil-water interface and the white line is the particle. (b-1) is the flow field when c(θ) = 1+0.1cosθ (symmetric concentration field) and (b-2) is the flow field when c(θ) = 1 + 0.2cosθ + 0.1sinθ (asymmetric concentration field). As a result of these process, there appears inhomogeneity of interfacial tention, and the droplet is driven by the Marangoni flow induced by the interfacial tension gradient. In a laboratoryframe, theparticle canfreelymove along theinterface. Without lossofgenerality, we may also consider a frame in which the particle and the center of mass of the droplet are fixed. Polar coordinates are introduced with the center of mass of the droplet regarded as the origin. The position of the particle is fixed at r = R and θ = 0. With the assumption of low Reynolds number, the velocity field is calculated using Stokes equation as follows, ∇p+η∇2v = 0, (1) − where v is the velocity, p is the pressure field and η is the viscosity coefficient. In our model, v satisfies the incompressive condition, ∇ v = 0. Due to the interfacial gradient, there is · the following stress jump at the interface [8]: 1 ∂γ (i) (o) σ = σ + , (2) rθ r=R rθ r=R R∂θ (cid:12) (cid:12) (cid:12) (cid:12) where σ is the stress tensor, the(cid:12) superscrip(cid:12)ts “(i)” and “(o)” correspond to the in- side and outside of the droplet, and γ(θ) is the interfacial tension between two flu- ids. Here, γ linearly depends on the interface concentration of the surfactant, c(θ) = ∞ c0 + (cc cosmθ +cs sinmθ), as γ(θ) = γ kc(θ). m=1 m m 0 − AsPsumingthatthesurfactantsuppliedfromtheparticleisdistributedonlyattheinterface 3 due to fast decomposition in the bulk, the time evolution of c(θ) can be described as ∂c u ∂c D ∂2c + = αc+παβδ(θ), (3) ∂t R∂θ R2∂θ2 − where D is the diffusion constant at the interface, α is the decomposition rate, and παβ is the rate of supply. In Eq. (3), the stationary flow field satisfying Eq. (1) is used as u, which is determined by c(θ). III. SOLUTION OF THE TWO DIMENSIONAL STOKES EQUATION In this section, v satisfying Eq. (1) is calculated under given concentration field at the interface, c(θ). In the polar coordinate, v is described as v = v (r, θ)e +v (r, θ)e , (4) r r θ θ where e and e are the unit vectors in the radial and angular directions, respectively. Due r θ to incompressibility, 1 ∂ 1∂v ∇ v = (rv )+ θ = 0; (5) r · r∂r r ∂θ therefore, using the stream function, Ψ, the solution of Eq. (1) can be described as 1∂Ψ v = , (6) r r ∂θ ∂Ψ v = . (7) θ − ∂r Taking rotation of both sides of Eq. (1) and substituting Eqs. (6) and (7), we can obtain ∇2∇2Ψ = 0, (8) where 1 ∂ ∂ 1 ∂2 ∇2 = r + . (9) r∂r ∂r r2∂θ2 (cid:18) (cid:19) 4 The solutions inside the droplet, Ψ(i), and the outside the droplet, Ψ(o), are r 2 Ψ(i)(r, θ) =B(i) 0 R (cid:16) (cid:17)r 3 r + Ac +Cc(i) cosθ 1 R 1 R (cid:18) (cid:19) (cid:16) (cid:17) r 3 r + As +Cs(i)) sinθ 1 R 1 R (cid:18) (cid:19) ∞ (cid:16) (cid:17) r m+2 r m + Ac +Cc cosmθ m R m R mX=2(cid:26)(cid:18) (cid:16) (cid:17) (cid:16) (cid:17) (cid:19) r m+2 r m + As +Cs sinmθ , (10) m R m R (cid:18) (cid:19) (cid:27) (cid:16) (cid:17) (cid:16) (cid:17) r 2 r Ψ(o)(r, θ) =B(o) +C ln 0 R 0 R (cid:16) r(cid:17) r r R + Bc ln +Cc(o) +Dc cosθ 1R R 1 R 1 r (cid:18) (cid:19) r r r R + Bs ln +Cs(o) +Ds sinθ 1R R 1 R 1 r (cid:18) (cid:19) ∞ m−2 m R R + Bc +Dc cosmθ m r m r ! m=2(cid:26) (cid:18) (cid:19) (cid:18) (cid:19) X m−2 m R R + Bs +Ds sinmθ . (11) m r m r ! (cid:18) (cid:19) (cid:18) (cid:19) (cid:27) (o) (i) The boundary conditions of flow field at the interface are v = v = 0 and r r r=R r=R (o) (i) (cid:12) (cid:12) v = v since the droplet shape is fixed. Since there is no(cid:12) external fo(cid:12)rce, a torque θ θ r=R r=R (cid:12) (cid:12) free(cid:12)condition,(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) T = dl n σ(o) r = 0, (12) · r=R × Z (cid:0) (cid:12) (cid:1) (cid:12) and a force free condition, f = dln σ(o) = 0, (13) · r=R Z (cid:12) (cid:12) have to be satisfied, where l is the distance along the interface. Using the conditions de- 5 scribed above and Eq. (2), the coefficients are calculated as (i) (o) B =B , 0 0 C =0, 0 Bc =Bs = 0, 1 1 κcs R Ac = Cc(i) = Cc(o) = Dc = m , 1 − 1 1 − 1 2 κcc R As = Cs(i) = Cs(o) = Ds = m , 1 − 1 1 − 1 − 2 κcs R Ac =Bc = Cc = Dc = m m m − m − m 2 (m = 0,1), 6 κcc R As =Bs = Cs = Ds = m m m − m − m − 2 (m = 0,1), (14) 6 where κ = k/ 2 η(o) +η(i) . From Eqs. (6), (7) and (14), flow profiles are obtained as (cid:8) (cid:0) (cid:1)(cid:9) ∞ κm R m−1 R m+1 v(o) = r − r r − h(cid:0) (cid:1) 2 (cid:0) (cid:1) i m=1 X (cc cosmθ+cs sinmθ), m m ∞ κm r m+1 r m−1 v(i) = R − R r − h(cid:0) (cid:1) 2 (cid:0) (cid:1) i m=1 X (cc cosmθ+cs sinmθ), m m (i) B v(o) = 2 0 r θ − R2 ∞ κ (m 2) R m−1 m R m+1 − r − r − h (cid:0) (cid:1) 2 (cid:0) (cid:1) i m=1 X (cc sinmθ cs cosmθ), m − m (i) B v(i) = 2 0 r θ − R2 ∞ κ (m+2) r m+1 +m r m−1 − R R − h (cid:0) (cid:1)2 (cid:0) (cid:1) i m=1 X (cc sinmθ cs cosmθ), (15) m − m 6 and the tangential flow at the interface (r = R), u(θ), is obtained as (i) ∞ B u(θ) = 2 0 R+κ (cc sinmθ cs cosmθ). (16) − R2 m − m m=1 X As the particle is a solid, the flow velocity on the particle is 0 (no-slip condition). As- suming that the particle size is small enough, this condition can be described as u(0) = 0. (i) Using this condition, B is determined as 0 ∞ κR2 cs B(i) = m. (17) 0 − 2 R m=1 X At infinity, v(o) = κ(cccosθ +cs sinθ)/2 and v(o) = 2B0(i)r κ( ccsinθ+cs cosθ)/2, r − 1 1 θ − R2 − − 1 1 whichmeansthedroplettranslationalvelocity(V )andtheangularvelocityinthelaboratory frame (Ω) are cc cs V =κ 1e +κ 1e , (18) x y 2 2 B(i) ∞ cs Ω =2 0 = κ m. (19) R2 − R m=1 X e and e in Eq. (18) are the unit vectors in the direction of x axis (θ = 0) and y axis x y (θ = π/2), respectively. The calculated flow fields are illustrated in the case of translational motion and in the case of rotational motion using Line Integral Convolution [9] in Fig. 1 (b). It is noted that the rotational speed of a droplet is proportional to the sum of the Fourier co- efficients corresponding to an antisymmetric concentration field about the anterior-posterior axis. IV. ANALYSIS OF THE MODEL Using Eq. (3), we calculate the velocity and the angular velocity of the droplet in this section. Since higher Fourier modes of c decay quickly due to diffusion, we neglect cc, and n cs for n 2. Letting cc, cc, and cs be replaced by X, Y and Z, Eq. (3) is described as n | | ≥ 0 1 1 dX κ = α(X X )+ Y2 (Y )2 +Z2 , (20) 0 0 dt − − R − dY κ (cid:8) (cid:9) =λ (Y Y ) Z2, (21) Y 0 dt − − R dZ κ =λ Z + (Y Y )Z. (22) Z 0 dt R − 7 Here, X = κβ2ρ2/R(ρ2 +1)+β/2, and Y = βρ2/(ρ2 +1) are the values of X and Y at the 0 0 fixed point at which Z = 0. Here, ρ = R/l is the radius normalized by the characteristic d length of diffusion, l = D/α. λ and λ are the eigenvalues for the linearized equations d Y Z for Y and Z, p 1 λ = α 1+ , (23) Y − ρ2 (cid:18) (cid:19) κ 1 αρ2 L 1 2 λ = Y α 1+ = 1+ , (24) Z R 0 − ρ2 ρ2 +1 ρ − ρ2 ( ) (cid:18) (cid:19) (cid:18) (cid:19) respectively, where L = l /l is the ratio between another characteristic length, l = a d a κβ/α, and l . l is the length of advective transportation during 1/α since u κcc d a ∼ 1 ∼ κβρ2/(ρ2 +1) κβ. If λ is positive, the fixed point, (X ,Y ,0), is unstable, which means Z 0 0 ∼ a finite angular velocity of a droplet obtained from Eq. (19). The phase diagramof the motion can be illustrated as shown in Fig. 2 (a) using two kinds of motion: translational motion and rotational motion. To realize rotation, it is necessary that the advection term in Eq. (3) is dominant; that is, L ' 1 is needed. In fact, λ can Z be positive only when L > 16√3/9 3.08. When R / l (ρ / 1), the concentration d ≈ field is almost uniform due to diffusion, so that flow is too weak to generate rotation. On the other hand, when R ' l (ρ ' L ), the uniform decomposition term in (3), αc, a − becomes dominant and the concentration field is almost uniform. Therefore, only when L is large enough and R is between Rc l and Rc l , can the droplet spontaneously l ≈ d u ≈ a rotate due to the coupling between the mirror symmetric concentration field (Y ) and mirror 0 antisymmetric concentration field (Z) through the rotational flow field. Since there is no other bifurcation, there are only a pair of stable fixed points with non-zero Z, which are calculated as X = 3αl /2κ, Y = Rλ /κ and Z = R√ λ λ /κ, when Rc < R < Rc. a − Y ± − Y Z l u Using Eq. (18) and Eq. (19), V and Ω in the steady state are calculated as ( λ l , 0) (R < Rc,R > Rc) V = − Y a l u , (25)   λYR, R√ λ λ (Rc R Rc)  − 2 ±2 − Y Z l ≤ ≤ u (cid:0)0 (cid:1) (R < Rc,R > Rc) l u  Ω =  √ λ λ . (26)  Y Z ∓ −    2 = α la 1+ ld 2 (Rc R Rc) ∓ R − R l ≤ ≤ u  r n (cid:0) (cid:1) o   8 (a) (b-1) (b-2) 10 8 rotational 6 motion L 4 translational 2 motion 0 0 2 4 6 8 10 r (c) (d) 4 1.5 2 1 0 | W W| 0.5 -2 -4 0 -10 -5 0 5 10 0 1 2 3 4 5 6 Lr 3/(r 2+1)-1 R FIG. 2. (Color Online) Bifurcation of the motion of a droplet. (a) Phase diagram of the motion. One-pointdashedline(blue)anddashedline(green) representR = l andR = l , respectively. (b) d a Schematic diagrams of the interface concentration field and the motion of a droplet in a laboratory frame. When the droplet moves straight (b-1) and rotates (b-2), the concentration field is mirror symmetric and asymmetric about the anterior-posterior axis, respectively. (c) Dependence of Ω˜ on Lρ3/ ρ2+1 2 1. (d) Dependence of Ω on R. α = 1, l = 1, and l = 5. With these values, d a − | | Rc =(cid:0)0.83, R(cid:1)c = 4.55 and Ω reaches a maximum value at R = 1.28. l u | | The interface concentrationandthemotionof thedroplet inalaboratoryframeareschemat- ically illustrated in Fig. 2 (b). When a droplet exhibits translational motion (R < Rc,R > l Rc), the interface concentration is symmetric about the anterior-posterior axis and the u anterior-posterior axis always corresponds to the direction of the motion. On the other hand, when a droplet rotates (Rc < R < Rc), the interface concentration is asymmetric l u about the anterior-posterior axis and the particle is always the inside of the trajectory. Con- sidering the normalized Ω, Ω˜ = Ω/λ , Ω˜ is the square root of the distance from the critical Y point, Lρ3/(ρ2 +1)2 1 , which corresponds to the stable solution of the normal form of | − | pitch-fork bifurcation [10], as shown in Fig. 2 (c). Figure 2 (d) shows the dependence of Ω | | on R when l , l and α are fixed. Ω has a maximum value at a radius between Rc and Rc. d a | | l u 9 (a) (b-1) (b-2) oleic acid 0 s 1 s 2 s sodium oleate 10 s 3 s 3 s 4 s 5 s 0 s 8 s aqueous phase FIG. 3. (Color Online) Outline of the experimental system. (a) Schematic diagram. A petri dish with a radius of 7.5 cm was filled with 100 ml water. A droplet of oleic acid together with a 3 mm particle of sodium oleate was placed on the water surface. The droplet moved spontaneously driven by the Marangoni effect [7]. (b) Trajectory of a droplet with a radius of 5.9 mm (300 µl). (b-1) Snapshots of a droplet per 1 s. The arrow indicates the particle fixed on the interface. (b-2) Trajectory of the center of mass of the droplet. Both scale bars represent 10 mm. V. COMPARISON WITH THE EXPERIMENT Using above theoretical results, we analyzed the corresponding experiment reported in [7] with additional data. The schematic diagram of the experimental setup is illustrated in Fig. 3 (a). One hundred milliliters of water, which was purified with a MilliQ filtering system (Millipore), was placed to a petri dish. A droplet of oleic acid (Wako Pure Chemical Industries; 159-00246) and a solid sodium oleate (soap) were floated on the aqueous phase in the petri dish. A solid column of sodium oleate (3 mm in length) was chosen from a commercial sample (Nacalai Tesque; 257-02). The movement of the solid/liquid composite was captured by a digital video camera at 30 frames per second at room temperature, and then analyzed using Image J (http://rsbweb.nih.gov/ij/docs/index.html). A time series of snapshots of a droplet, and the corresponding trajectory of the center of mass of the droplet are shown in Fig. 3 (b). Due to the fluctuation of rotational speed, the trajectory does not exhibit a closed circle. From the analysis of observed trajectories, the angular velocity, Ω, was measured. Time courses of Ω are shown in Fig. 4 (a). Five | | experiments were performed for each volume to yield the distribution of Ω , which is shown | | in Fig. 4 (b). To determine the maximum value of the distribution, the distribution was fit- ted with the function, exp b4Ω2(Ω2 a)+c , using the weighted nonlinear least-squares {− − } Levenberg-Marquardt algorithm, where the fitting parameters were a, b and c [11]. The number of events was taken as the weight for fitting. Using these parameters, the peak of the distribution, Ω , is calculated as 0 when a < 0, and a/2 when a 0. We used typ ≥ p 10