The application of the Young–Laplace equation to a solid–liquid interface is considered. Computer simulations show that the pressure inside a solid cluster of hard spheres is smaller than the external pressure of the liquid (both for small and large clusters). This would suggest a negative value for the interfacial free energy. We show that in a Gibbsian description of the thermodynamics of a curved solid–liquid interface in equilibrium, the choice of the thermodynamic (rather than mechanical) pressure is required, as suggested by Tolman for the liquid–gas scenario. With this definition, the interfacial free energy is positive, and the values obtained are in excellent agreement with previous results from nucleation studies. Although, for a curved fluid–fluid interface, there is no distinction between mechanical and thermal pressures (for a sufficiently large inner phase), in the solid–liquid interface, they do not coincide, as hypothesized by Gibbs.

## I. INTRODUCTION

Under certain conditions (i.e., constant number of particles *N*, volume *V*, and temperature *T*), it is possible to have a spherical phase in equilibrium with another phase around it. There, the Helmholtz free energy *F* is a local/global minimum representing a metastable/stable equilibrium state.^{1–8} This equilibrium implies that *T* and chemical potential *μ* are homogeneous. Thus, ∇*T*(**r**) = 0 and ∇*μ*(**r**) = 0, where **r** is the position vector. However, the number density *ρ*(**r**) = *dN*(**r**)/*dV*(**r**) and the pressure tensor **p**(*r*) are inhomogeneous.^{9,10} By taking the center of mass of the cluster (COM) as the origin and using spherical coordinates,

where **e**_{r}, **e**_{θ}, and **e**_{ϕ} are the unitary vectors. Then, the condition of mechanical equilibrium, ∇ · **p** = 0, implies^{9–11}

In the late 1970s and 1980s, Rusanov, Rowlinson, and Gubbins and co-workers pioneered the application of computer simulations to study fluid–fluid spherical interfaces at equilibrium.^{1,11,12} Lately, there has been a revival in the study of these systems,^{5,6,8,13–27} including also solid–fluid curved interfaces.^{4,7,28–31} In fact, we have recently shown that the stable equilibrium observed in the isochoric-isothermal (NVT) ensemble is an unstable equilibrium in the isobaric-isothermal (NpT) ensemble that corresponds to a maximum in the Gibbs free energy *G*. Thus, nucleation can be studied via both stable and unstable equilibrium as they are two sides of the same coin.^{7,8}

The best thermodynamic description of a system with a curved interface in equilibrium can be found in the book of Rowlinson and Widom.^{9,32} Following in the spirit of Gibbs, one assumes two macroscopic phases that are homogeneous up to the interface and accounts for an additional contribution due to the interface itself. Taking into account that *μ* is homogeneous,

where *γ* is the interfacial free energy, *R* is the radius of the spherical phase, and *p*_{int} and *p*_{ext} are the respective internal and external pressures.

At a molecular scale, there is some arbitrariness in determining *R*. Since *F*, *μ*, *p*_{int}, and *p*_{ext} are fixed, changing *R* also changes *γ*. There are two popular choices for *R*. The first is the Gibbs dividing surface, *R* = *R*_{e}, for which the number of excess particles is zero (meaning that particles belong either to the solid or to the liquid, but not to the interfacial region). The second is the surface of tension, *R* = *R*_{s}, for which *γ* is a minimum (*γ*_{s}).

Actually, by taking the notational derivative (i.e., an arbitrary change in *R* without any physical change in the system), one obtains^{9,32}

By definition, $d\gamma /dR=0$ when *R* = *R*_{s} leading to the celebrated Young–Laplace equation,

Since *γ*_{s} is positive, this equation shows that the pressure inside the spherical phase is higher than outside and the difference depends on *γ*_{s} and *R*_{s}. Simulation studies of fluid–fluid phases by Gubbins and co-workers^{11} and Vrabec and co-workers^{33} have confirmed the higher pressure of the internal phase. However, this equation has not been tested for a solid–fluid curved interface. This is the goal of this letter.

## II. METHODOLOGY

Recently, we simulated several solid clusters in equilibrium with a liquid^{7} via the pseudo-hard-sphere PHS continuous potential^{34} (hereafter, simply HS), which allows us to simulate with the standard molecular dynamics package GROMACS.^{35}

Here, for three selected clusters labeled IV, VII, and VIII in our previous work^{7} (see Ref. 7 for further details on the size of these solid clusters and the way they were obtained using NVT simulations), we launch new trajectories (in the NVT ensemble) saving configurations very often allowing us to compute the pressure tensor. Since the definition of the pressure tensor is locally arbitrary, we choose to use the Irving–Kirkwood^{36} convention in which the forces between two particles act in the line connecting them. Further details can be found in the supplementary material (SM). In addition, density profiles are provided.

The simulation details including the interaction potential, GROMACS setup, and order parameter to label particles as solid or liquid are exactly the same as in our previous work.^{7} We shall use here reduced units. The lengths are given in units of *σ* (i.e., the hard sphere diameter), the time unit is $\sigma m/(kT)$, the densities are given as *ρ* = *N*/*Vσ*^{3}, the pressures are given in units of *kT*/*σ*^{3}, the interfacial free energies are given in units of *kT*/*σ*^{2}, and the chemical potentials are given in units of *kT*. The pressure profiles are computed up to half of the simulation box *L*/2, whereas the density profiles, following Ref. 37, cover the whole system.

## III. RESULTS

The density profile and the normal and tangential components of the pressure tensor are presented in Fig. 1. The values of the densities of the solid and the liquid when they reach a plateau and that of the pressure (far from the interface) are presented in Table I. These are obtained by averaging the data from the corresponding plateaus.

Label . | R_{s}
. | ρ_{sol}
. | ρ_{liq}
. | p_{sol}
. | p_{liq}
. | Δp
. |
---|---|---|---|---|---|---|

IV | 10.791 | 1.0613 | 0.9619 | 12.6046 | 12.7437 | −0.1391 |

VII | 15.20 | 1.0548 | 0.9560 | 12.3053 | 12.4047 | −0.0994 |

VIII | 17.467 | 1.0529 | 0.9541 | 12.2199 | 12.3003 | −0.0804 |

Label . | R_{s}
. | ρ_{sol}
. | ρ_{liq}
. | p_{sol}
. | p_{liq}
. | Δp
. |
---|---|---|---|---|---|---|

IV | 10.791 | 1.0613 | 0.9619 | 12.6046 | 12.7437 | −0.1391 |

VII | 15.20 | 1.0548 | 0.9560 | 12.3053 | 12.4047 | −0.0994 |

VIII | 17.467 | 1.0529 | 0.9541 | 12.2199 | 12.3003 | −0.0804 |

Close to the interface, the normal and tangential components of the pressure tensor are different, albeit, far from it, both are identical. Surprisingly, the pressure inside (solid) is smaller than outside (liquid). This result, in principle, contradicts the Young–Laplace equation. Notice though that having a lower pressure for the solid phase does not violate the mechanical equilibrium condition, which only requires a certain relation between *p*_{N} and *p*_{T} [Eq. (2)].

This is opposite to the fluid–fluid curved interface. Actually, all previous studies on curved interfaces with fluid phases found higher pressure for the internal phase.^{1,11,12} Here, for a solid spherical cluster, we found lower pressure in the internal phase. One may think that this behavior is peculiar for HS, for which there are no attractive forces. However, recently, Gunawardana and Song^{31} have reported a similar behavior for a solid cluster of Lennard-Jones particles surrounded by liquid.

Interestingly, one can already learn the behavior of the curved interface by analyzing the behavior of the pressure tensor for the planar interface. It turns out that, in the interfacial region of a planar interface, *p*_{T} < *p*_{N} for fluid–fluid interfaces^{38} and *p*_{T} > *p*_{N} for solid–fluid interfaces.^{39} Thus, a simple analysis of the behavior of the pressure tensor for the planar interface is sufficient to know if one will have higher or lower pressure in the internal spherical phase [see Eq. (2)]. It is also interesting to point out that the pressure of the external phase, *p*_{ext}, is identical to the average pressure of the system, ⟨*p*⟩, obtained from the virial equation applied to the entire system provided that the normal and tangential components are identical at *L*/2, as demonstrated in the supplementary material.

Although one must accept the fact that the pressure inside a solid cluster is smaller than the pressure outside, the consequence of this appears to be dramatic as this would imply (apparently) from the Young–Laplace equation that *γ* is negative.^{40} The Young–Laplace equation is explained in any textbook of physics, and now, we have a problem about how to use it in the case of a solid cluster surrounded by liquid. How to reconcile the results of this work with the Young–Laplace equation? The key was provided by Tolman in his celebrated paper discussing the variation of *γ* with *R* in a droplet. In particular, there is a remark by Tolman,^{41} which we believe is highly important in this context. The remark is as follows (adapting his notation to this paper):

“*In applying Eqs. (2.2) and (2.3) to very small droplets, it is to be noted that p*_{int} *and ρ*_{int} *are to be taken as the pressure and density for a large mass of internal phase in a condition at the temperature of interest to give the same value of μ as that of the vapor (*cf. *Gibbs,* Ref. 1 *, p. 253).*”

Notice that Tolman was describing the equilibrium between a droplet of liquid and its vapor. However, it also applies for the solid–liquid interface as we are about to show. Therefore, when using the formalism of Gibbs,^{42} the pressure of the internal spherical phase should not be taken as its actual value but rather from that of a bulk having the same *μ* as the external phase. Similar reasoning was also used by ten Wolde and Frenkel.^{43}

Determining the exact value of *μ* of inhomogeneous systems of high density is very difficult^{44} (notice though some recent progress^{45}). Thus, we do not know the exact value of *μ* for the three systems considered in this work. However, to illustrate our main point, this is not crucial. We shall assume that the external liquid has bulk behavior so that *μ* in the system corresponds to that of a bulk liquid at *p*_{liq}. In the inset of Fig. 2(a), *μ*(*p*) for solid and liquid bulks are presented (obtained via thermodynamic integration^{46} from *p* = 11.648 that is the coexistence pressure^{47} where *μ* is identical in both phases).

Taking, for instance, system VIII where *p*_{liq} = 12.3003 and *p*_{sol} = 12.2199, it is possible to determine *μ* for a bulk liquid phase at this pressure and also the pressure of a solid that has the same value of *μ*, which we found to be $psol\mu =12.37$. The superscript *μ* reminds us that this pressure is not the mechanical pressure of the solid but rather the pressure of a bulk solid that has the same *μ* as that of a bulk liquid at *p*_{liq}. We shall denote this as the thermodynamical pressure (as opposed to the mechanical pressure *p*_{sol}). Notice that $psol\mu >pliq>psol$. This finding suggests that the cluster must be different from a perfect bulk at the same pressure *p*_{sol}; otherwise, *μ* could not be homogeneous and no equilibrium could be reached. In order to find differences, we followed the evolution of the closest particles to the COM finding that the solid cluster is a “living” structure that can melt in a certain region and grow in another while keeping the size approximately constant. As can be seen in Fig. 2(b), the selected particles ended up quite close to the interface and some of them changed their neighbor. This is likely due to the presence of vacancies in the cluster that lead to a relative diffusion. By computing such diffusion for the cluster as well as for solid bulks with and without vacancies at $psolVIII$, we could estimate the cluster to have about one vacancy per four thousand particles (1/4000). More relevant is the distribution of distances from a given particle to its 12th closest neighbor considering only particles that fulfill COM < *r* < 10*σ* in order to avoid surface effects (setting the upper limit in 7*σ* did not produce any difference). As shown in Fig. 2(c), such distribution is shifted to higher distances with respect to a bulk with and without vacancies indicating a tiny expansion in the cluster lattice. Further work is needed to completely understand this. Nevertheless, it is clear that the cluster is not identical to a bulk solid with or without vacancies at the same mechanical pressure.

Therefore, for the solid–liquid interface, the Young–Laplace equation must be written as

In Table II, the value of $psol\mu $ for the three systems considered in this work is presented along with the difference in pressure $\Delta p\mu =psol\mu \u2212pliq$.

Label . | R_{s}
. | $psol\mu $ . | Δp^{μ}
. | 2γ_{s}/R_{s}^{7}
. |
---|---|---|---|---|

IV | 10.791 | 12.8627 | 0.1190 | 0.1164 |

VII | 15.20 | 12.4846 | 0.0799 | 0.0793 |

VIII | 17.467 | 12.3700 | 0.0697 | 0.0694 |

Label . | R_{s}
. | $psol\mu $ . | Δp^{μ}
. | 2γ_{s}/R_{s}^{7}
. |
---|---|---|---|---|

IV | 10.791 | 12.8627 | 0.1190 | 0.1164 |

VII | 15.20 | 12.4846 | 0.0799 | 0.0793 |

VIII | 17.467 | 12.3700 | 0.0697 | 0.0694 |

As can be seen, Δ*p*^{μ} is positive. Thus, by using Tolman’s suggestion, one recovers a “normal” Young–Laplace equation. In previous work where the same clusters (IV, VII, and VIII) were studied, we obtained the values of *γ*_{s} from nucleation studies. Since, according to Eq. (6), Δ*p*^{μ} corresponds to 2*γ*_{s}/*R*_{s}, it is of interest to analyze whether our previously reported values of *γ*_{s} and *R*_{s} are consistent with this difference of pressures. As shown in Table II, results are fully consistent. Thus, the physical meaning of 2*γ*_{s}/*R*_{s} obtained for the values of *γ*_{s} and *R*_{s} from nucleation studies^{7,48} is now clear. Note that, for a fluid–fluid interface, there is no difference between $pint\mu $ and *p*_{int} (for a sufficiently large inner phase), whereas, in a solid–liquid system, we could not find such agreement even for very large clusters. We have plotted the difference in pressure between the solid and liquid as a function of 1/*R*_{s} in Fig. 2(a). As can be seen, there is no evidence that this difference can become positive for a certain value of *R*_{s}.

## IV. CONCLUSIONS

In summary, we have computed the pressure tensor for a HS system at constant *N*, *V*, and *T* where one has a stable solid cluster in contact with a liquid away from coexistence conditions. We found that the internal pressure (solid) is lower than the external one (liquid). This would lead to a negative *γ*. However, as suggested by Tolman (and insinuated by Gibbs), defining a thermal pressure for the inner phase, which corresponds to that of a solid with the same chemical potential as the external liquid, allows us to recover a normal Young–Laplace equation, where the pressure of the internal phase is higher, leading to a positive *γ*. The values of *γ* from this scheme are in excellent agreement with recent results from nucleation studies. Thus, for a solid–liquid interface, one should distinguish between the mechanical and the thermodynamic pressure. This distinction is not so necessary for a fluid–fluid curved interface as they are comparable.^{25,51} However, it is crucial in understanding the meaning of the Young–Laplace equation for a solid–fluid interface. Computer simulations have been the key to solve this subtle issue.

## SUPPLEMENTARY MATERIAL

See the supplementary material for a description of the system and its interaction potential, details on the pressure tensor calculation, values for the fitting parameters, and the demonstration that the average pressure in the system equals the external pressure. Additional figures and information can also be found.

## ACKNOWLEDGMENTS

This work was funded by Grant Nos. FIS2016-78117-P and PID2019-105898GB-C21 and FIS2017-89361-C3-2-P of the MEC, by Project No. UCM-GR17-910570 from UCM, and by the USA National Science Foundation under Award No. CBET-1855465. P.M.d.H. acknowledges financial support from the FPI Grant No. BES-2017-080074.

## DATA AVAILABILITY

The data that support the findings of this work are available within the article and its supplementary material.