Local structure of liquid/vapour interfaces approaching the critical point

Investigating the structure of fluid interfaces at high temperatures is a particularly delicate task that requires eﬀective ways of discriminating liquid from vapour and identifying the location of the liquid phase boundary, thereby allowing to distinguish intrinsic from capillary fluctuations. Several numerical approaches require introducing a coarse-graining length scale, often heuristically chosen to be the molecular size, to determine the location of the liquid phase boundary. Here, we propose an alternative rationale for choosing this coarse-graining length scale; we require the average position of the local liquid phase dividing surface to match its flat, macroscopic counterpart. We show that this approach provides additional insight into the structure of the liquid/vapour interface, suggesting the presence of another length scale beyond the bulk correlation one that plays an important role in determining the interface structure.


Introduction
A detailed understanding of liquid/vapour interfaces has proven to be an elusive objective because of the interplay between fluctuations in the bulk and surface capillary excitations. Theories dealing mainly with the effects of bulk fluctuations on the interface structure are exemplified by the approach of van der Waals, 1,2 whereas the theory of Buff, Lovett, and Stillinger 3 was the first to place the capillary fluctuation in focus. When approaching the critical point, fluctuations in the bulk phase are expected to become dominant, while far from it, one expects the intrinsic interface thickness to be of the order of the molecular size for simple liquids and dominated by entropic fluctuations. Several approaches tried to reconcile these two pictures. Weeks, 4 for example, considered the statistical mechanics of columns of fluid separated by the bulk correlation distance x. In this way, all capillary fluctuations with a wavelength larger than x are removed from the statistical average. Further, the local equimolar Gibbs dividing surface is used to define the location of the interface. Several alternative routes have been made since then to overcome the limitations of the capillary wave theory, for example, including a wave-vector dependency into the continuous Hamiltonian that describes the coarse-grained dynamics of the interface, 5-7 exploring nonlocal models [8][9][10] or, recently, extending density functional models to asymmetric interfaces. 11,12 Distinguishing capillary from bulk fluctuations has been proven difficult not only in density functional treatments 6 but also in computer simulations, as the abundance of different numerical approaches attests. For example, columnar averaging has been applied in simulations of polymer mixtures [13][14][15] and water. 16 Other computational strategies [17][18][19][20][21][22] have taken alternative routes to the definition of the liquid/vapour interface. The computational approaches aimed at locating the local boundary between liquid and gas are characterised by the presence of an adjustable parameter (typically, a characteristic length) that influences the determination of surface atoms and, as a consequence, the distinction between capillary fluctuations and corrugations happening below that scale (intrinsic fluctuations). So, for example, the intrinsic sampling method (ISM) of Chacon and Tarazona 17,18 involves a distance cutoff within which new pivot points are searched in each iteration as well as a short-wavelength cutoff for the surface modes; the ITIM and GITIM methods 20,22 introduce the cutoff via the probe-sphere size; the Willard-Chandler method 21 identifies the instantaneous surface with an isodensity surface determined by Gaussian kernels located at the atomic centres, and needs a choice of kernel width and isodensity threshold. These parameters are usually set initially to some physically meaningful quantity (typically, the molecular size) and refined according to some reasonable criterion, including, for example, stability of the results with respect to variation of the parameter itself.
Here, we propose a method to eliminate the ambiguity in selecting this parameter by ensuring consistency between the locally-defined and macroscopic Gibbs dividing surfaces. As it turns out, this characteristic length appears to scale with a different (smaller) exponent than that of the bulk correlation length, x, when approaching the critical temperature and is related to the density contrast between the liquid and the vapour phases. The appearance of another scaling might be surprising at first, but one should not forget that x is not the only fundamental length scale, another example being the distance between the equimolar and the zero excess energy surfaces. 2 The local position of the liquid surface so defined suggests that a new length scale, smaller than the bulk correlation one, might play an important role in determining the interfacial structure of fluids.
Let us stress that according to the commonly accepted theories, only the bulk correlation length x, the capillary broadening D cw , and, in case periodic boundary conditions are in place, the lateral cell size L x enter the description of fluid interfaces close to the critical point. However, in our attempt to define a practical protocol for determining an appropriate coarse-graining length, we found, against our expectations, evidence of a different scaling from that of x. This empirical finding has piqued our interest, and here we present the results of our investigation.
2 Tagging liquid phase, vapour phase, and surface atoms The first step in determining the interfacial atoms using algorithms such as GITIM or ISM is the separation between the liquid and the vapour phases. This classification is usually performed using clustering techniques (falling in the general category of unsupervised algorithms in the language of machine learning, where clustering algorithms play an important role). Far enough from the critical point, it can be safely assumed that almost all molecules are in the liquid phase and that only a few isolated molecules are found in the vapour phase. In these cases, the simplest solution is to perform a cluster analysis based on the presence of neighbours within a given cutoff radius. Molecules with no neighbours are tagged as belonging to the vapour phase, and usually, only one large cluster is found (the liquid phase). This approach can no longer distinguish between the two phases at higher temperatures because chains of molecules connecting the liquid and the vapour phase start appearing. In this case, density-based clustering approaches are more appropriate.
The Density-Based Spatial Clustering of Applications with Noise 23 (DBSCAN) is a popular density-based clustering algorithm. DBSCAN works by defining clusters as dense regions of data points separated by areas of lower density and requires as an input a search radius used to determine the local density and a threshold density value to distinguish points in highdensity regions (called core points) from points in low-density ones. Core points and points within a cutoff radius to a core point are clustered together. In the present case, the highdensity cluster represents the liquid phase.
DBSCAN has been used successfully in the characterisation of fluid phases. 24,25 Here, we use an automatic determination of the threshold density that assumes a bimodal local density distribution 25 and a search radius of 12 Å. Fig. 1 shows the result of the DBSCAN clustering using an ad-hoc generated system, built by filling a cubic simulation box with argon atoms taken from the middle of the liquid phase of an explicit coexistence simulation (details in the Methods section) equilibrated at 153 K, in which a spherical cavity is carved and filled with a patch of atoms taken from the middle of the vapour phase of the same simulation. Fig. 1 shows a cut-out slab passing through the middle of the simulation box with the sphere's location highlighted. The points represent the centres of the argon atoms, coloured according to the way they are clustered by DBSCAN.
Even though the densities of the two phases are pretty close and the curvature radius of the carved sphere is relatively small, the algorithm does an excellent job at distinguishing molecules sampled from the original liquid or vapour phases. Notice that 153 K is, concerning argon, the highest temperature investigated in the present work, where the distinction between vapour and liquid is the most difficult due to the proximity to the critical point. The DBSCAN clustering provides, in this sense, a meaningful and reliable way to determine to which phase atoms belong, on a local basis, when interfaces are present.
Once the liquid phase is determined, one must understand which molecules have to be considered interfacial to determine the capillary fluctuations. We have chosen to use the GITIM algorithm, 22 a generalisation for arbitrarily shaped interfaces of ITIM, 20 as implemented in the Pytim 26 package. In a nutshell, the GITIM algorithm performs the Delaunay triangulation of the atomic or molecular centres in the liquid phase. The algorithm tags as surface atoms those belonging to simplices that can host a sphere larger than a given value R p , called the probe-sphere radius. The algorithm considers the excluded volume of the atoms. Here we use the contact distance, namely, where the pair distribution function starts differing from zero, to determine the excluded volume. Fig. 2 and 3 show molecular dynamics simulation snapshots of the explicit coexistence simulations of argon and water at different temperatures. The simulation details are reported in the Methods section, and the snapshots were produced using the optimal coarse-graining length R p = R cg , as discussed later. Atoms in the liquid phase are shown in blue (the internal ones) or orange (the surface ones), whereas atoms in the vapour phase are shown in red. Hydrogen atoms are always white. The interfacial atoms appearing in the bulk liquid phase mark spontaneous cavitation effects. 27 Notice also that due to the attractive potential of the liquid phase, the equilibrium vapour distribution must show an increased density close to the interface.
Even though the density-based clustering approach is quite robust, it is important to ensure that this accumulation of vapour is not an artefact of the algorithm, especially at high temperatures. We investigated this point further by building a simple model in the convolution approximation. 28 In this approximation, the interface is described at the microscopic scale by the intrinsic profile and at the macroscopic scale by the capillary wave theory, neglecting the coupling between capillary fluctuations and the intrinsic structure. The apparent profile is then obtained by performing the thermal average of the intrinsic profile over the independent capillary modes, resulting in the convolution of the former with a Gaussian probability distribution with width s determined by the continuum Hamiltonian of the capillary wave theory. The effective profile r(z) of an interface centered at z = 0 with intrinsic profile r I (z) is then rðzÞ / ð e Àz 02 =ð2s 2 Þ r I ðz 0 Þdz 0 : Our simple model involves switches from a constant intrinsic density profile r L in the liquid phase (z o 0) to the mean-field solution in the vapour phase, r V e ÀU(z)/(k B T) We have chosen the potential U(z) = U mean (z) + 3U LJ (z) as the sum of the potential obtained from integrating the contribution of a homogeneous distribution of Lennard-Jones centres in the half-space occupied by the liquid, U mean , and the contribution at close range of the three nearest neighbouring molecules, each with a Lennard-Jones potential U LJ : In Fig. 4, we report the result of this empirical model (with effective width s used as a fit parameter) applied to the high-temperature regime, where our approximation is applicable. The results show that despite its qualitative nature, the model can reproduce the sampled density profiles relatively well, including the vapour accumulation at the interface.

Local liquid phase dividing surface
The thermodynamic, macroscopic definition of the equimolar Gibbs dividing surface is the plane at which the excess surface density is identically zero. This condition is often written (liquid phase located in the region z o z eq ) as and defines implicitly the location of the equimolar surface plane z = z eq along the surface normal z for a single component system with density profile r(z) and bulk values r L and r V in the Fig. 2 Snapshots of the argon system at different temperatures. The colours blue and orange mark atoms belonging to the largest cluster of liquid-like molecules, with orange denoting the surface ones as determined using the optimal probe-sphere radius (different for each temperature). The red colour marks the remaining atoms, regardless of whether they belong to smaller liquid droplets or low-density, vapour-like regions.
liquid and vapour regions, respectively. For a periodic simulation cell of edge (0, L) along z, the integration can be carried out explicitly. If the liquid phase of a slab configuration is placed in the middle of the box, the location of the equimolar surface is given simply by In an explicit coexistence simulation in the canonical ensemble, the average density r is known precisely. In contrast, the values of the bulk liquid and vapour densities must be determined by an appropriate fit. As Stillinger noted, 29 the boundary of the liquid phase (however defined) does not need to coincide with the equimolar Gibbs dividing surface. If the density profile r(z) is piecewise constant and equal to the two bulk densities in the respective phases, then z eq is located precisely between the two. When r L c r V , there should be no ambiguity in this statement. If the bulk values and the distribution on the liquid side remain unchanged, but the vapour phase develops an accumulation of material at the interface, as the simple model presented before shows, then r must increase. Consequently, z eq will move toward the vapour phase, and it will not coincide with the location of the liquid phase surface.
We aim to start from this intuitive picture to provide a formulation for the liquid phase dividing surface along the lines of the equimolar Gibbs dividing surface, using the information about the tagged liquid phase molecules. In other words, we want to define the equimolar surface of the liquid phase, which from now on, we will call the ''liquid phase dividing surface'''. To locate the liquid phase dividing surface, we start by computing the profile r L (z) of the liquid phase (which takes the value r L in the bulk). We define the position of the liquid phase dividing surface z L similarly to the equimolar one, eqn (5), through or, for a simulation box with periodic boundary conditions, where r L ¼ ð1=LÞ Ð r L ðzÞdz is the average density of molecules in the liquid phase over the whole box. This relation is valid only if the liquid density reaches zero far from the bulk liquid phase, which was always the case in the systems considered in this work. It is worth noticing that computing the liquid phase dividing surface z L does not require knowledge of liquid interface location, but just of the average profile r L (z). One can think of z L as the average position of a locally corrugated dividing surface z L (x, y) (the local liquid phase dividing surface) that contains all liquid molecules. The calculation of z L (x, y), unlike z L , does require the knowledge of the fluid phase surface.
It is worth reminding the reader that our goal is to match a local definition of the surface that divides the two phases with a Fig. 3 Snapshots of the water system at different temperatures. The colours are the same as in Fig. 2, except for hydrogen atoms that are always white. macroscopic one. The Gibbs equimolar dividing surface, z eq , would not be suitable because, as we discussed, it mixes liquid and vapour even in the case of a perfectly flat interface. 29 This is why we are bound to choose the liquid phase dividing surface. Again, note that separating the two phases (here, using the DBSCAN algorithm) still leaves the question open of where the microscopic liquid phase surface is located. By performing the thought experiment of carving out a portion of material from a bulk liquid, it is easy to understand that, to yield the equimolar condition eqn (7), the location of z L (x, y) must be shifted with respect to the location of the surface passing through the interfacial atomic centres (defined with some suitable interpolation scheme), z S (x, y;R p ), by the typical intermolecular distance in the liquid phase % r r ' 3 4pr L 1=3 ; hence, Fig. 5 reports a sketch to help visualise some of these quantites. As z S (x, y;R p ) depends implicitly on the probe sphere radius, we need an additional condition to find the optimal value R p = R cg , here named coarse-graining length, that satisfies eqn (10). A sensible condition requires that the average position along the surface normal is where the angular brackets denote a suitable average over the points of the surface; here, we simply use the location of the surface atoms. The meaning of eqn (12) is that R cg is the probesphere radius that makes the average position of the local liquid surface coincide with the one computed using the Gibbs criterion eqn (8), once the excluded volume of the atoms is taken into account. Note, again, that z L does not coincide with z eq . As remarked earlier, the Gibbs equimolar surface z eq applied to a capillary-wave-free ideal case would include, on the liquid side, also a number of vapour molecules and is therefore not suitable to identify the microscopic boundary between liquid and vapour.

Results
We tested these concepts on a series of molecular dynamics simulations of water and argon, performed as described in the methods section over a range of reduced temperatures 1 À T/T c from about 0.5 to 0.03, where T c indicates the critical temperature of the respective models. For the Lennard-Jones fluid, we have chosen 30 k B T c /e = 1.3365 (the interested reader could also refer to recent reviews on simulation and equation of state results 31,32 ). For the TIP4P-2005 water model we have used 33 T c = 657.5 K.
For each temperature, we computed the average position of the interfacial atoms for a set of values of the probe sphere radius R p . We then solved eqn (12) numerically for the optimal value R cg using a linear interpolation scheme. In Fig. 6, we report the distance dz(R p ) = z L À hz S (x, y;R p )i À % r between the liquid phase dividing surface z L and the mean local liquid surface position in argon, as a function of the probe-sphere radius, for different temperatures. We determined the optimal probe-sphere radius using a simple linear interpolation to find the zeroes of dz(R p ) In this way, we avoid using any assumption on the scaling of the probe-sphere radius, as the procedure relies only on the determination of the liquid phase (the DBSCAN search radius is set at 12 Å, more than twice as large as the largest R cg ). As we will show, the optimal probe-sphere radius appears to scale with an exponent smaller than that of the bulk correlation length; therefore, R cg can not be identified with the short-wavelength cutoff that is used to define capillary waves, intended as independent surface modes driven by the restoring force provided by the (macroscopic) surface tension.
In Fig. 7, we report the values of the coarse-graining sphere radius as a function of the reduced temperature. When approaching the critical temperature, the values of the coarsegraining radius appear to show a universal behaviour. This behaviour is apparently compatible with R cg B (1 À T/T c ) Àb , where b C 0.32643 is the critical exponent 34 for the order parameter The appearance of a length scale following a scaling different from the one set by the exponent n C 0.63 of the bulk Fig. 5 Sketch of some of the quantities involved in the determination of the optimal probe sphere radius, R cg . The same configuration is reported twice with atoms tagged as inner (blue) or surface (orange) liquid and vapour (red) according to the chosen probe sphere radius R p . The surface passing through interfacial atoms, z S (dashed line) and the local liquid phase dividing surface z L (solid line) are also reported. The liquid phase (blue and orange atoms) is determined by DBSCAN and does not depend on R p . Fig. 6 Distance dz(R p ) between the liquid phase dividing surface z L and the mean location of the liquid phase interface for argon, as a function of the probe-sphere radius R p . The temperature increases from 70 (leftmost curve) to 153 K (rightmost curve). The zeroes of dz yield the optimal probesphere radius, R cg . correlation length might initially seem surprising. Not all distances, however, need to scale with n, another example of a fundamental length being the distance between two dividing surfaces (e.g., the equimolar surface and the zero excess energy surface), which scales with exponent 2 2n À 1 À b. A scaling exponent Àb would imply R cg B 1/Dr, where an inverse square length constant must appear in front of the proportionality relation.
This result underlines the dependence of the coarsegraining length on the density contrast between the two phases. This makes intuitively sense, as R cg must be large enough to ''recognise'' the liquid surface, essentially by preventing the probe sphere from penetrating the liquid phase.
The reader might wonder if defining a length scale that is (at high enough temperatures) bound to be smaller than the bulk correlation length makes sense. Probably, the best way to convince oneself about the soundness of this approach is to look (a posteriori) at the intrinsic density profiles sampled as a function of the elevation z À z S on the liquid phase interface (here z 4 0 denotes the vapour phase). As shown in Fig. 8, even at the highest temperature considered here, one can distinguish liquid from vapour as long as the first peak is considered. At low temperatures, this requires considering a neighbourhood larger than x, which is smaller than the molecular radius, while at high temperatures, the neighbourhood can be smaller than x. This can be appreciated by computing x directly from the decay of the intrinsic density profile in the vapour phase, reported in Fig. 9.
The decay Dr(z À z S ) = r(z À z S ) À r bulk,vap towards the bulk vapour density r bulk,vap can be fitted with an exponential function A exp[ À (z À z 0 )/x], with A and z 0 fitting parameters, to obtain an estimate of the bulk correlation length x.
We report x measured in this way in Fig. 10, showing that it follows indeed the expected scaling (1 À T/T c ) Àn . Notice that we calculated the intrinsic profiles along the macroscopic interface normal. This might become an issue at high temperatures, where the interface undulation becomes more prominent, and the direction of the local normal should be taken into account, 35 but does not seem to be a problem in the present case, at least for the qualitative scaling reported in Fig. 10.
At high temperatures, the bulk correlation length reaches 8 Å, but despite this, close to the interface, it is still possible to Fig. 7 Coarse-graining length R cg as a function of the reduced temperature 1 À T/T c for water (circles) and argon (squares) as well as the scaling law (1 À T/T c ) Àb (black dashed line).   . 9 Decay of the intrinsic density profile of argon, Dr(z À z S ) = r(z À z S ) À r bulk,vap , semi-logarithmic scale, normalised to peak value. The dashed lines are fit to an exponential decay as described in the text. Fig. 10 The bulk correlation length obtained from the fit of the argon intrinsic density profile decay towards the bulk vapour value, and scaling law (1 À T/T c ) Àn . recognise in the intrinsic density profile (Fig. 8) the distinguishing features of liquid and vapour (a different local density), even below this scale. In this sense, we deem R cg a meaningful quantity to identify the location of the liquid phase surface. Finally, let us notice that the fact that R cg is smaller also than the cutoff radius used in DBSCAN is not a contradiction: the latter is chosen to be large enough to resolve better the local liquid and vapour distributions, but this does not influence the size that the probe sphere should have in order not to penetrate the liquid phase.
Nevertheless, the scaling of R cg should be taken with a pinch of salt as we do not have a theory to back it, and we observed it over roughly one decade only. We underline that the highest reported temperatures, especially for the Lennard-Jones fluid, are close to the limit where the system starts fragmenting into multiple droplets (a smaller system might be more stable, but we are constrained to simulate cells larger than the bulk correlation length). In literature, explicit coexistence simulation of simple liquids usually reaches at most 1 À T/T c C 0.01; see, for example, Heyes and Woodcock. 30 In addition, we checked for the influence of finite size effects by simulating a smaller box for argon (5299 atoms in a 44 Â 44 Â 400 Å box). The results of the smaller system are also reported in Fig. 7 and show that finite size effects provide similar results, albeit with a discrepancy from the scaling law that grows with the temperature. It is worth noticing that this finite size effect is not related to a small box size, as the side lengths are still many times larger than the bulk correlation length; rather, the number of atoms is so small that at the highest temperature, the liquid slab width starts being comparable with the bulk correlation length, and the system cannot attain a fully developed density plateau in the liquid phase. Once both the liquid and vapour phases meet the condition of being larger than the bulk correlation length at each temperature (as in the case of the larger system size), the plot of R g falls back to the (1 À T/T c ) Àb scaling.
Again, we would like to stress that the equimolar surface is not expected to coincide with the liquid phase dividing surface, even in simple model systems. As a byproduct of the coarse-graining length investigation, we can now quantify this difference. In fact, the separation between the two surfaces, z eq À z L , grows when the temperature increases, as reported in Fig. 11, for both argon and water. Although it is difficult to draw a clear picture in this case, the gap increase close to the critical point is at least not incompatible with the theoretical expectation of the separation between two dividing surfaces (like the equimolar and the surface of zero excess energy 2 ), with exponent 2n À 1 À b. Of course, one should also consider the possibility of this being a mere coincidence.
To investigate further the interfacial region, we computed the apparent interfacial widths w and w L from the profiles r(z) and r L (z), respectively. We performed nonlinear leastsquare fits to the complementary error function, e.g., . This functional form is the prediction of the convolution approximation 28 of the capillary wave theory with a step-like intrinsic profile and has also been recently proven to be the analytical solution for the threedimensional Ising model. 36 It is worth noting that there is an analytical solution for exponentially decaying intrinsic profiles, 37 even though, in our case, this does not apply on the liquid side until very close to the critical temperature. In passing, we also tried to use the hyperbolic tangent solution of the Cahn-Hilliard theory, 38 but it fits the profiles noticeably worse than the complementary error function. The scaling of w L and w, show in Fig. 12, is compatible with (1 À T/T c ) Àn . Here, one should not forget that the apparent interface width, even if only that of the liquid phase, emerges as the superposition of the intrinsic width, which scales indeed with the exponent Àn, and the capillary broadening D cw in the presence of finite lateral size L x , where a is the microscopic cutoff below which the continuum description of capillary waves loses meaning. In Fig. 12, we also Fig. 11 Separation between the equimolar (z eq ) and liquid phase (z L ) dividing surfaces as a function of the reduced temperature 1 À T/T c for water (circles) and argon (squares). report the scaling law obtained for R cg and for x, showing a crossover value for 1 À T/T c C 0.2, below which (i.e., at high temperature) x is always larger than R cg . In our opinion, the fact that R cg shows the hallmark of a universal behaviour is very suggestive. The physical meaning of R cg is that of a length scale that determines the roughness of the liquid phase surface, as a probe sphere of radius R cg should be large enough not to penetrate the liquid phase. Since the scaling exponent of R cg is smaller than that of the bulk correlation length x, this implies that sooner or later, at high temperatures, x must become larger than R cg ; the liquid phase surface will then be characterised, at distances larger than R cg but smaller than x, by modes that are not independent. This tells us that the features of the liquid surface can be distinguished based on the local density contrast, even if their size is smaller than the bulk correlation length. Capillary waves, intended as independent modes governed by the restoring force of the surface tension, will be observed only at a larger scale. At length scales between R cg and x, one can find surface modes that are not independent. Below R cg , it is not possible anymore to distinguish any surface feature, and everything is characterised by bulk fluctuations. Notice that nonlocal effects are expected to become dominant below the scale of the bulk correlation length but can still be appreciated, for example, for argon close (T C 90 K) to the triple point, at distances in the order of nanometers. 9

Methods
We performed molecular dynamics simulations of argon with the Lennard-Jones potential and water with the TIP4P-2005 potential 39 using the GROMACS package, release 2018.2. 40 We integrated the equations of motion in the canonical ensemble using the Verlet algorithm with a 2 fs integration time step, a Nosé-Hoover 41,42 thermostat with a time constant of 1 ps. We kept water molecules rigid using the SETTLE algorithm. 43 We computed the long-range part of Coulomb and dispersion interactions using the smooth version of the particle-mesh Ewald method 44 using a grid spacing of 1.2 Å, fourth order polynomial interpolation scheme, short-range cutoff of 12 Å, metallic boundary conditions and a relative interaction strength at the cutoff of 10 À5 and 10 À3 for the Coulomb and dispersion terms, respectively.
We simulated 20 133 argon atoms and 13 824 water molecules in rectangular unit cells of size 88 Â 88 Â 300 and 75 Â 75 Â 250 Å 3 , respectively. We started from a slab configuration with normal along the z axis and integrated the equations of motions for 10 ns, the last 5 of which we used for production. Except for the lowest temperatures, each simulation started from the last configuration of the lower-temperature run. We stored configurations to disk every ps for subsequent analysis using the MDAnalysis 45 and Pytim 26 analysis packages. We assigned molecules to the liquid or vapour phases according to the protocol reported first in our study on mixtures with high partial miscibility, 25 which we summarise here. The protocol uses the (typically bimodal) neighbours distribution within a distance of 12 Å to determine for every frame the density threshold for the DBSCAN algorithm. 23 DBSCAN was employed in the context of molecular liquids for highly miscible 25 and supercritical fluids. 24 In DBSCAN, points are tagged as highdensity ones if either the local density (number of neighbours within a given distance r n ) is higher than a threshold one (the so-called core points), or if they lie within r n from a core point (core-reachable points). If r n is large enough, the algorithm's outcome is stable and does not depend noticeably on the radius itself anymore. 25 Here, we define the liquid phase as the largest cluster of high-density molecules. We consider the smaller, disconnected clusters (droplets) as part of the vapour phase to the end of the present analysis. This choice is inconsequential for determining the surface of the largest liquid cluster. We identified each system's surface molecules with the GITIM algorithm in each frame using a series of different probe-sphere radii R p . We collected the histograms of their position, r S (z;R p ), as well as the histograms of the liquid-like and vapour-like molecules, r L (z) and r V (z), respectively. The function x L À z S (R p ) turned out to be always a monotonously decreasing function of R p . Eventually, we determine the zero of the function x L À z S (R p ) by linear interpolation between the two values closest (but with opposites sign) to zero.

Conclusions
The equimolar Gibbs dividing surface is one of the many possible choices to separate two coexisting phases and obtain helpful thermodynamic relations between bulk values and surface-excess quantities, like the adsorption equations. However, even when defined locally, the equimolar surface does not necessarily (and in most cases, it does not) separate all molecules in the liquid phase from those in the vapour one. Here, we investigated the implications of determining the surface separating liquid from vapour using the local density as an underlying criterion. In this context, a coarse-graining length seems to emerge naturally in the definition of the surface of the liquid phase and is found to follow a universal scaling law.
The procedure reported here provides solid ground for determining the otherwise free parameter (like the probesphere radius) in a class of algorithms for identifying interfacial molecules. While we applied this concept to the GITIM algorithm, nothing prevents following a similar procedure for different approaches, like the Willard-Chandler one, where the Gibbs condition could be used to choose the threshold density or the Gaussian kernel parameters. The coarse-graining length is set by requiring the average location of surface molecules to be compatible with the macroscopic liquid phase dividing surface. Strictly speaking, this is not a thermodynamic quantity like the Gibbs equimolar dividing surface, as the set of molecules belonging to the liquid phase cannot be determined by only thermodynamic means. However, this internal selfconsistency removes the need for an assumption on how the probe-sphere radius should scale with temperature and relies only on a recipe to separate liquid-like molecules from vapour-like ones, which we implemented using a density-based clustering approach. Given that the scaling relation's prefactor is unknown, it is necessary to find the optimal probe sphere radius for at least one temperature value when simulating a new liquid.
From a more general perspective, the coarse-graining length R cg sets the scale at which it makes sense to define the liquid surface. Below this coarse-graining length, fluctuations are not living anymore on the liquid/vapour interface but are penetrating the liquid phase. The structure of the liquid/vapour interface is only dictated by bulk fluctuations below that scale set by R cg . On length scales between R cg and x it is already possible to distinguish features of the liquid surface, even though these are composed of modes that are not independent. Eventually, at scales larger than x, independent capillary modes appear. In this sense, the present finding suggests the presence of another fundamental length scale, set implicitly by the order parameter Dr, which is intermediate between strongly correlated bulk fluctuations and independent capillary waves.

Conflicts of interest
There are no conflicts to declare.