Theory for the coalescence of viscous lenses

Drop coalescence occurs through the rapid growth of a liquid bridge that connects the two drops. At early times after contact, the bridge dynamics is typically self-similar, with details depending on the geometry and viscosity of the liquid. In this paper we analyse the coalescence of two-dimensional viscous drops that float on a quiescent deep pool; such drops are called liquid lenses. The analysis is based on the thin-sheet equations, which were recently shown to accurately capture experiments of liquid lens coalescence. It is found that the bridge dynamics follows a self-similar solution at leading order, but, depending on the large-scale boundary conditions on the drop, significant corrections may arise to this solution. This dynamics is studied in detail using numerical simulations and through matched asymptotics. We show that the liquid lens coalescence can involve a global translation of the drops, a feature that is confirmed experimentally.


I. INTRODUCTION
Coalescence of drops is one of the most common capillarity-driven phenomena which can be observed in multiphase fluid dynamics. The early-time dynamics of coalescence is dependent on both the viscosity of the drops and their geometry. Different power laws for the growth of the connecting structure (referred to as neck or bridge) have been found for viscous and inviscid freely suspended drops [1,7,10,20,24], as well as for sessile drops in the viscous and inviscid limit [8,15,18,19,21]. The study of coalescence phenomena is also relevant for many applications where the underlying substrate of the coalescing drops is a liquid. Some examples are wet-on-wet printing [13], emulsions [17,22], and lubricant impregnated substrates [2,23].
Here, we focus on liquid lenses [12], consisting of liquid drops floating on a quiescent pool of another liquid. This case was studied for Newtonian drops [5] and liquid crystals [6], where the authors analysed the growth of the bridge in top-view experiments. Recent work considered the coalescence of lenses using side-view experiments [14]. This perspective is sketched in figure 1, providing a quasi-two-dimensional view of the problem. The experiments revealed a self-similar dynamics of the bridge profiles, with scaling laws for the bridge height h 0 with time t that depend on the viscosity of the lenses [14]. At low viscosity, the dominant balance during coalescence is between surface tension and inertia, and it was found that h 0 ∼ t 2/3 . At high viscosity, the dominant balance between surface tension and viscosity leads to h 0 ∼ t. These scaling laws are the same as those described in the merging of liquid wedges [3]. Owing to the slender geometry of the drops -typically the contact angle θ in figure 1 is small -the coalescence of liquid lenses can be analysed using the thin-sheet equations [11,25]. Using a similarity analysis, the experimentally observed inertial and viscous scaling laws are recovered [14]. For example, the viscous coalescence speed was found to be dh 0 /dt ≈V 0 = 0.5525γθ 2 /η, where γ and η respectively are the drop surface tension and viscosity (for consistencȳ x z θh u R Figure 1. Side-view sketch of two coalescing lenses initial radius R and equilibrium contact angle θ. The drop profile is described byh(x,t), while the minimal height in the neck regionh0(t) =h(0,t). The coalescence velocityV is defined asV = dh0/dt.  (c) Comparison of the center of mass velocityŪ and the "bridge" coalescence velocityV = dh0/dt for liquid lenses, taken for three oils of different viscosity and nearly identical contact angles (blue: η =115 Pa·s and θ = 27 • , green: η =33 Pa·s and θ = 32 • , red: η =9 Pa·s and θ = 31 • , see Hack et al. [14] for experimental details). The center of mass velocityŪ is not a small effect, as it is comparable in magnitude in comparison to the coalescence velocityV .
we here use overbars to indicate dimensional variables). This prediction was found to be in very good agreement with experiments.
The viscous similarity analysis, however, contains a salient feature that remains to be explained: the obtained self-similar velocity profile does not decay at large distance from the thin bridge region, but reaches a finite value. This is rather unusual for problems involving coalescence (or drop breakup, cf. Eggers & Fontelos [9]). Namely, the similarity analysis is typically based on the assumption that the flow remains confined to the scale of the bridge -at large scale, i.e. the scale of the drop, the flow is usually assumed to vanish. Such is the case for the coalescence of sessile drops as illustrated in figure 2(a). It shows an experimental top-view sequence of merging drops that are in contact with a solid substrate. During the initial growth of the bridge the global features of the drops appear nearly stationary -away from the bridge region one observes only a minute spreading of the drops. Figure 2(b) shows the equivalent top-view sequence for liquid lenses, for which the situation is manifestly different. Clearly, these floating drops do not remain stationary, but their centers of mass exhibit an inward motion as soon as the drops establish contact. This inward motion is not a small effect. Figure 2(c) compares the center of mass velocityŪ (measured in top-view) to the bridge coalescence velocityV = dh 0 /dt (measured in side-view), showing that the two velocities are proportional. These velocities were obtained using the experimental method described in Hack et al. [14]. We remark that this inward motion does not at all arise for liquid lenses of very low viscosity -this is in line with the inviscid similarity solutions, whose velocity rapidly decays away from the bridge [14].
In this paper, we provide a detailed analysis of the coalescence of highly viscous lenses and elucidate the coupling between the inner "bridge" solution and the global dynamics of the drops. We treat the problem using the two-dimensional thin-sheet equations, reflecting an analysis along the cross-section shown in figure 1. Such a twodimensional approximation turned out successful for the geometry of spherical caps [14,15,21], since close to the bridge the length scale in the third dimension, O((Rh 0 /θ) 1/2 ), is much larger than the horizontal and vertical length scales of the bridge, O(h 0 /θ) and O(h 0 ).
We demonstrate that, in general, the coalescence velocity exhibits significant logarithmic corrections, O(1/ ln t), where t 1 is the (dimensionless) time after coalescence. On the other hand, the thin sheet equations admit an outer solution where the drop's center of mass can migrate freely, closely resembling the motion observed in figure 2(b). In this latter case, the corrections to the leading-order result are much smaller. The analysis is confirmed in detail by comparison to time-dependent numerical simulations of the thin-sheet equations.
The article is organized as follows: The governing equations for coalescing lenses are presented in §II. In §III, we study the inner scale dynamics, near the point of coalescence, followed in §IV by an analysis of the outer region, which is where the two different boundary conditions manifest themselves. We end with our conclusion and outlook in §V.

A. Formulation
Following the approach of Hack et al. [14], the process of coalescence is modeled by the two-dimensional viscous thin-sheet equations [11,25]. The underlying approximations are the following: (i) Similar to sessile drops [15,21] the flow in the bridge region is quasi-two-dimensional in the early stage of coalescence. (ii) The equilibrium contact angle θ is small, such that a slender body approximation can be employed. (iii) The influence of the bath on the dynamics is negligible, i.e., free slip boundary conditions can be employed at both interfaces of the 2D lenses. (iv) Due to negligible differences in surface tension between the bath and liquid lens and the liquid lens and air, the liquid lenses are assumed to be symmetrical with respect to the bath-air interface -though asymmetric surface tensions can actually be mapped to an "effective" symmetric surface tension (cf. Hack et al. [14], Supplementary Material). We note that these assumptions are in accordance to the previous experiments of where the ratio of viscosities of the lenses and the bath is more than a thousand and the corresponding surface tension asymmetry is only about ten percent. The former justifies (iii) and the latter (iv).
The resulting model equations for negligible inertia of the flow, in dimensionless variables, take the form: Here, h(x, t) and u(x, t), respectively, are the dimensionless interface height and horizontal velocity, and dots and primes indicate derivatives with respect to t and x, respectively. Dimensional variables (x,h,ū,t) are scaled as where R is the initial lens radius and θ is the contact angle (cf. figure 1). The surface tension is denoted as γ and η is the viscosity of the liquid inside the lenses. The thin-sheet equations (1) correspond to mass and (horizontal) momentum conservation. The latter gives the balance between capillary forces (first term) and viscous forces (second term), while inertia has here been neglected. The momentum equation (1b) can be readily integrated to give the horizontal force balance The terms on the left-hand side represent the horizontal force transmitted through the thin drop by pressure, surface tension, and viscous stresses, respectively. We have introduced the constant −1/2 on the right-hand side corresponding to the total force in the static solution, see (6) below. Therefore F (t) measures any additional horizontal force that arises during the coalescence dynamics.

B. Two-dimensional numerical simulations
In order to illustrate the interplay between the dynamics of the small bridge region and the large bulk of the drops, we simulate the coalescence of two-dimensional lenses numerically. Due to symmetry about the point of coalescence x = 0, we can impose the boundary conditions h = u = 0 at x = 0 and focus on the domain x ≥ 0. We consider two different sets of outer boundary conditions. The first case corresponds to the experimentally realised setup of two freely floating lenses. In this case, the length L(t) of each lens decreases with time (as can be seen in figure 2) from its initial value L(t = 0) = 2, and at the edge of the lens we impose the thickness h(L), the contact angle, and conservation of mass: Substitution into the momentum equation (3)  The second case we consider is one where the lenses do not move, due to a symmetry condition being imposed about the centers of the lenses, such as in a periodic array of simultaneously coalescing lenses, This yields three boundary conditions h (0) = u(0) = h (1) = 0 on the governing equations and a fourth condition u(1) = 0 that determines the unknown F (t).
The key finding is that the different outer boundary conditions will lead to different spatial and temporal dynamics also in the inner region, at the scale of the bridge. As we shall see, the leading-order solution in the two cases remains the same -however, the coalescence velocity can exhibit significant corrections at the next order. The theoretical initial condition is taken to be the static solution of (1) with non-dimensional radius 1 and contact angle 1 (corresponding to the dimensional values R and θ, respectively), but in the numerical simulations a small-scale perturbation (h 0i − x/2) exp(−x/2h 0i ) is added to initialize the coalescence. The perturbation profile is chosen to satisfy the symmetry boundary condition h (0) = 0 and have an initial bridge height h 0i = h 0 (t = 0), which is taken to be 10 −10 unless stated otherwise. We solve the governing equations (1a) and (3) numerically in Matlab using a finite difference method with a Crank-Nicolson time-stepping scheme. Both the spatial and temporal grids are non-uniform with a relative resolution of 1%, and the discretisation is second-order accurate in both space and time. For the periodic lenses (5), the computational domain is 0 ≤ x ≤ 1. For the free lenses (4), where the domain 0 ≤ x ≤ L(t) changes with time, we use a rescaled position variablex = 2x/L(t) and solve the rescaled equations on the fixed domain 0 ≤x ≤ 2 instead. As explained above, four boundary conditions are imposed on the third-order system of equations, in order to determine the additional unknownL(t) or F (t).
Some snapshots of height and velocity profiles from the simulations are shown in figure 3. The free-floating lenses immediately begin to move towards each other at a constant velocity, and eventually merge into one larger lens with double the volume (and hence √ 2 times the radius and height). The periodic lenses instead flatten out towards a uniform height of 1/3, determined by the initial volume in the lens. Importantly, and this is one of the central points of the paper, the velocity profiles are significantly different between the two cases, even at very early times, due to the different outer boundary conditions. We will show how this affects the coalescence velocity,ḣ 0 .

III. THE INNER REGION
In order to study the dynamics of the region close to the point of coalescence, we rescale the variables as with τ = t. In choosing the scaling used above, we are motivated by the fact that in the inner region the length scale of importance is the bridge height h 0 . The resulting governing equations become where V =ḣ 0 is the unknown coalescence velocity that we wish to determine. At ξ = 0, the definition h(0, t) = h 0 (t) and the symmetry about ξ = 0 yield The system of equations (8) has three spatial derivatives and accordingly three boundary conditions (9). However, the equations contain two unknowns V and F that are determined on matching to the outer solution. The matching necessitates that we evaluate H and U as ξ → ∞, and from (8) we obtain, on neglecting the time-derivative term, This asymptotics contain two degrees of freedom α and β that are determined during the solution process. 1 This matching will be performed explicitly in §IV, but below we already anticipate some of the results necessary to evaluate the inner solutions.
A. The leading-order similarity solution The leading-order solution to the foregoing equations is a steady self-similar solution, which was previously obtained by Hack et al. [14]. The leading-order quantities H 0 (ξ), U 0 (ξ) and V 0 are independent of τ , so that (8) reduces to where we have anticipated that F (τ ) 1 -this indeed is true as we see later. The boundary conditions required to evaluate these quantities are given by (9), complemented by H 0 = 1 as ξ → ∞ (α = 1) set by the leading-order outer solution given in (6). These equations (and also the higher-order equations) are solved numerically in Mathematica using a shooting method on the domain 0 ≤ ξ ≤ 10 5 . We find the leading-order coalescence and far-field translational velocities to be [14] V 0 = 0.5525, B. Next-order corrections The solutions H 0 and U 0 are plotted in figure 4(a,b), where they are compared to the numerical solutions. The two columns correspond to the two types of boundary conditions, free-floating lenses and periodic lenses. For both cases and excellent agreement is found at early times after coalescence. However, the velocity profiles in figure 3(a,b) revealed large differences between these two situations. While these differences do not turn up in the leading order inner solutions H 0 and U 0 , they impact the next order corrections. As we will see, this also has a strong effect on the evolution of the coalescence velocity V . We therefore now address the inner solution beyond the leading order.
In the case of free-floating lenses, for which F (τ ) ≡ 0, it turns out that we can set up a consistent expansion based on the expansion parameter h 0 (τ ). On substituting the expansions (where the dependence of H and U on ξ is understood) into equations (8)-(9) we obtain at O(h 0 ), The differential equations for H 1 and U 1 are of third order with three boundary conditions, so one further condition is required in order to determine the unknown V 1 . This is obtained from a matching with the outer region. As ξ → ∞, the generic solutions to (14) are quadratic, H 1 ∝ ξ 2 , and so will need to match the second derivative h s (0) = −1 (an O(h 0 ) quantity when expressed in inner variables) of the outer solution. We thus impose H 1 → −1 as ξ → ∞, and obtain the numerical solutions plotted in figure 4(c), with the coefficients

The case of non-zero F (τ )
When the horizontal force F (τ ) is non-zero, the next-order corrections are at O(F ). As we shall see in the next section, these corrections dominate the O(h 0 ) corrections that arise due to time evolution of the outer solution, since F , though unknown yet, is evaluated to be O(1/ ln(h 0 )). For this reason, we now expand the solution in terms of F , i.e.
After substituting into the governing equations, we obtain, at O(F ), Since the corrections to the outer solution at O(h 0 ) are to be neglected when matching the inner and outer solutions at O(F ), we now impose the condition H 1 → 0 as ξ → ∞. On solving these equations numerically, we obtain the solution profiles plotted in figure 4(d) and the coefficients As F = O(1/ ln h 0 ) will typically not be very small, it is useful to calculate the second-order O(F 2 ) correction, which includes contributions proportional to F 2 as well as contributions proportional to (h 0 /ḣ 0 )Ḟ . We perform this calculation in appendix A.

IV. THE OUTER REGION AND MATCHING
We now turn to the outer solution, and calculate the corrections to the initial static profile h s (6) that are generated by the coalescence in the bridge region. Time integration of the evolution equation (1a) from this initial condition yields The The coefficients B(t) and F (t) are determined by boundary conditions and by matching to the inner solution.
which reveals that the first-order correction simply represents a translation of the initial profile h s (x) by the steady leading-order velocity u o < 0. (In fact, it can be shown that to all orders in h 0 , the velocity profile is spatially uniform and hence the drop is undergoing pure translation, with deformation only occuring in the bridge region x = O(h 0 ).) In order to match with the bridge region, we substitute x = h 0 ξ into (21) and expand in powers of h 0 , making use of Comparing this to the definition of the inner variables (7), we now identify the appropriate far-field behaviour of the inner solution to be This yields the condition H 1 → −1 as ξ → ∞ that we anticipated in §III B 1, which was imposed to obtain the solution (15). The resulting expression for the coalescence velocity is, from (12) and (15), predicting a linear correction to the coalescence velocity. This prediction is tested quantitatively in comparison to numerical simulations in figure 5. The result (24) is shown as a dashed line labeled "free", while the solid lines are numerics that were initialised at three different initial bridge heights h 0i . After a short time, the numerical curves for free lenses rapidly converge to the predicted asymptotics (24). Note that in typical experimental conditions where h 0 ∼ 10 −3 · · · 10 −2 , the coalescence velocity is very close to the leading order value V 0 .  (24) and (28).

B. The case of periodic lenses
The periodic lenses come with the symmetry boundary condition (5) on the velocity profile (20). This yields B(t) ≡ 0 and hence This logarithmic velocity at small x can indeed be matched to the velocity at large ξ from the inner region calculated in §III B 2. Specifically, we equate (25) (taking the limit x 1) with the far-field inner velocity profile (10) (using u = V U and ξ = x/h 0 ), where it is noted that V and β are expanded in F themselves. Equation (26) finally enables us to express F in terms of h 0 , which to leading order gives the result This confirms the slow, logarithmic decay of the force term, suggesting strong corrections with respect to the leading order similarity solution. Most importantly, this leads us to the sought-after coalescence velocity for periodic drops: where V 0 = 0.5525, V 1 = −0.6227 andV 1 = 4V 0 β 0 V 1 = 0.7892. The O(ln(2/h 0 ) −2 ) correction to the coalescence velocity is calculated in appendix A. Once again, the asymptotic prediction is tested quantitatively in comparison to numerical simulations in figure 5. The result (28) is shown as a dashed line labeled "periodic", while the solid lines are numerics that were initialised at three different initial heights h 0i ; again the numerics are in excellent agreement with the asymptotics. We notice a dramatic difference between the coalescence velocities for periodic drops as compared to the free drops, owing to the logarithmic decay of F . At values where h 0 ∼ 10 −3 · · · 10 −2 , the coalescence velocity for periodic drops is significantly below the leading order value V 0 .

V. CONCLUSION
In the present work, we have analysed the coalescence dynamics of viscous liquid lenses using both matched asymptotics and numerical simulations. We restricted our attention to the case where the flow in the bath is negligible, which is a consistent approximation for lenses of high viscosity, and where the contact angles are small to enable a slender thin-sheet description. In addition, following a previously used assumption in coalescence, we treated the problem as quasi-two-dimensional. The common scenario in coalescence and pinch-off, is that the flow remains localised into the narrow bridge region, while the far field remains stationary. Here we have found that this is not the case for viscous lens coalescence, and demonstrated that the bridge region affects the global dynamics, using both matched asymptotics and numerical simulations.
For freely floating two-dimensional viscous lenses, as soon as the lenses start coalescing there is a motion of the outer contact line towards the point of coalescence. In dimensional form, the growth rate of the bridge height h 0 (t) =h(x = 0,t) (24), and the horizontal translation velocity of the drop's center of massū o (12,15) can be written asV Here, we have eliminated the lens radius R in favour of the second derivativeh s (0) = −θ/R of the initial condition (6), in order to highlight that the coalescence process depends only on the local shape of the lens near the bridge. For a periodic array of two-dimensional viscous lenses, the centers of mass of the drops do not move due to symmetry. Instead, a horizontal force is generated that resists the inward motion towards the bridge, which results in a logarithmic velocity profile in the lens. In dimensional form, the growth rate of the bridge (28) and the additional force generated (27) are then given bȳ A further, second-order, correction is calculated in appendix A. It is important to note that the corrections are logarithmic in time, so these can be significant even at the very early stages of coalescence. In fact, logarithmic corrections also arise for viscous coalescence of freely suspended drops [10,16]. The structure of the problem is, however, different. In the freely suspended case the dynamics is h ∼ t ln t, which in contrast to (30) does not exhibit a finite velocity at early times. It is of interest to compare these two-dimensional, slender predictions to the experiments shown in figure 2(c). In Hack et al. [14] it was already shown that the experimental vertical coalescence velocityV was in very good quantitative agreement with the leading-order prediction of (29). Indeed, for freely floating drops, the higher-order terms are expected to be negligible in the experimental range. The current theory predicts a ratio of horizontal to vertical velocityŪ /V ≈ 0.57/θ for the freely moving drops, which for the experimental contact angles amounts tō U /V ≈ 1.1. The experimentally measured ratio in figure 2(c) was foundŪ /V ≈ 1.7, which implies an even stronger center of mass motion than predicted. This quantitative disparity could be due to three-dimensional effects, or due to the fact that the contact angles are not very small. Still, our theory offers an explanation for the appearance of a center of mass motion for viscous lenses, and provides the relevant scaling laws. This center of mass motion is not observed for inertial drops, which is in accordance with the decaying velocity in the inertial similarity solution [14]. For future experiments, it might be of interest to study coalescence while the centers of the drop are prevented from translating towards each other (e.g. by attaching the drops to fixed capillaries). In that case, we expect the appearance of significant (logarithmic) corrections to the coalescence velocity.
Besides three-dimensional effects, it is also of interest to discuss the influence of gravity. Gravity becomes important if the (dimensional) radius R of the lenses is no longer small compared with the capillary length c = γ/∆ρ g, where ∆ρ is the difference in density between the lenses and either of the external fluids. In this case, we expect gravity to flatten the static lens profile h s (x) (6), but near the rim of the lens, when x c , gravity becomes negligible and the capillary thin-sheet equations (1) and contact angle θ are recovered, and hence the leading-order inner solution (12) will still hold. For (two-dimensional) free lenses, which can undergo uniform translation even in the presence of gravity, the O(h 0 ) correction (29) also holds, but with a modified value of h s (0) that simply follows from the static drop. For periodic lenses, the correction becomes more involved as the outer velocity profile (20) is affected by the change in h s .
More generally, the dynamical structure of the problem bears a strong similarity with drop spreading on a rigid substrate, where the motion of the contact line also induces a weak flow on the scale of the drop [4]. The spreading velocity exhibits a logarithmic dependence on the scale separation between drop size and the characteristic scale of the contact line in that case too. An important difference, however, is that for drop spreading the universal leading-order similarity solution for the inner problem captures the phenomenon, as it has only algebraically small corrections. Here we found that for viscous lens coalescence, the corrections are themselves only logarithmically small, and therefore can be significant.
We are grateful to J. Eggers for discussions. We would also like to thank the referees for all their very valuable comments and suggestions, and in particular for highlighting the different large-scale boundary conditions. W.
The second-order coefficientV 2 is coincidentally quite small whenḣ 0 is expanded in terms of ln(2/h 0 ), so including the second-order correction in figure 5 would not have a noticeable effect. However, whenḣ 0 is expanded in terms of e.g. ln(1/h 0 ) or ln(1/t), the second-order correction yields a significant improvement.