Quantifying interfacial tensions of surface nanobubbles: How far can Young's equation explain?

Abstract

Nanobubbles at solid–liquid interfaces play a key role in various physicochemical phenomena and it is crucial to understand their unique properties. However, little is known about their interfacial tensions due to the lack of reliable calculation methods. Based on mechanical and thermodynamic insights, we quantified for the first time the liquid–gas, solid–liquid, and solid–gas interfacial tensions of submicron-sized nitrogen bubbles at graphite–water interfaces using molecular dynamics (MD) analysis. It was revealed that Young's equation holds even for nanobubbles with different radii. We found that the liquid–gas and solid–liquid interfacial tensions were not largely affected by the gas density inside the nanobubbles. In contrast, the size effect on the solid–gas interfacial tension was observed, namely, the value dramatically decreased upon an increase in the gas density due to gas adsorption on the solid surface. However, our quantitative evaluation also revealed that the gas density effect on the contact angles is negligible when the footprint radius is larger than 50 nm, which is a typical range observed in experiments, and thus the flat shape and stabilization of submicron-sized surface bubbles observed in experiments cannot be explained only by the changes in interfacial tensions due to the van der Waals interaction-induced gas adsorption, namely by Young's equation without introducing the pinning effect. Based on our analysis, it was clarified that additional factors such as the differences in the studied systems are needed to explain the unresolved open issues – a satisfactory explanation for the nanobubbles in MD simulations being ultradense, non-flat, and stable without pinning.

---

Introduction

Bubbles at solid–liquid interfaces are ubiquitous and important both in nature and in industry. Accordingly, a rigorous understanding of their physics from the molecular scale is of fundamental importance for various physicochemical phenomena involving bubble nucleation, such as boiling, cavitation, chemical reactions, and water splitting. In particular, small bubbles at solid–liquid interfaces with a submicron-sized footprint radius and a height of tens of nanometers, which are called surface/interfacial nanobubbles,¹ were first reported a quarter century ago² and have attracted increasing attention because of their roles for a wide range of phenomena, such as long-range hydrophobic interactions,³ fluid slippage,⁴ and stabilization of water thin films in nanoscale spaces under vacuum.⁵ The utilization of surface nanobubbles is also a promising way for engineering applications including the wafer-scale transfer of graphene films,⁶ the onset of nucleate boiling at a low superheat,⁷ the drag reduction of fluid flow inside a microchannel,⁸ and inhibition of metal corrosion in acidic geothermal fluids.⁹

To date, a number of studies have reported the unique properties of surface nanobubbles, such as the unexpected day-order long-lifetime,^10–12 significantly larger contact angles than those observed on a macroscale,¹³ and superstability against the disturbance.^14,15 Subsequently, a variety of theories have been proposed to explain them.^1,16 So far, it has been speculated that pinning of their three-phase contact line and local gas oversaturation are essential factors to explain the characteristics,^17,18 which was first proposed by Liu and Zhang in 2013.¹⁹ Indeed, it has been experimentally confirmed that surface nanobubbles are strongly pinned even on atomically flat surfaces.^13,15,20,21 However, some have questioned the necessity of the pinning for the long-lifetime of surface nanobubbles. Chen et al. carried out molecular dynamics (MD) simulations of steady argon surface nanobubbles without pinning sites and argued that the pinning was not necessarily needed for their stabilization.²² Nag et al. experimentally observed a free growth/shrinkage of surface nanobubbles for longer than 1000 seconds in a nano-confined water film using an in situ liquid phase electron microscopy (LPEM) technique, which indicated the existence of non-pinned surface nanobubbles.²³ Thus, it is still under debate whether the pinning is necessary to explain the unique properties of surface nanobubbles.

In addition to the aforementioned properties, a new characteristic of surface nanobubbles was experimentally reported: an ultrahigh density of inner gas molecules. Zhou et al. revealed using scanning transmission X-ray microscopy (STXM) that oxygen surface nanobubbles with a footprint diameter as large as 500–800 nm had a density one to two orders of magnitude higher than that under ambient conditions, which cannot be explained by the ideal gas law and the Laplace pressure.²⁴ This new feature would lead to innovations such as effective oxygen transport, storage of hydrogen gases at room temperature and pressure, and improvement in chemical reactions. However, the cause for the ultrahigh density has not been understood yet. Moreover, Zhang et al. using attenuated total internal reflection-Fourier transform infrared (ATR-FTIR) spectroscopy reported that the average density inside the CO₂ surface nanobubbles was 44 ± 16 mol m⁻³ (i.e., 1.2 ± 0.4 kg m⁻³), equivalent to the gas density at atmospheric pressure,²⁵ which is in contrast to the aforementioned STXM results.

The ultrahigh inner density of surface nanobubbles also raises a new question whether their interfacial tensions change. Because interfacial tensions originate from the interaction among molecules at the interfaces, it is natural to expect that the values of ultradense nanobubbles deviate from those of macroscopic bubbles. Such a change in the interfacial tensions could rationalize the unique characteristics of surface nanobubbles without introducing pinning as an additional cause. For example, Zhang et al. proposed that the theoretical lifetime of surface nanobubbles would approach the timescale observed in experiments if the liquid–gas interfacial tension could be reduced due to the high inner pressure.²⁶ Petsev et al. also theoretically predicted that the adsorption of gas molecules onto the solid surface just under the surface nanobubbles could reduce the solid–gas interfacial tension, leading to their distinctive flat shape.²⁷

To reveal the mechanisms underlying the properties of surface nanobubbles, quantification of their interfacial tensions is highly desired. However, direct experimental measurement of the interfacial tensions of surface nanobubbles is extremely difficult. Specifically, the liquid–gas interfacial tension has been measured in a few studies, but the results are inconsistent: one study reported that the value was not changed,²⁸ while another reported that it was changed.²⁹ More importantly, quantitative values of the solid–gas and solid–liquid interfacial tensions of such nanobubbles have never been reported. To calculate the reliable interfacial tensions of surface nanobubbles, it is desirable to employ a theoretical approach applicable at the molecular scale, rather than experiments that inherently involve various uncertainties.

In this study, by integrating well-proven mechanical and thermodynamic insights^30–37 into MD simulations, we quantified the liquid–gas, solid–gas, and solid–liquid interfacial tensions of nitrogen surface nanobubbles at graphite/water interfaces. The resulting contact angles evaluated by substituting these three interfacial tensions into Young's equation agreed well with those measured from the apparent shapes of the surface nanobubbles. We further extended our method to estimate the contact angles not only for quasi-two-dimensional nanobubbles but also for semispherical submicron-sized bubbles, which enabled us to directly compare our simulation results with the previous experimental results. Based on this analysis, as well as the interfacial tensions, we shed light on the unresolved issues of surface nanobubbles mentioned above: the necessity of the pinning and the effect of the inner gas molecules on the flatness and long-lifetime.

Results and discussion

Density-dependence of the contact angles of surface nanobubbles

We have carried out MD simulations of quasi-two-dimensional hemi-cylindrical surface nanobubbles composed of nitrogen molecules (referred to as N₂ molecules hereafter) in water, which are confined between the top graphene and the bottom AB stacked four-layer graphene (referred to as graphite hereafter). Fig. 1 presents an example of the nanobubble simulation systems in the present study. Five cases were tested for the number of N₂ molecules in the system N_N2, as N_N2 = 400, 300, 250, 200, and 150 with a constant number of water molecules N_H2O = 11 000. The top graphene worked as a piston to attain the pressure in the liquid at the atmospheric pressure P₀ = 1 atm, as it has been successfully used in previous studies.^22,38,39 The Lennard-Jones (LJ) energy between the carbon atoms in the bottom graphite and N₂ molecules was changed by multiplying a coupling parameter η = 1.0 or 1.5 to investigate the influence of the solid–gas interaction. The details of the setup of these nanobubble simulations are described in the Methods section.

Fig. 1 Simulation system of an N₂ nanobubble in water confined between graphite and hydrophilic graphene with a number of N₂ molecules N_N2 = 400.

Fig. 2 shows (i) equilibrium MD snapshots and (ii) the time averaged density distributions of N₂ surface nanobubbles and a nanodroplet with only 5000 water molecules on the graphite surface. As the number of N₂ molecules decreased, the nanobubbles became smaller and the N₂ density in the bulk of the nanobubbles, , increased. The values were approximately 130–230 kg m⁻³ and 100–170 kg m⁻³ for η = 1.0 and 1.5, which are two orders of magnitude higher than the nitrogen gas density under room conditions (1.14 kg m⁻³) and reaches 12–28% of liquid nitrogen density (809 kg m⁻³). This ultrahigh density is caused by the high Laplace pressure generated by their small radii of curvature. The inner pressures of these nanobubbles P_gas, which were calculated from the normal stress working on the graphite underneath them, were indeed very high, ranging from 9.0 to 23.7 MPa (see ESI Fig. S1†). In addition, adsorbed layers with higher concentrations of N₂ molecules were observed near the solid surface underneath the nanobubbles, which is consistent with the experimental results.⁴⁰

Fig. 2 (i) MD snapshots of nanobubbles with a number of gas molecules N_N2 = 400, and (ii) density distributions of nanobubbles with N_N2 = 400, 250, and 150 (from the top) in the cases of (a) η = 1.0 and (b) η = 1.5, and (c) a water nanodroplet equilibrated on the same graphite. Time average was taken for 20 ns after 40 ns equilibration for the calculation of the density distribution.

Furthermore, we evaluated the apparent contact angle of the nanobubbles θ_app and that of the nanodroplet θ_droplet by least-squares circular fitting to a contour of the gas–liquid interface away from the wall.^30,33–35 Masuduzzaman and Kim reported the importance of the definition of the solid–liquid boundary position for the evaluation of the apparent contact angle.⁴¹ In this study, the solid–liquid interface was defined as the position where the fluid density becomes zero which gives a consistent mechanical balance near the three-phase contact lines.³⁰ The contact angle θ_app increased as the nanobubbles became smaller, i.e., as increased (89.2–98.9° for η = 1.0 and 109.8–121.0° for η = 1.5). They were much smaller than those observed in experiments (160–175°)¹ but clearly higher than the contact angle of the nanodroplet θ_droplet = 71.7° shown in Fig. 2(c-ii), which is within a range of that on graphite surfaces experimentally reported.⁴² Because the pinning effect can be negligible on an atomically smooth solid surface,³⁰ this density-dependent change in the contact angles indicates that the gas molecules inside the nanobubbles induced the change in their interfacial tensions.

Quantitative evaluation of interfacial tensions

To quantitatively investigate the effect of gas molecules on the interfacial tensions, we employed quasi-one-dimensional liquid–gas, solid–gas, and solid–liquid systems that mimic the interfaces of surface nanobubbles shown in Fig. 2. Specifically, the liquid–gas and solid–gas interfacial tensions were calculated based on Bakker's equation,^30–37 which describes the relationships between the fluid stress anisotropy and interfacial tensions. The solid–liquid interfacial tension was quantified by the thermodynamic integration (TI) method that calculates the minimum work needed to separate the liquid from the solid surface per area, which corresponds to the interfacial free energy.^30,34,35,43 These types of calculation methods of interfacial tensions are called the “mechanical route”³² and “thermodynamic route”,^30,34,35,43 respectively. The details of the methods used for the calculation of interfacial tensions are described in the Methods section.

The liquid–gas, solid–gas, and solid–liquid interfacial tensions of the nanobubbles calculated for various bulk N₂ densities are shown in Fig. 3(a). Although the value of the liquid–gas interfacial tension γ_LG decreased from that of the liquid–vapor interfacial tension γ_LV = 53.9 ± 0.4 × 10^–3 N m⁻¹ (95% confidence interval) with the increase of , the difference γ_LV − γ_LG was as small as about 4 × 10^–3 N m⁻¹ (7% of γ_LV). This may be because the strong Coulomb force among the water molecules governed the interfacial tension and the effect of the weak van der Waals force between the water and N₂ molecules was relatively small. The value of γ_LV obtained by the present SPC/E water model is consistent with other calculated values previously reported in the literature.^44,45 The calculated values of γ_LV obtained with several water models and at different temperatures are known to be underestimated compared to the established experimental values.⁴⁶ On the other hand, the solid–gas interfacial tension relative to the solid–vacuum interfacial tension, γ_SG – γ_S0, dramatically decreased in both cases of η = 1.0 and 1.5 with increasing , i.e., with increasing pressure in the nanobubbles P_gas. This is because the density inside the gas molecule-adsorbed layer increased as P_gas increased (shown in ESI Fig. S2†). Although such an adsorption of gas molecules was also observed for the liquid–gas interfaces (Fig. 2(ii)), the density increase was much smaller than that for the solid–gas interfaces, resulting in a small effect on γ_LG. These results imply that the change in the contact angles observed in Fig. 2 was attributed to the significant decrease in the solid–gas interfacial tension caused by the strong adsorption of N₂ molecules.

Fig. 3 (a) Interfacial tensions and (b) contact angles of θ_Y and θ_app estimated by Young's equation (eqn (1)) and apparent shape of the nanobubbles as a function of the bulk N₂ density . The pink and green dashed lines indicate the values of γ_LV = 53.9 ± 0.4 × 10⁻³ N m⁻¹ and –(γ_SL − γ_S0) = 17.3 ± 1.0 × 10⁻³ N m⁻¹, respectively.

From these results, we compared the apparent contact angle θ_app obtained from the density distribution as shown in Fig. 2(ii) and the contact angle θ_Y estimated by Young's equation:

as a function of in Fig. 3(b). These values calculated via different approaches were in good agreement for both cases of η = 1.0 and 1.5, which assures the reliability of our calculation results and proves the conclusion that the density-dependence of the contact angles observed in Fig. 2 is attributed to the change in γ_SG – γ_S0. This result also highlights the important point that Young's equation holds even for nanoscopic gas bubbles at solid–liquid interfaces. Young's equation broke down only when the nanobubbles became so small that there was no bulk N₂ density region inside them where the interfacial tensions between the two bulk phases could not be defined any more (refer to ESI Fig. S3†).

Are submicron-sized surface bubbles flattened by gas molecular adsorption?

At first glance, the density dependence of the contact angle induced by the adsorption of gas molecules appears to be a key factor to explain the flatness of surface nanobubbles. However, because the nanobubbles in our MD systems are quasi-two-dimensional and their footprint radius (∼5 nm) is much smaller than that observed in experiments (∼500 nm), a direct comparison between them is not possible. Although it is currently not feasible to directly analyze semispherical submicron-sized bubbles in MD simulations due to the very demanding computational cost involved, our method, which solves Young's equation by quantifying the interfacial tensions corresponding to the Laplace pressure of three-dimensional nanobubbles, makes it possible to calculate their contact angles including the effects of gas molecular adsorption (the detail is shown in ESI Fig. S4–S6†).

The relationship between the contact angle θ_Y and the footprint radius R_FP of semispherical submicron-sized nanobubbles calculated by our method is shown in Fig. 4(b). When R_FP > 50 nm, the change in the contact angles from one of nanodroplet θ_droplet obtained in Fig. 2(c-ii) was small, at most 4.1° and 15.3° for η = 1.0 and 1.5, respectively. When R_FP < 50 nm, there is a significant increase in the contact angle as the R_FP value decreases, but only when R_FP < 1 nm, it reaches 160–175°, which are typical values for nanobubbles on graphite in experiments.¹

Fig. 4 (a) Relative solid–gas interfacial tension γ_SG − γ_S0 as a function of the pressure P_gas. The red and blue lines are drawn by fitting the Szyszkowski eqn (3). The light pink and light blue bands were determined to cover all data points, which have the ranges of K^A_eq between 0.38 × 10⁻⁷ and 0.65 × 10⁻⁷ Pa⁻¹ and between 1.31 × 10⁻⁷ and 2.89 × 10⁻⁷ Pa⁻¹, respectively. (b) A semispherical nanobubble's contact angle versus footprint radius R_FP. The red and blue curves are drawn by using Petsev's model (eqn (2)) with the K^A_eq value obtained in (a). The green dashed lines indicate θ_droplet shown in Fig. 2(c). Enlarged view for a small R_FP range is shown as an inset.

Petsev et al. theoretically proposed that the flat nanobubble morphology can be explained by considering the effect of gas adsorption at solid–gas interfaces on their contact angles.²⁷ The equation used can be rewritten as

where k_B is the Boltzmann constant and b is the cross-sectional area of the adsorbed molecule, and K^A_eq is the equilibrium adsorption constant with a unit of inverse pressure, which represents the strength of adsorption. Their study used the values of K^A_eq = 1.0 × 10⁻⁴ to 1.0 × 10⁻⁶ Pa⁻¹ to explain the flatness of the nanobubbles, while these values are for materials functionalized for gas adsorption^47,48 and would be too high for an atomically-flat and non-functionalized graphite surface typically used for the experimental observation of surface nanobubbles. Here, we estimated the values of K^A_eq for graphite surfaces used in our MD simulation using the Szyszkowski equation,⁴⁹ which predicts the effect of gas adsorption onto a solid surface on γ_SG – γ_S0 by combining the Langmuir isotherm and the ideal gas model as

The plots of γ_SG – γ_S0 as a function of the pressure P_gas (Fig. 4(a)) were fitted using eqn (3) with K^A_eq as the fitting parameter. The trends between the MD results and the theoretical prediction presented good agreement and the resulting values were K^A_eq = 0.43 × 10⁻⁷ and 1.53 × 10⁻⁷ Pa⁻¹ for η = 1.0 and 1.5, which were significantly smaller than the K^A_eq values used in the Petsev's study.²⁷ We then evaluated the simulated contact angle θ_Y by substituting the obtained K^A_eq value to Petsev's model (eqn (2)) with b = 10.74 × 10⁻² nm² estimated from the LJ length parameter used in our study. As shown in Fig. 4(b), the curves obtained using eqn (2) followed the plots of θ_Y well. If we follow the value of K^A_eq = 1.0 × 10⁻⁴ to 1.0 × 10⁻⁶ Pa⁻¹ used in a previous study,²⁷ the value of η needs to be considerably large, which induces the unrealistic behavior of N₂ molecules on graphite. From these results, the present K^A_eq value obtained in our analysis is more reasonable than the ones of the previous study, which seems to have overestimated the effect of gas adsorption on the contact angles. Therefore, although the size effect appears when the nanobubbles are sufficiently small, the flat shape of semispherical nanobubbles observed in experiments cannot be explained only by the van der Waals interaction-induced gas molecular adsorption.

We note that these results are still limited to the ideal case, namely, the nanobubbles composed of non-polar gas molecules in ultrapure water. In more realistic systems, the gas molecular adsorption should be affected by the existence of ions in solutions and the electric charge at interfaces. Further quantitative characterization of them remains an interesting topic but is beyond the range of the present paper.

Are submicron-sized surface bubbles ultradense or not?

Zhou et al. reported by STXM observation²⁴ that oxygen nanobubbles with a footprint radius of hundreds of nanometers had a density of tens to hundreds kg m⁻³. On the other hand, Zhang et al. reported by ATR-FTIR measurements that the inner densities of surface nanobubbles are almost the same as that of the ideal gas at atmospheric pressure.²⁵ In our MD analysis, we found that the ideal gas law is applicable to surface nanobubbles with a gas density of lower than approximately 100 kg m⁻³ (shown in ESI Fig. S6†). Therefore, by considering the ideal gas law and the Laplace pressure, the bulk N₂ density inside the surface nanobubbles is given bywhere m is the mass of a gas molecule (here we used 4.65 × 10⁻²⁶ kg of N₂ molecule) and we used T = 298 K and γ = 72 × 10⁻³ N m⁻¹ instead of using γ_LG obtained in the MD simulation as γ. The graph of the bulk N₂ density inside the nanobubbles estimated by eqn (4) is shown in ESI Fig. S7.†

If we adopt the typical values of surface nanobubbles, namely a footprint radius of 300 nm and a contact angle of 160°, the density is estimated to be 3.0 kg m⁻³ from eqn (4). Although this is about 2.6 times higher than the N₂ gas density at atmospheric pressure (1.14 kg m⁻³), it is still two orders of magnitude lower than the STXM results, supporting the ATR-FTIR ones. There are two possible reasons for this result. One is the overestimation of the gas density inside the nanobubbles in STXM observations caused by the existence of an ultradense gas molecule-adsorbed layer underneath them. The formation of the dense adsorbed layers of gas molecules has been reported not only by MD simulations^22,24,40 but also by AFM^40,50–52 and LPEM^53,54 experiments. Our analysis also found that even when the pressure P_gas was reduced to 0.2 MPa corresponding to the nanobubbles with a footprint radius of 250 nm and a contact angle of 170°, the peak density inside the adsorbed layer was 37 and 241 kg m⁻³ for η = 1.0 and 1.5 (shown in ESI Fig. S2†), which is 1–2 orders of magnitude higher than 2.3 kg m⁻³ of the bulk N₂ density inside the nanobubbles. Such highly dense adsorbed layers would have some effect on the measured absorption spectra. The other is the nano-confined effect. While ATR-FTIR experiments were conducted in open systems, the nanobubbles in STXM observations were inside a nanoscopic water film sandwiched in two Si₃N₄ windows like those in LPEM observation. Under such nano-confined conditions, the overlap of the electric double layer formed at the two solid–liquid interfaces would deviate the properties of the nanobubbles from those predicted by the classical theory and MD simulations, implying the importance of the system conditions to the properties of surface nanobubbles discussed in the next subsection.

Eqn (4) suggests that surface nanobubbles with a large footprint radius in open systems like those observed by AFM or optical methods cannot have an ultrahigh density of tens to hundreds kg m⁻³ except for inside the adsorbed layers. In contrast, because the typical footprint radius of nanobubbles observed by LPEM is at most tens of nanometers,^54–56 the bulk density inside them can reach a hundred kg m⁻³. In addition, when a contact angle is 90°, the gas density estimated using eqn (4) corresponds to that inside the bulk nanobubbles in which the footprint radius R_FP is defined as the radius of curvature. We can estimate that the typical bulk N₂ nanobubble with a diameter of 100 nm would have a gas density of 34 kg m⁻³. Such a high density would make it possible to effectively supply gas molecules, and this is a possible explanation of the mechanism of the nanobubble solution-induced promotion of the growth of plants and living beings.⁵⁷

Is pinning required for the stabilization of surface nanobubbles?

Some experiments and simulations have reported that pinning is responsible for the stabilization and flattening of surface nanobubbles,^{13,15,19,20,38,58} while several MD simulations, including our study, showed that the nanobubbles can be stable without pinning.^22,59 Zhang et al. predicted that the lifetime of surface nanobubbles would be extended if the liquid–gas interfacial tension could be significantly reduced by the high inner gas density.²⁶ However, the effect of the gas density on the liquid–gas interfacial tension is small as shown in Fig. 3(a). Then, what caused the nanobubbles to stabilize with no pinning sites in the MD simulations? We propose that it is because of the abnormally high gas oversaturation in water due to the closed condition. In our analysis, N₂ molecules have dissolved in water more than 1 kg m⁻³ (shown in ESI Fig. S8†), which is considerably supersaturated compared to the nitrogen solubility of 0.017 kg m⁻³ of ambient water. Because the simulation box is closed due to the periodic boundary in the horizontal direction and the walls in the vertical direction, the decrease in gas supersaturation via diffusion does not occur. Therefore, further dissolution is not induced even if the nanobubbles are at such a high inner pressure, resulting in their metastability. Indeed, it was reported that the nanobubble nucleation in LJ fluids can be controlled by adjusting the solubility through LJ parameters,⁶⁰ which supports our interpretation.

Such highly saturated and closed conditions are also present in the experimental systems. STXM observations of a confined water film reported that the water surrounding oxygen nanobubbles were highly saturated with oxygen gas.²⁴ In LPEM experiments, hydrogen nanobubbles are nucleated in a nano-confined liquid film by the radiolysis of water molecules. Because the radiolysis induced by electron beam irradiation continues to occur during observation, the gas supersaturation in the solution should be extremely high. In such cases, surface nanobubbles would be stable without pinning. Nag et al. recently found an unpinned nanobubble that was stable for over 1000 seconds in LPEM observation.²³ Shin et al. also observed the cross-sections of surface nanobubbles with diameters of 5–15 nm using the vended graphene liquid cell, and found that their contact angles were 90–120°,⁵⁶ which are consistent with the values of the unpinned nanobubbles on graphite shown in Fig. 4(b). In contrast, in open systems, the pinning and local gas oversaturation near solid–liquid interfaces will be indispensable for nanobubble stability. These indicate that the dynamic properties of surface nanobubbles can be different in open systems (e.g., typical systems in AFM and optical microscopy observations) and in closed systems (e.g., typical systems in STXM, LPEM, and MD simulations). Thus, the results obtained in different systems should be carefully compared.

Our analysis showed that it would be necessary to consider the force that causes the contact angles deviate from those predicted by Young's equation, i.e., pinning, to explain the properties of surface nanobubbles observed in open systems. The remaining question in experiments is why the pinning occurs on an atomically flat surface such as graphite. It is conceivable that the adsorbed layers of molecules at solid–liquid interfaces work as pinning sites.^{50–52,61–63} It has been found so far that the adsorbed layers have a thickness of monoatomic to several nanometers,⁵⁰ appear preferentially at the three-phase contact line of surface nanobubbles,⁶² and pin the gas phases nucleated on them.^50,52 However, it is still under debate whether the constituent molecules are gases or airborne contaminants.^63,64 Moreover, the pinning of surface nanobubbles by gas molecule-adsorbed layers and their nanoscopic thickness have never been reproduced by MD simulations. We expect that further experiments and simulations considering not only the gas molecule-adsorbed layers but also the concentration of ions in solutions and the electric charge at interfaces will provide valuable insights and address the origin of nanoscale pinning, which is essential for the development of micro- and nano-fluidic systems.

Conclusions

We quantified the interfacial tensions working on nitrogen surface nanobubbles at graphite–water interfaces through mechanical and thermodynamic routes. It was found that Young's equation holds even on the nanobubbles with a footprint diameter of a few nanometers. Furthermore, although a size-dependent change in their contact angle was observed, it was revealed that the effect of inner gas molecules on their interfacial tensions is negligible for the submicron-sized semispherical surface bubbles typically observed by AFM measurements. Therefore, the van der Waals interaction-induced gas molecular adsorption cannot explain the experimentally observed stability and flatness of surface nanobubbles. Our theoretical prediction also shows that surface nanobubbles with a footprint diameter of hundreds of nanometers cannot have an ultrahigh density of tens to hundreds kg m⁻³ only by the van der Waals interaction. The unexpected stability of surface nanobubbles without pinning in MD simulations, which contradicts the previous AFM results, was explained by extremely high gas supersaturation induced by the closed conditions, implying the necessity of the pinning for the stability of surface nanobubbles in open systems.

Methods

The detailed setup of the nanobubble simulation systems

All MD simulations of surface nanobubbles have been performed using the LAMMPS⁶⁵ package with a timestep of 1.0 fs. Periodic boundary conditions were employed in the x- and y-directions for a rectangular calculation cell of 16.17 × 3.19 × 13.0 nm³. The SPC/E model⁶⁶ was adopted for water and a nitrogen molecule was modeled by a monoatomic 12-6 LJ particle.⁶⁷ Because the Young's modulus of graphite is on the order of 10 GPa, which is considerably higher than the value that causes non-negligible surface deformation (on the order of 100 kPa⁶⁸), the influence of surface deformation on three-phase contact lines is not considered. We have also confirmed in our previous work that the vibration of atoms constituting a solid surface does not influence the validity of Young's equation.³⁰ Thus, the positions of the carbon atoms in the graphite were fixed onto the space. The gap between the layers in the graphite is 0.335 nm. Because the cutoff distance for the LJ interaction between the carbon and water and between the carbon and N₂ was set at 1.0 nm, further increase in the number of the graphene layers does not affect the simulation results. The inter-molecular interactions among the species were modeled by the 12-6 LJ and coulombic potential aswhere r_ij is the distance between the interaction sites i and j, and ε, σ, q_i and ε₀ denote the LJ energy, length parameter, and point electric charge at the site i and vacuum permittivity, respectively.

The top graphene basically used as a barostat (piston) was set to be hydrophilic by adjusting the LJ energy between the constituent carbon atom and fluid molecules (N₂ and water) so that the surface nanobubbles could nucleate only on the bottom graphite. The hydrophilicity of the top graphene does not affect the properties of surface nanobubbles because there are sufficiently large regions of bulk water between the top hydrophilic graphene and the nanobubbles (shown in ESI Fig. S9†). The cross LJ parameters except for the top graphene were determined using the Lorentz–Berthelot combining rules,⁶⁹ in which the LJ parameters of the carbon atom were assumed to be σ = 0.34 nm and ε = 0.06575 kcal mol⁻¹.⁷⁰ The LJ parameters and charge values used in the present study are listed in Table 1.

Table 1 Lennard-Jones (LJ) parameters and charge values adopted for the simulations

Pair	q (e)	σ ij (nm)	ε ij (kcal mol^−1)
O–O	−0.8476	0.3166	0.1555
H–H	0.4238	0	0
N2–N2		0.3698	0.1888
N2–H2O		0.3285	0.1924
Wall(C, top) – H2O		0.3367	0.2
Wall(C, top) – N2		0.3549	0.01
Wall(C, bottom) – H2O		0.3367	0.1015
Wall(C, bottom) – N2		0.3549	0.1114 × η (η = 1.0, 1.5)

---

The fluid temperature was kept at 300 K using the Nosé–Hoover thermostat. The coulombic interactions were calculated by the particle–particle-particle–mesh method.⁶⁵ A space three times longer than the length of the simulation domain was inserted in the z-direction to avoid the electrostatic interactions from the dummy cells in the wall-normal direction.⁷¹ The water bond length was constrained at 0.1 nm by the SHAKE algorithm.⁷² The system was divided into bins of 0.1 × l_y × 0.1 nm³ for the calculation of the density distribution, where l_y is the system size in the y-direction.

The details of the methods for the calculation of interfacial tensions

In addition to the simulation systems of surface nanobubbles, we also employed quasi-one-dimensional systems to calculate the corresponding interfacial tensions. We note that the basic setting of molecular models, potential parameters, and the control temperature is the same as that of the nanobubble systems.

For the calculation of the liquid–gas interfacial tension, we used two equilibrium simulation systems with a flat-water film: with two liquid–vapor interfaces and with two liquid–gas interfaces (shown in ESI Fig. S10(a)†). Liquid–vapor or liquid–gas interfacial tensions γ_LV and γ_LG of the systems can be expressed as

where P_xx, P_yy and P_zz are the diagonal components of the average pressure tensor, and l_z is the system size in the z-direction. Note that the right-hand side of eqn (6) is divided by two considering that the systems have two interfaces.

Next, we employed the system in which N₂ molecules were sandwiched between the top graphene and the bottom graphite (shown in ESI Fig. S10(b)†). The solid–gas interfacial tension relative to the solid–vacuum interfacial tension defined as γ_SG − γ_S0 was calculated by the Bakker's equation extended for the solid–gas interface³⁰ given by

where τ_t(z) and τ_blk were the fluid stress tangential to the solid–gas interface and bulk stress, respectively. To calculate the stress distributions with sufficient accuracy, the system was divided into small bins of l_x × l_y × 0.01 nm³, where l_x is the system size in the x-direction. The tangential stress τ_t = τ_xx = τ_yy is a function of the distance z and it starts to have non-zero values from the solid–gas interface at z = z_sg, whereas the bulk stress τ_blk is a constant equivalent to the bulk gas pressure with inverted sign (−P_gas).

Finally, we calculated the work of adhesion at solid–liquid interfaces W_SL as the minimum work needed to separate the liquid from the solid surface per area, which corresponds to the interfacial free energy. Specifically, we employed the equilibrium MD simulation systems of a water–graphite interface shown in ESI Fig. S11(a)† to quasi-statically strip the liquid off from the solid surface using a “phantom wall” under the constant temperature conditions. This is an implementation of the TI method for the calculation of the solid–liquid interfacial free energy.^34,43

In the systems of the TI method, a phantom wall that interacts only with water molecules via the repulsive range of the 12-6 LJ model is at first placed sufficiently below the water film, z_w = z₁. For the present case, z₁ is positioned exactly at the middle of the two topmost graphene layers, i.e., approximately 0.4 nm below the lowest water molecules. Then, by quasi-statically moving the phantom wall step by step up to the place where the water film is completely stripped from the graphite surface (z_w = z₀), the work of adhesion was extracted through the TI along a reversible path. The LJ energy ε and length parameter σ of the phantom wall were 0.316 nm and 0.16 kcal mol⁻¹, respectively. The cutoff distance is set at 0.3546 nm to make the interaction repulsive only. The position of the phantom wall in the z-direction z_w(λ) was expressed with a coupling parameter λ as follows:

In our case, we calculated multiple equilibrium systems of graphene–water interfaces with discrete coupling parameters λ from 0 to 1.0. The solid–liquid interface at λ = 1.0 is separated to the solid–vacuum and liquid–phantom wall interfaces by changing the coupling parameter to λ = 0 as shown in the left and right panels in ESI Fig. S11(a).† The phantom wall was displaced by 0.01 nm in the z-direction until the force exerted by the wall was confirmed to be close to zero, i.e., the interaction between the liquid film and the graphene surface had become neglectable. A total of 46 simulations were performed, with a displacement of 0.45 nm, i.e., z₀ corresponds to the phantom wall positioned 0.3 nm above the topmost graphene layer. The variation of the force exerted by the wall on the liquid as a function of z_w is shown in ESI Fig. S11(b).†

Because the number of molecules N, volume V, and temperature T are constant in the systems, the difference in the Helmholtz free energy ΔF between the systems at z_w = z₁ and z_w = z₀ under a constant NVT is expressed by the interfacial free energies as follows:

where A denotes the surface area and the phantom wall is denoted by the subscript “W”, γ_LW and γ_S0 are the liquid–phantom wall and solid–vacuum interfacial energies per unit area. Note that γ_LW was substituted by the liquid–vapor interfacial tension γ_LV because the LJ parameters were set so that only the repulsive force worked between the water molecules and the phantom wall.

Using the NVT canonical ensemble, the difference in the Helmholtz free energy ΔF in eqn (9) is calculated through the TI as

where H is the Hamiltonian, i.e., the internal energy of the system, and f_w(z_w) is the sum of forces exerted on water molecules by the phantom wall. The ensemble average was substituted by the time average and was denoted by the angular brackets. By substituting eqn (10) and γ_LV calculated using eqn (6)–(9), γ_SL − γ_S0 was obtained to be −17.3 ± 1.0 × 10⁻³ N m⁻¹. Although the N₂ molecules were adsorbed onto the graphite–water interfaces of nanobubble systems (Fig. 2), they did not give considerable influence on the contact angles, i.e., the contact angles estimated by Young's equation did not deviate from the apparent contact angles, as shown in Fig. 3(b). This may be because the density of N₂ molecule-adsorbed layers at graphite–water interfaces (up to 50 kg m⁻³) were considerably lower than the adsorbed layer of water molecules (∼1800 kg m⁻³).

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Prof. Koji Takahashi at Kyushu University, Prof. Takeshi Omori at Osaka City University, and Haruki Oga and Minori Shintaku at Osaka University for fruitful discussion. This work was partially supported by JST CREST Grant No. JPMJCR18I1, JSPS KAKENHI Grant No. JP18K03978 and JP21K20405, and the Grant-in-Aid for JSPS Research Fellow No. JP20J01307 and JP20J20251.

References

1. D. Lohse and X. Zhang, Rev. Mod. Phys., 2015, 87, 981–1035.
2. J. L. Parker, P. M. Claesson and P. Attard, J. Phys. Chem., 1994, 98, 8468–8480.
3. E. E. Meyer, K. J. Rosenberg and J. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 15739–15746.
4. P. G. de Gennes, Langmuir, 2002, 18, 3413–3414.
5. Y. Tomo, A. Askounis, T. Ikuta, Y. Takata, K. Sefiane and K. Takahashi, Nano Lett., 2018, 18, 1869–1874.
6. L. Gao, G. X. Ni, Y. Liu, B. Liu, A. H. Castro Neto and K. P. Loh, Nature, 2014, 505, 190–194.
7. Y. Nam and Y. S. Ju, Appl. Phys. Lett., 2008, 93, 103115.
8. Y. Gao, J. Li, H. C. Shum and H. Chen, Langmuir, 2016, 32, 4815–4819.
9. A. Aikawa, A. Kioka, M. Nakagawa and S. Anzai, Geothermics, 2021, 89, 101962.
10. X. H. Zhang, A. Quinn and W. A. Ducker, Langmuir, 2008, 24, 4756–4764.
11. S. R. German, X. Wu, H. An, V. S. J. Craig, T. L. Mega and X. Zhang, ACS Nano, 2014, 8, 6193–6201.
12. H. An, G. Liu, R. Atkin and V. S. J. Craig, ACS Nano, 2015, 9, 7596–7607.
13. H. Teshima, T. Nishiyama and K. Takahashi, J. Chem. Phys., 2017, 146, 014708.
14. X. Zhang, H. Lhuissier, C. Sun and D. Lohse, Phys. Rev. Lett., 2014, 112, 144503.
15. A. Brotchie and X. H. Zhang, Soft Matter, 2011, 7, 265–269.
16. B. H. Tan, H. An and C. D. Ohl, Curr. Opin. Colloid Interface Sci., 2021, 53, 101428.
17. D. Lohse and X. Zhang, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 031003.
18. B. H. Tan, H. An and C. D. Ohl, Phys. Rev. Lett., 2019, 122, 134502.
19. Y. Liu and X. Zhang, J. Chem. Phys., 2013, 138, 014706.
20. X. Zhang, D. Y. C. Chan, D. Wang and N. Maeda, Langmuir, 2013, 29, 1017–1023.
21. B. H. Tan, H. An and C. D. Ohl, Phys. Rev. Lett., 2017, 118, 1–5.
22. Y. X. Chen, Y. L. Chen and T. H. Yen, Langmuir, 2018, 34, 15360–15369.
23. S. Nag, Y. Tomo, H. Teshima, K. Takahashi and M. Kohno, Phys. Chem. Chem. Phys., 2021, 23, 24652–24660.
24. L. Zhou, X. Wang, H. J. Shin, J. Wang, R. Tai, X. Zhang, H. Fang, W. Xiao, L. Wang, C. Wang, X. Gao, J. Hu and L. Zhang, J. Am. Chem. Soc., 2020, 142, 5583–5593.
25. X. H. Zhang, A. Khan and W. A. Ducker, Phys. Rev. Lett., 2007, 98, 136101.
26. L. J. Zhang, H. Chen, Z. X. Li, H. P. Fang and J. Hu, Sci. China, Ser. G: Phys., Mech. Astron., 2008, 51, 219–224.
27. N. D. Petsev, L. G. Leal and M. S. Shell, Phys. Rev. Lett., 2020, 125, 146101.
28. S. R. German, M. A. Edwards, Q. Chen and H. S. White, Nano Lett., 2016, 16, 6691–6694.
29. P. Attard, 2015, arXiv, arXiv:1505.02217, 1–16.
30. Y. Yamaguchi, H. Kusudo, D. Surblys, T. Omori and G. Kikugawa, J. Chem. Phys., 2019, 150, 044701.
31. G. Bakker, Kapillarität und Oberflächenspannung, Wien-Harms, 1928, vol. 6.
32. S. M. Thompson, K. E. Gubbins, J. P. R. B. Walton, R. A. R. Chantry and J. S. Rowlinson, J. Chem. Phys., 1984, 81, 530–542.
33. H. Kusudo, T. Omori and Y. Yamaguchi, J. Chem. Phys., 2019, 151, 154501.
34. C. Bistafa, D. Surblys, H. Kusudo and Y. Yamaguchi, J. Chem. Phys., 2021, 155, 064703.
35. D. Surblys, F. Leroy, Y. Yamaguchi and F. Müller-Plathe, J. Chem. Phys., 2018, 148, 134707.
36. S. Nishida, D. Surblys, Y. Yamaguchi, K. Kuroda, M. Kagawa, T. Nakajima and H. Fujimura, J. Chem. Phys., 2014, 140, 074707.
37. D. Surblys, Y. Yamaguchi, K. Kuroda, M. Kagawa, T. Nakajima and H. Fujimura, J. Chem. Phys., 2014, 140, 034505.
38. Y. Liu and X. Zhang, J. Chem. Phys., 2014, 141, 134702.
39. B. R. Novak, E. J. Maginn and M. J. McCready, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 085413.
40. S. Wang, L. Zhou, X. Wang, C. Wang, Y. Dong, Y. Zhang, Y. Gao, L. Zhang and J. Hu, Langmuir, 2019, 35, 2498–2505.
41. M. Masuduzzaman and B. H. Kim, Microfluid. Nanofluid., 2021, 25, 1–14.
42. D. Kim, E. Kim, S. Park, S. Kim, B. K. Min, H. J. Yoon, K. Kwak and M. Cho, Chem, 2021, 7, 1602–1614.
43. F. Leroy and F. Müller-Plathe, J. Chem. Phys., 2010, 133, 044110.
44. A. E. Ismail, G. S. Grest and M. J. Stevens, J. Chem. Phys., 2006, 125, 014702.
45. N. Giovambattista, A. B. Almeida, A. M. Alencar and S. V. Buldyrev, J. Phys. Chem. C, 2016, 120, 1597–1608.
46. F. Chen and P. E. Smith, J. Chem. Phys., 2007, 126, 221101.
47. Y. H. Ma, W. Sun, M. Bhandarkar and J. Wang, Sep. Technol., 1991, 1, 90–98.
48. H. Swenson and N. P. Stadie, Langmuir, 2019, 35, 5409–5426.
49. B. von Szyszkowski, Z. Phys. Chem., 1908, 64U, 385–414.
50. H. Teshima, Y. Takata and K. Takahashi, Appl. Phys. Lett., 2019, 115, 071603.
51. H. Teshima, Q. Y. Li, Y. Takata and K. Takahashi, Phys. Chem. Chem. Phys., 2020, 22, 13629–13636.
52. H. Teshima, N. Nakamura, Q. Y. Li, Y. Takata and K. Takahashi, RSC Adv., 2020, 10, 44854–44859.
53. W. H. Hsu and I. S. Hwang, Chem. Sci., 2021, 12, 2635–2645.
54. S. Nag, Y. Tomo, K. Takahashi and M. Kohno, Langmuir, 2021, 37, 874–881.
55. Y. Tomo, K. Takahashi, T. Nishiyama, T. Ikuta and Y. Takata, Int. J. Heat Mass Transfer, 2017, 108, 1460–1465.
56. D. Shin, J. B. Park, Y. J. Kim, S. J. Kim, J. H. Kang, B. Lee, S. P. Cho, B. H. Hong and K. S. Novoselov, Nat. Commun., 2015, 6, 1–6.
57. K. Ebina, K. Shi, M. Hirao, J. Hashimoto, Y. Kawato, S. Kaneshiro, T. Morimoto, K. Koizumi and H. Yoshikawa, PLoS One, 2013, 8, 1–7.
58. S. Maheshwari, M. Van Der Hoef, X. Zhang and D. Lohse, Langmuir, 2016, 32, 11116–11122.
59. T. H. Yen, C. H. Lin and Y. L. Chen, Langmuir, 2021, 37, 2759–2770.
60. J. H. Weijs, J. H. Snoeijer and D. Lohse, Phys. Rev. Lett., 2012, 108, 104501.
61. I.-S. Hwang, Y.-H. Lu, C.-W. Yang, C.-K. Fang and H.-C. Ko, Sci. Rep., 2014, 4, 1–8.
62. C. K. Fang, H. C. Ko, C. W. Yang, Y. H. Lu and I. S. Hwang, Sci. Rep., 2016, 6, 1–10.
63. I. Schlesinger and U. Sivan, J. Am. Chem. Soc., 2018, 140, 10473–10481.
64. M. R. Uhlig, S. Benaglia, R. Thakkar, J. Comer and R. Garcia, Nanoscale, 2021, 13, 5275–5283.
65. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
66. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
67. J. O. Hirschfelder, C. F. Curtiss and R. B. Bird, Molecular Theory of Gases and Liquids, Wiley, 1964.
68. R. Pericet-Cámara, A. Best, H. J. Butt and E. Bonaccurso, Langmuir, 2008, 24, 10565–10568.
69. G. H. Hudson and J. McCoubrey, Trans. Faraday Soc., 1959, 56, 761–766.
70. S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472–6486.
71. I. C. Yeh and M. L. Berkowitz, J. Chem. Phys., 1999, 111, 3155–3162.
72. J. P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23, 327–341.

---

Footnote

† Electronic supplementary information (ESI) available: Normal stress profiles under surface nanobubbles; peak density inside the N₂ molecule-adsorbed layers; breakdown of Young's equation for the case of N_N2 = 100 and η = 1.0; calculation of the contact angles of semispherical submicron-sized bubbles; bulk N₂ density inside the nanobubbles estimated using eqn (4); density distributions of N₂ molecules in the nanobubble systems; density distributions of water molecules in the nanobubble systems; and systems used for the calculation of interfacial tensions. See DOI: 10.1039/d1nr07428h

---