The pressure tensor across a liquid-vapour interface

Inhomogeneous systems are characterised by an interface, across which physical quantities vary sharply, over a length scale of a few molecular diameters.¹ Thermally induced fluctuations of the interface in the form of capillary waves smooth out the transition of properties across the interface. Although capillary wave behaviour has been observed experimentally,^2,3 analysis of the interfacial structure at the sharpest molecular resolution is difficult as a result of these fluctuations. Capillary wave theory (CWT), introduced by Buff et al.,^4,5 accounts for these by assuming the existence of an intrinsic surface $ξx$⁠, where x is a vector in the coordinate system parallel to the surface. This function acts as an instantaneous boundary between fluid phases, against which statistical averages are computed, e.g., the intrinsic density profile,

where $xi,zi=xi,yi,zi$ are the atom positions and A = L_xL_y is the transverse area. In its original form,⁴ CWT sees the interface at the largest wavelength as a step function, whose local position is described by a Gaussian distribution $Pξ∝exp(−(ξ−⟨ξ⟩)2/2w2)$⁠. The square width $w2$ characterises the correlation length of the surface fluctuations and is known to have the taxing property of, in the absence of an external field, diverging in the thermodynamic limit.^6,7 The connection with the singlet density $ρz$ is often accomplished by ignoring the coupling between capillary wave fluctuations and the interfacial structure, resulting in a congruent fluctuation of the intrinsic profile. The mean density profile is then obtained from the convolution over the capillary waves, $ρz=∫dξρ̃z−ξPξ$⁠.⁸ 

Among the different results reported previously, the interfacial Hamiltonian approach (e.g., Refs. 9–12) gives useful information on the wavenumber dependence of surface tension but relies on various phenomenological parameters which are not known beforehand. Furthermore recently, there has been a renewed interest in the analysis of the pressure tensor in inhomogeneous systems, notably by Sega et al. in Ref. 13 for a Lennard-Jones (LJ) fluid and in Ref. 14 for molecular liquids. In these studies, the authors performed a deconvolution of the pressure tensor from the fluctuations in the interface position, thereby obtaining its intrinsic counterpart. The intrinsic profile of single particle properties is trivially obtained by appropriate weighting of the sum in Eq. (1), $ρ̃w(z+ξ)=∑i⟨wiδ(z+ξxi−zi)⟩/A$⁠, where $wi$ is an arbitrary atomic weight. By contrast, the potential component of the pressure, P^c(r), is not easily defined due to the intermolecular interactions.¹⁵ For a pairwise potential, $Pc(r)=1/2∑i≠j⟨Fij∫Cijdlδ(l−r)⟩$⁠, where F_ij is the force between the pair of atoms and C_ij is the path of the line integral connecting atoms i and j.

There is no unique way of defining C_ij, and so there is also no unique way of determining which atom pairs contribute to the force across a given surface element.¹⁶ The Irving-Kirkwood (IK) convention¹⁷ stipulates that the force between two atoms is said to be across if the vector connecting the centres of mass of the two atoms intersects the surface element. The alternative Harasima (H) contour¹⁸ breaks the line connecting a pair of atoms into two mutually orthogonal paths. For planar interfaces, the H representation is seen as particularly convenient because, by construction, the transverse contribution of each pair of particles becomes path independent;¹³ on the other hand, there are always two possible H contours connecting a pair of particles, first transverse and then normal ↑○ →, and vice versa, → ○↑. In the presence of an interface, this distinction is important because it defines the position where the contribution to the total pressure is accounted for. This is usually mitigated by equally partitioning the contribution at each particle’s location, but it does not resolve the fundamental ambiguity.

Furthermore, it is also recognised that the normal component of this contour is known not to be physically meaningful.^13,19 This additional conceptual difficulty is usually reconciled by arguing that mechanical stability, ∇·P = 0, requires uniformity of the normal pressure across the system. As it stands, this assumption can only hold under equilibrium conditions. Gradients in mass, momentum, and energy will generate thermodynamic fluxes that change the liquid structure, vis-à-vis the interfacial properties.^20,21 Physical modeling of coupling between these is key to understanding the behaviour of fluids under non-equilibrium conditions. Examples include the suppression of thermal fluctuations by shear flow,^22,23 Faraday instability at the liquid-vapour interface,²⁴ and Marangoni flow.²⁵ 

The motivation of this work is thus to extend the IK definition of the pressure tensor as a force acting across a non-planar surface, fluctuating as a result of the thermal motions of the fluid. Starting from the continuity equations of hydrodynamics, expressed relative to a fluctuating reference frame, new expressions are derived which are valid for any inhomogeneous fluid far from equilibrium. But unlike the H or IK line integrals, the new contour of integration is a dynamic quantity, reflecting the instantaneous local curvature of the interface. Using the control volume approach²⁶ to the IK equations, we compute the pressure profile of a well-known liquid-vapour LJ fluid. We find the emergence of the surface tension acting exactly across the intrinsic surface, free from the smearing effects of the surface separating the liquid and vapour phases.

In its original form, the IK microscopic expressions for thermodynamic fluxes are written in terms of ensemble averages. Thermodynamic fluxes are, however, based on conservation laws, valid instantaneously over every member in the statistical ensemble. The Dirac delta functions become a convenient operator whose physical meaning only becomes apparent when averaging over the volume.²¹ We therefore drop the ensemble averaging of the density operators in our analysis following the procedure of Todd et al.^15,20 The mass density $ρr,t$ and momentum density $Jr,t$ are defined by

where m_i, r_i, and $ṙi$ are the particles’ mass, position, and total velocity, respectively. Their hydrodynamic formulation²⁶ allows the densities to be expressed in control volume form rendering them more analytically tractable. Given a region of interest, the integral of the local density,

is the total mass inside the control volume. The region interval is arbitrary as long as it is compact. Assuming that one can obtain the intrinsic surface $ξx$ over the xy-plane for a given configuration, the volume can be foliated into regular elements

where the limits of integration are defined as Δr = r⁺ − r⁻. Here, we are interested in the distribution of the pressure as a function of the distance to the interface. This can be achieved by a change of coordinate system and, in accordance with CWT assumptions, we adopt the view that surface fluctuations are independent of interfacial structures. The arbitrary surface $ξx$ is therefore taken to be independent of time, thereby avoiding the appearance of time dependent contributions to the pressure due to the dynamic motion of the surface. Integration of the operators $δr−ri$ over the control volume has the effect of isolating the particles inside the region of interest,²⁶ and Eq. (4) expands to

where the x and y integrals, in the last line, use the sifting property $∫−∞∞dαf(α)δ(α−αi)=f(αi)$⁠; the boxcar function ϑ results from the integral of the Dirac delta between finite limits, $𝜗αi[α−,α+]≡θα+−αi−θα−−αi$⁠, α = x, y, z, and $θα$ is the Heaviside step function. Assuming periodic boundary conditions in the transverse directions, $𝜗xi0,Lx=1∀xi∈(0,Lx)$ and $𝜗yi0,Ly=1∀yi∈(0,Ly)$⁠,

where $xi=xi,yi$⁠. The mean value theorem applied on the left-hand side of Eq. (7) gives a measure of the volume average density $ρ̃(z+ξ)≈∫Vdrρr,t/ΔV$⁠, z ∈ (z⁻, z⁺), with ΔV = A Δz. Using the definition of the Dirac delta function in the zero volume limit, Δz → 0, $δz≡limΔz→0θz+Δz/2−θz−Δz/2/Δz$⁠, the density at a given position relative to the surface, $z+ξx$ is given by

The above equation defines the instantaneous density relative to a given surface. Assuming that one can somehow compute $ξx$⁠, the average over the ensemble of configurations and corresponding surfaces is the intrinsic density profile, $⟨ρ̃z+ξ⟩$⁠.⁴ The above equation deserves particular emphasis because its derivation is based on a microscopic expression of a hydrodynamic quantity, without any assumption regarding the distribution function $frN,pN,t$⁠. It can therefore form the basis for a statistical mechanical derivation of the pressure tensor in a Lagrangian frame of reference relative to the instantaneous surface. Starting from the definition of momentum, Eq. (3), the same process can be followed to obtain an instantaneous momentum density relative to the surface,

The time derivative of Eq. (3) at a given point results in the momentum conservation equation where, for a system of pairwise interacting particles, the microscopic representation of the IK pressure tensor²¹ is given by

where $Pkr,t$ is the kinetic contribution to the pressure, obtained by decomposition of the total velocity $ṙi$ into a thermal component v_i and a streaming part $uri,t$⁠, taken to be zero in the current work, $ṙi=vi+uri,t$⁠. Akin to density [cf. Eq. (2)] and momentum [cf. Eq. (3)], $Pkr,t$ is localised at the atom positions by $δr−ri$⁠. The configurational part of the pressure tensor, $Pcr,t$⁠, is obtained by integrating along the line connecting a pair of atoms, where each point r_λ = r_i − λr_ij is sampled by $δr−rλ$⁠. Starting from the IK form of the pressure in Eq. (10) and following the same procedure used to obtain intrinsic density above, we aim to obtain the IK pressure tensor mapped onto a coordinate system defined by the instantaneous intrinsic surface. The integral of the pressure tensor over a control volume bounded by a surface $ξx$ is given by

where x_λ, y_λ, and z_λ are the components of a point in the vector connecting particles i and j. The integral of λ from 0 to 1 moves along the line and partitions the energy between the different volume elements across which the contour line passes through. This seemingly complicated looking function has the clear geometric meaning that the path between each pair of particles is not a straight line but a curve—unique to each pair of interacting particles—whose coordinates are dictated by the shape of the interface at the transversal coordinates x_λ and y_λ. Approximating the energy in Eq. (11) by its volume average, $P̃z+ξ≈∫VdrPr,t/ΔV$⁠, over the xy-plane, and using the Heaviside definition of the Dirac delta function when Δz → 0, the instantaneous pressure relative to the surface is given by

where $xλ=xλ,yλ$ are the coordinates of the line between a pair of particles transversal to the interface. In the particular case where the intrinsic surface is planar, $ξx=c$⁠, which can be absorbed into the coordinate position z, and the above equation gives the IK expression for the pressure tensor, Eq. (10), averaged over the xy-plane,

Equation (12) is the main result of the current work. Unlike the IK and H conventions, the contour line for each pair of particles is a unique path that reflects the local curvature of the instantaneous interface, setting the current analysis apart from previous studies.^13,14,27,28 Its average over an ensemble of configurations gives a precise definition of the intrinsic pressure tensor, $P̃z+ξ=⟨Pz+ξ⟩$⁠, free from the smearing effect of surface fluctuations. Because its foundations are hydrodynamic, it is not bound to any equilibrium distribution, thereby paving the way to the study of complex fluids under non-equilibrium conditions.

We apply Eq. (12) to the analysis of a liquid-vapour interface of a simple LJ fluid. Surface tension is the defining physical property characterizing the interfacial structure. From a mechanical standpoint, it arises from the anisotropy of the pressure tensor, reflecting the distinct way molecular interactions vary across the system,

where P_N, and P_T are the components of the pressure tensor normal and transverse to the surface (assuming cylindrical symmetry). The link between these two approaches is given by the principle of virtual work, summarised here for completeness. Given a fluid element dV = L_xL_ydz = Adz and a deformation tensor given by $ϵxx=ϵyy=1+ϵ$⁠, ϵ_zz = 1/(ϵ_xxϵ_yy), and ϵ_αβ = 0, α ≠ β and then to first order in the Cauchy strain ϵ, the change in area is δA = 2 ϵ A and the work done by the system is

By virtue of its definition, the transformation preserves the volume, so the free energy cost per unit area is δF/δA = dz (P_N(z) − P_T(z)), where P_N = P_z and P_T = (P_x + P_y)/2. This analysis merits some justification: (i) the approximation is only valid for infinitesimal deformations ϵ ≪ 1, hinting at the fact that surface tension is an equilibrium property, and (ii) any excess pressure in the transversal direction will cost energy wherever its position across the system. What is perhaps unnoticed from this argument is that gradients in density may break the symmetry of the pressure components and are bound to create tension in the direction transversal to the interface. Although this has been observed previously,¹³ one aim of the current study is to quantitatively analyze in detail this behaviour across the interfacial region. In accordance with CWT, taking the interface at the longest wavelength, we can use Eq. (12) to decouple the capillary wave fluctuations and use it to analyze the molecular origins of the surface tension at the interface,

Molecular dynamics (MD) simulations of a liquid-vapour interface are performed with a system of N = 6000 particles inside a simulation cell of dimensions L_x = L_y = 14.0 and L_z = 80.0, where unless otherwise indicated, reduced units are used throughout this work. The system size was chosen to be larger than the minimum value of N = 1000 particles as illustrated by Trokhymchuk et al.²⁹ The shorter dimension of the simulation box was chosen to be larger than L_x ≈ 10 to avoid spurious correlations in the interface fluctuations,³⁰ and a longer dimension L_z > 5L_x was chosen to avoid coupling between the two interfaces in the system. Atoms interact through a spherical truncated and shifted LJ potential with a cutoff radius of r_cut = 2.5,

where r is the distance between a pair of particles and the corresponding energy and length scales are defined by ϵ = 1 and σ = 1, respectively.

Two temperatures are considered, T = 0.6 and T = 0.7. These were chosen to emphasise the effects of a dense fluid with a sharp liquid-vapour interface. Although different values of the triple point temperature of a LJ fluid have been reported (e.g., Mastny and de Pablo give a value of T = 0.694,³¹ while Errington et al.³² report a value of T = 0.560 using the same cutoff as the current work), it has been shown previously³³ that the current pair radial distribution functions of the studied systems do not exhibit evidence of solid-phase structures. The equations of motion are integrated using the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) software package.³⁴ The temperature is controlled with a Nosé-Hoover^35,36 thermostat with mass m_η = 50.0. For each temperature, 200 statistically independent trajectories are equilibrated for 5 × 10⁵ steps using a time step Δt = 0.001. After equilibration, 50 configurations are sampled from each trajectory every 10⁴ time steps over further 5 × 10⁵ steps giving a total set of 10⁴ configurations over which statistical analysis of the interfacial properties is performed.

Analysis of the interfacial region at the molecular scale requires an operational definition of the collective coordinates of the intrinsic surface $ξx$ for a given molecular configuration. This step is accomplished through a percolation analysis to separate the liquid phase—identified by the spanning cluster—from the isolated clusters in the vapour phase. Originally put forward by Stillinger,³⁷ several approaches have been developed to this end.^38–45 In particular, the intrinsic sampling method (ISM) was the first implementation allowing a quantitative comparison between computer simulations and CWT^38,39 and is the method of choice in the current analysis. The ISM can be interpreted as an outlier detection method⁴⁶ that iteratively estimates the amplitudes of a Fourier series representation of the intrinsic surface,

where x = (x, y) is a vector parallel to the surface, $ξ^k$ are the amplitudes associated with each wavevector in the basis set, k = 2π(n_x/L_x, n_y/L_y), with n_x, n_y = 0, ±1, ±2, …, and k_u = 2π/λ_u is an upper wavevector cutoff that sets the lower resolution limit λ_u of the surface. Here the judicious choice within CWT is to set a cutoff wavelength commensurate with the molecular size λ_u ≈ σ.³⁸ A basic assumption behind ISM is that the molecular system consists of a set of surface atoms, or inliers, whose position is described by the assumed mathematical model, Eq. (18), and outliers—representing the atoms in the liquid and vapour phases—which are not accounted by the model. The spanning cluster is obtained by a graph traversal algorithm that identifies the largest set of atoms with at least 3 neighbours at distances below 1.4. From this set, an initial estimate of the surface inliers is obtained by dividing the xy-plane into a 3 × 3 mesh and finding the N_s = 9 atoms with the largest z coordinate in each mesh element. The intrinsic surface is then defined as the surface of least area passing through the selected atoms. This is accomplished by minimizing the objective function

where A = L_xL_y and ω is a constant that specifies the strength of the harmonic restraints tethering the surface to the atomic positions. The method then searches the percolating cluster for the closest atom to the minimal surface that is not part of the set of inliers. The procedure is iterated until an optimum value of the surface density n_s = N_s/A is reached. At the highest level of resolution, the surface atoms are represented in the density profile as a delta function n_sδ(z). Too low a value of n_s results in a unphysical shoulder between the surface and the first layer in the liquid phase of the intrinsic density profile $ρ̃z$⁠, and too large values are characterised by a degenerate mathematical surface that tries to fit all the surface atoms.⁴⁷ For the LJ system at the studied temperatures, it was found n_s = 0.8 to best describe the interface at the molecular scale with a clear peak at the origin followed by a significant oscillatory structure whose amplitude decays toward the liquid bulk density.

Figure 1 depicts the intrinsic profiles at T = 0.6 and T = 0.7, respectively, obtained using the largest wavevector cutoff k_u = 2π. Intrinsic quantities are plotted in a frame centred on the intrinsic surface, with the corresponding fixed grid quantities centred on the location of the appropriate Gibbs dividing surface. For convenience, both surfaces are denoted as being at z = 0 in Fig. 1, and this convention is followed for all other figures (Figs. 1–4) in this work. The removal of the blurring effect of thermal fluctuations results in the clear appearance of the molecular structure at the interface, given by the delta-function at z = 0 representing the surface layer, followed by a strong oscillatory structure whose amplitude decays to zero at the bulk liquid density, and an adsorption peak at z ≈ 1.0 in agreement with previous results.^47,48

The kinetic and configurational components of the intrinsic pressure profile [cf. Eq. (12)] at T = 0.6 are depicted in Fig. 2. The dyadic terms m_iv_iv_i in the kinetic component $P̃k$⁠, being a single particle property, are as a result directly proportional to the intrinsic density through the kinetic equation of state $P̃kz=ρ̃zkBT$⁠, where k_B is Boltzmann’s constant. At constant temperature, these are all equal, emphasising the well-known fact that, in the absence of temperature gradients, surface tension is a configurational property as attested by the configurational terms $P̃c$⁠.⁴⁹ On the liquid phase, away from the interface (viz., z ≲ −4.0), these are equivalent by a symmetry argument, i.e., the configurational pressure tensor is isotropic (the negative value reflects the cohesive energy of the LJ interactions).

In the layers closer to the interface, fluctuations in density break the symmetry of interactions. Figure 3 illustrates this difference in more detail for $P̃N=P̃z$ and $P̃T=(P̃x+P̃y)/2$⁠, at T = 0.6 and T = 0.7, respectively. The transverse pressure $P̃T$ follows closely the oscillations in density, supporting the notion that, at the peaks in the intrinsic density, the larger number of neighbours interacting with an atom in the directions parallel to the interface results in a higher proportion of short ranged repulsive interactions contributing to the extra stress.

On the other hand, in the normal direction, $P̃N$ exhibits a more steady profile up until the region between the first two layers in the liquid phase. While being zero at the peak locations, reflecting mechanical equilibrium of the liquid layers (cf. Fig. 3, viz., z ≈ −1 × 2^1/6 and z ≈ −2 × 2^1/6), there is a clear positive excess normal pressure between these at z ≈ 1.5. The physical argument behind this excess energy resides in the observation that an atom located at this position would experience the short range repulsive forces of the adjacent fluid layers. The resultant instability would force the atom to transit to one of the layers, thus preserving the liquid structure next to the interface.

We note that at larger distances from the interface, the mean normal pressure P_N is equal in the liquid and vapour phase, as can be seen in Fig. 3. Over the surface, the mean divergence of the pressure will still be equal to zero; i.e., mechanical equilibrium must be satisfied as the surface is stationary.

Figure 4 depicts the excess pressure of the fluid, defined by $P̃Nz−P̃Tz$⁠, whose integral over the entire domain, $γz=∫−∞zdz′PNz′−PTz′$⁠, gives the Kirkwood-Buff⁵⁰ definition of the surface tension [cf. Eqs. (15) and (16)]. At the highest molecular resolution, the intrinsic pressure tensor Eq. (12) allows us to ascertain the microscopic mechanism involved in giving rise to the interfacial energy.

We can observe that the surface tension has contributions from both the liquid and vapour facing sides of the (intrinsic) interface. Importantly, however, we note that the depletion region between the first layer of fluid at the surface layer is precisely where the onset of the surface tension takes place. We are also able to see that the adsorption peak—approximately one fluid layer away from the surface—is vital to the contribution of the other half of the surface tension. Crucially, this is truly an adsorption peak rather than a density oscillation around the bulk vapour density (in contrast to the liquid facing side of the interface), due to the attraction of the vapour particles to the liquid (of course, as they are the same atoms). As the density increases and returns to the bulk vapour density, there is a net additive contribution to surface tension.

The physical picture is thus as follows: the surface tension, unlike the traditional mechanical picture—defined as the excess pressure acting across a surface of zero width between the liquid and vapour phases—arises instead across a boundary region spanning the three key layers involved in the intrinsic surface—the adsorbed layer at z ≈ 1.0, the surface layer at z = 0, and the first layer in the liquid phase at z ≈ −1.0. This region is constructed by pairs of layers attracting each other in an interleaved manner, at a distance commensurate with the attractive tail of the Lennard-Jones potential—the adsorbed layer at z ≈ 1.0 with the first liquid layer at z ≈ −1.0 and the surface layer at z = 0 with the second liquid layer at z ≈ −2.0.

The results shown in Fig. 4 demonstrate that the surface tension force has similar contributions on either side of the interface, an observation first given in Ref. 13 but shown here with unprecedented clarity using the intrinsic pressure derived in this work.

We have proposed a generalisation of Irving and Kirkwood’s expression for the pressure tensor free from the smearing effect of thermal fluctuations of capillary waves. This is accomplished by letting the contour of integration become a dynamic quantity, unique to each pair of interacting particles, that reflects the instantaneous local curvature of the interface. Expressing hydrodynamic fields in the control volume format, we are able to compute the intrinsic density profile that matches the profiles reported previously.

Applying the same technique to pressure allows the derivation of an intrinsic pressure tensor, which gives unprecedented insight into the liquid-vapour interface. More importantly, however, the notion of a mathematical surface of infinitesimal width, central to the mechanical definition of the surface tension, looses meaning at the microscopic level where the discrete nature of matter becomes apparent. Instead, we observe the emergence of the surface tension over a finite range of three atomic layers, acting exactly across the intrinsic surface.

The surface tension and other quantities of interest are derived from the continuity equations of hydrodynamics—expressed in a Lagrangian frame of reference relative to the intrinsic surface—and therefore are valid arbitrarily far from equilibrium. This makes the derived expressions particularly useful to the study of non-equilibrium inhomogeneous systems—e.g., Marangoni-driven flows, hydrodynamic instabilities, and bubble nucleation—which are inherently non-local and unsteady. Of particular interest would also be extension of this study to systems in confinement including heterogeneous substrates and nanochannels (e.g., Refs. 51 and 52).

We acknowledge financial support from the European Research Council (ERC) through Advanced Grant No. 247031 and by the Engineering and Physical Sciences Research Council (EPSRC) of the UK through Grant No. EP/L020564.