# Molecular-sized bubbles in a liquid: free energy of formation beyond the capillarity approximation

Joël Puibasset

<figure><img src="image_1.png" /></figure>

## To cite this version:

Joël Puibasset. Molecular-sized bubbles in a liquid: free energy of formation beyond the capillarity approximation. The Journal of Chemical Physics, 2025, 163 (3), pp.034110. 10.1063/5.0272545. hal-05147414

HAL Id: hal-05147414

https://hal.science/hal-05147414v1

Submitted on 7 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Molecular-sized bubbles in a liquid: free energy of formation beyond the capillarity approximation

Joël Puibasset

Université d'Orléans, CNRS, ICMN, UMR 7374, Orléans, France.

(*Electronic mail: puibasset@cnrs-orleans.fr)

(Dated: 7 July 2025)

We investigate the transient bubbles that spontaneously appear in a simple liquid using molecular simulations. The objective is to deduce the free-energy of formation of the bubbles $W(s)$ from the bubble size distribution $p(s)$ through the hypothesis of a Boltzmann distribution: $W(s) = -kT \ln p(s)$. The bubbles are detected and characterized using a method based on a grid superimposed on the liquid, efficient for bubbles larger than the grid mesh. We first investigate how the results are affected by the mesh choice, and show that using several mesh values allows for the detection of bubbles in a wide range of sizes with minimal computing cost. The free-energy of formation of a bubble can then be deduced for a large range of sizes, with particular emphasis in the region of vanishing bubbles scarcely investigated in previous works. We first show that the usual Boltzmann relation has to be modified when the bubble size is characterized by its volume. In particular, the bubble volume distribution diverges for a vanishing bubble, which should be taken into account before calculating its free-energy of formation from the above formula. An analytical expansion, valid for any interacting spherical molecules, confirms this observation. We then show that the capillarity approximation fails for small bubbles: an extra contribution, linear with the bubble radius, has to be added to the usual quadratic (surface) and cubic (volume) contributions to the free-energy. This extra term most probably relates to the irregular shape of the tiny bubbles.

## I. INTRODUCTION

Among density fluctuations appearing spontaneously in a liquid, the voids or cavities (empty spaces between molecules) and bubbles (region filled with vapor) play an important role regarding its dynamic or thermodynamic properties, for instance its ability to solubilize molecules.¹⁻¹²

These fluctuations also play a crucial role in phase transitions, in particular during the liquid-to-vapor transition of a superheated or stretched liquid.¹³⁻²² In this context, the popular Classical Nucleation Theory (CNT) relies on the idea that bubbles can grow from seeds formed by spontaneous cavities.²³⁻²⁵ Quantitatively, bubbles can reach the size $s$ with a probability $p$ proportional to the Boltzmann factor $\exp[-W(\{s\})/kT]$, where $W(\{s\})$ is the free energy of formation of the bubble, $k$ is Boltzmann's constant, $T$ is the temperature, and $\{s\}$ is a set of variables that characterizes the bubble size and shape (volume, gyration radius, aspect ratio etc.).²,²⁶,²⁷ To simplify the notations, the braces will be dropped, and an expression like “bubble size $s$”, used in a general context, should be understood as including a description of the shape of the bubble. The main objective of this work is to assess the validity of such a relationship between $p$ and $W$ in the early stages of bubble formation. In particular, one has to determine the prefactor that fixes the origin of the free energy, and check that it is well behaved for a vanishing bubble.

In principle, the bubble can adopt any shape, possibly including strong distortions.²⁸⁻³⁶ However, it is generally admitted that, at least for bubbles large enough, the surface energy cost will favor spherical shapes. In CNT, the shape is assumed to be spherical for any size, which is parametrized by its radius $r$. In the capillarity approximation, the free energy of

formation contains surface and volume contributions:

$$W(r) = 4\pi r^2 \gamma + \frac{4}{3}\pi r^3 \Delta, \quad (1)$$

where $\gamma$ is the surface tension and $\Delta$ is for instance a pressure difference in the case of bubble nucleation in a liquid at constant pressure and temperature. While the volume contribution $\Delta$ changes sign when crossing the coexistence line, the surface contribution, which dominates for tiny bubbles, is always positive and generates an energetic cost. The equation then tells us that the first stage of bubble growth is always disfavored, and that tiny bubbles have similar probabilities to occur spontaneously on both sides of the equilibrium line, i.e. in the stable as well as in the metastable liquid, with a probability approaching one for vanishing cavities.³⁷,³⁸

However, nothing can be said about the number of spontaneous bubbles of size $s$ in the liquid, although, in many cases, one is rather interested in this quantity: for instance in CNT for the calculation of the nucleation rate, or in molecular simulations where the number of bubbles is directly measurable. In the context of condensation in a supersaturated vapor or crystallization in a supercooled liquid, the number $N_n$ of nuclei of size $s=n$, where $n$ is defined as its number of particles, has been shown, when clusters are rare, to verify²⁷,³⁸⁻⁴⁴

$$N_n = Ne^{-W(n)/kT} \quad (2)$$

where $N$ is the total number of molecules in the system and $W(n)$ is the free energy of formation of a nucleus of size $n$. Note that $N_n$ is extensive, while the corresponding $N_n/N$ has been referred to as an “intensive” probability.⁴⁰,⁴¹ In the case of bubble formation in a liquid, their size is a continuum, and one defines the density of bubbles of size $s$ as $p(s) = dN_s/ds$ where $dN_s$ is the infinitesimal number of bubbles of size between $s$ and $s+ds$.⁴⁵ $p(s)$ is again extensive, and one expects
<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematic representation of $p(s)$, the density of bubbles of size $s$ in a liquid, as given by Eqs 1 and 3. If $Q$ is constant or varies smoothly with $s$, $p(0)$ is finite (solid line). However, our simulations suggest that $p(s)$ may diverge around $s = 0$ when the bubble size is characterized by its volume (dotted line).</figcaption></figure>

that an equation similar to Eq 2 should hold for $p(s)$:

$$p(s) = Qe^{-W(s)/kT}. \qquad (3)$$

It is however unclear how to generalize the arguments developed in the context of discrete clusters (Eq 2) to the case of continuous ones (Eq 3).²⁷,⁴⁰ As previously, an intensive bubble size density can be defined. For instance, several authors introduced a normalization by the volume $V$ of the system: $\rho(s) = p(s)/V$.²⁹,⁴⁵,⁴⁶ To our knowledge, $Q$, or equivalently $\rho_0 = Q/V$, has generally been taken as a constant, and adjusted so that $W(0) = 0.29$ In the case of stretched water, Menzl et al.²⁹ found $\rho_0 = Q/V = 0.022 \text{ nm}^{-6}$. It is however unclear how to define the free energy of a bubble of size exactly zero ($s = 0$), and, therefore, it is generally assumed to be the limiting value of $W(s)$ for a vanishing bubble ($s \to 0$). This however supposes that the corresponding limit does exist.

Figure 1 gives a schematic representation of $p(s)$ as given by Eq 3 with the hypothesis that $Q$ is constant or varies smoothly around $s = 0$ (solid line). However, the simulation results presented in this work suggest that $Q$ may not be constant, and may even present, in some cases, a divergence at vanishing size (see Fig. 1, dotted line). Note however that, if the condition that bubbles do not overlap is satisfied, the total volume occupied by the bubbles $\int_0^\infty vp(v)dv$ has to be finite ($\le V$), which puts a limit on the degree of divergence of $Q$. Integrability of the distribution imposes an even stronger limit.

The last point to be mentioned is that the exact expression for $p(s)$ is expected to explicitly depend on the chosen parameter used to characterize the bubble size, even when there is a one to one correspondence between them. For example, a spherical bubble may either be characterized by its radius $r$ or its volume $v$. Since the number of bubbles between $s$ and $s + ds$ is the same for $s = r$ or $v$, one has $dN_s = p_r(r)dr = p_v(v)dv$ where we have introduced an index to differentiate the two mathematical expressions for $p$. Therefore,

$$p_r(r) = 4\pi r^2 p_v\left(\frac{4}{3}\pi r^3\right), \qquad (4)$$

and, necessarily, in Eq 3, $Q$ cannot be a constant in both representations. In literature, in most cases, one uses volume distributions to calculate free energies. This is the case in the context of condensation of vapors or crystallization,³³,³⁴,³⁹⁻⁴⁴,⁴⁷⁻⁵² where the natural definition of the nucleus size is its number of atoms, but also in the context of bubble nucleation.²⁹,³⁶,⁵³ On the other hand, the radial distribution is more scarcely used.⁹ Two questions then arise. Is it possible to determine $Q(s)$ for a given choice of size parameter? And, is there a particular choice for which $Q$ would be constant?

Experiments cannot help because the tiny spontaneous bubbles in the liquid are unobservable, but molecular simulations are reliable enough at this scale to provide an efficient alternative. The principle consists of performing a Molecular Dynamics or a Monte Carlo run of a stable or metastable liquid, and following the spontaneous growth and shrink of bubbles. Defining a bubble when this one is large is not difficult. However, when its size approaches atomic scale, it is more involved. There are different methods, all relying on an estimate of the local fluid density, but using quite different approaches, and giving a priori different results for small bubbles.³⁰,³⁶,⁵³,⁵⁴ In this paper, we focus on a popular method proposed by Wang et al.,⁵⁴ denoted as the W-method in the following. The principle consists of analyzing the liquid structure based on a geometric criterion and a grid superimposed on the molecular configuration. For a large grid mesh, the method is fast, but it is unable to capture accurately the shape of the bubbles, and completely misses the smaller ones. Conversely, the smaller the grid mesh, the more precise the characterization of the bubble shape and volume. But the computational cost increases prohibitively, in particular for large systems. In practice, the mesh is generally taken slightly smaller than the molecular size.

The objective of this work is to use the W-method to characterize the bubble size distribution $p(s)$, with a particular focus for $s \to 0$. The molecular model is presented, as well as a brief reminder of the W-method. A careful analysis of the influence of the grid mesh is presented, in particular for small bubbles. It is shown that the bubble size distribution converges toward a definite distribution when the grid mesh is chosen arbitrarily small, which allows one to extrapolate $p(s)$ to vanishing bubble sizes. An unexpected divergence is observed at $s = 0$ when the bubble size is characterized by its volume ($s = v$), which has been confirmed by an analytical calculation that fits very well with the numerical results for tiny bubbles. This divergence however disappears if the bubble size is characterized by its radius, allowing to define a free energy of formation. It is then shown that the capillarity approximation fails to reproduce the free energy of the smallest bubbles.

## II. NUMERICAL DETAILS

### A. Molecular model

The ideas developed in this study are generic and depend essentially on the atomistic structure of the liquid and the
method used to define the cavities or bubbles. We there-
fore focus on a system of particles interacting through the
(12-6) Lennard-Jones potential, depending on the two con-
stants $\epsilon$ and $\sigma$ corresponding to the depth of the potential well
and the molecular kinetic diameter respectively. The poten-
tial is truncated and quadratically shifted at a distance taken
equal to $3\sigma$.⁵⁵⁻⁵⁷ All quantities will be expressed in reduced
units, using $\sigma$ for distances, $\epsilon$ for energies, $\epsilon/k$ for temper-
atures, and $\epsilon/\sigma^3$ for pressures. This model is expected to
apply to simple liquids with spherical molecules. For more
complex molecules, highly anisotropic or flexible, more re-
fined definitions of bubbles or cavities would be required, with
probably different conclusions. For instance, in the case of
network-forming liquids like water, the grid method has to be
adapted.⁵³

The property we are interested in, i.e. the spontaneous for-
mation of cavities, is mainly controlled by thermal activation.
We therefore perform a Monte Carlo simulation, but molecu-
lar dynamics would give identical results. To enhance statis-
tics and favor density fluctuations, we work at the reduced
temperature kT/ε = 1.0 and at a pressure in the metastable
region slightly below the saturating pressure. However, since
the method we use to detect the cavities relies on a mapping of
the liquid with a grid, because of the periodic boundary condi-
tions, the grid mesh has to follow the fluctuations of the sim-
ulation box during an isobaric run. Since we precisely want
to study the impact of a particular choice for the grid mesh,
we prefer to avoid mesh fluctuations and work at constant
box volume in the Grand Canonical ensemble and impose
the chemical potential instead of the pressure. More specif-
ically, we work at ln (zσ³) = -3.20 where z = eμ/kT/Λ³ is
the activity, μ is the chemical potential and Λ is the thermal
de Broglie wavelength.⁵⁵,⁵⁶ Note that the activity is the natural
parameter entering the Grand Canonical Monte Carlo, while
the chemical potential would depend explicitly on the mass of
the molecules. In these conditions, the average reduced pres-
sure is pσ³/ε = 0.026, a value below the reduced saturating
pressure p sat σ³/ε = 0.055 at the working temperature, which
places the system in a slightly metastable state (superheated).
Although the formation of tiny bubbles is slightly favored by
the metastability, no liquid-to-vapor transition was observed
during the simulation runs: the metastability is too weak to
allow the system to overcome the energy barrier. It is also em-
phasized that we have not used any bias to sample the system
(except to acquire information for large bubbles, see below):
the observed bubbles are all spontaneous. A cubic system of
edge L = 18σ has been considered, with periodic boundary
conditions, containing on average 3300 particles.

B. Bubble identification and characterization

The bubbles appearing in the liquid are identified and char-
acterized by estimating the local fluid density. Is this work,
we focus on a method which consists of working on a discrete
grid, and evaluate the impact of a particular choice for the
grid mesh. We anticipate that using different mesh values will
produce different results for the smallest fluid structures, but

<figure><img src="image_4.png" /><figcaption>FIG. 2. Distributions in the vapor and liquid phases of the number of neighbors within the Stillinger sphere defined by the first minimum of the radial distribution function (insert) of the liquid in our thermodynamic conditions (kT/ε = 1 and ln (zσ³) = -3.20)</figcaption></figure>

the largest ones should be unchanged. In this study we have
chosen a four-steps algorithm which has proven to be very ef-
ficient, developed by Wang *et al*.⁵⁴ and denoted as W-method.
We proceed as follows:

(i) The first step consists of identifying liquid-like and vapor-like molecules in the system. For each molecule, one counts the number of neighbors closer than 1.625σ, corresponding to the first minimum of the radial distribution function in our thermodynamic conditions (Stillinger's criterion).⁵⁸ Molecules having five or less neighbors are labeled as vapor-like, while those having six or more neighbors are labeled as liquid-like.⁴⁰ This criterion comes from the fact that the distribution of the number of neighbors calculated in the bulk vapor phase and the bulk liquid phase exhibit two well separated lobes (below 5 and above 6 neighbors, see Fig. 2). In the next steps, the vapor-like molecules are discarded, and only liquid-like ones are considered.

(ii) The simulation box is then divided into cubic cells of edge $l_{cell}$ (see Fig. 3 for an illustration in two dimensions). In practice, $l_{cell}$ is chosen smaller that the particle size to gain insight in the fluid structure at the molecular level, but not too small to keep the procedure computationally fast. A good compromise commonly accepted is $l_{cell} = 0.5\sigma$.⁵⁴ In this work, we explicitly vary $l_{cell}$ to evaluate the influence of this parameter, starting with a significantly larger value $l_{cell} = 2\sigma$, and reducing it by factors of two down to $2^{-4}\sigma$. The volume of the cell is denoted as $\delta = l_{cell}^3$, ranging from $2^{-12}\sigma^3$ to $2^3\sigma^3$.

(iii) The grid cells in the vicinity of a liquid-like molecule (center-to-center distance < 1.625σ) are marked as liquid, while the others are marked as vapor (see Fig. 3).

(iv) We then perform a cluster analysis on the vapor cells to define the bubbles: two cells are considered to belong to the same bubble if they are connected through the faces, edges or vertices, i.e. their center-to-center distance is less or equal to $\sqrt{3}l_{cell}$ (In Fig. 3, the corresponding criterion in two dimensions is $\sqrt{2}l_{cell}$). The size of the bubble may either be defined as the number *n* of cells in the cluster or its total volume *v* = *nδ*.

Step (i) is the most important since it defines the fluid struc-
ture in terms of liquid-like and vapor-like molecules. Steps
<figure><img src="image_5.png" /></figure>

FIG. 3. Illustration of the W-method$^{54}$ in the two-dimensional case. The dark gray disks represent the liquid-like molecules (kinetic diameter $σ$), and the light gray corona are their Stillinger$^{58}$ disks (radius $R$ taken equal to $1.625σ$ like in the three-dimensional case). Panels (a), (b) and (c) correspond to three grids with meshes $l_{cell} = 2σ$, $σ$ and $0.5σ$ respectively. Cells whose centers (black dots) do not fall within any Stillinger sphere are marked as vapor. Bubbles are defined as clusters of vapor cells (highlighted regions with a thick boundary). For a vanishing mesh, this definition of bubbles converges toward the white regions between Stillinger disks.

(ii) to (iv) correspond to the fast discrete numerical procedure used to characterize the bubbles. In the vanishing limit $l_{cell} \to 0$, the discrete procedure (ii) to (iv) gives the exact volume and shape of the bubbles defined as one-piece regions of space out of the Stillinger spheres associated to the liquid-like molecules in step (i). Note that these regions can either contain vapor-like molecules, as expected for a bubble, or be completely empty, in which case it would be more appropriate to call them cavities or voids. However, we will not do the distinction in the following.

### C. Large bubbles and umbrella sampling

The bubbles forming spontaneously in a liquid are in general quite small. The probability of occurrence of larger bubbles decreases very fast, so a biased method is required to study their statistics. We have used the standard umbrella sampling scheme of Torrie and Valleau$^{59}$ to calculate the probability distribution $p(r)$ for $r \ge \sigma$. This method has been a breakthrough for exploring the nucleation barrier of various systems.$^{29,34,40,44,48–50,53,60–63}$ We have also taken advantage of a recently developed method$^{64,65}$ that cures approximations pointed out by Goswami and coworkers.$^{52}$ Although the method applies in principle to any size, we are limited in practice by the shortest dimension of the simulation box ($18σ$ in our case) and the periodic boundary conditions: bubbles

<figure><img src="image_6.png" /></figure>

FIG. 4. Schematic representation rendered with Ovito software$^{66}$ of the bubbles detected with the W-method for a particular molecular configuration. The grayed cells correspond to the vapor-like cells obtained in step (iii) of the method. Panels (a) to (f) correspond to the same molecular configuration analyzed with a grid mesh $l_{cell}$ decreasing from $2σ$ to $2^{-4}σ$ by factors of two.

larger than $r = 4σ$ have a significant probability of interacting with their own images, or of transforming irreversibly into cylindrical bubbles crossing the box. Using larger boxes being out of the scope of this work, the distribution will be limited to $r \le 4σ$.

## III. RESULTS

### A. Molecular configuration analysis

Figures 4 and 5 show the result of the application of the W-method on two particular molecular configurations chosen to illustrate different situations. In each case, the panels (a) to (f) correspond to the six grid meshes considered in this work, ranging from $2σ$ to $2^{-4}σ$ respectively. Tables I and II gather different quantities for the six panels of each figure: the mesh size $l_{cell}$, the number of bubbles $N_B$ in the molecular configuration, and the sizes of the three largest bubbles (when they exist) given in terms of their number of cells $n_i$ or their volume $v_i = n_iδ$.

For the largest mesh of the first example [Fig. 4(a)], one observes only one bubble of size one cell, corresponding to a volume $v_1 = 8σ^3$. Reducing $l_{cell}$ by a factor two [Fig. 4(b)] still gives one cluster, the size being now 12 cells and the volume $v_1 = 12σ^3$. The shape of the bubble looks irregular and wormlike. For $l_{cell} = 0.5σ$ [Fig. 4(c)], two clusters are observed. The largest one, made of 70 cells, is positioned at the same place as the bubble previously detected in panels (a) and (b), and obviously corresponds to the same bubble seen with
<table><caption>TABLE I. Grid mesh $l_{cell}$, number of bubbles $N_B$, and sizes of the largest (1), second largest (2) and third largest (3) bubble, given in terms of their number of cells $n_1$, $n_2$ and $n_3$ respectively, or in terms of their volume $v_i = n_i l_{cell}^3$, for the six grid meshes labeled (a) to (f), and for the configuration depicted in Fig. 4.</caption><thead><tr><th>label</th><th>a</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th></tr></thead><tbody><tr><td>l<sub>cell</sub>/σ</td><td>2.0</td><td>1.0</td><td>0.5</td><td>0.25</td><td>0.125</td><td>0.0625</td></tr><tr><td>N<sub>B</sub></td><td>1</td><td>1</td><td>2</td><td>7</td><td>11</td><td>17</td></tr><tr><td>n<sub>1</sub></td><td>1</td><td>12</td><td>70</td><td>567</td><td>4552</td><td>36193</td></tr><tr><td>n<sub>2</sub></td><td>-</td><td>-</td><td>1</td><td>16</td><td>88</td><td>695</td></tr><tr><td>n<sub>3</sub></td><td>-</td><td>-</td><td>-</td><td>3</td><td>28</td><td>189</td></tr><tr><td>v<sub>1</sub>/σ<sup>3</sup></td><td>8.0</td><td>12.0</td><td>8.75</td><td>8.859</td><td>8.891</td><td>8.836</td></tr><tr><td>v<sub>2</sub>/σ<sup>3</sup></td><td>-</td><td>-</td><td>0.125</td><td>0.25</td><td>0.172</td><td>0.170</td></tr><tr><td>v<sub>3</sub>/σ<sup>3</sup></td><td>-</td><td>-</td><td>-</td><td>0.047</td><td>0.055</td><td>0.046</td></tr></tbody></table>

<table><caption>TABLE II. Same as Table I for the configuration depicted in Fig. 5.</caption><thead><tr><th>label</th><th>a</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th></tr></thead><tbody><tr><td>l<sub>cell</sub>/σ</td><td>2.0</td><td>1.0</td><td>0.5</td><td>0.25</td><td>0.125</td><td>0.0625</td></tr><tr><td>N<sub>B</sub></td><td>2</td><td>6</td><td>7</td><td>16</td><td>9</td><td>16</td></tr><tr><td>n<sub>1</sub></td><td>2</td><td>9</td><td>39</td><td>496</td><td>4052</td><td>32332</td></tr><tr><td>n<sub>2</sub></td><td>1</td><td>3</td><td>24</td><td>49</td><td>390</td><td>3082</td></tr><tr><td>n<sub>3</sub></td><td>-</td><td>1</td><td>4</td><td>31</td><td>249</td><td>2004</td></tr><tr><td>v<sub>1</sub>/σ<sup>3</sup></td><td>16.0</td><td>9.0</td><td>4.875</td><td>7.75</td><td>7.91</td><td>7.89</td></tr><tr><td>v<sub>2</sub>/σ<sup>3</sup></td><td>8.0</td><td>3.0</td><td>3.0</td><td>0.766</td><td>0.76</td><td>0.752</td></tr><tr><td>v<sub>3</sub>/σ<sup>3</sup></td><td>-</td><td>1.0</td><td>0.5</td><td>0.484</td><td>0.486</td><td>0.489</td></tr></tbody></table>

a better definition. On the other hand, the second cluster at the bottom of the simulation box contains only one cell and corresponds to a smaller bubble previously undetected in panels (a) and (b). Reducing further the grid mesh allows for the detection of an increasing number of bubbles in the molecular configuration: for instance, $l_{cell} = 0.0625\sigma$ allows one to detect 17 bubbles or cavities [Fig. 4(f)]. Note however that the newly detected bubbles are quite small, with dimensions of the order of the mesh size that allow to detect them. It is also interesting to notice that, generally, a bubble detected for a given mesh will be detected for smaller meshes. This is in particular the case for the initially detected bubble close to the center of the simulation box [Fig. 4(a)] which appears in all subsequent panels with a better definition. Obviously, the number of cells in this cluster increases, but the size of the bubble measured by its absolute volume $v_1$ converges (see Table I). This is also true for the smaller bubbles detected in Fig. 4 (c) to (f). This leads to the conclusion that the newly detected bubbles in that molecular configuration are 50 times smaller than the largest one. It can also be seen that the shape of the largest bubble converges when the grid mesh decreases. This is also true for the smaller bubbles.

The second example given in Fig. 5 confirms the previous observations, with however some differences. For instance, two clusters are detected for the largest mesh $2\sigma$ [Fig. 5(a)]. The largest cluster made of two cells corresponds to the largest bubble detected in panel (f), with volume $v_1 = 7.9\sigma^3$, while the second cluster made of one cell in (a) corresponds to the second largest bubble in (f) with volume $v_2 = 0.75\sigma^3$. As can be seen, the large grid mesh severely misestimates the volume ratio between the two bubbles: ratio 2 in panel (a), while the

<figure><img src="image_7.png" /><figcaption>FIG. 5. Same as Fig. 4 for another molecular configuration.</figcaption></figure>

ratio is more likely to be 10.5 with a better precision in panel (f) [see Table II, columns (a) to (f)]. Choosing the standard value $l_{cell} = 0.5\sigma$ instead of $2\sigma$ does not cure the problem since the ratio in that case equals 1.6. Another interesting feature in this configuration is that the number $N_B$ of clusters is not necessarily always increasing when the grid mesh decreases, as can be seen in Table II, columns (d) and (e). This is due to the merger of clusters close to each other. In particular, one can notice that the two large clusters in panel (c) actually belong to the same cavity when detected with a better accuracy in panel (f). The last point to be noticed is that the single cell cavity pointed out by the arrow in panel (b) disappears in panel (c), and reappears in subsequent panels. This reveals that cavities of a given size may be undetected when analyzed on a grid with a mesh of comparable size, depending on the position and shape of the cavity with respect to the grid.

These two examples show that the number of cavities detected with the W-method strongly depends on the grid mesh. It is however noticed that the smaller the mesh, the better the accuracy. In particular, the two previous examples show that the largest bubbles are accurately described (in volume and shape) when using the smallest mesh. It looks reasonable to think that whatever the minimal bubble size we want to characterize, it is possible to choose a grid mesh small enough to detect them. In other words, reducing the mesh to zero should allow to characterize all bubbles. However, there remain two important questions.

(i) From a fundamental point of view, it would be desirable that the method somehow converges, in the sense that the characterization of each configuration should reach a definite limit for a grid mesh approaching zero. In particular, the number of detected cavities should converge toward a finite value.

(ii) From a practical point of view, it would be desirable that an accurate enough characterization can be reached with a “reasonable” value of the grid mesh, in particular in terms
<figure><img src="image_8.png" /><figcaption>FIG. 6. Normalized frequency histograms $h_{\delta}(n)$ of the bubble size measured as the number $n$ of cells in the cluster. Large gray bars: $l_{\text{cell}} = 0.5\sigma$; thin blue bars: $l_{\text{cell}} = 0.125\sigma$.</figcaption></figure>

of computing resources: as a matter of fact, reducing the grid mesh is highly demanding in memory and computing time, in particular for large simulation boxes. Taking this constrain into account, is the generally admitted $l_{\text{cell}} = 0.5\sigma$ the best choice?

To answer these questions, the bubble size distribution is carefully analyzed as a function of the grid mesh.

## B. Bubble size distribution

We first focus on the bubble size distribution with the size defined as the number $n$ of cells in the clusters. By construction, $1 \le n \le V/\delta$ where $V$ is the volume of the simulation box and $\delta = l_{\text{cell}}^3$ is the cell volume. One first calculates the frequency histogram, $H_{\delta}(n)$, cumulated over the configurations generated by the simulation and taken with an equal weight. This is equivalent to establishing the histogram based on the ensemble of all cavities detected during the simulation. The histogram is then normalized to one to get $h_{\delta}(n) = H_{\delta}(n)/\Sigma_i H_{\delta}(i)$ which can be interpreted has a probability distribution.

Since distant bubbles in the simulation box are uncorrelated, this ensemble average is equivalent to considering an infinite system, except for the fact that the histogram is truncated at the upper bound $V/\delta$ that explicitly depends on the system size. This limit can be pushed as far as allowed by computer capabilities, using large enough simulation boxes. On the other hand, the first bin $h_{\delta}(0)$ cannot be defined and remains inaccessible to numerical calculations since the smallest bubble detected by the W-method contains at least one cell.

The bubble size histogram $h_{\delta}$ obviously depends on the grid mesh, as can be seen on Fig. 6 that gives the result for the two grid meshes $l_{\text{cell}} = 0.5\sigma$ and $0.125\sigma$. In both cases, the probability decreases monotonically when the bubble size increases. The highest probability is for the first bin (single cell cavities), reaching 38% for $l_{\text{cell}} = 0.5\sigma$ and 19% for $0.125\sigma$, meaning that a large portion of detected cavities have dimensions close to the grid mesh.

Comparison between these normalized histograms requires to represent them as a function of the absolute bubble volume $v = n\delta$.

The result is shown in Fig. 7 for the six grid meshes considered in this study. Because of the large range covered by the bubble volume and probability, the histograms are given in a log-log scale. Whatever the grid mesh, all histograms exhibit a monotonically increasing probability toward smaller cavities. The smallest grid mesh [Fig. 7(f)] suggests that the probability should follow a power law $h(v) \sim v^{-\alpha}$ for small cavity sizes. This is a puzzling result to be discussed below. Visual inspection of the evolution of the histograms with grid mesh does not exhibit anything special with the particular value $l_{\text{cell}} = 0.5\sigma$ commonly adopted in literature.⁵³,⁵⁴ In particular, it is evident that a lot of information below that size is lost.

We are now in a position to calculate the density $p_v(v) = dN_v/dv$ using the finite difference approximation $p_v(v) \simeq \Delta N_v/\Delta v$ where $\Delta N_v$ is the number of bubbles with volume between $v$ and $v + \Delta v$. Taking the natural choice $\Delta v = \delta$, $\Delta N_v$ is given by $H_{\delta}(n)$ with $n = v/\delta$ and $p_v(n\delta) = H_{\delta}(n)/\delta$ for $1 \le n \le V/\delta$. The approximation should be better for bubble sizes much larger than the grid mesh, while it is obviously not appropriate for smaller bubbles. The result, given in Fig. 8, shows that all histograms fall on the same curve independent of $\delta$. It is however interesting to notice that each choice of a grid mesh gives an optimal volume interval for the calculation of $p_v(v)$, between few $\delta$ and hundreds of $\delta$. Below few $\delta$, the grid mesh is too rough to capture the exact volume of the bubble. This can be seen in Fig. 8 where the very first bins of each histogram depart from the main curve. Above hundreds of $\delta$, the bins are too thin, which induces increasing fluctuations that need to be smoothed out. The inset of Fig. 8 shows the data for the three smallest grids in linear scale, beyond the third bin where data are reliable.

## IV. DISCUSSION

### A. The vanishing bubble limit

How do the numerical results for $p_v(v)$ compare to Eq 3? At first glance, the noticeable point is that Eq 3 leads to a

<figure><img src="image_9.png" /><figcaption>FIG. 7. Normalized frequency histograms $h_{\delta}(v)$ given as a function of the reduced cavity volume $v/\sigma^3$ for the six grid meshes $l_{\text{cell}}$ given in the figure.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Symbols: reduced bubble size densities $\sigma^3 p_v(v)$ in log-log scale deduced from the histograms $h_\delta(v)$ shown in Fig. 7, and given as a function of the reduced bubble volume $v/\sigma^3$. The inset shows, in linear scale, the same data where the first three bins have been removed. The color code identifying each grid mesh is the same as in Fig. 7. The thin lines connecting the symbols are guides to the eye. The black dotted line is a fit with Eq 3 with constant $Q$; the black solid line is a fit with Eq 6 and the dashed lines are fits with Eq 7 including all data (black) or only the data for $l_{cell} = 0.5\sigma$ (red).</figcaption></figure>

finite value for $p_v(0)$, while we observe an apparent divergence in Fig. 8. This is highlighted in the inset where the data obtained for the three smallest grids are reported in linear scale. In order to establish the expected behavior of the bubble size density $p_v(v)$ for vanishing $v$, an analytical calculation has been attempted. We first notice that for a vanishing grid mesh, the voids detected with the W-method converge to the regions out of the Stillinger spheres (corresponding to the white regions in Fig. 3). Therefore, the limit of the distribution $p_v(n\delta) = H_\delta(n)/\delta$ for vanishing $\delta$ can be obtained by calculating the density $p_v(v) = dN_v/dv$ of the distribution of these empty regions (inversion of limits). As we are interested in the behavior of vanishing bubbles, we focus on the tiny empty voids that may appear between the Stillinger spheres. Figure 10 gives an example in two-dimensions where a small triangular void is delimited by the Stillinger spheres of three molecules. The Appendix gives the detailed calculation of the volume expansion around zero of $p_v(v)$, in one, two and three dimensions, and proposes a generalization in dimension $d$ (Eq A15). In three dimensions, one has:

$$p_v(v) = \frac{A_v}{v^{2/3}} + \frac{B_v}{v^{1/3}} + C_v + o(1) \quad (5)$$

The details of the fluid-fluid interactions are absorbed into the constants $A_v$, $B_v$ and $C_v$.

In order to make Eqs 1, 3 and 5 mutually compatible, $Q$ cannot be a constant and has to vary around $v=0$ according to $Q = A_v v^{-2/3}$. Discarding in Eq 1 the negligible volume contribution to the free energy compared to the surface term (we are considering small bubbles), the expression for $p_v(v)$ reads:

$$p_v(v) = \frac{A_v}{v^{2/3}} e^{-(36\pi v^2)^{1/3} \gamma / kT} \quad (6)$$

where we have made the hypothesis that the bubbles are spherical.

A best fit with the distributions calculated from simulations gives $A_v = 2.15\sigma^{-1}$, and $\gamma = 0.108\epsilon/\sigma^2$ (solid black line in Fig. 8). As can be seen, the agreement is considerably improved compared to a fit with Eq 3 where $Q$ is supposed to be constant (dotted line). This is particularly true on the leftmost side of the figure where the slope of $v^{-2/3}$ clearly accommodates the data.

On the other hand, for $v \ge 0.5\sigma^3$, the agreement is quite poor. The surface contribution in the exponent varies too rapidly to fit the data. Introducing an asphericity of the bubble would result in an effective surface tension without changing the fit.⁶⁷ Restoring the volume contribution in the exponent would not help since the variation would be even faster. Taking into account the term $B_v v^{-1/3}$ from Eq 5 into the prefactor of Eq 6 does not improve neither. The most efficient way to improve the fit is to introduce a slowly varying contribution $v^{1/3}$ in the exponent:⁶⁸⁻⁷⁰

$$p_v(v) = \frac{A_v}{v^{2/3}} e^{\left[ \frac{(36\pi v^2)^{1/3} \gamma + \kappa_v v^{1/3}}{kT} \right]} \quad (7)$$

A best fit of simulation data with Eq 7 gives $A_v = 5.8\sigma^{-1}$, $\gamma = 5 \times 10^{-3}\epsilon/\sigma^2$ and $\kappa_v = 2.2\epsilon/\sigma$ (black dashed line in Fig. 8). Note that the obtained value for the surface tension is quite small and the corresponding term plays actually a minor role. Omitting this surface tension term gives essentially the same result with indistinguishable fits in Fig 8: $A_v = 6.0\sigma^{-1}$ and $\kappa_v = 2.3\epsilon/\sigma$. As a matter of fact, the new term $\kappa_v v^{1/3}$ largely dominates for small bubbles with $v \le 50\sigma^3$. A physical interpretation of this term will be given below when larger bubbles will be considered. The last point to be mentioned is that performing the same fit with only the data for $l_{cell} = 0.5\sigma$ gives indistinguishable results (red dashed line in Fig 8). The reason is that the information gained with smaller meshes is entirely contained in the $v^{-2/3}$ dependence for vanishing bubbles. This is an important observation regarding the practical use of the W-method: combined with Eq 7, $l_{cell} = 0.5\sigma$ is a good compromise between precision and computational cost.

## B. Volume versus radius densities and free energy profile

In most cases, when the capillarity approximation is invoked and the bubble is supposed to be spherical, its size is characterized by its radius instead of its volume. From Eqs 4 and 7 it is seen that the prefactor in this radial representation is now expected to be a constant, and the new term $\kappa_v v^{1/3}$ is now linear with $r$. In the general case where bubbles are not spherical, one can introduce an effective radius $r_e = (3v/4\pi)^{1/3}$. The corresponding expression for the radius density reads:

$$p_r(r_e) = A_r e^{-[\kappa_r r_e + W(r_e)] / kT} \quad (8)$$

where $A_r = (36\pi)^{1/3}A_v$, $\kappa_r = (4\pi/3)^{1/3}\kappa_v$, and $W(r_e)$ is the capillarity approximation given by Eq 1 where $r$ is taken equal to $r_e$. This equation allows one to define an effective free energy

$$W_{\text{eff}}(r_e) = -kT \ln p_r(r_e) \quad (9)$$
<figure><img src="image_11.png" /></figure>

FIG. 9. Effective free energy profiles $W_{\text{eff}}(r_e)/kT$ given by Eq 9, as a function of the bubble effective radius $r_e/\sigma = (3v/4\pi)^{1/3}\sigma^{-1}$. The data (colored symbols) are deduced from Fig. 8 and Eq 4. The color code identifying each grid mesh is the same as in Fig. 8. The first three bins of each histogram have been omitted. In the inset, all data for $r_e \le 2\sigma$ are shown, while in the main panel, data for $l_{\text{cell}} > 0.5\sigma$ are omitted for clarity. The results from biased umbrella sampling for $r_e \ge \sigma$ are shown as squares in the main panel. The lines are best fits with Eq 8 (solid line), or without $W(r_e)$ (dashed line), or with $\kappa_r = 0$ (dotted line). Typical bubbles are shown for $r_e/\sigma = 1.3, 2.5$ and 4.0.

cleared from the spurious divergence for vanishing bubbles. This can be seen in the inset of Fig. 9 giving $W_{\text{eff}}/kT$ calculated from the simulation data. As previously, the first bins of the histograms are omitted since the corresponding densities deduced from finite difference are erroneous, see Fig. 8. As can be seen, all histograms fall on the same curve. But now, the noticeable point is that the plotted quantity clearly reaches a finite value at vanishing bubble size. It is possible to define and calculate the limit $W_{\text{eff}}(r_e \to 0)$ ($=-3.3kT$ in our case) which corresponds to $-kT \ln A_r$. This allows us to unambiguously define a free energy difference $W_{\text{eff}}(r_e) - W_{\text{eff}}(0)$ that can be interpreted as the free energy of formation of the bubble (that correctly reaches zero for a vanishing bubble). It is emphasized that this is not possible if one uses $p_v$ instead of $p_r$ in Eq 9.

### C. Free energy of formation of a bubble beyond the capillarity approximation

How does the previous definition of the free energy of formation of a bubble $W_{\text{eff}}(r_e) - W_{\text{eff}}(0)$ compare with the capillarity approximation? We have already seen that an extra term $\kappa_r r_e$ has to be introduced to fit the data for small bubbles. What happens for larger bubbles?

We have used the umbrella sampling method to calculate the probability distribution up to $r_e = 4\sigma$.⁵⁹ The results are given in Fig. 9 (squares), as well as the previous data for $l_{\text{cell}} \le 0.5\sigma$. The previous fit with Eq 7 and $\gamma = 0$ (dashed line in Fig 8) is also reported in Fig 9 (dashed line, corresponding to Eq 8 with $W(r_e)$ set to zero).

As can be seen, the simulation data for large bubbles ($r_e \ge 2\sigma$) now depart from the dashed straight line corresponding to the previous fit without the surface and volume contributions.

Re-introducing the surface tension allows one to fit nicely the data over the whole available range (solid line in Fig 9, with $\gamma = 0.02\epsilon/\sigma^2$; the volume term was unnecessary, but would undoubtedly contribute for larger bubbles). For comparison, a best fit with Eq 8 with $\kappa_r = 0$ (dotted line) has also been attempted. This would correspond to considering only the usual capillarity approximation. As can be seen, the agreement is quite poor.

What is the significance of the introduction of the linear term $\kappa_r r_e$? The first point to be mentioned is that the cavity distribution may be method-dependent, in particular for small bubbles. This could influence the magnitude of the linear term: further work is required to evaluate this point. However, in the case of the W-method, a linear term is required to fit the data down to the smallest bubbles. It might be tempting to interpret this term as a consequence of the dependence of the surface tension on the bubble curvature $-2/r$, according to $\gamma(r) = \gamma_{\infty}/(1 - 2\tau/r)$, where $\tau \approx -0.1\sigma$ is Tolman's length.¹⁵,¹⁹,⁷¹⁻⁷⁴ In this context, if one identifies the linear term $\kappa_r r$ with $8\pi\tau\gamma_{\infty}r$ appearing in the expansion for large bubbles of the surface energy $4\pi r^2 \gamma(r)$, one gets $\tau = \kappa_r/(8\pi\gamma_{\infty}) = 1.2\sigma$. This is at odds both in sign and magnitude with the expected value. As a matter of fact, this is not contradictory since Tolman's expansion is for large bubbles while $\kappa_r r$ is introduced for very small ones. From a physical point of view, $\kappa_r r$ is more likely interpreted as an effective term to take into account the fact that small bubbles are not of spherical shape (see snapshots in Fig 9). $\kappa_r$ would correspond to an additional linear energetic contribution associated with surface corrugation.

What is the consequence of neglecting the linear term $\kappa_r$ in the context of phase transitions, in particular for the calculation of nucleation barriers? By definition, the barrier is the free energy of formation of the critical nucleus, i.e. $W(r_c) - W(0)$ where $r_c$ is the critical radius. If the linear contribution is omitted in Eq 8 ($\kappa_r = 0$), i.e. $W_{\text{eff}}$ is replaced by the capillarity approximation Eq 1, a best fit can be performed on the largest bubbles where this approximation is expected to be better ($r_e \ge 3\sigma$ in our case). The result is given in Fig 9 (dotted line). As can be seen, this results in a significant difference for $W(0)$, $5.5kT$ in our case, that reports directly on the nucleation barrier. For comparison, the CNT predicts a reduced barrier around $20kT$ for $kT/\epsilon = 1$ and $\Delta P\sigma^3/\epsilon = -0.03$ corresponding to our thermodynamic conditions: the error then reaches 25%. For larger barriers, more likely to appear in realistic experimental situations, the relative error would be smaller, but not necessarily negligible.

### V. CONCLUSION

Spontaneous tiny bubbles are ubiquitous in liquids, in the stable as well as in the metastable state. They are particularly important in the context of the liquid-to-vapor phase transition as they play a fundamental role in the process of nucleation. In particular, the probability of finding bubbles of the critical size appears in the classical nucleation theory. More generally, the bubble size distribution is generally translated into a free
energy profile through Eq 3, with an associated free energy of formation that can be written in terms of surface and volume contributions in the capillarity approximation (Eq 1) valid for large bubbles. In this work we have focused on small bubbles or cavities, typically of the order of the molecular size. We have shown that:

(i) The W-method⁵⁴ is efficient to detect tiny bubbles if the mesh $l_{cell}$ is chosen small enough. In practice, it is enough to take $l_{cell} = 0.5\sigma$: large bubbles are well characterized, while the size distribution of small bubbles is well accounted for thanks to the following points (ii) and (iv).

(ii) The bubble volume distribution $p_v$ follows a $v^{-2/3}$ law in the vanishing limit, while the distribution $p_r$ calculated for the equivalent radius $r_e = (3v/4\pi)^{1/3}$ converges toward a constant in virtue of Eq 4.

(iii) One can define an effective free energy according to Eq 9. Thanks to point (ii), $W_{eff}(0)$ exists, and it is possible to define the free energy of formation of the bubble as $W_{eff}(r_e) - W_{eff}(0)$, which correctly reaches zero for a vanishing bubble. Despite the fact that, in most cases, Eq 9 seems to be correctly used in literature, it has never been clearly stated, and implications of Eq 4 have never been discussed. It is however emphasized that using $p_v$ instead of $p_r$ in Eq 9 does not allow to define properly a free energy of formation of the bubble.

(iv) Comparison of $W_{eff}$ with the capillarity approximation Eq 1 shows that it is necessary to introduce a linear term $\kappa_r r_e$, that can be interpreted as an effective correction to the fact that bubbles significantly depart from the spherical shape.

As a conclusion, gathering these recommendations in the context of nucleation should allow to improve the calculation of free energy barriers, in particular when the critical bubble is small and the corresponding barrier is low.

## ACKNOWLEDGMENTS

The author acknowledges fruitful discussions with E. Rolley and P. Porion, and is grateful to the CNRS interdisciplinary “Défi Needs” through its “MiPor” program (Project DARIUS) and the financial support of Agence Nationale de la Recherche through the project ANR-23-CE30-0028-03 (NanoCav).

## AUTHOR DECLARATIONS

### Conflict of interest

The author has no conflicts to disclose.

### Author contributions

Joël Puibasset: Conceptualization; Data curation; Writing – original draft; Writing – review & editing.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## Appendix A

This section is devoted to the calculation of the bubble size distribution in the vanishing limit. It is shown that a general relation can be obtained with minimal hypotheses. For clarity, the argument is first presented in details in dimensions 1 and 2. In dimension 3 we give only the important steps, leading to a generalization to dimension $d$.

### 1. Definitions

The first step of the W-method consists of identifying the liquid-like atoms with the Stillinger’s criterion,⁵⁸ which compares the number of neighboring molecules with a given threshold. The liquid-like atoms are then replaced by a sphere of radius corresponding to the first minimum of the radial distribution function (Stillinger’s sphere). These criteria obviously depend on the space dimension $d$: the radius of the Stillinger $d$-dimensional sphere will be denoted as $R_d$, where $R_3 = 1.625\sigma$ in our case of a Lennard-Jones fluid at $kT/\epsilon = 1$. The positions of the Stillinger spheres are denoted as $\mathbf{r}_i$, and the distance between spheres $i$ and $j$ is denoted as $r_{ij}$.

Let us now focus on the conditions to have a vanishing bubble, corresponding to the smallest voids that can be built between the Stillinger spheres. In one dimension, these voids consists in linear intervals between two adjacent and non-overlapping 1D-spheres. In two dimensions, they consist in curvilinear triangles formed by three disks which overlap in pairs, but not all three together. In three dimensions, four spheres are required, which overlap three by three, but not all four together. Let us denote $V_d(\mathbf{r}_{i_0}, \dots, \mathbf{r}_{i_d})$ the volume of the void between the $(d+1)$ $d$-dimensional spheres $i_0, \dots, i_d$ verifying the above mentioned conditions. The probability to have a void of volume equal to a given infinitesimal $d$-volume $v_d$ is proportional to the average number of such voids:

$$p_v(v_d) \propto \int e^{-H(\mathbf{r}_{1}, \dots, \mathbf{r}_{N})/kT} \sum_{i_0, \dots, i_d} \delta(V_d(\mathbf{r}_{i_0}, \dots, \mathbf{r}_{i_d}) - v_d) d\mathbf{r}_{1} \dots d\mathbf{r}_{N} \quad (A1)$$

where the sum runs over all $(d+1)$-uplet $(i_0, \dots, i_d)$ of spheres and $\delta$ is the Dirac distribution that selects those satisfying the condition to form a small void with the given volume $v_d$. The exponential is the Boltzmann weight associated with the statistical ensemble under consideration, chosen to be canonical in this example for simplicity. However, the argument can be generalized to any statistical ensemble. The integration runs over the coordinates of the $N$ spheres in the system.

The vanishing voids between Stillinger spheres are rare enough to be considered as spatially independent. Therefore,
their number distribution follows a Poisson law, and the average number of voids of a given size is proportional to the probability to have exactly one void of the given size. By symmetry, all the (d+1)-uplets of the sum contribute equally. Therefore, the calculation can be done only for the particular choice of spheres 1 to (d+1):

$$p_v(v_d) \propto \int e^{-H'(\mathbf{r}_1, \dots, \mathbf{r}_{d+1})/kT} \delta(V_d(\mathbf{r}_1, \dots, \mathbf{r}_{d+1}) - v_d) d\mathbf{r}_1 \dots d\mathbf{r}_{d+1} \quad (A2)$$

where the integration now runs only on the (d+1) positions of the Stillinger spheres that form the void. This is possible through a redefinition of the Boltzmann weight by integration on the remaining spheres:

$$e^{-H'(\mathbf{r}_1, \dots, \mathbf{r}_{d+1})/kT} = \int e^{-H(\mathbf{r}_1, \dots, \mathbf{r}_N)/kT} d\mathbf{r}_{d+2} \dots d\mathbf{r}_N \quad (A3)$$

## 2. In dimensions one and two

In one-dimensional space, the volume of the void between the two 1D-spheres is $V_1(\mathbf{r}_1, \mathbf{r}_2) = V_1(r_{12}) = r_{12} - 2R_1$ if the mutual distance $r_{12} \ge 2R_1$ and is zero otherwise. By translation invariance, $H'(\mathbf{r}_1, \mathbf{r}_2) = H'(r_{12})$, and can be expanded around $2R_1$: $H'(\mathbf{r}_1, \mathbf{r}_2) = H'(2R_1) + \alpha V_1 + o(V_1)$. Equation A2 then gives $p_v(v_1) \propto V e^{H'(2R_1)}(1 + \alpha v_1 + o(v_1))$, where the total volume $V$ comes from the global translation invariance. The probability is then finished around $v_1 = 0$.

In two-dimensional space, the expression for $V_2(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$ is more involved, but depends only on the mutual distances: $V_2(r_{12}, r_{13}, r_{23})$. The first condition, that the 2D-spheres overlap by pairs, is $r_{ij} \le 2R_2$. Such a situation is depicted in Fig. 10, where the three molecules at $\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3$ (materialized by the points $M_1, M_2$ and $M_3$) are represented by solid circles of diameter $\sigma$, and their corresponding Stillinger 2D-spheres are represented by the dotted circles. The second condition, that the three spheres do not overlap together, is equivalent to the existence of the 3-branch star delimited by the points $M_1A M_2 B M_3 C$, where A, B and C are at the intersection between the three Stillinger disks that fall inside the triangle $M_1M_2M_3$ (see Fig. 10). The area of this star is given by the area $A_{123}$ of the main triangle $M_1M_2M_3$, minus the sum of the areas of the three minor triangles $M_1A M_2, M_2B M_3$ and $M_1C M_3$ in gray in Fig. 10. Using Heron's formula, $A_{123} = 0.25\sqrt{(r_{12}^2 + r_{13}^2 + r_{23}^2)^2 - 2(r_{12}^4 + r_{13}^4 + r_{23}^4)}$ and the sum of the area of the minor triangles is given by $\sum_{ij} 0.25 r_{ij} \sqrt{4R_2^2 - r_{ij}^2}$ for $ij = 12, 23$ and $13$. When these conditions are fulfilled, the void ABC (in black in Fig. 10) exists and its area is calculated as follows. Starting from the area of the full triangle $M_1M_2M_3$, one first removes the area of the three circular sectors defined by the intersection of the disks 1, 2 and 3 with the triangle $M_1M_2M_3$. This area equals $0.5(\theta_1 + \theta_2 + \theta_3)R_2^2 = 0.5\pi R_2^2$ (see Fig. 10). However, each hatched area $A_{ij}^h$ in Fig. 10 delimited by the intersection between the two disks i and j and the main triangle

$M_1M_2M_3$ is removed twice, and has to be compensated. Using $A_{ij}^h = R_2^2 \arccos(0.5r_{ij}/R_2) - 0.25r_{ij} \sqrt{4R_2^2 - r_{ij}^2}$, the final result reads:

$$V_2(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = A_{123} - 0.5\pi R_2^2 + A_{12}^h + A_{23}^h + A_{13}^h. \quad (A4)$$

This formula can be used to calculate numerically the exact volume of voids during the simulation of a 2D-fluid without resorting to a grid method.

The analytical calculation of $p_v(v_2)$ using Eqs A2 and A4 can be done by noticing that $V_2$ is invariant by translation and rotation of the three atoms as a whole. One then introduces the relative positions $\mathbf{r}'_i = \mathbf{r}_i - \mathbf{r}_1$ for $i = 2$ and 3. The integration over $d\mathbf{r}'_2$ is done in polar coordinates with trivial integration on the angular coordinate: $\int d\mathbf{r}'_2$ is replaced by $2\pi \int r_{12} dr_{12}$. The remaining integration on $d\mathbf{r}'_3$ is done as follows. When $r_{12}$ is chosen, the point A in Fig. 10 is determined. We introduce polar coordinates centered on A, the angles being measured with reference to the perpendicular bisector of points $M_1$ and $M_2$, denoted as $\Delta$ in Fig. 11. The polar coordinates of the third atom $M_3$ are denoted as $l_3 = |AM_3|$ and $\omega$ (see Fig. 11). The void $V_2$ exists only if $l_3 \ge R_2$. Let us denote $\xi = l_3 - R_2 \ge 0$. One has: $V_2(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = V_2(r_{12}, \xi, \omega)$. For a vanishing void, $\xi \to 0$, and the loci of the points $M_3$ compatible with this condition fall on the circular arc of radius $R_2$, centered on A, and limited by the two straight lines $M_1A$ and $M_2A$ extended on the side of $M_3$, corresponding to the constraint that the Stillinger disk of the third atom has to cut the two other ones (see the solid circular arc in Fig. 11). For a finite $V_2 = v_2$, the loci of $M_3$ points can be determined numerically by calculating $\xi^*(v_2, r_{12}, \omega)$ such that $V_2(r_{12}, \xi^*, \omega) = v_2$ (see for example the dash-dotted line in

<figure><img src="image_12.png" /><figcaption>FIG. 10. Schematic representation in the two-dimensional case of three atoms (solid circles of diameter $\sigma$) at $\mathbf{r}_1, \mathbf{r}_2$ and $\mathbf{r}_3$ materialized by the points $M_1, M_2$ and $M_3$, and their corresponding Stillinger 2D-spheres represented by the dotted circles. A small void (central curvilinear triangle ABC in black) may exist between these spheres if they overlap by pairs, but do not overlap all three together. The hatched areas resulting from the intersection between disks i and j and the main triangle $M_1M_2M_3$ are denoted as $A_{ij}^h$.</figcaption></figure>
Fig. 11 that gives the loci of points $M_3$ giving the void volume corresponding to the depicted atomic configuration). Gathering everything:

$$p_v(v_2) \propto \int e^{-H'(r_{12}, \xi, \omega)/kT} \delta(V_2(r_{12}, \xi, \omega) - v_2)(\xi + R_2)r_{12}dr_{12}d\omega d\xi \quad (A5)$$

where the integral runs from 0 to $2R_2$ for $r_{12}$, from 0 to $R_2$ for $\xi$ and from $\theta - \pi/2$ to $\pi/2 - \theta$ for $\omega$. The calculation can be done numerically for any void $v_2$, but an analytical expansion can be found in the vanishing limit $v_2 \to 0$. Expanding Eq A4, it can easily be shown that

$$V_2(r_{12}, \xi, \omega) = 0.5\xi^2\{\tan(\theta + \omega) + \tan(\theta - \omega)\} + o(\xi^2) \quad (A6)$$

where $\theta$ is defined in Fig. 11. Integration of the delta function with respect to $\xi$ gives a factor $(dV_2/d\xi)^{-1} = 0.5\xi^*(v_2, r_{12}, \omega)/v_2$ where $\xi^*$ as been previously defined. One gets:

$$p_v(v_2) \propto \int e^{-H'(r_{12}, \xi^*, \omega)/kT} \{0.5\xi^*(\xi^* + R_2)/v_2\} r_{12} dr_{12} d\omega \quad (A7)$$

Noticing that, from Eq A6, $\xi^*(v_2, r_{12}, \omega) = f(r_{12}, \omega)\sqrt{v_2} + o(\sqrt{v_2})$, and expanding $H'(r_{12}, \xi^*, \omega) = H'(r_{12}, 0, \omega) + \alpha(r_{12}, \omega)\xi^* + o(\xi^*)$, the probability can be written as:

$$p_v(v_2) \propto \frac{A_2}{\sqrt{v_2}} + B_2 + o(1) \quad (A8)$$

where

$$A_2 = 0.5 \int e^{-H'(r_{12}, 0, \omega)/kT} f(r_{12}, \omega) R_2 r_{12} dr_{12} d\omega \quad (A9)$$

<figure><img src="image_13.png" /><figcaption>FIG. 11. Same atomic configuration as in Fig. 10. $\Delta$ is the perpendicular bisector of points $M_1$ and $M_2$. $l_3$ and $\omega$ are the polar coordinates of $M_3$ with respect to point A and axis $\Delta$. The solid circular arc is the locus of $M_3$ giving a vanishing void $V_2$, while the dashed-dotted curve is the locus of $M_3$ giving the finite void volume $V_2$ corresponding to the depicted atomic configuration.</figcaption></figure>

and $B_2 =$

$$0.5 \int e^{-H'(r_{12}, 0, \omega)/kT} f(r_{12}, \omega)^2 (1 + \alpha(r_{12}, \omega)R_2) r_{12} dr_{12} d\omega \quad (A10)$$

Figure 12 shows the results of a two-dimensional simulation of disks in a plane. We have chosen to ignore the interactions (random disks). This is possible since their influence is limited to the constants $A_2$, $B_2$, etc. appearing in the expansion. The main behavior $\sim 1/\sqrt{v_2}$ is thus independent of the interactions. As can be seen in Fig. 12 the agreement is excellent in the vanishing limit (straight line in log-log scale).

### 3. In dimensions three and higher

In three dimensions, the expression of $V_3$ is even more involved, and an exact expression has not been found. However, as previously, invariance by translation and rotation can be used, as well as an expansion of the void volume in the vanishing limit. The calculation in dimension three essentially proceeds in a way similar to the 2D case, with the following differences:

(i) The conditions to have a closed void between the four Stillinger spheres require that they overlap three by three. In the previous 2D case, we encountered the opposite condition that three disks do not overlap, requiring that the area of the main triangle $M_1M_2M_3$ in Fig. 10 be larger than the sum of the areas of the three minor triangles $M_1AM_2$, $M_2BM_3$ and $M_1CM_3$ in gray in Fig. 10. This can easily be translated into a condition of overlapping of three spheres by considering their 2D intersection with the plane defined by their centers. Similarly, the condition of non-overlapping of the four spheres imposes a condition on their relative distances that will define the integration limits (next points).

(ii) As previously the integration is done using radial and angular coordinates. The previous integration on $r_{12}$ in 2D polar coordinates is replaced by a 3D spherical integration with the measure $r_{12}^2 dr_{12}$.

(iii) The position of the third atom is measured with respect to the first two ones, in cylindrical coordinates $(r_3, \phi_3, z_3)$

<figure><img src="image_14.png" /><figcaption>FIG. 12. Symbols: cavity size density $p_v$ in log-log scale for random disks in a two-dimensional box (green triangles) or spheres in a three-dimensional box (black circles). The green and black lines correspond to the main term in Eqs A8 and A14 for the 2D and 3D cases respectively.</figcaption></figure>
around the z axis defined by the two first atoms. $r_3$ is the distance of atom 3 with respect to this axis. The corresponding integration measure is $r_3dr_3d\phi_3dz_3$. By rotation invariance, integration over $\phi_3$ is trivial.

(iv) The position of the fourth atom is measured in spherical coordinates with the origin at the intersection of the first three spheres. As in the 2D case, a void between the four spheres exists only if the radial coordinate $r_4$ is larger than $R_3$, the Stillinger radius in dimension three. We thus introduce $\xi = r_4 - R_3 \ge 0$. The orientation is measured with two angles denoted $\omega_4$ and $\phi_4$, and the integration measure is $(\xi + R_3)^2 d\xi \sin \omega_4 d\omega_4 d\phi_4$.

The probability can therefore be written as follows:

$$p_v(v_3) = \int e^{-H'(\xi, \Omega)/kT} \delta(V_3(\xi, \Omega) - v_3)(\xi + R_3)^2 d\xi d\Omega \quad (A11)$$

where $\Omega$ is a shortcut for $(r_{12}, r_3, z_3, \omega_4, \phi_4)$ and $d\Omega = \sin \omega_4 dr_{12} dr_3 dz_3 d\omega_4 d\phi_4$.

Simple geometric considerations give the expansion of $V_3$ in the vanishing limit:

$$V_3(\xi, \Omega) = \xi^3 g(\Omega) \quad (A12)$$

We introduce $\xi^*(v_3, \Omega)$ such that $V_3(\xi^*, \Omega) = v_3$. The delta function introduces the factor $(dV_3/d\xi)^{-1} = (1/3)\xi^*/v_3$, and one finally gets:

$$p_v(v_3) \propto \int e^{-H'(\xi^*, \Omega)/kT} \{(1/3)\xi^*(\xi^* + R_3)^2/v_3\} d\Omega \quad (A13)$$

As previously, expanding $H'(\xi, \Omega)$ around $\xi = 0$, one finally gets

$$p_v(v_3) \propto \frac{A_3}{v_3^{2/3}} + \frac{B_3}{v_3^{1/3}} + C_3 + o(1) \quad (A14)$$

where $A_3$, $B_3$ and $C_3$ follow from equations similar to Eqs A9 and A10.

Figure 12 shows the result of a three-dimensional simulation of spheres in space. As previously, the interactions are ignored since their influence is limited to the constants $A_3$, $B_3$, etc. As can be seen, the distribution in the vanishing limit follows the $v_3^{-2/3}$ law (straight line in log-log scale). It is worth emphasizing that this divergent prefactor still results in a valid probability distribution, in the sense that its definite integral from 0 to infinity converges.

The previous argument suggests a generalization to dimension $d$. In particular, we use the symmetries as well as the radial and angular coordinates. Dimensional argument is used to obtain the cavity expansion in the vanishing limit resulting in an equation similar to Eq A12. After integration, one gets:

$$p_v(v_d) \propto \frac{A_d}{v_d^{1-1/d}} + \frac{B_d}{v_d^{1-2/d}} + \frac{C_d}{v_d^{1-3/d}} + \dots \quad (A15)$$

where the details of the interatomic forces are absorbed into the constants $A_d$, $B_d$, etc. As can be seen, except in one-dimension, the probability diverges for vanishing voids.

In the context of the classical nucleation theory, the bubbles are supposed to be spherical, and the probability distribution of bubble sizes is generally given as a function of the sphere radius $r$. In this case, the bubble size density should verify $p_r(r)dr = p_v(v_d)dv_d$ (= the number of bubbles between $r$ and $r+dr$). Since $v_d \propto r^d$, one has:

$$p_r(r) = d^{\frac{V_d}{r}} p_v(v_d) \propto A'_d + B'_d r + C'_d r^2 + \dots \quad (A16)$$

In this radial representation, the bubble size density is expected to be finite for vanishing bubbles. This is an important point that allows one to define the free energy of formation of a bubble (see main text).

<footnote n="1">¹J. G. Kirkwood and F. P. Buff, “The statistical mechanical theory of solutions. I,” The Journal of Chemical Physics **19**, 774–777 (1951).</footnote>

<footnote n="2">²H. Reiss, H. L. Frisch, and J. L. Lebowitz, “Statistical mechanics of rigid spheres,” The Journal of Chemical Physics **31**, 369–380 (1959).</footnote>

<footnote n="3">³E. Helfand, H. Reiss, H. L. Frisch, and J. L. Lebowitz, “Scaled particle theory of fluids,” The Journal of Chemical Physics **33**, 1379–1385 (1960).</footnote>

<footnote n="4">⁴D. S. Corti, P. G. Debenedetti, S. Sastry, and F. H. Stillinger, “Constraints, metastability, and inherent structures in liquids,” Physical Review E **55**, 5522–5534 (1997).</footnote>

<footnote n="5">⁵P. J. in 't Veld, M. T. Stone, T. M. Truskett, and I. C. Sanchez, “Liquid structure via cavity size distributions,” The Journal of Physical Chemistry B **104**, 12028–12034 (2000).</footnote>

<footnote n="6">⁶J. M. Simon, P. Krüger, S. K. Schnell, T. J. H. Vlugt, S. Kjelstrup, and D. Bedeaux, “Kirkwood–Buff integrals: From fluctuations in finite volumes to the thermodynamic limit,” The Journal of Chemical Physics **157**, 130901 (2022).</footnote>

<footnote n="7">⁷N. Giri, M. G. Del Pópolo, G. Melaugh, R. L. Greenaway, K. Rätzke, T. Koschine, L. Pison, M. F. C. Gomes, A. I. Cooper, and S. L. James, “Liquids with permanent porosity,” Nature **527**, 216–220 (2015).</footnote>

<footnote n="8">⁸H. Reiss, H. L. Frisch, E. Helfand, and J. L. Lebowitz, “Aspects of the statistical thermodynamics of real fluids,” The Journal of Chemical Physics **32**, 119–124 (1960).</footnote>

<footnote n="9">⁹J. P. M. Postma, H. J. C. Berendsen, and J. R. Haak, “Thermodynamics of cavity formation in water. A molecular dynamics study,” Faraday Symposia of the Chemical Society **17**, 55–67 (1982).</footnote>

<footnote n="10">¹⁰A. Pohorille and L. R. Pratt, “Cavities in molecular liquids and the theory of hydrophobic solubilities,” Journal of the American Chemical Society **112**, 5066–5074 (1990).</footnote>

<footnote n="11">¹¹J. L. Prelesnik, J. L. Chen, K. Mahesh, and J. I. Siepmann, “Anomalous pressure–density relations and speed of sound in bubbly water systems,” The Journal of Chemical Physics **161**, 204501 (2024).</footnote>

<footnote n="12">¹²J. L. Chen, B. Xue, K. Mahesh, and J. I. Siepmann, “Molecular simulations probing the thermophysical properties of homogeneously stretched and bubbly water systems,” Journal of Chemical & Engineering Data **64**, 3755–3771 (2019).</footnote>

<footnote n="13">¹³A. Arvengas, E. Herbert, S. Cersoy, K. Davitt, and F. Caupin, “Cavitation in heavy water and other liquids,” The Journal of Physical Chemistry B **115**, 14240–14245 (2011).</footnote>

<footnote n="14">¹⁴V. G. Baidakov and K. S. Bobrov, “Spontaneous cavitation in a Lennard-Jones liquid at negative pressures,” The Journal of Chemical Physics **140**, 184506 (2014).</footnote>

<footnote n="15">¹⁵V. G. Baidakov and K. R. Protsenko, “Molecular dynamics simulation of cavitation in a Lennard-Jones liquid at negative pressures,” Chemical Physics Letters **760**, 138030 (2020).</footnote>

<footnote n="16">¹⁶V. Doebele, A. Benoit-Gonin, F. Souris, L. Cagnon, P. Spathis, P.-E. Wolf, A. Grosman, M. Bossert, I. Trimaille, and E. Rolley, “Direct observation of homogeneous cavitation in nanopores,” Physical Review Letters **125**, 255701 (2020).</footnote>

<footnote n="17">¹⁷J. Puibasset, “Cavitation in heterogeneous nanopores: The chemical ink-bottle,” AIP Advances **11**, 095311 (2021).</footnote>
¹⁸M. Bossert, A. Grosman, I. Trimaille, F. Souris, V. Doebele, A. Benoit-Gonin, L. Cagnon, P. Spathis, P.-E. Wolf, and E. Rolley, "Evaporation process in porous silicon: Cavitation vs pore blocking," Langmuir **37**, 14419–14428 (2021).

¹⁹M. Bossert, I. Trimaille, L. Cagnon, B. Chabaud, C. Gueneau, P. Spathis, P. E. Wolf, and E. Rolley, "Surface tension of cavitation bubbles," Proceedings of the National Academy of Sciences of the United States of America **120**, e2300499120 (2023).

²⁰K. M. Bal and E. C. Neyts, "Extending and validating bubble nucleation rate predictions in a Lennard-Jones fluid with enhanced sampling methods and transition state theory," The Journal of Chemical Physics **157**, 184113 (2022).

²¹K. R. Protsenko and V. G. Baidakov, "Classical nucleation theory and molecular dynamics simulation: Cavitation," Physics of Fluids **35**, 014111 (2023).

²²C. P. Lamas, E. Sanz, C. Vega, and E. G. Noya, "Estimation of bubble cavitation rates in a symmetrical Lennard-Jones mixture by NVT seeding simulations," The Journal of Chemical Physics **158**, 124109 (2023).

²³J. C. Fisher, "The fracture of liquids," Journal of Applied Physics **19**, 1062–1067 (1948).

²⁴M. Blander and J. L. Katz, "Bubble nucleation in liquids," American Institute of Chemical Engineers Journal **21**, 833–848 (1975).

²⁵P. G. Debenedetti, *Metastable Liquids: Concepts and Principles* (Princeton University Press, Princeton, NJ, 1996).

²⁶J. W. Gibbs, *The Scientific Papers of J Willard Gibbs* (Longmans Green, London, 1906).

²⁷H. Reiss and R. K. Bowles, "Some fundamental statistical mechanical relations concerning physical clusters of interest to nucleation theory," The Journal of Chemical Physics **111**, 7501–7504 (1999).

²⁸T. Kinjo and M. Matsumoto, "Cavitation processes and negative pressure," Fluid Phase Equilibria **144**, 343–350 (1998).

²⁹G. Menzl, M. A. Gonzalez, P. Geiger, F. Caupin, J. L. F. Abascal, C. Valeriani, and C. Dellago, "Molecular mechanism for cavitation in water under tension," Proceedings of the National Academy of Sciences USA **113**, 13582 (2016).

³⁰V. K. Shen and P. G. Debenedetti, "A computational study of homogeneous liquid-vapor nucleation in the Lennard-Jones fluid," The Journal of Chemical Physics **111**, 3581–3589 (1999).

³¹A. Vishnyakov, P. G. Debenedetti, and A. V. Neimark, "Statistical geometry of cavities in a metastable confined fluid," Physical Review E **62**, 538–544 (2000).

³²Y. W. Wu and C. Pan, "A molecular dynamics simulation of bubble nucleation in homogeneous liquid under heating with constant mean negative pressure," Microscale Thermophysical Engineering **7**, 137–151 (2003).

³³A. V. Neimark and A. Vishnyakov, "The birth of a bubble: a molecular simulation study," The Journal of Chemical Physics **122**, 054707 (2005).

³⁴J. Wedekind, G. Chkonia, J. Wölk, R. Strey, and D. Reguera, "Crossover from nucleation to spinodal decomposition in a condensing vapor," The Journal of Chemical Physics **131**, 114506 (2009).

³⁵H. Watanabe, M. Suzuki, and N. Ito, "Cumulative distribution functions associated with bubble-nucleation processes in cavitation," Physical Review E **82**, 051604 (2010).

³⁶S. L. Meadley and F. A. Escobedo, "Thermodynamics and kinetics of bubble nucleation: Simulation methodology," The Journal of Chemical Physics **137**, 074109 (2012).

³⁷J. Frenkel, "A general theory of heterophase fluctuations and pretransition phenomena," The Journal of Chemical Physics **7**, 538–547 (1939).

³⁸J. Frenkel, "Statistical theory of condensation phenomena," The Journal of Chemical Physics **7**, 200–201 (1939).

³⁹K. J. Oh and X. C. Zeng, "Formation free energy of clusters in vapor-liquid nucleation: A Monte Carlo simulation study," The Journal of Chemical Physics **110**, 4471–4476 (1999).

⁴⁰P. R. ten Wolde and D. Frenkel, "Computer simulation study of gas-liquid nucleation in a Lennard-Jones system," The Journal of Chemical Physics **109**, 9901–9918 (1998).

⁴¹I. Kusaka, Z. G. Wang, and J. H. Seinfeld, "Direct evaluation of the equilibrium distribution of physical clusters by a grand canonical Monte Carlo simulation," The Journal of Chemical Physics **108**, 3416–3423 (1998).

⁴²K. Yasuoka and M. Matsumoto, "Molecular dynamics of homogeneous nucleation in the vapor phase. I. Lennard-Jones fluid," The Journal of Chemical Physics **109**, 8451–8462 (1998).

⁴³K. Yasuoka and M. Matsumoto, "Molecular dynamics of homogeneous nucleation in the vapor phase. II. Water," The Journal of Chemical Physics **109**, 8463–8470 (1998).

⁴⁴S. Auer and D. Frenkel, "Numerical prediction of absolute crystallization rates in hard-sphere colloids," The Journal of Chemical Physics **120**, 3015–3029 (2004).

⁴⁵H. M. Ellerby, C. L. Weakliem, and H. Reiss, "Toward a molecular theory of vapor-phase nucleation. I. Identification of the average embryo," The Journal of Chemical Physics **95**, 9209–9218 (1991).

⁴⁶H. M. Ellerby and H. Reiss, "Toward a molecular theory of vapor-phase nucleation. II. Fundamental treatment of the cluster distribution," The Journal of Chemical Physics **97**, 5766–5772 (1992).

⁴⁷P.-R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, "Simulation of homogeneous crystal nucleation close to coexistence," Faraday Discussions **104**, 93–110 (1996).

⁴⁸A. C. Pan and D. Chandler, "Dynamics of nucleation in the Ising model," J. Phys. Chem. B **108**, 19681–19686 (2004).

⁴⁹I. Saika-Voivod, P. H. Poole, and R. K. Bowles, "Test of classical nucleation theory on deeply supercooled high-pressure simulated silica," The Journal of Chemical Physics **124**, 224709 (2006).

⁵⁰S. E. M. Lundrigan and I. Saika-Voivod, "Test of classical nucleation theory and mean first-passage time formalism on crystallization in the Lennard-Jones liquid," The Journal of Chemical Physics **131**, 104503 (2009).

⁵¹I. Saika-Voivod, R. K. Bowles, and P. H. Poole, "Crystal nucleation in a supercooled liquid with glassy dynamics," Physical Review Letters **103**, 225701 (2009).

⁵²Y. Goswami, V. V. Vasisht, D. Frenkel, P. G. Debenedetti, and S. Sasy, "Thermodynamics and kinetics of crystallization in deeply supercooled Stillinger-Weber silicon," J. Chem. Phys. **155**, 194502 (2021).

⁵³M. A. González, G. Menzl, J. L. Aragones, P. Geiger, F. Caupin, J. L. F. Abascal, C. Dellago, and C. Valeriani, "Detecting vapour bubbles in simulations of metastable water," The Journal of Chemical Physics **141**, 18C511 (2014).

⁵⁴Z.-J. Wang, C. Valeriani, and D. Frenkel, "Homogeneous bubble nucleation driven by local hot spots: A molecular dynamics study," The Journal of Physical Chemistry B **113**, 3776–3784 (2009).

⁵⁵M. P. Allen and D. J. Tildesley, *Computer Simulation of Liquids* (Clarendon Press, Oxford, 1987).

⁵⁶D. Frenkel and B. Smit, *Understanding Molecular Simulation* (Academic Press, London, 2002).

⁵⁷H. Watanabe, N. Ito, and C.-K. Hu, "Phase diagram and universality of the Lennard-Jones gas-liquid system," The Journal of Chemical Physics **136**, 204102 (2012).

⁵⁸J. Stillinger, Frank H., "Rigorous basis of the Frenkel-Band theory of association equilibrium," The Journal of Chemical Physics **38**, 1486–1494 (1963).

⁵⁹G. M. Torrie and J. P. Valleau, "Monte Carlo free energy estimates using non-Boltzmann sampling: Application to the sub-critical Lennard-Jones fluid," Chem. Phys. Lett. **28**, 578–581 (1974).

⁶⁰L. Maibaum, "Comment on "elucidating the mechanism of nucleation near the gas-liquid spinodal"," Physical Review Letters **101**, 019601 (2008).

⁶¹S. Chakrabarty, M. Santra, and B. Bagchi, "Chakrabarty, Santra, and Bagchi Reply," Physical Review Letters **101**, 019602 (2008).

⁶²J.-M. Leyssale, J. Delhommelle, and C. Millot, "Atomistic simulation of the homogeneous nucleation and of the growth of N₂ crystallites," The Journal of Chemical Physics **122**, 104510 (2005).

⁶³P. R. ten Wolde and D. Frenkel, "Enhancement of protein crystal nucleation by critical density fluctuations," Science **277**, 1975–1978 (1997).

⁶⁴J. Puibasset, "A general relation between the largest nucleus and all nuclei distributions for free energy calculations," The Journal of Chemical Physics **157**, 191102 (2022).

⁶⁵P. Porion and J. Puibasset, "A statistical analysis of the first stages of freezing and melting of Lennard-Jones particles: Number and size distributions of transient nuclei," The Journal of Chemical Physics **161**, 074501 (2024).

⁶⁶A. Stukowski, "Visualization and analysis of atomistic simulation data with OVITO—the Open Visualization Tool," Modelling and Simulation in Materials Science and Engineering **18**, 015012 (2010).

⁶⁷A. K. Sharma and F. A. Escobedo, "Nucleus-size pinning for determination of nucleation free-energy barriers and nucleus geometry," The Journal of
Chemical Physics **148**, 184104 (2018).

son, "Communications: Evidence for the role of fluctuations in the thermo-dynamics of nanoscale drops and the implications in computations of the surface tension," The Journal of Chemical Physics **132**, 141101 (2010).

<footnote n="68"><sup>68</sup>S. Prestipino, A. Laio, and E. Tosatti, "Systematic improvement of classical nucleation theory," Physical Review Letters **108**, 225701 (2012).</footnote>

<footnote n="69"><sup>69</sup>A. Aasen, Ø. Wilhelmsen, M. Hammer, and D. Reguera, "Free energy of critical droplets—from the binodal to the spinodal," The Journal of Chemical Physics **158**, 114108 (2023).</footnote>

<footnote n="70"><sup>70</sup>W. Gispen, J. R. Espinosa, E. Sanz, C. Vega, and M. Dijkstra, "Variational umbrella seeding for calculating nucleation barriers," The Journal of Chemical Physics **160**, 174501 (2024).</footnote>

<footnote n="71"><sup>71</sup>R. C. Tolman, "The effect of droplet size on surface tension," The Journal of Chemical Physics **17**, 333–337 (1949).</footnote>

<footnote n="72"><sup>72</sup>J. G. Sampayo, A. Malijevský, E. A. Müller, E. de Miguel, and G. Jack-</footnote>

<footnote n="73"><sup>73</sup>P. Rehner and J. Gross, "Surface tension of droplets and tolman lengths of real substances and mixtures from density functional theory," The Journal of Chemical Physics **148**, 164703 (2018).</footnote>

<footnote n="74"><sup>74</sup>J. Julin, I. Napari, J. Merikanto, and H. Vehkamäki, "A thermodynamically consistent determination of surface tension of small Lennard-Jones clusters from simulation and theory," The Journal of Chemical Physics **133**, 044704 (2010).</footnote>