A novel coarsening mechanism of droplets in immiscible fluid mixtures

In our daily lives, after shaking a salad dressing, we see the coarsening of oil droplets suspended in vinegar. Such a demixing process is observed everywhere in nature and also of technological importance. For a case of high droplet density, domain coarsening proceeds with interdroplet collisions and the resulting coalescence. This phenomenon has been explained primarily by the so-called Brownian coagulation mechanism: stochastic thermal forces exerted by molecules induce random motion of individual droplets, causing accidental collisions and subsequent interface-tension driven coalescence. Contrary to this, we demonstrate that the droplet motion is not random, but hydrodynamically driven by the composition Marangoni force due to an interfacial tension gradient produced in each droplet as a consequence of composition correlation among droplets. This alters our physical understanding of droplet coarsening in immiscible liquid mixtures on a fundamental level.

In our daily lives, after shaking a salad dressing, we see the coarsening of oil droplets suspended in vinegar. Such a demixing process is observed everywhere in nature and also of technological importance. For a case of high droplet density, domain coarsening proceeds with interdroplet collisions and the resulting coalescence. This phenomenon has been explained primarily by the so-called Brownian coagulation mechanism: stochastic thermal forces exerted by molecules induce random motion of individual droplets, causing accidental collisions and subsequent interface-tension driven coalescence. Contrary to this, we demonstrate that the droplet motion is not random, but hydrodynamically driven by the composition Marangoni force due to an interfacial tension gradient produced in each droplet as a consequence of composition correlation among droplets. This alters our physical understanding of droplet coarsening in immiscible liquid mixtures on a fundamental level.
Phase separation is one of the most fundamental physical phenomena that lead to pattern formation in various types of systems [1][2][3][4][5][6], including all kinds of classical condensed matter, biological systems (both ordinary [7] and active matter [8]), quantum systems such as liquid 3 He-4 He mixtures [9] and Bose-Einstein condensates [10,11], nuclear matter [12], and cosmology [13]. In many of these systems, momentum conservation plays a crucial role in phase separation. The simplest among such systems is a classical binary liquid mixture, where material can be transported by both diffusion and hydrodynamic flow. Thus, the relevant field variables are the composition of one of the components φ and the fluid flow velocity v. The process of phase separation is then described by the Ginzburg-Landau-type φ 4 free energy and a set of coupled dynamical equations of these two variables satisfying the composition and momentum conservation under the incompressiblity condition ∇ · v = 0 (see Methods for the explicit forms of the equations). This model is widely known as model H in the Hohenberg-Halperin classification of dynamical critical phenomena and phase separation kinetics to dynamic universality classes [14]. Phaseseparation behaviour in many of the above-mentioned systems can basically be described in a similar framework. Thus, revealing the fundamental physical mechanism of domain coarsening in model H is of significant importance for the general understanding of spinodal decomposition in these systems.
From a technological point of view, since phase separation starts at nanoscale and its characteristic size proceeds to grow indefinitely with time, we can control the domain size of materials intentionally by freezing the process at a certain time by using liquid-solid transitions such as glass transition and crystallization, by chemical * E-mail:tanaka@iis.u-tokyo.ac.jp reaction, or by emulsifying liquid droplets with surfactants. This strategy is widely used in industrial processes of various materials such as liquid mixtures, emulsions, polymer blends, metallic alloys, ceramics, cosmetics, and foods. Among various types of phase separation, liquid phase separation is quite important in materials science since a liquid state is most suitable for materials processing. Thus, the understanding of domain coarsening dynamics in the demixing process of liquid mixtures is not only of fundamental importance but also of crucial technological importance for controlling domain structures of materials.
The theoretical framework of model H has a very firm basis since it replies solely on the composition and momentum conservation, and there is little doubt on its validity. This phase-field model is widely used for numerical simulations not only in physics, but also in materials science and chemical engineering communities. However, the analytical theoretical analysis of model H is not an easy task due to its non-locality, non-linearity, complex non-local dynamical coupling between φ and v, and stochastic nature due to thermal noise. Nevertheless, the concept of dynamical scaling provides us with a powerful scaling argument for the self-similar domain growth in the late stage of phase separation after the formation of a sharp domain interface [1][2][3][4]15], on the basis of clear physical pictures on the elementary process of domain coarsening. Thus, the basic mechanisms of domain coarsening in immiscible liquid mixtures are now believed to be reasonably understood [1][2][3][4][5][6]. Recently, research interests have shifted towards other aspects of phase separation including aging dynamics during the coarsening process [16], geometrical features [17], and response functions of coarsening systems [18], inertia effects [15,19], and viscoelastic effects [4,20].
Here we summarize the current understanding of liquid phase separation. We consider a binary liquid mixture, whose components have the same viscosity η. A schematic phase diagram is shown in Fig. 1a, together with typical phase-separation structures (Fig. 1b,c). The relevant domain coarsening mechanism depends solely on the volume fraction of the minority phase Φ: For small Φ (case A), the evaporation-condensation (Lifshits-Slyozov-Wagner (LSW)) mechanism [21][22][23] is responsible for droplet coarsening: steady diffusion flux from smaller droplets to nearby larger droplets leads to the growth of the latter at the expense of the former. In this case, the average domain size R grows as R(t) ∼ = [(4/27)D 0 ξt] 1/3 [4], where t is the time, D 0 is the diffusion constant of a component molecule (or atom), and ξ is the correlation length of composition fluctuations (or the interface thickness). In this mechanism, translational diffusion of component molecules (or atoms) is the main transport process. Near a symmetric composition (Φ ∼ 1/2) (case B), bicontinuous spinodal decomposition takes place (Fig. 1b). There, the hydrodynamic coarsening (Siggia's) mechanism characteristic to a bicontinuous pattern leads to rapid coarsening [24]: [4], where σ is the interface tension between the two phases. In this mechanism, hydrodynamic flow is the main transport process. For intermediate Φ (case C), droplet spinodal decomposition takes place (Fig. 1c). It has been believed that the Browniancoagulation (Binder-Stauffer-Siggia (BSS)) mechanism is mainly responsible for the coarsening [24,25]: droplets grow by accidental collisions between droplets undergoing random Brownian motion, which leads to the coarsening law R(t) ∼ = [(6k B T Φ/(5πη))t] 1/3 [4,26], where k B T is the thermal energy. In this mechanism, thermal Brownian motion of droplets (rather than molecules) is the main transport process. The growth exponents of these coarsening laws have been confirmed by both experiments [27][28][29][30][31][32][33] and numerical simulations [19,[34][35][36].
The coarsening mechanisms for case A and B are rather well established, whereas that for case C is not so clear due to difficulties associated with many-body effects and hydrodynamic effects. The simplifications of case A and B come from the fact that for case A transport is dominated by diffusion and for case B by hydrodynamic flow. On the other hand, both may play an important role for case C in a complex manner. Although there has been a wide belief that the primary coarsening mechanism is the BSS mechanism, additional mechanisms have also been suggested. For example, Tanaka observed the process of droplet coarsening in a binary liquid mixtures confined between two glass plates with optical microscopy, and found [26,37] that droplet collision enhances subsequent collisions via hydrodynamic flow (phenomenon 1) and droplets sometimes move directionally, contrary to the expectation from random Brownian motion (phenomenon 2). Phenomenon 1 is enhanced for a quasi-twodimensional (quasi-2D) situation due to stronger hydrodynamic interactions [38], but even for 3D its importance has been pointed out when the droplet volume fraction is very high [36,39]. On phenomenon 2, some theories were proposed to explain directional motion of droplets by a In the (light yellowgreen) metastable region (i.e., the low Φ region (Φ <∼ 0.21)), nucleation-growth-type phase separation takes place and there droplets coarsen by the evaporation-condensation (or, LSW) mechanism: diffusional transport from smaller to larger droplets. Near the symmetric composition (the light blue region, or, Φ ∼ 1/2), bicontinuous spinodal decomposition (b) takes place and the domain coarsening is governed by the Siggia's hydrodynamic coarsening mechanism: hydrodynamic transport from thin to thick parts of bicontinuous tube structures. In the intermediate composition region (the green region), droplet spinodal decomposition (c) takes place and its coarsening mechanism has been believed to be the Brownian coagulation mechanism: random Brownian motion of droplets and the resulting collisions. Contrary to this, we propose a new coarsening mechanism, in which droplet motion is induced by the composition Marangoni effect due to the diffusional composition coupling among droplets, and not by random Brownian motion. The border between bicontinuous and droplet spinodal decomposition is located around Φ = 0.34, according to our simulation.
coupling between composition and velocity fields [40,41]; however, they turned out to be wrong even on a qualitative level, as will be described later. The difficulty of this problem stems from the non-local nature of diffusion and the many-body and long-range nature of its coupling with hydrodynamic degrees of freedom. Although directional motion of droplets was experimentally reported for such a quasi-two-dimensional situation [26,37], there has so far been no firm experimental evidence supporting such spontaneous motion of droplets for a bulk 3D system. Because of this lack of clear experimental evidence together with the theoretical difficulty mentioned above, phenomenon (ii) has not attracted much attention. Thus, it is now widely believed that droplet coarsening in immiscible liquid mixtures can be explained by the BSS mechanism [1][2][3][4][5][6].
In this Article, we will demonstrate by numerical simulations that, contrary to the BSS mechanism, droplet motion is not random, but hydrodynamically driven by the composition Marangoni force, which is induced by longrange diffusional composition correlation among droplets, indicating that a dynamical coupling between the composition and velocity field leads to a novel mechanism of efficient droplet coarsening in droplet spinodal decomposition.

Results
Phase diagram and the setting of basic parameters Before presenting results of phase-separation simulation, we summarize phase separation behaviour of a binary fluid mixture, focusing on how it depends on the composition φ and the temperature T (see the schematic phase diagram, Fig. 1a). Here we consider phase demixing induced by an instantaneous temperature quench at time t = 0 from an initial temperature in the one-phase region (white region) to a final temperature T in the twophase region. First the phase separation behaviour can be grouped into two types: nucleation&growth (NG)type phase separation in the metastable region (the light yellow region in Fig. 1a) and spinodal-decomposition (SD)-type one in the unstable region. Then, SD-type phase separation can further be classified into droplet and bicontinuous SD. Domain coarsening is mainly driven by the LSW mechanism for NG-type phase separation, which takes place for Φ 1, and by Siggia's mechanism for bicontinuous SD, which takes place for Φ ∼ 1/2. For droplet SD, which takes place for the intermediate Φ (roughly, ∼ 0.21 < Φ <∼ 0.34; the green region in Fig.  1a), it has been widely believed that the BSS mechanism is responsible for domain coarsening. We note that the boundary between NG-and SD-type phase separation is sharp only in the mean-field limit, and should be diffuse for ordinary binary liquid mixtures [42].
The composition region on which we focus here is the region of droplet SD. The thermodynamic state of a binary fluid mixture is characterized by the composition φ and the temperature difference between the critical temperature T c and the temperature T . Then the correlation length of the composition fluctuation ξ, which also characterizes the interface thickness, is described as ξ ∼ = a|(T − T c )/T c | −ν , where the critical exponent ν ∼ 0.63 for 3D Ising universality class, to which a binary mixture belongs [4,14]. On the other hand, the kinetics is characterized by the fluid viscosity η. Then the cooperative diffusion constant is expressed as D ξ = k B T /(6πηξ) and the characteristic lifetime of composition fluctuations is given by τ ξ = ξ 2 /D ξ = 6πηξ 3 /k B T . As shown in Methods, the phase separation behaviour of a binary fluid mixture can be described solely by four dimensionless parameters. One is the scaled compositionφ = φ/φ e , where φ e is the final equilibrium composition, or equivalently the volume fraction of the minority phase Φ. The others are the fluidity parameter A, which is a measure of the relative importance of the hydrodynamic to the diffusional transport [43], the thermal noise parameter B, which is the strength of the noise term (see Methods), and the Reynolds number Re, which characterizes the relative importance of the inertia term to the viscous term. A and B are defined respectively as A = 6π|γ|ξ 3 φ 2 e = 18πσξ 2 2 3/2 kBT , and B = 2 |γ|ξ 3 φ 2 e = 2 5/2 3 kBT σξ 2 [35,43]. We note that σ in these relations is the interface tension in an equilibrium two-liquid coexistence state (see Methods for its expression in the Ginzburg-Landau model). We set Re = 0 since in ordinary droplet phase separation we can safely neglect the inertia effects [24] (see Methods).
We note that near a critical point T c , it is known that the renormalized value of A is a universal constant, A = 4 ∼ 8 [35,43], since the equilibrium interface tension is given by σ = (0.2 ∼ 0.4)k B T /ξ 2 according to the twoscale-factor universality [4]. However, far from a critical point, this parameter can become very large for typical immiscible binary liquids (e.g., for a water/hexane mixture, A ∼ 100), whereas it becomes very small for polymer mixtures since A ∼ N −1 (N : the degree of polymerization of polymer) [4,35]. On the other hand, the estimation of the value of B needs some care, since this represents the strength of thermal fluctuations that directly lead to the break-down of the mean-field picture due to renormalization effects: A larger value of B shifts T c to a lower temperature. We estimate that B should be the order of 1 near a critical point, and it decreases when the temperature is much lower than T c . So, in our simulations, we chose A = 5 and B = 1 as typical parameters for a system near a critical point. We note that the choice of these values of A and B does not affect our basic conclusion on a qualitative level. As shown below, our new mechanism, where hydrodynamic transport plays a key role, becomes more important with an increase in the distance from a critical point, i.e., for larger A and smaller B. Note that A/B is a measure of the ratio of the hydrodynamic transport to the noisedriven diffusional transport. This means that our new mechanism should become more important for a liquid mixture far from a critical point, which is often the case for realistic applications, than for a critical mixture.
Finally, we note that our study concerns only the socalled late stage of phase separation after the formation of a sharp domain interface, where the domain size is large enough compared to the interface thickness: R ξ. Below we study Φ = 0.25, which is in the composition region of droplet SD (see Fig. 1a). Hereafter we use the scaled domain sizeR = R/ξ and the scaled timet = t/τ ξ to describe the coarsening dynamics.

Results without thermal noise
Firstly, we show a coarsening process during droplet phase separation in Fig. 2a, which is simulated by solving the dynamical equations of model H (see Methods) without thermal noise (see Methods). Because of the absence of noise (i.e., B = 0), the droplet interface is smooth. The time evolution of the average size R and the total number of droplets N for the fluidity parameter A = 5 (black) and A = 50 (red) are shown in Fig. 2b and c, respectively. In addition to A = 5, we use the large fluidity parameter A = 50 to access a large separation between the interfacial thickness and the droplet size as well as to see hydrodynamics effects clearly. Because of the composition asymmetry, the symmetry of initial composition fluctuations around the average composition, which holds in the initial linear regime, is broken once the non-linear terms in the free energy start to come into play, and thus many small droplets are formed by fragmentation in the early stage. As stated above, in this study we are concerned only with the late stage of coarsening after the formation of a sharp domain interface. In the late stage, we find R ∝t 1/3 and N ∝t −1 for both fluidity parameters (Fig. 2b,c). Even when thermal noise is absent, we can see spontaneous motion of droplets and subsequent collisions: the droplet coarsening proceeds even without thermal noise. This result can never be explained by the Brownian coagulation mechanism, in which thermal noise is the only cause of droplet motion and coarsening. This is the most direct evidence for the presence of an unknown coarsening mechanism other than the BSS mechanism, which is responsible for the spontaneous motion of droplets and their collisions. Because of the absence of thermal noise, the droplet coarsening process is perfectly deterministic in this case, and thus the initial composition noise introduced at t = 0 completely determines what takes place afterwards.
For A = 5, we can see both collisions of droplets and gradual disappearance of droplets due to evaporation, although the former is much more frequent than the latter. For A = 50, on the other hand, the growth of droplets proceeds almost solely through direct collisions between droplets and evaporation of droplets is very rare. This difference can be explained by the fact that A is a measure of the relative importance of hydrodynamic to diffusional transport. Figures 2e and f show respectively the average diffusion flux and the average magnitude of the velocity of droplets for A = 50 during a coarsening process. Each spike reflects interdroplet collision. We can see that, in the late stage, collisions take place after a complete decay of flow and thus are not affected by the preceding collision even for A = 50, clearly indicating that the migration of droplets and the following collision are due to spontaneous motions of droplets and 'not' induced directly by the preceding collision.

Results with thermal noise
Next we show a coarsening process during droplet phase separation with thermal noise in both composition and force fields. We note that the incorporation of both types of noise in the model H equations has been technically challenging [44][45][46]. However, this has been properly done for 2D in the framework of model H by Camley et al. [44]. We have employed their method to simulate a 3D system (see Methods). The merit of this method is that we can drop the inertia term.
The process of phase demixing is shown in Fig. 3a (see also Supplementary Movie 1). The simulation is made for the fluidity parameter A = 5 and the thermal noise parameter B = 1 (see the above for the choices of these values). We can clearly see the temporal growth of droplets, which are induced by spontaneous droplet motion and the resulting collision and coalescence, as observed in Fig.  2a. Unlike in Fig. 2a, however, the interface of droplets are rough due to thermal fluctuation effects.
As shown in Fig. 3b, the structure factor S(k) can be scaled by using the characteristic wavenumberk 1 (see Methods), indicating that droplet patterns grow self-similarly [1][2][3][4][5][6]. The time evolutions of the average droplet size R and the number of droplets N are shown in Fig. 2c. We find R ∼ t 1/3 and N ∼ t −1 . The exponents are consistent with the prediction of the BSS mechanism, where the coarsening rate is predicted as [4,26]. In our scaled units, this relation becomes R (t) 3 6t. This indicates that the droplet size R grows 1.9 times faster than the prediction of the BSS mechanism although the exponent of the power law is the same. We note that our simulation result is much consistent with the experimental finding of Perrot et al. under microgravity [31] that the droplet size grows twice faster than the prediction of the BSS mechanism. This supports the validity of our simulations with the above choices of A and B.

Spontaneous directional motion of individual droplets
To seek the cause of the non-Brownian nature of droplet coarsening in simulations without noise (Fig. 2) and the faster coarsening than the BSS prediction in simulations with noise ( Fig. 3), we analyse trajectories of droplets during the interval of two successive collisions for the case with thermal noise. First we focus on the nature of droplet motion, i.e., whether it is completely random as assumed in the BSS mechanism (Fig. 4a) or directional (Fig. 4b). We carefully choose the observation period to avoid the effects of flow induced by collisions of droplets and to see only spontaneous motion of the droplets. In Fig. 4c, a typical trajectory of a droplet during a certain time period are shown. We can clearly see that the droplet moves directionally (towards a neighbouring droplet). We confirm that this type of motion is common to most of droplets and not special (see, e.g., Fig. 4d). Below we will show that it is this directional motion that leads to the fast coarsening process. We note that this motion is not due to inertial effects [36]  must be some thermodynamic force acting on droplets that drives droplets to spontaneously move directionally.
For a droplet of radius R which is undergoing random Brownian motion freely, its mean-square displace-  Fig. 3. We note that during this time interval no collisions take place. Symbols with different colours are for different droplets in the simulation box. The blue straight line is R δr ∼ = 0.4t 1/2 (see Eq. (1)) and the red dashed line has the slope of 1. The displacements are always much larger than the prediction of Brownian motion (the blue straight line) for non-evaporating droplets. Note that droplets having displacements smaller than the prediction are evaporating ones.
ment, ∆r 2 , should be given by ∆r(t) 2 = δr(t) 2 = (k B T /(5πηR))t. In our scaled units, for A = 5 and B = 1, In order to compare the displacements of droplets with different radii, we plot R δr for all droplets in Fig. 4d. In this plot, the difference in the droplet size are scaled out; and, thus, the displacements of all droplets should agree to the prediction of Eq. (1) and thus fall on the blue solid line in Fig. 4d, if droplets are doing free random Brownian motion. However, we find that most of droplets have much larger displacements far beyond the prediction of Eq. (1) (or, the blue solid lines), strongly indicating the non-random directional nature of droplet motion. During this time interval of observation, the change in the droplet radius is within a few %, so this effect does not play any role. Hydrodynamic interactions between droplets and finite-size effects should slow down the self-diffusion of droplets. Furthermore, the displacement should still increase in proportion to √ t as long as the motion is random. Nonetheless, Fig. 4d clearly show that the displacements of droplets are much larger than the prediction of Eq. (1) and even increase almost linearly at later times. This cannot be explained by the scenario based on random Brownian motion.
Without noise, we should be able to see clearly the relation between the thermodynamic force and the resulting transport, diffusional and hydrodynamic, without suffering from noise that obscures the details. Thus, we study the chemical potential on droplet interfaces for the case without thermal noise, and the result at a certain time is shown in Fig. 5a. We can clearly see that the chemical potential on the interface of a droplet, or the interfacial composition profile, is not spherical symmetric but rather anisotropic, reflecting the spatial configuration of the neighbouring droplets. As will be shown below, this causes the composition Marangoni force, which induces spontaneous directional motion of droplets. Nevertheless, the droplets are nearly spherical because of the action of the Laplace pressure ∆p = 2σ/R , which induces a quick hydrodynamic relaxation to the stable spherical geometry with a characteristic relaxation time of τ v ∼ ηR/σ. A small variation of the interface tension, ∆σ, in a droplet does not cause shape deformation of the droplet, but causes a drastic change in the type of droplet motion from stochastic to deterministic nature. We note that the former is a relative effect, whereas the latter is an absolute effect. Figure 5b is an enlargement of the blue small droplet and its surrounding droplets shown in Fig. 5a. The arrows show the flow field at their interfaces. The flow field is caused by an imbalance between the Laplace pressure and the thermodynamic force under the incompressibility condition: the break down of the spherical symmetry of the interfacial chemical potential in a droplet and the resulting tangential force acting along the interface allow the directional flow field to be induced even under the condition of ∇·v = 0, unlike the case of an isolated spherical droplet. We can see that the smaller bluish droplet in a higher free energy state is slowly evaporating and the neighbouring droplets in lower energy states are growing. Thus the interface of a neighbouring droplet that is closer to the evaporating small droplet has smaller chemical potential (more bluish) than the opposite side. This anisotropy of the chemical potential on the interface of each droplet is the origin of the deterministic directional motion of the droplet. The variation of the chemical potential on a droplet interface is less than a few %, but the resulting composition Marangoni force is strong enough to overwhelm the thermal force acting on a droplet.
The emerging physical picture is as follows. An small (evaporating) droplet causes deterministic flow a b FIG. 5: Coupling between composition and velocity fields. a, A snapshot of droplets during phase separation without thermal noise. The colour reflects differences of the chemical potential on the interface. Reddish and bluish colour means higher and lower chemical potential, respectively. As shown in the colour bar, the colour is assigned for each droplet, which is bound by the minimum and maximum value of the chemical potential µ. This is because the amplitude of the change of µ inside each droplet interface is much smaller than the fluctuations of the absolute value of µ of each droplet. We note that the magnitude of difference relative to mean value is less than 3% even for a droplet which has the largest intradroplet variation. b, Enlargement of the blue small droplet in a and its surrounding droplets. The meaning of colour is the same as in a. The velocity field at the droplet interfaces is also shown by the arrows. A longer arrow means high velocity.
fields on its neighbouring larger droplets, which lead to motion of the neighbouring droplets towards it and eventually induce an interdroplet collision, if the droplet does not evaporate before collision. In other words, a larger droplet tends to eat a neighbouring smaller droplet by direct collision. This hydrodynamic eating speed is much faster than the diffusional one in the evaporation-condensation (LSW) mechanism, as long as hydrodynamic transport is faster than diffusional transport. More precisely, however, the composition correlation between neighbouring droplets induces both the hydrodynamic and diffusional transport. Thus, the crossover between the present mechanism and the evaporation-condensation mechanism should be controlled by the competition between hydrodynamic and diffusional transport. The former is always dominant as long as Φ is not so low. For lower Φ, however, the direct compositional correlation between droplets, which is the origin of the thermodynamic force inducing hydrodynamic droplet motion, becomes weaker due to the less overlap of the composition fields around droplets. For very low Φ, thus, droplets interact with each other via the mean-field matrix, as assumed in the original mean-field theory of the LSW mechanism [21,22]. In this regime, the LSW mechanism becomes dominant. We speculate that the crossover takes place around Φ ∼ 0.21 separating NG and SD (see Fig. 1) since by definition droplets are formed rather independently for NG whereas in a correlated manner for SD. But this crossover composition may also depend on the value of A. Although this is a very interesting problem, we leave it for future investigation.

Analysis of a pair of droplets
To elucidate the exact physical mechanism, we also study the motion of pairs of droplets (see Methods). Typical results are shown in Fig. 6a and b, where the chemical potential and the flow field on the interfaces are shown as in Fig. 5b. Figure 6c and d show the composition profile of two droplets along the line connecting the centres of mass of the two droplets for case a and b. Figure 6e and f are their enlargements, which clearly show the presence of the composition difference by ∆φ between the two sides of a droplet. We can see that such an asymmetric interfacial profile of a droplet induces a chemical potential gradient in the droplet and leads to spontaneous directional droplet motion. Finally, the temporal change in the centre-of-mass velocities of the two droplets are shown in Fig. 6g and h, respectively, for case a and b.
For the case of Fig. 6a, we calculate the interfacial tension σ for each interface by integrating F φ = −∇ 2φ ∇φ [40] across the interface along the radial direction. A part of this thermodynamic force is balanced by the Laplace pressure ∆p = 2σ/R. For the left droplet whose radius is 11.5 and velocity is -0.0045, the interfacial tension is 0.9286 for its left side and 0.9296 for the right side. This small variation of the interfacial tension is consistent with a nearly spherical shape of droplets, as discussed above. It is the interfacial tension gradient due to the composition gradient on the interface that leads to the motion of a droplet. This is known as the composition Marangoni effect [47,48]. We note that it is well known that this effect induces spontaneous droplet motion for droplets with weak surfactants [49] and those subjected to composition or temperature gradients [50,51]. We stress that unlike these well-known cases, for the present case the interfacial tension gradient in a droplet is produced by a non-trivial coupling between the composition fields around droplets. The analytical expression for the droplet velocity under a uniform composition gradient is give as V = −(R/(5η))(dσ/dx) [47]. By inserting to this relation the above values of σ estimated from the simulation, we estimate the droplet velocity as |V| ∼ 0.0050, which is much consistent with the directly measured one (|V| ∼ 0.0045) (see Fig. 6g). A slight difference may come from the difference in the geometry, i.e., a single particle under a uniform composition gradient (theory) vs. a droplet pair under diffusional coupling (our simulation), and the finite size of the simulation box. To avoid the finite-size effects, we need to perform larger size simulations, but unfortunately a large numerical cost of hydrodynamic simulations do not allow us to do so. The velocities of the droplets in Fig. 7b are also estimated from the same analysis, and they also agree well with the measured ones (see Fig. 7b and its caption). These results strongly support our mechanism. In the particular configuration in Fig. 6a, due to the perfect symmetry around the centre of mass of the two droplets and the absence of noise, coarsening does not proceed. We note that the same situation should be realized for the LSW mechanism. This is not the case for the asymmetric situation in Fig. 6b, where two particles approach each other and eventually collide. For real droplet phase separation, the randomness of droplet size and many-body interactions always lead to such asymmetric situations and their continuous presence is responsible for continuous coarsening. We emphasize that the randomness is a direct consequence of thermal fluctuations and thus intrinsic to any phase separation phenomena. Combining this result on pairs of droplets and the results shown in Fig. 5a and b, we can conclude that the composition Marangoni effect stemming from the randomness of droplet size is responsible for spontaneous directional motion of droplets and the resulting coalescence. The importance of randomness indicates that there should be a non-trivial feedback between the coarsening and the droplet size distribution. We find that the normalized droplet size distributions at various phaseseparation times can be superimposed on the universal curve (see Fig. 7), consistent with the self-similar nature of pattern evolution. We speculate that more frequent collisions between particles with large size difference may be responsible for the self-similarity, but this problem should be studied carefully in the future. Our mechanism cannot lead to domain coarsening for a monodisperse emulsion. However, the Brownian coagulation mechanism should work even in such a situation; and once a collision happens and produces the droplet size difference, our mechanism should become operative.

Roles of other non-trivial mechanisms
Finally we briefly consider other possibilities. In our previous papers [26,37], we suggested two other nontrivial droplet collision mechanisms: multiple-collision and collision-induced collision mechanisms, which was confirmed experimentally and numerically [38,39,52]. Here we show that neither of them can explain the abovedescribed phenomena. The former mechanism is a rather obvious one. When two droplets collide with each other, the shape of the droplet interface is strongly deformed and relaxes to a spherical shape, accompanied by the interfacial-tension driven flow. This flow induces hydro- We can see that the normalized droplet size distribution functions at various phase-separation times can be scaled onto a universal curve at least approximately, consistent with the self-similar nature of the domain growth (see Fig. 3b).
dynamic motion of the surrounding droplets, which leads to subsequent collisions. The latter mechanism is, on the other hand, a result of a coupling between droplet collision and composition diffusion. After droplet collision, strong diffusion flux towards the resulting merged droplet is created due to what we call the interface quench effect [43] (see Fig. 2d), which is the process of local equilibration of the droplet with the matrix. This strong diffusion affects the translational motion of other neighbouring droplets due to the dynamical coupling of the diffusion and flow fields, which leads to subsequent collisions. Both mechanisms are expected to enhance successive collisions, but are effective only immediately after previous collisions. To reveal whether the coarsening behaviour can be described by these mechanisms or not, we calculate the average diffusion flux and the average magnitude of the velocity during a coarsening process without noise respectively in Fig. 2d and e. Occasional interdroplet collisions accompany sharp spikes and subsequent relaxation for both quantities. The characteristic decay times of diffusion flux and flow are given by τ φ = R 2 /D 0 = 6πηR 2 a/(k B T ) and τ v = ηR/σ, respectively. As can be seen in Fig. 2d and e, in the late stage a collision takes place after the diffusional and hydrodynamic fluxes induced by the preceding collision have completely decayed, suggesting that it is not affected directly by the preceding collision: there are no chain collisions. This strongly indicates that the coarsening behaviour we observe in the late stage of phase separation can be explained by neither of the above-mentioned mechanisms.

Discussion
The dynamical scaling argument on domain coarsening in the late stage of phase separation has been so successful [1][2][3][4][5][6]. In this argument, it has been assumed that the compositions of the two phases already reach the final equilibrium one φ e in the late stage. For droplet phase separation, thus, it is implicitly assumed that the interfacial tension σ is uniform over a droplet's interface. In liquid mixtures, which are incompressible in ordinary conditions, the hydrodynamic flow can be induced only by the interfacial force, which acts normal to the interface in the above assumption. For a spherical droplet, this leads to the conclusion that the spherical symmetric force fields which act along the interface normal unit vector n of a droplet cannot cause any flow under the constraint of ∇ · v = 0 and accordingly the interfacial force should be balanced by the Laplace pressure ∆p = 2σ/R (σ being the interfacial tension) [20]. This is the basis for the long belief that droplet motion is not induced hydrodynamically but solely by its thermal diffusion.
In the above, we have shown that this physical picture is not correct and the interfacial tension is actually inhomogeneous over each droplet's surface due to diffusional composition correlation among droplets, which allows interfacial-tension-driven flow even for a spherical droplet. The key is the fact that droplets exchange the component molecules (or atoms) with neighbouring droplets via the surrounding matrix phase by diffusion. The time necessary for a component molecule (or atom) (size a) to diffuse over the same distance is τ s = L 2 /D 0 . On the other hand, the average time necessary for a droplet of radius R to diffuse over an average inter-droplet distance L is estimated as τ D = L 2 /D R , where D R is the diffusion constant of a droplet. Since D 0 and D R are inversely proportional to the size of a component molecule (or atom), a, and the domain size, R, respectively, the ratio τ D /τ s becomes very large in the late stage of droplet SD, where R a. Thus the composition correlation between neighbouring droplets develops via molecular (or atomic) diffusion long before an accidental collision between them by thermal Brownian motion takes place. It is this composition correlation among neighbouring droplets that is responsible for internal inhomogeneity of the interfacial tension over a droplet's interface and the resulting Marangoni force.
Here we note [24,26] that van der Waals attraction between droplets is too weak to affect the motion of droplets in the late stage. Furthermore, directional motion of droplets observed in our previous experiments in a quasi-2D situation, where neither gravitational force nor electrostatic interactions are relevant, can be explained by the present mechanism.
Here it is worth mentioning that the mechanism we find here is essentially different from mechanisms previously proposed by one of the authors (H.T.) [40] and Kumaran [41]. In the above, we reveal that the force acting on a droplet comes from the non-uniformity in the amount of the composition jump across the interface in the droplet and thus is localized on the interface, whose characteristic width is given by ξ. This non-uniformity of the interfacial tension due to the coupling of the composition fields around droplets causes the composition Marangoni flow, which leads to the spontaneous motion of droplets. In both of the previous theories, however, the non-uniformity of the Gibbs-Thomson condition on the droplet interface was completely ignored and only the composition gradient in the matrix phase was considered. The characteristic lengthscale of this inhomogeneous composition is the order of R. Thus, the strength of the composition Marangoni force in the new mechanism is much stronger than the composition-gradientinduced force in the previous ones since 1/ξ 1/R. Furthermore, there is another important, even qualitative, difference between them; for example, both of the previous theories [40,41] (wrongly) predict that for a pair of large (growing) droplets or for a pair of small (shrinking) droplets, the two droplets should come towards each other, whereas for a pair of large (growing) and small (shrinking) droplets the two droplets should move away from each other. We stress that this is opposite to what is shown in Fig. 6a and b.
The novel mechanism reported here can be explained on an intuitive level as follows. There is a difference in the local free energy between droplets, reflecting the difference in the size between them. Larger droplets are more stable and in a lower free energy state simply because the interface energy cost is smaller. This leads to the diffusion flux from small droplets to nearby larger droplets due to the chemical potential difference. This correlation of the composition field between droplets evolves much faster than any significant diffusional centre-of-mass motion of the droplets. This is exactly the same process as that of the LSW mechanism: if there are neither hydrodynamic degrees of freedom nor translational droplet diffusion, this leads to the coarsening of the LSW type. However, the presence of hydrodynamic degrees of freedom in a system entirely changes the kinetic route for lowering the free energy. An intradroplet gradient of the interfacial tension generates hydrodynamic motion of droplets towards the side with lower interfacial tension. Thus, larger droplets tend to eat nearby small droplets, as in the case of the LSW mechanism. However, the crucial difference comes from the fact that the transport process is dominantly controlled by hydrodynamic translational motion of droplets due to the composition Marangoni effect in our mechanism, whereas by translational diffusion of molecules (or atoms) in the LSW mechanism.
In a more general perspective, we may say that domain coarsening is under the influence of the free energy at any moment. This has a significant implication for the unified physical picture of domain coarsening: all coarsening mechanisms, including the evaporation-condensation mechanism (for Φ 1), the hydrodynamic coarsening mechanism ( for Φ ∼ 1/2), and the mechanism found here (for intermediate Φ), have a common principle of material transport: Transport can be either diffusion, flow, or by their coupling, but always takes place from smaller to larger domains: Domain coarsening obeys a rule that "big always wins over small". We note that there is no such a rule for the BSS mechanism, in which the process is purely random and stochastic.
Finally, we stress that our new mechanism may play an important role in phase separation of many binary mixtures in practical use, which takes place in a state far from a critical point. We note that our mechanism should be more operative for a mixture of lower viscosity, as mentioned earlier. We also stress that our mechanism can operate even near zero temperature, as shown in Fig.  2, unlike the Brownian-coagulation mechanism; and thus it should be relevant to spinodal decomposition in quantum systems, which often takes place at extremely low temperatures [9][10][11]. We note that in such a condition, quantum diffusion associated with tunnelling mass transfer can be operative but there is little thermal motion of droplets. Our study shows that the coupling between the composition and velocity fields alone can lead to domain coarsening even with little thermal noise. The relative importance of the velocity field to the composition field monotonically increases with an increase in the volume symmetry between the two phases; and thus the primary mechanism of coarsening changes from the LSW, to our mechanism, and to the Siggia's mechanism with an increase of Φ from 0 to 1/2. We hope that our mechanism will provide a novel perspective for spinodal decomposition of classical and quantum fluid systems in various fields.

Methods
Basic dynamic equations of model H. For our simulations, we employ the following kinetic equations known as a fluid model, or model H [1,4,14]: where φ is the composition, v is the fluid velocity, p is a part of pressure, ρ is the density, η is the viscosity, L is the kinetic coefficient, and β = 1/k B T (k B : the Boltzmann constant). The pressure p is determined to satisfy the incompressiblity condition ∇ · v = 0. Since our primary concern is the dynamical coupling between the composition and flow field, we do not care about the detailed properties of the mixture and thus adopt the following standard Ginzburg-Landau free energy functional: where γ = γ 0 (T − T c ) and γ 0 , u and K are positive constants. In Eq. (3), F φ is the thermodynamic force density acting on the fluid due to spatial inhomogeneity of the composition field and F φ = −φ∇µ = −∇π+k B T Kφ∇ 2 ∇φ (π: osmotic compressibility), where µ = δH/δφ is the chemical potential. Here θ and ζ are the thermal noise terms in composition and force, which satisfy the following fluctuation-dissipation relation, respectively [1,4]. and Here we rewrite Eqs. (2) and (3) by scaling the space and time unit by the correlation length of critical composition fluctuations ξ = (K/|γ|) 1/2 and the characteristic lifetime of fluctuations τ ξ = ξ 2 /(|γ|L). Then we define new variables,r = r/ξ (r being the position vector), t = t/τ ξ andṽ = (ξ/τ )v. The composition is normalized asφ = φ/φ e , where φ e = (|γ|/u) 1/2 is the equilibrium composition of the coexisting phase. The scaled equation which corresponds to Eq. (2) is then and Eq. (3) with the Stokes approximation becomes where T ij is the Oseen tensor given by θ andζ are scaled thermal noises which satisfy the following fluctuation-dissipation relations: θ (r,t)θ(r ,t ) = −B∇ 2 δ(r −r )δ(t −t ), and ζ i (r,t)ζ j (r ,t ) = −(B/A)δ ij ∇ 2 δ(r−r )δ(t−t ). (11) Here we note that σ/(k B T ) = (2 3/2 /3)|γ|ξφ 2 e , which is the interface tension in an equilibrium two-liquid coexistence state.
Simulations without thermal noise. We performed simulations without thermal noise to show that droplet coarsening takes place even without thermal noise (i.e., B = 0). We note that the Brownain coagulation mechanism fully relies on the presence of thermal force noise. Simulations without thermal noise are also useful to see the spatial distribution of physical quantities such as the chemical potential and the fluid velocity without suffering from thermal noise. For these simulations, we introduce Gaussian random noise intoφ with an intensity δφ = 0.01 around the averaged compositionφ 0 = −0.5 (i.e., Φ = 0.25) as an initial condition only att = 0, switch off noise fort > 0, and follow the pattern evolution process.
Numerical scheme to solve the above set of equations properly with both composition and force noises and the details of simulations. Here we explain how thermal composition and velocity fluctuations can be incorporated in our simulations. Although we have written Eqs. (7) and (8) as two separate equations, they can be unified into only one kinetic equation in the Stokes regime, since the velocity field is set solely by the composition field. Substituting Eq. (8) into Eq. (7), the dynamic equation of composition is given by [4,35] ∂φ(r,t) ∂t = dr L(r,r ) δH δφ(r ) + θ R (r,t), where the kinetic coefficient is given as L(r,r ) = ∇ 2 δ(r −r ) + ∇φ(r) · T(r −r ) · ∇φ(r ) and the random noise term θ R (r) = −v R (r) · ∇φ +θ satisfies the following fluctuation-dissipation relation θ R (r,t)θ R (r ,t ) = 2L(r,r )δ(t −t ).
Here the random velocity noise v R is given by v R (r,t) = dr T(r −r )ζ(r). Since the kinetic coefficient of the above equation depends upon the composition fieldφ, the noise term is multiplicative and should be treated via the Stratonovich interpretation [53]. Thus we adopted a Stratonovich scheme with semi-implicit terms developed by Camley and Brown [44] under a periodic boundary condition with a staggered grid. This scheme was originally conceived to study the phase separation process in a two-dimensional membrane. We modify their approach by using the Oseen tensor appropriate for a threedimensional (3D) Newtonian fluid in the Stokes regime.
Here it is worth noting that 3D simulations of model H with thermal noises were also performed by a Lattice-Boltzmann method [46]. The system size was set to 256 3 lattices, which is large enough to study the late stage of phase separation. We set the grid size to be ∆x = ∆y = ∆z = 1. The volume fraction of the minority phase Φ is set to 0.25 and the composition is set to spatially homogeneous att = 0.
Simulations of the motion of a pair of droplets. We place two droplets whose interfacial profile is described by a step function att = 0 and then follow the process of equilibration by solving the model H equations. In this simulation, the fluidity parameter A is 50 and the thermal fluctuations are not included (i.e., B = 0). Here the large fluidity is used to access a large separation between the interfacial thickness and the droplet size within a limited simulation time as well as to see hydrodynamics effects clearly. Note that the process is much slower without thermal noises. The simulation box size is 192 3 .
Analysis. We analyse the sizes and positions of droplets by counting the number of lattice points belonging to each droplet and calculating the positions of the centre of mass of the lattice points after binarization of the composition field. We also obtain the velocity of a droplet by calculating the average of velocity over lattice points inside the droplet. The structure factor S(k) is obtained by spherically averaging the power spectrum of the Fourier transformation ofφ(r). The characteristic wavenumber k 1 is calculated ask 1 = dkkS(k)/ dk S(k).