Mechanical Stability of Surface Nanobubbles.

Bubble cavitation is important in technologies such as noninvasive cancer treatment and diagnosis, surface cleaning, and waste-water treatment. The cavitation threshold is the critical external tensile pressure that induces unstable growth of the bubble. Surface nanobubbles have been previously shown experimentally to be stable down to -6 MPa, in disagreement with the Blake threshold, which is the classical cavitation model that predicts bulk bubbles with radii ∼100 nm should be unstable below -0.6 MPa. Here, we use molecular dynamics to simulate quasi-two-dimensional (2D) and three-dimensional (3D) nitrogen surface nanobubbles immersed in water, subject to a range of pressure drops until unstable growth is observed. We propose and assess new cavitation threshold models, derived from mechanical equilibrium analyses for both the quasi-2D and 3D cavitating bubbles. The discrepancies from the Blake threshold are attributed to the pinned contact line, within which the surface nanobubbles grow with constant lateral contact diameter, and consequently a reduced radius of curvature. We conclude with a critical discussion of previous experimental results on the cavitation of relatively large surface nanobubbles.


■ INTRODUCTION
Cavitation is a complex phenomenon in fluid mechanics where a bubble of dissolved gas or vapor forms within a liquid, usually as a result of a local drop in pressure. 1 Repeated formation and collapse of cavitating bubbles can be destructive to adjacent solid surfaces, such as in turbo-machinery, over prolonged periods of time. 1,2 However, there are also beneficial uses proposed for cavitation, including noninvasive cancer treatment, 3 surface cleaning, 4 and industrial waste-water treatment. 5 Heterogeneous nucleation is the formation of cavitation bubbles on the surface of an immersed substrate. 1 The acoustic cavitation threshold is the liquid pressure required to induce unstable growth of a bubble and has been modeled by the crevice theory for bubbles trapped in microscopic cracks of immersed solids. 6 This theory has been shown to predict the acoustic cavitation threshold in well-controlled, cylindrical pits down to 100 nm in diameter. 7 Surface nanobubbles (see Figure 1) rest on a solid substrate and have sizes of the order of 1−100 nm. They have recently attracted interest because their long lifetimes seemingly go against previous understanding of bubble stability. Due to the relatively stronger effects of surface tension at decreasing length scales, surface nanobubbles tend to form spherical cap shapes. Their lifetimes are up to 10 orders of magnitude longer than predicted for equivalently sized spherical nanobubbles. 8,9 It has been suggested previously that the diffusive stability of surface nanobubbles is due to an organic skin composed of amphiphilic molecules or other contaminants. 10,11 However, surface nanobubbles have been shown experimentally to exist stably without such contaminants. 12 Experimental, 13−15 simulation, 16−18 and analytical 19,20 investigations have all concluded that pinning of the three-phase contact line (CL) and a local supersaturation of dissolved gas in the surrounding bulk liquid are essential for both mechanically and diffusively stable equilibria of these bubbles.
The pressure inside a surface nanobubble (or any bubble) in steady-state mechanical equilibrium can be estimated by the Young−Laplace equation where P g is the internal gas bubble pressure, P v is the vapor pressure inside the bubble, P ∞ is the surrounding, far-field liquid pressure, γ is the liquid/gas interfacial surface tension, and R is the bubble radius of curvature. When P ∞ drops below a cavitation threshold, P ∞,c , it can be shown that the equality in eq 1 cannot hold for any bubble size, and the nanobubble experiences unstable, uncontrolled growth. The cavitation threshold of spherical macroscale bubbles can typically be predicted by the Blake threshold, which depends on R. 21 However, experiments have reported that surface nanobubbles are stable under impinging shock wave pressures of −6 MPa, despite the Blake threshold predicting a cavitation threshold of −0.55 MPa for R ∼ 100 nm. 22 So, what is the threshold pressure drop which induces unstable growth in surface nanobubbles? Our aim in this paper is to investigate the stable and unstable growths of surface nanobubbles using molecular dynamics (MD) and compare the observed threshold for stability with that predicted by considering the balance of pressures in eq 1 as the bubble grows.

■ SIMULATION METHODOLOGY
The LAMMPS MD software was used in this work to simulate surface nanobubble growth. 23 We performed simulations in quasi-two-dimensions (hereon referred to just as "2D") and three-dimensions (3D), with some variations in the corresponding theoretical models to account for the cylindrical and spherical cap shapes, respectively, and the different Laplace pressure contributions. An example of the MD setup for a 2D surface nanobubble is shown in Figure 1.
The internal gas and surrounding bulk liquid phases were chosen to be nitrogen (N 2 ) and water (H 2 O), respectively, given their relevance to most of the experimental literature. 3,5,9,13,25−27 All of the atoms comprising the fluid were contained between upper and lower solid walls. Each surface nanobubble was first equilibrated on the lower wall, which was textured with alternating patterns of hydrophobic (S o ) and hydrophilic (S i ) atom types with a regular "pattern wavelength", λ S = 3.14 nm. The equilibrium, gas-side contact angles of the pure S o and S i substrates were θ = 43°and θ = 99°, respectively, as measured from individual MD presimulations. 28 The substrate patterning allowed the contact line to expand to new pinning sites during bubble growth. 29 A neutrally charged, two-site nitrogen model 30 where ϵ ij and σ ij are the LJ potential well depth and characteristic length scale, respectively, between atoms i and j, separated by a distance r ij ; q i and q j are the charges of atoms i and j, respectively; ϵ 0 is the permittivity of free space. 23 Potentials were truncated and shifted, with relatively large cutoff radii (r c ) of 1.65 and 1.45 nm for the LJ and Coulombic potentials, respectively, to accurately model molecular effects in a multiphase system. 16 Interaction parameters between the H 2 O molecules and the S o and S i substrate atoms were previously determined to provide the desired wettabilities, as mentioned above. For the gas−solid energy parameters (S o −N; S i −N), we use similar values from refs 16, 32. All of the LJ and Coulombic parameters are listed in Table 1.
All of the simulations were run in an NVT ensemble with the Nose−Hoover thermostat applied to all of the water molecules at constant temperature, T = 300 K. The thermostat operated on y-component velocities in the 2D simulations. In the 3D simulations, the thermostat was applied to all velocity components, relative to the total water molecules' center of mass velocity. All liquid and gas molecules were equilibrated to 300 K initially, but the thermostat was removed from the nitrogen molecules during the pressure drop to allow the gas phase to expand in a more realistic manner without forcing isothermal or adiabatic behavior, which are commonly assumed in cavitation analyses. 1,7 A time step of Δt = 1 fs was used throughout; time integration was performed by the velocity Verlet algorithm. 23 The S o and S i substrate atoms were kept rigid.
All surface nanobubbles were initialized from a cap shape to reduce computational time and then allowed to equilibrate for 1.4 ns at an initial surrounding pressure of P ∞,0 = 10 MPa. The bulk liquid phase was doped with N 2 molecules to achieve supersaturation, otherwise the surface nanobubble would dissolve during equilibration. 16,19,20,33−35 The CL was pinned at a lateral diameter ϕ L = 14.11 nm (equivalent to 4.5λ S ) by a single S o patch for initialization of the system, then switched to a textured substrate for equilibration. The end state of the equilibrated bubble was then used as the starting point for all subsequent pressure drop simulations.
a LJ parameters for atoms in bold are assumed for pairs of like atoms. Any interaction pairs not given are equal to zero.
The drop in pressure was applied as a smooth hyperbolic tangent function over a period of 0.1 ns to prevent any separation of the water from the top channel wall; the period of the pressure drop was still small enough that it could be considered a step change. Care was taken not to run the simulations too long in order to ensure that growth was purely mechanical, i.e., due to the pressure imbalance in the gas, liquid, and curved surface tension components, not due to an increase in bubble mass by diffusion. Diffusive growth of surface nanobubbles has been demonstrated to occur over a time scale of microseconds, 8,9,20,36 so should not be a problem in this present setup.
We considered the vapor pressure, P v , inside the bubble to be negligible compared to the ∼MPa contribution of 2γ/R at the nanoscale. So, the tensile load applied during the surface nanobubble MD simulations was assumed to produce P v − P ∞ ≈ −P ∞ . 1 Periodic boundaries were applied in the x and y directions, whereas fixed boundaries were used in the z direction. During the pressure drop, the fluid domain was allowed to expand in the z direction by vertical motion of the upper channel wall, subjected to a force F = P ∞ A xy /N, where A xy is the area of the piston in the x, y plane and N is the number of atoms in the upper channel wall.
The 2D/3D steady-state surface nanobubble shapes resembled cylindrical/spherical caps. 16,19,25 To measure R, ϕ L , and θ from the MD simulations, cylindrical/spherical profiles were fitted to the 50% isodensity contour of the normalized fluid density (see Figure 1). 16,37 ■ RESULTS AND DISCUSSION Two-Dimensional Cavitation Threshold. We performed 13 simulation cases in which an equilibrated 2D surface nanobubble was subjected to a drop in pressure to final absolute values in the range +3 to −4.5 MPa. Negative values indicate tensile pressures in which the upper channel wall (piston) forcing was directed in the positive z direction. Each target pressure was sustained for 15 ns after the pressure drop, during which the bubble either reaches a mechanically stable size or exhibits unstable growth. A total of 43 000 H 2 O molecules and 1200 N 2 molecules were used in each of the 2D simulations.
All our simulations started from the same equilibrated surface nanobubble at P ∞,0 = 10 MPa. The radius of curvature and gas-side contact angle for this equilibrated bubble were R 0 = 10.52 nm and θ 0 = 43°, respectively, and the bubble was pinned to the substrate at a lateral contact diameter ϕ L,0 = 14.37 nm (see Figure 1). The initial gas pressure, P g,0 , is estimated from the 2D equivalent of eq 1, where the right hand side becomes γ/R and using P ∞ = P ∞,0 and R = R 0 .
During the simulations, the surface nanobubbles grew with a pinned contact line, i.e., constant contact radius (CCR) growth, until some time at which the contact line expanded out rapidly to the next pinning site. For macroscopic bubbles and droplets, this behavior is typically known as "stick-slip", although in our surface nanobubble simulations the expansion of the CL is so rapid compared to the time spent during CCR growth that we consider this "stick-jump" behavior. 16 The variations in the contact diameters in three bubble growth cases are shown in Figure 2. The lateral contact diameters ϕ L of the surface nanobubbles were observed to grow in finite increments of λ S , i.e., ϕ L,n = ϕ L,0 + nλ S , where n is an integer. The discrete pinning sites ϕ L,n are labeled in Figure 2 Surface nanobubbles experiencing a final applied pressure of −4 MPa and below were found to exhibit unstable growth and Langmuir XXXX, XXX, XXX−XXX outgrew the domain within 5 ns. All of the other cases reached stable steady states within around 3 ns. The variation in the steady-state bubble cross-sectional areas, A, with applied pressures, P ∞ , is shown in Figure 3. The lowest liquid pressure that sustained a stable surface nanobubble had a nominal applied pressure of −3.75 MPa, so we conclude the cavitation threshold to be between −3.75 and −4 MPa. We compare the results of the MD simulations to the 2D Blake threshold equation 21 (see the derivation in the Supporting Information) where P v is the vapor pressure at the bulk liquid temperature and γ = 57.35 mJ/m 2 from a preliminary MD simulation of a simple plane interface system under equilibrium conditions; k is the constant exponent in the polytropic gas equation, P g A k ∝ P g V k = const., which is estimated from the MD simulations (see Figure 3) to be k = 1.18 ± 0.14. 1,21 The lower bound of this estimate is close to 1, which is indicative of isothermal expansion; indeed, the stable bubble expansion cases were found to grow with near constant temperature, as shown in Figures S3a,b and S4. The unstable cases exhibited a more substantial decrease in temperature, as shown in Figures S3c  and S4, suggesting that past the stability threshold the bubble expansion cannot be considered isothermal. It is likely that bubbles well beyond the threshold tend to more adiabatic behavior, where the growth rate (as in Figure 2c) is too rapid for adequate heat transfer to the gas to occur. With these parameters, the Blake threshold pressure is calculated using eq 3 to be −0.78 MPa, which is approximately 3 MPa higher than the lowest sustained pressure from our MD simulations. This result is consistent with experiments, which have shown Blake's threshold to be inadequate for predicting the acoustic superstable cavitation threshold of surface nanobubbles. 22 The Blake prediction is also included in Figure  3. where is the initial crosssectional area of the 2D surface nanobubble. There is no analytical solution to eq 4 in terms of θ c , so we solve this numerically.
The maximum contact angle observed in all of the MD simulations was just less than 90°. It has been previously suggested that surface nanobubbles cannot be in stable diffusive equilibrium with gas-side contact angles greater than 90°, although they can still return to a point of diffusive stability with an angle less than 90°. 20,38 There is therefore no reason why they could not sustain a gas-side contact angle greater than 90°for short time scales in mechanical equilibrium. However, given how close this is to the equilibrium contact angle of the S i substrate, we assume that the angle at which the CL would unpin and jump to the next pinning site is closely linked to the S i surface's wettability. Here, to match the simulations, we assume that the maximum contact angle that can be sustained is 90°, which imposes a constraint on eq 4.
Once θ c is determined from eq 4, the threshold pressure for each pinning site can be calculated from  Equations 4 and 5 predict the critical size and threshold pressure, respectively, for a bubble on a single pinning site with lateral contact diameter ϕ L . The surface nanobubbles in our simulations expanded to new pinning sites as they grew, so the assumption of CCR growth breaks down in our model. However, the benefit of writing eq 4 in terms of P g,0 and A 0 (at the original pinning site, ϕ L,0 ) means that θ c can then be found for any of the succeeding pinning sites ϕ L,n . By substituting θ for θ c and P ∞ for P ∞,c in eq 5, the quasistatic growth path of a surface nanobubble can be found, relating the far-field liquid pressure to the equivalent mechanically stable shape at a particular pinning site. These growth paths are the solid black lines in Figure 3. The corrected threshold pressure from eq 5 is determined to be P ∞,c = −3.79 MPa, with the critical bubble angle θ c = 90°w hile pinned at ϕ L,1 . This threshold is the lowest point in the predicted pressure lines in Figure 3 and is in good agreement with the lowest stable pressure in our MD simulations. The growth paths predicted by our corrected threshold model agree well with the MD results too. The size of the bubble at the threshold, as predicted by our model, is termed the critical size and is also marked in Figure 3. Our model predicts the contact angle, radius of curvature, and pinned contact diameters, which agree well with the MD results in Figure 4a−c, respectively. Figures 3 and 4 show that the surface nanobubble jumps from its original pinning site ϕ L,0 to the next pinning site ϕ L,1 at around 84°. The nanobubble growth then continues to agree with our corrected threshold model for the new pinning site ϕ L,1 before reaching its critical size, this time at θ = 90°.
The growth path of surface nanobubbles at new pinning sites is predicted well by our model, even if it is not possible to determine at what contact angle they jump. The simulated surface nanobubble jumped from a contact angle of 84°at the original pinning site to 80°at the new site. Because this is only a small decrease in contact angle, it could be possible to model this stick-jump growth as constant contact angle for larger bubbles as λ S /ϕ L → 0. However, for the purposes of this study, we considered all growths to be stick-jump CCR acting at discrete pinning sites. 36 In Figure 4b, the radius of curvature is shown to decrease, as the bubble grows on a particular pinning site. This is expected for a cylindrical/spherical cap shape with constant ϕ L and θ ≤ 90°. However, the 2D Blake threshold, eq 3, does not capture this growth behavior accurately because it assumes bubble growth is a monotonically increasing function of R. From eq 1, the surface tension component of the pressure balance varies ∼γ/R, so a decrease in R during surface nanobubble growth means an increase in the surface tension contribution. This suppresses excessive growth and enables the surface nanobubble to withstand lower pressures than the Blake threshold for a spherically equivalent bubble. This is the main reason why the Blake threshold does not correctly predict the surface nanobubble cavitation threshold. Figure 4c shows the contact diameters of the surface nanobubbles do remain pinned to discrete positions during growth, as we expected from the textured S i −S o substrate patterning. The jump in the contact diameter is equal to the pattern wavelength, λ S = 3.14 nm. Given that our model has been shown to agree with MD simulations, we then used it to predict the cavitation thresholds for a range of surface nanobubble sizes (θ 0 and ϕ L,0 ), with other parameters the same as in the 2D MD simulations. These results are shown in Figure 5. Decreasing θ 0 and ϕ L,0 is equivalent to reducing A 0 and typically lowers the predicted cavitation threshold. Increasing k from 1 (isothermal expansion) to 1.4 (adiabatic) also tends to lower the predicted cavitation threshold, which is a similar trend to the Blake model. 1,21 The inset in Figure 5 shows that the threshold tends Langmuir XXXX, XXX, XXX−XXX to P v for large bubbles, which is similar to the Blake threshold in eq 3 for large bubbles.
Three-Dimensional Cavitation Threshold. To assess the validity of our corrected threshold model, a 3D surface nanoubble in MD was also equilibrated on a substrate, patterned with alternating S i and S o concentric rings, as shown in Figure 6. The patterning for this 3D case was chosen to be equivalent to the axisymmetric patterning used in the 2D simulations and likewise allowed the CL to be pinned at a particular diameter ϕ L for CCR growth. 37 A total of 230 000 H 2 O molecules and 3100 N 2 molecules were used in each of the 3D simulations. At an initial equilibrated pressure P ∞,0 = 10 MPa, the radius of curvature and contact angle of the 3D bubble were R 0 = 10.5 nm and θ 0 = 46°, respectively, pinned to the substrate with ϕ L,0 = 15.22 nm. The properties P v , γ, and k were assumed identical to the 2D case. Equation 1 is used to estimate P g,0 , using P ∞ = P ∞,0 and R = R 0 .
Due to the significant computational cost of the 3D simulations, only two pressure drop cases were investigated, with final applied absolute pressures of −7.5 and −10.5 MPa, based on conservative estimates of the cavitation threshold model extended to 3D (which is discussed below). Both simulations were run for a maximum of 5 ns: given that our 2D simulations reached either a mechanically stable state or exhibited unstable growth within the first 3 ns (see Figure 2), we deemed this sufficient simulation time. The transient expansions in the bubble volume for these two cases are shown in Figure 7. The case with −7.5 MPa applied pressure was stable, whereas the case with −10.5 MPa was unstable; the cavitation threshold for the 3D surface nanobubble must be between these two values. A video showing the bubble growths for the two cases can be found in the Supporting Information.
The 3D form of the Blake threshold (see the derivation in the Supporting Information) is and we use the parameters taken from the 3D MD simulation. Equation 6 predicts a cavitation threshold of −4.08 MPa,   is the initial volume of the surface nanobubble. As in the 2D case, there is no analytical solution for θ c , so we solve eq 7 numerically. Once obtained, θ c can then be used to find the threshold pressure for each particular pinning site, i.e., The θ c ≤ 90°constraint is applied again, as in the 2D case. One further modification for the 3D model is that ϕ L should now jump in increments of 2λ S , since the system is axisymmetric, i.e., ϕ L,n = ϕ L,0 + 2nλ S . This corrected threshold model predicts a cavitation threshold of −9.18 MPa, which is in good agreement with the range observed in the MD simulations. The predicted variation in the cavitation threshold with surface nanobubble size, using the same parameters as in the 3D MD simulations, is shown in Figure 8. Decreasing the initial contact angle, θ 0 , produces a lower cavitation threshold. This could be because the initial radius of curvature will be larger and hence the Laplace pressure would increase relatively more as the bubble grows, making it more stable to cavitation. Qualitatively, increasing the polytropic exponent, k, from an isothermal expansion (k = 1) to an adiabatic expansion (k = 1.4) tends to lower the cavitation threshold, which is similar in the Blake threshold prediction too. 21 The inset in Figure 8 shows how the threshold decreases as the bubble gets larger and tends to P v ; eq 6 shows that the Blake threshold also tends to P v for large bubbles.
Previous experiments have indicated that surface nanobubbles are stable down to at least −6 MPa. 22 Using our model, the predicted variation in the cavitation threshold with surface nanobubble size at experimental conditions, i.e., P ∞,0 = 0.1 MPa and γ = 71.69 mJ/m 2 , is shown in Figure 9a, with other model parameters taken from the previous 3D simulation. The cavitation threshold is strongly dependent on both the initial contact angle, θ 0 , and the initial contact diameter, ϕ L,0 . However, it is clear that typically decreasing θ 0 and ϕ L,0 decreases the cavitation threshold, making the nanobubbles more stable. For ϕ L,0 < 20 nm, we predict that most surface nanobubbles would be resistant to a −6 MPa pressure fluctuation, and for some cases, e.g., θ 0 = 10°and ϕ L,0 = 10 nm, the cavitation threshold is predicted to be as low as −28 MPa. For bubbles, the same initial size and shape as in the 3D MD simulations, i.e., θ 0 = 46°and ϕ L,0 = 15.22 nm, the cavitation threshold is predicted to be P ∞,c = −15 MPa under typical experimental conditions. Although seemingly low, this    Figure 9a shows that we predict most surface nanobubbles should cavitate by −6 MPa pressure, provided that ϕ L,0 > 50 nm. Yet the published experiment 22 reported surface nanobubbles of sizes up to ϕ L,0 = 300 nm that still did not undergo cavitation at −6 MPa.
As already noted, the Blake model gives a poor prediction of the cavitation threshold when using the nanobubble's radius of curvature, R 0 , in eq 6. However, consider instead a free spherical bubble of equal mass to the surface nanobubble, with a radius of curvature, R m ; combining eq 1, the polytropic gas law, and the volume of a spherical bubble enables us to obtain the radius of this equivalent spherical bubble numerically from Then by substituting R m for R 0 in eq 6, the Blake threshold is found to be many MPa lower than for a surface nanobubble of equal mass, as seen in Figure 9b. This means that if a surface nanobubble detached from the surface and formed a spherical bubble, while still remaining very close to the surface, there is a larger range of tensile pressures within which the bubble could still be stable. 10 We also acknowledge the possibility of contamination from polymer droplets in the cited experiments, 22 which could be mistaken for surface nanobubbles and lead to anomalous results, particularly for the larger bubbles. 9,40 ■ CONCLUSIONS The cavitation threshold for unstable growth of a surface nanobubble cannot be predicted using the classical Blake threshold equation. Here, we have proposed a corrected cavitation threshold model, which assumes that the surface nanobubble expands with a pinned lateral diameter, i.e., the constant contact radius mode of growth. The decreasing radius of curvature increases the surface tension component of the internal bubble pressure. The growth paths of 2D surface nanobubbles can be predicted by our model and are in good agreement with molecular dynamics simulations of air bubbles in water.
The cavitation threshold of 3D surface nanobubbles predicted by our model is also in agreement with molecular dynamics simulations. For experimental conditions, i.e., γ = 71.69 mJ/m 2 and P ∞,0 = 0.1 MPa, the model predicts a cavitation threshold of P ∞,c ≈ −15 MPa for small surface nanobubbles (ϕ L = 15 nm, θ 0 = 46°), which is many MPa lower than the classical Blake threshold prediction.
Our model also predicts that larger surface nanobubbles considered in experiments (i.e., 60 nm < ϕ L < 300 nm) should become unstable in the constant contact radius mode of growth. We hypothesized a possible mechanism for the experimentally observed stability, in which unstable surface nanobubbles detach from the surface during their unstable growth phase and form more stable bulk nanobubbles with reduced radius of curvature, R m . Whether or not this proposed detachment process is feasible requires further investigation.
The cavitation threshold model we have presented here assumes that the maximum contact angle a surface nanobubble can sustain is θ = 90°, which is the equilibrium contact angle on the hydrophilic region of the patterned substrate. Future work could vary the surface/liquid interatomic potentials to assess whether wettability has a significant effect on the cavitation threshold.