Leakiness of Pinned Neighboring Surface Nanobubbles Induced by Strong Gas-Surface Interaction.

The stability of two neighboring surface nanobubbles on a chemically heterogeneous surface is studied by molecular dynamics (MD) simulations of binary mixtures consisting of Lennard-Jones (LJ) particles. A diffusion equation-based stability analysis suggests that two nanobubbles sitting next to each other remain stable, provided the contact line is pinned, and that their radii of curvature are equal. However, many experimental observations seem to suggest some long-term kind of ripening or shrinking of the surface nanobubbles. In our MD simulations we find that the growth/dissolution of the nanobubbles can occur due to the transfer of gas particles from one nanobubble to another along the solid substrate. That is, if the interaction between the gas and the solid is strong enough, the solid-liquid interface can allow for the existence of a "tunnel" which connects the liquid-gas interfaces of the two nanobubbles to destabilize the system. The crucial role of the gas-solid interaction energy is a nanoscopic element that hitherto has not been considered in any macroscopic theory of surface nanobubbles and may help to explain experimental observations of the long-term ripening.

Understanding the stability and the dynamics of surface nanobubbles,[1−6] gaseous nanoscopic entities attached to hydrophobic surfaces, is relevant from both a fundamental and an applicational point of view. A single surface nanobubble will remain stable, provided a sufficient gas oversaturation in the bulk liquid and pinning of the three phase contact line.[7−13] Lohse and Zhang[7] showed that the equilibrium contact angle is not dictated by Young’s law, but by the expression:where L is the fixed (due to pinning) diameter of the footprint (assumed to be spherical) and L = 4γ/P0 is the capillary length scale, γ is the surface tension, and P0 is the ambient pressure. Furthermore, in expression 1, is the gas oversaturation, C∞ is gas concentration in the bulk liquid, and C is the saturation concentration at pressure P0. The stability conditions and expression 1 for the contact angle of a single surface nanobubble at equilibrium were also confirmed by MD simulations[14] and experiments.[15] In most practical scenarios, nanobubbles are however not single, but exist in the neighborhood of other nanobubbles.[1,16−20] A recent theoretical study[21] shows that two nanobubbles remain stable against Ostwald ripening (mathematically, a linear instability), due to the pinning of the contact line. In that calculation it was assumed that the distance between the two nanobubbles was large compared to the lateral dimension of the nanobubbles. However, as confirmed in recent numerical simulations based on the advection-diffusion equation,[22] nanobubbles can also remain stable even when the distance between them is comparable to the lateral length, provided that the radii of curvature of the nanobubbles are equal. In that case, the stability is the consequence of the absence of any concentration gradient of gas molecules in the liquid. However, interestingly, while for pinning, oversaturation ζ > 0, and equal radius of curvature R = L/sin θ = L/ζ, the macroscopic theory predicts linear stability of the two-nanobubble system, there have been many experimental observations which show some sort of ripening/shrinking of surface nanobubble populations.[11,18,23−25] This motivated us to look into the problem at a microscopic scale using molecular dynamics (MD) simulations, which is perfectly suited for the length scales and time scales involved in the problem.

In this paper, we have studied the dynamics of two nanobubbles sitting next to each other with their contact line pinned on a chemically heterogeneous surface. Our simulations suggest a possible explanation for the experimental observations of the ripening/shrinking of surface nanobubbles, where we observe that under certain conditions the system allows for transfer of gas particles from one nanobubble to a neighboring one along the solid substrate. We find that the ripening of surface nanobubbles strongly depends on the interaction energy between the gas molecules and the solid substrate relative to the interaction energy between the liquid molecules and the solid substrate.

Results and Discussion

We performed MD simulations of two surface nanobubbles sitting next to each other on a chemically heterogeneous surface, so that the contact line is pinned (see Figure ). At time t = 0, both nanobubbles have the same radius of curvature, lateral length, contact angle, and number of gas particles inside the bubble. We performed the simulations at constant pressure which means the solubility of gas particles is fixed and can be manipulated by varying the system pressure. Figure a shows that the number of gas particles in both nanobubbles remains constant with time. We found that this is also true for the contact angle, the lateral length, and the radii of both nanobubbles. In other words, the system seems to be in steady state (also see Supplementary Movie S1). Note that the time for which we run the simulation (∼550 ns) is longer than the diffusion time scale, τ ∼ L2/D = 400 ns, where L is the largest dimension of the box, and D is the diffusivity of gas particles in the bulk liquid.

Figure 1

A typical simulation box which consists of four kind of particles. Red particles (S) form the solid surface, yellow particles (SP) form pinning sites, blue particles (L), which are predominantly in the liquid phase, and green particles (G), which are predominantly in the gas phase. The nanobubbles shown in the figure are cylindrical in shape where the length of the cylinder is in the x-direction. Periodic boundaries are employed in all three directions.

Figure 2

Variation of number of gas particles, θ, L, and R of two nanobubbles with time for (a) ϵ/ϵ = 0.50 and (b) ϵ/ϵ = 1.11, and C = 17.3 × 10–5. Inset of top left plot shows the number of gas particles inside two nanobubbles up to 560 ns which confirms the steady state. For ϵ/ϵ = 1.11 and at t = 175 ns, one can clearly notice the “stick-jump” dissolution mode for the pinned case from the variation of θ, L, and R.

The stability of our system is consistent with the theoretical arguments[21] which suggest that the two-nanobubble state should be stable due to the contact line pinning. However, many experimental observations[11,18,23−25] show the ripening of surface nanobubbles in which one nanobubble grows at the expense of others. One can argue that from the free-energy point of view, it is always favorable to form one bubble instead of two because of the minimization of interface area. Yet to reach the energetically favorable state, one first has to overcome the contact line pinning barrier, which requires a (strong) concentration gradient, which is absent in the case of identical radii of curvature. Namely, according to Henry’s law, the concentration of the gas at the interface of a bubble is directly proportional to the pressure inside the bubble, which is a function of the radius of the bubble as given by Laplace’s law. So from the macroscopic point of view, no concentration gradient is present in the system that can facilitate the transfer of particles from one nanobubble to another.

Since there is no net flux of particles through the bulk liquid, the conclusion is that the only way gas molecules can get transferred from one nanobubble to another is along the solid substrate at the solid–liquid interface. Such a mechanism would require a film layer of gas particles on the solid surface, which can be achieved by increasing the interaction energy between gas and solid particles, relative to the other interactions. To this end, we have performed MD simulations in which we increased the ratio of the solid–gas to solid–liquid interaction (ϵ/ϵ) from 0.50 to 1.11 (for complete list of all the interactions, see Table ). We simulated the system for an increased ratio of ϵ/ϵ and observed that indeed after a few nanoseconds, one nanobubble grows at the expense of another, as can be seen in Figures b and 3 (also see Supplementary Movie S2). Figure shows the number of gas particles as a function of time in the two nanobubbles for various system pressures, or equivalently, saturation concentrations C. Note that C is measured exactly in the same manner as in ref (14). The results shown in Figure prove that the conclusion that one bubble absorbs the gas content of the other is fairly robust and occurs for a wide range of saturation concentrations.

Table 1

Value of Various LJ Parameters Used in the MD Simulations

i–j	σij, nm	ϵij, kJ/mol
S-L	0.34	1.8 and 2.0
SP-L	0.34	1.5
S-G	0.40	1.0 and 2.0
SP-G	0.40	5.0
L-G	0.40	1.55
G-G	0.46	0.8
L-L	0.34	3.0

Figure 3

Number of gas particles inside the two nanobubbles as a function of time when the contact lines of both bubbles are pinned for ϵ/ϵ = 1.11. Different lines indicate the simulation at various solubilities which is indicated by the color bar. For all solubilities that we studied, we find ripening of one bubble at the expense of the other. Note that there is no systematic dependence of the time of the onset of the instability on the solubility, and thus the color coding mainly serves to identify the two bubbles (one growing and one shrinking) which belong to one simulation.

In Figure b we focus on a system with one particular solubility (C = 17.3 × 10–5). It can be clearly seen that even while the stability criteria are met (the contact lines of both nanobubbles should remain pinned and the radius of both nanobubbles should be equal),[21] the system becomes unstable due to the gas leakage along the surface, as we will demonstrate later. At t = 175 ns, a sudden jump is observed in the values of the contact angle, the lateral length, and the radius of both nanobubbles. This sudden jump in the values shows the “stick-jump” motion of the contact line which was also observed during the dissolution of a single surface nanobubble[14] and the dissolution of a microdrop in another liquid.[26,27] The “stick-jump” motion of the contact line is a clear indication of contact line pinning, where the contact line prefers to stay on the chemical heterogeneities and shifts from one pinning site to another during the dissolution of a nanobubble.

As mentioned earlier, increasing the ratio of ϵ/ϵ leads to a formation of a gas layer near the solid surface. To show this more clearly, we plotted the normalized concentration of gas particles within a 5 particle diameter from the solid substrate, and indeed the concentration of gas particles has increased by more than 10 times, for an increased ratio ϵ/ϵ by a factor of around 2 times (see Figure ). This gas layer links the interface of the two nanobubbles like a tunnel that overcomes the contact line pinning barrier and aids the system to reach the energetically favorable state. We calculated the net movement of particles between two nanobubbles for ϵ/ϵ = 1.11 and for 0.50. Figure a,b shows the net movement of particles in the bulk liquid, along the surface and the total net movement of particles. It is clear that almost all particles are transferred along the substrate from one nanobubble to another when ϵ/ϵ = 1.11 (see Figure b). In contrast, Figure a shows that the net flux is almost zero when ϵ/ϵ = 0.50. Inset of Figure b shows that during the growth period, particles transferred at an exponential rate ∝ e (where τ is the time constant of the growth rate) from one nanobubble to other, which demonstrates that the system is highly unstable. The exponential dependence of rate on time shows the linear instability.

Figure 4

Variation of normalized concentration of gas particles within 5 particle diameter of distance from the solid surface for C = 17.3 × 10–5. Here the concentration is normalized with the concentration of gas particles in the bulk liquid far away from the solid surface.

Figure 5

Cumulative number of gas particles moving from one nanobubble to another as a function of time, when (a) ϵ/ϵ = 0.5 and (b) ϵ/ϵ = 1.11. Inset of (b) shows the cumulative transfer of gas particles during the growth period on a semilog scale. It shows the exponential growth of the nanobubble after destabilization with a time constant τ = 28.9 ns. The red line indicates the number of particles that are moving along the surface (that is, below a height of 5σ), while the green line indicates the number of particles moving in the bulk (that is, above a height of 5σ, where σ is the diameter of liquid particles).

In Figure we show the flux of the particles leaving the dissolving nanobubble near the contact line as a function of the contact angle, where it is clearly observed that the flux of the particles increases with decreasing contact angle. This increase in the flux of the gas particles from the dissolving nanobubble further demonstrates the instability of the system.

Figure 6

Variation of flux of particles near the contact line of dissolving nanobubble as a function of contact angle. Once the contact angle of the pinned nanobubble becomes smaller, the flux along the surface increases, signaling instability of the two bubble system.

Note that while we observe in Figure that the onset of growth or dissolution differs for the various gas solubilities, this is mainly due to the statistical fluctuations, that is, no specific trend with gas solubility would be observed. For different realizations of the same numerical experiment, but at different initial conditions, we also find that the moment of the onset t of the instability is fluctuating (see Figure ). While it would be interesting to measure the distribution of t, and its dependence on the solubility, this is far beyond our computational resources. In order to get statistically significant data one would have to perform thousands of simulations at different initial conditions for each solubility and then average, which would be extremely expensive computationally.

Figure 7

Variation of the number of gas particles inside the two nanobubbles with time. Different line colors indicate the simulation at various initial conditions at a saturation concentration C = 15.4 × 10–5 and ϵ/ϵ = 1.11.

Conclusions

To summarize, we showed that a two-nanobubble system can remain stable when the contact line is pinned, consistent with the theory, as long as the gas–solid interaction (ϵ) is low compared to the liquid–solid interaction (ϵ). We find that upon increasing ϵ relative to ϵ, a film layer of gas particles at the surface forms a channel which allows gas particles to transfer along the solid surface at the solid–liquid interface and not via the bulk liquid. Our simulations thus showed that the nanoscopic interaction with solid surface can play a crucial role in the stability of two-nanobubble systems, which hitherto has not been considered in the macroscopic theories of surface nanobubbles. We showed that the system is unstable for ϵ/ϵ = 1.1 and stable for ϵ/ϵ = 0.5. It will be interesting to know the rough estimate of ϵ/ϵ at which the transition occurs. Figure a shows the number of gas particles in both the nanobubbles as a function of time for various values of ϵ/ϵ. It can be observed that the curves diverge when ϵ/ϵ > 0.8, while the number of gas particles in each nanobubble remains constant when ϵ/ϵ < 0.8, at least for the time that we ran the simulation. We also plotted the concentration of gas particles near the solid surface (see Figure b) which grows exponentially with the increase in ϵ/ϵ. Although the gas layer near the solid surface enhances the transfer of gas particles from one nanobubble to another, the onset of this transfer varies quite a lot for different realizations (see Figures and 7). This stochastic variability refrains us from claiming ϵ/ϵ = 0.8 as an exact boundary for the instability, as for some microscopic conditions, the system may show stability or instability for the value of ϵ/ϵ close to 0.8. Nevertheless ϵ/ϵ ≈ 0.8 can still be considered as a rough estimate for the boundary of instability based on our simulations. It would require an extensive number of simulations performed at various initial conditions, if one wanted to obtain better statistical estimates of the exact value of ϵ/ϵ beyond which the system is unstable. This was not feasible for us, as these MD simulations were computationally very expensive. In real systems, surface heterogeneities along the “channel” between the two bubbles will moreover modify the onset value.

Figure 8

(a) Number of gas particles inside the two nanobubbles with time for various values of ϵ/ϵ. (b) Normalized concentration of gas particles within 5 particle diameter of distance from the solid surface as a function of ϵ/ϵ. The inset shows the same data on a log linear scale. It can be seen that the concentration of gas particles near the surface increases exponentially with ϵ/ϵ. The value of ϵ/ϵ ∼ 0.8 is estimated as a rough boundary between the region of a stable and unstable two nanobubble system.

For the sake of connecting our simulations with real systems, we calculated the value of ϵ/ϵ for some real gas–solid–liquid combinations with known interaction energy parameters and identified the possible combinations where one can expect stable and unstable two-nanobubble systems. For example, two nanobubbles of krypton gas in water on a graphite surface will give the value of ϵ/ϵ as 0.98, which means this combination will give an unstable two-nanobubble system. However, if we replace krypton with neon, then the value of ϵ/ϵ becomes 0.46, which makes the two-nanobubble system stable, according to our simulation results. The interaction energy parameters used for the comparison are taken from refs (28 and 29).

Simulation Methodology

Molecular dynamics (MD) simulations were performed to simulate the nanobubbles on a solid substrate for which we used the open source code GROMACS.[30] We used four types of particles or molecules in our simulations: two types of solid particles (S and SP), which remain fixed in a fcc lattice during the whole simulation, and two types of moving particles, liquid (L) and gas (G). The SP particles form the pinning sites and have different interaction strengths toward the two types of moving particles. The L particles form a bulk liquid phase (and hence we refer to these as “liquid particles”) as the system temperature and pressure are below the critical point of L particles, whereas the G particles form a bulk gaseous phase (to which we refer as “gas particles”) because the critical point for G particles is much below the thermodynamic conditions at which we performed our simulations.

The interaction between the particles is described by a Lennard-Jones potential:where ϵ is the interaction strength between particles i and j, and σ is the characteristic size of the particles. The potential is truncated at a relatively large cutoff radius (r) of 5σ. The time step for updating the particle velocities and positions was set at , where m is mass of the liquid particles, and ϵ is the Lennard-Jones interaction parameter for the liquid particles. The time step was chosen such that its value is sufficiently less than the shortest time scale available in the system.[31] Periodic boundary conditions were employed in all three directions, which implies that the same solid substrate is also present above the liquid layer.

The simulations were performed in an NPT ensemble where the temperature was fixed at 300 K, which is below the critical point for the Lennard-Jones parameters (σ, ϵ) that we set for the liquid particles. Semi-isotropic pressure coupling was used for maintaining constant pressure, which means that the simulation box can expand or contract only in the z-direction to keep the pressure constant. This was done to avoid the creation of gaps along the solid surface boundaries in the x and y directions. The complete set of Lennard-Jones parameters that we used in our simulations is given in Table . The typical system size was 5.6 × 60 × 22 nm3 in x, y, and z directions, respectively, where the length of the z dimension was changed during the simulation to maintain the pressure constant.

In the initial configuration, gas and liquid particles were arranged in an fcc lattice above the solid substrate. Initially, the liquid near the surface was highly oversaturated with gas particles in order to aid the bubble nucleation on the surface. This decreased the equilibration time, which is around 108 time steps (around ∼18 ns). For simulations at different pressures, we used the final configuration of the previous simulation as an initial configuration to save computation time. In Figure we show a typical profile of two nanobubbles sitting next to each other on a chemically heterogeneous surface.

The time-dependent average density field of liquid particles was then calculated, correcting for the center of mass motion in the lateral direction. Quantities like the radius of curvature of the bubble and the contact angle were obtained by fitting a circle to the isodensity contour of 0.5 of the normalized density field, ρ*(r), defined as , where ρ and ρ are the bulk vapor and liquid density, respectively. Since the liquid very near to the solid substrate is subject to layering, we excluded the density field in the range of 2σ from the substrate for the circular cap fitting.