Theoretical analysis for the optical deformation of emulsion droplets

We propose a theoretical framework to predict the three-dimensional shapes of optically deformed micron-sized emulsion droplets with ultra-low interfacial tension. The resulting shape and size of the droplet arises out of a balance between the interfacial tension and optical forces. Using an approximation of the laser field as a Gaussian beam, working within the Rayleigh-Gans regime and assuming isotropic surface energy at the oil-water interface, we numerically solve the resulting shape equations to elucidate the three-dimensional droplet geometry. We obtain a plethora of shapes as a function of the number of optical tweezers, their laser powers and positions, surface tension, initial droplet size and geometry. Experimentally, two-dimensional droplet silhouettes have been imaged from above, but their full side-on view has not been observed and reported for current optical configurations. This experimental limitation points to ambiguity in differentiating between droplets having the same two-dimensional projection but with disparate three-dimensional shapes. Our model elucidates and quantifies this difference for the first time. We also provide a dimensionless number that indicates the shape transformation (ellipsoidal to dumbbell) at a value $\approx 1.0$, obtained by balancing interfacial tension and laser forces, substantiated using a data collapse.


I. INTRODUCTION
Optical manipulation of fluid interfaces having ultra-low interfacial tension was realised recently in systems of micron-sized emulsion droplets [1,2]. Potential applications of such techniques lie in pumping fluids, performing reaction chemistry at the attolitre scale and understanding the behaviour of oil droplets in surfactant-enhanced oil recovery [3,4]. Similar techniques have been applied to deform cells, which have a non-zero bending energy of the lipid bilayer in addition to the usual interfacial tension [5][6][7][8].
The precision and non-destructive nature of manipulating particles by optical tweezers can be a benefit over standard mechanical techniques. In most physically realisable situations the radiation pressure exerted by the beam on the trapped particle is orders of magnitude weaker than the Young's modulus of the solid or the Laplace pressure of the confined fluid. As a result, particles do not deform in the trap: a solid particle retains its initial shape while a fluid droplet assumes a spherical geometry. Therefore, historically many applications of optical tweezing have been limited to exerting [9][10][11][12] and measuring [13][14][15] external forces on a rigid trapped object. The situation is, however, different for fluids with low interfacial tension, where the object itself can be deformed in response to the laser field.
Ashkin and Dziedzic reported small but measurable deformations on planar soft surfaces using a focused and pulsed laser beam [16]. Further, Casner and Delville showed that large deformations could be observed if the interfacial tension was lowered [17]. It is known that the interfacial tension can be lowered substantially by the addition of surfactants. Ward et al. showed that the interfacial tension of heptane droplets could be lowered to γ ≈ 10 −5 − 10 −6 N m −1 by the addition of suitably selected surfactants [1]. As a result, they were able to deform the emulsion droplets into various shapes using multiple continuous-wave lasers. Our theoretical work presented here is motivated by these experiments.
The problem is similar in spirit to that of lipid vesicles deformed via micro-manipulation techniques [18] where a bending free energy of the lipid bilayer is balanced against externally applied forces to give its resulting shape. The difference between the membrane and fluid geometry arises from the nature of the constraints imposed on the shape equations. In the vesicle case, the total area is conserved while for a fluid droplet the volume of enclosed fluid is constant throughout the deformation process. Furthermore, based on dimensional arguments one can show that so long as κ ξ 2 < γ, where κ is the bending modulus and ξ is a length scale based on the radius of curvature of the droplet surface, in the hydrodynamic limit the interfacial tension contribution dominates over the curvature energy contribution.
We propose a numerical model that predicts the three-dimensional steady-state shapes of droplets with ultra-low interfacial tension under the influence of one or more optical traps. Several authors have recently published models exploring different regimes of optical deformation of liquid droplets [19,20]. The key feature of the model we propose here is that it does not assume small, linear deformation of the droplet. Our model makes no assumptions on the final shape of the droplet and there is no restriction that the optical traps have to be focused at the centre of the droplet. Our model allows us to investigate the equilibrium shapes of droplets as a function of different parameters, such as laser power (P 0 ), numerical aperture (N A), initial droplet radius (R d ), interfacial tension (γ) and the parameter (ñ) which is related to the ratio of refractive indices asñ = 1 − n2 n1 , where n 1 and n 2 are the refractive indices of the droplet and external media, respectively. From these results we have defined a dimensionless deformation number N d which is a function of P 0 , N A, R d , γ andñ. Using this dimensionless number we are able to predict previously undescribed shape transformations of a droplet in a single optical trap.
The mathematical model is described in the Section II. This is followed by an outline of the numerical implementation of the model in Section III. We then present and discuss the results of the deformation of a droplet in single and multiple optical traps. Finally, we discuss future applications and extensions to the present model.

II. MATHEMATICAL MODEL FOR DROPLET DEFORMATION
In our model we describe the surface of the droplet (the interface between the liquid droplet and its host medium) in terms of a single-valued function R(θ, φ) in spherical polar coordinates, where R defines the distance between the interface and a preassigned fixed origin. In the absence of any external forces, a liquid droplet assumes a spherical shape, as a result of two antagonistic forces: the interfacial tension which tends to minimise the area, and the bulk pressure of the internal fluid which translates into a force acting along the local normal to an infinitesimal surface element dS. The energy function for an isolated droplet can be written as: Thus, the internal pressure for an isolated droplet at equilibrium turns out to be P int = 2γ R d . In the presence of one or more optical tweezers, the total pressure acting at a position r on an infinitesimal surface element dS at the oil-water interface can be written as: P (r) = P lap (r) + P opt (r) + P int (2) where P lap is the Laplace pressure, P opt is the optical pressure due to momentum transfer from the laser field to the interface, and P int is the internal pressure within the droplet, which we will treat as an unknown variable dependent on the specific experimental configuration, but uniform throughout the droplet since we are only considering equilibrium structures. The Laplace pressure P lap (r) acting on a surface element on the oil-water interface (interfacial tension γ), at a position r and with outward normaln(r), is calculated using the Young-Laplace equation: from using Gauss' theorem on the first integral in Eq. 1. Note that by convention this pressure acts inwards in the case of a convex droplet surface. The optical pressure P opt (r) across the oil-water interface is calculated by considering the momentum transferred to the interface as light is reflected and refracted at the interface. An exact solution for the light field in the presence of a dielectric object could be calculated in the form of a series expansion within the framework of T-matrix theory [21]. However, due to the excessive computational demands that this would impose, we instead adopt a localized ray-optics approximation similar to that described in [22]. We generalize the treatment to an arbitrary shape, but at the same time we note that the refractive index ratio m = n2 n1 between the oil and water phases is close to one, such that the phase shift of rays crossing the droplet is small compared to the wavelength of the trapping laser and we are close to the Rayleigh-Gans regime [23]. We will therefore make the approximation that the field is unperturbed by the presence of the droplet and calculate the localized Fresnel reflection and refraction for each surface element using the unperturbed laser field, without considering higher-order reflections.
If we consider a surface element on the oil-water interface, at a position r and with outward normaln(r), then for an incident beam with momentum density p 0 n 1ŝ , whereŝ is the Poynting vector, we can apply standard laws of reflection and refraction [24] and, after some algebraic manipulation, we find the optical pressure acting on the surface (in the direction of the outward normal) to be: where µ =ŝ ·n, F r and F t are the Fresnel power reflection and transmission coefficients for the angle of incidence θ inc = arccos(|µ|), and n 1 and n 2 are the refractive indices for the media in which the incoming and refracted beams are propagating. We also require a description for the momentum density p 0 n 1ŝ of the beam. The fact that our droplets are several wavelengths in diameter, significantly larger than the trapping beam waist, implies that a series expansion around the beam focus, such as that presented in [25] and used in surface optical stress calculations in [22], is not suitable here. However, observing that for the droplet radii we are interested in the interface is generally several Rayleigh lengths from the (tight) laser focus, it becomes apparent that a far-field scalar model of a tightly focused Gaussian beam is appropriate [26].
We can now return to the pressure acting on the surface of the droplet (Eq. 2). In order to obtain equilibrium droplet shapes we employ a relaxation dynamics scheme at each "spoke" of the coordinate system. Thus the time evolution of each point on the surface of the droplet is given by: where the (n ·r) term takes account of the fact that the interface normal is not collinear to the (fixed) radial direction along which R(θ, φ) is measured -see Figure 1, and hence R will in general vary faster than the rate of normal motion of the interface. The friction coefficient β dictates the inverse time scale for a droplet to relax to equilibrium. This, however, cannot be chosen to be arbitrarily large as it compromises the stability of the numerical scheme used to integrate Eq. 5. This equation of motion represents overdamped motion. Thus our model does not attempt to predict the detailed transient motion of the droplet surface over time as it approaches its final shape, since we have deliberately selected a simplified equation of motion that does not take complete account of the hydrodynamic environment of the droplet. Although a final solution to our model represents a physically correct steady-state solution for a given set of experimental parameters, we note that the "time" variable in our calculations does not represent real-world time, the trajectory taken through parameter space to attain the final solution is not physically rigorous, and if experimental parameters support multiple steady-state solutions then our model may not be able to predict for certain which solution represents the kinetic product from a given set of initial conditions. The equation of motion is integrated with respect to time in order to determine the physically correct steady-state solution, subject to a constant-volume constraint. The volume of the droplet is defined as: and hence the boundary condition can be expressed as the constraint: Substituting Eq. 5 into this, we obtain an expression for the free parameter P int representing the internal pressure of the droplet: For completeness we note that once the shape of the droplet has converged, as a result of these pressure contributions, by removing the presence of the optical tweezers from our model we return to the energy function of an isolated droplet (Eq. 1). As a result the deformed droplet converges back to a perfect sphere with radius R d .

III. NUMERICAL IMPLEMENTATION
In our model the continuous surface of the droplet is sampled at a certain number of discrete points. These points lie on a regular grid in spherical coordinates, with For each node a radius R ij is defined, each representing a point (R ij , θ i , φ j ) on the surface of the droplet (the oil-water interface), so that the radius of each node defines a "spoke" extending from the origin in the direction (θ, φ) -see Figure 1. We then solve the equation of motion (Eq. 5) on this discretized grid, using finite-difference expressions as an estimate for local derivatives on the surface of the droplet.
We can calculate the surface normaln as follows [27]: where r is the radial coordinate and R (θ, φ) is the radius of a "spoke" at a given (θ,φ) coordinate. At first glance it might appear convenient to implement Eq. 3 and 9 in a two-stage finite-difference scheme, however in practice this two-stage approach leads to a decoupling between odd and even points in the grid, which is difficult to eliminate. The solution is to derive a single expression for the Laplace pressure directly in terms of derivatives of the radius. This combined expression proves to be rather complicated in spherical coordinates: where: Now, by evaluating Eq. 4 and 10 with the help of second-order-accurate finite-difference expressions, and numerically integrating Eq. 8, we are able to numerically integrate the equation of motion (Eq. 5) with respect to time, using the Runge-Kutta-Fehlberg method, until motion ceases and the droplet has converged to a stable shape.
To improve computational efficiency we initially use a mesh size of N θ = 25 and N φ = 24 nodes. Once the droplet shape has converged using this coarse mesh, we refine to N θ = 51 and N φ = 48. Figure 2 shows a comparison between the converged shapes of a droplet with R d = 5.0 µm in an increasing number of optical traps. It can be seen that the refinement improves the shapes of the droplets in the regions of high curvature where the optical traps are positioned, whilst only minor improvements can be observed for the rest of the droplet. A mesh size of N θ = 13 and N φ = 12 is also shown in green in Figure 2. Such a coarse mesh is incapable of ensuring a smooth surface is obtained, while the initial mesh size of N θ = 25 and N φ = 24 nodes represents a good compromise between computational demands and accuracy, being able to give a reasonable initial estimate of the droplet shape that can later be refined. Furthermore, the calculation of the Laplace pressure would be erroneous with an insufficient number of mesh points. Our choice of N θ and N φ is also verified by considering the well-known fact that the integral of the mean curvature H over a closed surface S is equal to zero [28]: As the number of N θ and N φ nodes is increased in our model, this calculated quantity will tend to zero.

IV. RESULTS AND DISCUSSION
Continuous-wave, unpolarized, lasers having a wavelength λ = 1064 nm have been modelled for results reported in this paper. Measured values of the refractive index of a heptane droplet and an external water medium, n 1 = 1.38 and n 2 = 1.33 respectively, have been used in our calculations, unless otherwise stated.

A. Single Optical Trap
We begin our discussion with the deformation of a droplet in a single optical trap. The optical trap is positioned such that the centre of the droplet initially coincides with the focal point of the beam. The steady-state profiles are obtained when the "spoke" positions parameterising the surface do not change as a function of time. Figure 3 shows the converged three-dimensional shape of a deformed droplet with initial spherical geometry R d = 2.0 µm, and a surface tension γ = 10 −6 N m −1 , as a function of increasing laser power P 0 . The radius of the beam waist for each optical trap was kept constant at w 0 = 0.282 µm, which corresponds to a numerical aperture of N A = 1.20. As the laser power increases the droplet elongates to assume a lozenge form, with its long axis parallel to the direction of propagation of light (taken to be along the +z axis). For laser powers P 0 0.055 W the droplet has a dumbbell-like shape with an hour-glass connecting two spherical caps. For these configurations one of the principal curvatures is negative. Experimental observations of these deformations have so far been limited to their two-dimensional projections along the axis of laser propagation [1]. Note that the shapes are asymmetric with respect to z. If we consider the gradient force alone then we would expect a symmetric shape, reflecting the symmetry of the intensity distribution in the trapping laser field. However the scattering force (which acts in the local direction of propagation of the laser field) breaks this symmetry and pushes the interface, and hence the droplet, in the direction of propagation of the laser and leads to "bulging" of the droplet beyond the laser focus. The colour scheme used in Figure 3 indicates the variation in the calculated mean curvature over the surface of the droplet. For low laser powers the mean curvature is close to 0.5 µm −1 corresponding to that of a perfect sphere with radius R d = 2.0 µm, except at the "tips" of the droplet. As the laser power is increased we see a larger variation in the mean curvature, with the lowest mean curvature at the waist of the droplet, and the highest mean curvature at the two spherical caps. For laser powers greater than P 0 = 0.12 W the mean curvature at the waist becomes negative.
Since current experimental observations are limited to two-dimensional xy projections, the deformation of these droplets in a single optical trap is only known due to a decrease in the maximum projected radius [1]. The variation of the maximum projected droplet radius R xy as a function of the laser power P 0 for a constant value of interfacial tension γ = 10 −6 N m −1 and initial droplet size R d = 5.0 µm for different values of numerical aperture N A is shown in Figure 4a. R xy changes non-monotonically as a function of P 0 , with the minimum occurring at the transition to a dumbbell-like form, a transition that we will discuss in more detail below. The initial linear decrease in the slope of the R xy vs. P 0 curve is due to the elongation along the z-axis. Beyond the transition point, the equilibrated configurations show a narrowing of the neck region connecting the two spherical caps. A similar trend is observed for different values of the initial droplet size R d as shown in Figure 4b. The laser power at which a droplet undergoes a transition to a dumbbell-like shape, in a single optical trap, is found to be linearly dependent on both the initial radius of the droplet (R d ) and the interfacial tension (γ), and inversely proportional to the refractive index ratioñ. The relationship between the laser power and the numerical aperture is slightly more complicated. We find that the laser power P 0 ∝ exp (N A). These relationships allow us to define a dimensionless number which characterises when a droplet will deform into a dumbbell-like shape, for a given initial radius, surface tension, refractive index ratio and numerical aperture: When N d 1.0 the droplet deforms into a dumbbell-like shape, whereas for N d 1.0 the droplet only elongates in the direction of the propagation of light. Figure 5a shows this for a constant initial spherical radius R d = 5.0 µm whilst varying the numerical aperture for the optical tweezer. Figure 5b shows the variation of R xy vs. N d for different starting radii, whilst using a numerical aperture of N A = 1.10. To obtain a data collapse for the graphs in Figure 5b we rescale the observed R xy values as:  The transition to a dumbbell-like shape can be explained in terms of a balance between the forces due to the optical tweezer and the interfacial tension of the droplet. At values of N d 0.7 the interfacial tension energy is stronger than that of the optical tweezer, resulting in a linear deformation of the droplet. Between N d ≈ 0.7 and 1.0 the forces exerted by the laser field begin to overcome the interfacial tension of the droplet, until N d ≈ 1.0 at which the optical forces dominate. The droplet interface then conforms locally to the iso-intensity contours of the laser beam, which results in the observed dumbbell-like geometry.
Beyond a certain laser power our model is not capable of accurately predicting the converged shapes of droplets. This breakdown of the model is due to a single radial "spoke" intersecting the interface multiple times. This critical point for our numerical model will occur as the droplet elongates further and the dumbbell-like shape becomes more defined. Our model as described here is implemented using a spherical coordinate system, where each "spoke" radiates out from the origin of the system. Refining the number of "spokes" can allow the model to approach this critical point more closely (albeit at larger computational cost), but the model as presently formulated is unable to pass this critical point.

B. Multiple Optical Traps
In addition to calculating the shape of a droplet in a single optical trap, our model is able to successfully calculate how a droplet will deform in multiple optical traps. To our knowledge this is the first time such a model has been presented. Each optical tweezer is moved slowly outwards by 0.25 µm in the appropriate directions, and we ensure that the droplet shape has converged at each stage before the lasers are moved again.   Figure, moving from left to right, there is an increase in the separation between each optical trap. Initially, all traps are positioned at the centre of the droplet. The total laser power acting on the droplet is kept constant at P total = 0.12 W and is equally shared between the total number of lasers being modelled. The interfacial tension and numerical aperture were also kept constant at γ = 10 −6 N m −1 and N A = 1.20 respectively. The top row of each figure shows the two-dimensional xy projections of the deformed droplet, as might be seen in a brightfield microscope image. The chosen colour scheme indicates the variation in the calculated mean curvature; as the separation of each laser increases there is a larger variation in the mean curvature. For the largest separation between optical traps, presented in Figures 6c, 7c & 8c the mean curvature is highest at the "tips" or "corners" of the shapes formed, where the tightly focused lasers are positioned.
As mentioned above, a major benefit of these numerical results is that it is possible to view the three-dimensional geometries of the deformed droplet. These are shown in the bottom row of each of the Figures 6 -8. To date such observations have not been achieved experimentally. When all lasers are positioned at the centre of the droplet then P total = P 0 and we obtain a value of N d = 2.18 for the given R d , γ and N A. Hence this is beyond the transition to a dumbbell-like shape, as shown in Figures 6 -8a. The hour-glass connecting the two spherical caps has a negative mean curvature. As each individual laser is then moved outwards from the centre, the surface of the droplet at the "tips" or "corners" retains this concave shape. This observation is interesting since from just viewing the two-dimensional projections one might assume the surfaces to be convex along the axis of propagation of light.
It is worth mentioning that for suitably selected parameters we can observe a range of deformed droplets with a high proportion of negative mean curvature. For example, a droplet with initial radius R d = 5.0µm and interfacial tension γ = 10 −6 N m −1 , deformed in four optical traps each with a laser power of P 0 = 0.1 W and numerical aperture N A = 0.8 has negative mean curvature on the top and bottom faces, in addition to that seen on the side faces in Figure 8. This indicates that the deformation of droplets by optical forces is extremely sensitive to the selected parameters.
As an experimental validation of our model, we compare our calculated two-dimensional shapes of a deformed emulsion droplet to those presented in the work of Ward et al. [1]. To summarise the parameters used in their experiments: an approximate value of R d = 2.5 µm was taken for the initial spherical geometry of the droplet, and an interfacial tension of γ ≈ 10 −6 N m −1 was reported. For the optical traps, each laser had a numerical aperture of N A = 1.20 and a total combined laser power of P total = 24 mW was equally distributed between the total number of lasers. The experimentally observed shapes and steady-state shapes predicted by our model are presented in the first and second rows of Figure 9 respectively, for an increasing number of optical traps.
As mentioned above, the highest calculated mean curvature values occur near the focal positions of lasers. It can be seen from Figure 9 that using the experimental parameters of Ward et al. we obtain droplets with conical ends. The high electric field at the laser foci exerts a very strong force on the interface, and when combined with volume conservation and surface curvature considerations, this results in a locally very small radius of curvature at the tips of the shape, as previously modelled by Stone et al. [29] based on observations of "Taylor cones" in static electric fields [30].
The bottom row of Figure 9 show the three-dimensional geometries of these emulsion droplets. Unlike the structures presented in Figures 6-8, the "tips" or "corners" of each droplet remain convex, as expected from the value of N d = 0.35 deduced from the parameters reported by Ward et al. [1]. Hence when all the lasers are positioned at the centre of the droplet at the start of the experiment, the droplet remains approximately spherical, and does not take on the dumbbelllike structure discussed earlier. This in turn reduces the deformation of the droplet along the light propagation axis at the "tips"/"corners" as each laser is moved outwards.
There is a very strong agreement between our theoretical predictions of the two-dimensional xy projections of the deformed droplets and those obtained experimentally. We hope that this work will prompt experimental studies to visualise the three-dimensional configurations as a function of the physical parameters explored in this paper.

V. CONCLUSION
In conclusion we have developed a theoretical framework to compute three-dimensional equilibrium shapes of liquid droplets with ultra-low interfacial tension in one or more optical traps. Taking a cue from experiments, we assume an isotropic, temperature-independent interfacial tension coefficient γ. The optical traps were described using a far-field scalar model of a tightly focused Gaussian beam, within the Rayleigh-Gans regime.
The equilibrium droplet shape arises as a result of the interfacial tension trying to minimise the droplet surface area, the internal fluid pressure resisting such a change, and the external optical pressure causing local deformations, without any volumetric change. Using this model we numerically compute the droplet shapes as a function of both droplet and laser parameters e.g. interfacial tension, initial droplet size, laser power and numerical aperture.
We obtain droplet shapes similar to those obtained experimentally by Ward et al. for similar parameter values. The close agreement between theoretical predictions of two-dimensional projections in the xy plane and experiments for known geometries and parameter values gives us confidence in the correctness of our predicted three-dimensional droplet shapes. It is worth noting that these predictions of three-dimensional droplet shapes for large deformations (where the linear response breaks down) have not been reported in the literature. Experiments are currently under development to measure three-dimensional shapes for comparison with our model Further, the theoretical work revealed a few surprises along the way. For example, in a single optical trap as the laser power is increased the trapped droplet not only elongates in the direction of light propagation, but also deforms into a dumbbell-like shape. We have predicted the onset of this transformation in terms of a dimensionless deformation number N d , obtained using dimensional arguments balancing antagonistic surface tension and optical forces. The mathematical expression has been substantiated using a data collapse of the numerical solution of the shape equations.
We have also predicted droplet shapes having a total negative mean curvature along its faces for optical traps having a high intensity. Though for some configurations (e.g. four laser traps) this may be intuitive, the existence of droplet shapes having a total negative mean curvature in a single trap in not immediately obvious.