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 identify 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 bio-locomotion [Hu et al., 2003, Bush and Hu, 2006, Houghton et al., 2018], plastron respiration [Flynn and Bush, 2008], and in numerous technological processes ranging from efficient oil recovery [Nobakht et al., 2007] to organic electronics [Forrest and Thompson, 2007] to capillarity driven thermal management [Antao et al., 2016, Rath and Flynn, 2020].
To characterize the liquid-vapor interface, a number of approaches have been developed over the past century, from macro-scale treatment [De Gennes et al., 2013] to the sophisticated molecular dynamics (MD) simulations [Ghoufi et al., 2016], which have changed our perception of the interface quite dramatically [Tarazona et al., 2012].

Among various available MD techniques, the mechanical route to the surface tension measurement has gained popularity for its relative simplicity and accuracy. This method replaces the scalar pressure with a second-rank pressure tensor [Walton et al., 1983] 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] ([Irving and Kirkwood, 1950]) definition, the diagonal elements of the pressure tensor can be expressed as:

\begin{equation}\label{Eq:IKpressure}
\boldsymbol{P}(\boldsymbol{x},t) = \sum_{i=1}^{N} \frac{\boldsymbol {p}_i \boldsymbol{p}_i}{m_i} \delta(\boldsymbol{x}_i - \boldsymbol{x}) + \frac{1}{2} \sum_{i,j}^{N} \boldsymbol{f}_{ij}\boldsymbol{x}_{ij} O_{ij}\delta(\boldsymbol{x}_i -\boldsymbol{x})
\end{equation}

where, the momentum $\boldsymbol{p_i} = m_i (\boldsymbol{\dot{x}}_i -u)$, ${\boldsymbol{\dot{x}}}$ is the total particle velocity and $u$ is its streaming part, $\delta(\dots)$ is the Dirac delta function, $O_{ij}$ is the Irving Kirkwood operator [Todd et al., 1995], $\boldsymbol{f}_{ij}$ is the force exerted by atom $j$ on atom $i$, and $\boldsymbol{x}_{ij} = \boldsymbol{x}_i -\boldsymbol{x}_j$ 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 ($P_N$) and the tangential ($P_T$) components of the pressure, surface tension can be obtained from the Hulshof integral, [Hulshof, 1901] as in Eq. (<ref>):

\begin{equation}\label{Eq:surfTenMec}
\gamma  = \int_{-\infty}^{\infty} \left[ P_N (x) - P_T (x) \right] dx .
\end{equation}

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 neighbourhood of the interface and only at a location far from the interface, this becomes uniform and equal to the normal component. Although Eq. (<ref>) 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 away from the interface. We will further elucidate this in the sec:resultsndiscuss section of this manuscript.

Attaining a molecular/atomic description of the interface, which is essential to understand the microscopic origin of surface tension [Antonin et al., 2011], is challenging due to the inhomogeneous nature of the interface and the difficulties 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 [Malijevskỳ and Jackson, 2012].
In the past years, numerous efforts have been made, both experimentally [Kuzmin and Romanov, 1994, Mitrinović et al., 2000] and by computer simulations [Senapati and Berkowitz, 2001, Willard and Chandler, 2010, Höfling and Dietrich, 2015, Stephan et al., 2018] to advance our understanding of the interface. A comprehensive review can be found in Ghoufi et al., 2016 and Ghoufi and Malfreyt, 2019.
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 behaviours [Bresme et al., 2008, Braga et al., 2018, Willard and Chandler, 2010].
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 [Bresme et al., 2008].
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., 1965 In their capillary wave theory, [Buff et al., 1965] formalised the concept of the existence of an interface $\zeta (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., 1999 exploited this idea of decoupling the capillary-wave broadening from the intrinsic interface and showed that the surface tension measured from Eq. (<ref>) 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 [Chacón et al., 2006].
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 between water-oil surface [Bresme et al., 2008] and to define the thickness of an adsorbed liquid film on a substrate [Fernández et al., 2011].

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., thermo-capillary flows. [Jasnow and Vinals, 1996]
In the present study, we carefully examine the temperature effects on the intrinsic density and pressure profiles. Afterwards,
we apply fractal analysis and the concept of percolating networks to analyse the temperature dependency of the spatial correlation of the interatomic interactions at the surface layer.

§ METHODOLOGY

MD simulations with the Flowmol MD code [Smith, 2013] were used to model the liquid-vapor interface of a Lennard-Jones (LJ) [Jones, 1924] fluid at different temperatures.
The initial simulation domain was a cubic box of dimensions $L_x = 120.64$, $L_y =19.05$ and $L_z = 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-centred cubic (FCC) lattice with a starting density of $\rho=1$, then molecules were randomly deleted until the pre-set liquid ($\rho=0.5$) and vapor ($\rho=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 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 Nośe-Hoover thermostat, [Nosé, 1984, Hoover, 1985] for 50,000 time steps with $\Delta T = 0.005$. A shifted LJ potential was used with periodic boundary conditions applied in all three Cartesian directions and the Verlet Leapfrog [Verlet, 1967] integrator was used to integrate the equations of motion.

For surface tension calculation, a small cutoff radius, $r_c$ results in significant deviation from experimental data for argon [Nijmeijer et al., 1988], whereas, $r_c\geq 4$ is found to give better agreement [Smith et al., 2016]. Surface tension being one of our prime interests in this study, and as recommended by Ghoufi et al., 2016, $r_c = 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. [Ghoufi et al., 2016] 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 cut-off length, [Stillinger Jr, 1963] $r_d = 1.5$, and with the required criterion of each atom having more than three neighbours. [Smith and Braga, 2020]
Instead of assuming an average interaction contour, the functional form of the interface, $\xi (y,z,t)$ was refitted at every time step using the intrinsic sampling method (ISM), [Chacón and Tarazona, 2003]  which approximates the liquid-vapor interface by means of a Fourier series representation as in Eq. (<ref>).

\begin{equation}\label{eq:fourier}
\xi (y,z,t) = \sum_{\textbf{k}<k_u} \hat{\xi_\textbf{k}} (t) ~\mbox{exp}(2\pi i \textbf{k} \cdot \boldsymbol{r}_\parallel)
\end{equation}

Here, $\boldsymbol{k}$ is the wave vector that corresponds to the periodic boundary conditions, i.e., $k = 2 \pi (n_y/L_y, n_z/L_z)$ with $n_y, n_z = 0, \pm 1, \pm 2, ... ,~k_u$ is the modulus of the wave vector, and $\boldsymbol{r}_\parallel$ denotes the parallel surface components in $y$ and $z$ directions.
The fitted coefficients, $\hat{\xi_\textbf{k}}$ are expressed as  time dependent functions because they are refitted to the surface each time the position of the surface atoms change (i.e. every time step).
Fig. <ref> (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 Fig. <ref> (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.

[Color online] (a) Liquid-vapor coexistence (b) fitted liquid-vapor surface function $\xi$ ($yz$ plane) (c) mapped grids at the interface, $x=0$ (mapped coordinates are in red at the top and un-mapped 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, $\xi$ 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, 1950 definitions can be integrated over a volume where the surfaces in the $x$ direction follow the function given by Eq. (<ref>), with uniform grids in $y$ and $z$ directions.
The density in a volume moving with the interface is then,

\begin{equation}
\int_V \rho\, dV = \displaystyle\sum_{i=1}^N m_i \vartheta_i
\end{equation}

where $\vartheta_i$ is the product of Heaviside functions which is unity when an atom is inside a volume and zero if outside [Smith and Braga, 2020].
This can be expressed in terms of this product, $\vartheta_i = \Lambda_x(x_i) \Lambda_y(y_i) \Lambda_z(z_i) $, where, the difference between two Heaviside functions is known as a boxcar function and represented by $\Lambda$. 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,

\begin{align}
\Lambda_a (x) = 1_{a^- < x < a^+} (x) = \left\{
\begin{array}{ll}
1 & \quad a^- < x < a^+ \\
0 & \quad \mathrm{otherwise},
\end{array}
\right.
\label{Indicator_fn_boxcar}
\end{align}

which checks whether $x$ is between the two limits, $a^-$ and $a^+$.
In the $y$ and $z$ direction, these indicator functions $\Lambda_y(y_i)$ and $\Lambda_z(z_i)$ check whether the atom positions $\{y_i, z_i\}$ 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 centred on $x$ with size $\Delta x$ the limits are $x^{\pm} = x \pm \Delta x/2 + \xi(y_i, z_i)$.
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_i^{'} = x_i - \xi(y_i, z_i)$ and then simply bin as if using a uniform grid, i.e., ${\tt{bin} = round}(x^{'}/ \Delta x)$.
The grid is uniform in the other two directions so the binning in $y$ and $z$ is unchanged, see fig. <ref> (c).
Quantities such as momentum, temperature and pressure inside a volume [Cormier et al., 2001] 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 [Smith and Braga, 2020].
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 [Smith, 2021].
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,

\begin{equation}
\boldsymbol{P}_{\alpha}^\text{Total} = \boldsymbol{P}_{\alpha}^\text{kin.} + \boldsymbol{P}_{\alpha}^\text{SM.} + \boldsymbol{P}_{\alpha}^\text{config.}
\end{equation}

where $\alpha \in \{x,y,z\}$ denotes the three directions where each is a vector $\boldsymbol{P}_{\alpha}=[P_{x \alpha}, P_{y \alpha}, P_{z \alpha} ]^T$, with the kinetic contribution, $\boldsymbol{P}^\text{kin.}_{\alpha}$, which comes from the momentum transport of atoms crossing the interface [Berry, 1971], the configurational contribution, $\boldsymbol{P}^\text{config.}_{\alpha}$ that arises from the atomic interaction [Walton et al., 1983],  and $\boldsymbol{P}^\text{SM.}_{\alpha}$ $-$ 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, $\boldsymbol{n}_y=\boldsymbol{e}_y$ and $\boldsymbol{n}_z=\boldsymbol{e}_z$, while for the intrinsic surface, it is given by $ {\boldsymbol{n}_{x}}=\boldsymbol{\nabla}_{s}\left(\xi-x_{s}\right)/||\boldsymbol{\nabla}_{s}\left(\xi-x_{s}\right)||$ where the subscript, $s$, denotes the derivative taken at the point of crossing $x_s$.
By assuming the convective term to be zero, $\rho \boldsymbol{u} \boldsymbol{u} = 0$, the three pressure contributions discussed above can be obtained in a molecular simulation as follows,

\begin{align}
\int_{t_1}^{t_2}  \boldsymbol{P}^\text{kin.}_{\alpha}  dt
= & \frac{1}{\Delta S_{\alpha}}  \displaystyle\sum_{i=1}^N  m_i \boldsymbol{\dot{r}}_i \boldsymbol{r}_{i_{12}} \cdot
\frac{ {\boldsymbol{n}_{\alpha}} }{|\boldsymbol{r}_{12} \cdot {\boldsymbol{n}}_{\alpha} |} dS_{\alpha}
\nonumber  \\
\int_{t_1}^{t_2}  \boldsymbol{P}^\text{SM.}_{\alpha}  dt
= &  \frac{1}{\Delta S_{\alpha}} \displaystyle\sum_{i=1}^N  m_i \boldsymbol{\dot{r}}_i  \vartheta_t
\nonumber  \\
\boldsymbol{P}_{\alpha}^\text{config.}= & \frac{1}{2 \Delta S_{\alpha}}\sum_{i,j}^{N} \boldsymbol{f}_{ij}  \boldsymbol{{r}}_{ij} \cdot
\frac{ {\boldsymbol{n}_{\alpha}} }{|\boldsymbol{r}_{ij} \cdot {\boldsymbol{n}_{\alpha}}|}dS_{\alpha}
\label{P_equation}
\end{align}

where $\Delta S_{\alpha}$ is the surface area, $\boldsymbol{r_{i_{12}}} = \boldsymbol{r}_i(t_2) - \boldsymbol{r}_i(t_1)$ denotes the vector for the movement of an atom from time $t_1$ to $t_2$ while $\boldsymbol{r}_{ij}$ is the separation vector.
The $\vartheta_t$ 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 $\vartheta_t = \Lambda_{\xi} (x_i)\, \Lambda_{y} (y_i)\, \Lambda_{z} (z_i)$, the surface evolution from the start of a timestep $\xi^- = \xi(t)$ to the end $\xi^+=\xi(t+\Delta t)$ is multiplied by an indicator function, i.e.,  Eq. (<ref>).

The $dS_{\alpha}$ term in Eq.(<ref>) is the derivative of the $\vartheta$ function, with respect to $\alpha=\{x,y,z\}$, and is only non-zero if the separation vectors $\boldsymbol{r}_{i_{12}}$ or $\boldsymbol{r}_{ij}$ are crossing the surface of the volume in question.  Without loss of generality, we consider the expression for the surface and the separation vector, $\boldsymbol{r}_{ij}$ in the $x$ direction as

dS_x = ∑_k=1^N_roots
Λ_λ (λ_k)

where $\lambda$ parameterises the line between atoms $\boldsymbol{r}_{\lambda} = \boldsymbol{r}_i + \lambda \boldsymbol{r}_{ij}$ with $\lambda_k$ the value at a point of crossing on a surface.
The function $\Lambda_{\lambda}$ therefore checks if a $\lambda_k$ value is between $\boldsymbol{r}_i$ ($\lambda^- = 0$) and $\boldsymbol{r}_j$ ($\lambda^+ = 1$), with the remaining functions checking if the point of crossing $\{y_{\lambda}(\lambda_k), z_{\lambda}(\lambda_k) \}$ is on the volume surface between $y^-$ and $y^+$ and $z^-$ and $z^+$.
The form of $dS_y$ and $dS_z$ are similar and atomic motions are parameterised in the same way with $\boldsymbol{r}_{\lambda} = \boldsymbol{r}_1 + \lambda \boldsymbol{r}_{i_{12}}$.

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. (<ref>) is converted to a set of bilinear patches of the form

ξ^BL (y,z) = a_o + a_y y + a_z z + a_yz yz

The 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. (<ref>).
The procedure to fit the intrinsic interface and to choose the number of bilinear patches, as well as to calculate the intersect are described in previous work  [Smith and Braga, 2020, Smith, 2021].

§ RESULTS AND DISCUSSIONS

§.§ Temperature effects on density profiles

[Color online] (a) Density profile (blue circles) and derivative of density (black dashed line) for $T=0.70$ shows oscillations near the liquid-vapor interface (b) density profiles for a range of temperatures, oscillations are damped as temperature increases.

The intrinsic density profile at temperature $T=0.7$ is shown in Fig. <ref> (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 behaviour of the free surface [Chacón et al., 2001], 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.

As seen from Fig. <ref> (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\rho/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., 1999][Sides et al., 1999] for further details.
Fig. <ref> (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 dampens quite abruptly in the vapor side right after the interface ($x>0$), despite the existence of a small peak due to the adsorbed layer.

[Color online] (a) Oscillation of the intrinsic density profiles:  circles denote results from MD simulations, solid lines are exponential fitting, as in Eq. (<ref>), to the density profiles at different temperatures, (b) the damping coefficient (left $y$ axis,  ) increases with temperature, whereas, the frequency of oscillation (right $y$ axis, ) decreases. The solid lines are linear fits to the data.

The oscillations of the density profile at the liquid side for $x \in (-\infty, -1]$ can be approximated, as illustrated in Fig. <ref> (a), by an exponential decay function of the form:

\begin{equation}\label{eq:dampen}
\rho (x) = \rho_\mathrm{bulk} + e^ {bx} \cos (\mathrm{f}\, x)
\end{equation}
where, $\rho (x)$ denotes the local intrinsic density, $\rho_\mathrm{bulk}$ is the liquid density of the bulk phase far away from the interface, $b$ is the damping coefficient, and $\mathrm{f}$ is the frequency of oscillation.
Note that the negative sign of the exponent, i.e., the damping coefficient in Eq. (<ref>) 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 Fig. <ref> (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.51\,T - 0.81$ with $R^2 = 0.95$ and $\mathrm{f} = -0.7\pi\,T + 2.5\pi$ with $R^2 = 0.99$.
Thus the intrinsic density in Eq. (<ref>) can be expressed as a function of location, bulk density of the liquid and the temperature, i.e.,

ρ(x) = ρ_bulk + e^(2.51T-0.81)x [π(2.5-0.7T) x]  for  x ∈(-∞, -1]

Whilst 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 [Evans et al., 1993].

§.§ Temperature effects on pressure profiles

[Color online] Pressure profiles at and near the surface layer for T= 0.70. (a) Components of the normal pressure: () configurational part, $P^\mathrm{config.}$,   () kinetic part, $P^\mathrm{kin.}$, () contributions from surface movement, $P^\mathrm{SM.}$, and () the total normal pressure, $P^\mathrm{Total}$. (b) Normal pressure without () and with () the consideration of surface movement term, $l_s$ denote the effective fluctuating length, over which the surface movements are visible. (c) Linear increase of the effective fluctuating length, $l_s$ with temperature, the solid line is a linear fit to the data with a slope of $2.04$, an intercept of $0.23$, and $R^2= 0.97$)

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.
Fig. <ref> (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, 2020, and Smith, 2021 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, $P^\mathrm{SM.}$.
This, when considered along with the configurational and the kinetic contributions of the normal pressure, as shown in Fig. <ref> (a), exactly balances the momentum change to machine precision.

The essence of the surface movement contribution is illustrated in  Fig. <ref> (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 Fig. <ref> (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, $l_s$ (schematically shown in Fig. <ref> b).
More formally, we define $l_{s}$ as

\begin{equation}\label{eq:eff_fluc_len}
\end{equation}

which is the maximal distance between any two roots of the equation $\left|P(x)\right|=\epsilon$, for $0<\epsilon\ll1$ is the amplitude of the very small, but finite oscillations either side of the effective fluctuating length. In this study, we use a cut-off value $\epsilon=0.05$ in order to numerically determine $l_{s}$ from the fluctuations of $P(x)$.
As shown in Fig. <ref> (c), $l_s$ increases linearly with the temperature.
Despite that the interfacial thickness assumes approximately the size of an atom at a temperature away from the critical temperature [Lu and Hu, 2008], 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 [Goujon et al., 2016].

Although not straight forward, a plausible link between the effective length, $l_s$ and the shear viscosity can be inferred through the hydrodynamic description of the capillary waves. It is well established [Levich, 1962, Delgado-Buscalioni et al., 2008, Shen et al., 2018] that the damping rate, $\Gamma$ of an over-damped capillary wave mode satisfies the asymptotic relation, $\Gamma \sim q \gamma / \mu$ where $q$ is the wavenumber and $\mu$ is the dynamic shear viscosity. At critical damping, the real part of the complex wave frequency is zero, i.e., $\mathrm{Re} (\omega) =0$ and $q \sim \gamma \rho / \mu^2$. If the propagation velocity of the density fluctuations, $v_0$, is further assumed to be constant, it can then be postulated that $l_s$ scales inversely with $\Gamma$, i.e, $\Gamma \sim {v_0/l_s} \sim \gamma^2 \rho / \mu^3$. Crucially, this rearranges to give

\begin{equation}
l_s \sim \frac{v_0 \mu^3}{\gamma^2 \rho},\label{eq:l_s-mu}
\end{equation}

thereby provides 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. However, a detailed analysis to confirm this postulate is beyond the scope of this manuscript but this relationship between $l_s$ and $\mu$ certainly prompts further investigation.

[Color online] (a-d) Constituent parts of the () normal component of pressure, i.e.,  () the kinetic part and, () 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 () the tangential component of the pressure tensor and of () $P_N - P_T$  near the interface for $T = 0.70$ to $1.0$.

The temperature effect on the intrinsic pressure profile has further been examined in Fig. <ref>.
Though the shape of the normal pressure profile remains unaltered irrespective of temperature,  the oscillations of the configurational and the kinetic parts, in Fig. <ref> (a-d), are seen to smooth out as temperature increases.
The normal and the tangential components of the pressure, along with $(P_N-P_T)$, are plotted in the Fig. <ref> (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. (<ref>) is a particle property, and thereby proportional to the intrinsic density through $P^\mathrm{kin.} = \rho k_B T$, where $k_B$ is the Boltzman constant. As such these oscillations are directly comparable to those in the intrinsic density profiles, as in Figs. <ref> and <ref>. The equal and opposite of this kinetic component is the configurational pressure, as shown in Fig. 5 (a-d), which ensures equilibrium.
The observations made in Fig. <ref> (c) and Fig. <ref> allow to conclude that an increase in temperature dampens the oscillations of the intrinsic profiles, and at the same time, broadens the effective fluctuating length, $l_s$, over which surface movement effects are present.

[Color online] (a) Difference between the normal and the tangential pressure,  $P_N - P_T$ () near the interface. The cumulative integral of $P_N - P_T$, i.e., the surface tension, $\gamma(x) = \int_{-10}^x \left[P_N (x^\prime) - P_T (x^\prime)\right] dx^\prime$ is shown by the 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 modelling () 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 $R^2 = 0.99$) to the surface tension obtained from the present study.

The influence of the oscillations discussed, are confined to a small interfacial region.
Far from the liquid-vapor interface, both the normal 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 (both in the liquid and the vapor sides) contributes to the surface tension.
The term inside the integral of Eq. (<ref>) is plotted (black solid line) in Fig. <ref> (a) along with the cumulative integral, i.e, the surface tension, $\gamma$.
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 Fig. <ref> (a), where the large negative peak of the pressure difference at $x = 0$ contributes significantly to the integral of Eq. (<ref>).
We compare, in Fig. <ref> (b), the surface tension evaluated at different temperatures with available experimental data for liquid argon [Lemmon, 1998],
results from molecular simulations with truncated and shifted potential with a cut-off radius of 2.5 [Vrabec et al., 2006, Goujon et al., 2014], and 4.0 [Smith et al., 2016], and truncated potential with long range corrections and cut-off radius of 3.0 [Werth et al., 2013]- as seen in the figure, the comparison suggests good agreement.

§.§ Surface fractals

The preceding sections discuss the temperature effects 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 towards the surface tension, we examine the distribution of these stresses and ignore the kinetic components, at the intrinsic surface layer ($x=0$); i.e.,

P^config. (x=0) =  1/2 [ P_yy^config. (x=0) + P_zz^config. (x=0) ] .

Fig. <ref>(a) shows an instantaneous configurational stress distribution map at the interface. The filled black squares correspond to the atomic locations which the intrinsic interface is fitted to, at that particular instant.
In Fig. <ref>(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 Fig. S1, S2, movie S1 and related discussions).
These `non-trivial' networks and their variation with temperature are analysed by means of fractal analysis.

[Color online] Instantaneous stress map at the liquid-vapor interface for T = 0.80. (a) Color map of the configurational stress distribution with atomic positions overlaid (filled black 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 standard error of the mean of fractal dimension, SEM, the solid line is a linear fit with $R^2 = 0.95$.

To quantify the temperature dependency of the structural complexity of the stress networks, we measure their fractal dimensions, $\Phi$. This dimension is often used to quantify the space filling nature, heterogeneity or self-similarity of surfaces, clusters etc. [Avnir et al., 1984, Jelinek et al., 2006, Chen et al., 2015, Shen et al., 2020]
$\Phi$ is calculated through the Minkowski dimension [Hassan et al., 2012, Kempkes et al., 2019], which is often referred to as the box-counting dimension, whereby, for non-overlapping $N$ boxes with sides $\epsilon$,

\begin{equation}\label{eq:boxcount}
\Phi = \lim_{\epsilon \to 0} \frac{\log N (\epsilon)}{\log (1/\epsilon)}.
\end{equation}

The algorithm employed to calculate $\Phi$ consists of converting the network maps (as in Fig. <ref> b) into binary images such that the stress networks are depicted by black pixels on white background. For a given grid side, $\epsilon$ the number of grids, $N(\epsilon)$ 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, $\Phi$ is then obtained from the slope of $\mbox{log}(N)$ vs $\mbox{log}(\epsilon)$. See  supporting information for further details.

With the variation of temperature, we have observed consistent behaviour of the fractal dimension of the network comprised of the negative configurational stresses that lie within $\left[ P^\mathrm{config.}_\mathrm{n,m} - \sigma,  P^\mathrm{config.}_\mathrm{n,m}\right]$, where $P^\mathrm{config.}_\mathrm{n,m}$ denotes the mean of all the negative configurational stresses, and $\sigma$ denotes the standard deviation.
The fractal dimensions thus obtained are compared against the corresponding surface tensions in
Fig. <ref> (c).
For the temperature range considered here, $\Phi$ is seen to linearly correspond to the surface tension as $\gamma = 9.86\, \Phi\ - 13.86$ which is shown by the solid line.
This reflects the greater space filling nature [Jelinek et al., 2006] 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, $\gamma$ is often modelled 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 localised surface fractal stress approach to surface tension via the calculation of localised fractal dimensions of the interface stress distribution would improve on both of these issues.
Firstly, in a system where thermo- and soluto-Marangoni effects are in play, the localised fractal stress would take into account both of these two effects and their interactions with each other without any loss of generality. Secondly, the localised fractal stress can provide a standardised 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 localised fractal stress can be directly interpreted as a higher energy cost required to form the interface in that region, thereby resulting in a higher localised surface tension.
A complete hydrodynamical description of this interesting problem is not pursued in this study, but we anticipate numerous applications of this approach to the modelling of nano-fluidic interfacial phenomena with a highly variable localised surface tension which would be realised in a future contribution.

[Color online] 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.

§.§ 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 [Pathak et al., 2017, Sega et al., 2014, Stauffer and Aharony, 2018, Heyes and Melrose, 1988, de Souza and Harrowell, 2009, Dapp et al., 2012]. 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$ percent of the configurations, a system is described as a percolating system [Zarragoicoechea et al., 2005, Seaton and Glandt, 1987] and the corresponding stress is known as the percolating threshold [Bug et al., 1985, Heyes and Melrose, 1988, Heyes and Melrose, 1989, Heyes and Melrose, 1990, Hasmy et al., 2021].
Here, by configuration, we mean the instantaneous behaviour 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 neighbours algorithm. Fig. <ref> (a) and (b) 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.
Fig. <ref> (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.
Similar functional dependency between the percolating threshold and the fractal dimension of the stress networks can be seen in Fig. <ref> (d), which further echos the higher space filling nature of the surface clusters with lower stress heterogeneity at a lower temperature.

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

Not only the stress networks at distinct temperatures differ at the percolating threshold, but their ability of cluster formation also varies at any stress.
In Fig. <ref>, the average number 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 Fig. <ref> 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 $\mathrm{I}$ to $\mathrm{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: for identical range of stresses, it is more probable to obtain a larger cluster for a lower temperature surface, see Fig. S3 and related discussions.

From a thermodynamic point of view, a liquid-vapor interface at a higher temperature is (thermally) more energetic[Palmer, 1976] and thereby, favours a broader interface and weaker spatial correlation between atoms with larger nearest-neighbours distances causing a lower surface tension [Bhanvase and Barai, 2021].
This, however, is challenging to perceive via a mechanical route - the inherent difficulty lying in separating the actual surface layer from the capillary fluctuations [Sega et al., 2014].
Such 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' [Burnley, 2013, Alvarado et al., 2017] or higher heterogeneity of the stress networks. Stated differently, a high temperature surface (of which, surface tension is lower) can be mechanically described as a surface that is loosely inter-connected 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 [Mu et al., 2006, Menger and Rizvi, 2011], the stress network approach presented in this paper is quite straight forward. 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, $\rho = \rho(x,T)$ for $T\in [0.7,1.0]$ and $x\in(-\infty,-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.

§ CONFLICT OF INTEREST

The authors declare no conflict of interest.

§.§ Supporting Information

– Fractal Analysis, Percolation Network, Liquid-Vapor coexistence data (PDF).

– Temporal variation of fractal networks (Movie).

§.§ Acknowledgement

The authors wish to thank Prof. D.M. Heyes (Department of Mechanical Engineering, Imperial College London, UK) for insightful discussions.
M.R.R. acknowledges PhD studentship funding from Shell via the University Technology Centre for Fuels and Lubricants and the Beit Trust for the Beit Fellowship for Scientific Research. L.S. thanks the Engineering and Physical Sciences Research Council (EPSRC) for a Postdoctoral Fellowship (EP/V005073/1). J.P.E. was supported by the Royal Academy of Engineering through a Research Fellowship. D.D. thanks the EPSRC for an Established Career Fellowship (EP/N025954/1).

[Hu et al., 2003]
Hu, D. L.; Chan, B.; Bush, J. W. The Hydrodynamics of Water Strider Locomotion.
Nature 2003, 424, 663–666
[Bush and Hu, 2006]
Bush, J. W.; Hu, D. L. Walking on Water: Biolocomotion at the Interface.
Annual Review of Fluid Mechanics 2006, 38,
[Houghton et al., 2018]
Houghton, I. A.; Koseff, J. R.; Monismith, S. G.; Dabiri, J. O. Vertically
Migrating Swimmers Generate Aggregation-Scale Eddies in a Stratified Column.
Nature 2018, 556, 497–500
[Flynn and Bush, 2008]
Flynn, M. R.; Bush, J. W. Underwater Breathing: the Mechanics of Plastron
Respiration. Journal of Fluid Mechanics 2008, 608,
[Nobakht et al., 2007]
Nobakht, M.; Moghadam, S.; Gu, Y. Effects of Viscous and Capillary Forces on
C$O_2$ Enhanced Oil Recovery under Reservoir Conditions. Energy &
Fuels 2007, 21, 3469–3476
[Forrest and Thompson, 2007]
Forrest, S. R.; Thompson, M. E. Introduction: Organic Electronics and
Optoelectronics. Chemical Reviews 2007, 107,
[Antao et al., 2016]
Antao, D. S.; Adera, S.; Zhu, Y.; Farias, E.; Raj, R.; Wang, E. N. Dynamic
Evolution of the Evaporating Liquid–Vapor Interface in Micropillar Arrays.
Langmuir 2016, 32, 519–526
[Rath and Flynn, 2020]
Rath, A.; Flynn, M. R. Core Annular Flow Theory as Applied to the Adiabatic
Section of Heat Pipes. Physics of Fluids 2020, 32,
[De Gennes et al., 2013]
De Gennes, P.-G.; Brochard-Wyart, F.; Quéré, D. Capillarity and
Wetting Phenomena: Drops, Bubbles, Pearls, Waves; Springer Science &
Business Media, 2013
[Ghoufi et al., 2016]
Ghoufi, A.; Malfreyt, P.; Tildesley, D. J. Computer Modelling of the Surface
Tension of the Gas–Liquid and Liquid–Liquid Interface. Chemical
Society Reviews 2016, 45, 1387–1409
[Tarazona et al., 2012]
Tarazona, P.; Chacón, E.; Bresme, F. Intrinsic Profiles and the Structure
of Liquid Surfaces. Journal of Physics: Condensed Matter
2012, 24, 284123
[Walton et al., 1983]
Walton, J.; Tildesley, D.; Rowlinson, J.; Henderson, J. The Pressure Tensor at
the Planar Surface of a Liquid. Molecular Physics 1983,
48, 1357–1368
[Irving and Kirkwood, 1950]
Irving, J.; Kirkwood, J. G. The Statistical Mechanical Theory of Transport
Processes. IV. The Equations of Hydrodynamics. The Journal of Chemical
Physics 1950, 18, 817–829
[Todd et al., 1995]
Todd, B.; Evans, D. J.; Daivis, P. J. Pressure Tensor for Inhomogeneous Fluids.
Physical Review E 1995, 52, 1627
[Hulshof, 1901]
Hulshof, H. The Direct Deduction of the Capillarity Constant $\gamma$ as a
Surface Tension. Annalen der Physik 1901, 4,
[Antonin et al., 2011]
Antonin, M.; Joost, W.; Snoeijer, J.; Andreotti, B. Why is Surface Tension a
Force Parallel to the Interface? American Journal of Physics
2011, 999–1008
[Malijevskỳ and Jackson, 2012]
Malijevskỳ, A.; Jackson, G. A Perspective on the Interfacial Properties of
Nanoscopic Liquid Drops. Journal of Physics: Condensed Matter
2012, 24, 464121
[Kuzmin and Romanov, 1994]
Kuzmin, V.; Romanov, V. Influence of the Surface Profile on the Roughness
Contribution to the Ellipticity Coefficient. Physical Review E
1994, 49, 2949
[Mitrinović et al., 2000]
Mitrinović, D. M.; Tikhonov, A. M.; Li, M.; Huang, Z.; Schlossman, M. L.
Noncapillary-Wave Structure at the Water-Alkane Interface. Physical
Review Letters 2000, 85, 582
[Senapati and Berkowitz, 2001]
Senapati, S.; Berkowitz, M. L. Computer Simulation Study of the Interface Width
of the Liquid/Liquid Interface. Physical Review Letters 2001,
87, 176101
[Willard and Chandler, 2010]
Willard, A. P.; Chandler, D. Instantaneous Liquid Interfaces. The Journal
of Physical Chemistry B 2010, 114, 1954–1958
[Höfling and Dietrich, 2015]
Höfling, F.; Dietrich, S. Enhanced Wavelength-Dependent Surface Tension of
Liquid-Vapour Interfaces. Europhysics Letters 2015,
109, 46002
[Stephan et al., 2018]
Stephan, S.; Liu, J.; Langenbach, K.; Chapman, W. G.; Hasse, H. Vapor-Liquid
Interface of the Lennard-Jones Truncated and Shifted Fluid: Comparison of
Molecular Simulation, Density Gradient Theory, and Density Functional Theory.
The Journal of Physical Chemistry C 2018, 122,
[Ghoufi and Malfreyt, 2019]
Ghoufi, A.; Malfreyt, P. Calculation of the Surface Tension of Water: 40 Years
of Molecular Simulations. Molecular Simulation 2019,
45, 295–303
[Bresme et al., 2008]
Bresme, F.; Chacón, E.; Tarazona, P. Molecular Dynamics Investigation of
the Intrinsic Structure of Water–Fluid Interfaces via the Intrinsic Sampling
Method. Physical Chemistry Chemical Physics 2008, 10,
[Braga et al., 2018]
Braga, C.; Smith, E. R.; Nold, A.; Sibley, D. N.; Kalliadasis, S. The Pressure
Tensor across a Liquid-Vapour Interface. The Journal of Chemical
Physics 2018, 149, 044705
[Buff et al., 1965]
Buff, F.; Lovett, R.; Stillinger Jr, F. Interfacial Density Profile for Fluids
in the Critical Region. Physical Review Letters 1965,
15, 621
[Sides et al., 1999]
Sides, S. W.; Grest, G. S.; Lacasse, M.-D. Capillary Waves at Liquid-Vapor
Interfaces: A Molecular Dynamics Simulation. Physical Review E
1999, 60, 6708
[Chacón et al., 2006]
Chacón, E.; Tarazona, P.; Alejandre, J. The Intrinsic Structure of the
Water Surface. The Journal of Chemical Physics 2006,
125, 014709
[Fernández et al., 2011]
Fernández, E. M.; Chacón, E.; Tarazona, P. Thickness and Fluctuations
of Free and Adsorbed Liquid Films. Physical Review B 2011,
84, 205435
[Jasnow and Vinals, 1996]
Jasnow, D.; Vinals, J. Coarse-Grained Description of Thermo-Capillary Flow.
Physics of Fluids 1996, 8, 660–669
[Smith, 2013]
Smith, E. On the Coupling of Molecular Dynamics to Continuum Computational
Fluid Dynamics, PhD thesis, Imperial College London. 2013
[Jones, 1924]
Jones, J. E. On the Determination of Molecular Fields.—II. From the Equation
of State of a Gas. Proceedings of the Royal Society of London. Series
A, Containing Papers of a Mathematical and Physical Character 1924,
106, 463–477
[34]
SAT-TMMC: Liquid-Vapor Coexistence Properties - Linear-Force Shifted
Potential at $4.0\,\sigma$; National Institute of Standards and Technology,
accessed: 2022-04-16
[Nosé, 1984]
Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical
Ensemble. Molecular Physics 1984, 52, 255–268
[Hoover, 1985]
Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions.
Physical Review A 1985, 31, 1695
[Verlet, 1967]
Verlet, L. Computer Experiments on Classical Fluids. I. Thermodynamical
Properties of Lennard-Jones Molecules. Physical Review 1967,
159, 98
[Nijmeijer et al., 1988]
Nijmeijer, M.; Bakker, A.; Bruin, C.; Sikkenk, J. A Molecular Dynamics
Simulation of the Lennard-Jones Liquid–Vapor Interface. The Journal of
Chemical Physics 1988, 89, 3789–3792
[Smith et al., 2016]
Smith, E.; Müller, E.; Craster, R.; Matar, O. A Langevin Model for
Fluctuating Contact Angle Behaviour Parametrised using Molecular Dynamics.
Soft Matter 2016, 12, 9604–9615
[Stillinger Jr, 1963]
Stillinger Jr, F. H. Rigorous Basis of the Frenkel-Band Theory of Association
Equilibrium. The Journal of Chemical Physics 1963, 38,
[Smith and Braga, 2020]
Smith, E. R.; Braga, C. Hydrodynamics Across a Fluctuating Interface. The
Journal of Chemical Physics 2020, 153, 134705
[Chacón and Tarazona, 2003]
Chacón, E.; Tarazona, P. Intrinsic Profiles beyond the Capillary Wave
Theory: A Monte Carlo Study. Physical Review Letters 2003,
91, 166103
[Cormier et al., 2001]
Cormier, J.; Rickman, J.; Delph, T. Stress Calculation in Atomistic Simulations
of Perfect and Imperfect Solids. Journal of Applied Physics
2001, 89, 99–104
[Smith, 2021]
Smith, E. R. The Importance of Reference Frame for Pressure at the
Liquid–Vapour Interface. Molecular Simulation 2021,
48, 57–72
[Berry, 1971]
Berry, M. The Molecular Mechanism of Surface Tension. Physics Education
1971, 6, 79
[Chacón et al., 2001]
Chacón, E.; Reinaldo-Falagán, M.; Velasco, E.; Tarazona, P. Layering at
Free Liquid Surfaces. Physical Review Letters 2001,
87, 166101
[Evans et al., 1993]
Evans, R.; Henderson, J.; Hoyle, D.; Parry, A.; Sabeur, Z. Asymptotic Decay of
Liquid Structure: Oscillatory Liquid-Vapour Density Profiles and the
Fisher-Widom Line. Molecular Physics 1993, 80,
[Lu and Hu, 2008]
Lu, X.; Hu, Y. Molecular Thermodynamics of Complex Systems; Springer,
2008; Vol. 131
[Goujon et al., 2016]
Goujon, F.; Ghoufi, A.; Malfreyt, P.; Tildesley, D. J. Can We Approach the
Gas–Liquid Critical Point Using Slab Simulations of Two Coexisting Phases?
The Journal of Chemical Physics 2016, 145,
[Levich, 1962]
Levich, V. G. Physicochemical Hydrodynamics; Prentice-Hall Inc.,
[Delgado-Buscalioni et al., 2008]
Delgado-Buscalioni, R.; Chacon, E.; Tarazona, P. Hydrodynamics of Nanoscopic
Capillary Waves. Physical Review Letters 2008, 101,
[Shen et al., 2018]
Shen, L.; Denner, F.; Morgan, N.; van Wachem, B.; Dini, D. Capillary Waves with
Surface Viscosity. Journal of Fluid Mechanics 2018,
847, 644–663
[Lemmon, 1998]
Lemmon, E. W. Thermophysical Properties of Fluid Systems; Nat. Inst.
Stand. Tech. Gaithersburg, 1998
[Vrabec et al., 2006]
Vrabec, J.; Kedia, G. K.; Fuchs, G.; Hasse, H. Comprehensive Study of the
Vapour–Liquid Coexistence of the Truncated and Shifted Lennard–Jones Fluid
Including Planar and Spherical Interface Properties. Molecular Physics
2006, 104, 1509–1527
[Goujon et al., 2014]
Goujon, F.; Malfreyt, P.; Tildesley, D. J. The Gas–Liquid Surface Tension of
Argon: A Reconciliation between Experiment and Simulation. The Journal
of Chemical Physics 2014, 140, 244710
[Werth et al., 2013]
Werth, S.; Lishchuk, S. V.; Horsch, M.; Hasse, H. The Influence of the Liquid
Slab Thickness on the Planar Vapor–Liquid Interfacial Tension. Physica
A 2013, 392, 2359–2367
[Avnir et al., 1984]
Avnir, D.; Farin, D.; Pfeifer, P. Molecular Fractal Surfaces. Nature
1984, 308, 261–263
[Jelinek et al., 2006]
Jelinek, H. F.; Jones, C. L.; Warfel, M. D.; Lucas, C.; Depardieu, C.;
Aurel, G. Understanding Fractal analysis? The Case of Fractal Linguistics.
Complexus 2006, 3, 66–73
[Chen et al., 2015]
Chen, D. Z.; Shi, C. Y.; An, Q.; Zeng, Q.; Mao, W. L.; Goddard III, W. A.;
Greer, J. R. Fractal Atomic-Level Percolation in Metallic Glasses.
Science 2015, 349, 1306–1310
[Shen et al., 2020]
Shen, L.; Denner, F.; Morgan, N.; Van Wachem, B.; Dini, D. Transient Structures
in Rupturing Thin Films: Marangoni-induced Symmetry-breaking Pattern
Formation in Viscous Fluids. Science Advances 2020, 6,
[Hassan et al., 2012]
Hassan, N.; Soltero, A.; Pozzo, D.; Messina, P. V.; Ruso, J. M. Bioinspired
Templates for the Synthesis of Silica Nanostructures. Soft Matter
2012, 8, 9553–9562
[Kempkes et al., 2019]
Kempkes, S. N.; Slot, M. R.; Freeney, S. E.; Zevenhuizen, S. J.;
Vanmaekelbergh, D.; Swart, I.; Smith, C. M. Design and Characterization of
Electrons in a Fractal Geometry. Nature Physics 2019,
15, 127–131
[Pathak et al., 2017]
Pathak, S. N.; Esposito, V.; Coniglio, A.; Ciamarra, M. P. Force Percolation
Transition of Jammed Granular Systems. Physical Review E
2017, 96, 042901
[Sega et al., 2014]
Sega, M.; Horvai, G.; Jedlovszky, P. Two-Dimensional Percolation at the Free
Water Surface and its Relation with the Surface Tension Anomaly of Water.
The Journal of Chemical Physics 2014, 141,
[Stauffer and Aharony, 2018]
Stauffer, D.; Aharony, A. Introduction to Percolation Theory; CRC press,
[Heyes and Melrose, 1988]
Heyes, D.; Melrose, J. Percolation Thresholds of Simple Fluids. Journal
of Physics A: Mathematical and General 1988, 21,
[de Souza and Harrowell, 2009]
de Souza, V. K.; Harrowell, P. Rigidity Percolation and the Spatial
Heterogeneity of Soft Modes in Disordered Materials. Proceedings of the
National Academy of Sciences 2009, 106, 15136–15141
[Dapp et al., 2012]
Dapp, W. B.; Lücke, A.; Persson, B. N.; Müser, M. H. Self-affine
Elastic Contacts: Percolation and Leakage. Physical Review Letters
2012, 108, 244301
[Zarragoicoechea et al., 2005]
Zarragoicoechea, G. J.; Pugnaloni, L. A.; Lado, F.; Lomba, E.; Vericat, F.
Percolation of Clusters with a Residence Time in the Bond Definition:
Integral Equation Theory. Physical Review E 2005, 71,
[Seaton and Glandt, 1987]
Seaton, N.; Glandt, E. D. Aggregation and Percolation in a System of Adhesive
Spheres. The Journal of Chemical Physics 1987, 86,
[Bug et al., 1985]
Bug, A.; Safran, S.; Grest, G. S.; Webman, I. Do Interactions Raise or Lower a
Percolation Threshold? Physical Review Letters 1985,
55, 1896
[Heyes and Melrose, 1989]
Heyes, D. M.; Melrose, J. R. Microscopic Simulation of Rheology: Molecular
Dynamics Computations and Percolation Theory. Molecular Simulation
1989, 2, 281–300
[Heyes and Melrose, 1990]
Heyes, D.; Melrose, J. Continuum Percolation of 2D and 3D Simple Fluids.
Molecular Simulation 1990, 5, 329–343
[Hasmy et al., 2021]
Hasmy, A.; Ispas, S.; Hehlen, B. Percolation Transitions in Compressed Si$O_2$
Glasses. Nature 2021, 599, 62–66
[Palmer, 1976]
Palmer, S. The Effect of Temperature on Surface Tension. Physics
Education 1976, 11, 119
[Bhanvase and Barai, 2021]
Bhanvase, B., Barai, D., Eds. Nanofluids for Heat and Mass Transfer;
Academic Press, 2021; pp 101–166
[Burnley, 2013]
Burnley, P. The Importance of Stress Percolation Patterns in Rocks and other
Polycrystalline Materials. Nature Communications 2013,
4, 1–6
[Alvarado et al., 2017]
Alvarado, J.; Sheinman, M.; Sharma, A.; MacKintosh, F. C.; Koenderink, G. H.
Force Percolation of Contractile Active Gels. Soft Matter
2017, 13, 5624–5644
[Mu et al., 2006]
Mu, Q.; Lu, J.; Ma, Y.; Paz de Banez, M.; Robinson, K.; Armes, S.; Lewis, A.;
Thomas, R. Neutron Reflection Study of a Water-Soluble Biocompatible Diblock
Copolymer Adsorbed at the Air/Water Interface: The Effects of pH and Polymer
Concentration. Langmuir 2006, 22, 6153–6160
[Menger and Rizvi, 2011]
Menger, F. M.; Rizvi, S. A. Relationship between Surface Tension and Surface
Coverage. Langmuir 2011, 27, 13975–13977