What experiments on pinned nanobubbles can tell about the critical nucleus for bubble nucleation

The process of homogeneous bubble nucleation is almost impossible to probe experimentally, except near the critical point or for liquids under large negative tension. Elsewhere in the phase diagram, the bubble nucleation barrier is so high as to be effectively insurmountable. Consequently, there is a severe lack of experimental studies of homogenous bubble nucleation under conditions of practical importance (e.g., cavitation). Here we use a simple geometric relation to show that we can obtain information about the homogeneous nucleation process from Molecular Dynamics studies of bubble formation in solvophobic nanopores on a solid surface. The free energy of pinned nanobubbles has two extrema as a function of volume: one state corresponds to a free-energy maximum (“the critical nucleus”), the other corresponds to a free-energy minimum (the metastable, pinned nanobubble). Provided that the surface tension does not depend on nanobubble curvature, the radius of the curvature of the metastable surface nanobubble is independent of the radius of the pore and is equal to the radius of the critical nucleus in homogenous bubble nucleation. This observation opens the way to probe the parameters that determine homogeneous bubble nucleation under experimentally accessible conditions, e.g. with AFM studies of metastable nanobubbles. Our theoretical analysis also indicates that a surface with pores of different sizes can be used to determine the curvature corrections to the surface tension. Our conclusions are not limited to bubble nucleation but suggest that a similar approach could be used to probe the structure of critical nuclei in crystal nucleation.

Bubble nucleation is a very common process but the observation of truly homogenous bubble nucleation is very rare [1]. The reason is that the free-energy barrier for homogenous bubble nucleation tends to be very high (order 10 3 -10 4 eV), except near the critical point (typical T > 0.9T c ) where the surface tension is very low, or under conditions of extreme under-pressure [2]. For this reason, bubble nucleation under less extreme circumstances (e.g., bubble formation in boiling water) is always dominated by heterogeneous nucleation. The rate of heterogeneous nucleation is normally extremely sensitive to the properties of the surface, in particular to the value of the solid-liquid and solid-vapor surface free energies. Hence, normally, heterogeneous bubble nucleation does not provide direct information about the homogenous nucleation process.
Supplementary material in the form of a .pdf file available from the Journal web page at http://dx.doi.org/10.1140/epje/i2017-11604-7 a e-mail: df246@cam.ac.uk b e-mail: jd489@cam.ac.uk c e-mail: zhangxr@mail.buct.edu.cn The situation is very different in the case of welldefined hydrophobic (or, more generally, "solvophobic") nanopores in a hydrophilic ("solvophilic") surface. Inside these pores, vapor may nucleate easily, sometimes even at a pressure where the bulk liquid is still thermodynamically stable. As the pressure of the liquid is lowered, the vapor phase will bulge out from the entrance of the nanopores. However, if the remainder of the surface is solvophilic, these incipient surface nanobubbles [3][4][5][6][7][8] are pinned at the perimeter of the nanopore [9][10][11]. The radius of curvature of the metastable nanobubble is, as is shown here, the same as (or very close to) the one of the critical nucleus in homogeneous nucleation. This quantity is usually not accessible in bulk experiments.
Below, we show that if the classical nucleation theory (CNT) relation holds, the radii of pinned metastable nanobubbles and critical nuclei in the bulk are the same. The relation is not exact if the surface tension is curvature dependent [12,13]. However, for the geometry that we study, the lowest-order curvature effects vanish. We then present molecular dynamics (MD) simulations to show that metastable nanobubbles do indeed have the same curvature radius of curvature as the critical nuclei. The simulations also show that, as CNT would predict, the pressure difference is independent of the radius of the pore that pins the nanobubble, and is given by the same Young-Laplace relation that would apply to critical nuclei in homogenous bubble nucleation. Our findings therefore indicate that direct probe experiments on pinned bubbles on carefully patterned substrates should enable us to determine the elusive properties of the critical nucleus for homogenous bubble nucleation. Similarly, we can use simulations of pinned nanobubbles to study the properties of critical bubbles, even when using system sizes that are orders of magnitude smaller than the size of the bulk critical bubble.
Theoretical background. Within the classical nucleation theory (CNT) [1,[14][15][16], the free-energy cost for the formation of a bubble is With Δp the pressure difference across the interface, γ the vapor-liquid interface tension, V the volume of the bubble and A the area of the liquid-vapor interface. Note that for highly compressible bubbles, this expression is only meaningful if the bubble is in chemical equilibrium with the surrounding liquid. Critical nuclei are in chemical and mechanical equilibrium with the surrounding fluid while non-critical nuclei (or unstable nanobubbles) cannot be in mechanical and chemical equilibrium except in small systems with a fixed volume and number of particles. In such cases, the more modern approaches to compute the free energy [17][18][19] can be used. However, in the present work we focus on bubbles that are both in mechanical and chemical equilibrium with the surrounding fluid at constant pressure, in which case eq. (1) holds. In case of homogeneous nucleation, the change in free energy due to the formation of a spherical bubble of radius R is given by ΔG = −Δp 4πR 3 3 +4πR 2 γ and the critical radius R * is the radius at which the bubble is in equilibrium, i.e. ∂ΔG/∂R = 0, yielding the Laplace-Young equation The expression in eq. (2) depends on the geometry, e.g. for bubbles forming in a slitpore, the numerical factor on the right-hand side is 1 instead of 2. We express the excess free energy of a nanobubble pinned on a solid pore as the sum of two terms, a volume term Δp ΔV , where ΔV is the excess volume (i.e. the volume outside the nanopore) and a surface term γΔA, where ΔA is the difference between the area of the nanobubble and the cross-section of the pore As we consider the case where the contact line is pinned at the rim of the pore, the magnitudes of the solid-liquid and solid-vapor interfaces do not depend on the bubble size, and hence can be ignored in eq. (3). To estimate the height and the radius of curvature of a metastable bubble, we first write down the expression for the variation of the free energy: dΔG NB = −ΔpdΔV + d(γΔA) = −ΔpdΔV +γ 0 dΔA+d(ΔAγ ), where γ 0 is the part of surface tension that does not depend on the curvature and γ = γ − γ 0 . The nanobubble radius R is obtained by requiring ∂ΔGNB ∂R = 0. We will study pores with circular cross-section and extended slit pores. As derived in the Supplementary Material (SM), in both cases, the excess volume and excess area are related through where D is one for a straight pore and two for a circular pore. Taking into account this relation, the bubble radius is given by the following expression: The first observation from eq. (5) is that when the surface tension does not depend on the curvature, i.e. γ = 0, the radius of the nanobubble is exactly equal to the radius of the critical nucleus, R = R * . This is true both for circular and straight pores. We performed a stability analysis for the pinned nanobubbles. From ∂ΔGNB ∂R = 0, as shown in the SM, we demonstrate that there exist two types of nanobubbles pinned on a circle pore, one with a contact angle determined by sin θ = r/R (and thus θ < π/2) and the other having a contact angle of π − θ (see fig. 1(b)). Both have the same curvature radius that is equal to that of the critical nucleus for homogeneous nucleation ( fig. 1(a)). The stability of these nanobubbles, as also shown in the SM, can be inferred from the equation, 2 cos θ . Hence, for θ > π/2, cos θ < 0 and ∂ 2 ΔG ∂R 2 | R=RNB < 0. In this case, a local maximum of free energy indicates that the pinned bubbles correspond to the critical nuclei for pinned nanobubbles losing their stability, leading to a liquid-to-vapor phase transition under pinning condition (see fig. 1(b)). This kind of pinned nanobubbles, which are also in equilibrium with the environment, behave similarly as the critical bubble for homogeneous liquid-to-vapor nucleation: they are of the same mean curvature and both unstable. This differs from the critical nucleus for heterogeneous nucleation because its contact angle is determined by sin θ = r/R, rather than by the substrate chemistry. But for θ < π/2, cos θ > 0, and thus ∂ 2 G ∂R 2 | R=RNB > 0 the pinned nanobubbles correspond to a local free-energy minimum, and thus are thermodynamically (meta)stable (see fig. 1(b)). This indicates that the accessible range of the contact angles for stable nanobubbles should always be smaller than π/2.
For sufficiently small nanobubbles, the dependence of the surface tension on the curvature of the bubble cannot be ignored [12,13]. In that case, R depends on the pore size and is generally not the same as the radius of the critical nucleus. This observation is interesting, because it suggests that experiments could detect the presence of curvature effects in a simple way. To demonstrate this, we note that in eq. (5) we can determine the derivative by measuring the left-hand side of this equation for a range of R values with varying pore sizes. By integrating the expression, we can evaluate γ . In the simple case the surface tension can be written as with δ T the Tolman length and κ the surface (Helfrich) bending stiffness. By performing measurements on circular and straight pores with various sizes, it is possible to determine both curvature corrections δ T and κ. Simulation validation. The above analysis is based on CNT, which is a "macroscopic" theory. To validate our main result, i.e. that the radius of curvature of a pinned nanobubble in a solid pore is equal to the radius of a critical nucleus in homogeneous nucleation process, we performed MD simulations of bubble formation in a Lennard-Jones fluid in contact with a slit-like pore, using an open-source program (LAMMPS [20]). We employed a pseudo-2D simulation box with a slab-like geometry of 65.8 × 13.2 × H in units σ 3 , where σ is the Lennard-Jones diameter and H, the height of simulation box, fluctuates during the simulation process to fix the pressure ( fig. 2(a)). Periodic boundary conditions were used in x and y directions. This geometry is convenient, to generate "infinite" (i.e. periodically repeated), cylindrical nanobubbles. In principle, it is also possible to generate spherical bubbles on circular pores. However, the present approach yields better statistics, and therefore more accurate answers.
At the top and bottom of the simulation box we placed two substrates that consist of frozen solid atoms on the (100) surface of FCC lattice with a lattice size of 1.64 σ. On the bottom substrate of the pseudo-2D system, we placed a rectangular nanopore of width r and depth 16.4 σ, that is periodically repeated along the short dimension (13.2 σ) of the simulation box. This slit pore exerts a pinning force on the contact line of nanobubbles. For simplicity, we prepared pores that already contained nanobubbles -hence there was no need to make the inside walls of the pores solvophobic. The height of the simulation box could change, subject to an external force along the z-direction to maintain the imposed pressure.
For the intermolecular potentials, we used the standard Lennard-Jones (LJ) 12-6 potential, with well-depth ε. The potential was truncated and shifted at a cutoff distance of 3.2σ. In addition, the interaction between bottom solid and fluid atoms was adopted to a value of 0.55ε, generating a medium wettability so that a pure LJ droplet would have a contact angle about 89 • , as that of HOPG used experimentally [10]. We carried out isothermal-isostress (NP ZZ T ) simulations with a fixed number of fluid molecules, N = 31294. The equations of motion were integrated using the velocity-Verlet algorithm with a time step of 0.0023 in the unit of mσ 2 /ε. The fluid temperature was controlled at a reduced temperature k B T /ε = 0.96 using a Nosé-Hoover thermostat with a time constant 0.23 in reduced units. In our simulations, surface nanobubbles were produced as follow. First, we filled N fluid molecules into the region enclosed by the solid substrates, without a pore, followed by an equilibration run of 6 × 10 6 MD time steps, at a temperature k B T /ε = 0.96 and the corresponding coexistence pressure p ext = 0.029 (in units εσ −3 ). After equilibration, we reduced the pressure to the desired target value. Next, we introduced slit pores of varying width. For each pore size, we performed a simulation of 14×10 6 time steps, to produce stable surface nanobubbles (a typical sequence of simulation snapshots is shown in fig. 2(c)).
To obtain the radius of curvature of a stable nanobubble at a constant external pressure, time evolution of the height of simulation box at different pore radii, r, were recorded (see fig. 2(b)). Note that for pseudo-2D pores, the pore radius is equal to half the pore width. The box height increases with increasing the pore radius. For the pore radii between 7.5 and 15.0 σ, stable nanobubbles were observedand the box height H fluctuates strongly around the mean value (see fig. 2(b)). For r < 7.5 σ (e.g., r = 5.0 σ in fig. 2(b)) the height reaches an equilibrium value but fluctuations are strongly suppressed. The corresponding snapshots show that at this pore radius, the nanobubble disappeared and a Cassie state with the planar vaporliquid interface on the pore wall was formed. For pore radii larger than 16.5 σ, H increases rapidly, indicating unbounded expansion of the nanobubble. In this case, the nanobubble initially grew rather slowly until the contact line was de-pinned at about 12 × 10 6 MD steps followed by a rapid expansion of the simulation box ( fig. 2(b)).
The vapor-liquid interfaces for stable nanobubbles were determined through finding the locations at which the fluid density is equal to half of the liquid density. Note that the fluid density profile was obtained via averaging over 1000 configurations within 10 6 time steps. As expected [21][22][23][24] the vapor-liquid interface of the pinned nanobubble can be well described by a cylinder segment. This observation facilitates the calculation of the nanobubble curvature radius R and its contact angle θ, which are depicted in fig. 3. We note that for nanobubbles on various pore widths r, the radius of curvature R is almost constant, while the contact angle θ measured from the liquid side decreases with the increase of pore radius. This indicates that at a given external pressure, stable surface nanobubbles have an equilibrium curvature radius independent of the pore size and as the pore grows larger, the nanobubble contact angle decreases.
Next, we show that the curvature radius of nanobubbles observed in MD simulations is indeed equal to that of the critical nucleus of the homogeneous liquid-tovapor transition at the same thermodynamic conditions. Nanobubbles or nuclei in the bulk solution are thermodynamically unstable: For bubbles with radii larger than a critical value, the liquid-to-vapor phase transition occurs, while bubbles smaller than the critical size tend to shrink and disappear. This is the reason why observing the critical nucleus in experiments is challenging.
Here, we first determined the critical nucleus for the vapor-liquid transition by means of MD simulations. To obtain the size of critical nucleus for the homogeneous nucleation, independent MD simulation runs were performed with predefined bulk nanobubbles in liquid phase. The predefined nanobubbles were produced as follows. First, N = 34531 liquid molecules were randomly arranged in a simulation box with dimensions of 63.6 × 13.2 × 63.6σ 3 and periodic boundary conditions were applied in all three directions. Note that in order to directly compare to nanobubbles obtained above, we used again pseudo-2D simulation box here and we studied periodically continued cylindrical bubbles, rather than spherical bubbles. A simulation of 2 × 10 6 MD steps was carried out at the same reduced temperature as before. To produce candidate critical bubble nuclei in the bulk, we created a cylindrical cavity in a bulk liquid by removing molecules inside the cylindrical region. A second MD run of 5 × 10 5 MD steps was performed to equilibrate the nucleus. Finally, another 5 × 10 5 MD steps NP ZZ T -MD run was carried out to obtain the time evolution of the nanobubble radius in the bulk.
The time evolution of the nucleus size at different initial radii of nuclei is depicted in fig. 4. We observe that the nuclei either expand into a vapor phase or shrink and return to the metastable liquid phase. We can determine Fig. 4. Time evolution of the bubble radius at the external pressure of p ext = 0.012. We expect the spontaneous nucleation to take place for bubbles with radius R = 18.5 σ. Of course, nucleation is a stochastic event (e.g., some post-critical nuclei may still disappear), and therefore our estimate of the critical nucleus size from the observed nucleation events is subject to an appreciable statistical error, which we estimate to be at least 1.5 σ.
the size of the critical nucleus: it is within a range of 16.4-22 σ at a reduced external pressure of 0.012.
We further note that under the same thermodynamic conditions (temperature and pressure), the critical nuclei have the same radius as the metastable surface nanobubbles (see fig. 3). Furthermore, when the radius of the nanopore is close to or larger than the size of the critical nucleus, the surface nanobubble becomes unstable ( fig. 2(b)). This observation confirms that anchored nanobubbles can mimic the critical nucleus, and this finding offers a new way to study its properties of experimentally.
We also determined nanobubble internal pressure and surface tension. According to the Young-Laplace equation for the cylindrical case, p in − p out = γ lv R with p in the internal pressure of nanobubble and p out the pressure of the bulk phase. In our simulations, we calculated the internal pressure of nanobubbles with different radii. This allowed us to determine the surface tension for highly curved vapor-liquid interface. Using these calculations, we can determine the pressure difference across the nanobubble interface and show that the surface nanobubbles behave analogous to nanobubbles in the bulk solution (the critical nuclei for homogeneous nucleation).
To calculate the nanobubble internal pressure at different radii, we varied the external pressure and evaluated the density profile from the MD simulations. From the density profile averaged over 2 ns of equilibrium simulations, the local density of the portion that is far from interface and of uniform particle distribution was chosen for determining the internal pressure via MBWR equation of state proposed by Johnson et al. [25]. Because we used the isothermal, isostress (NP ZZ T ) ensemble, the liquid pressure far from the liquid-vapor interface should be equal to the imposed pressure. We confirmed this by performing independent MD simulations with a planar liquid film to calculate the pressure in z-direction at fixed temperature and density of the liquid.
Having obtained the pressure inside and outside the nanobubble, we can compute the pressure difference. The results are shown in fig. 5(a). From the figure, we can see that the internal pressure of the nanobubble is nearly independent of the bubble radius at least for the range of nanobubbles we studied (10 σ-45 σ). Such behavior is to be expected, because changing the external pressure from −0.0024 to 0.022 changes the chemical potential of the fluid molecules by less than 4%, and hence, we should expect the (almost ideal) vapor density inside the bubble to vary by about the same amount. Clearly, the nanobubble pressure difference increases with the decrease of bubble radius. Using the Laplace-Young equation, we can estimate the liquid-vapor surface tension for the nanobubble ( fig. 5(b)). The figure shows that the surface tension for surface nanobubble is also nearly independent on its radius as for bubbles in the bulk solution [26], and is close to the planar value of γ = 0.31 [27] which was calculated from independent MD runs. Note that different from the simple fluid studied here, the curvature effects on surface tension for a more complex fluid, such as water, are non-negligible for small nanobubbles [28]. In contrast, the Tolman length for a LJ fluid has been estimated to be of order 0.3 σ [29], which is why we cannot observe curvature effects for bubbles with radii from 10 σ to 45 σ.
In summary, we used MD simulations to probe the relationship between the critical nucleus size and the radius of metastable surface nanobubble. We showed that if the curvature effect is absent, the radius of a stable nanobubble is approximately equal to the critical nucleus of the homogeneous liquid-to-vapor transition. We find that the nanobubble surface tension and pressure difference across the vapor-liquid interface show the expected dependence on the radius of the bubbles. The radius of the surface nanobubble is always equal to that of a "critical" bubble. We further demonstrate that when the curvature effects are present, measuring a surface with many pores of different sizes can be used to determine the curvature correction to the surface tension. Thus, we propose that surface nanobubbles can be used to probe the properties of the critical nuclei that transiently appear in bulk nucleation. Hence, our findings suggest that the study of pinned nanobubbles will yield direct insight into the properties of the elusive critical nucleus in bulk nucleation.
We note that, the main idea behind this work can easily be carried over the study of crystal nuclei of other phases, e.g. crystals. This could be extremely useful as direct (AFM or TEM) studies of true critical nuclei are virtually impossible. Within CNT, the present approach makes it possible to create stable, partial crystal nuclei with the same radius as the true critical nuclei. In the presence of curvature effects, a study of the dependence of the radius of the metastable nanocrystals on the width of the pore would make it possible to quantify the curvature corrections to the solid-liquid interfacial free-energy density.
This work is supported by National Natural Science Foundation of China (Nos. 21276007 and 91434204) and by the European Union's Horizon 2020 research and innovation programme under grant 674979-NANOTRANS.

Author contribution statement
QX, YL, ZG and ZL performed the calculations and analyzed data; JD, DF and XZ analyzed data and wrote the manuscript.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.