## Abstract

Viscous liquids often exhibit flow slippage on solid walls. The occurrence of flow slippage has a large impact on the liquid transport and the resulting energy dissipation, which are crucial for many applications. It is natural to expect that slippage takes place to reduce the dissipation. However, (i) how the density fluctuation is affected by the presence of the wall and (ii) how slippage takes place through forming a gas layer remained elusive. Here, we report possible answers to these fundamental questions: (i) Density fluctuation is intrinsically enhanced near the wall even in a quiescent state irrespective of the property of wall, and (ii) it is the density dependence of the viscosity that destabilizes the system toward gas-layer formation under shear flow. Our scenario of shear-induced gas-phase formation provides a natural physical explanation for wall slippage of liquid flow, covering the slip length ranging from a microscopic (nanometers) to macroscopic (micrometers) scale.

## INTRODUCTION

The no-slip boundary condition of liquid flow is one of the most fundamental assumptions in fluid dynamics. However, there is no rigorous physical foundation for this condition. Therefore, its validation and violation have been discussed since the early days of the development of fluid mechanics (*1*–*8*). Recently, notable progress in experimental technology has unambiguously shown that the condition does not generally hold on a microscopic scale for a high enough shear rate. The widely used measure that characterizes the degree of the violation is the so-called slip length *b*, which is defined as the length at which *v*_{wall} = *b*(*∂v*/*∂y*)_{wall}, where *v* is the flow velocity, *y* axis is set to the direction normal to the wall, and the subscript “wall” means the value at the wall surface. For example, the slip length of water, which is the most typical and essential liquid and has thus been studied by many researchers, usually takes a value ranging from a few nanometers to 10 nm but sometimes reaches a few micrometers (*2*, *9*–*13*).

There have been a large number of studies on the physical mechanism behind such flow slippage. For example, the nature of fluid slip has been studied from a microscopic aspect by molecular dynamics (MD) simulations, which take into account microscopic interactions between the fluid and solid atoms (*14*, *15*). These studies have shown that a depletion layer of a few angstrom is formed between the fluid and the wall when the contact angle has a large value (e.g., a hydrophobic wall for water). Although a consistent scaling relationship has been established among the slip length, contact angle, and depletion layer, the physical mechanism behind flow slippage has remained elusive. For example, it is not clear how simple shear flow that is not coupled to the volume change can induce density fluctuations toward the formation of a gas layer. The shear rate used in MD simulations is typically 10^{9} to 10^{12} s^{−1}, which is higher than the one in experiments by many orders of magnitude. Furthermore, the slip length estimated in this way is far smaller than the typical experimental values, and thus, it is rather difficult to explain a slip length exceeding a few tens of nanometer. Although this approach is very transparent and attractive, thus, the relationship between the slippage phenomenon seen in microscopic MD simulations and that in macroscopic experiments is not so clear yet.

There is another type of theoretical approach, which assumed a priori that localized low-viscous gas layers (*16*–*20*) or nanobubbles (*21*–*23*) are formed on the wall. Once we accept this assumption, we may explain the apparent violation of the no-slip boundary condition, including the case of a large slip length. This scenario is rationalized by assuming the presence of nanobubbles on the wall (*8*, *20*, *24*). The existence of nanobubbles has been shown for a quiescent state by recent experimental studies [see (*8*) for review]. Consistent with this scenario, there has been experimental (*25*, *26*) and numerical (*27*) evidence that dissolved gases in liquids play a crucial role in flow slippage of simple liquids; for example, it was shown experimentally that the fluid slip sensitively depends on the amount and kind of a dissolved gas (*25*, *26*). We note that under ordinary experimental conditions, gas dissolution is inevitable. If the liquid is oversaturated with dissolved gas, gas bubbles will be spontaneously nucleated even in a quiescent condition if we wait for a sufficiently long time. That is, we may regard such a liquid to be metastable against gas formation.

The above scenario is attractive, and the surface enrichment of gas certainly takes place for hydrophobic surfaces. However, there have been debates on whether it is due to nanobubble formation or just the gas enrichment on the wall (*15*, *27*, *28*). Furthermore, there has so far been no satisfactory physical understanding of the impact of shear on nanobubble or gas-layer formation. For example, it is mostly unknown whether nanobubble formation under shear flow is the same as that in a quiescent state. The existence of stable nanobubbles was even paradoxical since a simple classical estimate suggests that they should dissolve quickly because of the large internal Laplace pressure coming from the large gas-liquid interface curvature (*8*). The stability of nanobubbles in a quiescent state has recently been explained by the pinning of their three-phase contact line. This stabilization mechanism suggests a possible fundamental difference between in a quiescent state and under shear. For example, it is not clear how the presence of a tiny amount of nanobubbles can lead to the formation of a rather thick gas layer on the wall. In addition, the formation of nanobubbles requires unique surface properties such as roughness and chemical inhomogeneity. This fact implies that the nucleation of the gas phase in bulk may suffer from a substantial barrier and is practically hard to take place. Then, the critical question is how shear flow affects the nucleation kinetics. This question is also related to another fundamental question of what determines the critical shear rate. The scenario based on nanobubbles seems to be difficult to explain the fact that flow slippage takes place even for walls of θ ≤ 90°.

Here, it is worth noting that slippage phenomena similar to those of simple liquids under flow were observed also for much more viscous fluids such as lubricants (*29*, *30*), for which direct experimental evidence for shear-induced cavitation (or bubble formation) has been shown. The results of these two types of systems, simple nonviscous liquids and highly viscous liquids, apparently look quite different; at the same time, however, there is a commonality that both phenomena involve gas-phase formation. Thus, it is natural to consider that there may be a common physical mechanism behind the two phenomena. This motivates us to explain the flow slippage phenomena based on shear-induced cavitation or gas-layer formation on solid surfaces in a unified manner.

Here, we propose a novel physical mechanism of flow slippage based on the shear-induced liquid-to-gas transition in a liquid metastable against gas formation, taking into account the boundary conditions, including the friendliness (wettability) of the gas phase to the wall. To this end, we numerically study the mechanism of gas generation on a flat solid wall in a viscous metastable liquid, which mimics a liquid oversaturated with a gas. We assume that a metastable liquid obeys the mass conservation and compressible Navier-Stokes equations with thermal fluctuations and satisfies the no-slip boundary condition on a flat solid wall. We note that, in this thermal system, the thermally activated nucleation of gas bubbles should occur for a metastable liquid even in a quiescent state if we wait for a sufficiently long time. This incubation time can be enormously long, and thus, nucleation practically does not happen besides in the form of nanobubbles. This fact implies that there should be some mechanism of shear-induced reduction of the nucleation barrier. We investigate how the nucleation barrier for gas-bubble formation is reduced by shear flow and how shear-induced gas-bubble formation is affected by the presence of the solid wall and its surface wettability. Our study demonstrates that the amplitude of density fluctuation increases near the wall, and the application of shear further enhances the fluctuation and helps the formation of a gas layer on the wall. This indicates that the presence of nanobubbles is a sufficient condition to have flow slippage for low-viscosity liquids, such as water, but may not be a necessary condition. For very viscous liquids, for example, nanobubbles may be unnecessary, and spontaneous shear-induced slippage can take place even at a low shear rate.

## RESULTS AND DISCUSSIONS

Here, we consider a van der Waals liquid, whose viscosity η weakly depends on the density ρ, near a gas-liquid transition as a model, and use a general coarse-grained model of compressible liquid in the framework of fluctuating hydrodynamics (*31*) (see Materials and Methods on the details). Unless otherwise stated, we set the homogeneous initial normalized density of the fluid as ψ_{0} = 0.480, which corresponds to a thermodynamically metastable state (see the Supplementary Materials for the van der Waals phase diagram). We introduce a solid wall, which has a contact angle θ with a liquid, i.e., the angle within the body of the liquid formed at the gas-liquid-solid interface, and apply simple shear, whose shear rate is

### Enhancement of density fluctuations near walls in a quiescent condition (at γ ̇ = 0 )

It is well known that the average density of ψ_{0}(*y*) near the wall increases or decreases compared to the average density in bulk, depending on the contact angle θ. Here, the average density ψ_{0}(*y*) is defined as* _{x}* and 〈⋯〉

*are the averaging operations in terms of the coordinate*

_{t}*x*(along the wall) and the time

*t*, respectively. Such behavior can be explained by the standard density functional theory (

*32*). We show in the inset of Fig. 1A the spatial profile of the average density of ψ

_{0}(

*y*) for various values of the contact angle θ. We can see that the density profile near the wall is crucially dependent on the contact angle θ, as expected. The liquid density near the wall increases (decreases) for θ ≤ 90° (θ ≥ 90°).

Besides the well-known dependence of the density profile near the wall on θ mentioned above, we find a nontrivial influence of the wall on the spatial profile of the amplitude of density fluctuation, δψ(*y*) in a quiescent condition, i.e., for the shear rate *y*) on the contact angle θ. Here, the density fluctuation δψ(*y*) are defined as

We can see that the amplitude of density fluctuation at the wall is larger than that in bulk for any contact angle. As shown in Fig. 1, this enhancement can be seen even when θ = 90°, i.e., for the case of no preferential wetting. This behavior is unexpected since there is no fluid-wall coupling in the free-energy functional for θ = 90°. This implies that the enhancement should arise solely from the fluid-solid boundary condition itself.

Here, we explain the physical mechanism behind this wall-induced enhancement of density fluctuation. The equilibrium amplitude of thermal density fluctuation is determined by the competition between the strength of thermal fluctuation and the compressive resilience against the fluctuation. Now, we note that the boundary condition sets the influx and outflux through the boundary wall to be zero. Thus, the boundary condition restricts the direction of the compressive resilience only along the wall, leading to the anisotropy in the compressibility. This anisotropy effectively decreases the compressive resilience, resulting in the enhancement of density fluctuation near the wall than in bulk.

Since this enhancement mechanism is just a consequence of the constraint that the density flux is zero through the wall, we can consider this problem in the framework of model B (*33*, *34*), where diffusion is the only transport mechanism, without losing the generality. By solving the one-dimensional problem, we can derive the density fluctuation near the wall analytically for θ = 90° as follows (see the Supplementary Materials for the derivation)_{bulk} is the equilibrium value of thermal fluctuation in bulk. Therefore, the amplitude of density fluctuation should be theoretically

Here, we emphasize that this enhancement of density fluctuation near the wall plays an essential role in shear-induced destabilization of a liquid against the gas-phase formation. This effect may be crucial to understand why the flow slippage takes place near the wall even for θ ≤ 90°.

### Nontrivial effects of shear flow on density fluctuation

Now, we consider the impact of this enhancement of density fluctuation near the wall on shear-induced cavitation. For a metastable liquid, we should eventually observe inhomogeneous nucleation of gas bubbles via a thermal activation process if we wait for a sufficiently long time. According to the classical nucleation theory, the probability of the gas-phase nucleation *p*_{N} is estimated as *p*_{N} ∝ exp [− Δ*G**/*k*_{B}*T*], where Δ*G** is the height of the free-energy barrier for nucleation. Considering the above boundary effect, which increases density fluctuation near the wall, it is natural to expect that the probability of inhomogeneous nucleation on the solid surface should be larger than the above bulk value. Here, we emphasize that this conclusion should hold even when θ ≤ 90°.

Keeping this in mind, we consider the effect of shear flow on density fluctuation. Usually, it is believed that shear flow is not coupled to density fluctuation simply because simple shear deformation does not cause any volume deformation. However, it was shown recently that this is limited to the case when the viscosity does not depend on the density (*31*, *35*). Since the viscosity of a liquid generally depends on its density, shear flow can enhance density fluctuation through modification of the effective compressibility along compression and extension directions of shear deformation (*31*). When the imposed shear rate *P* being the pressure), the sign of the effective compressibility becomes negative, and the liquid intrinsically becomes unstable against gas-phase formation. For *G**. Thus, inhomogeneous nucleation is accelerated by this flow-induced instability mechanism, and gas-bubble formation should occur on the wall much easier under shear flow than in a quiescent state. We expect that the application of shear can induce the nucleation of a much larger amount of the gas phase than nanobubbles formed in a quiescent state with the help of unique surface properties of the wall.

For a shear rate of

Here, we note that in this simulation, thermal fluctuation is switched on at *t* = 0. Thus, the time development of density fluctuation should be a superposition of the thermal equilibration process and the exponential growth process of density fluctuation. The former process can be obtained by the time evolution of the fluctuation *31*), one can show that the amplification factor should satisfy the relation

### Nucleation of gas bubbles accelerated by shear flow

As discussed above, even without shear flow, gas bubbles will be formed by the thermal activation process for a metastable state if we wait for a sufficiently long time. So, we need to consider nucleation-growth–type gas formation in addition to spinodal decomposition–type one. We show in Fig. 3A a typical gas-bubble generation process observed in a metastable state. We can see that the system remains spatially homogeneous in the metastable state until gas-bubble generation happens (see left panel). At a particular moment, a gas-bubble nucleus appears and rapidly grows in its size. When a gas bubble is nucleated, it often appears only on one of the two walls in our simulations since the full growth of a bubble takes place before the next bubble nucleation event takes place in our simulation box size. We define this birth time of the first gas bubble as the “incubation time.” The dependence of the incubation time on the shear rate and the contact angle is shown in Fig. 3B. We can see that the incubation time for bubble formation on the walls decreases with an increase in the contact angle and/or the shear rate.

Here, we also show in Fig. 3C how the incubation time depends on the shear rate of *t*_{inc} decays monotonically with an increase in

### Behavior of the gas phase under shear flow

If we stop shearing immediately after the nucleation of the gas phase and wait for a sufficient time to reach the thermodynamic equilibrium, the gas phase takes a stable equilibrium shape on the wall, depending on the contact angle. Under a steady shear flow, however, the structure of the gas phase is distorted from its equilibrium shape by the Kelvin-Helmholtz (KH) instability mechanism, which is the hydrodynamic instability occurring when two fluid layers with different densities move with the different horizontal speeds (*36*).

We investigate the temporal change of the shape of the gas phase under shear flow by performing the following simulations. We prepare the thermal equilibrium state of the gas phase coexisting with the liquid phase on the wall as the initial state. Then, we apply shear flow at *t* = 0 to investigate the temporal change of the gas-phase shape induced by shear flow.

*Time evolution of the gas-phase shape under shear*. First, we show an example of the time evolution of the gas-phase shape under shear flow in Fig. 4A for

Next, we show the case for the contact angle θ = 90° and

Last, we show the case for a large contact angle, θ = 150. 8° and

Thus, for a high shear rate regime, the gas-phase structure is strongly perturbed and detached from the wall by shear flow. For the realization of a large slip length, the presence of a sufficient amount of the gas phase on the wall is necessary. For a small contact angle, the detachment of the gas phase takes place easily, even for weak shear. We summarize the critical shear rate for detachment or large temporal fluctuations in Fig. 4F. The black dashed curve separating the steady and intermittent regions in Fig. 3C is drawn on the basis of the critical shear rate plotted in Fig. 4F. In the steady state, there is a well-defined stable slip length, whereas in the intermittent state, a slip length is not well defined and fluctuating with time as indicated in Fig. 2D.

*Theoretical estimation of the critical shear rate of the KH instability*. Here, we consider a mechanism by which a two-dimensional gas layer formed on a wall is destabilized under shear flow. This mechanism is known as the KH instability. We consider a situation that the gas and liquid phases coexist with the flat interface parallel to the wall located at *y* = 0. We also assume that the gas and liquid layers exist in the upper half plane and lower half plane, respectively. Then, we apply a simple shear flow, whose spatially averaged velocity field is given by *ik _{x}*(

*x*−

*ct*)].

This problem has already been discussed in several seminal papers (*37*–*40*). Here, we follow the discussion in (*39*) and estimate the critical shear rate _{LG} tends to suppress the vertical perturbation with the characteristic decay rate of ω_{decay} ∼ *O*(σ_{LG}*k _{x}*/η). These two competing effects are balanced when ω

_{growth}∼ ω

_{decay}. Thus, we obtain the critical shear rate for the KH instability as follows

Now, we calculate the actual value of the critical shear rate by substituting the parameters used in our simulations. The surface tension is estimated as σ_{LG} = 0.0029 by using the standard density functional theory (*41*). By setting the characteristic perturbation wave number as *k _{x}* = 2π/

*L*, the critical shear rate is estimated as

_{x}### Stability of gas bubbles formed on walls in a metastable state

*Destabilization of gas bubbles and the resulting intermittent behavior*. Destabilization of a gas phase by shear takes place not only for a flat gas layer but also for gas bubbles. For example, gas bubbles are deformed and split into small ones by shear. Thus, we consider the stability of gas bubbles generated in a metastable state under shear flow. If the generation and growth of a gas bubble occur in a quiescent condition, the bubble should be stable and exist as it is. This is simply because the final state is in thermal equilibrium. Under a steady shear flow, however, after the formation of gas bubbles, they are repeatedly broken up and coalesce unless the shear rate is very low

*Origin of the intermittent behavior*. To understand the mechanism of the intermittent behavior, we study the dynamical behavior of the gas phase under steady shear flow. We show an example of the time evolution of the density field for *34*, *42*), which leads to the emergence of tearing shear stress due to the difference in the viscosity between the liquid and gas, and the other is the advection due to the advective term in the Navier-Stokes equation. It is well known that the advective term stirs the system; however, this nonlinear term may not be so crucial because the Reynolds number is very small in our simulations. Once the tearing shear stress splits the gas bubble to pieces smaller than the critical gas-bubble nucleus size, they dissolve, and accordingly, the system goes back to a homogeneous state. Detachment of a gas bubble from the wall, which more easily happens for the case of lower gas wettability to the wall, also enhances their dissolution. It is the origin of the intermittent behavior. Here, we note that, in the intermittent state, the slip length *b* also fluctuates with time, as shown in Fig. 2D. When the wettability of the gas phase to the wall is high enough and the shear rate is low enough, the stabilizing force coming from the solid- and liquid-gas interfacial energy wins over the destabilizing force due to Taylor mechanism, leading to the emergence of a steady state. Thus, we can conclude that the formation of gas bubbles and their survival under shear is both easier under a lower shear rate and for a larger contact angle. In this steady state, the slip length *b* does not fluctuate largely with time.

### Slip length

Next, we estimate the slip length for shear-induced flow slippage described above.

*Slip in the steady state*. In Fig. 5A, we show the *y* dependence of the shear velocity *v _{x}*(

*y*) = 〈

*v*(

_{x}*x*,

*y*,

*t*

_{0})〉

*and the effective shear viscosity η(*

_{x}*y*) = 〈η(

*x*,

*y*,

*t*

_{0})〉

*at a certain moment*

_{x}*t*=

*t*

_{0}when the steady state is reached in the metastable region after gas-bubble generation and growth. The decrease in the averaged effective shear viscosity near the wall indicates the nonlinear shear velocity profile

*v*(

_{x}*y*) near the wall. We can see the apparent slip behavior by looking at the deviation of the profile from the red dashed line in Fig. 5A. We also show the dependence of the time-averaged apparent slip length on the shear rate and the contact angle in Fig. 5B. The slip length depends on both the shear rate and the contact angle, although its dependence is rather weak. Here, we take the time averaging operation only for the state with gas bubbles near the wall. We note that the way of this time average has a considerable effect on the estimation of slippage for the intermittent state. Although the slip length itself rather smoothly depends on the shear rate and the contact angle, the amplitude of the slip length fluctuations shows a distinct transition across the border between the steady and intermittent states (see Fig. 2C).

*Slip in the unstable regime*. Here, we note that in the unstable region, the flow slippage also takes place, but the slip length strongly fluctuates. We investigated the time dependence of the slip length and density distribution for the spinodal region. We select the shear rate and contact angle as

*System size effect on the slip length*. Until now, we have shown the results of the slip length for the system size of *L _{x}* =

*L*= 128. However, we find that when we change the system size while keeping the box shape, the apparent slip length is proportional to the size of the system. The reason is that the ratio of the total volume of generated gas bubbles localized near the wall to that of the system is determined by the lever rule in a nonequilibrium steady state. Therefore, even if the metastability of the liquid is weak, the slip length can become large for a system of large size. It may explain some experimental reports of large slip lengths (∼a few micrometers). The validity of our mechanism may be checked experimentally by examining the system size dependence of the slip length. It should be noted that the slip length should also depend on the amount of dissolved gas. This may explain a broad distribution of the slip lengths reported even for the same liquid-substrate pairs.

_{y}### Comparison of physical parameters between simulation and real systems

Here, we make a connection of our simulation results to experimental observation.

*Spinodal decomposition–type gas-phase formation*. As in many other simulations, there is a large gap in the time scale between the simulations and experiments. We estimate the actual physical values of the viscosity coefficient η_{c}, unit length *l*, and unit time τ used in our simulation by substituting typical parameters for liquids like water. In our simulation, we impose the following constraints: _{c}v_{0})/(5*k*_{B}*T*) between numerical parameters. Here, η_{c} is the viscosity, *m* is the molecular mass, *v*_{0} is the molecular volume, and *T* is the temperature. By using these relations and setting *T* = 300 K, *m* = 3.0 × 10^{−26} kg, and *v*_{0} = 3.0 × 10^{−29} m^{3}, we obtain

In particular, the viscosity coefficient η and the critical shear rate

On the other hand, it is well known that the experimental value of the viscosity coefficient for pure water at ambient condition is about 1 mPa·s. The value of *43*–*45*): ^{−1} at ambient condition.

There are not many experimental studies on the viscosity of water saturated with gas, but there are several studies on the aqueous solution of carbon dioxide (*46*–*48*). The value of ^{12} s^{−1} (*46*, *48*). For a low pressure range (from 0.7 to 5 MPa), on the other hand, the value of ^{10} s^{−1} (*47*), although large experimental error prevents accurate estimation. Anyway,

The time unit in our simulations is ∼12 ps, and the unit of the shear rate is thus ∼8 × 10^{10} s^{−1}. So, the shear rate in our simulations is quite high compared to experiments that reported flow slippage, which is comparable to MD simulations [see, e.g., (*15*)]. In any case, these estimations indicate that the spinodal decomposition–type gas-phase formation is not expected for low viscosity liquids such as water in a realistic experimental condition.

However, this is not necessarily the case for very viscous liquids near the glass transition (*31*, *49*). Now, we discuss the value of the unit time τ for highly viscous liquids. In simulations, we investigated the instability mechanism by solving the Navier-Stokes equation under a low Reynolds number condition. However, it is expected that the same conclusion should be obtained even under the Stokes approximation because the inertia term in the Navier-Stokes equation does not play an essential role in our instability mechanism. Under the Stokes approximation, it is unnecessary to assume the relation *29*, *30*). Surface effects on such instability is a fascinating problem for future experimental investigation. Flow slip in these viscous liquids should lead to a drastic reduction of energy dissipation in their hydrodynamic transport, which may have a considerable impact on various applications.

*Nucleation-growth–type gas-phase formation*. For nucleation-growth–type gas-phase formation, the story is quite different, since nanobubbles may exist even in a quiescent state without shear. For example, it is widely known that a liquid with dissolved gas, which is used in most of real experiments, tends to form nanobubbles on the wall if it is wettable to the gas phase (e.g., a hydrophobic wall for water) in a quiescent state (*8*). This fact seems to suggest that the gas layer may be easily formed by shear, which leads to the flow slippage [see, e.g., (*8*, *20*, *24*)].

The metastability should lead to the formation of such nanobubbles if we wait for a long enough time. Thus, the situation looks obvious, and no further consideration is necessary to understand the flow slippage. However, this is not the case. We stress that the existence of stable nanobubbles itself is not so obviously justified since a simple classical estimate suggests that they should dissolve quickly due to the large internal Laplace pressure coming from the large gas-liquid interface curvature (*8*). The stability of nanobubbles in a quiescent state has recently been explained by the pinning of their three-phase contact line. This stabilization mechanism suggests a possible fundamental difference between in a quiescent state and under shear. For example, it is not clear at all how the presence of a tiny amount of nanobubbles can lead to the formation of a rather thick gas layer on the wall. Thus, for example, we do not know the amount of the gas phase to be formed under shear. Furthermore, we should note that the formation of nanobubbles requires unique surface properties such as surface roughness and chemical inhomogeneity. This fact means that the actual free-energy barrier for gas bubble nucleation must be extremely high, and the nucleation does not happen in a realistic time scale besides nanobubbles. Even in such a situation, the shear flow should play a crucial role in gas-phase formation and the resulting slippage. Our scenario provides a simple physical picture on the impact of such shear effects on the formation of gas bubbles or gas layers.

As we show in Fig. 3 (B and C), there is a steep reduction of the incubation time of gas-bubble nucleation as a function of the shear rate, which is particularly the case for θ ∼ 90°. If we use the parameters for water, the time unit in this plot approximately corresponds to τ ∼ 10^{−12} s, and the unit of the shear rate is ∼10^{9} s^{−1} for Fig. 3 (B and C). Thus, for example, for θ ∼ 140°, we estimate the incubation time of spontaneous formation of gas bubbles on the walls to be about 0.01 to 0.1 s, although this estimation is very rough. For lower θ, this incubation time should become even much longer and may exceed the experimental time scale. The crucial point is that our scenario can explain the occurrence of the flow slippage even for θ < 90°, despite that the gas phase is less wettable for the wall in this condition. According to Fig. 3 (B and C), the dependence of the incubation time on the shear rate is very steep, particularly for a smaller contact angle. Thus, we expect that the shear-induced cavitation takes place even for a very small shear rate used in experiments within a reasonable observation period. It means that there may practically be an onset shear rate even for a metastable state [see, e.g., (*26*)]. Since the amount of a gas phase under shear is not constrained by the pinning of the three-phase contact line, it should simply follow the lever rule. Thus, it can even explain a macroscopic slip length if a sample size is large. Thus, we may say that the most natural scenario of flow slippage in a low-viscosity liquid may be nucleation-growth–type shear-induced cavitation. Here, we note that the density dependence of the viscosity η or the functional form of η(ρ), affects the strength of the shear effects on the gas-phase formation, as shown in the Supplementary Materials.

Furthermore, the state of the gas phase nucleated under shear flow is not necessarily in a steady state in the form of gas layers on the wall. The interaction of the gas phase formed on the wall with the shear flow can lead to complex behavior, as shown in Figs. 2 (C and D) and 4. The intermittent behavior should be observed in a broad region of the contact angle shear rate state diagram, which may be a cause of considerable fluctuations of the slip length observed in experiments (*12*, *13*).

## CONCLUSIONS

To summarize, we propose a novel physical mechanism of flow slippage originating from intrinsic enhancement of density fluctuation near solid surfaces due to the boundary condition and the coupling between shear stress and density through the density dependence of the viscosity η(ρ). Our theory provides essential information on how we can control the flow slippage, including the case of a metastable liquid saturated by gas, which we mainly discuss in this article. However, the relevance of our scenario to real experiments should be carefully examined. We note that the information on the amount of dissolved gas is crucial for such a comparison. The comparable shear rate required for the formation of slippage between our coarse-grained simulations and MD simulations indicates that, without dissolved gas, the experimentally unrealistic and exceptionally high shear rate is necessary for shear-induced gas-layer formation. We have also shown that this mechanism should also be essential for a viscous liquid near its glass transition point since the density dependence of the viscosity is very strong there (see the preceding section). In such a case, even spinodal decomposition–type shear-induced cavitation can take place preferentially on solid walls.

Our mechanism can naturally explain experimental reports of rather large slip lengths, which have so far been difficult to explain theoretically. For example, the shear rate–dependent slip lengths reported in (*11*) can be rather naturally understood by our mechanism, including complex, nonideal behaviors. It is known that the large-length slip experimental results largely fluctuate with time and depend on the spatial position of their measurements (*50*). Such features have been quite difficult to be explained by previous models. According to our results, the shape of gas bubbles near the wall is neither uniform nor steady in the metastable intermittent state and fluctuates spatiotemporally, even for a flat homogeneous solid wall. In the above, we have discussed the gas generation process on a completely flat wall. However, actual solid walls must have surface roughness and chemical inhomogeneity. Impact of such inhomogeneities on flow slippage phenomena and their application to the control of slippage (*50*–*53*) remains an interesting problem for the future investigation.

In this work, we focus on the slippage of a simple viscous fluid, in which the flow-induced instability stems from the density dependence of the viscosity. It has been shown that complex fluids also exhibit similar flow-induced instability in bulk (*54*), which resembles the shear-induced instability (*55*–*57*) of dynamically asymmetric complex fluids, which is a mixture of large (slow) and small (fast) components (*58*). The shear-induced instability in such complex fluids and simple viscoelastic liquids can be understood in a unified manner by using the mapping of *c* ↔ ρ and Π ↔ *p* (*31*), where *c* is the concentration and Π is the osmotic pressure in the complex fluids. However, rather a few studies have been conducted on the shear-induced instability of complex fluids influenced by a boundary wall (*59*–*67*). In this case, the liquid-rich dilute phase is to be formed near the walls for such a mixture to reduce viscous dissipation, as the gas phase is formed for a simple liquid. We believe that the instability mechanism described here may be applied commonly to shear-induced instability of viscoelastic materials and complex fluids under the influence of walls.

## MATERIALS AND METHODS

### Simulation method

Here, we consider a simple viscous liquid by using the following general model of fluid dynamics (*31*). The time evolution of density ρ(*r*, *t*) is given by the continuity equation*v _{i}* is the fluid velocity. The time evolution of fluid velocity is given by the momentum conservation equation, i.e., the Navier-Stokes equation

*are the viscous stress, thermal noise stress, and pressure tensors, respectively. The viscous stress tensor*

_{ij}Here, we assume that the ρ dependence of the shear viscosity η(ρ) and the bulk viscosity ζ(ρ) as η(ρ) = η* _{d}* · [1 + 6(ρ

*v*

_{0}/

*m*) + 36(ρ

*v*

_{0}/

*m*)

^{2}]/7 and ζ(ρ) = 8η(ρ)/15, respectively. To see the effect of the functional form of the density dependence of η(ρ), we also use another functional form (Supplementary Materials). The thermal noise stress

We use the Ginzburg-Landau–type free energy functional for the liquid, including the fluid-wall interaction

We use the van der Waals formula *f*(ρ) = − ε*v*_{0}(ρ/*m*)^{2} − (ρ/*m*)*k*_{B}*T* ln (*m*/ρ*v*_{0} − 1) as the bulk free-energy density *f*(ρ), which allows us to describe the gas-liquid phase transition properly. Here, *m*, *v*_{0}, and ε are the molecular mass, excluded volume, and strength of attractive interaction, respectively. We use the van der Waals formula as a simple model equation that can describe the transition between the gas and liquid phases. We set the fluid-wall surface interaction as ϕ_{S}(ρ) = *k*ρ, where *k* is a constant. The pressure tensor Π* _{ij}* is obtained by the functional derivative of ℱ. We use the normalized density, unit length, and unit time as ψ(

*r*,

*t*) = ρ(

*r*,

*t*)

*v*

_{0}/

*m*,

*l*= (

*m*

^{2}

*C*/

*v*

_{0}

*k*

_{B}

*T*)

^{1/2}, and τ = η

_{c}

*v*

_{0}/(5

*k*

_{B}

*T*), respectively. The velocity is scaled by

*l*/τ. Furthermore, we set

*31*).

We consider two-dimensional liquids (the spatial dimension *d* = 2) in this work, but the result should be essentially the same for *d* = 3. The system has a square shape whose horizontal and vertical lengths are *L _{x}* =

*L*= 128. The system is sandwiched by two parallel flat solid walls, and the flow direction is the

_{y}*x*direction. We impose a periodic boundary condition to

*x*direction. The contact angle θ between the gas phase, the liquid phase, and the solid wall is analytically calculated from the density functional theory (

*32*,

*41*). The velocity of the top and bottom walls are given by

*68*).

When we study the stability of gas bubbles formed on the wall under shear, we calculate the time evolution of the system while switching off the thermal noise after the gas bubble generation under thermal fluctuation. We detect the emergence of gas bubbles from a metastable liquid after a sufficiently long time, which is defined as *t* = 0. The incubation time of gas bubbles is defined as follows. First, we define the “gas phase” as lattices having a density less than ψ = 0.247. We calculate the time evolution of the system until the total area fraction of the gas phase exceeds 3%. If the area fraction exceeds 3%, then we go back in time and define the moment that corresponds to the gas area fraction equal to 0.03% as “gas bubble generation time” *t* = 0. We note that our results are not sensitive to the selection of the above thresholds, 3 and 0.03%.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/13/eaaz0504/DC1

Supplementary Text

Fig. S1. The van der Waals equation of state used in the main text.

Fig. S2. Normalized density dependence of the viscosity coefficient.

Fig. S3. Shear rate dependence of the incubation time for viscosity functions.

Movie S1. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0001_theta = 27.6).

Movie S2. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0016_theta = 90.0).

Movie S3. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0040_theta = 90.0).

Movie S4. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0016_theta = 150.8).

Movie S5. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0040_theta = 150.8).

Reference (*69*)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank A. Furukawa for valuable discussions and suggestions.

**Funding:**This work was partially supported by Grants-in-Aid for Specially Promoted Research (JP25000002) and Scientific Research (A) (JP18H03675) from the Japan Society of the Promotion of Science (JSPS) and by the Mitsubishi Foundation for support. Y.K. acknowledges the JSPS for financial support (JP15J00018).

**Author contributions:**H.T. conceived and supervised the project. Y.K. performed the simulations. Y.K. and H.T. discussed the results and wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).