Diffusive interaction of multiple surface nanobubbles and nanodroplets: shrinkage, growth, and coarsening

Surface nanobubbles are nanoscopic spherical-cap shaped gaseous domains on immersed substrates which are stable, even for days. After the stability of a {\it single} surface nanobubble has been theoretically explained, i.e. contact line pinning and gas oversaturation are required to stabilize them against diffusive dissolution [Lohse and Zhang, Phys.\ Rev.\ E 91, 031003 (R) (2015)], here we focus on the {\it collective} diffusive interaction of {\it multiple} nanobubbles. For that purpose we develop a finite difference scheme for the diffusion equation with the appropriate boundary conditions and with the immersed boundary method used to represent the growing or shrinking bubbles. After validation of the scheme against the exact results of Epstein and Plesset for a bulk bubble [J. Chem. Phys. 18, 1505 (1950)] and of Lohse and Zhang for a surface bubble, the framework of these simulations is used to describe the coarsening process of competitively growing nanobubbles. The coarsening process for such diffusively interacting nanobubbles slows down with advancing time and thus increasing bubble distance. The present results for surface nanobubbles are also applicable for immersed surface nanodroplets, for which better controlled experimental results of the coarsening process exist.


I. INTRODUCTION
Surface nanobubbles [1] -nanoscopic gaseous domains on immersed surfaces -were first speculated to exist about 20 years ago [2] and later found in atomic force microscopy (AFM) images [3][4][5]. While their long-time existence (often days) was first considered as puzzling [6] due to the supposedly large internal Laplace pressure, which should squeeze them out, it is now theoretically understood that they are stable thanks to a stable balance between the Laplace pressure inside the nanobubble and the gas overpressure from outside, which is enabled by pinning of the contact line [1,[7][8][9][10]. The equilibrium angle θ e (see figure   1 for a sketch of the surface nanobubble and the used notation) is determined by the gas oversaturation ζ = c ∞ /c s − 1, where c ∞ is the concentration far away and c s the solubility, and the contact diameter L by [10] sin θ e = ζ L L c , where L c = 4σ/P 0 = 2.84µm for air in water under ambient pressure P 0 = 1 bar and with its surface tension σ = 0.072N/m. Note that we have assumed a spherical-cap shape, which is well-justified theoretically and experimentally. The experimental confirmation of equation (1) through AFM experiments is difficult for various reasons [1], but it was confirmed in molecular dynamics (MD) simulations [11].
In this paper we will first add further numerical confirmation of the theory of Ref. [10] by directly solving the diffusion equation around a surface nanobubble, together with the appropriate boundary conditions, namely c ∞ far away from the bubble, no gas flux through the substrate, and a gas concentration given by Henry's law at the bubble-liquid interface, finding perfect agreement for the equilibrium contact angle θ e (Eq. (1)) (section III). Before, in section II, we will introduce the employed numerical method, namely a finite difference scheme coupled to an immersed boundary method [12][13][14].
Note that Eq. (1) implies that the Young-Laplace relation, which determines the contact angle on a macroscopic scale due to the mutual interfacial tensions, is irrelevant on the microscopic scale of the nanobubbles. This is in agreement with various experimental observations (see e.g. [1,15]) that the microscopic contact angle is constant and independent of the substrate and thus different from the macroscopic contact angle. According to Eq.
(1), the crossover from macroscopic to microscopic bubbles occurs at the length scale L c /ζ, below which the bubbles are small enough so that their Laplace pressure is large enough to counteract the gas influx by oversaturation.
The main focus of the present paper will however be on multiple surface bubbles which are diffusively interacting [7,10,16,17]. In general, no analytical solution is possible for this case. An exception is the case of two diffusively interacting surface bubbles far away from each other, i.e., with a distance d much larger than their surface contact diameter L.
For that case Dollet and Lohse [18] succeeded to analytically show that the pinning of the surface bubbles not only stabilizes each bubble against dissolution or growth, but that it also stabilizes the pair of surface bubbles against Ostwald ripening [19], i.e., the shrinkage of a bubble with smaller radius of curvature (corresponding to large Laplace pressure) to the benefit of a neighboring bubble with larger radius of curvature. Here we will numerically show that this stabilization of a pair of surface bubbles through pinning holds in general, i.e., is not limited to bubbles far away from each other. We will also show that the lack of pinning leads to Ostwald ripening (section III).
In section V we will extend the calculation to many surface nanobubbles in a row, studying their coarsening process. The coarsening of nanobubbles in principle can happen via Ostwald ripening or via coalescence. In Ref. [20] the analogous coarsening process of nanodroplets growing in an oversaturated solution was experimentally studied. There the nanodroplets also effectively sit in a row, namely at the rim of a spherical lense, and our assumption of periodic boundary conditions for the bubbles is justified. In that reference [20] it was speculated that the coarsening mainly happens via Ostwald ripening. Here within our model we will show under what conditions this indeed can be the case. We will moreover study the dynamics of the coarsening process and show that it slows down with advancing time and thus increasing distance between the bubbles, similar to other coarsening processes [21].
As mentioned above, our numerical scheme can not only be applied to diffusively interacting nanobubbles in a liquid, but equally well to diffusively interacting droplets in a liquid (see e.g. our own work on this subject, Refs. [20,22,23]) or in a gas, e.g., as they emerge in dew formation [24][25][26][27][28].
The paper ends with conclusions and an outlook (section VI).

BOUNDARY METHOD
We start by considering the diffusion equation where c is the concentration field, D the diffusion coefficient. In the immersed boundary boundary methods [13,14], the Eulerian source term s is used to mimic the effects of the boundaries of bubbles or droplets on the concentration.
The boundaries of bubbles or droplets are discretized into a series of Lagrangian points.
The Eulerian and Lagrangian sources are related to each other through a regularized delta function where x and X l are the position vectors of the Eulerian and Lagrangian points; S the Lagrangian source term; δ the delta function, respectively.
To enforce the prescribed concentration fields on the boundary, we define the Lagrangian concentration field. Using the regularized delta function again, this relation can be expressed as follows where C Γ is the Lagrangian concentration field which is prescribed, known beforehand, on the boundary. In the discretized form, the diffusion equation for the kth step is solved through the following procedures. First, an intermediate "guessed" concentration fieldc is calculated from the Eulerian source term of the last step s k−1 , with Here, the diffusion term ∇ 2 c is discretized by a second-order explicit scheme.
Next, we interpolate the intermediate concentration field from Eulerian (c) to Lagrangian (C) grid points through the discrete delta function δ h , i. e.
ApparentlyC does not satisfy the boundary condition C Γ . In order to achieve C Γ , from Eqn. 5 the Lagrangian source term S k for the current time step is derived as  Epstein-Plesset [29] results. However, our results deviate from that of a quasistatic approximation The next step is to spread the Lagrangian source term S k to the Eulerian counterpart s k through the discrete delta function δ h again, expressed as Finally, the concentration field with the Eulerian source term s k at kth step is solved Here the Crank-Nicolson scheme is adopted to ensure the stability of the code.
This ends one time step, after which the next time step is treated in the same way.
The regularized delta function used in the present study is defined as where φ in the present implementation is based on the four-point version of Peskin [13].

III. VALIDATION OF THE SCHEME FOR A SINGLE BULK BUBBLE AND A SINGLE SURFACE BUBBLE
We will now validate the scheme introduced in the previous section. We will assume two test cases: a spherical bubble in the bulk, whose growth or shrinkage behavior is analytically known since Epstein and Plesset [29] (section III.1), and a surface nanobubble, which in the pinned case has a stable equilibrium contact angle given by Eq. (1), and in the unpinned case either shrinks and then fully dissolves or grows and then finally detaches (section III.2).
All the simulations that are shown below are performed with nitrogen bubble, for which the material parameters are D = 2 × 10 −9 m 2 /s, ρ g = 1.165 kg/m 3 , and c s = 0.017 kg/m 3 .

A. The Epstein-Plesset bubble
In still liquid in an infinite domain, the mass loss or gain of a spherical bubble of radius R is given by the concentration gradient ∂c ∂r R at its surface and the diffusion constant D, equation together with Eq. (12) and the boundary condition far away from the bubble, c(r → ∞, t) = c ∞ , to obtain an ordinary differential equation (ODE) for the bubble radius Here the prescribed C Γ is calculated from Henry's law, taking the effects of surface tension into account, i.e., C Γ (R, t) = c s (1 + 2σ/R), where c s is the saturation concentration. Note that for small bubbles the Laplace pressure leads to an enhanced density, obtained from the ideal gas law, and this effect of the surface tension must also be taken into account. Equation (13) can be solved analytically to obtain R(t) [29]. Obviously, also in the simulations the bubble is assumed to keep its spherical shape during the diffusion process and equation (12) is used to update the bubble radius and the Lagrangian coordinate X during the simulation.
Our numerical results of the relation between the bubble radius and time based on the scheme developed in the previous section are shown in figure 3 and compared with the analytical results (or the results from Eq. (13)). Three cases are considered. In figure   3(a), the bubble surface concentration and gas density are kept constant, in figure 3(b), the density of the gas is kept constant and we use the Henry's law to calculate C Γ , and in figure   3(c), we vary the density of the bubble according to the ideal gas law and again the Henry's law is used to calculate C Γ . For all the cases, our simulations show excellent agreements with the predictions from Eq. (13).
We now come to dissolving or growing surface bubbles and droplets ("sessile droplets") [1,30]. For this axisymmetric case, Popov [31] could exactly solve the quasi-static case ∂ t c ≈ 0, i.e., the diffusion equation reduces to a Laplace question. For evaporating droplets as in the case of Popov, this in general is a very good approximation. Later the Popov model was also applied to surface nanobubbles [10]. Then the gas concentration at the interface is again given by Henry's law which for surface bubbles takes the form C Γ (R, t) = c s [1 + 4σ sin θ/(P 0 L)].
To check how important the time dependence of the concentration field is, we apply Popov's model for a dissolving bubble with a fixed contact angle of 90 • , written as One can see that the only difference between Eq. (13) and Eq. (14) is that in Eq. (14) the time dependent term in the right hand side of Eq. (13) is eliminated. It is observed from figure 3 that when Henry's law is used while the bubble density is kept constant, the quasi-static assumption of Popov's model leads to an overestimation of the bubble lifetime.
Therefore in the following, for appropriately simulating the diffusive dynamics of the bubbles, we do not use the quasi-static approximation, but employ the full diffusion equation with Henry's law for the bubble surface concentration and the ideal gas law for the bubble density. For nano-bubbles with pinned contact line, an ODE for the diffusive contact angle dynamics was derived in Ref. [10], namely with f (θ) = sin θ 1 + cos θ A stable nanobubble can therefore be formed with the condition of Eq. (1) where the bubble contact angle θ is a constant and stable.  We take the opportunity here to discuss the assumptions that lead to Eq. (15). Henry's law is used when deriving Eq. (15), however the gas density is assumed constant and the process is assumed quasi-steady. Let's first focus on the quasi-steady assumption. The typical diffusion time scale is t d = R 2 /D, while the evaporation/dissolution time scale . For a water droplet evaporation, t e /t d is of the order of 10 5 , thus Eq. (15) is a rather good approximation [32] without considering the time dependent term of the diffusion equation. However, for a gas bubble, t e /t d is of the order of 10 2 , thus the quasi-steady condition can not be valid anymore, as also shown in figure 3(b). Also the gas density might vary because of the Laplace pressure. However, it is easy to see from Eq. (15) and Eq. (3) that these considerations are only relevant for the time scale of the evolution towards the equilibrium contact angle θ e , not for the value of θ e itself.

PINNED CASE
We now move to the case of two bubbles, for which the general argument for nanobubble stability is not available anymore. One exception is the case where two bubbles are far away from each other. For this case Dollet and Lohse [18] theoretically show that pinning also suppresses the Ostwald ripening process between neighbouring surface nanobubbles. But this case is not given in most experiments, in which the nanobubbles sit very close to each other and nonetheless can remain stable for very long time [7]. In this paper we will now show with numerical simulations that this stabilization of a pair of surface bubbles through pinning is indeed not limited to bubbles far away from each other but also holds for bubbles that are close.  Here we can use the two-dimensional axis-symmetric simulations because the contact angle is 90 • . the coarsening process as seen in shaken compartimentalized granular matter [21].
We note that though here we give the results only for surface nanobubbles, corresponding results should also hold for surface nanodroplets. We also note that for the parameters of this study here, the dominant coarsening process is Ostwald ripening, i.e., mass exchange by diffusion, but for other parameters (e.g. larger oversaturation) the dominant process can also be bubble coalescence. To map out the parameter space when Ostwald ripening will be dominant and when bubble coalescence will be the subject of future work. Correspondingly, in future work we also want to extend this study from surface bubbles or surface droplets in a row to those in a two-dimensional array as experimentally done in e.g. Refs. [17,33] or to randomly distributed surface bubbles or droplets as in Ref. [22]. Future work can also address how heterogeneities on the gas-water interfaces through e.g. local surfactant accumulation can affect the overall dynamics of the bubble ensemble.
Finally, we caution the reader: Our results are based on continuum theory and hydrodynamic equations. However, at very short length scales the continuum approximation will break down. In very recent molecular dynamics (MD) simulations, Maheshwari et al.
[ 34] have revealed that in certain cases (very strong attraction between the dissolved gas molecules and the surface) surface nanobubbles very close to each other can communicate through a "new channel", namely diffusion of gas from one surface bubble to the other along the surface, and not through the bulk. If this is the case, some sort of ripening process of neighboring surface bubbles may be possible in spite of hydrodynamic stability against Ostwald ripening.

Conflicts of interest
There are no conflicts to declare.