The pressure tensor across a liquid-vapour interface

Inhomogeneous fluids exhibit physical properties that are neither uniform nor isotropic. The pressure tensor is a case in point, key to the mechanical description of the interfacial region. Kirkwood and Buff, and later Irving and Kirkwood, obtained a formal treatment based on the analysis of the pressure across a planar surface [J.G. Kirkwood and F.P. Buff, J. Chem. Phys., 17(3), (1949), J.H. Irving and J.G. Kirkwood, J. Chem. Phys. 18, 817 (1950)]. We propose a generalisation of Irving and Kirkwood's argument to fluctuating, non-planar surfaces and obtain an expression for the pressure tensor that is not smeared by thermal fluctuations at the molecular scale and corresponding capillary waves [F.P. Buff, R.A. Lovett, and F.H. Stillinger, Jr., Phys. Rev. Lett. 15, 621 (1965)]. We observe the emergence of surface tension, defined as an excess tangential stress, acting exactly across the dividing surface at the sharpest molecular resolution. The new statistical mechanical expressions extend current treatments to fluctuating inhomogeneous systems far from equilibrium.


I. INTRODUCTION
Inhomogeneous systems are characterised by an interface, across which physical quantities vary sharply, over a length scale of a few molecular diameters 1 . 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 (x i , z i ) = (x i , y i , z i ) are the atom positions and A = L x L y is the transverse area. In its original form 4 , CWT sees the interface at the largest wavelength as a step function, whose local position is described by a Gaussian distribution P (ξ) ∝ exp(−(ξ− ξ ) 2 /2w 2 ). The square width w 2 characterises the correlation length of the surface fluctuations and is known to have the taxing property of, in 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 interfacial structure, resulting in a congruent fluctuation of a) ccorreia@imperial.ac.uk b) edward.smith05@imperial.ac.uk c) andreas.nold09@imperial.ac.uk d) d.n.sibley@lboro.ac.uk e) s.kalliadasis@imperial.ac.uk the intrinsic profile. The mean density profile is then obtained from the convolution over the capillary waves, ρ (z) = dξρ (z − ξ) P (ξ) 8 . Among the different results reported previously, the interfacial Hamiltonian approach (e.g.  gives useful information on the wavenumber dependence of surface tension but relies on various phenomenological parameters which are not known before hand. Further recently, there has been a renewed interest in the analysis of the pressure tensor in inhomogeneous systems, notably by Sega et al. in 13 for a Lennard-Jones (LJ) fluid and in 14 for molecular liquids. In these works 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 w i δ(z + ξ (x i ) − z i ) /A, where w i is an arbitrary atomic weight. In contrast, the potential component of the pressure, P c (r), is not easily defined due to the intermolecular interactions 15 . For a pairwise potential, P c (r) = 1/2 i =j F ij Cij dl δ(l − r) where F ij is the force between the pair of atoms, and C ij is the path of the line integral connecting atom 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 16 . Irving-Kirkwood (IK) convention 17 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 18 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 13 ; 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 arXiv:1711.05986v3 [cond-mat.stat-mech] 27 Apr 2018 the position where the contribution to 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 to not 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-avis the interfacial properties 20,21 . Physical modelling 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 24 and Marangoni flow 25 .
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. Unlike 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 26 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.

II. INTRINSIC PRESSURE TENSOR
In its original form, 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 21 . 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 J (r, t) are defined by, where m i , r i andṙ i are the particles' mass, position and total velocity respectively. The hydrodynamic formulation 26 of the these 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 structure. 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 − r i ) over the control volume has the effect of isolating the particles inside the region of interest 26 , and Eq. (4) expands to where the x and y integrals, in the last line, use the sifting property ; the boxcar func-tion ϑ results from the integral of the Dirac delta between x, y, z and θ (α) is the Heaviside step function. Assuming periodic boundary conditions in the transverse directions, ϑ xi [0, The mean value theorem applied on the left hand side of Eq. (7) gives a measure of the volume average densityρ(z + ξ) ≈ V dr ρ (r, t) /∆V , z ∈ (z − , z + ) with ∆V = A ∆z. Using the definition of the Dirac delta function as the zero volume limit ∆z → 0, δ (z) ≡ lim ∆z→0 (θ (z + ∆z/2) − θ (z − ∆z/2)) /∆z, the density at given position relative to the surface, z + ξ (x) is given bỹ 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 + ξ) 4 . The above equation deserves a particular emphasis, because its derivation is based on a microscopic expression of a hydrodynamic quantity, without any assumption regarding the distribution function f r N , p N , 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 21 is given by where P k (r, t) is the kinetic contribution to the pressure, obtained by decomposition of the total velocitẏ r i into a thermal component v i and a streaming part u (r i , t), taken to be zero in the current work, Akin to density (cf. Eq. (2)) and momentum (cf. Eq. (3)), P k (r, t) is localised at the atom positions by δ (r − r i ). The configurational part of the pressure tensor, P c (r, t), is obtained by integrating along the line connecting a pair of atoms, where each point 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 + ξ) ≈ V dr P (r, 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 bỹ 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 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 + ξ) = P (z + ξ) , 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.

A. Intrinsic Pressure and Surface Tension
We apply Eq. (12) to the analysis of a liquid-vapour interface of a simple LJ fluid. Surface tension is the defin-ing 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 x L y dz = Adz and a deformation tensor given by xx = yy = (1 + ), zz = 1/( xx yy ) and αβ = 0 , α = β, 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. The analysis merits some justification: (i) the approximation is only valid for infinitesimal deformations 1, hinting at the fact 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 13 , one aim of the current study is to quantitatively analyse 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 analyse the molecular origins of the surface tension at the interface,

B. Molecular Dynamics Simulations
Molecular dynamics (MD) simulations of a liquidvapour 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. 29 . The shorter dimension of the simulation box was chosen to be larger than the L x ≈ 10 to avoid spurious correlations in the interface fluctuations 30 and 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 31 while Errington et al. 32 report a value of T = 0.560 using the same cutoff as the current work), it has been shown previously 33 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 software package 34 . 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 5 steps using a timestep ∆ t = 0.001. After equilibration, 50 configurations are sampled from each trajectory every 10 4 timesteps over a further 5 × 10 5 steps giving a total set of 10 4 configurations over which statistical analysis of the interfacial properties is performed.

C. Intrinsic Surface Analysis
Analysis of the interfacial region at the shaper 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 37 , several approaches have been developed to this end [38][39][40][41][42][43][44][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 46 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 wave vector 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 ≈ σ 38 . 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 the A = L x L 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 47 . 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 towards the liquid bulk density. Fig. 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, these are both denoted as being at z = 0 in all Figures. 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 mean density profiles were centred at their corresponding Gibbs dividing surfaces ξ .

IV. RESULTS
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 i v i v i in the kinetic componentP k , being a single particle property, are as result directly proportional to the intrinsic density through the kinetic equation of stateP k (z) =ρ (z) k B T , where k B is the 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 termsP c49 . 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. Fig. 3 illustrates this difference in more detail forP N =P z and P T = (P x +P y )/2, at T = 0.6 and T = 0.7 respectively. The transverse pressureP 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 argu- ment 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 pressure will still be equal to zero, i.e. mechanical equilibrium must be satisfied as the surface is stationary. Fig. 4 depicts the excess pressure of the fluid, defined byP N (z)−P T (z), whose integral over the entire domain, γ (z) = z −∞ dz (P N (z ) − P T (z )), gives the Kirkwood-Buff 50 definition of the surface tension (cf. Eq. (15) and Eq. (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 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 Top panel: components of the intrinsic pressure tensor tensor normal,P N =P z (red colour), and transversal to the interface, P T = (P x +P y )/2 (green colour), at T = 0.6. The mean profiles P N (dashed red line) and P T (dashed green line) are plotted for comparison. The intrinsic density, ρ (z) (grey), is plotted on the secondary vertical axis. Bottom panel: intrinsic and average pressure profiles normal and transversal to the interface at T = 0.7. The notation is the same as top panel.
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 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 13 , but shown here with unprecedented clarity using the intrinsic pressure derived in this work.

V. CONCLUSION
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 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 flow, 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). essential in a range of important industrial processes where the surface tension plays a role, such as Marangoni flows, hydrodynamic instabilities and bubble nucleation, which are inherently non-equilibrium, non-local and unsteady.