2. Molecular simulations of nanobubble motion in a temperature gradient

Figure 1. (a) A snapshot of a gas nanobubble in water between two solid plates in MD simulations. The top plate has a lower temperature $T_l$ , and the bottom plate has a higher temperature $T_h$ . (b) Sketch of a gas nanobubble in water between two solid plates. Here, $\gamma$ is the local surface tension of the bubble surface, and the bubble portion with a lower (higher) temperature has a higher (lower) surface tension, denoted as $\gamma +$ ( $\gamma -$ ). This creates a Marangoni flow around the bubble, resulting in a Marangoni force $F_{Ma}$ on the bubble.

Molecular dynamics (MD) simulations are used to simulate the motion of bulk gas nanobubbles in liquid between two solid plates whose temperatures are different. We adopt the open-source code LAMMPS (Plimpton Reference Plimpton1995). As shown in figure 1(a), the molecular system consists of water molecules (represented in orange), gas atoms (represented in green), atoms of the bottom solid (represented in red) with a higher temperature $T_h$ , and atoms of the top solid (represented in blue) with a lower temperature $T_l$ .

The mW water potential is adopted to model water. The mW water model is a monatomic water model proposed by Molinero & Moore (Reference Molinero and Moore2009). It uses the Stillinger–Weber potential:

(2.1) \begin{align} E&=\sum \limits _{i}{\sum \limits _{j\gt i}{{{\phi }_{2}}\left ( {{r}_{ij}} \right )}}+\sum \limits _{i}{\sum \limits _{j\ne i}{\sum \limits _{k\gt j}{{{\phi }_{3}}\left ( {{r}_{ij}},{{r}_{ik}},{{\theta }_{ijk}} \right )}}}, \nonumber \\ {{\phi }_{2}}\left ( r \right )&=A\varepsilon \left [ B{{\left ( \frac {\sigma }{r} \right )}^{p}}-{{\left ( \frac {\sigma }{r} \right )}^{q}} \right ]\exp \left ( \frac {\sigma }{r-a\sigma } \right ), \nonumber \\ {{\phi }_{3}}\left ( r,s,\theta \right )&=\lambda \varepsilon {{\left ( \cos \theta -\cos {{\theta }_{0}} \right )}^{2}}\exp \left ( \frac {\gamma \sigma }{r-a\sigma } \right )\exp \left ( \frac {\gamma \sigma }{s-a\sigma } \right ). \end{align}

Here, $A=7.049556277$ , $B=0.6022245584$ , $\varepsilon =6.189$ kcal mol $^{-1}$ , $\sigma =2.3925 {{\unicode{x00C5}}}$ , $a=1.80$ , $\lambda =23.15$ , $\gamma =1.20$ , $\cos \theta _0= -1/3$ , $p=4$ and $q=0$ . Compared to other water models such as TIP4P/2005 and SPC/E, the mW water model has a good accuracy of surface tension and a low computational cost (Molinero & Moore Reference Molinero and Moore2009). For example, a large time step can be chosen for the mW water model (Molinero & Moore Reference Molinero and Moore2009). In the current simulations, it is 8 fs.

Except for water itself, the intermolecular potentials $U$ between $i$ -type atoms and $j$ -type atoms are simulated with the standard Lennard-Jones 12–6 potential:

(2.2) \begin{equation} U({{r}_{ij}}) = \begin{cases} 4\varepsilon _{ij} \left [ {{\left ( \dfrac {\sigma _{ij} }{{{r}_{ij}}} \right )}^{12}}-{{\left ( \dfrac {\sigma _{ij} }{{{r}_{ij}}} \right )}^{6}} \right ] & \text {if} \,\,\,{{r}_{ij}}\leqslant {{r}_{c,ij}},\\ 0 & \text {if}\,\,\, {{r}_{ij}}\gt {{r}_{c,ij}}.\\ \end{cases} \end{equation}

Here, $r_{ij}$ , $\varepsilon _{ij}$ , $\sigma _{ij}$ and ${r}_{c,ij}$ are the pairwise distance, energy parameter, length parameter and cut-off distance, respectively. The complete list of parameters among water (W), gas (G) and solid (S) is given in table 1. The gas modelled by the standard 12–6 Lennard-Jones potential represents nitrogen, and it has density $\rho _{\infty }=1.147$ kg m⁻³ at 1 atm and 300 K. The molar mass of the gas is 28 g mol⁻¹. As noted earlier, clean bulk nanobubbles in saturated or undersaturated environments tend to dissolve in microseconds, as predicted by the Epstein–Plesset theory (Epstein & Plesset Reference Epstein and Plesset1950; Lohse & Zhang Reference Lohse and Zhang2015). We indeed find that in our molecular simulations, a nanobubble dissolves over time, leading to a time-dependent bubble radius. This obviously complicates the theoretical modelling and goes beyond the scope of current work to validate the YGB59 theory (where a constant bubble radius is assumed) to the nanoscale. To circumvent this problem, one can tune the energy parameter, distance parameter and cut-off distance between gas and water to obtain a small gas solubility. In current simulations, the cut-off distance between water and gas is set to be $2^{1/6}\sigma _{GW}=3.4$ nm so that only the repulsive force between gas and water is effective, leading to an extremely small gas solubility. The same simulation strategy was used when simulating nanobubble cavitation (Shekhar et al. Reference Shekhar, Nomura, Kalia, Nakano and Vashishta2013).

Table 1. Interaction parameters among water (W), gas (G) and solid (S).

The box has lateral size $L_x=19.2$ nm and $L_y=19.2$ nm; see figure 2(b). The height $L_z$ of this box is fixed during the bubble movement, but it varies for different cases with different bubble radii, as discussed below and shown in table 2. The height of each solid plate at the lower and upper ends of the domain is $L_p=0.96$ nm. Both the top and bottom plates are made of five layers of atoms with the same density as water for simplicity. Periodic boundary conditions are applied in the lateral directions. On the temperature conditions, the bottom hot plate is held at 350 K or 325 K, and the top cold plate is held at 300 K using the Langevin thermostat with damping factor 100 fs. A constant spring force 10 kcal mol⁻¹ ${{\unicode{x00C5}}}$ ⁻² is applied to each atom in the plates to tether them to their initial position. The positions and velocities of the fluid atoms (liquid and gas) are updated based on Newton’s second law without extra thermostating.

Table 2. Simulated three cases with different settings.

Figure 2. (a) Snapshots of the bubble movement in the MD simulations at three different times for case 1. (b) For the same simulation, the three-dimensional trace of the bubble motion is shown.

The tank is initially filled with 467 200 water atoms and is isothermal with temperature 300 K. Then $N$ water atoms next to each other centred at $h_0=28.8$ nm are switched into gas atoms, which forms a nanobubble. Notably, the nanobubble is initially tethered to its initial position through a spring force so that it can only expand or shrink without Brownian motions. During this process, the top plate can move vertically to maintain the ambient liquid pressure at 1 atm in response to the appearance of the nanobubble. Two different situations with $N=1196$ and $N=3316$ are simulated. The corresponding radii of the obtained nanobubbles at equilibrium are measured to be 3.1 nm and 5.0 nm, consistent with the prediction of the real gas law (Zhang & Lohse Reference Zhang and Lohse2023). At equilibrium, the position of the top plate barely moves, and $L_z=38.4$ nm for $R=3.1$ nm, and $L_z=39$ nm for $R=5.0$ nm, are measured and listed in table 2. After obtaining the initial configuration, the Langevin thermostat is applied to the two plates. The formed nanobubbles are still tethered to their initial positions for another 100 ns, during which a steady temperature gradient is produced between the two solid plates. After 100 ns, the nanobubbles are released so that they can move freely in the liquid. Due to the thermal motions of nanobubbles, for each case, several realizations (see table 2) have to be run to improve the statistics.

4. Extension of the YGB59 theory towards the nanoscale

Though Young et al. (Reference Young, Goldstein and Block1959) (YGB59 theory) derived the correct expression of bubble velocity, their work did not include the expression of Marangoni forces and drag forces. This has led to some approximate expressions in the literature (Morick & Woermann Reference Morick and Woermann1993; Lubetkin Reference Lubetkin2003; Zeng et al. Reference Zeng, Chong, Wang, Diddens, Li, Detert, Zandvliet and Lohse2021), which are valid only at the macroscale. Here, the YGB59 equation is revisited, and more general expressions for Marangoni forces and drag forces are obtained, which show important nanoscale effects, such as those of the thermal resistance and the slip at the interface.

The governing dimensionless numbers for both fluids (liquid and gas) are the Reynolds number $Re=\rho u_0 R/\mu$ and the Marangoni number $Ma=({{\rm d}\gamma }/{{\rm d}T})({{\rm d}T}/{{\rm d}z})({{R}^{2}}/ ( \mu \alpha ))$ . Here, $\rho$ , $\mu$ , $u_0$ , $R$ , $\gamma$ and $\alpha$ are the density, dynamic viscosity, velocity, radius of the nanobubble, surface tension and thermal diffusivity, respectively. Assuming that $Re$ and $Ma$ are much smaller than unity for both fluids, which is valid at the small scale (small $R$ ) of nanobubbles, the momentum equation and continuity equation for both fluids are

(4.1) \begin{equation} \mu\, {\boldsymbol {\nabla }^{2}} {\boldsymbol u}=\nabla p, \end{equation}

(4.2) \begin{equation} \boldsymbol {\nabla } \cdot {\boldsymbol u}=0, \end{equation}

respectively. In the momentum equation, the unsteady term $\rho\, \partial u / \partial t$ is dropped. This is justified because the ratio of the time scale required for the water velocity to adjust to the changes in the bubble position ( $\rho R^2/ \mu$ ) to the time scale of bubble movement ( $R/u_0$ ), i.e. $Re$ , is small. In fact, the largest Reynolds number in our simulations is evaluated to be 0.01 for case 1. Therefore, the water velocity field may be assumed to be quasi-static. The heat transfer equation for both fluids is

(4.3) \begin{equation} {{\nabla }^{2}}T=0. \end{equation}

The general solution to the temperature field for both fluids in spherical coordinates is

(4.4) \begin{equation} T={{A}_{1}}+\frac {{{A}_{2}}}{r}+{{A}_{3}}r\cos \theta +\frac {{{A}_{4}}}{{{r}^{2}}}\cos \theta . \end{equation}

Here, $A_i$ ( $i=1,2,3,4$ ) is the coefficient to be determined by the boundary conditions. At the far field, the temperature distribution $T^o$ for the outer fluid (i.e. water here – the superscript $^o$ is for the liquid, and the superscript $^i$ is for the inner fluid, i.e. gas hereafter) is

(4.5) \begin{equation} T^o=T_0+\Gamma r\cos {\theta }, \end{equation}

where $T_0$ is the temperature at the centre of the bubble, and $\Gamma$ is the temperature gradient. At the fluid–fluid interface, the continuity of the heat flux requires

(4.6) \begin{equation} k^o\frac {\partial T^o}{\partial r}= k^i\frac {\partial T^i}{\partial r}, \end{equation}

where we have assumed a spherical bubble with the coordinate system at its centre. For the temperature condition across the interface, motivated by the observation from our simulations, we adopt the temperature jump condition across the interface instead of the continuity of temperature, so that one can have

(4.7) \begin{equation} T^o-T^i=k^o\frac {\partial T^o}{\partial r}G, \end{equation}

where $G$ is the thermal resistance.

Using the above temperature boundary conditions, the temperature distribution can be obtained for both fluids:

(4.8) \begin{align} & {{T}^{o}}={{T}_{0}}+{\Gamma }r\cos \theta +\frac { ( {{k}^{o}}-{{k}^{i}} )R+{{k}^{o}}{{k}^{i}}G}{( 2{{k}^{o}}+{{k}^{i}} )R+2{{k}^{o}}{{k}^{i}}G}\,\frac {{\Gamma }{{R}^{3}}\cos \theta }{{{r}^{2}}}, \quad r\geqslant R,\nonumber \\ & {{T}^{i}}={{T}_{0}}+\frac {3{{k}^{o}}R}{ ( 2{{k}^{o}}+{{k}^{i}} )R+2{{k}^{o}}{{k}^{i}}G}{\Gamma }r\cos \theta, \quad r \lt R. \end{align}

In terms of the velocity field, the usual way to solve (4.1) and (4.2) is by introducing the Stokes stream function $\psi (r,\theta )$ (Bird, Stewai & Lightfoot Reference Bird, Stewai and Lightfoot2006)

(4.9) \begin{align} & {{u}_{r}}=-\frac {1}{{{r}^{2}}\sin \theta }\frac {\partial \psi }{\partial \theta },\nonumber \\ & {{u}_{\theta }}=\frac {1}{r\sin \theta }\frac {\partial \psi }{\partial r}. \end{align}

Equation (4.1) is then rewritten in terms of the Stokes stream function $\psi$ as

(4.10) \begin{equation} {{E}^{4}}\psi =0, \end{equation}

where ${{E}^{2}}=({\partial }/{\partial {{r}^{2}}})+({\sin \theta }/{{{r}^{2}}})({\partial }/{\partial \theta }) ( ({1}/{\sin \theta })({\partial }/{\partial \theta }) )$ . The general solution to (4.10) for both fluids is

(4.11) \begin{align} & {{\psi }^{o}}=\left ( {{C}_{1}}{{r}^{-1}}+{{C}_{2}}r+{{C}_{3}}{{r}^{2}}+{{C}_{4}}{{r}^{4}} \right ){{\sin }^{2}}\theta, \quad r\geqslant R, \end{align}

(4.12) \begin{align} & {{\psi }^{i}}=\left ( {{D}_{1}}{{r}^{-1}}+{{D}_{2}}r+{{D}_{3}}{{r}^{2}}+{{D}_{4}}{{r}^{4}} \right ){{\sin }^{2}}\theta,\quad r\lt R. \end{align}

The resulting expressions for the velocities are

(4.13) \begin{align} & u_{r}^{o}=-2\cos \theta \left ( \frac {{{C}_{1}}}{{{r}^{3}}}+\frac {{{C}_{2}}}{r}+{{C}_{3}}+{{C}_{4}}{{r}^{2}} \right ), \end{align}

(4.14) \begin{align} & u_{\theta }^{o}=\sin \theta \left ( -\frac {{C}_{1}}{{r}^{3}}+\frac {{{C}_{2}}}{r}+2{{C}_{3}}+4{{C}_{4}}{{r}^{2}} \right ), \end{align}

(4.15) \begin{align} & u_{r}^{i}=-2\cos \theta \left ( \frac {{{D}_{1}}}{{{r}^{3}}}+\frac {{{D}_{2}}}{r}+{{D}_{3}}+{{D}_{4}}{{r}^{2}} \right ), \end{align}

(4.16) \begin{align} & u_{\theta }^{i}=\sin \theta \left ( -\frac {{D}_{1}}{{r}^{3}}+\frac {{{D}_{2}}}{r}+2{{D}_{3}}+4{{D}_{4}}{{r}^{2}} \right ) . \end{align}

These factors $C$ and $D$ are determined by the velocity boundary conditions. We assume that the bubble has velocity $U\equiv -u_\infty$ ( $u_\infty \gt 0$ ). Taking the reference frame with the bubble, the far-field liquid velocity will be $u_r=u_\infty \cos {\theta }$ and $u_\theta =-u_\infty \sin {\theta }$ . Thus the far-field stream function will be

(4.17) \begin{equation} \psi ^{o}=-\tfrac {1}{2}u_\infty r^2 \sin ^2 \theta, \quad r\rightarrow \infty, \end{equation}

so $C_4=0$ and $C_3=- ({1}/{2})u_\infty$ using (4.11). For the flow inside the nanobubble, the velocity at $r=0$ has to be finite, so that $D_1=0$ and $D_2=0$ using (4.12). Furthermore, there are no masses across the bubble surface, and no bubble shape changes, so that $u_{r}^{o}=0$ and $u_{r}^{i}=0$ at $r=R$ . This leads to $C_1= (u_{\infty }R-2C_2 )R^2/2$ and $D_3=-R^2D_4$ .

For boundary conditions at the liquid–gas interface, the shear stress continuity is obeyed (Leal Reference Leal2007):

(4.18) \begin{equation} \tau _{r\theta }^{o}-\tau _{r\theta }^{i}=-{{\nabla }_{s}}\gamma. \end{equation}

Here, one can obtain $-{{\nabla }_{s}}\gamma =-({1}/{R})({\partial \gamma }/{\partial \theta })=- ({1}/{R}) ({\partial \gamma }/{\partial T}) ({\partial T}/{\partial \theta })= ({\partial \gamma }/{\partial T})$${(3{{k}^{o}}R+3k^ok^iG)}/{ (( 2{{k}^{o}}+{{k}^{i}} )R+2{{k}^{o})}{{k}^{i}}G}{\Gamma }\sin \theta$ , knowing the temperature distribution by (4.8). Later, we define $A=({\partial \gamma }/{\partial T}) {(3{{k}^{o}}R+3k^ok^iG)}/{ (( 2{{k}^{o}}+{{k}^{i}} )R+ 2{{k}^{o}}{{k}^{i}}G}{\Gamma )}$ such that we can simply write $-{{\nabla }_{s}}\gamma =A \sin \theta$ .

A number of recent works have shown that there exists slip at a fluid–fluid interface (Poesio et al. Reference Poesio, Damone and Matar2017; Telari et al. Reference Telari, Tinti and Giacomello2022; Hilaire et al. Reference Hilaire, Siboulet, Charton and Dufrêche2023). As seen from the simulation of the two-phase Couette flow in figure 6, this is indeed the case for the liquid–gas interface in our system. Thus it is natural to include slip in the theoretical model for the bubble velocity. However, how to model the fluid–fluid slip condition is not a well-explored topic. Even for the liquid–solid slip, the frequently used Navier’s slip condition (Zhang Reference Zhang2024), which assumes that the slip velocity is proportional to shear stress at the wall, is empirical (Leal Reference Leal2007). Thus at first, we still use the classic no-slip condition to derive the theoretical bubble velocity. Later, we then try to use empirical fluid–fluid slip models to see any improvements in the predictions.

The no-slip condition at the fluid–fluid interface is

(4.19) \begin{equation} u_{\theta }^o=u_{\theta }^i. \end{equation}

Using the shear stress condition (4.18) and the no-slip condition (4.19), one can obtain

(4.20) \begin{equation} {{D}_{4}}=\frac {4{{C}_{2}}-3{{u}_{\infty }}R}{4{{R}^{3}}}, \end{equation}

(4.21) \begin{equation} {{C}_{2}}=\frac {6{{\mu }^{o}}{{u}_{\infty }}{R}^{2}+9{{\mu }^{i}}{{u}_{\infty }}{{R}^{2}}-2A {{R}^{3}}}{12{{\mu }^{o}}R+12{{\mu }^{i}}R}. \end{equation}

Now the Marangoni force exerted on the bubble can be evaluated by integrating the shear stress:

(4.22) \begin{align} {{F}_{Ma}}&=-2\pi {{R}^{2}}{{\int _{0}^{\pi }{\tau _{r\theta } }}}\sin^{2}\theta\, {\rm d}\theta =-8\pi {{\mu }^{o}}\left ({{u}_{\infty }R-2{{C}_{2}}} \right ) \nonumber \\ & =-\frac {4\pi }{3}{{\mu }^{o}}\frac {2A{{R}^{2}}-3{{\mu }^{i}}{{u}_{\infty }}{{R}}}{{{\mu }^{o}}+{{\mu }^{i}}}. \end{align}

Similarly, by integrating the normal stress $\sigma _{rr}=-p+2\mu ({\partial u_r}/{ \partial r})$ , the pressure drag force experienced by the bubble is

(4.23) \begin{align} {{F}_{P}} &=2\pi {{R}^{2}}\int _{0}^{\pi }{\sigma _{rr} \sin \theta }\cos \theta\, {\rm d}\theta =8\pi \mu ^o \left ( {{u}_{\infty }}R-{{C}_{2}} \right ) \nonumber \\ & =2\pi \mu ^o \frac {2AR^2+6{{\mu }^{o}}{{u}_{\infty }} R+3{{\mu }^{i}}{{u}_{\infty }}{{R}}}{3{{\mu }^{o}}+3{{\mu }^{i}}}. \end{align}

One can see that the eventual expression of the pressure drag depends on $C_2$ , which again depends on the boundary condition at the liquid–gas interface. For a clean bubble with a negligible gas viscosity, $\tau _{r\theta }=0$ so that $C_2=u_{\infty }R/2$ and $F_P=4\pi \mu ^o u_{\infty }R$ , which recovers the classic drag on a conventional bubble. For our case, $C_2$ is different. By letting $F_{Ma}+F_{P}=0$ , one finds that $C_2=0$ . Using (C2), the velocity of the bubble ( $U\equiv -u_\infty$ ) is obtained as

(4.24) \begin{equation} U=-\frac {2R}{2{{\mu }^{o}}+3{{\mu }^{i}}}\,\frac {{{k}^{o}}R+k^ok^iG}{\left ( 2{{k}^{o}}+{{k}^{i}} \right )R+2{{k}^{o}}{{k}^{i}}G}\,\frac {\partial \gamma }{\partial T}\,{\Gamma }. \end{equation}

If there is no temperature jump ( $ G=0 $ ), then (4.24) reduces to the result of YGB59 theory (Young et al. Reference Young, Goldstein and Block1959):

(4.25) \begin{equation} U=-\frac {2R}{2{{\mu }^{o}}+3{{\mu }^{i}}}\,\frac {{{k}^{o}}}{2{{k}^{o}}+{{k}^{i}}}\,\frac {\partial \gamma }{\partial T}\,{\Gamma }. \end{equation}

In this case, the Marangoni force and drag force are

(4.26) \begin{equation} {{F}_{Ma}}=-{4\pi }{{\mu }^{o}}\frac {{{R}^{2}}\dfrac {\partial \gamma }{\partial T}\dfrac {2{{k}^{o}}}{2{{k}^{o}}+{{k}^{i}}}{{\Gamma }}-{{\mu }^{i}}{{u}_{\infty }}R}{{{\mu }^{o}}+{{\mu }^{i}}}, \end{equation}

(4.27) \begin{equation} {{F}_{P}}=2\pi {{\mu }^{o}}\frac {{{R}^{2}}\dfrac {\partial \gamma }{\partial T}\dfrac {2{{k}^{o}}}{2{{k}^{o}}+{{k}^{i}}}{{\Gamma }}+2{{\mu }^{o}}{{u}_{\infty }}R+{{\mu }^{i}}{{u}_{\infty }}R}{{{\mu }^{o}}+{{\mu }^{i}}}. \end{equation}

For macroscale gas bubbles whose viscosity and thermal conductivity are usually ignored, the velocity becomes

(4.28) \begin{equation} U=-\frac {R}{2{{\mu }^{o}}}\frac {\partial \gamma }{\partial T}{{\Gamma }}, \end{equation}

and the forces are ${{F}_{Ma}}=-4\pi ({\partial \gamma }/{\partial T}){{\Gamma }}{{R}^{2}}$ and ${{F}_{P}}=4\pi {{\mu }^{o}}{{u}_{\infty }}R+2\pi ({\partial \gamma }/{\partial T}){{\Gamma }}{{R}^{2}}$ . Thus the drag force is increased by the presence of the surface tension gradient. This extra term $2\pi ({\partial \gamma }/{\partial T}){{\Gamma }}{{R}^{2}}$ was not included in previous works about bubble motions in temperature gradients on the macroscale (Morick & Woermann Reference Morick and Woermann1993; Lubetkin Reference Lubetkin2003; Zeng et al. Reference Zeng, Chong, Wang, Diddens, Li, Detert, Zandvliet and Lohse2021). Reassuringly, if one investigates the total force ${{F}_{Ma}}+{{F}_{P}}$ on the bubble, then it is $4\pi {{\mu }^{o}}{{u}_{\infty }}R-2\pi ({\partial \gamma }/{\partial T}){{\Gamma }}{{R}^{2}}$ , which is the same for the previous study and the current study.

With all the measured transport properties of liquid, gas, and their interfaces, as summarized in table 3, now the theoretical velocity of nanobubbles without slip $U_{Nano-NS}$ (see (4.24)) can be evaluated. As shown in table 3, $U_{Nano-NS}$ values for the three cases agree reasonably well with the simulation results $U_{MD}$ . This highlights the consideration of the enhanced nanoscopic effects such as gas thermal conductivity, gas viscosity and liquid–gas thermal resistance for accurate descriptions of the nanobubble motion in a temperature gradient.

However, will the consideration of the slip improve the prediction even more? To answer this question, we have to use ‘empirical’ slip conditions for the fluid–fluid interface to include slip in the theoretical model of bubble velocity. In analogy to Navier’s slip condition for a liquid–solid interface, the slip condition for a fluid–fluid interface is proposed to be (Hilaire et al. Reference Hilaire, Siboulet, Charton and Dufrêche2023)

(4.29) \begin{equation} {{u}_{1}}-{u}_{2}=\frac {{{b}_{1}}}{{{\mu }_{1}}}\tau _{1}=\frac {{{b}_{2}}}{{{\mu }_{2}}}\tau _{2 }, \end{equation}

where ${u}_{1}$ , $b_1$ and $\tau _1$ ( ${u}_{2}$ , $b_2$ and $\tau _2$ ) are the interfacial velocity, slip length and interfacial shear stress of fluid 1 (fluid 2). For a clean surface where $\tau _1 = \tau _2$ , (4.29) means $b_{1}/b_{2}=\mu _{1}/\mu _{2}$ . However, in our case, $\tau _1 \neq \tau _2$ due to the extra Marangoni stress at the interface. Thus to account for this inconsistency, an extended slip model for the fluid–fluid interface may be (see Appendix B)

(4.30) \begin{equation} {{u_\theta }^{o}}-{{u_\theta }^{i}}=\frac {{{b}^{o}}}{{{\mu }^{o}}}\tau _{r\theta }^{o}+\frac {{{b}^{i}}}{{{\mu }^{i}}}\tau _{r\theta }^{i}. \end{equation}

Note that the definition of slip length in (4.30) is different from (4.29).

Using this boundary condition, the obtained bubble velocity with finite slip $U_{Nano-FS}$ is (see Appendix C)

(4.31) \begin{equation} U=-\frac {2R}{2{{\mu }^{0}}+3{{\mu }^{i}}\dfrac {1+2{{b}^{o}}/R}{1+3{{b}^{i}}/R}}\,\frac {{{k}^{o}}R+k^ok^iG}{\left ( 2{{k}^{o}}+{{k}^{i}} \right )R+2{{k}^{o}}{{k}^{i}}G}\,\frac {\partial \gamma }{\partial T}\,{\Gamma }. \end{equation}

One can see that the slip from the inner phase $b^i$ may speed up the bubble movement since the inner circulation is reduced by the slip. In contrast, the slip from the outer phase may slow down the bubble since the also reduced outer flow is the reason for bubble movement. However, using the slip lengths $b^i$ and $b^o$ measured from independent MD simulations (see table 3), the consideration of slip does not change the theoretical prediction too much in our case (see $U_{Nano-FS}$ in table 3) because $b/R$ is small. But a large liquid–liquid slip may also exist in other systems, such as at the interface between different types of polymer melts (Himmel & Wagner Reference Himmel and Wagner2013), so the consideration of slip may then become necessary.

Notably, the small difference between $U_{Nano}$ and $U_{MD}$ may be attributed to several reasons. First, the water viscosity shows a weak temperature dependence in the temperature range from 300 K to 350 K. However, we have used the mean viscosity in this temperature range for simplicity. Second, the Brownian motions of bubbles are strong in MD simulations, and only a limited number of realizations have been performed due to the high computational costs. Third, the solid plates may affect the flow field, though we have used the data of bubble displacements away from the plates.