# Effective diffusivity of Brownian particles in a two dimensional square
lattice of hard disks

M. Mangeat  Univ. Bordeaux, CNRS, Laboratoire Ondes et Matière d’Aquitaine
(LOMA), UMR 5798, F-33405 Talence, France  Center for Biophysics & Department
for Theoretical Physics, Saarland University, D-66123 Saarbrücken, Germany  T.
Guérin  Univ. Bordeaux, CNRS, Laboratoire Ondes et Matière d’Aquitaine (LOMA),
UMR 5798, F-33405 Talence, France  D. S. Dean  Univ. Bordeaux, CNRS,
Laboratoire Ondes et Matière d’Aquitaine (LOMA), UMR 5798, F-33405 Talence,
France  author to whom correspondence should be addressed:
<EMAIL_ADDRESS>
###### Abstract

We revisit the classic problem of the effective diffusion constant of a
Brownian particle in a square lattice of reflecting impenetrable hard disks.
This diffusion constant is also related to the effective conductivity of non-
conducting and infinitely conductive disks in the same geometry. We show how a
recently derived Green’s function for the periodic lattice can be exploited to
derive a series expansion of the diffusion constant in terms of the disk’s
volume fraction $\varphi$. Secondly we propose a variant of the Fick-Jacobs
approximation to study the large volume fraction limit. This combination of
analytical results is shown to describe the behavior of the diffusion constant
for all volume fractions.

###### pacs:

47.57.J-, 66.10.C-, 87.15.Vv, 87.16.dp

##  I Introduction

Determining the transport properties of tracer particles in complex media such
as colloidal crystals and suspensions, porous media or living cells is an
important and complex physical problem Bressloff and Newby (2013); Burada et
al. (2009, 2008); Höfling and Franosch (2013); Dentz et al. (2011); Dean et
al. (2007). In the case of particles diffusing in crowded or confined
environments, the presence of obstacles (which can be seen as entropic traps
or barriers) hinders their motion, leading to a decreased effective
diffusivity Malgaretti et al. (2013); Burada et al. (2009, 2008); Jacobs
(1967). The relation between this diffusivity and the geometry of the
confinement is known to be non-trivial and has interested physicists over the
last few decades Yang et al. (2017); Burada et al. (2008); Reguera et al.
(2006); Kalinay and Percus (2006); Malgaretti et al. (2013, 2016).

Figure 1: Equivalence of the dispersion of Brownian particles inside an array
of obstacles (a) and the dispersion in a periodic channel (b). From the
symmetries of the periodic array of disks, the dashed lines represented in (a)
can be replaced by ghost reflecting boundaries which correspond to the
boundaries of the channel. The same numerical realization of Brownian motion
is shown in both domains.

A standard, and extensively studied model of dispersion in a complex crowded
system is that of a Brownian tracer particle which diffuses in an environment
of fixed hard spheres (Fig. 1). The late time diffusivity $D_{e}$ of a
Brownian tracer system is related to the effective conductivity of a
heterogenous medium of fixed non-conducting spheres embedded in a conducting
medium (we recall this result below), and results for this problem can be
traced back to Maxwell Maxwell (1881) and Rayleigh Rayleigh (1892). For
example, in the dilute limit where the volume fraction $\varphi$ of obstacles
tends to zero, the calculation of $D_{e}$ can be formulated as a Laplace
equation in infinite space, leading to the well known formula Maxwell (1881);
Jeffrey (1973); Batchelor (1976)

$\frac{D_{e}}{D_{0}}=1-\frac{\varphi}{d-1}+{\cal O}(\varphi^{2})$ (1)

where $d$ is the spatial dimension, $D_{0}$ the microscopic (local)
diffusivity and $\varphi$ the volume fraction of obstacles. For higher values
of $\varphi$, the above formula is not accurate, and a partially renormalized
formula has been proposed by Kalnin et al. Kalnin et al. (2001) who obtained

$\frac{D_{e}}{D_{0}}\simeq\frac{d-1}{d-1+\varphi}.$ (2)

The above expression turns out to be a much better approximation for the
effective diffusivity in the dilute regime than the leading order expression
(1). This resummation is equivalent to that proposed by Maxwell for the
conductivity problem. However, its mathematical status is rather unclear and
only the term first order term in $\varphi$ is strictly exact. It is obtained
via an effective medium approach, assuming spherical symmetry around a single
obstacle, with an effective diffusivity at infinity determined self-
consistently and representing the effect of other obstacles.

In this paper, we show that, for regular square arrays of obstacles, the
renormalized expression (2) is exact up to (and including) the order
$\varphi^{4}$. This explains why the renormalized formula is accurate up to
relatively large volume fractions, around $\varphi\simeq 0.5$ in 2D, for such
volume fractions it is rather surprising that a result obtained in the dilute
limit should be relevant. In fact, our approach uses strong localized
perturbation theory Ward and Keller (1993) and enables us to obtain exact
analytical expressions for the coefficients of $\varphi$ at all orders. We
also show that the asymptotic expansion of $D_{e}(\varphi)$ can be extracted
from the literature as the result of the multipole expansion method, it agrees
with our results (derived with a new point of view). These considerations are
the first contribution of the present paper.

Next, in 2D, we will also consider the opposite limit where $R$ approaches the
value $R_{c}$ (half the square lattice period) at which the obstacles begin to
touch each other, leaving only small openings for the tracer particle to pass
to the next pore. In this limit, the diffusivity vanishes as
$\sqrt{|R_{c}-R|}$. In fact, one can see the 2D problem as equivalent to the
problem of a tracer particle in a channel Dagdug et al. (2012) (figure 1b), by
adding a ghost reflecting boundary condition identical to the zero flux
condition of the boundary of the channel. For channel geometries, we can apply
well known Fick-Jacobs approximation, which maps the 2D problem onto the 1D
problem of a diffusing particle in an effective, entropic, potential. It is
knownDagdug et al. (2012) that this approximation can be used to estimate the
effective dispersion in the limit $R\to R_{c}$, and even improved by
introducing effective local diffusion constant $D(x)$ for the effective one
dimensional dynamicsBerezhkovskii and Szabo (2011). The second contribution of
this paper is to propose an alternative expression for the local effective
diffusivity $D(x)$ in the effective 1D description:

$\displaystyle D(x)\simeq
D_{0}\\{1-[\hat{\mathbf{e}}_{x}\cdot\mathbf{n}(x)]^{2}/3\\},$ (3)

where $\mathbf{n}$ is the normal vector to the boundary and
$\hat{\mathbf{e}}_{x}$ the unit vector in the direction parallel to the
channel. The above formula can be seen as a natural renormalization of the
effective diffusivity tensor when studying diffusion between two hypersurfaces
of arbitrary dimensions, and is exact at the same order as other expressions
of $D(x)$ proposed in the literature Zwanzig (1992); Reguera and Rubi (2001);
Kalinay and Percus (2006). In the limit of nearly touching obstacles, we find
that the use of the local diffusivity (3) leads to an closed form analytical
estimate [Eq. (50)] of the effective diffusivity which turns out to be more
accurate than with the use of other choices of $D(x)$ (and which more over do
not have simple closed form expressions).

Taken together, the analytical formulas presented in this work, in the dilute
or crowded regime, provide an accurate description of $D_{e}$ in 2D for almost
all values of the volume fraction. The outline of this paper is as follows: in
Section II we present a general formalism to investigate dispersion in
heterogeneous media. Then, we investigate the dilute limit (Section III)
employing an exact perturbative approach. Finally we focus on the limit of
nearly touching obstacles (Section IV) within the context of the improved
Fick-Jacobs approximation.

##  II Dispersion in a regular array of spherical obstacles: general formalism

We consider a periodic array of hard spheres with each sphere is taken to be
of radius $R$ and at the centre of an (hyper)-cube of length $L$. The spheres
are impenetrable and we denote by

$\displaystyle\Omega=L^{d}(1-\varphi)$ (4)

the volume accessible to the tracer in each hypercube. The overall system is
constructed by repeating the underlying cell structure periodically over all
space. We consider the motion of small tracer particles inside this array,
with a microscopic diffusivity $D_{0}$ (see figure 1a) and we define the late
time effective diffusivity, say in the $x$ direction, via

$\displaystyle
D_{e}=D_{xx}=\lim_{t\to\infty}\frac{\left\langle[x(t)-x(0)]^{2}\right\rangle}{2t},$
(5)

where $\langle...\rangle$ denotes ensemble average. It is known that this
diffusivity can be obtained by solving a partial differential equation problem
for an auxiliary function $f$, see for example the macro-transport theory of
Brenner and Edwards Brenner and Edwards (1993), or Refs. Guérin and Dean
(2015a); Guérin and Dean (2015b); Mangeat et al. (2017a) for the case of
arbitrary (possibly out-of-equilibrium) heterogeneous media. Here we use the
notations of Ref. Mangeat et al. (2017a) and express the diffusivity as

$\frac{D_{e}}{D_{0}}=1-\frac{1}{\Omega}\int_{\partial\Omega}dS\
\mathbf{n}\cdot\mathbf{e}_{x}f({\bf x}),$ (6)

where $\mathbf{n}$ is the unit normal vector, oriented towards the interior of
the obstacles, $\mathbf{e}_{x}$ is the unit vector in the $x$ direction and
the integration is taken over the surface of a single obstacle. The auxiliary
function $f$ is periodic in both $x$ and $y$, is harmonic in the available
volume

$\displaystyle\nabla^{2}f=0$ (7)

and satisfies the boundary condition

$\displaystyle{\bf n}\cdot\nabla f({\bf x})={\bf n}\cdot\mathbf{e}_{x}$
$\displaystyle(\mathbf{x}\in\partial\Omega).$ (8)

where $\partial\Omega$ denotes the boundary of the obstacles. Note that $f$ is
defined up to an unimportant additive constant. The above set of equations are
not in general analytically tractable, in particular the periodic boundary
conditions for $f$ on $\Omega$ present the main obstacle to analytical
progress. The above set of equations can however be easily solved numerically
using standard partial differential equation solvers, eliminating the need for
stochastic simulations of the diffusion process itself.

Here, it is useful to make an explicit connection with an electrostatic
problem. If we set $V=f-x$, we see that $V$ is harmonic and satisfies Neumann
boundary conditions $\partial_{n}V=0$ at the surface of obstacles.
Furthermore, the periodicity of $f$ implies that $V(x+L)-V(x)=-L$. This means
that $V$ is the electrostatic potential in a medium of non-conducting spheres
embedded in a medium of conductivity $\sigma_{0}=1$ submitted to a constant
average electric field $\overline{E}=1$, where the over-bar denotes uniform
spatial average. The effective conductivity is defined by the relation
$\sigma_{e}=\overline{\sigma E}/\overline{E}$, and from this we find

$\displaystyle\frac{\sigma_{e}}{\sigma_{0}}=-\frac{\int_{\Omega}d\mathbf{r}\
\partial_{x}V}{L^{d}\overline{E}}=\frac{1}{L^{d}}\left(\Omega-\int_{\Omega}d\mathbf{r}\
\partial_{x}f\right).$ (9)

Using the divergence theorem and Eq. (6) we obtain

$\displaystyle\frac{\sigma_{e}}{\sigma_{0}}=\frac{\Omega}{L^{d}}\frac{D_{e}}{D_{0}}=(1-\varphi)\frac{D_{e}}{D_{0}}.$
(10)

This means that, up to the multiplicative factor $(1-\varphi)$, the effective
diffusivity of a periodic medium with reflecting obstacles is exactly the same
as the effective conductivity of a medium containing non-conducting elements.
This connection is well known Dean et al. (2007, 2004); Batchelor (1974);
Kalnin et al. (2001) but the above formula is useful to extract formulas on
diffusivity from the literature on electrostatics in heterogeneous media.

From now on, we focus on the case of two dimensions.

##  III Exact asymptotic expansions of the effective diffusivity in the dilute
limit $\varphi\to 0$ in 2D

We now consider dispersion in the dilute limit, $\varphi\to 0$. In this
Section, we choose units of length such that $L=1$. As mentioned above, the
main obstacle to analytical progress is the periodicity condition on $f$.
However, in the dilute limit, it is convenient to look at the behavior of $f$
at the vicinity of the obstacles and to solve the equations by requiring that
$f$ vanishes at infinity (meaning that that periodic conditions are
approximately satisfied). The solution then reads

$\displaystyle f\simeq-\cos\theta R^{2}/r,$ (11)

where $R$ is the radius of the obstacles, and $(r,\theta)$ are polar
coordinates when the origin is set at the center of the obstacles. The above
expression is the solution of Eqs. (7,8) in infinite space, and inserting it
into (6) yields

$\displaystyle D_{e}=1+\frac{R}{\Omega}\int_{0}^{\pi}d\theta\cos\theta
f(R,\theta)=1-R^{2}\int_{0}^{\pi}d\theta\cos^{2}\theta=1-\varphi+o(\varphi),$
(12)

where $\varphi=\pi R^{2}/L^{2}$ is the volume fraction. This is the Maxwell
result (1) at first order. To go to next order, we consider the problem as a
strong localized perturbation problem (as described in Ref. Ward and Keller
(1993)), in the spirit of boundary layer theory. The function $f$ is assumed
to have components at scale $R$, such that its behavior when $R\to 0$ at fixed
$r/R$ is given by

$\displaystyle
f(\mathbf{x})=R\left[\Phi_{0}(r/R,\theta)+R^{2}\Phi_{1}(r/R,\theta)+R^{4}\Phi_{2}(r/R,\theta)+...\right],$
(13)

with $\theta$ the angle with respect to $\mathbf{e}_{x}$, the usual convention
for polar coordinates. The first term of this expansion comes from the
comparison with Eq. (11), which enables us to identify $\Phi_{0}$ as

$\displaystyle\Phi_{0}(\tilde{r},\theta)=-\cos\theta/\tilde{r},$ (14)

where $\tilde{r}=r/R$ is the renormalized radius. The magnitude of the next
order terms in the expansion (13) at this stage is not justified, but it will
become clear that this expansion is in fact correct. For the outer solution,
obtained by considering the limit $R\to 0$ when $r$ is fixed, we assume the
following ansatz

$\displaystyle\Phi(\mathbf{x})=R^{2}[f_{1}(x,y)+R^{2}f_{2}(x,y)+R^{4}f_{3}(x,y)+...],$
(15)

where one has anticipated forthcoming calculations by assuming that the
underlying expansion parameter is $R^{2}$. The leading order term is however
determined from the fact that the outer solution and the inner solution must
match in the regime $R\ll r\ll 1$, so that the function $f_{1}$ must behave as

$\displaystyle f_{1}(r,\theta)\underset{r\to 0}{\sim}-\cos\theta/r,$ (16)

where the right hand side is identified by considering Eq. (11). The function
$f_{1}$ must be harmonic and periodic in both $x$ and $y$ with period $1$, and
can thus be expressed as

$\displaystyle f_{1}(x,y)=-2\pi\ \mathbf{e}_{x}\cdot\nabla G(\mathbf{x}),$
(17)

where $G$ is the pseudo-Green’s function Barton (1989) for the unit-square
with periodic boundary conditions, i.e. the solution of

$\displaystyle-\nabla^{2}G(\mathbf{x})=\delta(\mathbf{x})-1$ (18)

with periodic conditions in $x$ and $y$. It turns out that this pseudo-Green
function is known in closed form Mamode (2014), and is given by

$\displaystyle
G(x,y)=\frac{1}{4\pi}\ln\left\\{e^{-\pi(y^{2}+x^{2})}|\theta_{1}(\pi(x+iy)|i)\theta_{1}(\pi(y+ix)|i)|\right\\}+C,$
(19)

where, in the above, $G$ is defined up to an unimportant additive constant
$C$, $i^{2}=-1$, and $\theta_{1}$ is the elliptic Jacobi theta function,
defined as

$\displaystyle\theta_{1}(z|i)=\theta_{1}(z,q=e^{-\pi})=2e^{-\pi/4}\sum_{n=0}^{\infty}(-1)^{n}e^{-\pi
n(n+1)}\sin[(2n+1)z].$ (20)

Using Eq. (19), the behavior of $G$ near the origin reads

$\displaystyle G(\mathbf{x})\underset{r\to 0}{\sim}\frac{\ln
r}{2\pi}-\frac{r^{2}}{4}+C^{\prime}+A_{4}r^{4}\cos(4\theta)+\mathcal{O}(r^{8}\cos(8\theta)),$
(21)

where $C^{\prime}$ is another unimportant constant, and the coefficient of the
fourth order term reads

$\displaystyle
A_{4}=-\frac{\pi}{16}+\frac{\pi^{3}\sum_{n=0}^{\infty}(-1)^{n}e^{-\pi(n+1/2)^{2}}(2n+1)^{5}}{240\sum_{n=0}^{\infty}(-1)^{n}e^{-\pi(n+1/2)^{2}}(2n+1)}\simeq-0.1253...$
(22)

Using (17), we see that the behavior of $f_{1}$ near the origin is

$\displaystyle f_{1}(\mathbf{x})\underset{r\to
0}{\sim}-\frac{\cos\theta}{r}+\pi r\cos\theta+C^{\prime}-8\pi
A_{4}r^{3}\cos(3\theta)+\mathcal{O}(r^{7}).$ (23)

We now proceed to calculate the inner solution at the next order $\Phi_{1}$.
The boundary condition for $\Phi_{1}$ is deduced from the expansion in $R$ of
Eq. (8),

$\displaystyle\partial_{\tilde{r}}\Phi_{1}|_{\tilde{r}=1}=0$ (24)

Furthermore, $\Phi_{1}$ is harmonic and its behavior at infinity must match
with the behavior of the outer solution (23) near the origin, so that

$\displaystyle\Phi_{1}(\tilde{r},\theta)\underset{\tilde{r}\to\infty}{\sim}\pi\cos\theta\tilde{r}.$
(25)

The solution for $\Phi_{1}$ can therefore be written as

$\displaystyle\Phi_{1}(\tilde{r},\theta)=\pi\cos\theta\left(\tilde{r}+\frac{1}{\tilde{r}}\right).$
(26)

Inserting this expression into the inner expansion (13) and the expression (6)
for $D_{e}$ leads to

$\displaystyle D_{e}=1-\varphi+\varphi^{2}+\mathcal{O}(\varphi^{3})$ (27)

which means that the renormalized formula for $D_{e}$ is exact at second
order. The solution at next order is obtained as follows: the outer function
$f_{2}$ must diverge as $\pi\cos\theta/r$ for small $r$ to match with the
corresponding term in Eq. (26), meaning that $f_{2}=-\pi f_{1}$, which admits
the small $r$ behavior

$\displaystyle f_{2}(r,\theta)\underset{r\to
0}{\sim}\pi\frac{\cos\theta}{r}-\pi^{2}r\cos\theta+...$ (28)

The inner solution $\Phi_{2}$ is

$\displaystyle\Phi_{2}(\tilde{r},\theta)=-8\pi
A_{4}\left(\tilde{r}^{3}+\frac{1}{\tilde{r}^{3}}\right)+\pi\cos\theta\left(\tilde{r}+\frac{1}{\tilde{r}}\right)$
(29)

which is found by requiring that $\Phi_{2}$ is harmonic, with Neumann
condition at $\tilde{r}=1$ and that the cubic term in $\tilde{r}^{3}$ matches
with the corresponding term in the expression of $f_{1}$ [Eq. (23)] and that
the linear term corresponds to the linear term in the expression of $f_{2}$
[Eq. (28)]. Inserting this value into Eq. (6) leads to

$\displaystyle
D_{e}=1-\varphi+\varphi^{2}-\varphi^{3}+\mathcal{O}(\varphi^{4})$ (30)

so that the renormalized expression (2) is still valid at this order, even
though the non-trivial coefficient $A_{4}$ is already present in the inner
solution $\Phi_{2}$.

This calculation can be extended iteratively to higher orders. Since all
functions $f_{n}$ can be constructed with the Green’s function $G$, the
resulting expression for $D_{e}$ therefore depends on the coefficients
appearing in the expansion of the Green’s function. We describe in Appendix A
how to find $D_{e}$ to all orders. The main outcome of this calculation is
that the renormalized result (2) is exact up to (and including) the fourth
order in $\varphi$, i.e.

$\frac{D_{e}}{D_{0}}=\frac{1}{1+\varphi}+\mathcal{O}(\varphi^{5}).$ (31)

and this explains why the above formula is so precise up to seemingly
unreasonably high values of $\varphi$ (up to 0.5). The expression of the first
correction to the renormalized formula can be expressed in terms of the
coefficient $A_{4}$ appearing in the expression of the Green’s function:

$\displaystyle\frac{D_{e}}{D_{0}}=\frac{1}{1+\varphi}-\frac{384\
A_{4}}{\pi^{2}}(\varphi^{5}-\varphi^{6})+\mathcal{O}(\varphi^{7})\simeq\frac{1}{1+\varphi}-0.61165(\varphi^{5}-\varphi^{6})+\mathcal{O}(\varphi^{7})$
(32)

Finally, we also check that the present exact analytical approach is in fact
equivalent to the multipole expansion approach. Using the above mentioned
connection between diffusivity and conductivity, and an analytical series for
the effective conductivity for non-conducting inclusions proposed by Perrins
et al. Perrins et al. (1979) (correcting a result derived by Rayleigh), we see
that the multipole expansion approach for our problem is

$\displaystyle\frac{D_{e}}{D_{0}}=\frac{1}{1-\varphi}\left(1-\frac{2\varphi}{1+\varphi-\frac{0.305827\varphi^{4}}{1-1.402958\varphi^{8}}-0.013362\varphi^{8}}\right).$
(33)

Direct comparison of the above formula with the expansion obtained within our
formalism shows that Eq. (33) is valid up to the order $\varphi^{16}$.

##  IV The limit of nearly touching obstacles in 2D

We now consider the limit where the space between obstacles vanishes ($R\to
L/2$). For the 2D problem, the effective diffusivity vanishes in this limit,
$D_{e}\to 0$. As remarked in Ref. Dagdug et al. (2012), and mentioned above,
adding a ghost boundary condition identical to the zero flux condition at the
boundary of an effective channel does not change the effective diffusivity in
the $x$ direction, which can thus be obtained by considering Brownian motion
inside a channel of local height $h(x)$

$h(x)=\frac{L}{2}-\sqrt{R^{2}-x^{2}}\ \Theta(R-|x|),$ (34)

where $\Theta$ corresponds to the Heaviside function. The above expression is
valid for $-L/2<x<L/2$ and is then repeated with period $L$.

When $R\to R_{c}=L/2$, the effective diffusivity becomes limited by the
crossing of the narrow regions between the obstacles, separated by a length
$a$ at the neck. In these regions, the typical channel height $h$ is small
compared to the longitudinal dimensions $L$. In this limit, it is thus
appropriate to consider an adiabatic type of approximation, where one assumes
a much smaller relaxation time in the lateral direction than in the
longitudinal one. This approximation is known as the Fick-Jacobs
approximation, and, upon its application, the problem is mapped onto a one
dimensional dynamics for the coordinate $x(t)$. The effects of the geometry
are incorporated via an effective potential

$\Phi(x)=-k_{B}T\ln h(x),$ (35)

acting on the particle at the lowest order, and by an effective one
dimensional diffusivity $D(x)$ which is modified at higher orders. The
resulting effective one dimensional diffusion equation, for the marginal
probability distribution along the channel, takes the form

$\partial_{t}p^{*}(x,t)=\partial_{x}\\{D(x)[\partial_{x}p^{*}(x,t)+\beta
p^{*}(x,t)\partial_{x}\Phi(x)]\\}.$ (36)

with $\beta=1/k_{B}T$. We notice that the potential must always take the form
given in Eq. (35) as for a finite system the marginal equilibrium probability
distribution for the variable $x$ is given by

$p^{*}_{\mathrm{eq}}(x)=\frac{h(x)}{\int_{0}^{L}h(x^{\prime})dx^{\prime}}$
(37)

as the, joint, full two dimensional equilibrium distribution is uniform.

The effective diffusivity for the resulting one dimensional diffusion equation
can be computed using the exact Lifson Jackson formula Lifson and Jackson
(1962)

$D_{e}=\frac{1}{\langle h(x)\rangle\langle[D(x)h(x)]^{-1}\rangle},$ (38)

where here $\langle...\rangle$ denotes the uniform average over one period.

At leading order in the limit $h/L\to 0$, the height varies smoothly and the
effective local diffusivity is $D(x)\simeq D_{0}$. This approximation is
called the basic Fick-Jacobs (FJ) approximation. In this case, the averages in
Eq. (38) have been performed analytically by Dagdug et al. Dagdug et al.
(2012), leading to

$\frac{D_{e}}{D_{0}}\simeq\left[(1-\varphi)\left(1-\lambda-\frac{\pi}{2}+\frac{2}{\sqrt{1-\lambda^{2}}}\arctan\sqrt{\frac{1+\lambda}{1-\lambda}}\right)\right]^{-1},$
(39)

where the parameter

$\displaystyle\lambda=2R/L$ (40)

becomes unity when the obstacles touch. This expression recovers the
asymptotic result for $\lambda\rightarrow 1$, first obtained by Keller Keller
(1963),

$\frac{D_{e}}{D_{0}}\simeq\frac{\sqrt{2(1-\lambda)}}{\pi(1-\varphi)}.$ (41)

The next order correction to the Fick-Jacobs approximation (39) in terms of a
one-dimensional description leads to an effective local diffusivity depending
on the local height, the parameter $\varepsilon=h^{\prime}$ is assumed to be
small and the correction to (39) to leading order leads to

$\displaystyle D(x)=D_{0}\left[1-\frac{h^{\prime}(x)^{2}}{3}\right]$ (42)

where the local slope of the channel is considered as a small parameter. The
above formula is exact at order $\varepsilon^{2}=h^{\prime 2}$ (included).
However for the effective height profile generated by a square array of disks
we clearly see that $h^{\prime}(x)$ diverges for $x=\pm R$. (Note that when
the height profile is actually discontinuous, the standard perturbative
approaches to computing the diffusion constant beyond the FJ approximation
have to be modified Mangeat et al. (2018).)

Several partially resummed formulas have been proposed for $D(x)$ that are
compatible with the above expression,

$\displaystyle D(x)=\begin{cases}D_{0}/[1+h^{\prime}(x)^{2}/3]&(\text{Zw})\\\
D_{0}/[1+h^{\prime}(x)^{2}]^{1/3}&(\text{RR})\\\ D_{0}\
\text{arctan}(h^{\prime}(x))/h^{\prime}(x)&(\text{KP})\end{cases}$ (43)

which were proposed by Zwanzig (Zw) Zwanzig (1992), Reguera and Rubí (RR)
Reguera and Rubi (2001) and Kalinay and Percus (KP) Kalinay and Percus (2006),
respectively.

Here we propose an alternative form for $D(x)$, based on the argument that the
generalization of Eq. (42) for the diffusion between two hyper-surfaces of
location $y=\pm h(x_{1},...,x_{n-1})$, where $n$ is the spatial dimension. The
argument is given in the Appendix B. In this case we find that the effective
local diffusion tensor is given by

$\displaystyle D_{ij}(x)\simeq
D_{0}\left[1-\frac{1}{3}(\partial_{x_{i}}h)(\partial_{x_{j}}h)\right].$ (44)

We note that the approximate formulas given in Eq. (43) have the peculiar
properties that the effective diffusion constant tends to zero where
$h^{\prime}(x)$ diverges - that is to say if the channel profile has a kink.
This is clearly unphysical as only a drastic narrowing of the channel can
reduce the diffusion constant. To remedy this problem we recall that the local
normal of the surface has components

$n_{i}(\mathbf{x})=\frac{\partial_{x_{i}}h}{(1+[\nabla h]^{2})^{\frac{1}{2}}}$
(45)

and so to order $\varepsilon^{2}$, Eq. (44) can be written as

$\displaystyle D_{ij}(x)\simeq
D_{0}\left(1-\frac{1}{3}n_{i}(\mathbf{x})n_{j}(\mathbf{x})\right)$ (46)

where $n_{i}(x)$ is the $i^{\text{th}}$ coordinate of the normal to the
surface. Coming back to our 2D problem, this leads us to propose an
alternative form for the effective local diffusion constant $D(x)$ reads

$\displaystyle D(x)\simeq D_{0}\\{1-[n_{x}(x)]^{2}/3\\}$ (47)

where $n_{x}$ is the $x-$ component of the normal to the obstacle,

$\displaystyle
n_{x}(x)=-\frac{h^{\prime}(x)}{\\{1+[h^{\prime}(x)]^{2}\\}^{1/2}}$ (48)

An additional reason to postulate that the use of a local diffusion constant
(47) may lead to more accurate results than its non-renormalized form (42) is
that  the normal to the surface is the natural relevant quantity that appears
in the general equations (6), (8) for the dispersion (whereas the local
quantity $h^{\prime}(x)$ does not appear). Note that with Eq. (47), the local
diffusivity is always finite and positive.

Since $n_{x}(x)=-x/R$ for the disk, we obtain the, everywhere finite,
expression

$\frac{D(x)}{D_{0}}\simeq
1-\frac{1}{3}n_{x}(x)^{2}=1-\frac{x^{2}}{3R^{2}}\Theta(R-|x|).$ (49)

Using Eq. (38) with this local diffusion constant we find the, fully
analytical, result

$\frac{D_{e}}{D_{0}}\simeq\left[(1-\varphi)\left(1-\lambda-\frac{\sqrt{3}\lambda}{2}\frac{\sqrt{2}\pi\lambda-\ln(2+\sqrt{3})}{2\lambda^{2}+1}+\frac{6\lambda^{2}}{(2\lambda^{2}+1)\sqrt{1-\lambda^{2}}}\arctan\sqrt{\frac{1+\lambda}{1-\lambda}}\right)\right]^{-1}.$
(50)

The above expression is the main result of this Section.

Figure 2: Effective diffusivity of Brownian particles in two-dimensional
periodic array of non-attractive obstacles of volume fraction $\varphi$. The
circles represent the numerical solution obtained from Eq. (6-8 ) with a
partial differential equation solver (FreeFem++). (a) Dashed line: leading
order expression (1) in the dilute regime. Continuous line: renormalized
expression (2). The renormalized result is surprisingly valid up to a volume
fraction of $\varphi\sim 0.4$. (b) In the crowded limit, the Fick-Jacobs’
result (39) is represented with dashed lines while other results of the
literature shown by Reguera and Rubí Reguera and Rubi (2001) (RR) and Kalinay
and Percus (KP) Kalinay and Percus (2006) are plotted respectively using
dashed-dotted and dotted lines. Our result, valid only for the periodic array
of obstacles, Eq. (50) is shown using a solid line. The effective diffusivity
is then well-approximated for a volume fraction $\varphi\gtrsim 0.68$ while
existing results in the literature are only accurate down to $\varphi\gtrsim
0.77$.

We now turn to the numerical resolution of Eqs. (6-8) to determine the
accuracy of this new approximate formula, and to compare it with the basic
Fick-Jacobs approximation shown in Eq. (39), and the improvements proposed in
Eqs. (43).

We expect that the basic Fick-Jacobs approximation should work for narrow
channels and so in the region near $\varphi=\pi/4\simeq 0.785$ where the
obstacles touch each other and above which the diffusion constant vanishes. In
Fig.2b we see that Eq. (50) is valid up to $1\%$ precision all the way down to
$\varphi\approx 0.68$, whereas the basic FJ approximation is valid only down
to $\varphi\approx 0.77$ (for the same precision).

Note also that Eq. (50) is more accurate than the leading order asymptotics
(41). Further more, it is also more accurate than all of the modifications
proposed in Eq. (43) for the local diffusivity (and for which we recall that
the averages cannot be performed analytically).

We thus see that the combination of the results (2) and (50) perfectly
describe the effective diffusivity in, respectively, the dilute and the
crowded limits in 2D.

##  V Conclusion

We have revisited the old but important problem of the effective diffusivity
of a Brownian particle diffusing on a square lattice of hard reflecting disks.
This problem is also equivalent to computing the effective conductivity of a
medium of nonconducting spheres embedded in a conducting background as well as
a number of other problems such as flow in porous media and effective
dielectric properties Dean et al. (2007). Using a recently derived Green’s
function for the periodic square lattice Mamode (2014), a perturbative
expansion in the volume fraction of the disks was formulated. This formalism
allows a systematic almost algorithmic determination of the coefficients of
the expansion. We confirm that the first four terms in this expansion have a
particularly simple form which explains the the accuracy of Maxwell’s
effective medium approximation, which is in principle only correct to first
order but actually turns out to capture all terms up to fourth order.

The limit of high volume fraction was treated by mapping the problem onto that
of diffusion in a one dimensional periodic channel as proposed in Dagdug et
al. (2012). Here we proposed a new variant of the second order Fick-Jacobs
approximation where the height derivative terms in the effective diffusivity
are replaced, to the same order, by the expression for the surface normal.
This approximation, in common with a number of other suggestions in the
literature Zwanzig (1992); Reguera and Rubi (2001); Kalinay and Percus (2006),
avoids the appearance of a negative effective local diffusivity. However it
also eliminates the possibility that the effective diffusivity vanishes.
Furthermore, comparison with exact numerical calculations of the effective
diffusivity shows that it describes the behavior of the diffusivity with
better accuracy and down to much lower volume fractions that existing
approximation schemes.

A number of extensions of this work remain open, notably the treatment of
different lattices for hard disks in two dimensions, in both the dilute and
Fick-Jacobs limits. It would be particularly interesting to examine the case
of random arrangements of hard disks. In addition, extensions to the three
dimensional systems with spherical obstacles clearly merits further
exploration in both limits.

##  VI Acknowledgements

This work was partially performed with financial support from the Deutsche
Forschungsgemeinschaft (DFG) within the Collaborative Research Center SFB
1027. The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Calculation of $D_{e}$ at all orders

Here we describe how to calculate $D_{e}$ at all orders in the density of
obstacles. The notation used is the same as Section III. We write the pseudo-
Green’s function $G$ asMamode (2014)

$\displaystyle
G(x,y)=-r^{2}/4+\frac{1}{8\pi}\ln[\psi(z)\psi(i\overline{z})\psi(\overline{z})\psi(-iz)]$
(51)

with

$\displaystyle\psi(z)=\theta_{1}(\pi
z,e^{-\pi})=2e^{-\pi/4}\sum_{n=0}^{\infty}(-1)^{n}e^{-\pi
n(n+1)}\sin[\pi(2n+1)z]$ (52)

We can define coefficients $a_{p}$ which appear in the expansion of $\psi(z)$
near the origin

$\displaystyle\psi(z)=\sum_{p=0}^{\infty}a_{p}z^{p},$ (53) $\displaystyle
a_{2m}=0,\hskip
56.9055pta_{2m-1}=\frac{2(-1)^{m+1}\pi^{2m-1}}{(2m-1)!}e^{-\pi/4}\sum_{n=0}^{\infty}(-1)^{n}e^{-\pi
n(n+1)}(2n+1)^{2m-1}$ (54)

Our goal is to set up an iterative procedure to calculate the expansion of
$D_{e}$ in terms of the coefficients $a_{m}$, which can be straightforwardly
calculated with the series above.

We note that all functions $\Phi_{k}$ (with $k\geq 1$) appearing in the inner
expansion are harmonic and have Neumann boundary conditions, so that general
expression is

$\displaystyle\Phi_{k}(\tilde{r},\theta)=\sum_{q=1}^{k}\gamma_{k,q}\cos((2q-1)\theta)\left(\frac{1}{r^{2q-1}}+r^{2q-1}\right)$
(55)

where the maximal number of terms $k$ follows from the remark that one
additional divergence at most can appear at each iteration.

Next, we construct a set of functions $g_{n}$ which satisfy the following
requirements: (i) $g_{n}$ is harmonic on the unit square with periodic
boundary conditions, and (ii) the behavior of $g_{n}$ near the origin is

$\displaystyle g_{n}(r,\theta)\underset{r\to
0}{\sim}\frac{\cos((2n-1)\theta)}{r^{2n-1}}+\sum_{k=1}^{\infty}\alpha_{n,k}\cos((2k-1)\theta)r^{2k-1}$
(56)

This means that $g_{n}$ contains only one singularity of order $1/r^{2n-1}$
near the origin. Next, we remark that such functions $g_{n}$ can be
constructed by using the pseudo-Green’s function for the unit square with
periodic conditions by using

$\displaystyle g_{1}=2\pi\mathbf{e}_{x}\cdot\nabla G$ (57) $\displaystyle
g_{n+1}=\frac{1}{(2n-1)(2n)}(\mathbf{e}_{x}\cdot\nabla)^{2}g_{n}$
$\displaystyle(n\geq 1)$ (58)

and this recurrence relation implies that

$\displaystyle\alpha_{n+1,k}=\frac{(k+2)(k+1)}{(2n-1)(2n)}\alpha_{n,k+2}.$
(59)

The coefficients $\alpha_{1,l}$ can then be determined from the expansion near
the origin of $G(\mathbf{x})$ by using Eq. (57).

We now express the outer solution as a sum of such functions,

$\displaystyle f_{n}(r,\theta)=\sum_{q\geq 1}A_{n,q}g_{q}(r,\theta)$ (60)

Note that the coefficient of $g_{1}$ must match the $1/r$ singularity in
$\Phi_{n-1}$, the coefficient of $g_{2}$ must match the $1/r^{3}$ singularity
of $\Phi_{n-2}$, and we repeat this process until we find the last function
$\Phi_{n-k}$ that does not contain a $1/r^{2k+1}$ singularity, i.e. the
maximal value of $k$ is the integer part of $n/2$. We thus write

$\displaystyle A_{1,1}=-1$ $\displaystyle(n=1),$ (61) $\displaystyle
A_{n,q}=\gamma_{n-q,q}$ $\displaystyle(n\geq 2,q\leq\text{floor}(n/2)),$ (62)

where the particular case $n=1$ follows from the identification of $\Phi_{1}$.

Next, the coefficients of the inner solutions are determined by the fact that
the behavior of $\Phi_{n}$ must match the linear term $r$ in $f_{n}$, the
cubic term $r^{3}$ in $f_{n-1}$, and more generally the term $r^{2k-1}$ in
$f_{n-k+1}$. Since $f_{n-k+1}$ is a sum over the functions $g_{q}$, we obtain

$\displaystyle\gamma_{n,k}=\sum_{q=1}^{q_{m}(n,k)}A_{n-k+1,q}\alpha_{q,k}$
(63) $\displaystyle q_{m}(n,k)\equiv\text{max}[1,\text{floor}(n-k+1)/2]$ (64)

These formulas enable us to evaluate $D_{e}(\varphi)$ in terms of a series
whose coefficients are calculated iteratively. The first terms are

$\displaystyle D_{e}=1-$
$\displaystyle\varphi+\varphi^{2}-\varphi^{3}+\varphi^{4}-1.611(\varphi^{5}-\varphi^{6})-2.223(\varphi^{7}-\varphi^{8})-3.049\
\varphi^{9}+3.236\ \varphi^{10}$ $\displaystyle-4.248\ \varphi^{11}-4.622\
\varphi^{12}-6.754\ \varphi^{13}+7.446\ \varphi^{14}-9.953\
\varphi^{15}+11.077\ \varphi^{16}$ $\displaystyle-15.201\ \varphi^{17}+17.465\
\varphi^{18}+\mathcal{O}(\varphi^{19}).$ (65)

##  Appendix B Derivation of Equation (44)

Here we derive the approximate relation Eq. (44) which represents the first
correction to the basic Fick-Jacobs approximation. We consider the stochastic
trajectories of a particle $(\mathbf{x}(t),z(t))$ between the surface of
height $z=h(\mathbf{x})$ and the surface $z=0$. Here $z$ denotes the direction
perpendicular to the channel. We start by recalling an exact result for the
marginal probability density function $p^{*}({\bf x},t)$

$p^{*}({\bf x},t)=\int_{0}^{h(\bf x)}dz\ p({\bf x},z,t)$ (66)

for the position ${\bf x}$ along the channel. In Ref. Mangeat et al. (2017b)
it was shown that

$\partial_{t}p({\bf x},t)=D_{0}\nabla\cdot\left(\nabla p^{*}({\bf x},t)-p({\bf
x},h({\bf x}),t)\nabla h({\bf x})\right),$ (67)

where $\nabla$ above indicates the gradient operator parallel to the channel,
that is to say on the coordinates ${\bf x}$. However, the above is not a
closed equation for $p^{*}({\bf x},t)$ due to the presence of the surface term
$p({\bf x},h({\bf x})$ in the advection term.

We now proceed with a perturbative expansion for the full probability density
function $p({\bf x},z)$ based on the fact that $h({\bf x})$, and thus
$z<h({\bf x})$, is small. One writes,

$p({\bf x},z,t)=p_{0}({\bf x},t)+z^{2}p_{2}({\bf x},t)+{\cal O}(z^{4}),.$ (68)

the choice of the expansion in $z^{2}$ being justified by the no-flux boundary
condition at the flat lower surface $z=0$ $\partial_{z}p({\bf
x},z,t)|_{z=0}=0$. The no-flux boundary condition at the upper surface can be
written as

$\left[\nabla p({\bf x},z,t)\cdot\nabla h({\bf x})-\partial_{z}p({\bf
x},z)\right]_{z=h({\bf x})}=0.$ (69)

Inserting Eq. (68) into the above, we find that

$p_{2}({\bf x},t)=\frac{1}{2h({\bf x})}\nabla p_{0}({\bf x},t)\cdot\nabla
h({\bf x}).$ (70)

Integrating Eq. (68) to derive the marginal probability distribution gives

$p^{*}({\bf x},t)=h({\bf x})p_{0}({\bf x},t)+\frac{h^{3}({\bf
x})}{3}p_{2}({\bf x},t)+{\cal O}(h^{5}),$ (71)

and using Eq. (70) then gives to leading order

$p^{*}({\bf x},t)=h({\bf x})p_{0}({\bf x},t)+\frac{h^{2}({\bf x})}{6}\nabla
p_{0}({\bf x},t)\cdot\nabla h({\bf x}).$ (72)

The first order solution to the above is given by

$p_{0}({\bf x},t)=\frac{p^{*}({\bf x},t)}{h({\bf x})},$ (73)

and using this one recovers the basic Fick-Jacobs approximation. Solving Eq.
(72) to next order by iteration gives

$p_{0}({\bf x},t)=\frac{p^{*}({\bf x},t)}{h({\bf x})}-\frac{h({\bf
x})}{6}\nabla\frac{p^{*}({\bf x},t)}{h({\bf x})}\cdot\nabla h({\bf x}).$ (74)

Now using Eqs. (68),(70),(74), we find that the surface term is, to leading
order, given by

$\displaystyle p({\bf x},h({\bf x}),t)=\frac{p^{*}({\bf x},t)}{h({\bf
x})}+\frac{h({\bf x})}{3}\nabla\frac{p^{*}({\bf x},t)}{h({\bf x})}\cdot\nabla
h({\bf x}).$ (75)

Substituting this into Eq. (67) then leads to the closed form equation for
$p^{*}(\mathbf{x},t)$:

$\partial_{t}p^{*}({\bf x},t)=\partial_{x_{i}}\left[D_{ij}({\bf
x})\left(\partial_{x_{j}}p^{*}({\bf x},t)-p^{*}({\bf
x},t)\partial_{x_{j}}\ln(h({\bf x}))\right)\right],$ (76)

where $D_{ij}({\bf x})$ is given by Eq. (44).

## References

  * Bressloff and Newby (2013) P. C. Bressloff and J. M. Newby, Rev. Mod. Phys. 85, 135 (2013).
  * Burada et al. (2009) P. S. Burada, P. Hänggi, F. Marchesoni, G. Schmid, and P. Talkner, ChemPhysChem 10, 45 (2009).
  * Burada et al. (2008) P. S. Burada, G. Schmid, P. Talkner, P. Hänggi, D. Reguera, and J. M. Rubi, BioSystems 93, 16 (2008).
  * Höfling and Franosch (2013) F. Höfling and T. Franosch, Reports on Progress in Physics 76, 046602 (2013).
  * Dentz et al. (2011) M. Dentz, T. Le Borgne, A. Englert, and B. Bijeljic, J. Contam. Hydrol. 120, 1 (2011).
  * Dean et al. (2007) D. S. Dean, I. Drummond, and R. Horgan, J. Stat. Mech. Theory Exp 2007, P07013 (2007).
  * Malgaretti et al. (2013) P. Malgaretti, I. Pagonabarraga, and M. Rubi, Frontiers in Physics 1, 21 (2013).
  * Jacobs (1967) M. Jacobs, _Diffusion processes_ (Springer, New-York, 1967).
  * Yang et al. (2017) X. Yang, C. Liu, Y. Li, F. Marchesoni, P. Hänggi, and H. Zhang, Proc. Natl. Acad. Sci. U. S. A. p. 201707815 (2017).
  * Reguera et al. (2006) D. Reguera, G. Schmid, P. S. Burada, J. M. Rubi, P. Reimann, and P. Hänggi, Phys. Rev. Lett. 96, 130603 (2006).
  * Kalinay and Percus (2006) P. Kalinay and J. Percus, Phys. Rev. E 74, 041203 (2006).
  * Malgaretti et al. (2016) P. Malgaretti, I. Pagonabarraga, and J. Miguel Rubi, J. Chem. Phys. 144, 034901 (2016).
  * Maxwell (1881) J. C. Maxwell, _A Treatise on Electricity and Magnetism_ (Clarendon Press, 1881).
  * Rayleigh (1892) L. Rayleigh, Lond. Edinb. Dublin Philos. Mag. J. Sci. 34, 481 (1892), ISSN 1941-5982.
  * Jeffrey (1973) D. J. Jeffrey, Proc. R. Soc. Lond. A 335, 355 (1973), ISSN 0080-4630, 2053-9169.
  * Batchelor (1976) G. K. Batchelor, J. Fluid Mech. 74, 1 (1976), ISSN 1469-7645, 0022-1120.
  * Kalnin et al. (2001) J. R. Kalnin, E. A. Kotomin, and J. Maier, J. Phys. Chem. Solids 63, 449 (2001).
  * Ward and Keller (1993) M. J. Ward and J. B. Keller, SIA M J. Appl Math 53, 770 (1993).
  * Dagdug et al. (2012) L. Dagdug, M.-V. Vazquez, A. M. Berezhkovskii, V. Y. Zitserman, and S. M. Bezrukov, J. Chem. Phys. 136, 204106 (2012).
  * Berezhkovskii and Szabo (2011) A. Berezhkovskii and A. Szabo, J. Chem. Phys. 135, 074108 (2011).
  * Zwanzig (1992) R. Zwanzig, J Phys. Chem. 96, 3926 (1992).
  * Reguera and Rubi (2001) D. Reguera and J. Rubi, Phys. Rev. E 64, 061106 (2001).
  * Brenner and Edwards (1993) H. Brenner and D. A. Edwards, _Macrotransport theory_ (Butterworth-Heinemann, Boston, 1993).
  * Guérin and Dean (2015a) T. Guérin and D. S. Dean, Phys. Rev. Lett. 115, 020601 (2015a).
  * Guérin and Dean (2015b) T. Guérin and D. S. Dean, Phys. Rev. E 92, 062103 (2015b).
  * Mangeat et al. (2017a) M. Mangeat, T. Guérin, and D. S. Dean, Europhys. Lett. 118, 40004 (2017a).
  * Dean et al. (2004) D. S. Dean, I. Drummond, R. Horgan, and A. Lefevre, J. Phys. A: Math. Gen. 37, 10459 (2004).
  * Batchelor (1974) G. Batchelor, Annual Review of Fluid Mechanics 6, 227 (1974).
  * Barton (1989) G. Barton, _Elements of Green’s functions and propagation_ (Clarendon Press, Oxford, 1989).
  * Mamode (2014) M. Mamode, Boundary Value Problems 2014, 221 (2014).
  * Perrins et al. (1979) W. Perrins, D. R. McKenzie, and R. McPhedran, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 369, 207 (1979).
  * Lifson and Jackson (1962) S. Lifson and J. L. Jackson, J. Chem. Phys. 36, 2410 (1962).
  * Keller (1963) J. B. Keller, J. Appl. Phys. 34, 991 (1963).
  * Mangeat et al. (2018) M. Mangeat, T. Guérin, and D. Dean, J. Chem. Phys 149, 124105 (2018).
  * Mangeat et al. (2017b) M. Mangeat, T. Guérin, and D. S. Dean, J. Stat. Mech. Theory Exp p. 123205 (2017b).


# Bounced Model of Droplet on Moving Substrate

Chengwu Liu<EMAIL_ADDRESS>https://orcid.org/0000-0001-9067-1892 School of Physics, Shandong University,
Jinan 250100, China.  School of Physics, Shandong University, Jinan 250100,
China.

###### Abstract

Firstly, we get the completely bouncing criteria Cr for droplet on moving
substrate. The bouncing without splashing condition is $\mathrm{Cr}>1$. Then,
we mainly research the effect of wind field for droplet, and get the
completely bouncing criteria $\mathrm{Cr}_{\mathrm{wind}}$ for droplet with
wind.

Lastly, we get the contact angle of droplet on the moving substrate and
calculate the Time Independent Reynolds Equation with $rho$ and $\mu$ are
constant.

##  I Introduction

The questions of droplet on a surface are related to the interaction of
interface.

There is a micrometer-size gas film in the interface between liquid and solid.
This gas film was firstly observed by the way of snapshoot [8]. The
evolutionary process of gas film was firstly observed by X-Ray technology [9]
at the moment of contacting. They found that the gas film evolve to a bubble
with spending on microsecond-size time. E. Sawaguchi [10] found that the
distribution of thickness of droplet on a moving surface is similar to saddle
surface. In addition, the hydrophobicity of droplet on a moving surface is
enhanced and is similar to Leidenfrost effect [11]. Therefore, the interaction
between liquid and solid would be affected by the motion of surface. In this
paper, we will talk about how these parameters affect the hydrophobicity in
section 2.

Ted Mao [3] assumed a critical bounced state to deal with the question of
bounced on motionless surface and got a critical bounced criteria
$E_{\mathrm{ERE}^{*}}$. Actually, the gas film on motionless surface is
different from gas film on moving surface. So the interaction of liquid and
solid on motionless surface is also a little different from one on moving
surface. In this paper, we will talk about this question in section 3.

Droplet also might splash on a solid surface. We have some models to describe
splashing [12][13][14][7][15]. In this paper, we will talk about how extra
wind affects the splashing and bounced of droplet in section 4.

##  II Droplet on the moving substrate

The final states of droplet after it impacts on moving substrate are various.
We name bouncing without splashing as completely bouncing, bouncing and
splashing as partially bouncing. We find out many physical phenomenons about
droplet on the moving substrate with lots of experiments(The experiments
conditions are in the supplemental material[1]. ). For instance, there is the
critical speed which is the shift of bouncing and retention. We will deeply
research the completely bouncing condition and critical speed below.

###  II.1 Completely Bouncing Criteria for droplet on Moving Substrate

We will research the phenomenons without splashing in this part.

Droplet will spread, retract, bounce/retain after impacting on the moving
substrate. According to the different interaction modes on the solid-liquid
interface, the interface can be divided into three areas. As the figure 2(a)
shown, the first one is gas film(thickness $h\sim 10\mu\mathrm{m}$). The
cohesive force of gas film acts as normal capillary force to inhibit bouncing.
The second one is solid-liquid interface. The interaction of liquid and solid
is Van der Waals force. The third one is molecule area. The interaction of
this part mainly include the Derjaguin Disjoining Pressure
$\Pi\left(h\right)=\Pi_{\mathrm{vdW}}+\Pi_{\mathrm{EL}}+\Pi_{\mathrm{struc.}}+\Pi_{\mathrm{steric}}$.
In these areas, inhibition and promotion effects for bouncing should be
researched in detial.

####  II.1.1 Inhibition Effects

Firstly, the maximal viscous dissipation of droplet could be approxed is
$D\leq V_{G}=mgh\sim 10^{-3}\mathrm{J}$ with enery conservation(we order the
surface of substrate is the zero gravitational potential surface.) Nextly, the
dissipation of droplet also due to the microstructure on the surface of
substrate. The roughness could be evaluated by the size of roll angle in
general. We could ignore the dissipation of this part with considering a small
roll angle(smooth surface). Acturally, we use the substrate with small roll
angle in our experiments(more detials in supplemental material[1]). Then, we
will consider the interaction between liquid and solid on the gas film and
solid-liquid interface. The attraction is Van der Waals force mainly. And it
is the potential force and plays an obvious role on micrometer scale justly.
Therefore, the vdW force between solid and liquid has no effect on bouncing.
However, some gas are carried in the gas film during the droplet falling
process. This gas film will go away from the system constituted by droplet and
solid substrate during the bouncing process. And the capillary force of gas
film will dissipate the energy of droplet. So, the cohesive force of gas film
will inhibit bouncing. Nextly, we will explain it detially.

First of all, this process is different from the previous section. In this
section, the pressure of gas film $p$ changes continuously over time, and it’s
relatively large. The kinetic equation could be given through analysing the
droplet as the figure 2(b) shown.

$\left(p-p_{\ominus}\right)S=m\frac{\mathrm{d^{2}}H}{\mathrm{d}t^{2}}+mg$ (1)

where $p_{\ominus}$ is atmospheric pressure, $H$ is the vertical positionsl
coordinaten of droplet, $S$ is the contact area of liquid and solid.
Apparently, we could find out $H$ and $S$ change over time through observing
experiments results. So $p$ also should change over time as figure 2 shown.
Then the dissipation due to the cohesive force of gas film is be considered
below.

The gas film is divided into two parts $g_{1},\>g_{2}$ with $h$ thickness as
figure 2(d) shown. And two parts at a distance of $d\sim
10^{-9}\mathrm{m}$($d$ is the distance of two molecule). Hence, the per unit
area vdW potential between $g_{1}$ and $g_{2}$ is

$w_{g}=-\frac{A}{12\pi}\left[\frac{1}{d^{2}}+\frac{1}{\left(d+2h\right)^{2}}-\frac{2}{\left(d+h\right)^{2}}\right]\xlongequal{h\gg
d}-\frac{A}{12\pi d^{2}}$ (2)

where $A\sim 10^{-19}\mathrm{J}$ is Hamker constant. We assume that the
thickness of gas film changes from $h_{0}$ to $h$ from the state of maximum
spreading to bouncing exactly, $\tau^{\prime}$ is the time interval of this
process. Then the energy dissipation during retraction process is

$Q_{\mathrm{cap}}=\left|\int_{h_{0}}^{h}-\frac{\partial w_{g}}{\partial
d}S\mathrm{d}h\right|=\left|\int_{0}^{\tau^{\prime}}\frac{AS}{6\pi
d^{3}}\frac{\mathrm{d}h}{\mathrm{d}t}\mathrm{d}t\right|$ (3)

Then we expand $h$ to linear term and consider the ideal gas hypothesis.

$\frac{\mathrm{d}h}{\mathrm{d}t}=\frac{h-h_{0}}{t},\>\>d^{3}=\frac{kT}{p}$ (4)

We could give the energy dissipation due to capillary force with equation 4.

$\begin{split}Q_{\mathrm{cap}}&\approx\left|\int_{0}^{\tau^{\prime}}-\frac{ApS\left(h-h_{0}\right)}{6\pi
d^{3}t}S\mathrm{d}t\right|\\\
&\xlongequal[Sh_{0}=V_{0}]{Sh=V}\left|\int_{0}^{\tau^{\prime}}\frac{Ap\left(V-V_{0}\right)}{6\pi
kTt}\mathrm{d}t\right|\\\ &\xlongequal[p_{0}V_{0}=\nu RT]{pV=\nu RT}\frac{A\nu
R}{6\pi
k}\left|\int_{0}^{\tau^{\prime}}\frac{1-\left(p/p_{0}\right)}{t}\mathrm{d}t\right|\end{split}$
(5)

This is a improper integral. In order for the improper integral is convergent.
We have

$\lim_{t\to 0+}\frac{1-\left(p/p_{0}\right)}{t}=\lim_{t\to
0+}-\frac{1}{p_{0}}\frac{\mathrm{d}p}{\mathrm{d}t}=a\Rightarrow\frac{\mathrm{d}p}{\mathrm{d}t}=-ap_{0}$
(6)

Therefore, we have the pressure $p$ changes linearly with time below the above
approximation. Then bring equation 4 to equation 5.

$\begin{split}Q_{\mathrm{cap}}&=\frac{A\nu R}{6\pi
k}\left|\int_{p_{0}}^{p^{\ominus}}\frac{\mathrm{d}p}{p_{0}}\right|=\frac{A\nu
N_{\mathrm{A}}}{6\pi}\frac{p_{0}-p^{\ominus}}{p_{0}}=\frac{A\nu
N_{\mathrm{A}}}{6\pi}a\tau^{\prime}\sim 10^{-3}\mathrm{J}\end{split}$ (7)

where $\nu$ is the amount of substrate of gas film, $N_{\mathrm{A}}$ is
Avogadro constant. We could conclude that the order of energy dissipation due
to capillary force is same as the gravitational potential energy. So it
couldn’t be ignored.

Lastly, let’s consider the viscous dissipation. It’s so difficult to calculate
that we have to estimate it. Let’s consider an undemanding toy model.

Firstly, the viscous dissipation should relateto spreading factor with a
certain function. The viscous dissipation which droplet impacts on a smooth
motionless substrate relate to $\left(d_{m}/D\right)^{2}$[2] during spreading
process, and relate to $\left(d_{m}/D\right)^{2.3}$[3] during retraction
process. Hence, we associate the spring oscillator model of two degrees of
freedom with this question(Two springs placed perpendicular to each other
horizontally on a smooth substrate). We assume that the viscous dissipation is
$E_{\mathrm{diss}}\sim QE_{\mathrm{p\>max}}$, where $E_{\mathrm{p\>max}}$ is
the maximum spreading "elastic potential energy" for the first time, $Q$ is
similar to quality factor.

$E_{\mathrm{diss}}\sim\alpha
k_{n}\left(D_{\mathrm{n\>max}}-D_{0}\right)^{2}+\beta
k_{t}\left(D_{\mathrm{t\>max}}-D_{0}\right)^{2}$ (8)

Considering that the effect of tangential mainly due to surface tension and
viscous shear stress $T_{\nu}$, the effect of normal mainly due to surface
tension. We could give the $k_{n}$ and $k_{t}$ using above model.

$\displaystyle k_{\mathrm{t}}$ $\displaystyle\sim\frac{a\gamma
D_{\mathrm{t\>max}}+bT_{\nu}}{D_{\mathrm{t\>max}}-D_{0}}$ (9) $\displaystyle
k_{\mathrm{n}}$ $\displaystyle\sim\frac{a^{\prime}\gamma
D_{\mathrm{n\>max}}}{D_{\mathrm{n\>max}}-D_{0}}$ (10)

where $\alpha,\>\beta,\>a,\>b,\>a^{\prime}$ is constant, $T_{\nu}\sim\eta
VD_{\mathrm{n\>max}}D_{\mathrm{t\>max}}/\delta$ is the viscous shear
stress[4], $\delta$ is the thickness of gas film estimated by LLD’s law[5],
$V$ is the speed of substrate, the maximum tangential spreading diameter on
moving substrate[5] and the maximum spreading diameter on motionless
substrate[6] are respectively
$D_{\mathrm{t\>max}}/D_{0}\sim\mathrm{We}^{1/4}\mathrm{Ca}^{1/6}$ and
$D_{\mathrm{max}}/D_{0}\sim\mathrm{We}^{1/4}$. And the substrate speed has no
effect on normal maximum spreading diameter. Therefore,
$D_{\mathrm{n\>max}}\sim\mathrm{We}^{1/4}$. Then, we could estimate the
viscous dissipation with scaling law.

Figure 1: (a)The interface can be divided into three areas. The first one is
gas film(Zone 1). The second one is solid-liquid interface(Zone 2). The third
one is precursor film(Zone 3). The interaction of this part mainly include the
Derjaguin Disjoining Pressure. In addition to these, energy
dissipation(viscous dissipation and energy dissipation due to capillary force)
also play an important role on bouncing. (b)Pressure isn’t a constant during
the retraction process. (c)Effect of wind field caused by moving substrate on
droplet. This lift force could be ignored when the speed of substrate is
little.(d)The gas film is divided into two parts $g_{1},\>g_{2}$ with $h$
thickness. (e)Completely Bouncing Criteria for droplet on Moving Substrate
$\mathrm{Cr}$. The droplet will retain on the substrate if $\mathrm{Cr}\leq
1$, and bounce on the substrate if $\mathrm{Cr}>1$. (f)Completely Bouncing
Criteria for droplet on Moving Substrate $\mathrm{Cr_{V}}$, which is only due
to the speed of substrate. The arrow points to the critical speed.

####  II.1.2 promotion effects

Firstly, one of the promotion effects is initial kinetic energy
$E_{\mathrm{total}}=\frac{1}{2}\rho\left(\frac{1}{6}\pi
D_{0}^{3}\right)U^{2}$, where $U$ is the speed which droplet impacts the
moving substrate. Next, as the figure 2(c) shown. The movment of solid
substrate caused the flow of air beacuse air is the viscous fluid. Hence, the
pressure around substrate will be decreased by wind field. So, it has a lift
force $F_{L}\sim\rho_{g}U_{t}^{2}h$ on droplet, where
$\rho_{g}=1.185\mathrm{kg/m^{3}}$ is the density of air(25℃),
$\rho=997\mathrm{kg/m^{3}}$ is the density of water(25℃), $U_{t}\sim
1\mathrm{m/s}$ is wind speed around the substrate, $h\sim 10^{-3}\mathrm{m}$
is the maximum spreading thickness of droplet, $D\sim 10^{-3}\mathrm{m}$ is
the initial diameter of droplet. So, we have

$\frac{\rho DU^{2}}{\rho_{g}U_{t}^{2}h}\sim 10^{3}$ (11)

So lift force could be ignored when the speed of substrate is little($U\sim
1\mathrm{m/s}$). The situation of big wind speed will be discussed in section
4.

####  II.1.3 Completely Bouncing Criteria for droplet on Moving Substrate

Hence, we could conclude that the initial kinetic energy $E_{\mathrm{total}}$
promote to bounce, viscous dissipation $E_{\mathrm{diss}}$ and energy
dissipation due to capillary force $Q_{\mathrm{cap}}$ inhibit bouncing.

Considering a imaginary state which droplet bounce exactly. Then, we could get
the condition of bouncing using energy conservarion. We have

$E_{\mathrm{total}}+mg\frac{D}{2}>E_{\mathrm{Diss}}+Q_{\mathrm{cap}}+mg\frac{D}{2}$
(12)

So, we can get the completely bouncing criteria for droplet on moving
substrate.

$\mathrm{Cr}=\frac{\pi\mathrm{We}^{3/2}\sqrt{\gamma/D_{0}\rho}}{\frac{2A\nu
N_{\mathrm{A}}}{\pi
D_{0}^{2}\gamma}\left(\frac{p_{0}-p^{\ominus}}{p_{0}}\right)+12\left[\left(\beta\mathrm{We}^{1/4}\mathrm{Ca}^{1/6}+\frac{\beta^{\prime}D_{0}\mathrm{We}^{1/2}\mathrm{Ca}^{1/2}}{l\mathrm{c}}\right)\left(\mathrm{We}^{1/4}\mathrm{Ca}^{1/6}-1\right)+\alpha\mathrm{We}^{1/4}\left(\mathrm{We}^{1/4}-1\right)\right]}$
(13)

where $l_{c}=\left(\gamma/\rho g\right)^{-1/2}$ is the Capillary length of
water, $rho$ is the density of liquid, $\gamma$ is the liquid-gas surface
tension, $\alpha,\>\beta,\>\beta^{\prime}$ are the constants to be determined.
The numerator is shown that the effect of initial kinetic energy. The left
term of denominator is shown that the effects of capillary force in the gas
film. The right term of denominator is shown that the effects of viscous
dissipation. The bouncing without splashing condition is

$\mathrm{Cr}>1$ (14)

And completely bouncing criteria for droplet on moving substrate
$\mathrm{Cr_{V}}$, which is only due to the speed of substrate.

$\mathrm{Cr_{V}}=\frac{\mathrm{A}}{\mathrm{B}+\left(\mathrm{C}V^{1/6}+\mathrm{C^{\prime}}V^{1/2}\right)\left(\mathrm{D}V^{1/6}-1\right)}$
(15)

where A,B,C,C’,D are indepent to substrate speed $V$. The droplet will retain
on the substrate if $\mathrm{Cr_{V}}\leq 1$, and bounce on the substrate if
$\mathrm{Cr_{V}}>1$.

In a word, the completely bouncing criteria for droplet on moving substrate is
relate to capillary number $\mathrm{Ca}$ and Weber number $\mathrm{We}$ which
both are the initial state parameters of the droplet and moving substrate, as
figure 2(e) shown. And we also explain the experimental phenomenons which the
final state of droplet will shift with changing substrate speed as figure 2(f)
shown. We also could find out the critical speed which is the shift of
bouncing and retention as figure 2(f) shown. These results show that droplet
might undergo the transformation of "Retention to Bounced" or "Retention to
Bounced to Retention". But we didn’t consider the dissipation due to the
microstructure on the surface of substrate. It’s also a complex question

##  III The Effect of Wind Field for Droplet

In the discussion of previous section, we find that a lift force due to the
wind field couldn’t be ignored with high substrate speed. The lift force
promote to bounce. We even find out the splashingofdroplet in further
experiments. So we design the experiment(more detials are in supplemental
material[1]) to illustrate the importance of wind field for the final state of
droplet as figure 3(a)(b)shown.

Figure 2: (a)The arrow represents the direction in which time is passing. It
shows the droplet state with wind speed is 0, 3.8m/s and 10m/s. (b)The finial
state figure of dropplet. Three points which have various states respectively
represent that splashing+partially bouncing and completely
bouncing(coordinates are (7,5) and (9,8)), various states(coordinate is
(5,7)). (c)The spreading of droplet which inital height is 5cm in wind field.
The red points are the initial position of the center of mass(d)Considering
the loops $C_{\mathrm{droplet}}$ and $C_{\mathrm{air}}$. (e)R&G model with
wind field. $F_{\mathrm{wind,\>L}}$ is the extra effect on liquid finger.

We conclude that the final state of droplet is affected by both the initial
height of droplet(We) and the speed of wind(Ca) from figure 3(b). And the
state of droplet will shift to splashing and partially bouncing with creasing
the wind speed and the initial height. The state of droplet will shift to
splashing with only creasing the initial height. The state of droplet with
small initial height will shift to splashing and partially bouncing if only
creasing the wind speed. And state also could shift from retention to
completely bouncing/from completely bouncing to retention. We conclude that
the spreading factor et. al. could change with the wind speed from figure
3(c). Next, we will research the causes of these phenomenons.

###  III.1 The Transition Between Retention and Completely Bouncing

Firstly, we assume that (1)all flows could be ragard as the isentropic flows.
(2)The boundary conditions between liquid and gas obey Navier boundary
conditions, i.e., $\vec{V}_{\mathrm{droplet}}=\vec{V}_{\mathrm{air}}$(3)Wind
field is 2D incompressible laminar flow, i.e.,
$\nabla\cdot\vec{V}=0,\>V_{\mathrm{z}}=\mathrm{const.}$.

We take two circuits $C_{\mathrm{droplet}}$ and $C_{\mathrm{air}}$ around the
interface between liquid and gas. The velocity circulations respectively are

$\Gamma_{\mathrm{droplet}}=\oint\limits_{C_{\mathrm{droplet}}}\vec{V}_{\mathrm{droplet}}\cdot\mathrm{d}\vec{l}=\Gamma_{\mathrm{air}}=\oint\limits_{C_{\mathrm{air}}}\vec{V}_{\mathrm{air}}\cdot\mathrm{d}\vec{l}=\Gamma$
(16)

These velocity circulations $\Gamma$ are constant because of assumation 2. And
the interaction of wind field on the droplet could be divided into vertical
and horizontal direction. The bouncing state mainly depends on the vertical
interaction.

$\vec{L}=\rho_{a}\vec{V}_{\infty}\times\vec{\Gamma}_{\mathrm{air}}=-\rho_{a}\Gamma\vec{V}_{\infty}\times\vec{k}$
(17)

It can be seen from assumation 3

$\nabla\times\vec{L}=-\rho_{a}\Gamma\left[\left(\vec{k}\cdot\nabla\right)\vec{V}_{\infty}+\left(\nabla\cdot\vec{k}\right)\vec{V}_{\infty}-\left(\vec{V}_{\infty}\cdot{\nabla}\right)\vec{k}-\left(\nabla\cdot\vec{V}_{\infty}\right)\vec{k}\right]=0$
(18)

So the lift force $\vec{L}$ is potential force. We assmue the initial height
of droplet is $h_{0}$ and the lift force potential of initial position is 0.
Then the lift force potential is

$V_{\mathrm{L}}\left(y\right)=-\int-\rho_{a}\Gamma\vec{V}_{\infty}\times\vec{k}\cdot\vec{j}\mathrm{d}y=-\rho_{a}\Gamma\int
V_{\infty\mathrm{x}}\mathrm{d}y=\rho_{a}\Gamma\int_{h_{0}}^{y}V_{\infty\mathrm{x}}\mathrm{d}y$
(19)

So the completely bouncing criteria for droplet is

$\mathrm{Cr_{wind}}=\frac{\rho g\pi
D_{0}^{3}\left(h_{0}-0.5D_{0}\right)/6}{E_{\mathrm{Diss}}+Q_{\mathrm{cap}}+\rho_{a}\Gamma\int_{h_{0}}^{D_{0}/2}V_{\infty\mathrm{x}}\mathrm{d}y}$
(20)

The bouncing without splashing condition is

$\mathrm{Cr_{wind}}>1$ (21)

The final state of droplt could shift between bouncing and retention beacuse
$\Gamma\int_{h_{0}}^{D_{0}/2}V_{\infty\mathrm{x}}\mathrm{d}y$ could
bigger/smaller than 0.

###  III.2 The Transition Between Splashing and without Splashing

Then, considering the transition between splashing and without splashing. The
front of droplet maybe generate the liquid finger during the spreading
process. Liquid finger will be not only affected by the lubrication force of
bottom gas and the attraction of top
gas[7]$F_{L}=K_{l}\nu_{g}V_{t}+K_{u}\rho_{g}V_{t}^{2}H_{t}$, but also affected
by the wind field as figure 3(e) shown. The effect of extra wind field is
shown as lift force $F_{\mathrm{wind,L}}\left(V\right)$, where $V$ is the
speed of wind speed. Hence, we could introduce the
$F_{\mathrm{wind,L}}\left(V\right)$ to R&G model.

$\beta^{2}_{\mathrm{Wind}}=\frac{F_{\mathrm{L}}+F_{\mathrm{wind,L}}}{2\gamma}$
(22)

So, the state of droplet maybe have the transition between splashing and
without splashing with the wind field.

##  IV The Balance of Droplet on the Moving Substrate

###  IV.1 The Contact Angle of Droplet

The bottom of droplet will generate a very thin gas film when it impact on a
substrate[9]. The gas film will be saddle shape stably on moving
substrate[10]. The pressure of gas film is $p\sim 10\mathrm{Pa}$[10], the gas
film thickness is $\delta\sim 10\mathrm{\mu m}$. We also could approx the
Knudsen number of gas film is $K_{n}=\frac{\lambda}{h}\sim 10^{-4}$. So the
gas film could be regard as the continuons flow. If we assume that pressure is
a constance on the direction perpendicular to moving substrate, the thickness
and pressure of lubrication gas obey the Reylond Equation:

$\frac{\partial}{\partial x}(\frac{h^{3}}{\mu}\frac{\partial p}{\partial
x})+\frac{\partial}{\partial y}(\frac{h^{3}}{\mu}\frac{\partial p}{\partial
y})=6U\frac{\partial h}{\partial x}$ (23)

where $U$ is the speed of moving substrate, $\mu$ is the dynamic viscosity of
lubrication gas. This equation shows that the relation of distribution between
thickness and pressure. We could give the distribution of pressure with
measuring the distribution of thickness.

Figure 3: The diagrammatic map of gas film.

Then we research the effect of contact angle on moving surface. Considered the
difference form gas film on motionless and moving substrate originates from
the moving substrate. So, we focus on this element, then analyse this question
with minimum energy principle.

We order that the area of gas film is D and O is the center of area D, as
shown in the figure1, we research the infinitesimal area $\mathrm{D_{i}}$ with
infinitesimal angle and radial length $R_{i}$. It is full of air and saturated
vapour in area $\mathrm{D_{i}}$, the pressure of vapour obeys the Clapeyron
Equation $p=p_{0}\exp(-L_{\mathrm{v,m}}/RT)$. And the molecular number and
pressure of two components(saturated vapour and air) obey that

$1=\frac{N_{\mathrm{air}}}{N}+\frac{N_{\mathrm{H_{2}O}}}{N},\>1=\frac{p_{\mathrm{air}}}{p}+\frac{p_{\mathrm{H_{2}O}}}{p}$
(24)

The mean kinetic energy of two components can be given that
$\bar{e}_{\mathrm{air}}=\frac{5}{2}kT,\>\bar{e}_{\mathrm{H_{2}O}}=3kT$, if air
is regarded as diatomic molecule. So we could give the internal energy of gas
film

$E_{k}=\sum_{i=1}^{n}\iint\limits_{D_{i}}nh\bar{e}\mathrm{d}\sigma=\sum_{i=1}^{n}\iint\limits_{D_{i}}\frac{5}{2}ph+h\frac{p_{0}}{2}\exp(-\frac{L_{v,m}}{RT})\mathrm{d}\sigma$
(25)

Hence, assumed that the front of droplet has a infinitesimal virtual
displacement $\delta R_{i}$. Approximately, the interfacial energy between
solid and gas remain unchanged because of the gas film in the solid-gas
interface. So the variation in energy of system is $\delta E_{i}=(\Delta
L_{i}\cos\theta_{L}\gamma_{LG}+\Delta L_{i}\gamma_{SL})\delta R_{i}+\delta
E_{ki}$. And a stabilized system must obey that

$\frac{\delta E_{i}}{\delta R_{i}}=0$ (26)

A combination of equation 25 and equation 26 leads to

$\cos\theta_{Li}=\cos\theta_{0}-\frac{1}{\gamma_{LG}}\left[\frac{5}{2}ph+h\frac{p_{0}}{2}\exp(-\frac{L_{v,m}}{RT})\right]-\frac{\gamma_{SG}}{\gamma_{LG}}$
(27)

where $p=p(R_{i}\cos\theta_{\Delta L_{i}},R_{i}\sin\theta_{\Delta
L_{i}}),h=h(R_{i}\cos\theta_{\Delta L_{i}},R_{i}\sin\theta_{\Delta L_{i}})$,
i.e., $p$ and $h$ are the pressure and thickness of soild-liquid interface
boundary respectively.

Hence, we can get the mean contact angle with intergrating contact angle
$\cos\theta_{Li}$ along the boundary.

$\cos\theta_{L}=\cos\theta_{0}-\frac{1}{L\gamma_{LG}}\oint_{L}\left[\frac{5}{2}ph+h\frac{p_{0}}{2}\exp(-\frac{L_{v,m}}{RT})\right]\mathrm{d}l-\frac{\gamma_{SG}}{\gamma_{LG}}$
(28)

where $\theta_{0}$ is the contact angle obeying the Young Equation, $L$ is the
circumference of solid-liquid interface boundary, $L_{\mathrm{v,\>m}}$ is the
latent heat of phase transition from liquid to gas phase. $T$ is the
temperature of gas film. The element on the right side of equation 28 is the
influence of gas film. The element on the middle of equation 28 is the
influence of moving substrate.

We could conclude that the contact angle on moving substrate is 8∘ bigger than
the one obeyed Young Equation in the room tempurature approximately. So the
hydrophobicity will be reinforced on the moving substrate. And the contact
time[16][5], spreading factor[5], bouncing et.al. will change with the change
of hydrophobicity between the droplet and substrate. Then we will elucidate
them in detial.

###  IV.2 The analytical solution for Reynolds Equation

Firstly, we could get another equation which describes the gas film on the
moving surface from Reynolds transport equation:

$\frac{\partial h}{\partial
t}+\nabla\cdot\left(h\mathbf{u}\right)=0,\>\frac{\partial h}{\partial t}=0$
(29)

where $\mathbf{u}$ could be seen as the surface speed
$U\mathbf{i}+V\mathbf{j}$ because of the Navier Boundary Conditions.

Then, the problem is solving the differential equations:

$\left\\{\begin{matrix}h\frac{\partial^{2}p}{\partial
x^{2}}+h\frac{\partial^{2}p}{\partial y^{2}}+3\left(\frac{\partial h}{\partial
x}\frac{\partial p}{\partial x}+\frac{\partial h}{\partial y}\frac{\partial
p}{\partial y}\right)=0\\\ U\frac{\partial h}{\partial x}+V\frac{\partial
h}{\partial y}=0\par\end{matrix}\right.$ (30)

Then, we order that
$h\left(x,y\right)=h_{X}\left(x\right)h_{Y}\left(y\right)$,
$p\left(x,y\right)=p_{X}\left(x\right)p_{Y}\left(y\right)$. We could get
equation 31 through bringing these to equation 30.

$\frac{1}{p_{X}^{2}}\frac{\mathrm{d}^{2}p_{X}}{\mathrm{d}x^{2}}+\frac{1}{p_{Y}^{2}}\frac{\mathrm{d}^{2}p_{Y}}{\mathrm{d}x^{2}}+\frac{3}{p_{X}^{2}p_{Y}h_{X}}\frac{\mathrm{d}h_{X}}{\mathrm{d}x}\frac{\mathrm{d}p_{X}}{\mathrm{d}x}+\frac{3}{p_{Y}^{2}p_{X}h_{Y}}\frac{\mathrm{d}h_{Y}}{\mathrm{d}y}\frac{\mathrm{d}p_{X}}{\mathrm{d}y}=0$
(31)

Then, finding the derivative of equation 31 with respect to x and y in turn.
We can get

$\frac{\mathrm{d}^{2}p_{X}}{\mathrm{d}x^{2}}-C\frac{\mathrm{d}p_{X}}{\mathrm{d}x}+C^{\prime}p_{X}^{2}=0$
(32)
$\frac{\mathrm{d}p_{Y}}{\mathrm{d}y}\frac{\mathrm{d}h_{Y}}{\mathrm{d}y}-\frac{Cp_{Y}^{2}h_{Y}}{3}=C^{\prime}\frac{p_{Y}h_{Y}}{3}$
(33)

$C,\>C^{\prime}$are constant. And we could esaily find the solution of
equarion 30.

$h_{X}=C_{h_{1}}\exp\left(-\frac{\lambda}{U}x\right),\>h_{Y}=C_{h_{2}}\exp\left(\frac{\lambda}{V}y\right)$
(34)

bring them to equation 33, we can get

$\int\frac{\mathrm{d}p_{Y}}{C_{2}^{\prime}p_{Y}+C_{2}p_{Y}}=y$ (35)

So the $p_{Y}$ is

$p_{Y}=\frac{C_{2}}{-C_{2}^{\prime}+\exp\left(\frac{-y+C_{2}^{\prime\prime}}{C_{2}}\right)}$
(36)

Then, we solve the equation 32 with series method. Considering the series
solution $p_{X}=\sum_{n=0}^{\infty}a_{n}x^{n}$. Then, we can get

$a_{n+2}\left(n+2\right)\left(n+1\right)-Ca_{n+1}\left(n+1\right)+2C^{\prime}\left(a_{0}a_{n}+a_{1}+a_{n-1}+\cdots+a_{n/2}a_{n/2}\right)=0,\>\mathrm{n\>is\>an\>even.}$
(37)
$a_{n+2}\left(n+2\right)\left(n+1\right)-Ca_{n+1}\left(n+1\right)+2C^{\prime}\left(a_{0}a_{n}+a_{1}+a_{n-1}+\cdots+a_{\left(n-1\right)/2}a_{\left(n+1\right)/2}\right)=0,\>\mathrm{n\>is\>an\>odd.}$
(38)

And the radius of convergence $R$ obey that

$\lim_{n\to\infty}\left[n+2+\frac{2C^{\prime}}{n+1}\left(a_{0}R^{2}+a_{1}R^{3}+\cdots+a_{n/2}R^{n/2+2}\right)\right]-CR=0,\>\mathrm{n\>is\>an\>even.}$
(39)
$\lim_{n\to\infty}\left[n+2+\frac{2C^{\prime}}{n+1}\left(a_{0}R^{2}+a_{1}R^{3}+\cdots+a_{\left(n-1\right)/2}R^{\left(n+3\right)/2}\right)\right]-CR=0,\>\mathrm{n\>is\>an\>odd.}$
(40)

In addition, we can get equation 41 when $R=1$.

$\begin{split}\lim_{n\to\infty}\left[n+2+\frac{2C^{\prime}q\left(n/2+1\right)}{n+1}-C\right]\leq&\lim_{n\to\infty}\left[n+2+\frac{2C^{\prime}}{n+1}\left(a_{0}R^{2}+a_{1}R^{3}+\cdots+a_{n/2}R^{n/2+2}\right)\right]-CR\\\
&\leq\lim_{n\to\infty}\left[n+2+\frac{2C^{\prime}q^{\prime}\left(n/2+1\right)}{n+1}-C\right],\>\mathrm{n\>is\>an\>even.}\\\
\lim_{n\to\infty}\left[n+2+\frac{2C^{\prime}w\left(\left(n+1\right)/2\right)}{n+1}-C\right]\leq&\lim_{n\to\infty}\left[n+2+\frac{2C^{\prime}}{n+1}\left(a_{0}R^{2}+a_{1}R^{3}+\cdots+a_{\left(n-1\right)/2}R^{\left(n+3\right)/2}\right)\right]-CR\\\
&\leq\lim_{n\to\infty}\left[n+2+\frac{2C^{\prime}w^{\prime}\left(\left(n+1\right)/2\right)}{n+1}-C\right],\>\mathrm{n\>is\>an\>odd.}\\\
\end{split}$ (41)

where
$q=\min\left\\{a_{0},\>a_{1},\cdots,\>a_{n/2}\right\\},\>q^{\prime}=\max\left\\{a_{0},\>a_{1},\cdots,\>a_{n/2}\right\\}$
and
$w=\min\left\\{a_{0},\>a_{1},\cdots,\>a_{\left(n-1\right)/2}\right\\},\>w^{\prime}=\max\left\\{a_{0},\>a_{1},\cdots,\>a_{\left(n-1\right)/2}\right\\}$.

If the equation 41 is right, the equation 39 and 40 would be wrong. So the $R$
is either $\infty$ or $1<R<\infty$.

So $p\left(x,y\right)$ is

$p=\left(a_{0}+a_{1}x+\frac{Ca_{1}-C^{\prime}a_{0}^{2}}{2}x^{2}+\cdots\right)\left(\frac{C_{2}}{-C_{2}^{\prime}+\exp\left(\frac{-y+C_{2}^{\prime\prime}}{C_{2}}\right)}\right),\>1<x\leq\infty$
(42)

So $h\left(x,\>y\right)$ is

$h\left(x,\>y\right)=h_{X}\cdot
h_{Y}=C_{h}\exp\left(-\frac{\lambda}{U}x+\frac{\lambda}{V}y\right)$ (43)

##  V Conclusion

In section 2, we find that how the moving substrate affect the hydrophobicity
of droplet and we discuss the analytical solution for Reynolds Equation.
Therefore, we would analytically get the contact angle on moving substrate.
But we must have some boundary conditions such as
$h\left(x,\>y\right)|_{\mathrm{droplet\>boundary}}=H\left(x,\>y\right),\>p\left(x,\>y\right)|_{\mathrm{droplet\>boundary}}=P\left(x,\>y\right)$
and so on to get the whole solution. In section 3, we find out some promotion
and inhibition effects for bouncing question. Finally, we get a completely
bouncing criteria $\mathrm{Cr}$ for droplet on moving substrate. Some
phenomenons could be pridicted by using this criteria. In section 4, we
research the effect of extra wind field for droplet. we find that extra wind
field could change the final states of droplet. In addition, we get a
completely bouncing criteria $\mathrm{Cr_{Wind}}$ for droplet on moving
substrate with extra wind by the way that introduce the lift force potential.
We also get the splashing criteria $\beta^{2}_{\mathrm{Wind}}$ using the R&G
Model.

But we don’t analytically get a criteria because of the complex viscous
dissipation. We just use a simple model to calculate it. We also don’t
consider the energy dissipation due to the roughness of substrate. It is so
important that cound’t be ignored in some substrate with big roughness.

###### Acknowledgements.

Thanks for Shangqian Sun, Hongwang Lu, Jingcheng Hao and Ying Ma ’s support
for this work.

## References

  * Chandra and Avedisian [1991a] S. Chandra and C. T. Avedisian, On the collision of a droplet with a solid surface, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 432, 13 (1991a), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1991.0002 .
  * Chandra and Avedisian [1991b] S. Chandra and C. T. Avedisian, On the collision of a droplet with a solid surface, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 432, 13 (1991b), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1991.0002 .
  * Mao  _et al._ [1997] T. Mao, D. C. S. Kuhn, and H. Tran, Spread and rebound of liquid droplets upon impact on flat surfaces, AIChE Journal 43, 2169 (1997), https://aiche.onlinelibrary.wiley.com/doi/pdf/10.1002/aic.690430903 .
  * Tran  _et al._ [2013] T. Tran, H. J. J. Staat, A. Susarrey-Arce, T. C. Foertsch, A. van Houselt, H. J. G. E. Gardeniers, A. Prosperetti, D. Lohse, and C. Sun, Droplet impact on superheated micro-structured surfaces, Soft Matter 9, 3272 (2013).
  * Zhan  _et al._ [2021] H. Zhan, C. Lu, C. Liu, Z. Wang, C. Lv, and Y. Liu, Horizontal motion of a superhydrophobic substrate affects the drop bouncing dynamics, Phys. Rev. Lett. 126, 234503 (2021).
  * CLANET  _et al._ [2004] C. CLANET, C. BÉGUIN, D. RICHARD, and D. QUÉRÉ, Maximal deformation of an impacting drop, Journal of Fluid Mechanics 517, 199–208 (2004).
  * Riboux and Gordillo [2014] G. Riboux and J. M. Gordillo, Experiments of drops impacting a smooth solid surface: A model of the critical impact speed for drop splashing, Phys. Rev. Lett. 113, 024507 (2014).
  * Chandra and Avedisian [1991c] S. Chandra and C. T. Avedisian, On the collision of a droplet with a solid surface, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 432, 13 (1991c), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1991.0002 .
  * Lee  _et al._ [2012] J. S. Lee, B. M. Weon, J. H. Je, and K. Fezzaa, How does an air film evolve into a bubble during drop impact?, Phys. Rev. Lett. 109, 204501 (2012).
  * Sawaguchi  _et al._ [2019] E. Sawaguchi, A. Matsuda, K. Hama, M. Saito, and Y. Tagawa, Droplet levitation over a moving wall with a steady air film, Journal of Fluid Mechanics 862, 10.1017/jfm.2018.952 (2019).
  * Gauthier  _et al._ [2016] A. Gauthier, J. C. Bird, C. Clanet, and D. Quéré, Aerodynamic leidenfrost effect, Phys. Rev. Fluids 1, 084002 (2016).
  * Mundo  _et al._ [1995] C. Mundo, M. Sommerfeld, and C. Tropea, Droplet-wall collisions: Experimental studies of the deformation and breakup process, International Journal of Multiphase Flow 21, 151 (1995).
  * Xu [2007] L. Xu, Liquid drop splashing on smooth, rough, and textured surfaces, Phys. Rev. E 75, 056316 (2007).
  * Mandre  _et al._ [2009] S. Mandre, M. Mani, and M. P. Brenner, Precursors to splashing of liquid droplets on a solid surface, Phys. Rev. Lett. 102, 134502 (2009).
  * Hao and Green [2017] J. Hao and S. I. Green, Splash threshold of a droplet impacting a moving substrate, Physics of Fluids 29, 012103 (2017), https://doi.org/10.1063/1.4972976 .
  * Zhang  _et al._ [2020] X. Zhang, Z. Zhu, C. Zhang, and C. Yang, Reduced contact time of a droplet impacting on a moving superhydrophobic surface, Applied Physics Letters 117, 151602 (2020), https://doi.org/10.1063/5.0023896 .