Interpretation of Young’s equation for a liquid droplet on a flat and smooth solid surface: Mechanical and thermodynamic routes with a simple Lennard-Jones liquid

In this study, we carried out molecular dynamics simulations of a cylindrical Lennard-Jones droplet on a flat and smooth solid surface and showed that Young’s equation as the relation among solid-liquid, solid-vapor, and liquid-vapor interfacial tensions γSL, γSV, and γLV, respectively, was applicable only under a very restricted condition. Using the fluid stress-tensor distribution, we examined the force balance in the surface-lateral direction exerted on a rectangular control volume set around the contact line. As the mechanical route, the fluid stress integrals along the two control surfaces normal to the solid-fluid interface were theoretically connected with γSL and γSV relative to the solid-vacuum interfacial tension γS0 by Bakker’s equation extended to solid-related interfaces via a thought experiment, for which the position of the solid-fluid interface plane was defined at the limit that the fluid molecules could reach. On the other hand, the fluid stress integral along the control surface lateral to the solid-fluid interface was connected with γLV by the Young-Laplace equation. Through this connection, we showed that Young’s equation was valid for a system in which the net lateral force exerted on the fluid molecules from the solid surface was zero around the contact line. Furthermore, we compared γSL − γS0 and γSV − γS0 obtained by the mechanical route with the solidliquid and solid-vapor works of adhesion obtained by the dry-surface method as one of the thermodynamic routes and showed that both routes resulted in a good agreement. In addition, the contact angle predicted by Young’s equation with these interfacial tensions corresponded well to the apparent droplet contact angle determined by using the previously defined position of the solid-fluid interface plane; however, our theoretical derivation indicated that this correspondence was achieved because the zero-lateral force condition was satisfied in the present system with a flat and smooth solid surface. These results indicated that the contact angle should be predicted not only by the interfacial tensions but also by the pinning force exerted around the contact line. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5053881 J. Chem. Phys. 150, 044701 (2019); doi: 10.1063/1.5053881 150, 044701-1


I. INTRODUCTION
The behavior of the contact line, where a liquid-vapor interface meets a solid surface, has long been a topic of interest in various science and engineering fields because it plays a key role in the wetting properties. [1][2][3] By introducing the concept of interfacial tensions and contact angle θ, Young's equation 4 is given by where γ SL , γ SV , and γ LV denote solid-liquid, solid-vapor, and liquid-vapor interfacial tensions, respectively. The contact angle is a common measure of wettability at the macroscopic scale. Historically, Young's equation (1) was first proposed based on the horizontal force balance of interfacial tensions exerted on the contact line in 1805 before the establishment of thermodynamics. 5 Instead of using the concept of force balance, Young's equation is often re-defined from a thermodynamic point of view. 1 Various models have been further put forward to capture the details of the contact line, such as introducing microscopic contact angle, 6 adding line tension term to Eq. (1), 7,8 and dealing with precursor films. 1,9 However, it is difficult to experimentally validate these models mainly because measuring the interfacial tensions γ SL and γ SV , which include the solid phase, is not trivial. 10,11 It should also be noted that recent micro-structuring techniques can produce a heterogeneous solid surface with a well-defined boundary between areas having different wettability, and it has been shown that this boundary played a key role in the wetting behavior. 5 On the other hand, from a microscopic point of view, pioneering studies by Kirkwood and Buff 12 provided the theoretical framework of surface tension based on the statistical mechanics, and molecular dynamics (MD) analysis has been applied for the microscopic understanding of wetting.  Early study by Nijmeijer et al. 13,14 using mono-atomic Lennard-Jones (LJ) fluid film on a solid surface indicated that the balance in Eq. (1) was applicable. In their work, the liquid-vapor and solid-liquid surfaces were simulated and corresponding interfacial tensions were obtained from the integration of the difference between the pressure tensor components in the surface-normal and surface-lateral directions for the corresponding interfaces. This type of extraction of the interfacial tensions based on Bakker's equation 3,34 is called the mechanical route 35 and has been adopted for several following studies, and most of them indicated that Young's macroscopic model was also applicable to simple mono-atomic liquids at the microscale when the apparent contact angle of an equilibrium droplet on a planer solid surface was used. [15][16][17][18][19][20] In addition to the mechanical route, a number of new techniques are proposed for the calculation of the solid-related interfacial tensions, including the test-area perturbation method, 21,22 free-energy-based methods using transition matrix Monte Carlo simulations, 23,24 and the thermodynamic integration (TI) by the phantom-wall [25][26][27] and dry-surface (DS). [28][29][30] Calculation of interfacial tensions of realistic molecular liquids, e.g., water or hexane, on a complex solid surface was made more approachable by using these methods because the extraction of local distributions of the stress tensor was no longer needed. Interestingly, Young's macroscopic model using the above calculated interfacial tensions sometimes did not give a reasonable estimate of the apparent droplet contact angle for these complex systems, although the reason of this discrepancy was not clear. [31][32][33] Related to this feature, a very long relaxation time to obtain an equilibrium droplet contact angle was also reported. 33 Two of the present authors have also carried out MD simulations of an argon droplet 20 and a water-alcohol mixture droplet 27,36 on a flat crystal, where we obtained the spatial distribution of the stress tensor in the fluids and extracted the solid-liquid, solid-vapor, and liquid-vapor interfacial tensions γ SL , γ SV , and γ LV , respectively. Although these studies followed the mechanical route, we also obtained the spatial distributions of the stress field instead of merely using the resulting integrals. General questions arose from these studies: whether the mechanical route and other energy-based routes give the same solid-related interfacial tensions γ SL and γ LV , and what kind of condition must be provided for Young's equation to hold.
To make clear these questions, in this study, we carried out at first MD simulations of a hemi-cylindrical Lennard-Jones droplet on a flat and smooth solid surface with focusing on the following three points regarding the mechanical route: • The original Bakker's equation is intended to be applied only for the liquid-vapor or liquid-gas interface, and whether this equation is really applicable to the solidrelated interfacial tensions γ SL and γ SV . • Regarding the integration interval, where the lower and upper limits of the integration should be. • We have a variation of choices of the stress tensor: including the liquid-solid component or excluding it, for surface-normal and surface-tangential components in the calculation of the stress tensor close to the solid surface, and which choice gives the proper description of the interfacial tensions. For these purposes, we examined in detail the surface-lateral force balance exerted on a control volume set around the contact line using the distributions of the fluid stress-tensor components. Then, we related the three interfacial tensions γ LV , γ SL , and γ SV with the fluid stress integral exerted on each face of this control volume with an aid of extended thought experiment for the solid-related interfacial tensions. At this stage, we made the above-mentioned three points clear. Then, we showed that Young's equation was applicable to the present system because the lateral force exerted from the solid surface on the fluid molecules was negligibly small. A proper definition of the droplet contact angle was also provided. Finally, we compared the solid-liquid and solid-vapor interfacial tensions obtained by the mechanical route with the solid-liquid and solid-vapor works of adhesion obtained by the dry-surface method as one of the thermodynamic routes.

A. Potential model
Generic particles interacting through a 12-6 Lennard-Jones potential were adopted as the fluid molecules for ease of physical understanding as in our previous study. 20 The 12-6 LJ potential given by was used for the interaction between fluid molecules, where r ij was the distance between the molecules i at position r i and j at r j , while and σ denoted the LJ energy and length parameters, respectively. This LJ interaction was truncated at a cutoff distance of r c = 3.5σ, and quadratic functions were added so that the potential and interaction force smoothly vanished at r c with the Heaviside step function Θ. The constant values of c 2 and c 0 are shown in our previous study. 20 Hereafter, fluid and wall molecules are denoted by "f" and "w," respectively, and the corresponding combinations are indicated by subscripts.
A face-centered cubic (FCC) crystal with an exposed (001) face was used as the solid wall in contact with the fluid, and the interaction potential between wall atoms was expressed by the harmonic potential for the nearest neighbors with an equilibrium distance r 0 . The solid-fluid interaction was also expressed by the LJ potential equivalent to Eq. (2), where the length parameter σ fw was given by the Lorentz mixing rule, and the energy parameter fw was changed in a parametric manner by multiplying a fluid-wall interaction coefficient η to the base value as η 0 fw , where the base value 0 fw was given by the Berthelot mixing rule as 0 fw = √ ff ww . The energy parameter fw was further multiplied by a coupling parameter in the SL-DS and SV-DS systems as described later.
In addition to these intermolecular potentials, we employed a one-dimensional field Φ p for the piston in the SL-DS systems given by as a function of the distance z i ≡ z p − z i between a z-normal plane at z = z p and fluid molecule i at a height z = z i . This potential field mimicked a mean potential field created by a single layer of solid atoms with a uniform area number density ρ n = r −2 0 , and by setting the cut-off distance at σ fw , only repulsive force was exerted on the fluids by this field. We also used a similar one-dimensional field Φ b given by as a function of the distance z i ≡ z b − z i between a z-normal plane fixed at z = z b on the top of the simulation cell and fluid molecule i at a height z = z i . This potential field was used as the particle bath in the SV-DS systems. Further details of the piston and particle bath are described in Sec. II B.

B. Simulation systems
For the droplet system, periodic boundary conditions were set in the x-and y-directions, and a hemi-cylindrically shaped droplet was formed on the solid surface with the droplet axis parallel to the y-axis as shown in Fig. 1(a) so that the effect of line tension 7,8,17,26,37 was neglected. A mirror boundary condition was employed at the top boundary in the z-direction, whereas a solid wall consisting of the FCC crystal was located on the bottom of the calculation cell which directed its (001) plane normal to the z-direction. This solid wall had eight layers so that possible minimum distance between the argon molecule and wall atom in the bottom layer was longer than the cut-off distance of the fluid-wall interaction potential.
The position of the wall atoms in the bottom layer of the base crystal was fixed, and the temperature of those in the second layer from the bottom, which were sufficiently far from the solid-liquid interface, was controlled by using the standard Langevin thermostat in all directions at a control temperature T c of 85 K with a Debye temperature of 240 K. 20,38 On the other hand, no direct thermostatting was imposed on the fluid molecules, and the fluid temperature was maintained only through the heat conduction with the solid so that the effects of thermostatting attached to the second-bottom layer of the solid would be negligible on the present equilibrium wetting behavior.
A cylindrical liquid droplet was first equilibrated away from the solid surface, and after its automatic adsorption onto a solid surface followed by a further equilibration run of 20 ns, an initial equilibrium hemi-cylindrical droplet was obtained. The average of 40 ns thereafter was used for the analysis, where we obtained two-dimensional distributions of density and stress components in the frame of reference relative to the center of mass of the droplet, considering that the droplet showed a random Brownian motion without noticeable pinning on the homogeneous and smooth solid surface in this study.
For the SL-DS systems shown in Fig. 1(b), a solid-liquid interface was formed between the liquid and bottom solid wall with the same crystal structure, exposed face, and corresponding wettability parameter η as the droplet system. A periodic boundary condition was employed in lateral x-and ydirections, and the placement and temperature control of the solid surface were the same as in the droplet system, whereas the system width L x was set smaller since we only dealt with solid-liquid two-phase interface. In addition, we set a piston above the liquid to attain a constant pressure. This piston interacted with the fluid molecules with the potential Φ p given in Eq. (3), and its height z p obeyed the following equation of motion: where m p , A, p, and p set denoted the mass of the piston, surface area (=L x L y ), system, and control pressures, respectively. The system pressure p was obtained from the total normal force exerted on the piston from argon molecules given by where N f was the number of fluid molecules. By allocating a sufficient number of fluid molecules and by setting the pressure p set above the vapor pressure, a liquid bulk with a constant density was formed between the solid wall and piston. We also applied the Nosé-Hoover thermostat 39 only to the fluid molecules to maintain the system temperature at T c , which was the same as the control temperature for the solid, to maintain a constant fluid temperature even in case the fluid and solid were almost non-interacting during the DS procedure described in Sec. II D. The relaxation time of the thermostat was 1 ps. Note that the effects of this thermostat on the DS procedure were negligibly small because the present SL-DS system was in equilibrium. The interaction energy parameter fw was multiplied by the coupling parameter 1 − λ in this system for the dry-surface scheme, and we obtained multiple equilibrium systems with various λ values with 0 ≤ λ < 1 to numerically calculate the thermodynamic integration, as described in Sec. II D. Each system was obtained after a preliminary equilibration over 5 ns, and the time average of 20 ns was used for the analysis.
For the SV-DS systems, we investigated the interfacial energy between saturated vapor and corresponding solid surface set on the bottom of the simulation cell by placing an additional particle bath on the top, as shown in Fig. 1(c). The setup regarding the periodic boundary conditions employed in lateral x-and y-directions, temperature control and placement conditions for the solid surface, and additional Nosé-Hoover thermostat for the fluid molecules were the same as the SL-DS system, whereas the particle bath was kept in place by the potential field Φ b in Eq. (4) at a fixed height sufficiently far from the solid surface. Note that the potential field mimicked a completely wettable surface with an equilibrium contact angle of zero with the present potential parameters, i.e., a liquid film was formed on the particle bath. With this setting, a solid-vapor interface with the same density distribution as that in a droplet system was achieved. We formed multiple equilibrium systems with various values of the coupling parameter λ with the same recipe as the SL-DS systems.
The velocity Verlet method was applied for the integration of the Newtonian equation of motion with a time increment δt of 5 fs for all systems. Values of the simulation parameters are summarized in Table I with the corresponding nondimensional ones, which are normalized by the corresponding standard values. Note that the η value ranged up to 0.5 for the droplet systems because the system with η = 0.6 showed complete wetting and no hemi-cylindrical droplet was formed.

C. Calculation of fluid stress distribution
For the droplet system, we extracted the distribution of the two-dimensional fluid stress tensor averaged in the axial direction. The method of plane (MoP) was adopted in the present study instead of using the volume average (VA) used in the previous study 20,27 because the exact balance satisfied for an arbitrary control volume bounded by a closed surface was one of the most important prerequisites of the formulation in this study. Note that we extracted the fluid stress as an internal force while considering all solid contributions to be an external force field 14,40,41 in this study. The local fluid stress tensor τ(x, z) was calculated by dividing the system into x-and z-normal flat bin faces with a width of 0.08 nm for both.
3.92 nm 11.5 L y 3.92 nm 11.5 L z (droplet) 16 The fluid stress tensor component τ αβ , which expresses the stress in the β-direction exerted on a surface element with an outward normal in the α-direction, is given by kinetic term τ kin αβ and inter-molecular interaction term τ int αβ as In the MoP, the kinetic term on an α-normal bin face with an area A α was calculated by where m i and v i denoted the mass and velocity vector of the ith fluid molecule, and e α and e β were the unit normal vectors in αand β-directions, respectively. The angular brackets denoted the time average, and the summation acrossAα i,δt was taken for every fluid molecule i passing through the bin face within a time interval of δt, which was equal to the time increment for the numerical integration. A switching function, 2Θ(v i · e α ) − 1, where Θ was the Heaviside step function, giving ±1 depending on the sign of v i · e α was included through the Heaviside step function Θ in the RHS of Eq. (8). Strictly speaking, the velocity v i in Eq. (8) should be a relative one to the average; however, we did not use the relative velocity considering that the average velocity of the droplet under random Brownian motion was sufficiently small.
On the other hand, the intermolecular interaction term in Eq. (7) was given by where r ij and f ij denoted the relative position vector r j − r i and force vector exerted on molecule j at position r j from molecule i at r i , respectively, and the summation acrossAα (i,j)∈fluid was taken for all line segments between r i and r j which crossed the bin face. Note that technically the fluid-solid interaction can also be included as i − j pair in the summation acrossAα (i,j) in Eq. (9), but only the fluid-fluid interaction was taken into account as the fluid stress as the internal force, and fluidsolid contribution was considered as an external force field in this study, as mentioned above. 14,40,41 This will be discussed in detail in Sec. III A. Note that Eqs. (8) and (9) have a negative sign because they express stress, whose diagonal component is pressure with inverted sign. For example, Eq. (8) always has a non-positive diagonal stress value, i.e., positive pressure value, because molecules passing always transfer moment, which results in pressure force directed towards the control surface from either side. On the other hand, Eq. (9) can be both negative and positive, where a negative diagonal stress, i.e., positive pressure, indicates that two molecules on opposing sides are repulsing each other, while the opposing sign indicates that the molecules are attracting each other.

D. Works of adhesion of solid-liquid and solid-vapor interfaces
Leroy and Müller-Plathe 28 proposed the Dry-Surface method in which the coupling parameter for the thermodynamic integration (TI) 42 was embedded in the fluid-wall interaction parameter. They formed a solid-liquid interface at first, and then by turning off the attractive part of the fluidwall interaction potential through the coupling parameter, the work of adhesion was extracted through the thermodynamic integration along a reversible path. In the present study, we adopted the DS method for the SL interface and for the SV interface as shown in Figs. 1(b) and 1(c), respectively. A coupling parameter λ was embedded in the fluid-wall interaction potential for both systems as where Φ LJ fw (r ij ) denoted the LJ interaction in Eq. (2). For the case of SL-DS in Fig. 1(b), we obtained equilibrium solid-liquid interfaces with discrete coupling parameter λ from 0 to 0.999. With the setup described in Subsection II B, the number of molecules N, pressure p, and temperature T were kept constant, i.e., we formed constant NpT systems. Note that the maximum value of λ was set slightly below 1 to keep the solid-fluid interaction to be effectively only repulsive. This value is denoted by 1 − hereafter. The SL interface at λ = 0 with Φ DS fw (r ij , λ) = Φ LJ fw (r ij ) was separated into solidvacuum and liquid-vacuum interfaces by changing the coupling parameter to λ = 1 − as shown in left and right panels The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp of Fig. 1(b), as no fluid molecules were adsorbed onto the solid surface because the fluid-wall interaction was almost without attraction at λ = 1 − . Hence, the difference of the Gibbs free energy ∆G ≡ G | λ=1 − − G| λ=0 between systems at λ = 0 and λ = 1 − under constant NpT was related to the difference in the surface interfacial energies as where the vacuum phase was denoted by subscript "0" and γ S0 and γ L0 were the solid-vacuum and liquid-vacuum interfacial energies per unit area. Note that γ L0 was substituted by the liquid-vapor interfacial tension γ LV in the final approximation considering that the vapor density was negligibly small. The work of adhesion W SL was defined by the minimum work needed to strip the liquid from the solid surface per area under constant NpT.
Using the NpT canonical ensemble associated with the Gibbs free energy G, the difference of the Gibbs free energy ∆G in Eq. (11) was calculated through the following TI: where H and V are the Hamiltonian, i.e., internal energy of the system and system volume, respectively, and N w is the number of wall molecules. The ensemble average was substituted by the time average in the simulation and was denoted by the angular brackets.
On the other hand, for the case of SV-DS in Fig. 1(c), we simulated constant-temperature equilibrium systems of a fixed volume with a particle bath located on the top for coupling parameter λ ranging from 0 to 1 − . The SV interface in the SV-DS system at λ = 0 in the left panel of Fig. 1(c) was in equilibrium with a saturated vapor at this temperature and this approximately represented the SV interface away from the contact line in the droplet system with the corresponding η value. Similar to the SL-DS systems, the SV interface at λ = 0 was divided into S0 and V0 interfaces at λ = 1 − as shown in Fig. 1(c), while the SV-DS systems were under constant NVT. Thus, the solid-vapor work of adhesion W SV was given by the difference of the Helmholtz free energy ∆F per unit area and was related to the difference in the surface interfacial energy as where γ V0 was set zero in the final approximation.
Using the NVT canonical ensemble, ∆F in Eq. (13) was calculated through the TI as

III. RESULTS AND DISCUSSION
A. Force balance on a control volume around a contact line Figure 2 displays the time-averaged density distribution of an equilibrium argon droplet in the case of η = 0.4, where the origin of the x-axis was set relative to the droplet center of mass. As well-known, multiple adsorption layers with a thickness about 1-2 nm were formed at the solid-liquid interface, and a hemi-cylindrical liquid-vapor interface with a uniform curvature was observed above. This indicated that the liquid-vapor interfacial tension was uniform except around the contact line. We supposed a rectangular control volume set around the contact line with its bottom and top faces parallel to the solid surface and its vertical left face in the center of mass of the droplet as shown in Fig. 2 and examined the horizontal force balance on this control volume through the fluid stress tensor field τ(x, z). Before examining the details of the control volume, we will look at the fluid stress on x = 0. Figure 3 shows the distributions of the fluid stress components τ xx (0, z) and τ zz (0, z) in the center of the droplet, i.e., τ(x, z)| x=0 , superimposed with the distribution of density ρ(0, z). With respect to the relation between ρ and τ xx , the former fluctuated near the solid due to the multi-layered structure, and τ xx also fluctuated corresponding to the density: negative in the high density layers and positive in between. These were because the fluid was compressed in the lateral directions inside the adsorption layers and attractive force lines acting between fluid molecules in neighbouring adsorption layers crossed the x-normal bin faces in between. Above these adsorption layers, ρ and τ were both constant, indicating that a liquid bulk was formed. Note that the negative value of τ xx was due to the Laplace pressure. In contrast to the gradual decrease of density from liquid to vapor, τ xx showed a remarkable peak structure in the transition layer between liquid and vapor, which corresponded to the surface tension. A saturated vapor bulk with constant density and τ xx existed above the liquid-vapor interface.
On the other hand, τ zz was only equal to τ xx in the liquid and vapor bulks, i.e., the fluid stress was isotropic. In contrast to τ xx which fluctuated in the adsorption layers, τ zz was almost constant. This was because the strong external force from the solid was exerted only near the solid surface, and the force balance normal to the solid surface was satisfied only with τ zz since the contribution from τ xz cancelled out due to the symmetric feature of this area. Indeed, based on linear momentum conservation in a steady state system without any flow perpendicular to the wall, ∂τ zz /∂z = 0 is satisfied in case without external force. An equivalent explanation for this is that the oscillations in the kinetic and interaction stress profiles cancel out each other for τ zz . This was not the case near the solid surface where the external force from the solid was included in the force balance. 41 In other words, the pressure exerted on the fluids there from the liquid bulk, i.e., negative of surfacenormal fluid stress, balanced the external force from the solid. Throughout the liquid-vapor transition layer, τ zz showed a simple change between constant bulk values, which was also in contrast to τ xx .
Considering this feature of density and fluid stress distributions, we defined the position of the bottom and top faces of the control volume in the z-direction in Fig. 2 as follows: the former denoted by z SA was set at the boundary of solid and first adsorption layer. More specifically, it was at the position where the liquid density showed steep rise from zero as shown in Fig. 3, corresponding to the nearest limit to the solid wall the fluid molecules could reach. On the other hand, we set the top face position z AL above the top adsorption layer and below the liquid bulk, as displayed in Fig. 3. Different from z SA , z AL cannot be strictly determined. Actually, we only need strict definition of z SA , whereas z AL may shift into the liquid bulk. This will be mentioned in Sec. III C with the mechanical definition of the contact angle related to the extended Bakker's equation. In addition, we set z LT between the liquid bulk and LV-transition layer. The position of the right face of the control volume was set as a parameter and was denoted by x R . Now, we suppose an arbitrary volume V containing only fluid molecules bounded by a closed surface S in the present quasi-two-dimensional droplet system. Exact force balance must be satisfied for the present quasi-static system, and it follows for the surface integral of the fluid stress and volume integral of the external force that where the Einstein notation is used with a dummy index j, and τ ji (x, z), n j , and f ext i (x, z) are the fluid stress tensor, unit normal vector of the surface element dS, and external force vector per volume from the solid, respectively. By applying Eq. (15) to the rectangular control volume shown in Fig. 2 with bottom, top, and left positions fixed at z = z SA , z = z AL , and x = 0, respectively, the horizontal component of the LHS of Eq. (15) was given as a function of only x R , and the first term is written as where the three terms in the RHS denoted the fluid stress surface integrals on the right, left, and top faces, respectively. Note that the integral ∫ x R 0 τ zx (x, z SA )dx was omitted because τ zx (x, z SA ) was zero in the entire region of the bottom face because no fluid molecules existed under this plane to contribute to the sum in Eqs. (8) and (9). For short notation, we define the horizontal components of the fluid stress surface integrals in Eq. (16) by and as a function of x, and rewrite Eq. (16) as   (20) was satisfied in the entire range of x R even though T xx (x R ) and T zx (x R ) changed depending on x R . By inserting this result into Eq. (15), it followed for the second term of the LHS that meaning that the horizontal component of the external force was zero irrespective of x R , i.e., no horizontal force was exerted on the fluid from the solid surface with a densely packed crystal structure everywhere as a time-average for the present droplet system under random Brownian motion even though the fluid density was not uniform. Note that the LHS of Eq. (20) was not exactly zero, and this seemed to be due to the slight roughness of the present solid surface. In the following discussion, we set and for simplicity.
With respect to T zx (x R ), its value was zero up to around x R = 1 nm, and this indicated that τ zx (x, z AL ) was equal to zero because the fluid stress in the liquid was isotropic up to the liquid-vapor transition layer. After the decrease of T zx (x R ) with the increase of x R , it took a constant value above about x R = 4 nm because the fluid stress was isotropic with τ zx (x, z AL ) = 0 in the vapor bulk far from the transition layer.
Along with the change of T zx (x R ), the fluid stress integral T xx (x R ) was constant except in the liquid-vapor transition layer while satisfying Eq. (22) in the entire region as mentioned above.
Considering this feature, we defined the contact line region as x LC ≤ x R ≤ x CV between solid-liquid and solid-vapor interfaces in which T zx (x R ) and T xx (x R ) changed, as shown in

B. Bakker's equation extended to solid-related interfaces
As displayed in Fig. 3, the fluid stress tensor is not isotropic at phase interfaces. Bakker's equation describes the relation between the fluid stress anisotropy and surface tension through a thought experiment shown in Fig. 5(a). 34 In this thought experiment, one piston is set normal to a flat

ARTICLE
scitation.org/journal/jcp liquid-vapor interface and it covers from z = z L to z = z V at liquid and vapor bulk regions, respectively, across the plane of the liquid-vapor interface. Another piston parallel to the interface is set in vapor bulk far from the interface. Through simultaneous virtual infinitesimal displacements of the pistons, only the interface area can be changed without changing the vapor and liquid volumes. Let l be the depth normal to the xz-plane, and −δV and δx be the infinitesimal volume change given by the downward displacement of the top piston and corresponding displacement of the side piston, respectively, it follows that Assuming that this displacement is done quasi-statically under constant temperature, the minimum mechanical work δW required for this change is associated with the change in the Helmholtz energy F given by where τ zz and τ xx are the surface normal fluid stress components at the surfaces of top and side pistons, respectively. The former is constant in the entire region because of the force balance in the z-direction. We denote this constant value by τ bulk , which is equal to the saturated vapor pressure with its sign inverted. On the other hand, the latter is a function of z, while it is also equal to τ bulk in the liquid and vapor bulks. By equating δF with surface free energy of area lδx, and eliminating δV, l, and δx with Eq. (24), it follows for the surface tension γ LV that where A LV is the area of the liquid-vapor interface. The integration range in Eq. (26) is sometimes expressed by ∫ ∞ −∞ ; however, the necessary condition for this range is that z L and z V cover the entire range of the interface with anisotropic fluid stress, i.e., between liquid and vapor bulks.
We extend this thought experiment to the solid-liquid interface as shown in Fig. 5(b), where the top piston parallel to the interface is set away from the interface, whereas the side piston is set not across the plane of the interface but just on it so that the piston face can only contact the liquid. Then, one can achieve similar simultaneous infinitesimal displacements of the pistons without changing the liquid volume. However, the resulting change in the interface is not a simple increase of the SL interface area but a replacement of the solid-vacuum interface with the SL interface. Hence, the relations among the minimum mechanical work, change in the free energy and interfacial tensions are expressed by Note that the lower limit of the integration range should be defined physically at the limit nearest to the solid which the fluid can reach, i.e., exactly the same as z SA microscopically defined in Fig. 3.
If we consider another thought experiment system of the solid-vapor interface similar to Fig. 5(b), the following relation is derived: The difference from γ S0 described in Eqs. (27) and (28) was called "relative" interfacial tensions in our previous studies. 20,27 C. Relation between Young's equation and extended Bakker's equation Now we go back again to the horizontal force balance exerted on the control volume around the contact line. From   Fig. 4, it was known that T xx (x R ) and T zx (x R ) were both constant for 0 ≤ x R ≤ x LC and x R ≥ x CV away from the contact line. Thus, from Eqs. (17)- (19) and (22), it follows that In order to think about the third term of Eq. (29), at first, we suppose an imaginary equilibrium cylindrical liquid pillar having the same radius of curvature R as the LV interface of the droplet on the solid surface as shown in Fig. 6 and consider the horizontal force balance on a rectangular control volume whose bottom face is on the x-axis at z = 0 depicted in blue. Then, the shear stress on the bottom face τ zx (x, 0) is zero because of the symmetry there, and it follows that where the 1st and 2nd terms are the horizontal force per unit depth exerted on the left and right faces expressed by constant external and internal pressures p ext and p int , whereas the 3rd term is that on the top faces. Note that p int and p ext are equal to the uniform normal fluid stress −τ xx in the liquid and The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp vapor bulk regions, respectively. Let θ be the angle between the top face and LV interface shown in Fig. 6, for which the geometric relation is given by By applying Eq. (31) and Young-Laplace equation it follows for Eq. (30) that indicating that the integral of the fluid shear stress on the top face is equivalent to the lateral force exerted by the LV surface tension. Now, we go back to Eq. (29). Let τ bulk V and τ bulk L be the isotropic fluid stress values in the vapor and liquid bulks equal to −p ext and −p int , respectively, and by once removing them from the integrand of the first and second terms of the LHS of Eq. (29), and by adding their integrals, it follows that Note that when Eq. (23) is satisfied, Eq. (29) would hold even if the lower limit of the integral for the 1st and 2nd term of LHS was taken at an arbitrary position below z SA because τ xx (x, z) = 0 for z < z SA . Even in this case, Eq. (34) is also recovered, as shown in Appendix A. By substituting the 1st and 2nd term with the left-most hand sides of Eqs. (28) and (27), respectively, Eq. (34) is rewritten as indicating that the interfacial tensions are balanced on the control volume including the Laplace pressure. Finally, let θ be the angle between the bottom face and extended cylindrical plane of the LV interface satisfying as shown in Fig. 6, and by using Eq. (31) and Young-Laplace equation (32), Eq. (35) recovers Young's equation (1).
The following three points used to derive this consequence should be rephrased here. One point is that this result is derived because Eq. (23) is satisfied on the present flat and smooth solid surface as a prerequisite. In other words, this equivalence of Young's equation to the pure mechanical balance is true only when Eq. (23) holds, and it is not the case for a contact line subject to pinning due to physical or chemical inhomogeneity, e.g., surface roughness, edge structures, impurities, or surface deformation also seen in our case in Fig. 4. This condition of Eq. (23) would easily be violated for a water droplet on an inhomogeneous polar solid surface because the average dipole direction of water molecules at the contact line would be biased from the surface normal direction. 22,24,[31][32][33]36 Note that Kanduč and Netz 31 showed a good agreement for contact angles of water on a homogeneous polar surface calculated from their interface potential method and droplet simulation. Extracting the lateral force around the contact line could provide further insight into the matter and will be the target of future research. The second point is that Eq. (29) holds even though the lower limit of the integral for the 1st and 2nd term of LHS is taken below z SA on condition that Eq. (23) is satisfied, whereas z SA must be strictly given by Eqs. (27) and (28), which determines γ SL and γ SV , respectively. Thus, the position of the fluid-wall interface to define the contact angle θ is also determined. This is because different values of isotropic fluid stress values τ bulk V and τ bulk L , whose difference corresponds to the Laplace pressure, are included to derive γ SL and γ SV in Eq. (34). This is not the case for a contact line with the flat LV interface for which the fluid-wall interface does not have to be given to define the contact angle. We show in Appendix B that Eq. (33) holds also in this case with θ = θ. Note also that the fluid-wall interface position z SA may include ambiguity in case the surface is not flat. The third point is that the top face position set at z AL in Eq. (29) is not an essential prerequisite, and it may be varied between z AL and z LT . In other words, Young's equation should not be considered as the horizontal force balance on a point as in Eq. (1), but on a finite volume as in Eq. (35). 43

D. Relation between extended Bakker's equations and work of adhesion
We have extended Bakker's equation to SL and SV interfacial tensions and examined Young's equation from a mechanical point of view above. In Fig. 7, we compared these interfacial tensions obtained from a mechanical route with the corresponding work of adhesion W SL and W SV obtained from a thermodynamic route using the TI. Note that the error-bar for W SV became larger for larger η because of the fluctuating behavior of the solid-vapor adsorption layer. Considering the difference of the definitions, relative interfacial tensions γ SL − γ S0 and γ SV − γ S0 are shown as −(γ SL − γ S0 ) + γ LV and −(γ SV − γ S0 ), respectively, where γ LV = 11.3 × 10 −3 N/m is added to −(γ SL − γ S0 ). The value of γ LV was obtained from a standard simulation system with a planer liquid-vapor interface. 20,27 The interfacial tensions −(γ SL − γ S0 ) and −(γ SV − γ S0 ) showed the same dependence on η as the work of adhesion W SL and W SV , respectively, indicating that the interfacial tensions obtained by the mechanical and thermodynamic routes corresponded well for the present case of the flat and smooth fluid-wall interface on a densely-packed crystal surface. The good correspondence between −(γ SL − γ S0 ) and W SL also indicated that these two calculated by two routes were differed exactly by γ LV . Note that this correspondence was also because the Laplace pressure in the droplet was comparatively small so that the pressure dependence of γ SL was negligibly small. In addition, for wettable cases of η > 0.4, −(γ SV − γ S0 ) or W SV had a non-negligible value, indicating that the solid-vapor interfacial tension should be considered to properly estimate the contact angle.  7. Comparison between relative interfacial tensions obtained from the mechanical route and work of adhesion for flat SL and SV interfaces obtained by the DS method as a thermodynamic route. Considering the difference of the definitions, relative interfacial tensions γ SL − γ S0 and γ SV − γ S0 are shown as −(γ SL − γ S0 ) + γ LV and −(γ SV − γ S0 ), respectively, where γ LV = 11.3 × 10 −3 N/m is added to −(γ SL − γ S0 ). The value of γ LV was obtained from a standard simulation system with a planer liquid-vapor interface. 20,27 The error bars were obtained from the standard deviation. Figure 8 shows the comparison among three contact angles (1) estimated by using the interfacial tensions γ SL and γ SV obtained through the mechanical route, (2) estimated by

FIG. 8.
Comparison between the apparent contact angle obtained from the droplet system and expected contact angle from the interfacial tensions obtained through mechanical and thermodynamic routes.
using the SL and SV works of adhesion W SL and W SV , respectively, and (3) directly measured from the apparent shape in the droplet system in which the angle was defined by the angle between the least-squares fitting circle to a density contour line of ρ = 400 kg/m 3 for z > z AL and the plane of z = z SA at the intersection, as indicated in Fig. 6. 20 Note that we could give a strict choice of the density value for the contour as well, 44 but the result was not sensitive to the choice because the LV-interface was relatively thin. The value of γ LV = 11.3 × 10 3 N/m 20 was also used as well. It was shown that the three contact angles agreed quantitatively well for the present case. The correspondence between the mechanical route and apparent contact angle was achieved because the solid-fluid interface was properly given, which was used for the range of the lateral fluid stress integration. At present, we have not identified a clear reason for the difference between the mechanical route and droplet shape for a lower η of 0.1 yet, and further analysis on the following possible causes is needed: (1) the rather arbitrary choice of the density value of ρ = 400 kg/m 3 for the shape fitting, (2) spatial resolution of the stress calculation which directly affects the stress integral in Eq. (17), and (3) assumption of γ LV to be constant irrespective of the curvature.

E. Discussion
The surface tensions calculated from fluid stress tensors were well converged in the present work for the LJ system on a rather hard, flat, and smooth solid crystal surface. A question arises whether this may or may not be the case for molecular solids because the residual stress and internal dynamics in the solid usually cause problem. 45,46 The effects of the change of solid volume were explicitly taken into account in the TI scheme Kanduč and Netz. 31 For the mechanical path, since Eq. (15) is exact, it follows without assuming the zero-external force condition in Eq. (23) and using Eq. (33) that where V denotes the range of x LC ≤ x ≤ x CV and z bot ≤ z ≤ z AL , as long as we take the height of the bottom face of the control volume z bot well below the solid-fluid interface so that is satisfied. Then, the effects of the residual stress and the internal dynamics as well as the surface deformation 45,46 can be embedded into the external force in Eq. (39). We do not expect difficulties with the convergence of the fluid stress integrals and volume integral of the external force in Eq. (39). Ambiguity could arise when we transform the fluid stress integrals in the LHS of Eq. (39) by Eqs. (27) and (28), in which the The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp bottom position must be defined for τ bulk 0. We think that the definition becomes difficult for relatively soft molecular solids, where we would have to introduce different z SL and z SV due to the solid deformation, for instance. Regarding the mechanical path, we have only examined the force balance tangential to the surface, but it is also technically trivial to deal with the surface-normal force balance using the present fluid stress. Further analysis on this surface normal force would give valuable information about γ LV sin θ which is related to the disjoining pressure. [45][46][47] With respect to the comparison between mechanical and thermodynamic paths, since the SL-DS simulations were performed in the NPT ensemble, the resulting W SL was the work of adhesion at the prescribed pressure at 1 MPa rather than at the coexistence condition. On the other hand, the SV-DS simulations were under coexistence condition with the particle bath. These two mean that both the resulting W SL and W SV are under conditions slightly different from the coexistence condition in the droplet simulations. However, W SL was not sensitive to the pressure even within the present range of Laplace pressure of up to a few MPa, as shown in Fig. 3. For W SV , its absolute value was small compared to W SL , and the effect of vapor pressure on the evaluation of Young's equation was considered to be negligible. These are because we have used a relatively hard solid crystal surface without polarity in the present study. Indeed, non-negligible effect of vapor pressure was shown by Kanduč and Netz 31 for a system with water on a polar solid surface. Related to these points, the solid-vacuum interfacial tension γ s0 is included in Eqs. (11) and (13), and the pressure effect on γ s0 should be included for a system with a soft solid subject to considerable deformation due to pressure; however, this is not the case for the present system.
Finally, a further analysis on the system size dependence would give new insight with regard to the finite-size scaling analysis on the wetting transition. 18 Although we did not test the size dependence of the solid-related interfacial tensions, the curvature dependence of the liquid-vapor interfacial tension γ LV was small for the size range of the present droplets. 20,44

IV. CONCLUSIONS
We carried out MD simulations of a cylindrical Lennard-Jones droplet on a flat and smooth solid surface. We examined in detail the surface-lateral force balance on a control volume set around the contact line and theoretically connected the fluid stress integral along each face of the control volume with the corresponding interfacial tension. Through this connection, we showed that Young's equation was only applicable to a system in which the net lateral force exerted on the fluid molecules from the solid surface was zero around the contact line. This condition would be easily violated for systems with realistic molecular fluids and/or the solid surface terminated with functional groups. Furthermore, we showed that the solid-liquid and solid-vapor interfacial tensions relative to the vacuum-solid one obtained by mechanical and thermodynamic approaches through an extended Bakker's equation and DS method, respectively, were in good agreement. In addition, the contact angle predicted by Young's equation with these interfacial tensions corresponded well to the apparent droplet contact angle determined by using the previously defined position of the solidfluid interface plane; however, our theoretical derivation indicated that this correspondence was achieved because the zero-lateral force condition was satisfied in the present system with a flat and smooth solid surface. These results indicated that the pinning force density should be additionally included in Young's equation to predict the contact angle.