Influence of Surface Heterogeneity on Morphology of Interfacial Nanobubble

Gaseous domains formed on solid–liquid interface have attracted scientists’ attentions in recent 2 decades, and the existence of interfacial nanobubble (INB) has been basically confirmed. However, an overall understanding on INB is still lacking. This research studied the influence of surface chemical heterogeneity on the morphology of INB by molecular dynamics simulations technique. The results showed that the gaseous domains could not nucleate on the hydrophilic substrate, while only dense gas layer (DGL) could be observed from the time-averaged density map for homogeneously hydrophobic substrate due to the random moving of INB. If there was a hydrophobic patch on the hydrophilic surface, INB could form on the hydrophobic patch with contact line being pinned at the boundary of the patch. In this case, the contact angle (gas-side) increased with the gas oversaturation degree and decreased with surface hydrophobicity of the patch. For the case that there existed a more hydrophobic patch/site on the hydrophobic surface, the INB could have moved randomly along the hydrophobic surface, but its receding contact line was pinned by the more hydrophobic patch/site. Hence, the INB could only move in the vicinity of this pinning patch/site, so that an INB profile instead of a DGL formed due to the pinning effect, and the apparent contact angle of the INB is significantly lower than the actual one. Throughout this study, the apparent INB we observed from experiments may be different from its instantaneous state and is significantly affected by surface heterogeneity.


INTRODUCTION
Interfacial nanobubble (INB) was first proposed by Parker et al. (1994) in 1994 due to the steps and the long-range attractive force in the force curve between two approaching hydrophobic surfaces. Until 2000, the first image of solid-water interfacial nanobubble was independently obtained through tapping mode of the atomic force microscopy by Lou et al. (2000) and Ishida et al. (2000). Since then, research on INBs have been extensively focused on due to their potential applications in many fields including froth flotation, protein adsorption, drag reduction, catalysis, and electrolysis (Lohse and Zhang, 2015).
INBs are usually referred to as gas aggregates that are confined in a spherical cap with height less than 100 nm and contact line diameter less than 1 μm (Alheshibri et al., 2016). According to Young-Laplace equation, the internal pressure of INB could be described as where P is the internal pressure of INB, γ is the gas-liquid interfacial tension, R is the curvature radius of the INB, P 0 is the ambient pressure, L is the width of INB, and θ is the contact angle of the INB. For an INB with L = 100 nm and θ 30°(from gas side) in pure water, the calculated internal pressure of INB is p = 1.47 MPa = 14.7 atm, where γ 73mN/m(20℃). The huge internal pressure of INB would lead to the diffusive outflux of gas molecules, and the INB would dissolve at a time scale of τ R 2 /D 50μs according to the diffusion equation proposed by Epstein and Plesset (1950). Therefore, the research on INB had always been companied with the query about its existence or not (Ljunggren and Eriksson, 1997;Mao et al., 2004;McKee and Ducker, 2005;Takata et al., 2006). However, extensive studies in the last 2 decades demonstrated that INBs with long life-time indeed exist through various methods including total internal reflection fluorescence microscopy (Chan and Ohl, 2012;Tan et al., 2017), neutron reflectometry (Steitz et al., 2003), quartz crystal microbalance (Seo et al., 2007), attenuated total reflectance infrared spectroscopy , and surface plasmon resonance (Martinez and Stroeve, 2007). In addition, INBs could exist on various substrates, including hydrophobic surfaces, such as HOPG (Fang et al., 2016;Lu et al., 2014), polystyrene (Li et al., 2016), molybdenite (Wang et al., 2019), OTS-coated Si (Zhang et al., 2006), and decanethiol-coated gold (Zhang, 2008), as well as hydrophilic surfaces, such as mica (Lou et al., 2000;Wang et al., 2019), gold (Holmberg et al., 2003), and metals (Cavicchi and Avedisian, 2007). Throughout these studies, INBs are found to have abnormal contact angle, which is much lower than the contact angle derived from Young's equation. INBs usually form on hydrophobic surface, but it is beyond understanding why they could also form on the hydrophilic substrates, as well as why they have various contact angles even on the same surface (Lohse and Zhang, 2015). To data, the nature of INB has not been fully understood yet. Molecular dynamics (MD) simulation, as a state-of-the-art technology to simulate the nanoscale phenomenon, is an ideal method to research the nucleation and dynamics of INBs from the molecular-level insights (Pu et al., 2020). By using MD simulation, many discoveries that are hardly acquired from experimental techniques were made. Dammer and Lohse (2006) discovered the gas enrichment phenomenon at liquid-solid interface. Xiao et al. (2017) found that the MD method could simulate the process of solvent exchange that have been widely used to produce INBs in experiments. Weijs et al. (2012) researched the formation of INBs on the homogeneous surface and found that the INB could form under a low gas solubility. Maheshwari et al. (2016), Maheshwari et al. (2018), Maheshwari et al. (2020) studied the stability of INBs and found that the INBs could be pinned by the more hydrophobic part, and the INB nucleation process is affected by the surface morphology. Recently, Yen et al. (2021) demonstrated that the gas molecules might migrate from gas enrichment layer to the inside of INB through the study in potential of mean force. Wang et al. (2022) researched the interfacial phenomenon in cement fluidity and found the interfacial structure of water is influenced by the shearing fluid, which is significant in many interfacial systems such as boundary slip. However, seldom researches have been conducted to study the formation and morphology of INBs on different types of substrates.
In this study, we employed MD simulations to study how INBs are influenced by the surface heterogeneity through the formation of INB on various types of substrates including homogeneous surface with different hydrophobicity, hydrophilic surface  In this case, 300 gas particles and 4,000 liquid particles were randomly filled into the simulation box with a typical size of 10.20 × 2.75 × 11.00 nm 3 in x-, y-, and z-direction. The blue and green beads represent the liquid and gas particles, respectively, and the silver, orange, and red beads represent the fixed solid particles of S L , S M and S H . patterned with hydrophobic patch, and medium hydrophobic surface patterned with more hydrophobic patch/site. The current work facilitates understanding the interfacial phenomenon of surface nanobubble or nanodroplet and how it affects the morphology we observed by analytical instruments.

SIMULATION METHODOLOGY
MD simulations were performed with the open source code GROMACS 2019.6 package (Abraham et al., 2015) and were visualized by VMD software (Humphrey et al., 1996). Like many previous studies (Lu, 2019;Maheshwari et al., 2016;Maheshwari et al., 2020), this work employed Lennard-Jones beads to represent the liquid (L), solid (S), and gas (G) particles. Solid particles were fixed in a fcc lattice, while liquid and gas particles could move freely. The expression of LJ potential used to describe the interaction between different particles is shown as follows: where ε ij , σ ij , and r ij is the interaction strength, characteristic size, and separation distance between particles i and j. The LJ parameters of various particles, shown in Table 1, were determined by referring to previous reports (Lu, 2019;Maheshwari et al., 2020), where S L , S M , and S H represent low, medium, and high hydrophobic surfaces, respectively. A cutoff radius of 1.2 nm was used to compute the LJ potential. The leap-frog algorithm was used to integrate the equation of motion, with a time step of 2 fs. Periodic boundary conditions were  In this case, 300 gas particles and 4,000 liquid particles were randomly filled into the simulation box with a typical size of 10.20 × 2.75 × 11.00 nm 3 in x-, y-, and z-direction. The blue and green beads represent the liquid and gas particles, respectively, and the silver and orange beads represent the fixed solid particles of S L and S M . A movie for the process is available as Supplementary Video S1. used in all three directions. The simulations were conducted in NVT ensemble for 1 ns where the system volume was unchanged and the temperature was kept at 300 K, and then were conducted in the NP z T ensemble where the pressure maintained at 1 bar. Semi-isotropic pressure coupling was used, and the simulation box could scale in Z direction to maintain constant pressure. The substrate was divided into two parts, as shown in red/orange and silver colors in Figure 1. The upper part of the solid was the surface we studied, while the bottom part was always hydrophilic particles (S L ). Supplementary Figure S1 shows the typical initial configurations of Figures 2, 8. In the initial configuration, the liquid and gas particles were randomly dispersed above the solid substrate, and the gas particles were mainly arranged near the surface to aid the bubble nucleation for decreasing the equilibrium time. Supplementary Figure S2 shows the evolution of potential, temperature, and pressure of the system in Figure 2 with time, which clearly indicates that the simulation system already reaches equilibrium before 10 ns. In this study, coordinates were output every 1,000 steps, and the typical simulation time was 120 ns, unless otherwise specified. After the simulation, both the average densities of gas particles and of liquid particles in each bin of 0.1 × 2.75 × 0.1 nm 3 within a certain simulation time range were calculated. The obtained density distribution of gas particles was presented as a color map such as Figure 3. The density distribution of liquid particles was analyzed to calculate the contact angle. The Z coordinate of mass center of the top layer of solid plus σ SL /2 was considered as the solid-liquid interface. The gas-liquid interface was defined as the position where the liquid density equals half of the bulk liquid density. The height, radius, and contact angle of the INBs were obtained by fitting the gas-liquid interface to a circle function.

RESULTS AND DISCUSSION
Gaseous Domains on Homogeneous Surfaces Figure 1 shows the simulation snapshots of S L , S M , and S H . We could see an ordered layer of liquid particles formed adjacent to the S L surface. Bulk nanobubble is nucleated due to the gas oversaturation, but they could not adhere to the substrate because the hydrated layer on the hydrophilic surface is difficult to rupture. For S H surface, a dense gas layer (DGL) rather than an INB forms at the solid-liquid interface, which could also be seen from the averaged density map of gas molecules shown in Figure 3B. Only on a medium hydrophobic surface (S M substrate) could INB be formed.
When we inspected the trajectory of simulation for S M surface, we found that the INB moved along the solid-liquid interface randomly all the time. Figures 2A-C show the snapshots of INB formed on S M surface at three different simulation moments. The movie for the process is also available as Supplementary Video S1. The INB located at different positions at different moments. A similar phenomenon was also reported by Wu et al. (2017) who proposed that this random movement of INB was caused by the thermal fluctuation, and the movement frequency increased with the surface hydrophobicity. We believe Brownian movement might be also responsible for this phenomenon. When we averaged the density distribution of gas molecules over a time range of 100 ns, as shown in Figure 3A, we can see only a DGL profile instead of an INB on S M surface due to the rapid movement of INB, which is a visual effect from time-average Frontiers in Materials | www.frontiersin.org March 2022 | Volume 8 | Article 824125 5 rather than the instantaneous situation shown in Figure 2. Considering that the movement of INB is much faster than the temporal resolution of experimental instruments as that the former is at nanosecond magnitude while the latter is usually larger than millisecond for imaging INBs, this movement of INB is expected to influence the morphology and the calculated contact angle of the INB in experiments.

Hydrophilic Substrate Patterned by Hydrophobic Patch
The ideal homogeneous surface without pinning site rarely exists, and almost all the substrates coexist with heterogeneous pinning sites more or less. For hydrophilic substrate, it is usually hard to totally avoid the existence of hydrophobic patches on the surface. Figures 4, 5 show the gas density distribution maps for the hydrophilic substrates patterned by medium hydrophobic and strong hydrophobic patches, respectively, under different degrees of gas oversaturation. These results indicate that INBs could form on the hydrophobic patches patterned on the hydrophilic substrates. As shown in Figures 4, 5, the height and contact angle of INB increase with increasing gas particle number, while the width is unchanged due to contact line pinning, indicating that the height and contact angle of INBs are positively correlated with the degrees of gas oversaturation for INBs formed on the hydrophobic patches patterned on the hydrophilic substrate. However, if the degree of gas oversaturation is not high enough, such as 150 gas particles in Figure 4A or 200 gas particles in Figure 5A, a DGL instead of an INB would form, which is due to the rapid movement of small gas domain along hydrophobic patch just like what Figure 2 shows. In addition, the minimal gas particle numbers required to form an INB for S H patch is higher than that for S M patch, indicating that the minimal degree of gas oversaturation required to form INBs increases with the increase of hydrophobicity of the surface below INB.  Figure 6A, the intensities of the peaks in these curves increase with increasing gas particle number, while the peak intensities are almost unchanged with the increase of gas number for high hydrophobic pattern as shown in Figure 6B. It indicates that the gas density distribution next to the surface is affected by the degree of gas oversaturation for medium hydrophobic surface, but it is not affected by the gas oversaturation for strong hydrophobic surface. Figure 6C shows the contact angles of INBs formed on S M and S H patches when using different gas number. We can see that the contact angle increases with the increase of gas number, and the contact angles for S H are lower than those for S M , indicating that the contact angle of INBs increases with the gas oversaturation degree for a given surface hydrophobicity and decreases with the surface hydrophobicity for a given gas oversaturation. Figure 6D shows the number density of gas particles in the bulk increases with the increase of INB curvature, which might be due to the higher Laplace pressure for INB with a larger curvature.

Hydrophobic Substrate Patterned by Stronger Hydrophobic Patch
For a hydrophobic substrate, there might exist some patches with stronger hydrophobicity. Figure 7 shows the INBs There are 9,000 liquid molecules in these simulations, and the initial box size is 21.19 × 2.75 × 13.00 nm 3 . The yellow beads represent S M particles, the red beads represent S H particles, and the silver beads represent S L particles. formed on S M substrate patterned by S H patch. We could see that a cap-shaped INB forms on S H patterned S M surface, but the contact line stands on the S M surface rather than being pinned at boundary between S H and S M , which is not like the case of Figure 5. Although the contact line locates on the same surface (S M ), contact angle is found increasing with the increase of gas oversaturation rather than keeping constant. The contact angles of Figures 7B-D (Zhang et al., 2012). The gas domain of Figure 7A is more like a DGL, which is because that the gas oversaturation degree is too low to form an INB. The gas enrichment layer is also found existing outside of the INB, which is also reported by previous researches (Yen, 2020) and is considered to account for the stability of INB (Tortora et al., 2020). When we examined the simulation trajectory, we found that the advancing contact line of INB was not pinned by S H patch, but the receding contact line was always pinned. Hence, the INB could only move around the S H pattern, producing an apparent INB with contact line on S M surface and center on S H pattern. This phenomenon is more obvious when the S H pattern is diminished to a site, just as shown in Figure 8 where an INB forms on the S M surface patterned by S H site. It is apparent that the receding contact line of INB was always pinned by S H site, but the advancing contact line could not be pinned, so that the INB moved in a round area centered on S H site. These moving trails formed an apparent INB shown in Figure 9A. In this case, the apparent INB we observed is just a shadow rather than the actual nanobubble, and the time-averaged FIGURE 10 | Time-averaged gas density distribution map for substrate with two patches under gas density number of (A) 600 and (B) 1,000. Both the liquid particle number is 15,000. A movie for the process of (B) is available as Supplementary Video S3.
Frontiers in Materials | www.frontiersin.org March 2022 | Volume 8 | Article 824125 8 contact angle is also not the actual contact angle. Figure 9B shows a sketch for the influence of pinning site on contact angle of the INB. In this case, the white spherical cap represents the actual INB, while the red solid line represents the apparent INB. Obviously, the apparent contact angle we could see (θ 2 ) is much lower than the true value (θ 1 ) due to the pinning effect.
For the case with several patches, there is similar phenomenon. Figures 10A,B show the INB formed on the S M substrate with two S H patches where gas particle number are 600 and 1,000, respectively. For the case of 600 gas particles, the INB was not large enough to occupy two patches. Hence, two INBs were formed on the two patches, respectively ( Figure 10A). When the gas particle number increased to 1,000, a large INB crossing two patches was formed ( Figure 10B). This INB kept moving randomly around the two patches, with receding contact line always being pinned by the patch (Supplementary Video S3). In this case, the role of the two small patches is similar to a large patch just like what Figure 7 shows, and the influence of the two patches on morphology of the INB is also like what Figure 9B illustrates.

CONCLUSION
In this study, the influence of surface heterogeneity on interfacial nanobubble (INB) was studied by using MD method. The following conclusions were obtained: (1) For smoothly homogeneous ideal-surface, the gaseous domain could not nucleate on the hydrophilic surface, and only dense gas layer (DGL) forms on strong hydrophobic surface. For medium hydrophobic surface, although INB could nucleate on the substrate, it moves freely along the solid-liquid interface, so that the time-averaged density distribution map of gas particles exhibits a DGL shape.
(2) For hydrophilic substrate patterned by hydrophobic patch, the INB could form on the patch with the contact line being pinned at the boundary. The contact angle of INB increases with gas oversaturation degree, and decreases with the surface hydrophobicity of the patch. In addition, the minimal gas oversaturation degree required to form an INB increases with the surface hydrophobicity of the patch. If the degree of gas oversaturation is lower than this required minimal value, a DGL instead of an INB would form.
(3) For hydrophobic substrate patterned by stronger hydrophobic patch/site, the INB could also form. We found that the receding contact line rather than advancing contact line was always pinned by the boundary of the patch, so that the INB randomly moved around the patch. Hence, a cap-shaped INB centered on the patch was seen from the averaged gas density map. For this case, the width of the apparent INB was much larger than that of instantaneous state, and thus the contact angle of the apparent INB was lower than the actual one. If the INB is large enough to occupy several patches, these patches behave like a whole large patch.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HY conceived and preformed experiments, analyzed the data, and wrote the article. FZ analyzed the data. YX conceived the experiments. XG and YC wrote and edited the article.