Investigation of thermocapillary migration of nanodroplets using molecular dynamics

Molecular dynamics is used to investigate the thermocapillary motion of a water nanodroplet suspended in benzene subjected to a constant temperature gradient. This framework lets us identify the average behavior of the ﬂuid particles by revealing their mean evolution. We connect such statistics to the behavior of the temporally evolving nanodroplet, thereby providing a microphysical foundation to existing macroscopic models that rely on the assumption of continuum. It is shown that, despite the signiﬁcant Brownian effects, the droplet exhibits the macrophysical expected behavior, i


I. INTRODUCTION
Surface-tension-driven flows, also called Marangoni convection after the Italian physicist Carlo Giuseppe Matteo Marangoni, 1 refer to convective effects generated by a nonuniform distribution of surface tension along the boundary separating two fluids (e.g., a liquid and a gas or two immiscible or partially miscible liquids, see, e.g., Gaponenko et al. 2 and Shevtsova et al. 3 ).This nonuniformity can be caused by different possible root causes, such as active compounds at the interface (e.g., surfactants 4,5 ), solvent evaporation (leading to a variation of solute concentration at the interface), 6 or temperature gradients. 7,8n such context, the study of droplet-based systems represents a subject of great interest given the remarkable connections that it has with a variety of technological fields.Relevant examples range from the synthesis of metal alloys (where droplets are typically present as an immiscible or partially miscible liquid phase dispersed in an external liquid [9][10][11][12] ), microfluidics, [13][14][15][16][17][18] manipulation of particles, 19,20 combustion, 21 and nanosystems 22 to the production of drugs and medicines (droplets being employed as microreactors for the nucleation of proteins or other macromolecular substances) and other relevant biotechnological applications (where droplets are generally used to encapsulate biological entities into synthesized particles [23][24][25][26][27][28][29][30][31] ).
Droplet-based systems have enjoyed widespread attention over recent years owing to another important property.In some circumstances (in situations where the above-mentioned Marangoni effect plays a significant role) they can be turned into self-propelled entities.The spontaneous tendency of droplets to behave as a kind of machinery, potentially usable to transport and deposit quantities of interest in different points of a liquid matrix, is indeed a well-known phenomenon.A vast amount of literature exists where it has been shown that such dynamics can be produced as a result of the rupture of translational-invariance in space.In turn, such symmetry can be broken either as a result of a bifurcation, i.e., a process able to make an initially isotropic state unstable, 32,33 or due to the presence of an imposed ("external") gradient of concentration or temperature.
The latter case (imposed difference of temperature) has attracted much attention due to its remarkable property of enabling the "vectorial" (i.e., unidirectional) migration of droplets along a welldefined direction (droplets typically move from a cold area toward a warmer region along the direction of the prevailing temperature gradient).A first milestone work on which several other analyses would successively rely on was the analysis by Young, Goldstein, and Block. 7hese authors demonstrated experimentally that a gas bubble suspended in a heavier fluid (a liquid) could even migrate against buoyancy if the temperature difference applied to the host fluid (liquid) was properly tuned.The same authors derived an analytical solution of the problem in the limiting conditions of negligible convective transport of energy and momentum.
Even though, historically, the first analytical solution of the problem under consideration has been attributed to Young, Goldstein, and Block, 7 in more recent years, work from the Soviet scientist A. I. Fedosov, 34 originally published in 1956 in Zhurnal Fizicheskoi Khimii, has emerged.In his work, Fedosov 34 tackled the thermocapillary motion of a drop considering a zero-gravity environment using a perturbation expansion approach and obtained an expression for the droplet velocity.It was only a few years later that Young, Goldstein, and Block 7 derived an expression for the drop migration velocity independently.It is worth highlighting that contrarily to Fedosov's work, 34 in the work of Young, Goldstein, and Block 7 gravity was considered and experiments were performed.In the following, we will refer to this theory as F-YGB theory, to give credit to both works.
With the advent of orbiting platforms, other experimental studies succeeded in revealing the intimate nature of this process by filtering out the influence of buoyancy forces.][40] Given the limited opportunities for (and high costs being associated with) microgravity research, however, most of the existing knowledge on this subject has been produced by resorting to alternate strategies heavily relying on mathematical models and proper techniques for the solution of the related equations.As the number of these studies is extremely large, however, in the following, we will content ourselves with giving only some flavor of the most important findings.
In this regard, it is worth starting from the simple remark that, using the landmark analytical solutions of Young, Goldstein, and Block 7 as a basis (determined in the limit of vanishing Reynolds and Marangoni number and a perfectly spherical drop), other authors applied asymptotic-expansion methods to account for new effects (formally introducing them as small perturbations superimposed on the basic creeping flow).Relevant successful attempts along these lines are those by Subramanian 41 who took into account the convective transport of energy for Ma < 1, and Balasubramaniam and Chai, 42 Haj-Hariri, Nadim, and Borhan, 43 and Ford and Nadim, 44 who addressed the opposite circumstances where the convective transport of heat is negligible (Ma !0) while small inertial effects are retained in the momentum equation (i.e., for finite Re).Interestingly, Haj-Hariri, Nadim, and Borhan 43 and Ford and Nadim 44 also calculated the correction to the migration velocity caused by the shape deformation.Some year later, Balasubramaniam and Subramanian 45 analyzed the opposite asymptotic limit in which Re ! 1 and Ma ! 1.For this case, the migration velocity of the drop was obtained in the framework of a potential-flow theory by equating the rate at which work is done by the thermocapillary stress to the rate of viscous dissipation of energy.Under these assumptions, these authors found the velocity of an isolated drop to scale with the square of the temperature gradient and the cube of the radius of the drop (as opposed to the F-YGB's formula valid in the limit Re !0 and Ma !0, where the corresponding dependencies are linear).
Continuing with the development of a historical perspective on the subject, it should also be pointed out that typical difficulties relating to the analytical solution of the Navier-Stokes equations (and the concurrent increase in the performances of available computer hardware) have spurred over the years an increasing use of computational procedures.Among other things, these computational approaches have proven to have the potential of allowing insight into the processes that control droplet migration with no need to introduce limiting assumptions or simplifications (like those typically needed by the analytical approaches outlined before).This explains why understanding of the flow physics underlying these problems has relied to an unusual level on insights gained through methods such as the finite-difference Volume of Fluid (VOF) or level-set techniques, [46][47][48][49][50][51] or the lattice Boltzmann method. 52,53The level of sophistication in how the interface is treated generally underpins the computational cost of such numerical schemes (the reader interested in a more thorough discussion should consult the relevant references cited above).Here we limit ourselves to stating that, by providing an alternative cost-effective means of simulating real processes, CFD (computational fluid dynamics) has extensively been used as a relevant means to test theoretical advances for conditions often unavailable experimentally or having a prohibitive cost.
5][56][57][58] In principle, they can even be used to perform a variety of "digital fluidic" operations, which can be made programable and reconfigurable. 23][61][62][63] Despite this valuable way of thinking and the solid knowledge that has been developed in terms of droplet "deterministic behavior" (dynamics at the macroscopic scale as a result of the application of well-defined stimuli), however, the early promise of nanodroplets as advanced nanoengineered carriers has been tempered by a lack of consensus and clarity in identifying the relative importance at microscale of surface-tension induced and other (nondeterministic) effects.
Put simply, further developments in these fields have been hindered primarily because of our limited understanding of how these systems behave at submicron scale.With the descend of scale, other influential factors start to be dominating compared to purely surfacetension induced convective effects and complex interplay of these together contributes to significant challenges in predicting the effective dynamics at various scales.
This begs the question: how do thermocapillary flows behave in such miniaturized systems where viscous effects and Brownian dynamics dominate?
To the best of our knowledge, the only study combining Marangoni flows with a nanofluidic scenario is an investigation by Sumith and Maroo, 64 in which a thin film of liquid argon was placed in a vacuum, next to a solid wall.A hot-spot was then included at the center of the solid plain, creating a temperature gradient.The results indicated the presence of surface-tension-driven flows, as the argon started forming circulating regions.
In a similar way, in the present paper, molecular dynamics (MD) simulations are used as a viable approach to support the microscopic interpretation of the Marangoni effect in relation to the aforementioned droplet vectorial migration phenomenon.In particular, we focus on the motion of a water nanodroplet suspended in benzene, and confined between two solid, monatomic walls.The droplet velocity is measured for different temperature gradients and is compared against available (macroscopic) analytical models.This is made with the twofold intent to determine on the one hand the peculiarities of this phenomenon on a very small scale and, on the other hand, to support the elaboration of relevant microphysical reasoning through critical comparison with the information that can be obtained using macroscopic frameworks.
Particularly in the case of water, there are several forcefields that can be used for MD simulations.Water models can be categorized based on the number of point charges (sites) used to model the molecule.Three-site models, such as the simple point charge (SPC), 65 the extended SPC (SPC/E), 66 and the transferable intermolecular potential with 3 points (TIP3P), as well as the classic four-site model, TIP4P, have been studied for over four decades now.It is well established that each model has its strengths and weaknesses, and there is no one model that surpasses all others at everything. 67Subsequently, the fivesite model TIP5P was developed 68,69 which generally resulted in higher accuracy, at the expense, however, of computational efficiency.Later, a four-site model, named TIP4P/2005, 70 attempted to match the accuracy of TIP5P, while maintaining the computational cost of four-site models.A number of investigations comparing different water models favored the TIP4P/2005 model, particularly in fluid properties relevant to this study, such as viscosity, and diffusivity. 71Thus for this study, we chose the TIP4P/2005 water model.
The benzene molecules were modeled using the optimized potentials for liquid simulations model-All atom (OPLSAA) model, 72 which was specifically designed for organic molecules in the liquid phase and has been used in MD simulations studying similar organic molecules. 73 significant part of the present investigation is dedicated to the characterization of the thermal conductivity, viscosity, and surface tension of water (using the TIP4P/2005) and benzene (using OPLSAA).
Although not the primary objective of the paper, this practice is intentionally used to assess the strengths and weaknesses of these computational models in capturing the actual physical properties of the considered fluids.

II. PROBLEM STATEMENT
The system under investigation is constituted by a small droplet suspended in an otherwise quiescent fluid in which a temperature gradient, r 1 T, is maintained by external means, as indicated in Fig. 1.It is assumed that, as a result of the ensuing gradient of interfacial tension (the surface tension being a monotonic decreasing function of the temperature), the thermally induced capillary stresses generate flow inside and around the droplet.In such circumstances, droplet displacement is produced to counter-balance, through viscous stresses, the mismatch of the capillary stresses at the interface separating the two immiscible fluids.As a result, the droplet migrates following the direction of the imposed temperature gradient.
From a macroscopic point of view this balance of forces per unit area acting at the interface reads 74 where n is the unit vector normal to the droplet interface (directed toward the external liquid), I is the second-order unit tensor, and c T ¼ À@c=@T is the rate of change of surface tension with temperature.Quantities indicated with/without a tilde are referring to the droplet/surrounding phase.Moreover, s d is the so-called macroscopic viscous stress tensor related to the overall stress tensor by the relationship s ¼ ÀpI þ s d .Neglecting intermolecular forces and taking into account only forces that represent the diffusion of fluid particles in space (due to random motion) and/or result from the effective collision of molecules, s can formally be expressed as 75 where m i is the mass of a generic molecule, dðrÞ is the Dirac function, C i represents the aforementioned random motion of molecules, C i C i indicates a dyadic product, r i the molecule position vector, and the symbol hqi is used to indicate formally an "average" value from a stochastic standpoint (practically obtained by integrating the generic quantity q multiplied by a proper distribution function 76 ).Though this representation of the stress tensor should be regarded as an approximation of reality (it does not take into account stresses due to intermolecular forces), given its simplicity, it is used in the present work to support some of the explanations elaborated in the results section (with regard to the role played by Brownian effects).For Newtonian homogeneous and isotropic fluids, the macroscopic relationship emerging from it [via application of the distribution function and ensuing space integration of Eq. ( 2)] has proven to provide very good agreement with empirical observations.It can formally be cast in condensed form as s d ¼ 2lð$VÞ s , where l is the fluid dynamic viscosity and V is the macroscopic fluid velocity (being related to the total velocity of each molecule or atom v i by the rela- Therefore, the macroscopic stress balance condition (valid on the level of continuum) can be rewritten as In line with previous works (cf.Capobianchi, Lappa, and Oliveira, 8 Capobianchi et al. 77 and references therein), we define the Reynolds number as, Re ¼ qU T a=l, and the Marangoni number, Ma ¼ U T a=a, where U T is a characteristic velocity of the thermocapillary flow, a is the drop radius, q is the density, and a is the thermal diffusivity.The characteristic velocity is defined as U T ¼ c T r 1 Ta=l, where r 1 T is the modulus of the imposed temperature gradient.
The Marangoni number provides a measure of the relative importance of convective transport of energy with respect to its diffusive counterpart.In the limiting condition Ma !0, therefore, thermal diffusion is preponderant with respect to convection.This simplification allows a straightforward solution of the energy equation.Similarly, in the limit as Re !0, convective transport of momentum is also negligible which greatly simplify the solution of the Navier-Stokes equations (both simplifications were originally used by Young, Goldstein, and Block 7 to determine an analytical expression for the migration velocity, the reader being referred to Sec.IV B below for additional details).
In the present context, we consider a water droplet suspended in benzene, confined by two Lennard-Jones (LJ) solid walls that are perpendicular to the z direction.The simulation box has dimensions L x ¼ L y ¼ 60:5 A ˚, and L z ¼ 228 A ˚, in the x, y, and z directions, respectively (see Fig. 2).Moreover, the characteristic dimension of the water droplet is a ¼ Oð10 À9 Þm.The applied temperature gradient is r 1 T ¼ Oð10 10 Þ K/m while the variation of surface tension with temperature and viscosity of the surrounding phase determined from the MD simulations (cf.Sec.III B) are c T ¼ Oð10 À4 Þ N/(mK) and l ¼ Oð10 À3 Þ Pa s, respectively.Finally, the kinematic viscosity and the thermal diffusivity are ¼ Oð10 À6 Þm 2 =s; a ¼ Oð10 À7 Þm 2 =s (we recall that the Prandtl number for water is, Pr ¼ =a % 10, showing the good accuracy of the present numerical modeling).Using these values we find Re ¼ Oð10 À3 Þ and Ma ¼ Oð10 À2 Þ, indicating that convective transport of both momentum and energy can be considered negligible (though not totally absent).
We wish to remark that in these conditions buoyant effects can also be considered negligible.This could easily be verified by evaluating another typical characteristic parameters for these problems, namely, the so-called dynamic Bond number.This number, which measures the relative importance of buoyancy with respect to the surface tension force, can be defined as 78 G ¼ gaðq À qÞ=3c T r 1 T, where g is the acceleration due to gravity.Given the extremely small radius of the drop, and the comparable densities of the two considered liquids, the value taken by this number for the conditions examined in the present work is G ¼ Oð10 À13 Þ, which confirms the validity of our assumption.
Before moving to Sec.III, we wish to clarify an important aspect regarding the exceptionally large temperature gradient adopted in our simulations.To this end, evaluating a characteristic convective timescale as t c ¼ a=U T is beneficial.We wish also to remark that in our simulations we do not have control over the material properties, which are naturally determined by the MD model once a fluid pair is assigned.Considering a realistic temperature gradient, for instance, rT ¼ Oð10 3 Þ Km À1 , together with the values of c T , a and reported above, and observing that q ¼ Oð10 3 Þ kgm À3 , the above-mentioned timescale would be t ¼ Oð10 À2 Þ s.Following up on the previous point, it becomes clear that simulating the problem using MD for a time as large as this value, would be impracticable (note that our simulations have resolved the flow for a time which was on the order of nanoseconds and each run lasted for several weeks).In order to compare with the continuum F-YGB theory, we have paid more attention to the characteristic numbers which are obtained when all these parameters are combined to obtain nondimensional groups, namely, the Marangoni and Reynolds numbers.Our Ma and Re based on a temperature gradient as large as rT ¼ Oð10 11 Þ Km À1 are sufficiently small to assume that the continuum model remains applicable within a reasonable approximation.Thus, our decision to consider such a large temperature gradient should be regarded as a practical compromise between the need to to make the computational effort reasonably contained and the additional requirement of maintaining "sufficiently" small Marangoni and Reynolds numbers to obtain results that align with the F-YGB theory.

III. METHODOLOGY
This section describes the molecular models used for the thermocapillary migration of the droplet, as well as those used to calculate the thermal conductivity, viscosity, and interfacial tension of benzene and water.
The scripts for creating the initial molecular configuration, as well as example scripts for each type of MD simulations, are publicly available at https://github.com/michael-frank-research/marangoni_paper.All our models use open-source software, and the interested reader should be able to fully reproduce our results, the only potential constraint being the computational resources [the simulations for the droplet migration will require high-performance computing (HPC) facilities].

A. Molecular model
The initial molecular configuration was prepared with the help of the molecular builder "moltemplate." 79ach wall is made out of two atomic planes, with the atoms placed on the sites of a face-centered cubic (FCC) lattice.The lattice parameter of the FCC lattice is 4.05 A ˚, corresponding to the density of aluminum, q w ¼ 2700 kg/m 3 .
The benzene molecules are initially arranged on a cubic lattice with lattice parameters a ¼ 5:9, b ¼ 5:49, and c ¼ 4:39 A ˚in the x, y, and z directions, respectively.The initial configuration is relatively arbitrary and attempts to realize the density of benzene, q B ¼ 880 kg/ m 3 , while limiting any overlap between molecules that often causes the solution to blow up. 81imilarly to benzene, the atoms of the water droplet are also initially arranged on a primitive cubic lattice with lattice constants a ¼ b ¼ c ¼ 3:45 A ˚in the x, y, and z directions, respectively.The center of the cube is placed at the center of the xy plane (i.e., x ¼ y ¼ 30:25 Å), at a distance of %54 A ˚from the wall.The distance from the wall aims to prevent any potential interactions between the droplet and the walls.This initial setup results in significant overlap between water and benzene molecules.To circumvent this problem, we delete every benzene molecule which has at least one of its atoms within a distance of 0.6 A ˚from an atom in a water molecule.
A number of potentials have been used to model intermolecular interactions.The wall atoms were attached onto their initial lattice sites using spring potentials.Furthermore, the intermolecular interactions between the wall atoms was modeled using basic LJ potentials.We modeled the water molecules using the TIP4P/2005 model. 70In general, there is no water model that outperforms its alternatives in predicting all fluid properties.Overall TIP4P/2005 results in more accurate calculations when it comes to most properties relevant to this study, such as viscosity, surface tension, and diffusivity.While another viable option would be the SPC/E model, we favored the TIP4P/2005 model, as we believe it has been studied more extensively, providing us with a solid body of numerical work against which we could validate our model.Finally, the benzene molecules were modeled using the OPLSAA model, 72 which was specifically designed for organic molecules in the liquid phase and has been used in similar MD simulations. 73he LJ parameters, e and r, which correspond to the depth of the potential well and the molecular diameter, are presented in Table I, while the charges for the water and benzene molecules are presented in Table II.For benzene and water, the parameters between unlike atoms, either between molecules of the same kind or between unlike molecules (also shown in Table I), are calculated empirically using the geometric mixing rule where the indices i and j correspond to different atoms.We retrospectively mention that the Coulomb forces between atoms across water (c) the initial and final configuration demonstrating that the droplet has moved from the cold to the hot wall.As a visual aid, the benzene atoms in (a) and (c) where made transparent in the region surrounding the water droplet.The results were visualized using OVITO. 80nd benzene molecules were removed.The reason, which we detail in Sec.IV A, is that the inclusion of the Coulomb potential between the two types of molecules resulted in an almost vanishing gradient of the surface tension with temperature, which is the driving force of thermocapillary flows.
The LJ potentials between the walls and liquid atoms are relatively arbitrary, as long as they meet certain requirements.We generally want weak interactions between solid and water atoms, in order to minimize the effect of electrostatic forces on the hydrodynamics of the droplet.On the contrary, the interactions between the wall and benzene atoms should be relatively strong, thus minimizing temperature jumps at the solid-liquid interface. 82Yet we avoid very strong interactions as this would introduce liquid stratification across a large percentage of the channel, which might influence the droplet dynamics.
Note that we use the same potential parameters across all our simulations, for all temperatures and pressures.These parameters determine the properties of our fluids (e.g., thermal conductivity and viscosity) under different thermodynamic states.While previous MD simulations have already calculated such properties for water and benzene, they have not considered the full range of temperatures and pressures relevant to this study.Therefore, in addition to validation purposes, we spend a significant part of this paper, i.e., Sec.IV A, in evaluating these properties under different conditions.
The temperature of the channel was controlled exclusively through Langevin thermostats, applied only at the solid walls.No artifacts where introduced in the fluid, as thermostating the liquid can potentially result in unphysical results. 83Heat transfer across the solid-liquid interface was achieved naturally, i.e., through conduction.
The temperature of the lower wall was kept constant at T cold ¼ 275 K across all cases.Since the temperature difference is a key parameter of the study, the temperature at the upper wall varied between different cases, with T hot ¼ 325, 375, 425, 475 K corresponding to temperature differences of DT ¼ 50; 100; 150; 200 K, respectively.Note that, considering the small length scales characterizing our system, these temperature differences correspond to extremely large temperature gradients.These, however, are required in order to obtain Marangoni numbers that will allow the simulations to complete within a reasonable time.
Note, that the temperature at the hot wall, goes beyond the boiling temperature of water at atmospheric conditions.However, since the fluid is confined by walls, and considering that we always start with the same number of atoms, and therefore density, in our system, our fluids always remain in liquid form, and the pressure increases between different cases, reaching values of up to 2000 bars.
Starting from the initial configuration of the system, in which all atoms were arranged on a lattice, a number of short simulations were performed in order to develop a steady temperature gradient, and for the water molecules, initially arranged in a cube, to form a droplet.During this phase, the momentum of the center of mass of the water is zeroed at every time step.This allows the droplet to form and interact with the surrounding benzene naturally but ensures that the droplet will start from the same position across all cases.
The first step to this sequence was a short, energy-minimization process that smoothly and iteratively adjusts atomic positions through a backtracking optimization algorithm, leaving the system in a configuration that will not result in violent dynamics.Subsequently, the system is integrated in time using the Verlet algorithm, 81 with a time step Ds ¼ 1fs.For a short number of timesteps (2 Â 10 4 Ds), we added a Langevin thermostat, set to the average temperature of the system (ðT cold þ T hot Þ=2), on all liquid atoms.While this phase is not necessary, it can potentially reduce the time required for the system to reach a steady state.We then remove the Langevin thermostat from the liquid, and run the simulation for a further 2 Â 10 6 Ds, after which the system is ready, i.e., the temperature gradient is steady and linear, and the water molecules have formed a relatively spherical droplet.
Afterwards, we removed the restriction on the droplet's momentum, allowing it to move freely across the channel.We stress that the only artifacts at this point are the Langevin thermostats applied to the walls.
Once we remove the restriction on the droplet's momentum, we start measuring its position and velocity.We consider the simulation to be finalized when the droplet is in the proximity of the heated wall; specifically when it passes the point z ¼ 180 Å.As this depends on the velocity of the droplet, the total simulation time differs significantly between cases.
All simulations (as well as many of the calculations) were performed using the molecular simulator LAMMPS Plimpton 84 (http:// lammps.sandia.gov).

B. Calculation of thermodynamic properties
For the calculation of the thermal conductivity and viscosity of water and benzene, we considered each liquid individually in a simulation box with periodic boundary conditions along all three dimensions.We initially equilibrated the system within the isobaric, isothermal ensemble (NPT).Once we converged to the correct size for the simulation box (i.e., corresponding to the desired pressure and  temperature), the thermal conductivity and viscosity were calculated using the red-Kubo relations within the canonical ensemble (NVT).Specifically, we calculated the thermal conductivity as where j is the thermal conductivity; V is the volume of the system; k B is the Boltzmann constant; D is the number of dimensions of the system; the angled brackets indicate the ensemble average of the autocorrelation function of the microscopic heat flux, J, which is given by where v i is the velocity of atom i; E i its total energy (kinetic and potential); r ij is the interatomic distance between atom i and j; and F ij the pairwise force between the two atoms.Similarly, we calculated the viscosity, l using where now we consider the autocorrelation function of the pressure, i.e., the diagonal components of the stress tensor.
For the calculation of the interfacial tension between water and benzene, we considered a system in which the two liquids are stratified.If the interface is in the xy plane (i.e., normal to the z axes), we calculate the interfacial tension as where H Z is the height of the simulation box in the dimension normal to the interface, P N is the component of the total pressure of the system normal to the interface, and is the tangential component.The factor of 1 2 in Eq. ( 9) accounts for the periodic boundaries (i.e., there are two interfaces).
The exact same procedure is used for the calculation of the surface tension, with the exception that we remove one of the liquids.Note that introducing a vacuum over the liquid should initiate evaporation.Indeed, we have observed sporadic molecules passing through the interface and into the empty space.However, due to the small length scales that we consider, only a tiny fraction of free molecules can be seen at any time.Therefore, these are not expected to affect the calculations in any way.

IV. RESULTS AND DISCUSSION A. Characterization and validation of the MD models
In this section, we present the calculations of the thermodynamic properties that are relevant to this study.Specifically, we investigate the thermal conductivity, viscosity, and surface/interfacial tension of our liquids, which are relevant for thermocapillary flows.As water has been investigated by a number of MD studies, comparison of our results with the literature will serve as validation.On the other hand, MD studies on benzene are sparse and, to the best of our knowledge, there are no studies that consider benzene using the OPLSAA model.
Our calculations, therefore, are not only important to the objectives of this paper but will also provide a good reference point for future MD investigations of benzene.
We first consider the effect of temperature on the thermal conductivity of water for pressures 1, 1000, and 2000 bars [see Fig. 3(a)].At 1 bar our MD results overestimate the experimental thermal conductivity 85 by approximately 50% [cf.blue solid line to blue circles on the top left graph of Fig. 3(a)].This error is a known shortcoming of the TIP4P water model (as well as most other water models).With this in mind, our results are in good agreement with the thermal conductivity calculated by previous nonequilibrium MD simulations that also used the TIPI4P water model 86 [blue diamonds on top left plot in Fig. 3(a)].
The thermal conductivity of water increases with increasing pressure [cf.blue-solid, purple-dashed, and red-dotted lines in Fig. 3(a)].Its relationship with temperature, however, is a little more complex.For low temperatures, the thermal conductivity increases quite linearly.For high temperatures, i.e., >425 K, the thermal conductivity plateaus and even decreases for pressures 1 and 1000 bars.Arguably, this decrease is due to the temperature being significantly greater than the boiling point of water, resulting in phase change.
The MD results for the thermal conductivity of benzene are in excellent agreement with experimental results, 87 at least at a pressure of 1 bar [cf.blue solid line with blue, left-pointing triangle, markers in Fig. 3(b)].As with water, we observe that the thermal conductivity of benzene increases with pressure.Unlike water, however, increasing the temperature results in a decreasing thermal conductivity for the pressure range that we consider here.
The viscosity of the TIP4P model is in good agreement with the experimental data, 88 as well as with previous MD studies that have used the TIP4P model 89 [cf.blue-solid line, up-pointing triangle, and right-pointing triangle in Fig. 3(c)].As expected, the viscosity decreases exponentially with temperature while pressure has a minor impact.
The viscosity of benzene also seems to exhibit a sharp decrease with temperature.Yet pressure, unrealistically, seems to play a much greater role, with greater pressure resulting in greater viscosity.The MD calculations at 1 bar have errors of approximately 25% error compared to experimental data. 90he calculated surface tension for water is in good agreement with previous MD studies that use TIP4P 91 [cf.blue solid line with blue diamonds in Fig. 4(a)].The linearly decreasing relationship between the surface tension and temperature is expected, as it characterizes most liquids. 92ote a 15% difference between the MD results and experimental data 91 [cf.blue solid line with blue circles in Fig. 4(a)].Vega and De Miguel 91 show that this difference is accounted for by the tail of the LJ potential, i.e., the portion of the potential beyond the cutoff distance that we consider.Instead of increasing the cutoff distance, which would significantly increase the computational time, such effects can be postprocessed into the results by adding empirical terms into the energy and pressure. 93These adjustments, however, do not reflect the dynamic behavior of the system.In this study, we use the calculated thermodynamic properties to predict the thermocapillary motion of the droplet and, as such, we do not consider these correcting terms.
Similarly to water, the surface tension of benzene decreases linearly with temperature [see purple-dotted line in Fig. 4(a)].Although the OPLSAA model seems to underestimate the experimental data 94 by approximately 30%, the gradient of the surface tension with temperature is very similar [cf.red-dotted line with up-pointing triangles on the left plot of Fig. 4(a)].
Finally, we look at the interfacial tension between water and benzene [see Fig. 4(b)].When our potentials include the Coulombic interactions between the water and benzene atoms, we get, on average, relatively good agreement with experimental data for a temperature range 300 K T 375 K (Ref.95) [cf.dashed-dotted line and blue circles in Fig. 4(b)].Yet the gradient of the interfacial tension with temperature is practically nonexistent, a behavior that was also observed in previous MD studies. 96As thermocapillary flows are driven by the gradients of the surface tension with the temperature, such behavior is unacceptable for this investigation.
If we ignore the Coulombic interactions between the water and benzene molecules, the error with the experimental data is almost 100% across the temperature range that we consider.Yet the slope of the curve, which we re-emphasize is more relevant to this study, is within 5% of the slope of the experimental results.
We also note that the interfacial tension seems to increase significantly with increasing pressure.This behavior is inconsistent with experimental results that show no significant correlation between pressure and interfacial tension, at least for pressure of up to 800 bars. 95

B. Thermocapillary migration of the nanodroplet
Before starting to deal with the outcomes of the MD simulations, we wish to recall the analytical formula originally derived by Young et al. 7 with which we will compare the droplet migration velocity determined by means of MD.This formula is not valid in circumstances for which the Reynolds and/or the Marangoni number take finite appreciable values.Strictly speaking, it is relevant only to situations in which particles have a size sufficiently large to be considered non-Brownian, but inertia and departure from a thermally diffusive field are negligible (a situation generally known as the "Stokesian limit").Notably, in this limit, the momentum and energy transport (macroscopic) equations become linear.
Under the assumption of ðRe; MaÞ !0, considering steady-state conditions and absence of gravity, i.e., g ¼ 0, and taking advantage of the linearity of the governing equations, the analytical solution for the nondimensional terminal velocity of the drop can be expressed as 78 As already explained to a certain extent in the introduction, Eq. ( 10) is also affected by other limitations or constraints.It relates to perfectly spherical droplets evolving in a perfectly immiscible external matrix [this being the macroscopic equilibrium shape in the limit as ðRe; MaÞ !0] and relies on the assumption of an unbounded domain.Moreover, it assumes that all material properties are unaffected by temperature variations (except for the interfacial tension).In the present context, however, material properties (density, viscosity, and thermal conductivity) stem from the properties of and interaction between molecules.Therefore, the MD model naturally accounts for temperature variations of fluid properties that, in turn, may bring significant differences between the present approach and the continuum model.
For the conditions considered in the present work [i.e., Re ¼ Oð10 À3 Þ and Ma ¼ Oð10 À2 Þ], the inertia and convective transport of heat can still be considered negligible and the Stokesian limit in which the equations of fluid dynamics become linear is also still applicable.Nevertheless, some of the assumptions on which the F-YGB formula relies are no longer satisfied (we will come back to the implications of all these observations later).
The outcomes of the present simulations can be gathered from Fig. 5.
Once the system is initialized and all the constraints on the droplet are removed, the droplet starts moving toward the hot end of the wall [see Fig. 5(a)].Furthermore, on increasing the temperature difference between the two walls, the droplet moves faster and reaches the other side of the channel in a shorter amount of time.This is the qualitative behavior that we would expect to see macroscopically [as formalized by Eq. ( 10)].
As already demonstrated in Fig. 4(b), when the forces between benzene and water included the Coulombic contribution, the gradient of the interfacial tension with temperature vanishes.The purple, dashed line in Fig. 5(a) indicates, that in the absence of such a gradient, the droplet only exhibits Brownian motion, despite the temperature difference of DT ¼ 100 K.This suggests that other thermophoretic forces are absent.
Among other things, the migration of the droplet in this study is interesting, as it illustrates how the statistical treatment of the many particles of a fluid, with a key set of assumptions, can be used to derive in a relatively straightforward way the governing laws of continuum (in other words, through MD simulations, it is possible to recover for the droplet dynamics the same functional dependencies envisaged by continuum-based models).These can be obtained as the natural consequence of the ensemble behavior of many particles interacting with each other, that is, the macrophysical properties of fluids and related hydrodynamic limits for stochastic particle systems can be obtained starting from microscopic statistics [a concept that also applies to the fluid stress tensor, as implicitly indicated by Eq. ( 2)].
Remarkably, however, these simulations have also led to another very important conclusion or realization.They have revealed that if the problem is considered on such small scales, different initial microscopic configurations (i.e., a different randomly allocated velocities distribution) can yield significantly different results even though the problem parameters and the considered boundary conditions remain unchanged.For example, for three different runs of the same case with DT ¼ 100 K, the droplet crosses the channel at different times: %0:7 Â 10 8 ; % 1 Â 10 8 , and %1:3 Â 10 8 fs [see Fig. 5(b)].This suggests that the stochastic, Brownian dynamics are influencing the timeaveraged behavior of our system causing a significant loss in the predictability (in a deterministic sense) of its behavior.Among other things, in order to place this interpretation in a proper mathematical context, it would suffice to consider that the random motion of molecules directly affects the dynamics via the impact it has on the viscous stress [as formalized by Eq. ( 2), the stress tensor can be obtained as a stochastic measure of the exchange of microscopic momentum].This stress in turn directly enters the balance of forces responsible for the migration of the droplet [Eq.( 1)].
It should also be noted that in Fig. 5(a) we selected a single run for each different DT.The curves were selected in order to demonstrate the general trend that increasing the temperature difference increases the velocity.As demonstrated in Fig. 5(b), these are not unique, and Brownian dynamics result in large variations.
Following up on the previous point, a simple estimate of the relative importance of Brownian motion and thermal convection can be provided by evaluating the P eclet number, Pe ¼ U T a=D 0 , where D 0 represents the Einstein-Stokes diffusion coefficient.Recalling the expression of the "thermal" characteristic velocity U T ¼ c T r 1 Ta=l and taking into account that in our conditions, we have c  between 300 and 400 K. Using all the information gathered so far the P eclet number falls in the interval Pe ¼ Oð1Þ to Pe ¼ Oð10 1 Þ; this clearly indicates that the convective transport of momentum is comparable to Brownian diffusion.By contrast, assuming a droplet radius on the order of the millimeter, i.e., a ¼ Oð10 À3 Þ m and a reasonable temperature gradient of the order of r 1 T ¼ Oð10 3 Þ Km À1 and basing the evaluation on the same material properties, i.e., the same values of c T and l considered before, the P eclet number would read Pe ¼ Oð10 11 Þ which, on the one hand, shows that at macroscopic length scales Brownian motion is insignificant and, on the other hand, can be used to justify the peculiar dynamics determined in the current flow conditions.
In addition to the P eclet number, the scale between the expected Marangoni velocity and that due to Brownian motion have also been compared.Considering the thermodynamic properties calculated in Sec.IV A at 300 K, on the basis of Eq. ( 10), one should expect a steady-state velocity $0:03 ms À1 .On the other hand, calculating a characteristic velocity for the Brownian motion by considering the equipartition theorem where D is the dimensionality of the problem; and m is the mass of the droplet, for a 1D problem, and considering the density of water and the radius of the droplet, one should get v RMS $ 15 ms À1 .Indeed, we have found the instantaneous velocity of the droplet to fluctuate rapidly about zero, with values jvj $ 10 ms À1 [see Fig. 5(c)].This is very similar to the Brownian velocity, i.e., approximately three orders of magnitude greater than the expected thermocapillary velocity.To mitigate this inherent noise, we have used time-averaging; the effects of Brownian motion should indeed vanish when a sufficiently large time window is considered.
Taking a moving average of the velocities for each case, using a window of 20 ns (i.e., 20 Â 10 6 timesteps), we have observed mostly positive velocities which are much closer to our macroscopic expectations [see Fig. 5(d)].Yet we still observe significant velocity variations within each case, as well as between different runs of the same case.
To reduce this uncertainty as much as possible, we calculated the velocity for each case based on the entire duration of the simulation (i.e., based on the total time it took the droplet to cross the channel).We further considered five runs of the same case and calculated the average and standard error for each DT.
The outcomes of these calculations are qualitatively and quantitatively substantiated in Fig. 6, where the MD results and the prediction of the analytical solution [Eq. ( 10)] are plotted at the same time.Most notably, both curves display a similar trend: the velocity increases in a relatively linear fashion with a similar gradient.We note that the MD results for DT ¼ 150 K seem to deviate from this linear trend, but we ascribe this to the large standard error associated with that point, as indicated by the error bar.
Interestingly, in terms of scale, the MD results overestimate the analytical solution by up to approximately 100% for the lower temperature differences and down to 50% for the higher temperatures.
We argue that an explanation for these discrepancies should be sought in a number of key factors, which can be summarized as follows.
A first explanation for these differences can be elaborated in its simplest form taking into account that for the analytical model, we use the thermodynamic properties corresponding to the temperature at the center of the channel.This would be accurate if the thermodynamic properties varied linearly with temperature, which (as witnessed by Fig. 3) is not the case.
A further understanding of the observed discrepancies can be gained by stressing that some of the assumptions on which the analytical solution [Eq.(10)] relies do not hold here.The continuum model assumes a perfectly spherical, immiscible droplet, whereas our nanodroplet constantly deforms and is partly miscible.To elucidate further the significance of this observation, one should keep in mind that with the MD approach no boundary (in a mathematical sense) exists that physically separates the two different fluids considered.Their separation is essentially a result of the attractive forces among molecules.Put differently, it is an "outcome" of the simulation, not an intrinsic property of the model (as it would be for methods such as the VOF or the level-set mentioned in the introduction or analytical approaches such as that implemented by Ref. 7).
The analytical solution also requires that the domain is unbounded, whereas the height of our channel has dimensions comparable to the droplet's diameter.
Another reason could be rooted in the specific approach implemented for the calculation of the thermodynamic properties.In Sec.IV A, we considered each fluid individually, with periodic boundary conditions emulating a perpetual continuation of the domain.8][99] Furthermore, the water nanodroplet, less than 2nm in radius, is so small that all its molecules are affected by the surrounding inhomogeneities.This can have an effect on the thermodynamic properties of the liquid, compared to its bulk form.As a final remark, we wish to mention that the standard error of the MD results seems to increase with increasing temperature (see error bars in Fig. 6).The reason for such a behavior can be traced back to the overall time required by the droplet to cross the channel.For large temperature gradients, the droplet moves at a higher velocity than for smaller temperature gradients, resulting in a shorter duration over which the velocity is averaged.As the time average of Brownian effects converges to zero with increasing time, the uncertainty in the velocity is less for smaller temperature gradients.For DT ¼ 150K and DT ¼ 200K, the error bars have increased to the point that they overlap.This indicates that the velocity of the droplet is not consistently greater for DT ¼ 200K compared to DT ¼ 150K, i.e., some of the runs for the case where DT ¼ 150K result in a higher droplet velocity compared to some of the runs for DT ¼ 200K.Yet, on average, the trend persists and greater temperature gradients lead to a greater droplet velocity.

V. CONCLUSIONS
In this study, the thermocapillary motion of a water nanodroplet suspended in benzene has been investigated in the framework of an MD approach.Results have been compared with the F-YGB model by Fedosov, 7,34 an analytical formula relying on the assumption of continuum, which can adequately describe the steady-state, thermocapillary velocity of a spherical droplet for vanishingly small Reynolds and Marangoni number, i.e., in the limit as inertial and convective thermal effects can be neglected.The present analysis should therefore be regarded as a first attempt aimed to discern the genesis of the so-called (thermally induced) vectorial migration of droplets at a very small scale and "at a fundamental level" (to the best of our knowledge no other efforts of such a kind can be found in the existing literature).
It has been motivated essentially by the recent explosion of interest in the field of "microfabrication," "digital fluidic operations," and related aspects.These concepts have opened the door to new groundbreaking ideas such as the use of "deterministic hydrodynamics" for the design of new technological applications [100][101][102] by which distinct (liquid-gas, solid-liquid, or liquid-liquid) phases can be controlled on the basis of innovative strategies taking advantage of the intrinsic fluid-dynamic properties of the considered systems.We have shown that the applicability of these concepts on a very small (nano) scale might be hindered by factors, which in general are not present when micrometric (or with larger size) droplets are considered.
As made evident by the present results, with the descending scale, the Brownian motion can significantly interfere with the droplet migration process causing significant variations in its velocity with respect to the predictions of available macroscopic frameworks.On the one hand, the microscopic agitation of molecules can cause quantitative variations in the velocity; as such, it reduces the predictability of the solutions.As expected, the latter essentially manifests as a loss of determinism (making the connection between the initial conditions and the ensuing dynamics somehow elusive).
Nevertheless, by means of this framework, we could successfully determine the average behavior of particles and reveal their mean evolution.Such statistics, directly connected to the (emerging) macroscopic behavior of the temporally evolving droplet, have been compared with the functional relationships envisaged by the abovementioned analytical formula.This modus operandi has led to extend the validity of those functional dependencies to the nanoscale.Moreover, as a natural consequence of the strategy used to address the problem, we have also built a set of new solutions valid on this scale.
The discrepancies between the MD results for the nanodroplet and the corresponding ones predicted by the macroscopic (analytical) model have been justified on the basis of the triadic relationship established among the ability of Brownian effects to make the shape of the droplet extremely irregular, the effective confinement experienced by the droplet (as opposed to the assumption of infinite domain) and the intrinsic simplifications at the basis of the analytical formula (strictly speaking, valid for Re ¼ 0 and Ma ¼ 0 and applied in the present work using the thermodynamic properties evaluated in the center of the channel).
Finally, a necessary requirement for the calculation of analytical solution was an investigation of the thermodynamic properties (thermal conductivity, viscosity, and surface/interfacial tension) of water and benzene at different pressures and temperatures.To the best of our knowledge, a lot of these data are currently unavailable in the literature (this excludes investigations on the TIP4P water model, which has been studied extensively and has a number of excellent references available).Therefore, as a side contribution, we believe that this paper will also serve as a good academic reference for the modeling binary mixtures of water and benzene (and perhaps other hydrocarbons) using MD.

FIG. 1 .
FIG. 1. Schematic representation of a droplet suspended in an unbounded, otherwise quiescent fluid migrating according to the direction of the imposed temperature gradient r 1 T.

FIG. 2 .
FIG. 2. Schematic representation of the system, showing (a) the full three-dimensional molecular model; (b) the direction of the temperature gradient imposed on the system;(c) the initial and final configuration demonstrating that the droplet has moved from the cold to the hot wall.As a visual aid, the benzene atoms in (a) and (c) where made transparent in the region surrounding the water droplet.The results were visualized using OVITO.80

FIG. 3 .
FIG. 3. Thermal conductivity (first row) and viscosity (second row) of water (first column) and benzene (second column).The lines correspond to our MD calculations of the thermal properties as a function of temperature, for different pressures.The error bars indicate the standard error for five different runs of each case.The markers correspond to computational (filled symbols) or experimental (empty symbols) data from the literature.(a) Water thermal conductivity; (b) benzene thermal conductivity; (c) water viscosity; and (d) benzene viscosity.

FIG. 4 .
FIG. 4. (a) The surface tension (i.e., interfacial tension with vacuum) of water and benzene at 1 bar.Our MD calculations are indicated by the lines with error bars.The markers correspond to computational (filled symbols) or experimental (empty symbols) data from the literature.(b) Interfacial tension between water and benzene as a function of temperature at different pressures.The markers correspond to experimental data.
and l ¼ Oð10 À3 Þ Pa s, we find U T ¼ Oð1Þ ms À1 .Moreover, the Einstein-Stokes diffusivity can be expressed as D 0 ¼ k B T=6pal, where k B is the Boltzmann constant and T is the absolute temperature, which in our setup ranges approximately

FIG. 5 .
FIG. 5. (a) Position of the center of mass of the droplet as a function of time, subjected to different temperature gradients.The purple, dashed line, labeled as dc=dT ¼ 0, corresponds to a case with DT ¼ 100 K in which the Coulombic forces between water and benzene were considered, leading to a vanishing gradient of the surface tension with respect to temperature; (b) position of the center of mass of the droplet as a function of time.The three curves correspond to different runs of the case where DT ¼ 100K, each starting from a different microstate.Differently colored regions correspond to averaging windows [see subplot (d)].(c) Instantaneous velocity of the droplet as a function of time; (d) a moving average of the velocity for the cases outlined in (b), i.e., the figure above.The time window over which the average was taken by differently shaded areas in both graphs.

FIG. 6 .
FIG. 6. Velocity of the droplet's center of mass as calculated by our MD simulations and by Young's formula [Eq.(10)].

TABLE I .
Lennard-Jones parameters used in our molecular model.

TABLE II .
Electrical charges of the atoms in our system.