text:
# From C*algebra extensions to CQMS, 𝑆𝑈_𝑞(2), Podles sphere and other examples

## 1. Introduction

In noncommutative geometry, the natural way to specify a metric is by a “Lipschitz seminorm”. This idea was first suggested by Connes (), and developed further in . Connes pointed out (, ) that from a Lipschitz seminorm one obtains in a simple way an ordinary metric on the state space of a $`C^{}`$-algebra. A natural question in this context is when does this metric topology coincides with the weak\* topology. In his search for an answer to this question, Rieffel (, , ) has identified a larger class of spaces, namely order unit spaces on which one can answer these questions. He has introduced the concept of Compact Quantum Metric Spaces (CQMS) as a generalization of compact metric spaces, and used () this new concept for rigorous study of convergence questions of algebras much in the spirit of Gromov-Hausdorff convergence. One natural question in this regard is are there plenty of CQMS floating around? Rieffel (, ) has given some general principles for the construction of CQMS. In  we exploited one of his principles to construct CQMS. In fact Rieffel has shown () that indeed there are many examples. But in concrete C\*-algebras one would like to have more explicit description of these structures. Our objective here is construction of CQMS out of quantum SU(2) and Podles spheres. To achieve that we make a slightly general construction and produce CQMS starting from C\*-algebra extensions. Organization of the paper is as follows. In the next section we recall the basics of CQMS. In section 3 the basic construction has been described. In the final section we employ the principle developed in section 3 to special cases.

## 2.  Compact Quantum Metric Space: Preliminaries

We recall some of the definitions from .

###### Definition 2.1.

An order unit space is a real partially ordered vector space $`A`$ with a distinguished element $`e`$, the order unit satisfying
(i) ( Order Unit property ) For each $`aA`$ there is an $`r`$ such that $`are`$.
(ii) ( The Archimidean property ) If $`aA`$ and if $`are`$ for all $`r`$ with $`r0`$, then $`a0`$.

###### Remark 2.2.

The following prescription defines a norm on an order unit space.

$$a=inf\{r|reare\}$$

###### Definition 2.3.

By a state of an order unit space $`(A,e)`$ we mean a $`\mu A^{}`$, the dual of $`(A,)`$ such that $`\mu (e)=1=\mu ^{}`$. Here $`^{}`$ stands for the dual norm on $`A^{}`$. Collection of states on $`(A,e)`$ is denoted by $`S(A)`$.

###### Remark 2.4.

States are automatically positive.

###### Example 2.5.

Motivating example of the above concept is the real subspace of selfadjoint elements in a C\*-algebra with the order structure inherited from the C\*-algebra.

###### Definition 2.6.

Let $`(A,e)`$ be an order unit space. By a Lip norm on $`A`$ we mean a seminorm $`L,`$ on $`A`$ such that
(i) For $`aA`$, we have $`L(a)=0`$ iff $`ae`$;
(ii) The topology on $`S(A)`$ coming from the metric

$$\rho _L(\mu ,\nu )=sup\{|\mu (a)\nu (a)|:L(a)1\}$$

is the $`w^{}`$ topology.

###### Definition 2.7.

A compact quantum metric space is a pair $`(A,L)`$ consisting of an order unit space $`A`$ and a Lip norm $`L`$ defined on it.

The following theorem of Rieffel will be of crucial importance.

###### Theorem 2.8 (Theorem 4.5 of ).

Let $`L`$ be a seminorm on the order unit space $`A`$ such that $`L(a)=0`$ iff $`ae`$. Then $`\rho _L`$ gives $`S(A)`$ the $`w^{}`$-topology exactly if
(i) $`(A,L)`$ has finite radius, i.e, $`\rho _L(\mu ,\nu )C\text{ for all }\mu ,\nu S(A)`$ for some constant $`C`$, and
(ii) $`_1=\{a|L(a)1,\text{ and }a1\}`$ is totally bounded in $`A`$ for $``$.

## 3. Extensions to CQMS

In this section we describe the general principle of construction of CQMS from certain $`C^{}`$-algebra extensions. Let $`𝒜`$ be a unital $`C^{}`$-algebra. Fix a faithful representation $`𝒜()`$. Suppose we have a dense order unit space $`Lip(𝒜)𝒜_{s.a}`$, where $`𝒜_{s.a}`$ denotes the real partially ordered subset of selfadjoint elements in $`𝒜`$. Let $`L`$ be a Lip norm on $`Lip(𝒜)`$ such that $`((Lip(𝒜),I),L)`$ is a CQMS. Let $`\nu `$ be a state on $`𝒜`$, then define $`\stackrel{~}{𝒜_\nu }`$ to be the collection of $`((a_{ij}))𝒦(l^2())𝒜`$ such that (i) $`a_{ij}Lip(𝒜)`$, (ii) $`a_{ij}=a_{ji}`$, (iii) $`sup_{i1,j1}(i+j)^k(L(a_{ij})+|\nu (a_{ij})|)<\mathrm{}k`$. Clearly $`𝒜_\nu :=\stackrel{~}{𝒜_\nu }I`$, where $`I`$ is the identity on $`(ł^2())`$ is an order unit space. Define $`L_k:𝒜_\nu _+`$ by $`L_k(I)=0,`$

$$L_k((a_{ij}))=sup_{i1,j1}(i+j)^k(L(a_{ij})+|\nu (a_{ij})|).$$

###### Lemma 3.1.

Let $`d=`$ diameter of $`((Lip(𝒜),I),L)`$. Then for a “Lipschitz function” $`aLip(𝒜)`$ one has $`a(L(a)+|\nu (a)|)(1+d).`$

Proof: Let $`\mu `$ be an arbitrary state on $`𝒜`$. Then using $`sup\{|\mu (a)\nu (a)|:L(a)1\}d`$ we get,

$`|\mu (a)|`$ $``$ $`|\mu (a)\nu (a)|+|\nu (a)|`$
$``$ $`L(a)d+|\nu (a)|`$
$``$ $`(L(a)+|\nu (a)|)(1+d).`$
$`\mathrm{}`$

###### Lemma 3.2.

There exists a constant $`C>0`$ such that for $`((a_{ij}))\stackrel{~}{𝒜_\nu },`$

$$((a_{ij}))CL_2((a_{ij})).$$

Proof: Let $`\{e_i\}_{i1}`$ be the canonical orthonormal basis for $`l^2()`$. Let $`\lambda _ie_iu_i`$ and $`\mu _ie_iv_i`$ be two generic elements in $`l^2()`$. Here $`u_i,v_i`$ are unit vectors. Then clearly $`\lambda _ie_iu_i^2=|\lambda _i|^2,\mu _ie_iu_i^2=|\mu _i|^2.`$ Now observe that

$`|{\displaystyle \lambda _ie_iu_i},((a_{ij})){\displaystyle \mu _je_jv_j}|`$ $``$ $`{\displaystyle |\lambda _i||\mu _j||u_i,a_{ij}v_j|}`$
$``$ $`{\displaystyle |\lambda _i||\mu _j|(L(a_{ij})+|\nu (a_{ij})|)(1+d)}`$
$``$ $`(1+d){\displaystyle |\lambda _i||\mu _j|\frac{L_2((a_{ij}))}{ij}}`$
$``$ $`L_2((a_{ij}))(1+d){\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n^2}}\sqrt{{\displaystyle |\lambda _i|^2}}\sqrt{{\displaystyle |\mu _i|^2}}.`$
This proves the lemma with $`C=(1+d)_{n=1}^{\mathrm{}}\frac{1}{n^2}`$.  $`\mathrm{}`$

###### Lemma 3.3.

Let $`_1=\{a𝒜_\nu |L_k(a)1,a1\}.`$ Then $`_1`$ is totally bounded in norm for $`k>2`$.

Proof: Let $`ϵ>0`$ be given. Choose $`N`$ such that $`(\frac{1}{N})^{k2}<ϵ`$. For $`G=((g_{ij}))𝒜_\nu `$ let $`P_N(G)𝒦(l^2())𝒜`$ be the element given by

$$P_N(G)_{ij}=\{\begin{array}{cc}g_{ij}\hfill & \text{ for }i,jN,\hfill \\ 0\hfill & \text{ otherwise.}\hfill \end{array}$$

Now observe that

$`L_k(GP_N(G))`$ $`=`$ $`sup_{iN\text{o}rjN}(i+j)^k(L(g_{ij})+|\nu (g_{ij})|)`$
$``$ $`N^{k2}sup_{iN\text{o}rjN}(i+j)^2(L(g_{ij})+|\nu (g_{ij})|)`$
$`=`$ $`N^{k2}L_2(GP_N(G)).`$
Note that for $`G_1,L_k(GP_N(G))1`$, therefore

$`GP_N(G)`$ $``$ $`CL_2(GP_N(G))`$
$``$ $`CN^{(k2)}L_k(GP_N(G))<Cϵ.`$
Here the constant $`C`$ is the one obtained in the previous lemma. Note $`C`$ does not depend on $`N`$. By theorem 2.8 there exists $`N\times N`$ matrices $`((a_{ij}^{(r)}))M_N(𝒜),`$ for $`r=1,\mathrm{},l`$ such that for any $`N\times N`$ matrix $`((a_{ij}))_1,`$ there exists $`r`$ satisfying $`((a_{ij}))((a_{ij}^{(r)}))<ϵ`$. Now for $`G_1`$, get $`((a_{ij}^{(r)}))`$ such that $`P_N(G)((a_{ij}^{(r)}))<ϵ`$. Then,

$$G((a_{ij}^{(r)}))GP_N(G)+ϵ(1+C)ϵ.$$

This completes the proof.  $`\mathrm{}`$

###### Theorem 3.4.

$`((𝒜_\nu ,I),L_k)`$ is a compact quantum metric space for $`k>2.`$

Proof: In view of theorem 2.8 and the previous lemma we only have to show that $`(𝒜_\nu ,L_k)`$ has finite radius. Let $`\mu _1,\mu _2S(𝒜_\nu ),a𝒜_\nu `$ with $`L_k(a)1`$. By lemma 3.2 $`aC`$, because $`L_2(a)L_k(a)`$. Hence $`|\mu _1(a)\mu _2(a)|2C`$, that is $`diam(𝒜_\nu ,L_k)2C.`$  $`\mathrm{}`$

###### Proposition 3.5.

Let

$$0A_0\stackrel{i}{}A_1\stackrel{\pi }{}A_20$$

be a short exact sequence of $`C^{}`$-algebras, with $`A_1,A_2`$ unital and a positive linear splitting $`\sigma :A_2A_1`$. Let $`\varphi :A_1^{}A_0^{}A_2^{},\psi :A_0^{}A_2^{}A_1^{}`$ be the bounded linear maps given by

$$\varphi (\mu )=(\mu _1,\mu _2),\mu _1=\mu |_{i(A_0)},\mu _2=\mu \sigma $$

$$\psi (\mu _1,\mu _2)=\mu ,\mu (a)=\mu _2(\pi (a))+\mu _1(a\sigma \pi (a))$$

Then $`\mu _1,\mu _2`$ are inverse to each other.

Proof: Let $`\varphi (\mu )=(\mu _1,\mu _2),\psi (\mu _1,\mu _2)=\mu ^{}`$. Then

$`\mu ^{}(a)`$ $`=`$ $`\mu _2(\pi (a))+\mu _1(a\sigma \pi (a))`$
$`=`$ $`\mu (\sigma \pi (a))+\mu (a\sigma \pi (a))`$
$`=`$ $`\mu (a).`$
Therefore $`\psi \varphi =Id_{A_1^{}}`$. Similarly one can show that the other composition is also identity.  $`\mathrm{}`$
Let $`𝒜,Lip(𝒜),L`$ be as above. Suppose we have a short exact sequence of $`C^{}`$-algebras

$$0𝒦𝒜\stackrel{i}{}\stackrel{~}{𝒜_1}\stackrel{\pi }{}\stackrel{~}{𝒜_2}0$$

with $`\stackrel{~}{𝒜_1},\stackrel{~}{𝒜_2}`$ unital and a positive unital linear splitting $`\sigma :\stackrel{~}{𝒜_2}\stackrel{~}{𝒜_1}`$. Let $`(𝒜_2,L_2)`$ be a compact quantum metric space with $`𝒜_2`$ a dense subspace of selfadjoint elements of $`\stackrel{~}{𝒜_2}`$. Define $`𝒜_1=i(\stackrel{~}{𝒜_\nu })\sigma (𝒜_2)`$. Then we have

###### Theorem 3.6.

In the above set up $`L_1:𝒜_1_+`$, given by

$$L_1(a)=L_2(\pi ((a))+L_k(a\sigma \pi (a))$$

is a Lip norm for $`k>2`$.

Proof: We break the proof in several steps.
Step (i) $`L_1(a)=0`$ iff $`aId_{A_1}`$: If part is obvious for the only if part note $`L_1(a)=0`$ gives $`\pi (a)=\lambda Id_{A_2}`$ for some $`\lambda `$ and $`L_0(a\lambda Id_{A_1})=0`$. Hence $`a=\lambda Id_{A_1}`$.
Step (ii) $`(A_1,L_1)`$ has finite radius: Let $`\mu ,\lambda S(𝒜_1)`$ and $`(\mu _1,\mu _2)=\varphi (\mu ),(\lambda _1,\lambda _2)=\varphi (\lambda )`$, where $`\varphi `$ is as in proposition 3.5. Then from the norm estimate of $`\varphi `$ obtained in proposition 3.5 we get $`\mu _i,\lambda _i(1+\sigma ),`$ for $`i=1,2`$ and positivity of $`\sigma `$ implies $`\mu _2=\lambda _2=1`$. Let $`x𝒜_1`$ with $`L(x)1`$, then

$`|\mu (x)\lambda (x)|`$ $`=`$ $`|\mu _2(\pi (x))+\mu _1(x\sigma \pi (x))\lambda _2(\pi (x))\lambda _1(x\sigma \pi (x))|`$
$``$ $`|\mu _2(\pi (x))\lambda _2(\pi (x))|+|\mu _1(x\sigma \pi (x))\lambda _1(x\sigma \pi (x))|`$
$``$ $`diam(𝒜_2,L_2)+2(1+\sigma )C`$
where $`C`$ is the constant obtained in lemma 3.2. This proves $`(𝒜_1,L_1)`$ has finite radius.
Step (iii) In view of theorem 2.8 it suffices to show that $`_1=\{a𝒜_1:a1,L_1(a)1\}`$ is totally bounded. Since $`(𝒜_\nu ,L_k)`$ and $`(𝒜_2,L_2)`$ are compact quantum metric spaces it follows that if we have a sequence $`a_n_1`$, then there exists a subsequence $`a_{n_k}`$ such that both $`\pi (a_{n_k})`$ and $`a_{n_k}\sigma \pi (a_{n_k})`$ converges in norm. Hence $`a_{n_k}`$ converges in norm implying the totally boundedness.  $`\mathrm{}`$

## 4. Examples

###### Example 4.1.

Let $`\mathrm{\Omega }`$ be a strongly pseudoconvex domain in $`^n`$. Let $`H^2(\mathrm{\Omega })`$ be the closure in $`L^2(\mathrm{\Omega })`$ of boundary values of holomorphic functions that can be continuously extended to $`\overline{\mathrm{\Omega }}`$. For $`fC(\mathrm{\Omega })`$ let $`T_f`$ be the associated Toeplitz operator, that is the compression of the multiplication operator $`M_f`$ on $`L^2(\mathrm{\Omega })`$ on $`H^2(\mathrm{\Omega })`$. Let $`𝔗(\mathrm{\Omega })`$ be the associated Toeplitz extension, that is the $`C^{}`$-algebra generated by the operators $`T_f`$ along with the compacts. Then we have a short exact sequence of $`C^{}`$-algebras

$$0𝒦(H^2(\mathrm{\Omega }))\stackrel{i}{}𝔗(\mathrm{\Omega })\stackrel{\pi }{}C(\mathrm{\Omega })0$$

Since this sequence admits a positive unital splitting by the previous theorem we get CQMS structure on $`𝔗(\mathrm{\Omega })`$.

###### Example 4.2.

The $`C^{}`$-algebra of continuous functions on the quantum $`SU(2)`$, to be denoted by $`C(SU_q(2))`$, is the universal $`C^{}`$-algebra generated by two elements $`\alpha `$ and $`\beta `$ satisfying the following relations:

$`\alpha ^{}\alpha +\beta ^{}\beta =I,`$ $`\alpha \alpha ^{}+q^2\beta \beta ^{}=I,`$
$`\alpha \beta q\beta \alpha =0,`$ $`\alpha \beta ^{}q\beta ^{}\alpha =0,`$
$`\beta ^{}\beta `$ $`=`$ $`\beta \beta ^{}.`$
The $`C^{}`$-algebra $`C(SU_q(2))`$ can be described more concretely as follows. Let $`\{e_i\}_{i0}`$ and $`\{e_i\}_i`$ be the canonical orthonormal bases for $`L_2(_0)`$ and $`L_2()`$ respectively. We denote by the same symbol $`N`$ the operator $`e_kke_k`$, $`k0`$, on $`L_2(_0)`$ and $`e_kke_k`$, $`k`$, on $`L_2()`$. Similarly, denote by the same symbol $`\mathrm{}`$ the operator $`e_ke_{k1}`$, $`k1`$, $`e_00`$ on $`L_2(_0)`$ and the operator $`e_ke_{k1}`$, $`k`$ on $`L_2()`$. Now take $``$ to be the Hilbert space $`L_2(_0)L_2()`$, and define $`\pi `$ to be the following representation of $`C(SU_q(2))`$ on $``$:

$$\pi (\alpha )=\mathrm{}\sqrt{Iq^{2N}}I,\pi (\beta )=q^N\mathrm{}.$$

Then $`\pi `$ is a faithful representation of $`C(SU_q(2))`$, so that one can identify $`C(SU_q(2))`$ with the $`C^{}`$-subalgebra of $`()`$ generated by $`\pi (\alpha )`$ and $`\pi (\beta )`$. Image of $`\pi `$ contains $`𝒦C(𝕋)`$ as an ideal with $`C(𝕋)`$ as the quotient algebra, that is we have a useful short exact sequence

(4.1)
$$0𝒦C(𝕋)\stackrel{i}{}𝒜\stackrel{\sigma }{}C(𝕋)0.$$

The homomorphism $`\sigma `$ is explicitly given by $`\sigma (\alpha )=\mathrm{},\sigma (\beta )=0`$. It is easy to see that the above short exact sequence admits a positive splitting taking $`z^nC(𝕋)`$ to $`\mathrm{}^nI`$, for all $`n0`$. Hence we get a compact quantum metric space structure on $`C(SU_q(2))`$.

###### Example 4.3.

Quantum sphere was introduced by Podles in . This is the universal $`C^{}`$-algebra denoted by $`C(S_{qc}^2)`$, generated by two elements $`A`$ and $`B`$ subject to the following relations:

$`A^{}=A,`$ $`B^{}B=AA^2+cI,`$
$`BA=q^2AB,`$ $`BB^{}=q^2Aq^4+cI.`$
Here the deformation parameters $`q,c`$ satisfy $`|q|<1,c>0`$. For later purpose we also note down two irreducible representations such that the representation given by the direct sum of these two is faithful. Let $`_+=l^2(_0),_{}=_+.`$ Define $`\pi _\pm (A),\pi _\pm (B):_\pm _\pm `$ by

$`\pi _\pm (A)(e_n)=\lambda _\pm q^{2n}e_n`$ $`\mathrm{where}`$ $`\lambda _\pm ={\displaystyle \frac{1}{2}}\pm (c+{\displaystyle \frac{1}{4}})^{1/2}`$
$`\pi _\pm (B)(e_n)=c_\pm (n)^{1/2}e_{n1}`$ $`\mathrm{where}`$ $`c_\pm (n)=\lambda _\pm q^{2n}(\lambda _\pm q^{2n})^2+c,\text{and}e_1=0.`$
Since $`\pi =\pi _+\pi _{}`$ is a faithful representation we have (),
(i)$`C(S_{qc}^2)C^{}(𝔗)_\sigma C^{}(𝔗):=\{(x,y):x,yC^{}(𝔗),\sigma (x)=\sigma (y)\}`$ where $`C^{}(𝔗)`$ is the Toeplitz algebra and $`\sigma :C^{}(𝔗)C(𝕋)`$ is the symbol homomorphism.
(ii) We have a short exact sequence

(4.2) $`0𝒦\stackrel{i}{}C(S_{qc}^2)\stackrel{\alpha }{}C^{}(𝔗)0`$
As in the earlier case this short exact sequence is also split exact. Here a positive splitting is given by $`\mathrm{}C^{}(𝔗)(\mathrm{},\mathrm{})`$. Now to apply the basic theorem note that by the earlier example on Toeplitz extensions we already have a Lip norm on a dense subspace of $`C^{}(𝔗)`$.


format: text

text:
# MOLECULAR SCALE CONTACT LINE HYDRODYNAMICS OF IMMISCIBLE FLOWS

## I Introduction

Immiscible two-phase flow in the vicinity of the contact line (CL), where the fluid-fluid interface intersects the solid wall, is a classical problem that falls beyond the framework of conventional hydrodynamics . In particular, molecular dynamics (MD) studies have shown relative slipping between the fluids and the wall, in violation of the no-slip boundary condition . There have been numerous ad-hoc models  to address this phenomenon, but none was able to give a quantitative account of the MD slip velocity profile in the molecular-scale vicinity of the CL. While away from the moving CL the small amount of relative slipping was found to follow the Navier boundary condition (NBC) , i.e., relative slipping proportional to the tangential viscous stress, in the molecular-scale vicinity of the CL the NBC failed totally to account for the near-complete slip. This failure casts doubts on the general applicability of the NBC to immiscible flows and hinders a continuum formulation of the hydrodynamics in the CL region. In particular, a (possible) breakdown in the hydrodynamic description for the molecular-scale CL region has been suggested . In another approach , it was shown the MD results can be reproduced by continuum finite element simulations, provided the slip profile extracted from MD is used as input. This work demonstrated the feasibility of the hybrid algorithm, but left unresolved the problem concerning the boundary condition governing the CL motion. Without a continuum hydrodynamic formulation, it becomes difficult or impossible to have realistic simulations of micro- or nanofluidics, or of immiscible flows in porous media where the relative wetting characteristics, the moving CL dissipation, and behavior over undulating solid surfaces may have macroscopic implications.

From MD simulations on immiscible two-phase flows, we report the finding that the generalized Navier boundary condition (GNBC) applies for all boundary regions, whereby the relative slipping is proportional to the sum of tangential viscous stress and the uncompensated Young stress. The latter arises from the deviation of the fluid-fluid interface from its static configuration . By combining GNBC with the Cahn-Hilliard (CH) hydrodynamic formulation of two-phase flow , we obtained a consistent, continuum description of immiscible flow with material parameters (such as viscosity, interfacial tension, etc) directly obtainable from MD simulations. The convective-diffusive dynamics in the vicinity of the interface and the moving CL also means the introduction of two phenomenological dynamic parameters whose values can be fixed by comparison with one MD flow profile. Once the parameter values are determined from MD simulations, our continuum hydrodynamics can yield predictions matching those from MD simulations (for different Couette and Poiseuille flows). Our findings suggest the no-slip boundary condition to be an approximation to the GNBC, accurate for most macroscopic flows but failing in immiscible flows. These results open the door to efficient simulations of nano- or microfluidics involving immiscible components, as well as to macroscopic immiscible flow calculations, e.g., in porous media, that are physically meaningful at the molecular level . The latter is possible, for example, by employing the adaptive method based on the iterative grid redistribution introduced in Ref. . This method has demonstrated the capability of resolving, at the same time, both the global behavior of a partial differential equation solution with coarse mesh, and a strong singularity in a localized region with a refined local mesh of over $`10^5`$ ratio to the coarse mesh.

## II Molecular dynamics simulations

The MD simulations were performed for both the static and dynamic configurations in Couette and Poiseuille flows. The two immiscible fluids were confined between two parallel walls separated along the $`z`$ direction, with the fluid-solid boundaries defined by $`z=0,H`$ (see Fig. 1 for Couette geometry). Interaction between the fluid molecules was modeled by a modified Lennard-Jones (LJ) potential $`U_{ff}=4ϵ\left[\left(\sigma /r\right)^{12}\delta _{ff}\left(\sigma /r\right)^6\right]`$, where $`r`$ is the distance between the molecules, $`ϵ`$ and $`\sigma `$ are the energy scale and range of interaction, respectively, and $`\delta _{ff}=1`$ for like molecules and $`\delta _{ff}=1`$ for molecules of different species. Each of the two walls was constructed by two (or more)  planes of an fcc lattice (see Appendix A), with each wall molecule attached to the lattice site by a harmonic spring. The mean-squared displacement of wall molecules was controlled to obey the Lindemann criterion. The wall-fluid interaction was also modeled by a LJ potential $`U_{wf}`$, with energy and range parameters $`ϵ_{wf}=1.16ϵ`$ and $`\sigma _{wf}=1.04\sigma `$, and a $`\delta _{wf}`$ for specifying the wetting property of the fluid. Both $`U_{ff}`$ and $`U_{wf}`$ were cut off at $`r_c=2.5\sigma `$. The mass of the wall molecule was set equal to that of the fluid molecule $`m`$, and the average number densities for the fluids and wall were set at $`\rho =0.81/\sigma ^3`$ and $`\rho _w=1.86/\sigma ^3`$, respectively. The temperature was controlled at $`2.8ϵ/k_B`$, where $`k_B`$ is Boltzmann’s constant. Moving the top and bottom walls at a constant speed $`V`$ in the $`\pm x`$ directions, respectively, induced the Couette flow . Applying a body force $`m𝐠_{ext}`$ to each fluid molecule in the $`x`$ direction induced the Poiseuille flow . Periodic boundary conditions were imposed on the $`x`$ and $`y`$ boundaries of the sample. Most of our MD simulations were carried out on samples consisting 6144 atoms for each fluid and 2880 atoms for each wall. The sample is $`163.5\sigma `$ by $`6.8\sigma `$ along the $`x`$ and $`y`$, respectively, and $`H=13.6\sigma `$. Our MD results represent time averages over 20 to 40 million time steps. For technical details of our MD simulations, we followed those described in Ref. .

Two different cases were considered in our simulations. The symmetric case refers to identical wall-fluid interactions for the two fluids (both $`\delta _{wf}=1`$), which leads to a flat static interface in the $`yz`$ plane with a $`90^{}`$ contact angle. The asymmetric case refers to different wall-fluid interactions, with $`\delta _{wf}=1`$ for one and $`\delta _{wf}=0.7`$ for the other. The resulting static interface is a circular arc with a $`64^{}`$ contact angle. We measured six quantities in the Couette-flow steady states of $`V=0.25(ϵ/m)^{1/2}`$, $`H=13.6\sigma `$ for the symmetric case and $`V=0.2(ϵ/m)^{1/2}`$, $`H=13.6\sigma `$ for the asymmetric case: $`v_x^{slip}`$, the slip velocity relative to the moving wall; $`G_x^w`$, the tangential force per unit area exerted by the wall; the $`\sigma _{xx}`$, $`\sigma _{nx}`$ components of the fluid stress tensor ($`n`$ denotes the outward surface normal), and $`v_x`$, $`v_z`$.

We denote the region within $`0.85\sigma =z_0`$ of the wall the boundary layer (BL). It must be thin enough to render sufficient precision for measuring $`v_x^{slip}`$, while thick enough to fully account for the tangential wall-fluid interaction force, due to the finite range of the LJ interaction. Thus it is not possible to do MD measurements strictly at the fluid-solid boundary, not only because of poor statistics, but also because of this intrinsic limitation. The wall force can be identified by separating the force on each fluid molecule into wall-fluid and fluid-fluid components. For $`0<zz_0`$ the fluid molecules can detect the atomic structure of the wall. When coupled with kinetic collisions with the wall molecules, there arises a nonzero tangential wall force that varies along the $`z`$ direction and saturates at $`zz_0`$. $`G_x^w`$ is the saturated total tangential wall force per unit wall area (Fig. 2). In Appendix A we give account of our MD results on both the tangential and normal components of the wall force, plus the effect(s) of increasing the wall thickness in our simulations from $`2`$ layers of wall molecules to $`4`$ layers and to infinite layers (by using the continuum approximation beyond the $`4`$ layers).

Spatial resolution along the $`x`$ and $`z`$ directions was achieved by evenly dividing the sampling region into bins, each $`\mathrm{\Delta }x=0.425\sigma `$ by $`\mathrm{\Delta }z=0.85\sigma `$ in size. $`v_x^{slip}`$ was obtained as the time average of fluid molecules’ velocities inside the BL, measured with respect to the moving wall; $`G_x^w`$ was obtained from the time average of the total tangential wall force experienced by the fluid molecules in the BL, divided by the bin area in the $`xy`$ plane; $`\sigma _{xx(nx)}`$ was obtained from the time averages of the kinetic momentum transfer plus the fluid-fluid interaction forces across the constant-$`x(z)`$ bin surfaces, and $`v_{x(z)}`$ was measured as the time-averaged velocity component(s) within each bin. For the contribution of intermolecular forces to the stress, we have directly measured the fluid-fluid interaction forces across bin surfaces instead of using the Irving-Kirkwood expression , whose validity was noted to be not justified at the fluid-fluid or the wall-fluid interface (see the paragraph following equation (5.15) in the above reference). In Appendix B we give some details on our MD stress measurements. As reference quantities, we also measured $`G_x^{w0}`$, $`\sigma _{xx}^0`$, $`\sigma _{nx}^0`$ in the static ($`V=0`$) configuration. In addition, we measured in both the static and dynamic configurations the average molecular densities $`\rho _1`$ and $`\rho _2`$ for the two fluid species in each bin to determine the interface profile. The shear viscosity $`\eta =1.95\sqrt{ϵm}/\sigma ^2`$ and the interfacial tension $`\gamma =5.5ϵ/\sigma ^2`$ were also determined.

We have also measured the interface and velocity profiles for the Poiseuille flow in the asymmetric case, as well as for Couette flows with different $`V`$ and $`H`$ in the symmetric case.

## III Generalized Navier boundary condition

In the presence of a fluid-fluid interface, the static fluid stress tensor $`𝝈^0`$ contains static Young stress (surface tension) as well as those stresses arising from wall-fluid interaction. For the consideration of moving CL, we will be concerned with the part of the fluid stress tensor which is purely dynamic in origin, i.e., arising purely from the hydrodynamic motion of the fluid (and the CL). In the notations below, the over tilde will denote the quantity to be the difference between that quantity and its static part. Thus if $`𝝈`$ is the total stress, we will be concerned only with the hydrodynamic part, denoted by $`\stackrel{~}{𝝈}=𝝈𝝈^0`$. We note that in the absence of body forces, force equilibrium in bulk fluid is governed by the relation $`\stackrel{~}{𝝈}=0`$. As shown in Fig. 2, this relation is altered in the BL, where the fluid-wall interaction means the existence of a dynamic, tangential force density $`\stackrel{~}{g}_x^w`$ such that $`\stackrel{~}{𝝈}+\stackrel{~}{g}_x^w\widehat{𝐱}=0`$ inside the BL. The total tangential force exerted by the wall on the fluid is given by $`\stackrel{~}{G}_x^w=_0^{z_0}𝑑z\stackrel{~}{g}_x^w`$ per unit wall area. In steady state, this wall force is necessarily balanced by the tangential fluid force $`\stackrel{~}{G}_x^f=_0^{z_0}𝑑z\left(_x\stackrel{~}{\sigma }_{xx}+_z\stackrel{~}{\sigma }_{zx}\right)`$. Here $`_{x,z,n}`$ means taking partial derivative with respect to $`x`$, $`z`$, or surface normal.

We now present evidences to show that everywhere on the boundaries, relative slipping is proportional to $`\stackrel{~}{G}_x^f`$ (the GNBC, see also equation (3) below):

$$\stackrel{~}{G}_x^f=\beta v_x^{slip},$$
(1)
where $`\beta `$ is the slip coefficient and $`\stackrel{~}{G}_x^f`$ can be written as

$$\stackrel{~}{G}_x^f=_x_0^{z_0}𝑑z\stackrel{~}{\sigma }_{xx}(z)\stackrel{~}{\sigma }_{nx}(z_0),$$
(2)
where we have used the fact that $`\stackrel{~}{\sigma }_{zx}(0)=0`$. (More strictly, $`\stackrel{~}{\sigma }_{zx}(0^{})=0`$ because there is no fluid below $`z=0`$ and hence no momentum transport across $`z=0`$) Here the $`z`$ coordinate is for the lower fluid-solid boundary (same below), with the understanding that the same physics holds at the upper boundary, and $`_n=_z`$ for the lower boundary.

We have verified the steady state force balance $`\stackrel{~}{G}_x^w+\stackrel{~}{G}_x^f=0`$ on the two boundaries (inset to Fig. 3) . This relation reflects the fact that at steady state, the frictional force exerted by the solid wall on the moving/slipping fluid is fully accounted for in $`\stackrel{~}{G}_x^f`$. Thus the GNBC (or NBC) can be expressed in either $`\stackrel{~}{G}_x^f`$ or $`\stackrel{~}{G}_x^w`$, but not both. In Fig. 3 we show the measured MD data for the symmetric and asymmetric cases in the Couette geometry. The symbols represent the values of $`\stackrel{~}{G}_x^f`$ measured in the BL. The solid lines represent the values of $`\stackrel{~}{G}_x^f`$ calculated from $`\beta v_x^{slip}`$ by using $`\beta =\beta _1=\beta _2=1.2\sqrt{ϵm}/\sigma ^3`$ for the symmetric case and $`\beta _1=1.2\sqrt{ϵm}/\sigma ^3`$, $`\beta _2=0.532\sqrt{ϵm}/\sigma ^3`$ for the asymmetric case away from the CL region (straight line segments in Fig. 3), and $`\beta =(\beta _1\rho _1+\beta _2\rho _2)/(\rho _1+\rho _2)`$ in the CL region , with $`v_x^{slip}`$ and $`\rho _{1,2}`$ obtained from MD simulations. It is seen that for the lower boundary (upper right panel), the MD data agree well with the predictions of Eq. (1). For the upper boundary (lower left panel) the straight line segments also agree well with Eq. (1). However, there is some discrepancy in the interfacial region of the upper boundary that seems to arise from a “shear thinning” effect of decreasing $`\beta `$ at very large tangential stresses .

The fact that the wall force density is distributed inside a thin BL and vanishes beyond the BL necessitates the form of $`\stackrel{~}{G}_x^f`$ as defined by Eq. (2). However, it is intuitively obvious that the fluids would experience almost the identical physical effect(s) from a wall force density $`\stackrel{~}{G}_x^w\delta (z)`$, concentrated strictly at the fluid-solid boundary with the same total wall force per unit area. In the inset to Fig. 2, it is shown that the MD-measured wall force density is a sharply peaked function. The sharp boundary limit involves the approximation of replacing this peaked function by $`\delta (z)`$. The replacement of a diffuse boundary by a sharp boundary can considerably simplify the form of the GNBC, because local force balance along $`x`$ then requires $`_x\stackrel{~}{\sigma }_{xx}+_z\stackrel{~}{\sigma }_{zx}=0`$ away from the boundary $`z=0`$. Integration of this relation from $`0^+`$ to $`z_0`$ yields $`_x_0^{z_0}𝑑z\stackrel{~}{\sigma }_{xx}(z)+\stackrel{~}{\sigma }_{zx}(z_0)\stackrel{~}{\sigma }_{zx}(0^+)=0`$ and as a consequence (by comparing with Eq. (2)) $`\stackrel{~}{G}_x^f=\stackrel{~}{\sigma }_{nx}(0^+)`$. Therefore, $`\stackrel{~}{\sigma }_{zx}`$ changes from $`\stackrel{~}{\sigma }_{zx}(0^{})=0`$ to $`\stackrel{~}{\sigma }_{zx}(0^+)=\stackrel{~}{G}_x^f`$ at $`z=0`$, leading to $`\stackrel{~}{𝝈}=\stackrel{~}{G}_x^f\delta (z)`$. Comparing with the diffuse boundary, we see that the form of the equation remains the same, but the BL is now from $`0^{}`$ to $`0^+`$. Thus the GNBC (1) becomes $`\stackrel{~}{\sigma }_{nx}(0)=\beta v_x^{slip}`$ in the sharp boundary limit.

The tangential stress $`\stackrel{~}{\sigma }_{nx}`$ can be decomposed into a viscous component and a non-viscous component: $`\stackrel{~}{\sigma }_{nx}(z)=\sigma _{nx}^v(z)+\stackrel{~}{\sigma }_{nx}^Y(z)`$. In Fig. 4 we show that away from the interfacial region the tangential viscous stress $`\sigma _{nx}^v(z)=\eta (_nv_x+_xv_n)(z)`$ is the only nonzero component, but in the interfacial region $`\stackrel{~}{\sigma }_{nx}^Y=\sigma _{nx}\sigma _{nx}^v\sigma _{nx}^0=\sigma _{nx}^Y\sigma _{nx}^0`$ is dominant, thereby accounting for the failure of NBC to describe the CL motion. Therefore away from the CL region the NBC is valid, but in the interfacial region the NBC clearly fails to describe the CL motion. We wish to clarify the origin of $`\sigma _{nx}^Y`$ and $`\sigma _{nx}^0`$ as the dynamic and static Young stresses, respectively, so that $`\stackrel{~}{\sigma }_{nx}^Y=\sigma _{nx}^Y\sigma _{nx}^0`$ is the uncompensated Young stress. As shown in the inset to Fig. 4, the integrals (across the interface) of $`\sigma _{nx}^Y`$ ($`=\sigma _{nx}\sigma _{nx}^v`$, calculated by subtracting the viscous component $`\eta (_nv_x+_xv_n)`$ from the total tangential stress $`\sigma _{nx}`$) and $`\sigma _{nx}^0`$ are equal to $`\gamma \mathrm{cos}\theta _d`$ and $`\gamma \mathrm{cos}\theta _s`$, respectively, at different values of $`z`$, i.e., $`_{\mathrm{int}}𝑑x\sigma _{nx}^Y(z)=\gamma \mathrm{cos}\theta _d(z)`$ and $`_{\mathrm{int}}𝑑x\sigma _{nx}^0(z)=\gamma \mathrm{cos}\theta _s(z)`$, where $`\theta _d(z)`$ and $`\theta _s(z)`$ are respectively the dynamic and static interfacial angles at $`z`$ . Here $`_{\mathrm{int}}𝑑x`$ denotes the integration across the fluid-fluid interface along $`x`$. These results clearly show the origin of the extra tangential stress in the interfacial region to be the interfacial (uncompensated) Young stress. Thus the GNBC is given by

$$\beta v_x^{slip}=\stackrel{~}{\sigma }_{nx}(0)=\left[\eta _nv_x\right](0)\stackrel{~}{\sigma }_{nx}^Y(0).$$
(3)
Here only one component of the viscous stress is nonzero, due to $`v_n=0`$ at the boundary; and $`\stackrel{~}{\sigma }_{nx}^Y(0)`$ is denoted the uncompensated Young stress, satisfying $`_{\mathrm{int}}\stackrel{~}{\sigma }_{nx}^Y(0)𝑑x=\gamma (\mathrm{cos}\theta _d^{surf}\mathrm{cos}\theta _s^{surf})`$, with $`\theta _{d(s)}^{surf}`$ being a microscopic dynamic (static) contact angle at the fluid-solid boundary. The fact that $`\stackrel{~}{\sigma }_{nx}^Y(0)0`$ away from the CL shows that the GNBC implies NBC for single phase flows.

Due to the diffuse nature of the BL in MD simulations, the contact angle $`\theta _{d(s)}^{surf}`$ cannot be directly measured. Nevertheless, they are obtainable through extrapolation by using the integrated interfacial curvature within the BL. That is, in the sharp boundary limit the force balance in the fluids is expressed by $`_x\stackrel{~}{\sigma }_{xx}+_n\stackrel{~}{\sigma }_{nx}=0`$. Integration in $`z`$ across the BL gives

$$_x_0^{z_0}𝑑z\stackrel{~}{\sigma }_{xx}(z)\sigma _{nx}^v(z_0)+\sigma _{nx}^v(0)\stackrel{~}{\sigma }_{nx}^Y(z_0)+\stackrel{~}{\sigma }_{nx}^Y(0)=0.$$
(4)
Integration (of Eq. (4) along $`x`$) across the fluid-fluid interface then yields

$$\mathrm{\Delta }\left[_0^{z_0}𝑑z\stackrel{~}{\sigma }_{xx}(z)\right]_{\mathrm{int}}𝑑x\sigma _{nx}^v(z_0)+_{\mathrm{int}}𝑑x\sigma _{nx}^v(0)+\gamma 𝒦_d\gamma 𝒦_s=0.$$
(5)
where $`\mathrm{\Delta }\left[_0^{z_0}𝑑z\stackrel{~}{\sigma }_{xx}(z)\right]`$ is the change of the $`z`$-integrated $`\stackrel{~}{\sigma }_{xx}`$ across the interface, $`𝒦_d`$ and $`𝒦_s`$ denote the dynamic and static $`z`$-integrated interfacial curvatures:

$$𝒦_d=\mathrm{cos}\theta _d(z_0)\mathrm{cos}\theta _d^{surf},$$

and

$$𝒦_s=\mathrm{cos}\theta _s(z_0)\mathrm{cos}\theta _s^{surf}.$$

Here $`\mathrm{\Delta }\left[_0^{z_0}𝑑z\stackrel{~}{\sigma }_{xx}(z)\right]`$, $`\sigma _{nx}^v(z_0)`$, $`\theta _d(z_0)`$, and $`\theta _s(z_0)`$ are obtainable from MD simulations, $`𝒦_s\pm 2z_0\mathrm{cos}\theta _s^{surf}/H`$ for the circular static interfaces, while $`\sigma _{nx}^v(0)=\eta [_nv_x](0)`$ may be obtained by extrapolating from the values of tangential viscous stress at $`z=z_0`$, $`2z_0`$, and $`3z_0`$. Therefore the microscopic dynamic contact angle $`\theta _d^{surf}`$ can be obtained from Eq. (5). In Appendix B 3 we give a more detailed account of the relationship between the MD measured stresses and the stress components in the continuum hydrodynamics. The above extrapolation is based on this correspondence.

We have measured the $`z`$-integrated $`\stackrel{~}{\sigma }_{xx}=\sigma _{xx}\sigma _{xx}^0`$ in the BL. The dominant behavior is a sharp drop across the interface, as shown in Fig. 5 for both the symmetric and asymmetric cases. The value of $`\theta _d^{surf}`$ obtained is $`88\pm 0.5^{}`$ for the symmetric case and $`63\pm 0.5^{}`$ for the asymmetric case at the lower boundary, and $`64.5\pm 0.5^{}`$ at the upper boundary. These values are noted to be very close to $`\theta _s^{surf}`$. Yet the small difference between the dynamic and static (microscopic) contact angles is essential in accounting for the near complete slip in the CL region.

In essence, our results show that in the vicinity of the CL, the tangential viscous stress $`\sigma _{nx}^v`$ as postulated by the NBC can not give rise to the near-complete CL slip without taking into account the tangential Young stress $`\sigma _{nx}^Y`$ in combination with the gradient of the (BL-integrated) normal stress $`\sigma _{xx}`$. For the static configuration, the normal stress gradient is balanced by the Young stress, leading to the Young’s equation. It is only for a moving CL that there is a component of the Young stress which is no longer balanced by the normal stress gradient, and this uncompensated Young stress is precisely the additional component captured by the GNBC but missed by the NBC.

## IV Continuum hydrodynamic formulation

For Eq. (3) to serve as a boundary condition in hydrodynamic calculations, we need to derive the local value of the uncompensated Young stress $`\stackrel{~}{\sigma }_{nx}^Y(0)`$ from a continuum formulation of the immiscible flow hydrodynamics. Such a formulation is important for studying the macroscopic implications of moving CL’s under scenarios beyond the capability of MD simulations. As a first-order approximation, we formulate a hydrodynamic model based on the GNBC and the CH free energy functional  that has been successful in the calculations of fluid-fluid interfacial phenomena:

$$F[\varphi ]=𝑑𝐫\left[\frac{1}{2}K\left(\varphi \right)^2+f(\varphi )\right],$$
(6)
where $`\varphi =(\rho _2\rho _1)/(\rho _2+\rho _1)`$, $`f(\varphi )=\frac{1}{2}r\varphi ^2+\frac{1}{4}u\varphi ^4`$, and $`K`$, $`r`$, $`u`$ are parameters which can be directly obtained from MD simulations through the interface profile thickness $`\xi =\sqrt{K/r}`$ , the interfacial tension $`\gamma =2\sqrt{2}r^2\xi /3u`$, and the two homogeneous equilibrium phases given by the condition of $`f/\varphi =0`$, yielding $`\varphi _\pm =\pm \sqrt{r/u}`$ ($`=\pm 1`$ in our case).

To derive the effects of the CH free energy $`F`$ on immiscible flow hydrodynamics, let us consider a composition field $`\varphi (𝐫)`$. A displacement of the molecules from $`𝐫`$ to $`𝐫^{}=𝐫+𝐮(𝐫)`$ induces a local change of $`\varphi `$, $`\delta \varphi =𝐮\varphi `$, to the first order in $`𝐮`$. The associated change in $`F`$ is given by the sum of a body term and a surface term:

$$\delta F=𝑑𝐫\left[𝐠𝐮\right]+𝑑s\left[\sigma _{ni}^Yu_i\right],$$
(7)
where $`𝐠=\mu \varphi `$ is the capillary force density, with $`\mu =\delta F/\delta \varphi =K^2\varphi r\varphi +u\varphi ^3`$ being the chemical potential, and

$$\sigma _{ni}^Y=K_n\varphi _i\varphi ,$$
(8)
is the tangential Young stress due to the spatial variation of $`\varphi `$ at the fluid-solid boundary ($`in`$). Hence the two coupled equations of motion are the Navier-Stokes equation (with the addition of the capillary force density) and the convection-diffusion equation for $`\varphi (𝐫)`$:

$$\rho _m\left[\frac{𝐯}{t}+\left(𝐯\right)𝐯\right]=p+𝝈^v+\mu \varphi +\rho _m𝐠_{ext},$$
(9)
$$\frac{\varphi }{t}+𝐯\varphi =M^2\mu ,$$
(10)
together with the incompressibility condition $`𝐯=0`$. Here $`\rho _m`$ is the fluid mass density, $`p`$ is the pressure, $`𝝈^v`$ denotes the viscous part of the stress tensor, $`\rho _m𝐠_{ext}`$ is the external body force density (for Poiseuille flows), and $`M`$ is the phenomenological mobility coefficient.

Four boundary conditions are required to solve Eqs. (9) and (10). Two are given by the impermeability condition, i.e., the normal components of the fluid velocity and diffusive flux are zero: $`v_n=0`$ and $`_n\mu =0`$. The form of the other two differential boundary conditions may be obtained from the total free energy

$$F_{tot}[\varphi ]=F[\varphi ]+𝑑s\gamma _{wf}(\varphi ),$$
(11)
plus our knowledge of GNBC. Here $`\gamma _{wf}(\varphi )`$ is the interfacial free energy per unit area at the fluid-solid boundary. We use $`\gamma _{wf}(\varphi )=(\mathrm{\Delta }\gamma _{wf}/2)\mathrm{sin}(\pi \varphi /2)`$ to denote a smooth interpolation between $`\pm \mathrm{\Delta }\gamma _{wf}/2`$, with $`\mathrm{\Delta }\gamma _{wf}=\gamma _{wf}(\varphi _+)\gamma _{wf}(\varphi _{})`$ given by $`\gamma \mathrm{cos}\theta _s^{surf}`$ (Young’s equation). It should be noted that the form of the smooth interpolation has very little effect on the final results. Hence we have chosen a simple interpolation function. Similar to Eq. (7), the change in $`F_{tot}`$ due to the displacement of the molecules from $`𝐫`$ to $`𝐫^{}=𝐫+𝐮(𝐫)`$ is given by

$$\delta F_{tot}=𝑑𝐫\left[𝐠𝐮\right]+𝑑s\left[\stackrel{~}{\sigma }_{ni}^Yu_i\right],$$
(12)
where

$$\stackrel{~}{\sigma }_{ni}^Y=\left[K_n\varphi +\frac{\gamma _{wf}(\varphi )}{\varphi }\right]_i\varphi ,$$
(13)
is the uncompensated Young stress  (see below). The continuum (differential) form of the GNBC (3) is therefore given by

$$\beta v_x^{slip}=\stackrel{~}{\sigma }_{nx}(0)=\eta \left[_nv_x\right](0)+\left[L(\varphi )_x\varphi \right](0),$$
(14)
where $`\left[L(\varphi )_x\varphi \right](0)`$, with $`L(\varphi )=K_n\varphi +\gamma _{wf}(\varphi )/\varphi `$, is the differential expression for $`\stackrel{~}{\sigma }_{nx}^Y(0)=\sigma _{nx}^Y(0)+\sigma _{nx}^0(0)`$ in equation (3). Here $`K_n\varphi _x\varphi `$ is $`\sigma _{nx}^Y(0)`$ as seen in Eq. (8), and $`[\gamma _{wf}(\varphi )/\varphi ]_x\varphi =_x\gamma _{wf}(\varphi )`$  equals to $`\sigma _{nx}^0(0)`$, in accordance with the static force balance relation $`_x\gamma _{wf}(\varphi )\sigma _{nx}^0(0)=0`$. From $`_{\mathrm{int}}𝑑x[K_n\varphi _x\varphi ](0)=\gamma \mathrm{cos}\theta _d^{surf}`$  and $`_{\mathrm{int}}𝑑x_x\gamma _{wf}=\gamma \mathrm{cos}\theta _s^{surf}`$, we see that

$$_{\mathrm{int}}𝑑x\left[L(\varphi )_x\varphi \right](0)=\gamma (\mathrm{cos}\theta _d^{surf}\mathrm{cos}\theta _s^{surf}),$$

in agreement with $`\left[L(\varphi )_x\varphi \right](0)`$ being the uncompensated Young stress.

Another boundary condition may be inferred from the fact that $`L(\varphi )=0`$ is the Euler-Lagrange equation at the fluid-solid boundary for minimizing the total free energy $`F_{tot}[\varphi ]`$. That is, $`L(\varphi )=0`$ corresponds with the equilibrium (static) condition where $`\varphi /t+𝐯\varphi =0`$. The boundary relaxation dynamics of $`\varphi `$ is plausibly assumed to be the first-order extension of that correspondence for a nonzero $`L(\varphi )`$:

$$\frac{\varphi }{t}+𝐯\varphi =\mathrm{\Gamma }\left[L(\varphi )\right],$$
(15)
where $`\mathrm{\Gamma }`$ is a (positive) phenomenological parameter.

## V Comparison of MD and continuum hydrodynamics results

Motivated by the methods presented in , a second order scheme is designed to solve the CH hydrodynamic model, comprising the dynamic equations and the four boundary conditions. Details of the numerical algorithm are presented in Appendix C. Besides those parameters which can be directly obtained from MD simulations, $`M`$ and $`\mathrm{\Gamma }`$ are treated as fitting parameters, determined by comparison with MD results (values given in the caption to Fig. 6). In Figs. 1 and 6 we show that the continuum model can indeed quantitatively reproduce the interface and velocity profiles from MD simulations, including the near-complete slip ($`v_x0`$ in Fig. 6) of the CL, the fine features in the molecular-scale vicinity of the CL, and the fast pressure variation in the BL (inset to Fig. 6), with its implied large interfacial curvature. We wish to emphasize that for the comparison with the symmetric case, the parameters in the continuum model, including those in the GNBC, are directly obtained from the MD simulations, whose velocity profiles are then fitted by those from the hydrodynamic calculations with optimized $`M`$ and $`\mathrm{\Gamma }`$ values. Thus the comparison with the asymmetric (Couette) case, with $`\beta _2`$ directly evaluated from MD simulation data, is without adjustable parameters. We have also obtained $`\theta _d^{surf}=88.1^{}`$ and $`62.8^{}`$ for the symmetric and asymmetric (the lower boundary) cases shown in Figs. 1a and 1b, respectively. Both are in excellent agreement with their extrapolated values in MD simulations. For the upper boundary in the asymmetric case, our calculated $`\theta _d^{surf}=65.2^{}`$, which differs somewhat from the MD extrapolation value of $`64.5\pm 0.5^{}`$. This difference is a reflection of the discrepancy seen in Fig. 3. However, it is noteworthy that the difference in the dynamic contact angles does not show up in the velocity profiles, which agree well.

To further verify that the boundary conditions and the parameter values are local properties and hence applicable to flows with different macroscopic conditions, we have varied the wall speed $`V`$, the system size $`H`$, and the flow geometry to check that the same set of parameters plus the GNBC are valid for reproducing (a) the velocity profiles from a different set of Couette-flow simulations in the symmetric configurations, shown in Fig. 7, as well as (b) the velocity profiles of the Poiseuille flow simulations in the asymmetric case, shown in Fig. 8. The remarkable overall agreement in all cases affirms the validity of GNBC and the hydrodynamic model , as well as justifies the replacement of the diffuse fluid-solid boundary (force density) by a sharp boundary.

Another comparison is the dissipation incurred by the moving CL in the Couette flow geometry. Rate of heat generation per unit wall area is given by $`|\stackrel{~}{G}_x^w|V=\beta |v_x^{slip}|V`$. From this we have to subtract a small but constant relative slipping away from the interface, $`v_0^{slip}=2Vl_s/(H+2l_s)`$, where $`l_s=\eta /\beta `$ is a slip length for fluid 1 if $`\varphi <0`$ and for fluid 2 if $`\varphi >0`$. The resulting heat generation rate due to the CL is $`\beta V^2W_sL`$ (for one wall), where $`L`$ is the length of the CL and $`W_s`$ defines the width of the CL region:

$$W_s=\frac{1}{V}\left(|v_x^{slip}|v_0^{slip}\right)𝑑x.$$
(16)
Thus CL dissipation is equivalent to a segment, $`H(W_s/l_s)`$, of dissipation by single phase flow. Figure 9 shows the variation of $`W_s`$ as a function of capillary number $`Ca=\eta V/\gamma `$ for the symmetric case of Couette flow. Close to $`Ca0.1`$ the value of $`W_s`$ increases rapidly, in good agreement with the MD results, and beyond which the continuum model failed to converge. This corresponds to the breakup of the interface observed in MD simulations .

## VI Concluding remarks

In summary, we have found for the first time the boundary condition that yields near-complete slipping of the CL, in good agreement with MD results on the molecular scale. It should also be noted, however, that the present continuum formulation can not calculate fluctuation effects that are important in MD simulations. Long range interactions, e.g., that due to van der Waals interaction, have also been ignored. The latter is potentially important in the calculations involving wetting layers.

## acknowledgments

Partial support from HKUST’s EHIA funding and Hong Kong RGC grants HKUST 6176/99P and 6143/01P is hereby acknowledged.

## A wall-fluid interactions

We have measured both the tangential and normal components of the wall force exerted on the fluids. Both components vary along the $`z`$ direction and saturate somewhere away from the fluid-solid boundary. The tangential component saturates (by $`99.8\%`$) at $`z=z_0`$, which is well inside the wall-fluid interaction range ($`z_0=0.85\sigma `$, smaller than the cut-off distance $`r_c=2.5\sigma `$ for the wall-fluid interaction potential $`U_{wf}`$). On the other hand, the normal component is $`87\%`$ of the saturation value at $`z=z_0`$ and $`99.8\%`$ at $`z=2z_0`$. The different saturation ranges of the tangential and normal components may be understood as follows.

For a fluid molecule close to the solid wall, the interaction with one particular (the closest) wall molecule can be much stronger than all the others. As this fluid molecule moves laterally but still remaining its close proximity to the wall, it would thus experience a strong periodic modulation in its interaction with the wall. This lateral inhomogeneity is an important source for the tangential component of the wall force. Away from the fluid-solid boundary, each fluid molecule can interact with many wall molecules on a nearly equal basis. Thus the modulation amplitude of the wall potential would clearly decrease with increasing distance from the wall. Hence the tangential wall force tends to saturate at the relatively short range of $`zz_0`$. On the contrary, the normal wall force directly arises from the wall-fluid interaction, independent of whether the wall potential is “rough” or not. Consequently, the normal wall force saturates much slower than the tangential component.

The MD results presented in this paper were obtained from simulations using solid walls constructed by two  planes of an fcc lattice. We have also carried out MD simulations using thicker confining walls. First we changed the number of molecular layers ( planes of fcc lattice) from two to four in constructing each of the two walls. The wall-fluid interaction potential $`U_{wf}`$ were still cut off at $`r_c=2.5\sigma `$. It turned out that neither component of the wall force shows any noticeable change. The reason is that for the tangential component, the two outer planes are too distant to contribute to the roughness of the wall potential, while for the normal component, the fluid molecules closest to the wall are separated from the two outer planes by a distance $`r_c`$. Consequently, both the interface and velocity profiles do not show any noticeable change.

Additional wall layers do not contribute to the perceived modulation of the wall potential by the fluid molecules. Nevertheless, they can still affect the tangential wall force by modifying the organization of the fluid molecules near the wall. Such organization is governed by the wall-fluid interaction and can be greatly influenced by the normal wall force. To see the effects of normal wall force due to additional wall layers, we used four  planes of an fcc lattice plus a half-space continuum in constructing a wall. The first four solid layers show the atomic structure detectable by the fluid molecules while the half-space continuum models the deeper solid layers. The wall-fluid interaction was modelled as follows. For an in-range pair of fluid and wall molecules separated by a distance $`r<r_c`$, the interaction potential is still $`U_{wf}`$. Here the wall molecule must be from one of the four solid layers. In addition to this short-range interaction, the fluid molecules can also experience the long-range interaction potential due to (1) the distant wall molecules in the four solid layers and (2) the continuum. For (1) we integrated the $`1/r^6`$ term in $`U_{wf}`$ over the out-of-range ($`r>r_c`$) area of the solid layers while for (2) we integrated the same term over the half-space continuum. According to this model, only the in-range ($`r<r_c`$) part of the solid wall shows atomic structure to a fluid molecule while the out-of-range ($`r>r_c`$) part is effectively a half-space continuum.

We found that the effect of the long-range normal wall force (for $`\delta _{wf}>0`$) is to attract the fluid molecules to the wall. In fact the average number density in the BL can increase by $`34\%`$ once the long-range force is included. As a result, the slip coefficient $`\beta _{1(2)}`$ increases by $`515\%`$. This results in small but visible changes in the interface and velocity profiles.

These tests have convinced us that by using two  planes of an fcc lattice to model the solid wall, we have captured the dominant wall-fluid interaction. In fact, using two molecular layers to model the solid wall has been extensively practiced in the past MD simulations , although in some instances more molecular layers have also been used , where the accurate modeling of the normal component of the wall-fluid interaction force is important.

## B stress measurements in MD simulations

### 1 Microscopic formula of Irving and Kirkwood

Irving and Kirkwood  have shown that in the hydrodynamic equation of motion (momentum transport), the stress tensor (flux of momentum) may be expressed in terms of molecular variables as

$$𝝈(𝐫,t)=𝝈_K(𝐫,t)+𝝈_U(𝐫,t),$$
(B1)
where $`𝝈_K`$ is the kinetic contribution to the stress tensor, given by

$$𝝈_K(𝐫,t)=\underset{i}{}m_i\left[\frac{𝐩_i}{m_i}𝐕(𝐫,t)\right]\left[\frac{𝐩_i}{m_i}𝐕(𝐫,t)\right]\delta (𝐱_i𝐫),f,$$
(B2)
and $`𝝈_U`$ is the contribution of intermolecular forces to the stress tensor, given by

$$𝝈_U(𝐫,t)=\frac{1}{2}\underset{i}{}\underset{ji}{}(𝐱_i𝐱_j)𝐅_{ij}\delta (𝐱_i𝐫),f.$$
(B3)
Here $`m_i`$, $`𝐩_i`$, and $`𝐱_i`$ are respectively the mass, momentum, and position of molecule $`i`$, $`𝐕(𝐫,t)`$ is the local average velocity, $`𝐅_{ij}`$ is the force on molecule $`i`$ due to molecule $`j`$, $`f`$ is the probability distribution function

$$f(𝐱_1,\mathrm{},𝐱_N,𝐩_1,\mathrm{},𝐩_N,t),$$

which satisfies the normalization condition

$$d𝐱_1\mathrm{}d𝐱_Nd𝐩_1\mathrm{}d𝐩_Nf=1,$$

and the Liouville equation

$$\frac{f}{t}=\underset{i}{}\left[\frac{𝐩_i}{m_i}\frac{f}{𝐱_i}\mathbf{}_{𝐱_i}U\frac{f}{𝐩_i}\right],$$

with $`U`$ being the potential energy of the system, and $`\mathrm{},f`$ means taking the average for a probability distribution function $`f`$.

Although widely employed in the stress measurements in MD simulations, the above expression for $`𝝈_U`$ (Eq. (B3)) represents only the leading term in an asymptotic expansion, accurate when the interaction range is small compared to the range of hydrodynamic variation . This can be seen as follows.

Consider that all the molecules interact via a pair potential $`U_{\mathrm{pair}}(R)`$ such that the intermolecular force $`𝐅_{ij}=(𝐑/R)U_{\mathrm{pair}}^{}(R)`$ for $`𝐱_j=𝐱_i+𝐑`$. Accordingly, Eq. (B3) can be rewritten as

$$𝝈_U(𝐫,t)=\frac{1}{2}𝑑𝐑\frac{\mathrm{𝐑𝐑}}{R}U_{\mathrm{pair}}^{}(R)\rho ^{(2)}(𝐫,𝐫+𝐑,t),$$
(B4)
where $`\rho ^{(2)}`$ is the pair density defined by

$$\rho ^{(2)}(𝐫_1,𝐫_2,t)=\underset{ij}{}\delta (𝐫_i𝐫_1)\delta (𝐫_j𝐫_2),f.$$

It has been shown (see the appendix in Ref. ) that according to the definition that $`d𝐒𝝈_U`$ is the force acting across $`dS`$, the full expression for $`𝝈_U`$ is given by

$$𝝈_U(𝐫,t)=\frac{1}{2}𝑑𝐑\frac{\mathrm{𝐑𝐑}}{R}U_{\mathrm{pair}}^{}(R)\left[_0^1𝑑\alpha \rho ^{(2)}(𝐫\alpha 𝐑,𝐫\alpha 𝐑+𝐑,t)\right].$$
(B5)
It is readily seen that Eq. (B4) may be obtained from Eq. (B5) by keeping only the lowest order term in a Taylor’s series in $`\alpha `$, i.e.,

$$\rho ^{(2)}(𝐫\alpha 𝐑,𝐫\alpha 𝐑+𝐑,t)\rho ^{(2)}(𝐫,𝐫+𝐑,t).$$

That means $`𝐑\mathbf{}_𝐫\rho ^{(2)}(𝐫,𝐫+𝐑,t)`$ must be negligible compared with $`\rho ^{(2)}(𝐫,𝐫+𝐑,t)`$. Here $`𝐑`$ is on the order of the range of intermolecular force. This approximation, however, can not be justified at the fluid-fluid or the wall-fluid interface, where $`𝐑\mathbf{}_𝐫\rho ^{(2)}(𝐫,𝐫+𝐑,t)`$ can be comparable in magnitude to $`\rho ^{(2)}`$.

### 2 Stress measurement in the boundary layer

In the study of moving CL, it is of great importance to obtain the correct information about stress distributions at both the fluid-fluid and the wall-fluid interfaces. Therefore, we have directly measured the $`x`$ component of fluid-fluid interaction forces acting across the $`x(z)`$ bin surfaces, in order to obtain the $`xx(zx)`$ component of $`𝝈_U`$. For example, in measuring $`\sigma _{Uzx}`$ at a given $`z`$-direction bin surface, we recorded all the pairs of molecules interacting across that surface. Here “interacting across” means that the line connecting a pair of molecules intersects the bin surface. For those pairs, we then computed $`\sigma _{Uzx}`$ at the given bin surface using

$$\sigma _{Uzx}=\frac{1}{\delta s_z}\underset{(i,j)}{}F_{ijx},$$

where $`\delta s_z`$ is the area of $`z`$-direction bin surface, $`(i,j)`$ indicate all possible pairs of molecules interacting across the bin surface, with molecule $`i`$ being “inside of $`\widehat{𝐳}\delta s_z`$” and molecule $`j`$ being “outside of $`\widehat{𝐳}\delta s_z`$” (molecule $`i`$ is below molecule $`j`$), and $`F_{ijx}`$ is the $`x`$ component of the force on molecule $`i`$ due to molecule $`j`$. A schematic illustration is shown in Fig. 10.

For comparison, we have measured the $`xx`$ and $`zx`$ components of $`𝝈_U`$ using the discrete version of Irving-Kirkwood expression (B3):

$$𝝈_U=\frac{1}{2\delta v}\underset{i}{}\underset{ji}{}(𝐱_i𝐱_j)𝐅_{ij},$$

where $`\delta v`$ is the volume of sampling bin, $`i`$ runs over fluid molecules in the sampling bin, $`j`$ runs over fluid molecules in interaction with molecule $`i`$, and $`\mathrm{}`$ means taking the time average. We found that far from the the fluid-fluid and the wall-fluid interfaces, the results based on the Irving-Kirkwood expression agree well with those from direct force measurement, whereas near the fluid-fluid or the wall-fluid interface, the two results show appreciable differences (up to $`50\%`$), especially for the $`zx`$ component at the fluid-fluid interface.

### 3 Relation of MD-measured stresses <br>to the continuum hydrodynamic stress components

We want to note the correspondence between the MD-measured stresses and the continuum hydrodynamic stress components. This correspondence is essential to obtaining the microscopic contact angle $`\theta _d^{surf}`$, defined in the continuum hydrodynamic model but not directly measurable in MD simulations.

The GNBC for the diffuse BL is given by

$$\stackrel{~}{G}_x^f=\beta v_x^{\mathrm{slip}}=\frac{}{x}_0^{z_0}𝑑z\left[\sigma _{xx}(z)\sigma _{xx}^0(z)\right]+\left[\sigma _{zx}(z_0)\sigma _{zx}^0(z_0)\right],$$
(B6)
which involves only MD measurable quantities. To obtain the contact angle $`\theta _d^{surf}`$ from MD results, we need to interpret the MD-measured quantities in terms of the various continuum variables in the hydrodynamic model. In so doing it is essential to note the following.

(1) $`\sigma _{xx}`$ can be decomposed into a molecular component and a hydrodynamic component: $`\sigma _{xx}=T_{xx}+\sigma _{xx}^{HD}`$. Meanwhile, $`\sigma _{xx}^0`$ can be composed into the same molecular component and a hydrostatic component: $`\sigma _{xx}^0=T_{xx}+\sigma _{xx}^{HS}`$. The molecular component $`T_{xx}`$ exists even if there is no hydrodynamic fluid motion or fluid-fluid interfacial curvature. In particular, $`T_{xx}`$ in the BL depends on the wall-fluid interactions. The change of the BL-integrated $`T_{xx}`$ across the fluid-fluid interface equals the change in the wall-fluid interfacial free energy, i.e., $`_{\mathrm{int}}𝑑x{\displaystyle \frac{}{x}}\left[{\displaystyle _0^{z_0}}𝑑zT_{xx}(z)\right]=\mathrm{\Delta }\gamma _{wf}=\gamma _{wf}(\varphi _+)\gamma _{wf}(\varphi _{})`$. On the other hand, the hydrodynamic component $`\sigma _{xx}^{HD}`$ in $`\sigma _{xx}`$ results from the hydrodynamic fluid motion and fluid-fluid interfacial curvature. In the static ($`V=0`$ or $`g_{ext}=0`$) configuration, $`\sigma _{xx}^{HD}`$ becomes the hydrostatic component $`\sigma _{xx}^{HS}`$ in $`\sigma _{xx}^0`$.

(2) $`\sigma _{zx}(z_0)`$ can be decomposed into a viscous component plus a Young component: $`\sigma _{zx}(z_0)=\sigma _{zx}^v(z_0)+\sigma _{zx}^Y(z_0)`$ with $`\sigma _{zx}^v=\eta (_zv_x+_xv_z)`$ and $`_{\mathrm{int}}𝑑x\sigma _{zx}^Y(z_0)=\gamma \mathrm{cos}\theta _d(z_0)`$.

(3) $`\sigma _{zx}^0(z_0)`$ is the static Young stress: i.e., $`_{\mathrm{int}}𝑑x\sigma _{zx}^0(z_0)=\gamma \mathrm{cos}\theta _s(z_0)`$.

With the help of the above relations, integration of Eq. (B6) across the fluid-fluid interface yields

$$_{\mathrm{int}}𝑑x\beta v_x^{\mathrm{slip}}=\mathrm{\Delta }\left[_0^{z_0}𝑑z\sigma _{xx}^{HD}\right]+_{\mathrm{int}}𝑑x\sigma _{zx}^v(z_0)+\gamma \mathrm{cos}\theta _d(z_0)\mathrm{\Delta }\left[_0^{z_0}𝑑z\sigma _{xx}^{HS}\right]\gamma \mathrm{cos}\theta _s(z_0),$$
(B7)
where $`\mathrm{\Delta }\left[_0^{z_0}𝑑z\sigma _{xx}^{HD(HS)}\right]`$ is the change of the $`z`$-integrated $`\sigma _{xx}^{HD(HS)}`$ across the interface. According to the Laplace’s equation, the change of the hydrostatic $`z`$-integrated normal stress is directly related to the static $`z`$-integrated curvature $`𝒦_s`$:

$$\mathrm{\Delta }_0^{z_0}𝑑z\sigma _{xx}^{HS}=\gamma 𝒦_s=\gamma \left[\mathrm{cos}\theta _s(z_0)\mathrm{cos}\theta _s^{surf}\right].$$
(B8)
Note that $`𝒦_s`$ vanishes in the symmetric case. Substituting Eq. (B8) into Eq. (B7) then yields

$$_{\mathrm{int}}𝑑x\beta v_x^{\mathrm{slip}}=\mathrm{\Delta }_0^{z_0}𝑑z\sigma _{xx}^{HD}+_{\mathrm{int}}𝑑x\sigma _{zx}^v(z_0)+\gamma \mathrm{cos}\theta _d(z_0)\gamma \mathrm{cos}\theta _s^{surf}.$$
(B9)
If interpreted in the continuum hydrodynamic formulation with a sharp fluid-solid boundary, the last term in the right-hand side of Eq. (B9), $`\gamma \mathrm{cos}\theta _s^{surf}`$, is the net wall force along $`x`$ arising from the wall-fluid interfacial free energy jump across the fluid-fluid interface, in accordance with the Young’s equation $`\gamma \mathrm{cos}\theta _s^{surf}=\mathrm{\Delta }\gamma _{wf}`$. On the other hand, the sum of the first three terms on the right-hand side of Eq. (B9) is the net fluid force along $`x`$ exerted on the three fluid sides of a BL fluid element in the interfacial region, due to the hydrodynamic motion of the fluids.

To obtain an extrapolated value for the contact angle $`\theta _d^{surf}`$ from Eq. (B9), we turn to the Stokes equation in the BL:

$$_xp+_x\sigma _{xx}^v+_z\sigma _{zx}^v+\mu _x\varphi =0,$$
(B10)
obtained from the $`x`$-component of Eq. (9) by dropping the inertial and external forces. Integration in $`z`$ of Eq. (B10) across the BL, together with the integration along $`x`$ across the fluid-fluid interface yields

$$\mathrm{\Delta }\left[_0^{z_0}𝑑z\left(p+\sigma _{xx}^v\right)\right]+_{\mathrm{int}}𝑑x\sigma _{zx}^v(z_0)+\gamma \mathrm{cos}\theta _d(z_0)_{\mathrm{int}}𝑑x\sigma _{zx}^v(0)\gamma \mathrm{cos}\theta _d^{surf}=0.$$
(B11)
Here we have made use of two relations: (1) $`\mu _x\varphi \gamma \kappa \delta (xx_{\mathrm{int}})`$ in the sharp interface limit , with $`\kappa `$ being the interfacial curvature and $`x_{\mathrm{int}}`$ the location of the interface along $`x`$. (2) $`_0^{z_0}𝑑z\kappa `$ is the dynamic $`z`$-integrated curvature $`𝒦_d=\mathrm{cos}\theta _d(z_0)\mathrm{cos}\theta _d^{surf}`$. The local force balance along $`x`$ is expressed by Eq. (B10). Accordingly, the force balance along $`x`$ for the BL fluids in the integration region is expressed by Eq. (B11), where $`\mathrm{\Delta }\left[_0^{z_0}𝑑z\left(p+\sigma _{xx}^v\right)\right]`$ is the net force on the left and right (constant-$`x`$) surfaces, $`_{\mathrm{int}}𝑑x\sigma _{zx}^v(z_0)+\gamma \mathrm{cos}\theta _d(z_0)`$ is the tangential force on the $`z=z_0`$ surface, and $`_{\mathrm{int}}𝑑x\sigma _{zx}^v(0)\gamma \mathrm{cos}\theta _d^{surf}`$ is the tangential force on the $`z=0`$ surface. Substituting Eq. (B11) into Eq. (B9) and identifying the normal stress $`p+\sigma _{xx}^v`$ with $`\sigma _{xx}^{HD}`$, we obtain

$$_{\mathrm{int}}𝑑x\beta v_x^{\mathrm{slip}}=_{\mathrm{int}}𝑑x\sigma _{zx}^v(0)+\gamma \mathrm{cos}\theta _d^{surf}\gamma \mathrm{cos}\theta _s^{surf},$$
(B12)
which is identical to the integration of the continuum GNBC (Eq. (14)) along $`x`$ across the fluid-fluid interface.

In summary, to obtain Eq. (B12) from Eq. (B7), we have used both $`_x\sigma _{xx}^{HS}+_z\sigma _{zx}^0=0`$ and $`_x\sigma _{xx}^{HD}+_z\sigma _{zx}=0`$, whose integrated expressions are given by Eq. (B8) and Eq. (B11), respectively. We note that $`_x(\sigma _{xx}^{HD}\sigma _{xx}^{HS})+_z(\sigma _{zx}\sigma _{zx}^0)=0`$ is equivalent to the relation $`_x\stackrel{~}{\sigma }_{xx}+_n\stackrel{~}{\sigma }_{nx}=0`$ (integrated expressions given by Eqs. (4) and (5)), which has been used to obtain $`\theta _d^{surf}`$ through extrapolation in Sec. III.

## C numerical algorithm

We present our numerical algorithm for solving the continuum hydrodynamic model, comprising the dynamic equations (9) and (10) and the four boundary conditions $`v_n=0`$, $`_n\mu =0`$, plus Eqs. (14) and (15). We pay special attention to the application of boundary conditions, and restrict our analysis to the Couette flow because the generalization to Poiseuille flow is straightforward.

### 1 Dimensionless hydrodynamic equations

To obtain a set of dimensionless equations suitable for numerical computations, we scale $`\varphi `$ by $`|\varphi _\pm |=\sqrt{r/u}`$, length by $`\xi =\sqrt{K/r}`$, velocity by the wall speed $`V`$, time by $`\xi /V`$, and pressure/stress by $`\eta V/\xi `$. In dimensionless forms, the convection-diffusion equation reads

$$\frac{\varphi }{t}+𝐯\varphi =_d^2(^2\varphi \varphi +\varphi ^3),$$
(C1)
the Navier-Stokes equation reads

$$\left[\frac{𝐯}{t}+\left(𝐯\right)𝐯\right]=p+^2𝐯+(^2\varphi \varphi +\varphi ^3)\varphi ,$$
(C2)
the relaxation of $`\varphi `$ at the fluid-solid boundary is governed by

$$\frac{\varphi }{t}+v_x_x\varphi =𝒱_s\left[_n\varphi \frac{\sqrt{2}}{3}\mathrm{cos}\theta _s^0s_\gamma (\varphi )\right],$$
(C3)
and the GNBC becomes

$$\left[_s(\varphi )\right]^1v_x^{slip}=\left[_n\varphi \frac{\sqrt{2}}{3}\mathrm{cos}\theta _s^0s_\gamma (\varphi )\right]_x\varphi _nv_x.$$
(C4)
Here $`\beta (\varphi )=(1\varphi )\beta _1/2+(1+\varphi )\beta _2/2`$ and $`s_\gamma (\varphi )=(\pi /2)\mathrm{cos}(\pi \varphi /2)`$. Five dimensionless parameters appear in the above equations. They are (1) $`_d=Mr/V\xi `$, which is the ratio of a diffusion length $`Mr/V`$ to $`\xi `$, (2) $`=\rho V\xi /\eta `$, (3) $`=r^2\xi /u\eta V=3\gamma /2\sqrt{2}\eta V`$, which is inversely proportional to the capillary number $`Ca=\eta V/\gamma `$, (4) $`𝒱_s=K\mathrm{\Gamma }/V`$, which is the ratio of $`K\mathrm{\Gamma }`$ (of velocity dimension) to $`V`$, and (5) $`_s(\varphi )=\eta /\beta (\varphi )\xi `$, which is the ratio of the slip length $`l_s(\varphi )=\eta /\beta (\varphi )`$ to $`\xi `$.

### 2 Finite-difference scheme

For immiscible Couette flows, there are four variables $`\varphi `$, $`v_x`$, $`v_z`$, and $`p`$ to be solved in a two-dimensional (2D) system (in the $`xz`$ plane). We want to solve the convection-diffusion equation and the Navier-Stokes equation in a 2D system of length $`L_x`$ (along $`x`$) and height $`L_z`$ (along $`z`$). Here $`L_x`$ must be large enough to allow the single phase flows (far from the fluid-fluid interface) to approach uniform shear flows. A finite-difference scheme is employed as follows. (1) $`N_x`$ and $`N_z`$ equally spaced levels are introduced in the $`x`$ and $`z`$ directions, respectively. Grid size is given by $`\mathrm{\Delta }_x=L_x/(N_x1)`$ and $`\mathrm{\Delta }_z=L_z/(N_z1)`$ along $`x`$ and $`z`$, respectively. (2) Each variable ($`q`$) is defined at $`N_x\times N_z`$ sites distributed from $`x=L_x/2`$ to $`L_x/2`$ and from $`z=L_z/2`$ to $`L_z/2`$, represented by the array $`q_{i,j}`$, with $`i=1,\mathrm{},N_x`$ and $`j=1,\mathrm{},N_z`$. Here $`q_{i,j}q(x_i,z_j)`$, with $`x_i=(i1)L_x/(N_x1)L_x/2`$ and $`z_j=(j1)L_z/(N_z1)L_z/2`$. (3) In applying the various boundary conditions, “ghost” sites outside the system, i.e., $`i=0`$, $`i=N_x+1`$, $`j=0`$, or $`j=N_z+1`$, may appear in the discretization scheme. The values of the variables at the ghost sites are determined separately from the various boundary conditions, detailed below. (4) First and second spatial derivatives along $`\zeta `$ (=$`x`$ or $`z`$) are represented by $`_\zeta q(\zeta _k)=[q(\zeta _{k+1})q(\zeta _{k1})]/2\mathrm{\Delta }_\zeta `$ and $`_\zeta ^2q(\zeta _k)=[q(\zeta _{k+1})+q(\zeta _{k1})2q(\zeta _k)]/\mathrm{\Delta }_\zeta ^2`$.

### 3 Convection-diffusion equation

With the chemical potential $`\mu _{i,j}`$ given by

$$\mu _{i,j}=\left[\frac{\varphi _{i+1,j}2\varphi _{i,j}+\varphi _{i1,j}}{\mathrm{\Delta }_x^2}+\frac{\varphi _{i,j+1}2\varphi _{i,j}+\varphi _{i,j1}}{\mathrm{\Delta }_z^2}\right]\varphi _{i,j}+\varphi _{i,j}^3,$$
(C5)
the discretized convection-diffusion equation is

$$\frac{}{t}\varphi _{i,j}+\left[𝐯\varphi \right]_{i,j}=_d\left[\frac{\mu _{i+1,j}2\mu _{i,j}+\mu _{i1,j}}{\mathrm{\Delta }_x^2}+\frac{\mu _{i,j+1}2\mu _{i,j}+\mu _{i,j1}}{\mathrm{\Delta }_z^2}\right],$$
(C6)
with

$$\left[𝐯\mathbf{}\varphi \right]_{i,j}=v_{xi,j}\frac{\varphi _{i+1,j}\varphi _{i1,j}}{2\mathrm{\Delta }_x}+v_{zi,j}\frac{\varphi _{i,j+1}\varphi _{i,j1}}{2\mathrm{\Delta }_z}.$$
(C7)
The boundary conditions at $`x=\pm L_x/2`$ can be easily applied using $`\varphi =\pm 1`$ and

$$𝐯(z)=\frac{L_z}{L_z+2L_s}\frac{2z}{L_z}\widehat{𝐱},$$
(C8)
for single-phase uniform shear flows. Here we focus on the boundary conditions at $`z=\pm L_z/2`$: $`_n\mu =0`$ and Eq. (C3). We spell out the numerics for the lower boundary $`j=1`$, with the understanding that the same can be applied to the upper boundary.

To solve the discretized convection-diffusion equation (C6) at the lower boundary $`j=1`$, we need the values of $`\mu _{i,j}`$ at $`j=1`$ and $`j=0`$. We also need the values of $`\mu _{i,j}`$ at $`j=1`$ to solve the same equation at $`j=2`$. According to Eq. (C5), $`\mu _{i,j}`$ at $`j=1`$ and $`j=0`$ can not be directly evaluated from $`\varphi _{i,j}`$ with $`i=1,\mathrm{},N_x`$ and $`j=1,\mathrm{},N_z`$. But they can still be determined from the boundary conditions at $`z=L_z/2`$. $`\mu _{i,j}`$ at $`j=0`$ is obtained from $`_n\mu =0`$ at $`j=1`$ as

$$\mu _{i,j1=0}=\mu _{i,j+1=2}.$$
(C9)
To obtain $`\mu _{i,j}`$ at $`j=1`$, we need to determine $`\varphi _{i,j}`$ at $`j=0`$. This can be done by requiring that Eqs. (C1) and (C3) yield the same $`\varphi /t`$ at $`z=L_z/2`$. The discretized convection-diffusion equation is given by Eq. (C6) while the discretized relaxation equation for $`\varphi `$ at the boundary $`j=1`$ is given by

$$\frac{}{t}\varphi _{i,j}+\left[𝐯\varphi \right]_{i,j}=𝒱_s\left[\frac{\varphi _{i,j1}\varphi _{i,j+1}}{2\mathrm{\Delta }_z}\frac{\sqrt{2}}{3}\mathrm{cos}\theta _s^0s_\gamma (\varphi _{i,j})\right].$$
(C10)
Equating the right-hand side of Eq. (C6) at $`j=1`$ (with $`\mu _{i,0}`$ fixed by Eq. (C9) and other $`\mu `$’s given by Eq. (C5)) with that of Eq. (C10) leads to a tridiagonal system of linear equations for $`\varphi _{i,j}`$ ($`\varphi _{i,j}`$ coupled with $`\varphi _{i1,j}`$ and $`\varphi _{i+1,j}`$) at $`j=0`$. Solving this tridiagonal system determines $`\varphi _{i,j}`$ at $`j=0`$, from which we obtain $`\mu _{i,j}`$ at $`j=1`$ by using Eq. (C5).

### 4 Navier-Stokes equation

We now turn to the Navier-Stokes equation (C2) with the incompressibility condition $`𝐯=0`$. The difficulty in solving the Navier-Stokes equation is the lack of a time evolution equation for the pressure $`p`$. In the following, we will introduce a numerical method based on the pressure Poisson equation .

#### a Pressure Poisson equation

Taking the divergence of momentum equation (C2) and applying the incompressibility condition, we obtain the pressure Poisson equation

$$^2p=\left[\left(𝐯\right)𝐯\right]+[(^2\varphi \varphi +\varphi ^3)\varphi ].$$
(C11)
Dotting the momentum equation (C2) with the surface normal at the fluid-solid boundary and using $`v_n=0`$, we obtain for Eq. (C11) the boundary condition

$$_np=^2v_n+(^2\varphi \varphi +\varphi ^3)_n\varphi ,$$
(C12)
at $`z=\pm L_z/2`$. In addition, we use $`^2p=0`$ and $`_xp=0`$ for the values of $`^2p`$ and $`_np`$ at the boundaries $`x=\pm L_x/2`$. This reflects the single phase flow given by Eq. (C8).

From the momentum equation (C2) and the pressure Poisson equation (C11), we derive a diffusion equation

$$\frac{(𝐯)}{t}=^2(𝐯),$$

for $`𝐯`$. With $`𝐯=0`$ given at time $`t=0`$, and in order to ensure that $`𝐯`$ remains divergence-free at $`t>0`$, we must impose the additional boundary condition $`𝐯=0`$ at all times $`t0`$. We will show that this boundary condition is needed in solving for $`p`$ in a finite-difference scheme.

In order to solve the pressure Poisson equation, we need to evaluate $`[^2p]_{i,j}`$ for $`i=1,\mathrm{},N_x`$ and $`j=1,\mathrm{},N_z`$, $`[_xp]_{i,j}`$ for $`i=1,N_x`$ and $`j=1,\mathrm{},N_z`$, and $`[_zp]_{i,j}`$ for $`i=1,\mathrm{},N_x`$ and $`j=1,N_z`$. For $`^2p`$, we have

$$[^2p]_{i,j}=0,$$

for $`i=1,N_x`$ and $`j=1,\mathrm{},N_z`$;

$$\begin{array}{cc}[^2p]_{i,j}=\hfill & 2\left[\frac{v_{xi+1,j}v_{xi1,j}}{2\mathrm{\Delta }_x}\frac{v_{zi,j+1}v_{zi,j1}}{2\mathrm{\Delta }_z}\frac{v_{zi+1,j}v_{zi1,j}}{2\mathrm{\Delta }_x}\frac{v_{xi,j+1}v_{xi,j1}}{2\mathrm{\Delta }_z}\right]\hfill \\ & +\mu _{i,j}\left(\frac{\varphi _{i+1,j}2\varphi _{i,j}+\varphi _{i1,j}}{\mathrm{\Delta }_x^2}+\frac{\varphi _{i,j+1}2\varphi _{i,j}+\varphi _{i,j1}}{\mathrm{\Delta }_z^2}\right)\hfill \\ & +\left(\frac{\mu _{i+1,j}\mu _{i1,j}}{2\mathrm{\Delta }_x}\frac{\varphi _{i+1,j}\varphi _{i1,j}}{2\mathrm{\Delta }_x}+\frac{\mu _{i,j+1}\mu _{i,j1}}{2\mathrm{\Delta }_z}\frac{\varphi _{i,j+1}\varphi _{i,j1}}{2\mathrm{\Delta }_z}\right),\hfill \end{array}$$

for $`i=2,\mathrm{},N_x1`$ and $`j=2,\mathrm{},N_z1`$; and

$$\begin{array}{cc}[^2p]_{i,j}=\hfill & 2\frac{v_{xi+1,j}v_{xi1,j}}{2\mathrm{\Delta }_x}\frac{v_{zi,j+1}v_{zi,j1}}{2\mathrm{\Delta }_z}\hfill \\ & +\mu _{i,j}\left(\frac{\varphi _{i+1,j}2\varphi _{i,j}+\varphi _{i1,j}}{\mathrm{\Delta }_x^2}+\frac{\varphi _{i,j+1}2\varphi _{i,j}+\varphi _{i,j1}}{\mathrm{\Delta }_z^2}\right)\hfill \\ & +\frac{\mu _{i+1,j}\mu _{i1,j}}{2\mathrm{\Delta }_x}\frac{\varphi _{i+1,j}\varphi _{i1,j}}{2\mathrm{\Delta }_x},\hfill \end{array}$$

for $`i=2,\mathrm{},N_x1`$ and $`j=1,N_z`$ (where $`v_z=0`$ and $`_z\mu =0`$). We see that $`\varphi `$ and $`v_z`$ at ghost sites of $`j=0,N_z+1`$ appear in the last expression. The ghost $`\varphi `$’s have already been determined in solving the convection-diffusion equation, while the ghost $`v_z`$’s are determined through the additional boundary condition $`𝐯=0`$:

$$\frac{v_{xi+1,j}v_{xi1,j}}{2\mathrm{\Delta }_x}+\frac{v_{zi,j+1}v_{zi,j1}}{2\mathrm{\Delta }_z}=0,$$

for $`i=2,\mathrm{},N_x1`$, and $`j=1,N_z`$. For $`_np`$, we have

$$[_xp]_{i,j}=0$$

for $`i=1,N_x`$ and $`j=1,\mathrm{},N_z`$;

$$[_zp]_{i,j}=0$$

for $`i=1,N_x`$ and $`j=1,N_z`$; and

$$[_zp]_{i,j}=\frac{v_{zi,j+1}+v_{zi,j1}}{\mathrm{\Delta }_z^2}+\mu _{i,j}\frac{\varphi _{i,j+1}\varphi _{i,j1}}{2\mathrm{\Delta }_z}$$

for $`i=2,\mathrm{},N_x1`$ and $`j=1,N_z`$ (where $`v_z=0`$). The last expression involves the ghost $`\varphi `$’s and $`v_z`$’s at $`j=0,N_z+1`$. Given the above values of $`[^2p]_{i,j}`$ and $`[_np]_{i,j}`$, we apply a 2D Fast Fourier Transformation to solve $`p_{i,j}(0)`$ (up to a constant) for $`i=1,\mathrm{},N_x`$ and $`j=1,\mathrm{},N_z`$.

#### b Slip boundary condition

The discretized Navier-Stokes equation is given by

$$\begin{array}{cc}\frac{v_{xi,j}}{t}=\hfill & v_{xi,j}\frac{v_{xi+1,j}v_{xi1,j}}{2\mathrm{\Delta }_x}v_{zi,j}\frac{v_{xi,j+1}v_{xi,j1}}{2\mathrm{\Delta }_z}\frac{1}{}\frac{p_{i+1,j}p_{i1,j}}{2\mathrm{\Delta }_x}\hfill \\ & +\frac{1}{}\left(\frac{v_{xi+1,j}2v_{xi,j}+v_{xi1,j}}{\mathrm{\Delta }_x^2}+\frac{v_{xi,j+1}2v_{xi,j}+v_{xi,j1}}{\mathrm{\Delta }_z^2}\right)\hfill \\ & +\frac{}{}\mu _{i,j}\frac{\varphi _{i+1,j}\varphi _{i1,j}}{2\mathrm{\Delta }_x},\hfill \end{array}$$
(C13)
for $`i=2,\mathrm{},N_x1`$ and $`j=1,\mathrm{},N_z`$, and

$$\begin{array}{cc}\frac{v_{zi,j}}{t}=\hfill & v_{xi,j}\frac{v_{zi+1,j}v_{zi1,j}}{2\mathrm{\Delta }_x}v_{zi,j}\frac{v_{zi,j+1}v_{zi,j1}}{2\mathrm{\Delta }_z}\frac{1}{}\frac{p_{i,j+1}p_{i,j1}}{2\mathrm{\Delta }_z}\hfill \\ & +\frac{1}{}\left(\frac{v_{zi+1,j}2v_{zi,j}+v_{zi1,j}}{\mathrm{\Delta }_x^2}+\frac{v_{zi,j+1}2v_{zi,j}+v_{zi,j1}}{\mathrm{\Delta }_z^2}\right)\hfill \\ & +\frac{}{}\mu _{i,j}\frac{\varphi _{i,j+1}\varphi _{i,j1}}{2\mathrm{\Delta }_z}.\hfill \end{array}$$
(C14)
for $`i=2,\mathrm{},N_x1`$ and $`j=2,\mathrm{},N_z1`$, together with the boundary conditions that $`v_{zi,j}=0`$ at $`j=1,N_z`$ and $`𝐯`$ given by Eq. (C8) at $`i=1,N_x`$. Equation (C13) at $`j=1,N_z`$ involves $`\varphi `$ and $`v_x`$ at ghost sites of $`j=0,N_z+1`$. The ghost $`\varphi `$’s come from $`\mu _{i,j}`$ at $`j=1,N_z`$, and have already been determined. The ghost $`v_x`$’s are determined from the discretized GNBC

$$\left[_s(\varphi _{i,j})\right]^1v_{xi,j}^{slip}=\left[\frac{\varphi _{i,j1}\varphi _{i,j+1}}{2\mathrm{\Delta }_z}\frac{\sqrt{2}}{3}\mathrm{cos}\theta _s^0s_\gamma (\varphi _{i,j})\right]\frac{\varphi _{i+1,j}\varphi _{i1,j}}{2\mathrm{\Delta }_x}\frac{v_{xi,j1}v_{xi,j+1}}{2\mathrm{\Delta }_z},$$
(C15)
at the lower boundary $`j=1`$ with $`v_{xi,j}^{slip}=v_{xi,j}V`$, and

$$\left[_s(\varphi _{i,j})\right]^1v_{xi,j}^{slip}=\left[\frac{\varphi _{i,j+1}\varphi _{i,j1}}{2\mathrm{\Delta }_z}\frac{\sqrt{2}}{3}\mathrm{cos}\theta _s^0s_\gamma (\varphi _{i,j})\right]\frac{\varphi _{i+1,j}\varphi _{i1,j}}{2\mathrm{\Delta }_x}\frac{v_{xi,j+1}v_{xi,j1}}{2\mathrm{\Delta }_z},$$
(C16)
at the upper boundary $`j=N_z`$ with $`v_{xi,j}^{slip}=v_{xi,j}+V`$.

In summary, to solve the dynamic equations (9) and (10), we need to use $`\varphi =\pm 1`$ and Eq. (C8) at $`x=\pm L_x/2`$, with $`v_n=0`$, $`_n\mu =0`$, plus Eqs. (14) and (15) at $`z=\pm L_z/2`$. In particular, in applying the boundary conditions at $`z=\pm L_z/2`$, values of $`\varphi `$, $`v_x`$, and $`v_z`$ at ghost sites have to be introduced and solved for.

### 5 Time integration

We outline the scheme for time discretization and integration. For simplicity we only describe the forward Euler time stepping. In the following a superscript $`n`$ denotes consecutive time instants and $`\mathrm{\Delta }t`$ is the time interval.

Time Stepping: Given $`\left\{\varphi _{i,j}^n\right\}`$ and $`\left\{𝐯_{i,j}^n\right\}`$ at all the sites ($`i=1,\mathrm{},N_x`$ and $`j=1,\mathrm{},N_z`$) in the system:

Step 1: Determine $`\left\{\mu _{i,j}^n\right\}`$, $`\left\{\varphi _{i,j}^n\right\}`$, and $`\left\{𝐯_{i,j}^n\right\}`$ at the ghost sites from the various boundary conditions, as described in Secs. C 3, C 4 a, and C 4 b.

Step 2: Solve $`\left\{p_{i,j}^n\right\}`$ at all the interior sites ($`i=1,\mathrm{},N_x`$ and $`j=1,\mathrm{},N_z`$) from Eq. (C11) with appropriate boundary conditions for $`_np`$, as described in Sec. C 4 a.

Step 3: Compute $`\left\{\varphi _{i,j}^{n+1}\right\}`$ and $`\left\{𝐯_{i,j}^{n+1}\right\}`$ at all the interior sites (except those fixed by the boundary conditions at all times) using

$$\frac{\varphi ^{n+1}\varphi ^n}{\mathrm{\Delta }t}=𝐯^n\varphi ^n+_d^2\mu ^n,$$

and

$$\frac{𝐯^{n+1}𝐯^n}{\mathrm{\Delta }t}=\left(𝐯^n\right)𝐯^np^n+^2𝐯^n+\mu ^n\varphi ^n,$$

according to Eqs. (C6), (C13), and (C14) in discretized time. Here the ghost $`\left\{\mu _{i,j}^n\right\}`$, $`\left\{\varphi _{i,j}^n\right\}`$, and $`\left\{𝐯_{i,j}^n\right\}`$ determined in Step 1 and $`\left\{p_{i,j}^n\right\}`$ solved in Step 2 are needed.