The Intrinsic Fragility of the Liquid-Vapor Interface: A Stress Network Perspective.

The evolution of the liquid-vapor interface of a Lennard-Jones fluid is examined with molecular dynamics simulations using the intrinsic sampling method. Results suggest clear damping of the intrinsic profiles with increasing temperature. Investigating the surface stress distribution, we have identified a linear variation of the space-filling nature (fractal dimension) of the stress clusters at the intrinsic surface with increasing surface tension or, equivalently, with decreasing temperature. A percolation analysis of these stress networks indicates that the stress field is more disjointed at higher temperatures. This leads to more fragile (or poorly connected) interfaces which result in a reduction in surface tension.

Introduction

The intricate dynamics of the liquid–vapor interface have always attracted researchers from diverse research fields for its prevalence in myriad natural phenomena, such as biolocomotion[1−3] and plastron respiration,[4] and in numerous technological processes ranging from efficient oil recovery[5] to organic electronics[6] to capillarity driven thermal management.[7,8] To characterize the liquid–vapor interface, a number of approaches have been developed over the past century, from macroscale treatment[9] to sophisticated molecular dynamics (MD) simulations,[10] which have changed our perception of the interface quite dramatically.[11]

Among various available MD techniques, the mechanical route to surface tension measurement has gained popularity for its relative simplicity and accuracy. This method replaces the scalar pressure with a second-rank pressure tensor[12] and, thereby, gives access to the atomic/molecular level detail of the interfacial region through the microscopic definition of the pressure. Following the Irving and Kirkwood (1950) definition,[13] the diagonal elements of the pressure tensor can be expressed aswhere the momentum = m( – u), is the total particle velocity and u is its streaming part, δ(···) is the Dirac delta function, O is the Irving–Kirkwood operator,[14] is the force exerted by atom j on atom i, and = – is the central distance between the two interacting atoms. The average values of the off-diagonal elements in the pressure tensor are zero due to the axial symmetry about the normal direction and the translational symmetry about the plane parallel to the interface. With the normal (PN) and the tangential (PT) components of the pressure, surface tension can be obtained from the Hulshof integral,[15] as in eq :

The mechanical stability of the interface requires the normal component of the pressure tensor to be uniform and constant throughout the bulk phases and at the interface. On the contrary, the tangential component strongly depends on the position vector in the neighborhood of the interface, and only at a location far from the interface, this becomes uniform and equal to the normal component. Although eq requires integration over the entire domain of the two-phase system, the term inside the integration soon becomes zero as one moves a few atomic diameters from the interface. We will further elucidate this in Results and Discussion.

Attaining a molecular/atomic description of the interface, which is essential to understand the microscopic origin of surface tension,[16] is challenging due to the inhomogeneous nature of the interface and the difficulties in quantifying the surface forces. The non-uniformity of the particle number density across the interface gives rise to mathematical complexity in uniquely identifying the system and its thermodynamic properties.[17] In past years, numerous efforts have been made, both experimentally[18,19] and by computer simulations,[20−23] to advance our understanding of the interface. A comprehensive review can be found in Ghoufi et al.[10] and Ghoufi and Malfreyt.[24] Most MD simulation studies in the literature focused on the average profiles of the interface where the thermal fluctuation effects obscure the identification of the intrinsic behaviors.[21,25,26] The presence of thermal fluctuations, in the form of capillary waves, blurs the interfacial properties and makes it extremely difficult to extract the true nature of the interface.[25] Therefore, it is necessary to decouple the effects of capillary wave from the surface layer in order to investigate the intrinsic nature of the interface—an idea first introduced by Buff et al.[27] In their capillary wave theory, Buff et al. formalized the concept of the existence of an interface ζ(y, z, t) which acts as an instantaneous atomic border between the two phases, where the vector quantity, y, and z are parallel to the interface. Sides et al.[28] exploited this idea of decoupling the capillary wave broadening from the intrinsic interface and showed that the surface tension measured from eq shows good agreement with that measured from the interface width. However, to achieve a meaningful perspective of the intrinsic interface, adequate care must be taken in selecting the appropriate system size and the transverse resolution to avoid any blurring of the footprints of the intrinsic layer by the capillary wave fluctuations.[29] Following its development, the intrinsic sampling method has been applied to understand the interface layer in a number of studies, i.e., to investigate the intrinsic gap in the water–oil surface[25] and to define the thickness of an adsorbed liquid film on a substrate.[30]

Given the strong dependency of the surface fluctuations on temperature, it is surprising that more attention has not been devoted to understand the effect of temperature variation on the intrinsic profiles of the liquid–vapor interface. Not only that temperature modifies the surface properties—often it becomes one of the key driving factors for a number of interfacial events, e.g., thermocapillary flows.[31] In the present study, we carefully examine the effects of temperature on the intrinsic density and pressure profiles. Afterward, we apply fractal analysis and the concept of percolating networks to analyze the temperature dependency of the spatial correlation of the interatomic interactions at the surface layer.

Methodology

MD simulations with the Flowmol MD code[32] were used to model the liquid–vapor interface of a Lennard-Jones (LJ)[33] fluid at different temperatures. The initial simulation domain was a cubic box of dimensions L = 120.64, L = 19.05, and L = 19.05 in reduced LJ units, containing a total of 14 827 particles. The middle 40% of the box was initialized with LJ particles in liquid phase, and the remaining was designated as vapor. The initial state was created from a face-centered-cubic (fcc) lattice with a starting density of ρ = 1, and then molecules were randomly deleted until the preset liquid (ρ = 0.5) and vapor (ρ = 0.005) densities were obtained. Note, these initial units were somewhat arbitrary and used to set up the system only; a sufficient equilibration time was then allowed so that the system tended to the expected coexistence (see Table S1 and ref (34)). The equilibration phase was run in the canonical (NVT) ensemble, using a Nosé–Hoover thermostat,[35,36] for 50 000 time steps with ΔT = 0.005. A shifted LJ potential was used with periodic boundary conditions applied in all three Cartesian directions, and the Verlet Leapfrog[37] integrator was used to integrate the equations of motion.

For surface tension calculation, a small cutoff radius, rc, results in significant deviation from experimental data for argon,[38] whereas rc ≥ 4 is found to give better agreement.[39] Surface tension being one of our prime interests in this study, and as recommended by Ghoufi et al.,[10]rc = 4.5 was used to improve accuracy. Further agreement with experimental data would require a more complex interaction potential, such as those which explicitly include three-body interactions.[10] The final state of the initial NVT calculations was taken as an initial condition for production runs in the microcanonical (NVE) ensemble. We used three independent simulations with different initial conditions to improve the statistics.

The intrinsic interface was fitted to the outermost layer of the liquid by cluster analysis with the Stillinger cutoff length,[40]rd = 1.5, and with the required criterion of each atom having more than three neighbors.[41] Instead of assuming an average interaction contour, the functional form of the interface, ξ(y, z, t), was refitted at every time step by using the intrinsic sampling method (ISM),[42] which approximates the liquid–vapor interface by means of a Fourier series representation as in eq .

Here, k is the wave vector that corresponds to the periodic boundary conditions, i.e., k = 2π(n/L, n/L) with n, n = 0, ±1, ±2, ..., k is the modulus of the wave vector, and ∥ denotes the parallel surface components in the y and z directions. The fitted coefficients, are expressed as time dependent functions because they are refitted to the surface each time the positions of the surface atoms change (i.e., every time step). Figure a is representative of the coexisting system, where an intrinsic interface was fitted to the outermost liquid layer; the interface and the fluctuations are illustrated in Figure b. All results presented hereafter correspond to the interface at the right side of the domain; x > 0 corresponds to the vapor side and x < 0 corresponds to the liquid side.

Figure 1

(a) Liquid–vapor coexistence. (b) Fitted liquid–vapor surface function ξ (yz plane). (c) Mapped grids at the interface, x = 0 (mapped coordinates are in red at the top, and unmapped coordinates are in black at the bottom). The atoms are colored by their respective velocities; the bin resolution is reduced here for visualization purpose. The mapping shifts every time the interface moves in the x direction.

Once we have a mathematical form of the intrinsic interface, the density and pressure tensor can be obtained in a reference frame that moves with the surface, ξ, which has a different value at every point in y and z updated at every time step, t. To obtain quantities that move with the interface, the Irving and Kirkwood[13] definitions can be integrated over a volume where the surfaces in the x direction follow the function given by eq , with uniform grids in the y and z directions. The density in a volume moving with the interface is thenwhere ϑ is the product of Heaviside functions, which is unity when an atom is inside a volume and zero if outside.[41] This can be expressed in terms of the product ϑ = Λ(x) Λ(y) Λ(z), where the difference between two Heaviside functions is known as a boxcar function and is represented by Λ. This is formally obtained from integrating the Dirac delta function between finite limits but can be simply expressed in terms of an indicator function as follows:which checks whether x is between the two limits, a– and a+. In the y and z direction, the indicator functions Λ(y) and Λ(z) check whether the atom positions {y, z} are between the spatial extents of a surface bounded by y– and y+ and z– and z+, respectively.

For the x direction, the position of the interface is included in these limits, so for a box centered on x with size Δx the limits are x± = x ± Δx/2 + ξ(y, z). The easiest way to obtain the density relative to the intrinsic surface is to map the atomic positions based on the intrinsic surface at the same y and z location, x′ = x – ξ(y,z) and then simply bin as if using a uniform grid, i.e., bin = round(x′/Δx). The grid is uniform in the other two directions, so the binning in y and z is unchanged; see Figure c. Quantities such as momentum, temperature, and pressure inside a volume[43] can also be obtained using the same mapping approach, where the latter requires mapping of the line of interaction between atoms to get the configurational term.[41] However, only the pressure tensor obtained from surface fluxes, taken here over the surfaces of a binning volume, can be shown to satisfy the mechanical equilibrium condition; i.e., the normal pressure is constant when moving through the interface.[44] This form of pressure also includes a term for the movement of the intrinsic surface itself, which must be accounted for to ensure balance of the coarse-grained equations. The surfaces of a volume are flat in the y and z directions and follow the interface in the x direction. The surface pressure can be written as the sum of three contributions:where α ∈ {x, y, z} denotes the three directions where each is a vector α = [P, P, P]T, with the kinetic contribution αkin. that comes from the momentum transport of atoms crossing the interface,[45] the configurational contribution αconfig. that arises from the atomic interaction,[12] and αSM., a term due to the surface movement in time.

Introducing a surface normal vector, which for the flat surfaces in y and z directions are simply the standard basis unit vectors = and = . For the intrinsic surface, it is given by = ∇(ξ – x)/∥∇(ξ – x)∥, where the subscript s denotes the derivative taken at the point of crossing x. By assuming the convective term to be zero, ρ = 0, the three pressure contributions discussed above can be obtained in a molecular simulation as follows:where ΔSα is the surface area, = (t2) – (t1) denotes the vector for the movement of an atom from time t1 to t2, and is the separation vector. The ϑ term captures the movement of the surface by counting all atoms which enter or leave a volume as the intrinsic interface itself moves in time. Defined as ϑ = Λξ(x) Λ(y) Λ(z), the surface evolution from the start of a time step ξ– = ξ(t) to the end ξ+ = ξ(t + Δt) is multiplied by an indicator function, i.e., eq .

The dSα term in eq is the derivative of the ϑ function, with respect to α = {x, y, z}, and is only nonzero if the separation vector, or , is crossing the surface of the volume in question. Without loss of generality, we consider the expression for the surface and the separation vector in the x direction aswhere λ parametrizes the line between atoms λ = + λ with λ the value at a point of crossing on a surface. The function Λλ therefore checks if a λ value is between (λ– = 0) and (λ+ = 1), with the remaining functions checking if the point of crossing {yλ(λ), zλ(λ)} is on the volume surface between y– and y+ and z– and z+. The forms of dS and dS are similar, and atomic motions are parametrized in the same way with λ = 1 + λ.

The calculation of the pressure tensor, therefore, becomes a ray-tracing problem, namely, getting all the intersections of the surfaces of a volume due to interatomic interactions and atomic crossings. In order to accelerate the process of getting intersections on an intrinsic surface, for each binning volume, the Fourier function of eq is converted to a set of bilinear patches of the formThe intersection of a line and the bilinear patch is a local operation, which is much quicker than the root finding process on a full Fourier surface of eq . The procedures to fit the intrinsic interface and to choose the number of bilinear patches, as well as to calculate the intersect, are described in previous works.[41,44]

Results and Discussion

Temperature Effects on Density Profiles

The intrinsic density profile at temperature T = 0.7 is shown in Figure a. Considerable oscillations are evident near the intrinsic surface layer extending at least five atomic diameters into the liquid phase (density is averaged over time, as well as in the y and z spatial directions). Although such layering is a universal behavior of the free surface,[46] these oscillations are often smoothed out for an LJ fluid if the reference frame is static. On the other hand, a moving frame of reference makes the layering evident. This layering directly determines the stress that would be measured, both on the interface itself and in the bulk fluid. In the next sections, we focus on the effects of such layering on the surface stress and, thereby, on the surface tension.

Figure 2

(a) Density profile (blue circles) and derivative of density (black dashed line) for T = 0.70 showing oscillations near the liquid–vapor interface. (b) Density profiles for a range of temperatures; oscillations are damped as temperature increases.

As seen from Figure a, the amplitude of the oscillation increases as the surface is approached from the liquid side; the highest peak is attained at the intrinsic surface and quickly dampens within a distance of less than half an atomic diameter into the vapor side. The first derivative of the density profile, dρ/dx, further illustrates these fluctuations, and the zero crossing highlights the position of the intrinsic layer. Such derivatives are also instructive of the interfacial widths; readers are referred to Sides et al.[28] for further details. Figure b, where the intrinsic density profiles are plotted for a range of temperatures, shows increased damping at higher temperatures. Evidently, the fluctuations only exist at the liquid side and dampen quite abruptly in the vapor side right after the interface (x > 0), despite the existence of a small peak due to the adsorbed layer.

The oscillations of the density profile at the liquid side for x ∈ (−∞, −1] can be approximated, as illustrated in Figure a, by an exponential decay function of the formwhere ρ(x) denotes the local intrinsic density, ρbulk is the liquid density of the bulk phase far from the interface, b is the damping coefficient, and f is the frequency of oscillation. Note that the negative sign of the exponent, i.e., the damping coefficient in eq , is discarded as the location, x, in the liquid side is represented with negative numbers. The damping coefficient linearly increases, whereas the frequency of oscillation decreases as temperature is increased from T = 0.70 to 1.0; see Figure b. Both the damping coefficient, b, and the frequency of the oscillation, f, show linear fits with temperature within the considered range, i.e., b = 2.51T – 0.81 with R2 = 0.95 and f = −0.7πT + 2.5π with R2 = 0.99. Thus, the intrinsic density in eq can be expressed as a function of location, bulk density of the liquid, and the temperature, i.e.

Figure 3

(a) Oscillation of intrinsic density profiles. Circles denote results from MD simulations; solid lines are exponential fitting, as in eq , to the density profiles at different temperatures. (b) The damping coefficient (left y axis, black filled circles) increases with temperature, whereas the frequency of oscillation (right y axis, blue filled squares) decreases. The solid lines are linear fits to the data.

While a liquid–vapor interface is not directly comparable to a solid–liquid interface, one should obtain similar damping for the same bulk liquid in contact with a solid wall.[47]

Temperature Effects on Pressure Profiles

The oscillatory nature of the intrinsic profiles can be further realized from the pressure profiles. Both the normal and the tangential components of the pressure tensor can be decomposed into their constituent parts. Figure a shows these various components (averaged over time and yz plane) of normal pressure for T = 0.70. To track the movement of the interface and its temporal evolution at a resolution of atomic spacing, Smith and Braga[41] and Smith[44] discarded the otherwise applied concept of an average interaction contour and introduced an instantaneous frame of reference which evolves with the interface and hence describes the pressure tensor in a purely mechanical manner. The consideration of the surface movement contributes to an additional corrective term, PSM. This, when considered along with the configurational and the kinetic contributions of the normal pressure, as shown in Figure a, exactly balances the momentum change to machine precision.

Figure 4

Pressure profiles at and near the surface layer for T = 0.70. (a) Components of the normal pressure: (blue filled squares with lines) configurational part, Pconfig; (black filled squares with lines) kinetic part, Pkin.; (green filled squares with lines) contributions from surface movement, PSM; and (red line) total normal pressure, PTotal. (b) Normal pressure without (blue dashed line) and with (red solid line) the consideration of surface movement term; ls, denotes the effective fluctuating length, over which the surface movements are visible. (c) Linear increase of the effective fluctuating length, ls, with temperature. The solid line is a linear fit to the data with a slope of 2.04, an intercept of 0.23, and R2 = 0.97.

The essence of the surface movement contribution is illustrated in Figure b, where the normal pressure is plotted without and with the consideration of the pressure correction. Only when the surface movement effects are accounted for, the kinetic and the configurational components can precisely balance each other resulting in a perfectly flat profile for the normal pressure which signifies that the liquid–vapor interface is mechanically stable. It is apparent from Figure b that the pressure contribution due to surface movements fluctuates only over a length of few atomic diameters (where the dashed line shows fluctuations) and remains constant elsewhere. We denote this length as the effective fluctuating length, ls (schematically shown in Figure b). More formally, we define ls aswhich is the maximal distance between any two roots of the equation |P(x)| = ϵ, 0 < ϵ ≪ 1 is the amplitude of the very small, but finite oscillations at either side of the effective fluctuating length. In this study, we use a cutoff value ϵ = 0.05 in order to numerically determine ls from the fluctuations of P(x). As shown in Figure c, ls increases linearly with temperature. Although the interfacial thickness assumes approximately the size of an atom at a temperature away from the critical temperature,[48] the effective fluctuating length demonstrates that the kinetic contribution spans over at least a few atomic diameters. This is reasonably justified in light of the fact that a higher temperature corresponds to greater surface movements and, thereby, a (kinetically) thicker interface.[49]

Although not straightforward, a plausible link between the effective length, ls, and the shear viscosity can be inferred through the hydrodynamic description of the capillary waves. It is well-established[50−52] that the damping rate Γ of an overdamped capillary wave mode satisfies the asymptotic relation Γ ∼ qγ/μ, where q is the wavenumber and μ is the dynamic shear viscosity. At critical damping, the real part of the complex wave frequency is zero, i.e., Re(ω) = 0 and q ∼ γρ/μ2. If the propagation velocity of the density fluctuations, v0, is further assumed to be constant, it can then be postulated that ls scales inversely with Γ, i.e, Γ ∼ v0/ls ∼ γ2ρ/μ3. Crucially, this rearranges to givethereby providing a link between the dynamic shear viscosity and the effective fluctuating length. This is an interesting interpretation of the effect of viscosity on the fluctuations of the interface at the atomic scale. A detailed analysis to confirm this postulate is beyond the scope of this paper, but this relationship between ls and μ certainly prompts further investigation.

The effect of temperature on the intrinsic pressure profile has further been examined in Figure . Though the shape of the normal pressure profile remains unaltered irrespective of temperature, the oscillations of the configurational and the kinetic parts, in Figure a–d, are seen to smooth out as temperature increases. The normal and the tangential components of the pressure, along with (PN – PT), are plotted in Figure e–h for different temperatures. Interestingly, the tangential pressure profiles become less corrugated at higher temperatures; that is, the spatial oscillations or the expansion and contraction of the surface become less prominent. Such behavior can be ascribed to the more energetic interface at higher temperature, making the interface tracking difficult. The dyadic term in the kinetic pressure of eq is a particle property and thereby proportional to the intrinsic density through Pkin. = ρkBT, where kB is the Boltzmann constant. As such, these oscillations are directly comparable to those in the intrinsic density profiles, as in Figures and 3. The equal and opposite of this kinetic component is the configurational pressure, as shown in Figure a–d, which ensures equilibrium. The observations made in Figures c and 5 allow concluding that an increase in temperature dampens the oscillations of the intrinsic profiles and, at the same time, broadens the effective fluctuating length, ls, over which surface movement effects are present.

Figure 5

(a–d) Constituent parts of the (red solid line) normal component of pressure, i.e., (black filled squares with lines) the kinetic part and (blue filled squares with lines) the configurational part for (a) T = 0.70, (b) T = 0.80, (c) T = 0.90, and (d) T = 1.0. (e–h) Variation of (blue solid line) the tangential component of the pressure tensor and of (black solid line) PN – PT near the interface for T = 0.70–1.0.

The influences of the oscillations discussed are confined to a small interfacial region. Far from the liquid–vapor interface, both the normal pressure and the tangential pressure become equal and uniform. Hence, it is not surprising that only the region of a few atomic diameters from the interface (in both the liquid and the vapor sides) contributes to the surface tension. The term inside the integral of eq is plotted (black solid line) in Figure a along with the cumulative integral, i.e, the surface tension, γ. The magnitude of the surface tension is seen to reach a plateau right after the peak, within an atomic diameter in the vapor side. A careful examination of the intrinsic surface (i.e., x = 0) thus portrays its outright significance in determining the surface properties, as seen in Figure a, where the large negative peak of the pressure difference at x = 0 contributes significantly to the integral of eq . We compare, in Figure b, the surface tension evaluated at different temperatures with available experimental data for liquid argon,[53] results from molecular simulations with truncated and shifted potential with cutoff radii of 2.5[54,55] and 4.0,[39] and truncated potential with long-range corrections and cutoff radius of 3.0.[56] As seen in Figure b, the comparison suggests good agreement.

Figure 6

(a) Difference between normal and tangential pressures, PN – PT (black solid line), near the interface. The cumulative integral of PN – PT, i.e., the surface tension, γ(x) = ∫–10 [PN(x′) – PT(x′)] dx′, is shown by blue filled circles, where the lower limit of integration is kept fixed at x = −10 and the upper limit is gradually increased up to x = 2. (b) Surface tension from present modeling (black filled circles) shows good agreement with results from previous experiments and simulations. The error bars denote the standard deviation of surface tension as measured from three separate ensembles. The solid line is a linear fit (with slope = −2.26, intercept = 2.63, and R2 = 0.99) to the surface tension obtained from the present study.

Surface Fractals

The preceding sections discuss the effects of temperature on the intrinsic profiles, and on the corresponding values of the surface tension. How such alterations take place, however, remain unanswered until now. To investigate this, and because the configurational stresses are the sole contributors toward the surface tension, we examine the distribution of these stresses ignoring the kinetic components, at the intrinsic surface layer (x = 0), i.e.

Figure a shows an instantaneous configurational stress distribution map at the interface. The black filled squares correspond to the atomic locations which the intrinsic interface is fitted to, at that particular instant. In Figure b, only stresses that contribute to the tension of the surface (negative stresses) are shown, where a cluster spreading over the surface becomes apparent (also see Figures S1 and S2, Movie S1, and related discussions). These “nontrivial” networks and their variation with temperature are analyzed by means of fractal analysis.

Figure 7

Instantaneous stress map at the liquid–vapor interface for T = 0.80. (a) Color map of the configurational stress distribution with atomic positions overlaid (black filled squares). (b) Only negative configurational stresses. (c) Surface tension as a function of the fractal dimension of stress network at the surface layer over the range of temperatures considered in this study. Error bars (with a magnification factor of 2 for visualization purpose) denote the standard error of the mean (SEM) of fractal dimension; the solid line is a linear fit with R2 = 0.95.

To quantify the temperature dependency of the structural complexity of the stress networks, we measure their fractal dimensions, Φ. This dimension is often used to quantify the space-filling nature, heterogeneity, or self-similarity of surfaces, clusters, etc.[57−60] Φ is calculated through the Minkowski dimension,[61,62] which is often referred to as the box-counting dimension, whereby for nonoverlapping N boxes with sides ϵ,

The algorithm employed to calculate Φ consists of converting the network maps (as in Figure b) into binary images such that the stress networks are depicted by black pixels on a white background. For a given grid side, ϵ, the number of grids, N(ϵ), required to fill the projected surface area of the aggregate is counted and the grid is made increasingly finer at each subsequent iteration. The fractal dimension, Φ, is then obtained from the slope of log(N) vs log(ϵ). See the Supporting Information for further details.

With the variation of temperature, we have observed consistent behavior of the fractal dimension of the network comprised of the negative configurational stresses that lie within [Pconfig – σ, Pconfig], where Pconfig denotes the mean of all the negative configurational stresses and σ denotes the standard deviation. The fractal dimensions thus obtained are compared against the corresponding surface tensions in Figure c. For the temperature range considered here, Φ is seen to linearly correspond to the surface tension as γ = 9.86 Φ – 13.86, which is shown by the solid line. This reflects the greater space-filling nature[58] of the (surface) stress network at a lower temperature. At the same time, the fractal dimension of the stress clusters only at the outermost atomic layer proves to be predictive of the surface tension.

In hydrodynamics, the surface tension, γ, is often modeled as an equation of state in terms of temperature, surfactant concentration, etc. This seemingly mechanistic approach does not capture the subtleties of the mutual interactions between the various effects, nor does it directly model the fundamental thermodynamic nature of surface tension, that is, the free energy it would cost to form an interface. On the other hand, a localized surface fractal stress approach to surface tension via the calculation of localized fractal dimensions of the interface stress distribution would improve on both of these issues. First, in a system where thermo- and soluto-Marangoni effects are in play, the localized fractal stress would take into account both of these effects and their interactions with each other without any loss of generality. Second, the localized fractal stress can provide a standardized platform upon which we examine all possible effects on surface tension (whether local or global) in an uniform and consistent way; in this sense, a higher localized fractal stress can be directly interpreted as a higher energy cost required to form the interface in that region, thereby resulting in a higher localized surface tension. A complete hydrodynamic description of this interesting problem is not pursued in this study, but we anticipate numerous applications of this approach to the modeling of nanofluidic interfacial phenomena with a highly variable localized surface tension which would be realized in a future contribution.

Surface Percolating Clusters

The accurate identification of the intrinsic surface layer allows us to further the analysis of the atomic interaction network by applying the concept of percolation.[63−68] This method examines the “connectedness” of the different sites (or cells) on the interface that experience stress lower than a threshold. For a critical stress, a connected network of sites or a spanning cluster is formed that spans from left to right and from top to bottom of the lattice. If such spanning clusters form in at least 50% of the configurations, a system is described as a percolating system[69,70] and the corresponding stress is known as the percolating threshold.[66,71−74] Here, by configuration, we mean the instantaneous behavior of the system captured at a single time step. We determine the stress percolating threshold at the intrinsic interface and assess how temperature effects this threshold.

To quantify the threshold, we consider percolation in both directions (left to right and top to bottom) and apply next-nearest-neighbors algorithm. Panels a and b of Figure show two randomly selected instances, respectively, for T = 0.8 and 1.0, where spanning clusters form for the case of T = 0.8 but no such cluster is seen for T = 1.0. Such networks can be interpreted as a manifestation of the stress heterogeneity at the surface, which essentially increases with temperature resulting in a lower surface tension. Indeed, at lower temperature, the majority of the surface experiences higher negative stress, and thereby a spanning network can form at a relatively larger (negative) stress threshold. On the contrary, as temperature increases, the heterogeneity in stress distribution increases, too, and hence the formation of a spanning network requires the inclusion of smaller clusters. Figure c illustrates this, where the surface tension is plotted as a (linear) function of the percolating threshold. A lower value of the surface tension (which corresponds to a higher temperature) is seen to be associated with a lower (negative) percolating threshold, whereas a higher surface tension (lower temperature) is related to a larger magnitude of the threshold. A similar functional dependency between the percolating threshold and the fractal dimension of the stress networks can be seen in Figure d, which further echoes the higher space-filling nature of the surface clusters with lower stress heterogeneity at a lower temperature.

Figure 8

Percolation network, with periodic images shown faintly, at a random instant at temperature (a) T = 0.8 and (b) T = 1.0. The largest cluster is colored red. Panel a shows a spanning cluster, whereas in panel b no cluster is spanning throughout the surface for the same stress threshold. (c) Surface tension and (d) fractal dimension as functions of the percolation threshold over the range of temperatures considered in this study. Lines are linear fits.

Not only the stress networks at distinct temperatures differ at the percolating threshold, but their abilities for cluster formation also vary at any stress. In Figure , the average numbers of clusters are plotted for a range of stresses. It is evident that, for any particular stress, the (average) number of clusters at a lower temperature is less than that at a higher temperature, meaning that the low temperature clusters are larger and less heterogeneous, in terms of stress, than their high temperature counterparts. For instance, the top panel in the right-hand side of Figure shows three instantaneous networks for (negative) stresses with magnitudes greater than 8, 6, and 2, with the largest clusters colored in red. It is seen that the clusters grow in size as the field becomes more inclusive of stresses from panel I to panel III. Similar growth of the cluster size is seen for a lower temperature surface, i.e., for T = 0.75 as in the bottom panel. However, what differs between the networks at the two temperatures of interest here is that, for an identical range of stresses, it is more probable to obtain a larger cluster for a lower temperature surface; see Figure S3 and related discussion.

Figure 9

Average number of clusters at different stresses for T = 0.75 (red filled triangles), 0.9 (blue filled squares), and 1.0 (black filled circles). The top panel in the right-hand side corresponds to three instantaneous networks formed with stresses I (−∞, −8], II (−∞, −6], and III (−∞, −2] for T = 1.0. The bottom panel shows similar networks but for T = 0.75. The average numbers of clusters associated with I, II, and III are circled in the left panel.

From a thermodynamic point of view, a liquid–vapor interface at a higher temperature is (thermally) more energetic[75] and, thereby, favors a broader interface and weaker spatial correlation between atoms with larger nearest-neighbor distances causing a lower surface tension.[76] This, however, is challenging to perceive via a mechanical route, with the inherent difficulty lying in separating the actual surface layer from the capillary fluctuations.[64] Such a complication is circumvented in this study by detecting the interface using the intrinsic sampling method and collecting stresses in a reference frame moving with the surface which establishes the ground for an authentic surface-network analysis. Through this route, a thermodynamically weaker surface (due to the weaker spatial correlation between atoms) is mechanically represented by the lower “connectivity”[77,78] or higher heterogeneity of the stress networks. Stated differently, a high temperature surface (of which the surface tension is lower) can be mechanically described as a surface that is loosely interconnected by disjointed networks of stresses. Although the sensitivity of the analysis is limited by the probabilistic nature of the quantities investigated, the sole analysis of the sharpest atomic description of the intrinsic interface proves to be useful to unravel the variation of surface properties with temperature.

Whereas the intriguing notion of interpreting surface properties through surface coverage remains challenging,[79,80] the stress network approach presented in this paper is quite straightforward. Simplification of the interface by modeling an LJ fluid is a key limitation of this study, but the approaches developed here could be extended to molecular fluids in future work.

Conclusion

From the results of MD simulations of the liquid–vapor interface of an LJ fluid, we present a mechanical interpretation of surface tension and its variation with temperature. Intrinsic sampling method is used to define a moving frame of reference, and an equation for the interfacial density, ρ = ρ(x, T), for T ∈ [0.7, 1.0] and x ∈ (−∞, −1] is presented. We have identified stress clusters at the intrinsic surface and analyzed their non-uniform spatial correlations. The atomic interactions at the intrinsic surface layer can be thought of as a network of stresses holding the surface altogether. At an elevated temperature, the network becomes more disjointed and thus less stable.

Both the fractal dimension and the percolating threshold of the stress network can be correlated with the surface tension. These observations suggest that the pattern formation and the connectivity of the stress networks at the intrinsic interface are good indicators of surface tension. Importantly, the analysis of the stress network only at the outermost atomic layer suffices to provide a consistent prediction over a range of temperatures. The surface tension acquired from the simulation agrees well with previous MD simulations using different methods and experimental data for liquid argon, which advocates for the extrapolation of the findings of this study to liquid–vapor interfaces of molecular systems of interest.