Relaxation of surface tension in the freesurface boundary layer of simple LennardJones liquids.
Abstract
In this paper we use molecular dynamics (MD) to answer a classical question: how does the surface tension on a liquid/gas interface appear? After defining surface tension from the first principles and performing several consistency checks, we perform a dynamic experiment with a single simple liquid nanodroplet. At time zero, we remove all molecules of the interfacial layer, creating a fresh bare interface with the bulk arrangement of molecules. After that the system evolves towards equilibrium, and the expected surface tension is reestablished. We found that the system relaxation consists of three distinct stages. First, the mechanical balance is quickly reestablished. During this process the notion of surface tension is meaningless. In the second stage, the surface tension equilibrates, and the density profile broadens to a value which we call ”intrinsic” interfacial width. During the third stage, the density profile continues to broaden due to capillary wave excitations, which does not however affect the surface tension. We have observed this scenario for monatomic LennardJones (LJ) liquid as well as for binary LJ mixtures at different temperatures, monitoring a wide range of physical observables.
I Introduction
The phenomenon of capillarity and capillary flows is at the heart of numerous natural processes and technological applications. They are ranging from coating devices, polymer films and emerging technologies, such as micro and nanofluidics, to biological and medical applications of fluid dynamics. There is a large class of flows involving liquidgas interfaces which are controlled by surface tension. In many of them the free surface area undergoes significant changes over a period of time comparable to the characteristic diffusion time across the interfacial layer. The phenomena associated with this class of freesurface flows can be found in many applications and include such processes as capillary pinchoff, coalescence and generation of drops, formation of cusp regions and collapse of bubbles [1–7]. For example, the characteristic time scale of a pinchoff process is of the order of , where , and , are liquid density, viscosity and equilibrium surface tension respectively [1–3]. As an example, this characteristic time for water is at room temperature. On the other hand, typical relaxation time associated with a diffusion process in the watervapour interfacial layer is , taking characteristic width of the watervapour interface at and the coefficient of selfdiffusion of water molecules at room temperature and bulk conditions at [8, 9]. At this rate of change of the surface area, the interface will be at nonequilibrium conditions even for simple liquids. Apparently, this will result in transient nonequilibrium density profiles and as a consequence to dynamic surface tension.
The effect of dynamic surface tension at a liquidgas interface has been well studied for situations involving interfacially active molecules, or surfactants, when the interfacial tension is directly connected with the surfactant concentration at the interface [10, 11]. But as a matter of fact little is known from the first microscopic principles about dynamic surface tension effects at the interfaces of simple liquids in the absence of surface active molecules and the topic is the subject of strong debates [12–24].
The difference between two kinds of interfacial layers, with and without surface active molecules, is fundamental. When surfactant molecules are present, the dynamic surface tension is controlled by diffusion from the bulk area in case of solvable surfactants or by surface diffusion in case of unsolvable surfactants. In both cases relevant characteristic length scale of the diffusion process is of the order of the whole system size, , which is usually much larger than the interfacial thickness . As a consequence characteristic relaxation time of surface tension, , is found to be much larger than the diffusion time on the length scale of the interfacial layer , and could be well of the order of seconds [10, 11]. Note that interfacially active molecules reduce the value of equilibrium surface tension, therefore surface tension always relaxes from higher to lower values.
The situation is different in the case of liquidgas interfaces of simple liquids in the absence of surfactants. First of all, the relaxation time can be only associated with a diffusion process on the length scale of the interfacial layer, . This length scale suggests a very short relaxation time, of the order of . Secondly, the value of surface tension at a nonequilibrium interface of simple fluids could be in principle smaller than the equilibrium value. As a result, the relaxation process would have the character of an increasing function of time, in contrast to the case when surfactant molecules are present. Apart from simple estimates, however, at present we have no direct reliable information about the properties of fresh bare nonequilibrium liquidgas interfacial layers and the mechanism of relaxation process in that kind of systems. At the same time, experiments and macroscopic analysis of liquid flows with forming interfaces have already demonstrated that even relaxation process on this short time scale, , should have substantial impact on the character of the whole flow [7, 15–17].
For example, theoretical analysis of the freesurface flow breakup has shown that asymptotic behaviour and scaling of solutions to the NavierStokes equation close to the point of pinchoff depends on the assumptions made about the relaxation time of the surface phase [2, 3, 16, 17]. The first similarity solution to the pinchoff problem has been proposed by Eggers (1993, 1997) in the framework of a standard hydrodynamic model [2, 3]. In this approach, the surface tension relaxation time is simply and the surface tension is always equal to its equilibrium value independent of the surface area rate of change. A comparison of the freesurface profiles calculated on the basis of the similarity solution with experiments has shown very good agreement, see recent results [18]. However, the integrity of the standard approach to the pinchoff problem has been questioned by Shikhmurzaev (2005, 2007) [16, 17] on the basis of apparent inconsistencies between parametric dependencies predicted by the standard theory and those found in experiments. This concerns, for example, the minimal diameter of the microthread, at which the pinchoff process is triggered. In the standard approach (), the diameter is a function of liquid viscosity. At the same time, it has been found in experiments on simple fluids and binary mixtures by Kowalewski (1996) [1] that the microthread diameter is independent of viscosity of the liquid. Similar conclusions have been drawn later by Shikhmurzaev (2005, 2007) [16, 17] on the basis of an interface formation theory, where the surface tension relaxation time is intrinsically finite .
Another specific feature that follows from the assumption is a singularity of the axial velocity of the liquid jet, , at the point of pinchoff . At the same time molecular dynamics simulations of liquid nanojets rupture have shown rather smooth transition up to the moment of the breakup [25, 26]. One way to hinder the singularity of the axial velocity, as is shown in [16, 17], is to reduce the surface tension at the neck of the liquid thread, where the rate of change of the surface area is maximal.
In summary, for the rigorous macroscopic description of free surface flows of simple liquids with substantial rate of change of the surface area, it is imperative to determine the relaxation time of the surface phase and the fundamental mechanism of the relaxation process. A similar conclusion has been drawn on the basis of experimental studies of fluid necks rupture in the presence of surfactants [11]. The main effect of surfactants on the pinchoff dynamics is an increase of the local interfacial tension at the location of the minimum neck radius when the point of pinchoff is approached. This trend, as one can see, is opposite to what is expected during the pinchoff of simple fluids without surfactants [16, 17].
Despite the fact that the experimental evidence points out that the results of macroscopic analysis obtained in the assumption of finite relaxation time seems to be in a better agreement with the observations and molecular dynamics simulations, the issue of surface tension relaxation time is far from being resolved. There is still no consensus and the topic of surface tension relaxation time in simple fluids is the subject of disputes in the scientific community.
In the current study, we are going to directly establish properties of nonequilibrium liquidgas interfacial layers and focus on the relaxation phenomena at liquidgas interfaces of simple fluids without surface active molecules. To achieve this, we turn to MD simulations. We will set up clearcut model systems with fresh and sharp interfacial layers having bulk local structure. This way, we can observe for the first time the process of recreation of a liquidgas interface from an initially sharp surface and measure the corresponding timedependent surface tension. Our main objective is to address the following questions:

What is the value of surface tension at the bare fresh nonequlibrium liquidgas interfacial layer of simple fluids?

What is the characteristic relaxation time of an infinitely sharp nonequilibrium density profile and the surface tension associated with it?

What are the mechanisms associated with this relaxation process?

When can we assume that the surface tension of an evolving liquidgas interfacial layer of simple liquids (or miscible mixtures) is at equilibrium?
Ii Mathematical model and methodology
We focus on a liquidgas interface of simple monatomic and binary liquids consisting of particles interacting by means of 126 LennardJones (LJ) potential with the cutoff
where is the distance between particles and , and are the energy and the length scale parameters of LJ potential. In what follows all length scales are normalised by , energy and temperature by , pressure or stress by , surface tension by , viscosity by and time by , where is the particle mass. To understand the actual parameter range, we note that in the case of liquid argon the best choice of LJ parameters with is , and (where denotes the Boltzmann constant) [27]. This choice of parameters yields exact agreement between the calculated by MD critical point density and temperature and the values determined from the measurements, such as . This results in the following estimates , , and .
The system we study is a sufficiently large liquid drop consisting of up to 75000 particles giving the maximum droplet radius up to . The droplets have been prepared by means of an equilibration process at constant temperature during with the time integration step , which is used in MD simulations throughout this study. The sizes of the droplets in our study have been chosen to be sufficiently large in comparison to both interfacial thickness and the Tolman length so that one can work in the planar limit and neglect the curvature dependence of surface tension [28–33]. The curvature dependence is proportional to , that is, neglecting higher order terms, with Tolman length found in between depending on the droplet size [32, 33]. The temperature has been controlled by means of a DPD thermostat with friction to preserve liquid motion. Note that the friction due to collisions is , so that DPD friction introduces only negligible corrections. At low temperatures, the drop had been first equilibrated at a higher temperature value during and then brought to the lower temperature, by the thermostat, and equilibrated again during . The computational box had reflecting boundary conditions and the size, which was usually larger by than the characteristic drop size to allow for the vapour phase to settle in without total evaporation of the droplet.
In the framework of this numerical approach, static properties of liquidgas/liquid interfacial layers have been intensively studied by many groups over the last three decades, see for example [8, 34–42]. In the next two sections, we will discuss main procedures used in this study to analyse static and transient density profiles and to calculate the surface tension in liquid drops by MD simulations in equilibrium and nonequilibrium conditions. We will briefly summarize the results obtained in this area so far, while at the same time making a comparison with our own simulations to establish a robust connection between previous observations of interfaces in static conditions and our new results obtained in nonequilibrium conditions. In section 3, we will consider evolution of nonequilibrium interfacial density profiles and surface tension in monatomic LJ liquid drops. In section 4, we will extend our analysis of relaxation to binary LJ liquid drops. All nonequilibrium sharp density profiles in this study have been created by a cut off procedure, that is by removing all particles of the surface phase, see Fig. 1 for illustration. We will return to this cut off procedure for detailed consideration later and now consider density profile in a drop.
ii.1 The density profile in monatomic LJ liquid drops
Typical equilibrium density profile in a monatomic LJ drop consisting of 75000 identical particles, obtained as a result of averaging over at , is shown in Fig. 2. The temperature is slightly above the value at the triple point of LJ liquids [43]. The observed values of the bulk and the vapour densities are in a good agreement with the values observed in similar conditions [37, 41]. As is seen in Fig. 2, the calculated density profiles can be accurately approximated by means of an error function
(1) 
where parameter represents the interfacial thickness (half of the actual ”visible” size of the interface) and is the average position of the interface (equimolar surface in terms of excess density). The origin of this approximation derives from the capillary wave theory, [44–50], and is used in this study for parametrisation of the density profiles. Alternatively, the density profile can be also approximated by a hyperbolic tangent fit
(2) 
which is based on a meanfield approach, [30]. Although it has been noted in [8] that the fit given by eq. (2) may produce a less accurate approximation than that given by eq. (1), in our case we have not observed essential differences, Fig. 2. That is, in a liquid drop consisting of particles at , the fourparameter fit (1) produces , , , and the fourparameter fit (2) gives , , , . The only difference observed is the value of the gas density . But, as one can see from the magnified view of the tail region (insets in Fig. 2), it is the error function fit which produces much better approximation of this region. On the other hand, if parameter is fixed in the fitting procedure (2), the other three parameters have shown only very weak variations with within . For example, at , the fourparameter fit (2) gives , , . So the observed difference can only indicate that the accuracy of (2) in the tail region is lower than that of (1). We believe that the observed difference in approximation is not important, but the preference in the current study to use the error function fit (1) in the data profiling comes actually from the observed dependence of the interfacial width on the droplet size , Fig. 3, which is the characteristic feature according to the capillary wave theory, see further discussion. We will briefly return to this issue of the fitting procedure later in the manuscript, in the part dedicated to nonequilibrium situations for comparison. We will see, that the evolution curves of the density profile width are identical within the approximation error and only characteristic time scales are weakly dependant of the choice of the fitting function.
According to the capillary wave theory, the observed interfacial density profile is apparent and conditioned by the widening of an ”intrinsic” density profile by a spectrum of thermally excited capillary waves, that is by the fluctuations of density on the length scale larger than the width of the interface. Under the assumption that the density fluctuations on the large length scale are decoupled from the short length scale density fluctuations and neglecting effects of curvature (assuming ) the density profile can be represented as the convolution of and the probability to find the interface at a position , , centred at ,
The particular form of the density distribution (1) can be obtained exactly in the case of the Gaussian probability distribution function
and infinitely sharp intrinsic profile
where is the Heaviside step function.
In general, though, when is not just a step function, the apparent interfacial width defined as
consists of two components
provided that , that is is a symmetric function around , . Here and are used to designate the first derivatives and . The first contribution,
can be interpreted as the ”intrinsic width” of the interfacial layer, which is a manifestation of the thermal fluctuations on the length scale shorter or comparable with , and the second contribution,
is due to the capillary waves [36, 45, 49].
At equilibrium, the interfacial width has a weak logarithmic dependence on the size of the droplet with the coefficient depending on the equilibrium value of the surface tension [36]
(3) 
In eq. (3), the summation is over all spherical harmonics with the mode number running in between , with the frequency
(4) 
and the mean square amplitude
(5) 
This effect of the capillary wave motion has been clearly observed and studied in MD simulations and in the experiments on simple fluids [8, 36, 48, 49, 5155]. In particular, MD simulations of nanoscale capillary waves (with amplitude particle diameters) in liquid argon drops have decisively demonstrated that capillary wave motion closely follows predictions of continuum hydrodynamics (frequency and damping coefficient), the observed temperature dependence of the density profile width is in a very good agreement with experiments on liquid argon [51, 52]. The width of the density profile was found to be within for for a drop of . At , this width is slightly smaller than the one observed in our MD simulations () due to larger cutoff distance of LJ potential introduced in [51] to reproduce real argon. This effect of the LJ potential cutoff has been clearly observed in [41], due to, in particular, larger surface tension values leading to smaller contribution from the capillary wave widening. The effects of capillary waves at a watervapour interface have been studied in detail in [8], confidently establishing that at the interface width is and the intrinsic width , that is about and of the water molecule diameter () respectively. A similar study of capillary waves at interfaces between coexisting phases in polymer blends has been conducted in [49].
0.5  

0.6  
0.7  
0.8 
The effect of the capillary waves is also clearly observed in our MD simulations, as is illustrated in Fig. 3, where the interface width is shown as a function of drop size . A logarithmic twoparameter fit to that dependence (Fig. 3) given by eq. (3) results in, for example at , and . This value of the surface tension obtained from the fitting procedure, eq. (3), is very close to , directly obtained in [40] and in this study using MD simulations in similar conditions at and . The accuracy of calculation appears to be lower than that of directly obtained values of , though the obtained values are essentially the same. But, one needs to note that the latter does not probably include all approximation errors coming from calculation of the density profiles from MD simulations, while the former probably does account for that. For consistency, we will use directly obtained values of in what follows. Note, that in our simulations, we have been able to go below the temperature at the triple point of monatomic LJ liquids, , without crystallisation to expand the range of parameters, such as viscosity of the liquid. To ensure that crystallisation does not occur in the simulations, the state of the liquid has been monitored by means of the mean square displacement of particles both in the whole drop and in the equivalent conditions in the bulk phase with periodic boundary conditions
(6) 
where is the total number of particles.
The obtained values of ”intrinsic width” and surface tension at different temperatures are presented in Table 1. One can see that in general the equilibrium intrinsic width is of the order of less than one particle diameter and is monotonically increasing function of temperature. The value of the surface tension calculated from the logarithmic fit deviates in general by less than 10% from the values obtained independently and in this study by direct methods. In what follows, we will use this approach derived from the capillary wave theory in the analysis of relaxation of the interfacial phase, although with caution, since the separation of the density fluctuations into two parts is artificial and the density profile and the state of the surface phase are the result of the full spectrum of density fluctuations [47].
ii.2 Calculation of surface tension in the drops
Instead of using empirical fit eq. (3), the value of the surface tension generated in the interfacial layer of a drop is calculated in our MD experiments in the planar interface limit directly from the microscopic stress tensor by means of
(7) 
where and are the tangential and normal components of the stress tensor respectively, is the surface of tension . For example, if the interface is flat and is perpendicular to axis, and . The microscopic stress tensor components have been calculated through the procedure developed in [29, 37, 56], the details of the method can be found in the Appendix.
In practice, the actual integration in eq. (7) takes place in the interval , that is between the points where the integrand vanishes, see Fig. 4.
The definition of the surface tension through the integral (7) results from a mechanical argument, which is thermodynamically welldefined in the macroscopic limit () and valid in case of sufficiently slow motion, when, should we neglect the inertia of the surface layer, the general condition of mechanical equilibrium is fulfilled
In case of a spherically symmetric geometry, this condition is reduced to
(8) 
Equation (8) is the necessary condition for the correct evaluation of surface tension in drops, which should be always monitored under nonequilibrium conditions.
The distribution of stresses in the radial direction can be accurately approximated by a Gaussian function
(9) 
where is the point of maximum of the distribution , Fig. 4.
The fit, eq. (9), gives for a drop consisting of 75000 particles at equilibrium at (this case is identical to the one shown in Fig. 2) , and . As one can see from Fig. 4, the maximum of the distribution of is shifted into the drop interior by , which is a typical value in our simulations. That distance is close to the values of Tolman length reported in [32, 33] and could be interpreted as the distance between equimolar surface and the surface of tension, though, as is known, the later is not welldefined thermodynamically and path dependant [29, 30, 56]. The obtained width of the profile of , , if we were to assume that it is also a result of capillary wave widening, allows to calculate intrinsic contribution . At , this intrinsic width is obtained from eq. (3). So, in general, the intrinsic width of the stress difference distribution is on the same length scale as the density profile intrinsic width , Tables 1 and 2. Note however that unlike the distribution of density, distribution of stresses has no direct physical meaning, since it depends on the choice of the contour connecting interacting particles [29, 30, 56], see Appendix for details. So those values of can only be used to check the consistency of the results. On the other hand, the integral (7) of that distribution is invariant of the choice of contour in the leading order of the limit of a planar interface , as expected, since in macroscopic limit it represents the measurable quantity  surface tension. The equilibrium values of the distribution parameters, given by eq. (9) at different temperatures are presented in Table 2 for liquid drops consisting of particles.
0.5  

0.6  
0.7  
0.8 
In the next part, we consider evolution and recreation of interfaces from an initially infinitely sharp density profile and apply described density profile fitting and surface tension calculation techniques at nonequilibrium conditions.
Iii Interfacial dynamics of monatomic liquids
To study dynamic properties of liquidgas interfacial layers, we create a fresh sharp interfacial surface by removing all particles at (Fig. 4) in an initially equilibrated drop consisting of 75000 particles, so that after the cut off the drop size is reduced to approximately 40000 particles, Fig. 1a,b. We can then observe the recreation process of a new interfacial layer and measure its properties directly after the cut off as functions of time to reveal relevant relaxation times.
In this study, we are primarily interested in the relaxation process in isothermal conditions. But, even that the DPD thermostat friction introduces only negligible corrections, we have done control runs for the calculations presented in the paper at nonequilibrium conditions with the thermostat switched off to understand if the results would be qualitatively different. We have found no substantial evidence of any influence from the thermostat on the data so far, which can give rise to qualitatively different results. Although, a detailed study of the relaxation process in nonisothermal conditions wil be done later.
There are several characteristic time scales that one can anticipate here. First of all, it is the stress relaxation time , which is usually associated with NonNewtonian behaviour of liquids in classical hydrodynamics. This characteristic time can be calculated in the bulk from the offdiagonal () stressstress correlation function, see for example [57],
The integral of this correlation function also provides the value of the zero shear rate viscosity in the liquid at the bulk conditions. The values of the longest stress relaxation time in monatomic LJ liquids and corresponding viscosities at different temperatures are summarised in Table. 3. We note that the stress relaxation time is smaller than in all cases.
If we were to assume that the density fluctuations on the large and small length scales at the interface are decoupled, then one can expect another two separate time scales: one associated with the intrinsic width of the interfacial density profile and another one associated with the excitation of capillary waves. In the later case we may expect a full spectrum of capillary waves, which depends on the size of the droplet and is quantified by the mode number , eq. (4). The largest contribution to has to come from the mode with , which has the maximum amplitude at equilibrium, eq. (5).
The excitation time of the capillary waves is defined for each mode by the wave damping coefficient [58],
In the underdamped regime, ,
(10) 
and in the overdamped regime ,
(11) 
That is, from (10) and (11), for the mode with the dominant contribution into
(12) 
and
(13) 
which is in the limit of strong damping
The later implies that has different scaling with viscosity in the underdamped and overdamped regimes.
For monatomic liquid drops, the mode is either in the underdamped or overdamped regime. For example at for a liquid drop consisting of particles (), the characteristic frequency is , and thus , but at , . The surface waves relaxation time then will be found in the range depending on the temperature of the liquid, see Table 3. The estimated value of the damping coefficient at is close to actually observed in MD experiments on capillary waves in monatomic LJ nano droplets [51]. Note that in the case of binary liquids studied later in the second part of our work, the capillary wave motion will be in the overdamped regime.
Consider now results of MD simulations at nonequilibrium conditions. Evolution of surface tension and the density profiles after the cut off (at , which is always the case unless otherwise stated) is shown in Figs. 5–6 at different temperatures. The observed values of surface tension have been averaged over 20 time steps () and over approximately 150–200 statistically independent experiments to reduce the noise, which is inherently present in calculations of stresses in molecular dynamics simulations. As is seen in Figs. 5–6, the evolution of the surface tension and the width of the density profile has three distinguishable and separated stages, , and .
The first stage takes place immediately after the cut off during , which is on the brink between macroscopic and molecular time scales (). The appearance of this characteristic time scale seems to be solely due to the artificial nature of our experiment. Once the particles of the equilibrium interfacial layer have been removed, the particles in the newly formed interfacial layer immediately experience strong uncompensated force acting from the bulk particles. Strictly speaking, the notion of macroscopic stress tensor is not welldefined during this stage and a definition of the stress tensor could only be given in terms of the ensemble average. If we were to calculate surface tension using eq. (7) and check the condition of mechanical equilibrium (8), we would see that the condition is not fulfilled. Indeed, if we present eq. (8) in an integral normalised form
(14) 
we will see that condition (14) is fulfilled at equilibrium within the accuracy of , Fig. 7a. But after the cut off at , not surpsingly, condition (14) is violated and the system is out of mechanical balance at the interfacial layer, Fig 7b. This simply implies that the definition and the measurement of the surface tension during this initial phase have no direct macroscopic physical meaning and may be described as a microscopic relaxation phase. The mechanical balance should be of course restituted either on the time scale of the shear stress relaxation time or on the time scale of the pressure wave generated by the uncompensated force in the normal to the interface direction.
This relaxation process of the restitution of the mechanical equilibrium is indeed observed. For example, if we consider the evolution of condition (14) in a drop at , then one can see that at the mechanical balance has been almost restored, Fig. 7b, and at the balance is observed with almost the same accuracy as at equilibrium conditions. The separation between the first two stages is slightly smeared at higher temperatures, for example at , Fig. 5, as the relaxation process seems to take longer due to lower viscosity.
We can conclude that the observation of the surface tension that would have any macroscopic meaning should start after the first transitional phase of relaxation. The observed minimum value of surface tension (Fig. 5) at , when the interface is still relatively sharp ( for all temperatures, Fig. 6) is about of the equilibrium value.
In the second stage, the mechanical equilibrium is observed. During this stage, the surface tension relaxes to its equilibrium value, while, at the same time, the density profile is still in transition to equilibrium, Fig. 6. This indicates that we observe two separate time scales; the first time scale is associated with the relaxation of the surface tension, while the second time scale is associated with the relaxation of the density profile to its equilibrium shape. The surface tension relaxation time can be determined by fitting the time dependence with
(15) 
where is introduced to account for the transition from the first, fast relaxation stage. This particular choice of the fitting function is motivated by simplicity and serves to roughly reveal temperature dependence of the characteristic time. It may not reflect the nature of the surface tension as a derivative of the density distribution of particles, for example. The fit is applied after the first stage is completed () as illustrated in Fig. 8. Note that the accuracy of that definition due to nonmonotonic and noisy character of function is low and leaves some room for uncertainty. While the errors bars for in Table. 2 may look good, we think they are too optimistic, since the overall deviation of the fit from the fitting data points in this case was in between and . As a function of temperature surface tension relaxation time is found to be in the range , Table. 2. As is seen, this characteristic time is always larger (much larger at high temperatures) than the stress relaxation time , and has a clear temperature dependence, which is in anticorrelation with that of , so that one can completely rule out their possible connection.
0.5  
0.6  
0.7  
0.8 
To quantify the time scales relevant to the relaxation process at the interface, consider evolution of , Fig. 6, in more detail. These evolution curves have monotonic time dependence and have substantially lower level of noise and thus represent less formidable task for interpolation and time scale analysis than the evolution curves of the surface tension. To extract characteristic time scales, we fit function by
(16) 
We have found that it was sufficient to use just three characteristic time scales, , in eq. (16) for accurate approximation of , see Fig. 9 and Table. 3, where the results of the fitting procedure are summarised for different temperatures. For comparison, the same evolution curve has been reproduced using hyperbolic tangent fit (2), Fig. 9. One can see that the results are essentially identical within the approximation errors. The characteristic time scales extracted using (16) are sufficiently close, considering sensitivity of the exponential fits and the errors of approximation.
The first relaxation time, , as expected from the evolution of the surface tension, is in the range and corresponds to the first, microscopic stage of surface tension relaxation. This characteristic time is very short and just on the verge of the resolution . The second relaxation time is found to be in the range , which is close to the range of , , in the same temperature interval. The temperature dependence of , Table 3, is in a good qualitative agreement with the temperature dependence of . While the characteristic time scales are not very much apart, considering the approximation errors involved, exact quantitative comparison is difficult due to the difficulty in defining from nonmonotonic and noisy curves. Anyway, the characteristic values and the trend with temperature of both and suggest that these characteristic time scales correspond to the second, macroscopic stage of the surface tension and intrinsic density profile relaxation in the system.
The third fitting parameter , on the other hand, is apparently related with the longest relaxation time associated with the excitation of capillary waves . This may be evidenced from the scaling of both and with temperature, see Table 3. One can observe a sufficiently good quantitative correlation between those two characteristic times. The appearance of capillary waves can be also illustrated if we consider evolution of the offdiagonal components of the gyration tensor
with time, for example , which is an indicator of deviations from sphericity. From this dependence, Fig. 10, one can see that the shape of the drop is very close to a sphere just after the cut off, (some small deviation is actually expected due to fluctuations in the bulk area) and then the asymmetry develops.
So far we have established that the first, shortest relaxation time of the surface phase after the cut off is related to microscopic processes, while the third, the longest relaxation time should be attributed to the excitation of capillary waves. To understand the nature of the second time scale, , we consider transport properties of LJ liquids.
iii.1 Particle transport in the liquid and at the interface.
To understand the characteristic time scales of particle transport, consider the mean square displacement , eq. (6), in the bulk. The mean square displacement is used to calculate coefficient of selfdiffusion in the bulk
as illustrated in Fig. 11a. This, in turn, is used to estimate the characteristic diffusion time on the length scale of intrinsic equilibrium width
Alternatively, the effective time , which a particle needs to travel the distance in the direction across the interface (in the bulk conditions) can be obtained directly via
(17) 
assuming isotropic motion. The values of the selfdiffusion coefficient in the bulk, characteristic diffusion times and at different temperatures are summarized in Table 3. One can see that the estimate is close to , while . This is anticipated since if we were to calculate the characteristic time directly from (17), we would find that particle motion on this time scale is close to but not yet in the fully developed diffusive regime, which is reached at approximately (though this value depends, of course, on temperature), Fig. 11a. At the same time, for monatomic liquids in the range of temperatures , see Table 3 and Fig. 11b for illustration. The temperature dependences of and reflect the fact that both coefficient of selfdiffusion and intrinsic width of the density profile are monotonically increasing functions of temperature. Two effects almost compensate each other and both and demonstrate rather weak dependence on temperature, weaker than it is suggested by the dependence of , for example.
One needs to note that while those estimations may be practical and informative, they do not tell, strictly speaking, the whole story. Firstly, the particle motion is not in a fully developed diffusive regime. Secondly, the particles in the first layer at the interface are more mobile than those in the bulk, though this can be moderated by the uncompensated force, which keeps the interface width finite. Thus, the characteristic times and may only provide an upper limit on the contribution of the diffusion into the relaxation at the interface. On the other hand, a direct calculation of the residence time in the interfacial layer may be even less accurate. The interface is widened by capillary wave motion. Given that on average , it would be difficult to exactly locate particles belonging to the ”intrinsic” interface. To estimate characteristic residence time of particles at the interface, we calculate the mean square displacement in the radial direction of particles initially contained in a layer of thickness at the interface located at , that is ,
(18) 
where is the number of particles in the layer. The result is shown in Fig. 12 in comparison with the mean square displacement obtained in the radial direction in the bulk. The corresponding time calculated via
is listed in Table 3. One can see that in general is close to and as expected . So on average particles at the interface are more mobile, though it is difficult to instrument an exact measure of their mobility due to the interfacial ”roughness”, that is uncertainty in their positions.
Compare now the obtained values of the second relaxation time and its temperature dependence with those of and . One can observe that all of them, , and are in a very good quantitative agreement in the temperature range used in our MD experiments. This correlation between , , and is a clear indication of their connection. If we now compare the width of the interface at , , then we will see that it is only slightly above , that is . This is also illustrated in Fig. 9. All this implies that during the second relaxation stage, the interfacial width has reached equilibrium intrinsic width , while the surface tension has relaxed to its equilibrium value . To elaborate on this statement, consider snapshots of the profiles of distribution of stresses, Fig. 13, and density, Fig. 14. The snapshots of the profiles are shown before the cut off and at the subsequent times after the cut off. Corresponding values of the fitting parameters (eq. (9) ) are given in Table 2. One can see that at (the value at the beginning of the interval ) the amplitude of the profile (which is the value of the surface tension) has just returned to its equilibrium value or is close to it at high temperatures. At the same time the width is only slightly larger than , which is expected if we consider simultaneous relaxation to and excitation of the capillary waves. Thus, not only the density profile has relaxed to its intrinsic value, but also the distribution of . Therefore Fig. 13d and Fig. 14d show the snapshots of the virgin interface without apparent widening by the capillary waves. We note that we did not clearly observe the oscillating intrinsic interface profiles as in [59], though at after the cutoff, Fig. 14b, the density profile shows some signs of the fine structure similar to the one observed in simulations [59]. The reason for that ”discrepancy” might be that the effect was weak and simply below our resolution in dynamic simulations. This will need further studies.
One can then directly relate the relaxation of the surface tension and the relaxation of the density profile to its intrinsic value, and simply define the surface tension relaxation time in monatomic liquids as
So from this one can conclude that, first of all, the liquidgas interface in purely monatomic LJ liquids without apparent widening by the capillary waves is basically very thin and structureless. Secondly, the relaxation time of this density profile, , is directly related to a diffusionlike process on the length scale of intrinsic width of the density profile . Thirdly, the surface tension relaxation time is directly connected with the relaxation time of the density profile and therefore is also defined by the diffusionlike process on the length scale of the intrinsic width of the density profile . And finally, the process of surface tension relaxation is separated from the process of excitation of capillary waves, described by the third relaxation time of the density profile, which is much larger than .
We will further clarify these statements in the next section dedicated to similar dynamical studies but in binary LJ liquids, where we would be able to go into slightly different parameter range, especially with different viscosities, so that we can further separate the time scales discovered in monatomic liquids.
Iv Interfacial dynamics of a binary liquid
In this section, we turn our attention to a drop of a binary LJ liquid to test generality of our findings. The binary liquid consists of two particles, types A and B, such that their number density ratio in the mixture is . The particles interact with each other by means of a LJ potential with parameters , . This choice of LJ parameters corresponds to the wellknown KobAnderson model, which is less susceptible to crystallisation and has different range of viscosities in the larger operational window of temperatures in comparison to the monatomic LJ liquids [60]. Macroscopic parameters of the binary liquid used in this study, namely density, viscosity, coefficients of selfdiffusion for both components and , surface tension and the stress relaxation time calculated in a drop consisting of particles (or by the bulk MD simulations with periodic boundary conditions at the equivalent density) are listed in Table 4 at different temperatures. The minimum temperature of the binary liquid used in our simulations was above the critical temperature , when the coefficient of selfdiffusion for both particles goes to zero with [60].
0.45  

0.50 