Pressure Tensor of Nanoscopic Liquid Drops

This study describes the structure of an inhomogeneous fluid of one or several components that forms a spherical interface. Using the stress tensor of Percus–Romero, which depends on the density of one particle and the intermolecular potential, it provides an analytical development leading to the microscopic expressions of the pressure differences and the interfacial properties of both systems. The results are compared with a previous study and agree with the description of the mean field.


Introduction
Calculating the properties of curved interfaces plays a central role in statistical physics of inhomogeneous fluids. The relevance of this issue lies in the different situations in which the phenomenon occurs, for example in the nucleation of nanoscale drops or bubbles, wetting and drying, technology design for environmental pollution, and problems with recovering oil from porous rock [1][2][3][4][5]. Here is where we found interest in the knowledge of the details and first principles of the structural properties of this system. On the other hand, the spherical interface represents the most important case of curved interfaces, and its description has been a real challenge for a long time and still lacks full knowledge of their properties [4,[6][7][8][9][10][11][12]. In this paper, an analytical derivation of the microscopic expressions of the properties for a system of one or more components is performed. Interfaces that are formed in simple fluids owe their existence to the external field, defining a boundary and separating the homogeneous phase. Although this region is of molecular size, for practical purposes, it can be modeled as a two-dimensional surface. In this region, the density is a function through its position, which adds a mathematical complexity to the description of the system. A key quantity to perform an analysis of the first principles is the stress tensor (negative pressure tensor), which captures the details of the non-locality and the geometry of the system according to the density profile. The stress tensor in this system arises from a conservation equation and is related to the most important properties of the system [2,10,[13][14][15][16], the free energy of the system and the surface tension. From its definition, this amount is not the only quantity, because it introduces a gauge freedom. There is a discussion about whether one can avoid or not avoid the ambiguous character of this quantity [11,[17][18][19]. For the purposes of this paper, it is irrelevant if the stress tensor is unique or not, because the measurable physical properties of the system, such as the pressure difference between the two homogeneous phases and the free energy of the system, are independent of any arbitrary parameter. There are several routes to calculate the interfacial properties, depending on whether the focus is mechanical, thermodynamic or statistical mechanical. The first successful mechanical description was performed independently by Young and Laplace [6,7], who established that the condition of mechanical equilibrium must satisfy a liquid drop in equilibrium with its vapor. The result is the well-known expression relating the pressure difference between both phases δP and γ with the surface tension and the radius R of the drop.
where P l and P v correspond to the pressures of homogeneous bulk phases. The radius of the drop is considered large compared with the range of the intermolecular potential; this ensures that we can guarantee that it is homogeneous inside. The analysis of this system in the context of thermodynamics was initiated by Gibbs [8], introducing the ingenious concept of the dividing surface, which bears his name today. The Gibbs analysis was extended later by Tolman [9], who was able to derive the first-order correction because of the curvature, which is known as Tolman length δ. According to these studies, the equation for the pressure difference must be written as: where γ ∞ is the planar surface tension, and the Tolman length can be determined without ambiguity. Following this thermodynamic approach, we can obtain correction terms of higher order in development. However, there is no certainty that they have a physical meaning, because it is not obvious that purely thermodynamic arguments accurately capture information regarding molecular sizes. On the other hand, if the drop is quite small, in the order of a few molecular diameters, a homogeneous value of the pressure inside the droplet cannot be ensured and the thermodynamic analysis is ambiguous. In this sense, the correct path to describe the system regardless of the size of the drop is the statistical mechanical point of view. In following this approach, Helfrich has shown in a general form for surfaces of arbitrary shapes with two radii of curvature R 1 and R 2 that the pressure difference between two homogeneous phases can be expressed in terms of the mean curvature J = 1/R 1 + 1/R 2 and Gaussian curvature K = 1/R 1 R 2 , as [11,20,21]: where γ(J, K) = γ ∞ + κC 0 J + 1 2 κJ 2 +κK, with γ ∞ the surface tension of the planar surface, C 0 is the spontaneous curvature and κ andκ are constants for rigidity. The calculation of the microscopic expressions of the coefficients introduced in this model is a difficult task in the statistical physics of inhomogeneous fluids; although it has been made from different routes, it is still a subject of analysis and discussion [5,12,22].
We emphasize in this study calculating the microscopic expressions of the pressure difference between the two phases of a single and multicomponent system using the stress tensor in the mean field derived by Romero-Percus [23,24]. In this context, we consider that the most relevant previous studies are those by Blokhouis-Bedox [25] and Varea-Robledo [2,26]. Using the Irvin-Kirkwood pressure tensor in terms of the pair density of the spherical surface in the first one and the pressure derived from the tensor Laplacian model square in the second, both expressions of the difference pressure in terms of the density of the spherical surface were calculated. The essential differences of this study with regard to the above are the following: (1) the use of a general expression of the stress tensor in the mean field, which depends on the density of a particle at two different points in the system; (2) the final results are expressed in terms of the profile of the plane with more correction terms; (3) and the extension of the analysis to systems with an arbitrary number of components. To develop our analysis, we assume that the fluid drop is large compared with the range of the interaction potential and use the fact that the difference in the pressure between the two phases of the system is invariant under a shift of the Gibbs dividing surface. We develop the study as follows: the theoretical background on which this paper is based is presented in Section 2. In Sections 3 and 4, the main analysis for single and multiple system components is performed. Section 5 presents an illustrative selection of comparisons of the results of the present theory with other studies. A summary of the results is presented in Section 6.

Mechanical Equilibrium
In a fluid state of liquid-vapor coexistence, the pressure in the homogeneous bulk phases is a scalar quantity; whereby the pressure tensor is equal to the identity tensor I by a scalar. However, the interfacial region of the force acting across a unit area changes in different directions. This is because the interfacial region is of finite thickness and density is a function of the position, whereby the pressure is a tensor. The pressure tensor can be separated into two contributions: an isotropic part and a configurational part. The first contribution is well defined, but not the second, which introduces an element of arbitrariness that is manifested through a gauge freedom. As the system is in thermodynamic equilibrium, the condition of mechanical equilibrium is satisfied, which, in the absence of external forces, is expressed as ∇ · p = 0, [10,24]. For a spherical interface, based on symmetry arguments, we can say that this depends only on the radial distance and basically has two components, a normal p n (r) and a tangential p t (r) [2,10], which is: In this study, we carry out the formal derivation of this result for the van der Waals model, and we prove that the tangential components are the same. The phenomenological theory affirms that once the normal and tangential components of the pressure tensor have been obtained, all of the thermodynamic information of the system can be calculated. In what follows are shown some relations that may be established for these purposes [2,10]. The equilibrium condition implies the relation: where r denotes the radial direction. In this equation, the cases k = 0, 1, 2, 3 lead to quite different relations useful in the analysis of the interfacial region. By integrating the former equation with k = 0, it is obtained: This expression gives the value of one of the most important measurable physical properties of thermodynamics of an inhomogeneous fluid. To calculate this quantity, we need to know the microscopic expressions of the tangential and normal pressure tensor components. These are provided in the next section.

Results of Density Functional Theory
The pressure tensor can be derived from conservation principles of classical physics, which must satisfy a closed system and are governed by Noether's theorem. The expression obtained for the van der Waals model of the interfacial region is given by [16,23,24]: where α and β indicate the Cartesian components and λ is a parameter indicating gauge freedom in the stress tensor.w(| r |) is the direct correlation function of the fluid in the uniform phases and, in this approximation, is the attractive part of the interaction potential assumed to be of short range ξ B [1,10]. This is the most general expression for the pressure tensor for non-local system, from which can be derived the corresponding gradient and the Laplacian square approaches [1,10]. Now, we express the components of the pressure tensor in spherical coordinates and use the fact that, for this geometry, the density profile depends only on the magnitude of the argument, that is, ρ( r ) = ρ(| r|). We chose the base of vector r, which is the same for the operators involved; we express the vector r in the basis of vector (r,θ,φ). In general, if (r, θ, φ) represents the coordinates of the position vector r and (r , θ , φ ) indicates the position corresponding to the vector r , the transformation is given by: r = r cos θ cos θ + sin θ sin θ cos(φ − φ) r + r sin θ cos θ cos(φ − φ) − cos θ sin θ θ +r sin θ sin(φ − φ)φ.
The angle between the vectors r and r depends on θ, θ , φ, φ ; this dependence is captured in the argument of the density profile, which makes the method of calculating the stress tensor components quite tedious. However, it is possible to make a simplification of the calculation, considering the symmetry properties of the spherical surface, because it has the particularity that any vector that locates a point on the surface is normal. To simplify the calculations without loss of generality, we can choose the vector r in the z direction, which has transformation vector r to another basis: Now, the index represents the components of the involved quantities in spherical coordinates; summation over the repeated index is assumed. This applies to both vectors r and r , although the gradient operator acts on vector r. With these simplifications, we find that the normal component is given by: Following the same procedure with the tangential components: It is observed that the σ φφ and σ φφ components are apparently different; the variable φ has different functions. However, the density profile does not depend on this variable; the integral can be performed explicitly, and the following is obtained: Now, we analyze one of the off-diagonal components, particularly p rθ , The first term vanishes, because the density profile does not depend on θ, and the second term can be explicitly rewritten as: Because The same procedure that applies to all off-diagonal components is p rθ = p θr = p rφ = p φr = p θφ = p φθ = 0. Analytically, this shows that the stress tensor in spherical coordinate system is diagonal and the tangential components are identical, as assumed in the phenomenological theory [2,10,11].
For the system of many components that we are considering, we follow a similar process for obtaining the normal and tangential stress tensor components. We use the fact that the system is in equilibrium state and that the equilibrium condition satisfies ∇ · p = 0. The difference between this case and the previous one is that this object depends currently on several components [16], and the result is obtained by performing the operations. For the normal component: For the tangential components: Again, it is seen that the tangential components are explicitly equal, when the volume element is integrated into the φ variable. The expressions obtained are consistent with the case of a single component. In both cases in the system of single and arbitrary components, the pressure tensor is dependent on the λ parameter, which connects two points of the interface and reflects the fact that the stress tensor is not unique.

Microscopic Derivation of the Laplace Equation
Now, we use the results of the previous section and see that the calculated physical quantity is unique in the sense that it does not depend on any arbitrary parameter. We calculated the pressure difference between the homogeneous phases of Equation (3), considering in detail the case of a single-component system and using the previous expressions for p rr = p n and p θθ = p φφ = p t given in Equations (9) and (10).
We may notice that this expression contains the vectors r and r in the arguments of the density profile. The magnitude of the first one | r | indicates the range of the interaction potential, which is short range, whereas | r| is the order of the radius of the drop. According to this | r| | r |, in using this fact, we develop the argument of the density profile in powers of | r|/| r |.
where s = cos θ . Making λ r → −(1 − λ) r , in the relation also to calculate | r + (1 − λ) r |, among other quantities: These relations are contained in the argument of the density profile and enable us to approximate the density profile, which, so far, is accurate. With the above relations, we can develop the profile around r + λsr or around r − (1 − λ)sr in powers of r /r, as the case. For example, for ρ (| r − (1 − λ) r |), we develop to second order: The resulting calculation of the complete density profile is simple, but tedious, as it involves using all previous developments because the argument of this amount depends on the difference of the vectors r and r . Performing operations and simplifications, the approximate expression is obtained for each of the terms involved in the component pressure tensor difference Equation (17). Substituting these relations in Equation (6), an approximate expression for the pressure difference depends on the arbitrary parameter λ and the density profile, and the interaction potential between the particles is obtained: where δP (2) is the contribution of the second order, which is shown in Appendix A.
In the above expressions, we can observe that the integrals over r go from zero to infinity. Previously, it was assumed that | r| | r |, and the expressions above may seem as a contradiction. However, this can be done without problems, because the kernel ω(| r |) vanishes for r ξ B ; therefore, all of the primed variables cannot became greater than ξ B .
Although the density profile has been approximated by expansion in powers of the r /r, this amount is still the density of the spherical surface, and it contains all of the information of the geometry of the system. The dependence of the λ parameter is an indication of the non-local behavior of the system, because this parameter relates two points on the surface. Now, we make the choice of the Gibbs dividing surface at r = R. However, because the normal coordinate r varies across the interfacial region, to capture the physics of this region, we make r = R+z, which indicates that we take all of the neighboring regions to the Gibbs dividing surface. With this change, the density profile and its derivatives can be expanded in powers of the inverse of the radius of the curvature [28], as follows: where ρ i (z), with i = 1, 2, · · · , t being correction terms because of the curvature, whereas ρ 0 is the profile of the flat surface. In this development, it is assumed that the radii of the fluid drops are large compared with the range of the interaction potential, so that, locally, the density profile can be expressed as the density of the plane in more correction terms. Similarly, it will develop the other densities and its derivatives. The inverse of r can be expanded in powers of the quotient z/R [28], as follows: Similarly, we develop the other radii with larger powers. Incorporating the effects of these changes and making developments term to term, the result is obtained: In this result for the pressure difference, it can be noted that the density profile depends explicitly on the λ parameter, which expresses the fact that the pressure tensor is not unique. However, the relevant physical properties, such as free energy, the pressure difference between the phases, must be independent from this value. To remove this parameter, a strategic shift of the dividing surface of the Gibbs dividing surface [25] is necessary; in this case, it is z → z − λsr . The shift can alter the pressure tensor components, but not the pressure difference, because the choice election positing of the Gibbs dividing surface is a trick to define the properties of the interface relative to a surface; clearly, the physical quantities are invariant to the choice or shifting of the same. In this manipulation, we should add the terms arising from integrating extrapolated bulk regions [25]. However, this is not necessary, because all terms in the integrand contain derivatives of the uniform profile of the homogeneous phases and are canceled. In other words, our result is clearly invariant from its origin, because the pressure tensor only captures information used in the interfacial region [14]. After this shift, we integrate with respect to λ, performing integrations by parts with respect to z with some manipulations and simplifications. The result is: To perform a more detailed simplification and to obtain a simpler expression, we make a separation in Cartesian coordinates; this separation is possible because we have already done a profile development in powers of the inverse of the radius of the curvature. This indicates that the radii of fluid drops are large compared with the range of the interaction potential, and the correction terms indicated are measured from the plane. With this purpose, we identify r s = z , r 2 (1 − s 2 ) = R 2 = x 2 + y 2 and 2πr 2 dr = d r = d R 2 dz , where R is a vector in two dimensions. Furthermore, introduce the auxiliary functions W (R 2 + z 2 ), which also depend on the range of the potential and are related toω(| r |), as: z ω(R 2 + z 2 ) = d/dz W (R 2 + z 2 ) and also try the next relation: Integrating by parts with respect to z and considering the definitions and identifications previously mentioned, the final result is obtained: Comparing with Equation (3), for this geometry J = 2/R and γ(J.K) = γ(R), the result is: The previous result can be written as γ(R) = γ ∞ (1 − 2δ/R + · · · ). To identify explicit microscopic expressions of the planar surface tension and Tolman length, it is noted that the term of lowest order of δp is in agreement with the prediction of the microscopic expression for the surface tension of a flat surface [24], whose expression is: The following term provides the microscopic expression of Tolman length, The previous microscopic expressions for the surface tension and the first term of Tolman length are identical to that predicted using the density functional theory [29]. The difference in Tolman length would be due to R not being the equimolar radius. The main result Equation (29) is expected to be in agreement with the predictions of other points of view in the same level of approximation, because the pressure difference is a measurable physical quantity and is not altered by the ambiguity of the pressure tensor. Below, we make a comparison of the results with equivalent prediction.

Results for Mixtures
Now, we perform the analysis for a fluid mixture of an arbitrary number t of components. The difference with a single component system is that, currently, there are multiple coexistence regions in this system, so we cannot directly apply the analysis as in the previous section, except in a particular case. At high temperatures (compared with the temperature of mixing in the liquid state), the system is in a mixed vapor phase if the temperature is decreased in a controlled manner, and the system separates into two phases: mixed liquid-vapor or a rich or poor phase in all of the system components. In this particular case, there is a single region of coexistence [16,27], in which we can associate a Gibbs dividing surface. However, if the temperature continues to decrease, multiple coexistence regions arise, and it is not possible to apply the analysis of the previous section to the entire system. In this case, we must apply the analysis in part for each coexistence region, although at the end, we can calculate a net pressure difference because of all of the coexistence regions.
We initially consider the case when the system defines a single region of liquid-vapor coexistence; in this case, the Gibbs dividing surface is the same for all species. In developing this surface, we introduce ρ i (R + z + λsr ) = ρ j (R + z + λsr ). Substituting the normal and tangential pressure tensor component Equations (33) and (34) in the expression for the pressure difference, Equation (6), we perform the developments, approaches and manipulations similar to the one component case. The final result is obtained for the mixtures: From here, using Equation (3), we can identify the corresponding expression to the surface tension γ(R).
When the system has several coexistence regions, Equation (6) should be considered for every coexistence region; in particular, for the coexistence regions with density ρ i ( r), we have: where r(i) indicates that the associated radius of the surface, which defines the species i, p i n (r), p i t (r) is the component of the associated pressure tensor to balance this coexistence region. It is important to mention that these components only partially satisfy a balance equation and do not contain the information of the whole system [16]. The explicit expressions are as follows: For the tangential components, when substituting these components in Equation (32), we must take into account that the profiles are developed as: Performing all theoretical developments and simplifications, the result is obtained for the more general case: Clearly, this result is the most general expression for the surface tension γ(R i ) obtained, which is fully self-consistent with the previous results derived in this and the previous section. It is important to note that this expression captures the pressure difference of the interfacial region defined by the ρ i species alone. However, the system can simultaneously have other regions of coexistence, which is reflected in the density ρ j and radii R j . To account for the difference in pressure of all interfacial regions, we must add on index i, which involves using the information in the equation of the balance of the entire system.
The result is consistent with the single-component case, which is confirmed when considering t = 1 and also shows that the microscopic result is independent of arbitrary parameters, because it should be for any physical property. However, the physics contained in this expression is wider. With this result, we can perform numerical analysis in which the components of one region of coexistence can change, while the other components remain fixed; the segregation phenomena, among other aspects of the system, can also be studied.

Comparison with Previous Studies
In this section, a comparison is made with a previous study. For this purpose, we believe the most relevant are: Robledo-Varea [2,26] and Blokhouis-Bedeux [25], the first one with the square Laplacian model (SML) and the second using the Irvin-Kirkwood tensor. Below, we proceed to the predictions of these views.

Pressure Tensor in a Squared Laplacian Model
The free energy functional from this model is derived from the mean field model [1], supposing that gradients on density are small, such that ρ( r 2 ) can be developed around ρ( r 1 ). This gives a free energy function, depending on the density information at only one point. Considering the development level at the square of the Laplacian, we obtain the following expressions: where quantities f 0 , A(ρ), B(ρ) have expressions: and, B(ρ( r)) = 2kT 5! d r r 4 c( r ; ρ( r)).
In this expression, λ is a Broglie thermal wave length, c( r ; ρ( r)) is the correlation function for homogeneous fluids and isotropic with density ρ( r), f 0 is the free energy density for a uniform fluid and coefficients A and B are the second and four moments on c( r ; ρ( r)). The pressure tensor components at this level of approximation [21] are given by: and: The expression obtained for the pressure difference between the two phases is as follows: where P int = P l and P ext = P v .
To compare with our result, in these expressions, we make the changes that we have introduced for ρ( r) in Equation (22) and r in Equation (23), and the result obtained is as follows: from here: The corresponding expressions for the planar surface tension and the Tolman length are: When this quantity is compared with our result Equation (29), we noted that they differ due to the fact that in the SLM result, the densities depend on only one coordinate, while in our case, these quantities depend on the coordinates of two different points. This means that in spite of the approximations introduced, the physical content of the system has not been lost due to the fact that this level still captures the non-locality in the correct form. We should consider that the approximation of free energy that comes from the pressure tensor that we have used is superior to the SLM energy [16,24]. Evidently, our result can be reduced to one similar to those of the SLM, developing the profile of density around one of the coordinates.

Irvin-Kirkwood Pressure Tensor
The explicit microscopic expression for the Irvin-Kirkwood pressure tensor is given by: where r 12 = | r 12 | and ρ (2) (r 1 , r 2 ) is the pair density of the spherical interface using this relation, and Equation (17) gives the following result for the pressure difference: To make a real comparison, the pair density of a spherical interface ρ (2) s (z 1 , z 2 , r) is developed around the pair density of a flat interface, ρ (2) f (z 1 , z 2 , r) [30]: Substituting the terms of the development, we have the expression for the pressure difference δP in terms of the contribution of the plane and corrections due to the curvature. These are: where surface tension γ(R) is identified: It is observed that the corresponding expressions for the planar surface tension and the Tolman length are: In contrast with the previous case with the SLM, these results have a higher approximation level expressed in terms of the pair density of the flat surface and their corrections due to the curvature. It is noted that the superficial tension only depends on the pair density of the flat surface, in agreement with our result. Similarly, Tolman length depends on the pair density of flat surface and the first-order pair density correction of the spherical surface in correspondence with our result. Although the microscopic expressions of these quantities are not identical with our findings, they are consistent in their profile density dependence, so the agreement is satisfactory.

Conclusions
We have calculated the pressure difference between the two phases based on the results of the phenomenological theory and the use of the pressure tensor of Percus-Romero. Although this amount is not unique, because the nature of their origin introduces a gauge freedom, the physical properties are not affected by this ambiguity. The microscopic expressions obtained for the difference of the pressure of the system of one or several components reflect this fact. Although it was not the objective, in the development of the study, we have shown analytically that in the mean field, the pressure tensor is diagonal, and the tangential components are equal for the systems of one or several components, which was supposed by symmetry arguments [10]. On the other hand, despite the fact that the development is performed in a series of powers of the ratio r /R, the expressions of the stress tensor components continue, depending on the information from two different points of the interfacial region and the λ parameter. This means that they do not lose the effect of non-locality as in the approximation of square Laplacian Equations (42) and (43). Therefore, our results are expected to be much more general than those predicted by the SLM, which was confirmed in the previous section. We have also found that the result is in agreement with the mean field models [25]. However, the results are not identical, because in one case, the results are expressed in terms of the pair density, whereas in the other case, it is expressed in terms of the density of a particle, which depends on a couple of points of the interface. On the other hand, the results for systems of one or several components are self-consistent because the prediction for the mixture is reduced to the case of a single component when you take t = 1, and the results with multiple Gibbs dividing surfaces are reduced to the case of a single dividing surface when the choice is made on R i = R j = R. Finally, we emphasize that in this study, we have not only obtained consistent results with those in the literature, but we have gone beyond, considering the fluid mixture with an arbitrary number of components. However, the description of the behavior of the spherical interface in systems of one or several components is far from being a closed issue, because the numerical determination of its properties is not free of ambiguities [4,5,12,22]. We will address this task in the near future. Now, we introduce explicitly the approaches on the density profile ρ and radii r Equations (22) and (23). The result has been separated in two terms δP (2) = δP for a more comfortable manipulation: dθ sin θ dλω(| r |) × r 2 (1 − s 2 )λ z 2 ρ 0 (z − (1 − λ)sr )ρ 0 (z + λsr ) − zρ 0 (z − (1 − λ)sr )ρ 1 (z + λsr ) + ρ 0 (z − (1 − λ)sr )ρ 2 (z + λsr ) − zρ 1 (z − (1 − λ)sr )ρ 0 (z + λsr ) + ρ 1 (z − (1 − λ)sr )ρ 1 (z + λsr ) + ρ 2 (z − (1 − λ)sr )ρ 0 (z + λsr ) + (1 − s 2 ) 3 r 6 ρ 0 (z − sr )ρ 0 (z) + 1 3 s(1 − s 2 )r 3 ρ 0 (z − sr )ρ 0 (z) .
Now, we use the relation Equation (26) to identify the final terms, for example the term: After some algebra, we obtain the second order contribution to Equation (28).

Conflicts of Interest
The authors declare no conflict of interest.