PHYSICAL REVIEW E 103, 042607 (2021)

Brownian motion of a charged colloid in restricted confinement
Yael Avni ,1,2 Shigeyuki Komura ,2 and David Andelman
11Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Ramat Aviv 69978, Tel Aviv, Israel
2
Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University, Tokyo 192-0397, Japan(Received 13 December 2020; accepted 23 March 2021; published 19 April 2021)
We study the Brownian motion of a charged colloid, confined between two charged walls, for small separation
between the colloid and the walls. The system is embedded in an ionic solution. The combined effect ofelectrostatic repulsion and reduced diffusion due to hydrodynamic forces results in a specific motion in the
direction perpendicular to the confining walls. The apparent diffusion coefficient at short times as well as the
diffusion characteristic time are shown to follow a sigmoid curve as a function of a dimensionless parameter.
This parameter depends on the electrostatic properties and can be controlled by tuning the solution ionic strength.
At low ionic strength, the colloid moves faster and is localized, while at high ionic strength it moves slower andexplores a wider region between the walls, resulting in a larger diffusion characteristic time.
DOI: 10.1103/PhysRevE.103.042607

I. INTRODUCTION

Understanding the Brownian motion of colloids under confinement has been a great challenge in recent decades [1].
Such motion is present in the segregation and transport ofparticles through biocellular membranes [2], in microfluidic devices [3], and in particle trapping and tracking
technologies [4,5].
A simple example of confinement is observed when acolloid is placed in a confined space bounded by rigid surfaces. In this case, the colloid motion differs from its free
medium one (no boundaries) due to hydrodynamic forces between the colloid and the confining surfaces. At low Reynolds
numbers, this effect reduces to a change in the drag coefficient, which in turn leads to a position-dependent diffusion
coefficient (PDDC) of the colloid. For simple geometries
such as a spherical colloid near a single planar wall or in
between two flat walls, the PDDC was calculated by Brennerand co-workers [6,7] and Ganatos et al. [8,9]. These results
were verified experimentally for micro- or nanometer-sizedparticles, by different techniques of light scattering [10–13],
video microscopy [14], optical tweezers [15] and total internalreflection velocimetry [16].
In polar solvents, the colloid and the confining surfaces
often carry an electric charge, and the interaction between
them confines the colloid even further to the vicinity ofsome potential minimum [17]. This results in an intricate
motion which is affected both by the PDDC and the confining
potential.
Consequently, the colloid motion was examined experimentally in the lateral direction [18,19]. However, its motionin the direction perpendicular to the confining surfaces, where
the coupling between the PDDC and the interaction potential
is pronounced, was only examined when the diffusion coefficient does not change throughout the colloid motion [20,21].
Although this simplification is valid at times, a rigorous de2470-0045/2021/103(4)/042607(9)scription of the motion in which the PDDC and interaction
potential are fully coupled to one another is still missing.
Moreover, in previous studies, the electrostatic interaction
was approximated by the simplified form, U ∼ exp (−d/λD ),where λD is the Debye screening length and d is the separation
distance between the surface of the colloid and the wall.
While this is valid in the large separation limit d  λD , theinteresting regime where d is of the order of λD has yet to be
explored.
In this paper, we study the motion of a spherical charged
colloid, of radius a, confined between two charged walls,
and focus on the motion in the direction perpendicular to
the walls. We consider the large sphere limit, i.e., a  λD , d,for which the calculation greatly simplifies and analyticalresults are obtained. We make an additional simplification
by approximating the interaction potential to be a harmonic
one around the equilibrium position. We obtain a dynamical
equation that depends on a single dimensionless parameter,α, that quantifies the interplay between the walls and theinteraction potential.
In the small or large α limits, the calculated probability distribution function agrees with previous results, whilefor α ∼ 1, the motion deviates substantially from these two
known cases. The difference is quantified by the mean-squaredisplacement (MSD) behavior as a function of time, wherean analytic expression for the short-time behavior is derived,
while the long-time behavior is studied numerically. We also
show the dependence of α on the electrostatic properties of thesystem. In particular, we demonstrate how it can be tuned bychanging the salt concentration, causing the colloid motion
to cross over between the PCCD dominated motion, whichis slow and explores a wide region, and the electrostaticallydominated one, which is faster and more localized.
The outline of the paper is as follows. In Sec. II, we
calculate the interaction potential. In Sec. III, we calculate
the PDDC and derive the dynamical equation of the colloid

042607-1©2021 American Physical Society

AVNI, KOMURA, AND ANDELMANPHYSICAL REVIEW E 103, 042607 (2021)

FIG. 1. Schematic drawing of the system. A colloid of surface
charge density σc and radius a is positioned in between two wallswith distance L between them, each having a surface charge density
σw . Both the colloid and walls are embedded in an ionic solution
with dielectric constant ε, viscosity η, temperature T , and screening
length λD .

kB is the Boltzmann constant. The equation is solved together
with the boundary conditions of fixed charge on the colloid
surface and on the walls.
The DH equation does not have an analytic solution for thegeometry considered here. However, by assuming that a is thelargest length scale in the system, a  d , dr , λD (the large
sphere limit), Eq. (1) can be solved by using the Derjaguin
approximation. Note that the large sphere limit combines two
different and independent physical limits: a  d , dr (narrowconfinement) and a  λD (thin double layer).
In Appendix A, we calculate the electrostatic potential ψ inthe large sphere limit. In addition, the Appendix contains thedetailed derivation of the interaction potential, U (x), between
the colloid and the charged walls, in the same limit. For
convenience, we repeat here only the final expression relating
U (x) and the osmotic pressure, , d0 +x  ∞U (x) = − 2π adhdl (l )+h ∞

 d0 −xdhdl(l ) ,(2)hmotion. In Sec. IV, we analyze the MSD of the motion fordifferent values of the parameter α, and show the dependenceof α on the electrostatic properties of the system. In Sec. V,
we conclude with several general observations and commentson the validity of our results.
II. ELECTROSTATIC POTENTIAL

We consider a spherical colloid of radius a and charge Q,
embedded in a dilute ionic solution. The charge of the colloid
is assumed to be distributed homogeneously on its surface,with surface charge density σc = Q/(4π a2 ), where σc canbe either positive or negative. The ionic solution is characterized by a dimensionless dielectric constant ε, viscosity η,
temperature T , and bulk concentration of monovalent salt
n. The colloid is confined between two charged walls, as
shown in Fig. 1. The distance between the walls is L, and theposition of the center of the colloid relative to the midplane
between the two walls is x. The distance between the colloid
and the left (right) wall is denoted by d (dr ), such that d =L/2 − a + x and dr = L/2 − a − x. We define d0 = L/2 − aas the distance between the colloid and the walls when thecolloid is placed on the midplane (x = 0). We also assume
that the walls have a fixed charge density, σw , which can beeither positive or negative.
The electrostatic interaction between a spherical colloid
and a flat wall, embedded in an ionic solution, is described by
the Poisson-Boltzmann theory [22,23]. Assuming that the colloid and wall are not highly charged [see the exact conditionin Appendix A after Eq. (A2)], the theory can be linearized
and reduces to the Debye-Hückel (DH) theory. In the DH
theory, the electrostatic potential, ψ, at each point in spacer = (x, y, z) is given by the linearized equation
 2

∇ − κD2 ψ (r) = 0,(1)
where κD is the inverse screening length, κD =
2
1/2
λ−1
, λD is the Debye screening length,
D = [2e n/(ε0 εkB T )]e is the elementary charge, ε0 is the vacuum permittivity, andwhere the osmotic pressure, [23], is related to ψ by
εε0 −ψ 2 + κD2 ψ 2 ,=(3)
2
and ψ  is the first derivative along the x direction (note that
U (x) is defined up to a constant).
Equations (2) and (3) (see details in Appendix A) yield thefollowing interaction potential:
cosh(κD d0 ) + cosh(κD x)
πa
2σ
σ
ln
U (x) =
c w
cosh(κD d0 ) − cosh(κD x)
εε0 κD2 2cosh(2κD d0 )2.
+ σc + σw ln
cosh(2κD d0 ) − cosh(2κD x)(4)
Due to the symmetry, U  (x) = 0 at x = 0. However, from
Eq. (4) it follows that this equilibrium is stable (minimum ofU ) only if
σw
σw< −eκD d0 or> −e−κD d0 .
(5)
σc
σc
This is exactly the condition that the walls repel the colloid
placed at the midplane position x = 0, instead of attracting
it [24]. In what follows, we focus on the case where thecolloid and walls have the same charge sign, and the stability
condition, Eq. (5), is thus always satisfied.
Expanding U (x) to second order around the midplane,x = 0, we obtain the harmonic potential approximation
U ≈ 21 Kx 2 ,(6)where K = U  |x=0 is the effective spring constant. From
Eq. (4), this approximation leads to
2π a 2σc σw cosh(κD d0 ) + σc2 + σw2
K=
.
(7)
ε0 ε
sinh2 (κD d0 )The harmonic interaction potential in Eq. (6) will be usedas the potential energy in the following analysis. The approximation underestimates the strength of the interaction potential

042607-2

BROWNIAN MOTION OF A CHARGED COLLOID IN …

PHYSICAL REVIEW E 103, 042607 (2021)

with λ being∞

λ(ζ ) =

FIG. 2. Contour plot of the normalized effective spring constant,K̃ ≡ K/(π aσc2 /ε0 ε) where K = U  (x)|x=0 [see Eqs. (6) and (7)] as afunction of κD d0 and σw /σc . The color code is associated with ln K̃.
Large K̃ is obtained in the high σw /σc and small κD d0 limits. It decays
to zero as κD d0 is increased.

close to the walls, but this only has a minor effect on the
dynamics, as will be discussed in depth in Sec. V .
In Fig. 2, K̃ ≡ K/(π aσc2 /ε0 ε) is plotted for different κD d0
and σw /σc > 0 values. For fixed κD d0 and σc , K diverges as
σw → ∞, and monotonically decreases as σw decreases. For
given σc and σw , K decreases rapidly when κD d0 is increased.
In Appendix B we derive the corresponding K for fixedsurface potential on the walls, rather than fixed surface charge.

n(n + 1)
4sinh ζ
3(2n − 1)(2n + 3)n=1
2 sinh [(2n + 1)ζ ] + (2n + 1) sinh(2ζ )−1.
×
4 sinh2 [(n + 1/2)ζ ] − (2n + 1)2 sinh2 ζ(9)Here ζ = cosh−1 (1 + d/a), and F⊥ and v⊥ are the force andvelocity in the perpendicular direction, respectively. In thed  a limit, we have λ = 1 and the regular Stokes’ law isrecovered. In the da limit, which is our main interest,
Eq. (9) becomes, to leading order, λ ≈ a/d.
When the colloid is placed in between two walls, the modified Stokes’ law has a much more complicated expression,
relying on extensive numerical computations [8]. However,
several approximations have been proposed [11,27,28]. The
simplest one is the linear superposition approximation thatwas shown to agree fairly well with experiments [29],F⊥ ≈ 6π ηa(λ + λr − 1)v⊥ ,(10)where λ and λr are calculated from Eq. (9), with d = dand d = dr , respectively. Note that the negative third term in
Eq. (10) guarantees that we recover F⊥ = 6π aηv⊥ when thewalls are very far apart.
In the narrow confinement limit, d , dra, the perpendicular diffusion coefficient D⊥ can be approximated as


x2
d dr
D⊥ (x) ≈ D∞
(11)
= D0 1 − 2 ,a(d + dr )
d0
where

III. CONFINED BROWNIAN DYNAMICS

In a free medium (no boundaries), the colloid performs a
simple Brownian motion, described by a diffusion coefficient
D. By the Einstein relation, we have D = μkB T , where μ
is the mobility, defined by the ratio of the colloid terminal
drift velocity to an applied force, μ = v/F . If the colloid isneutral, and at low Reynolds number, the mobility follows the
Stokes’ law, μ = 1/(6π ηa), and the diffusion coefficient of
this motion is denoted as D∞ = kB T /(6π ηa).
The mobility of a charged colloid in an ionic solution
is reduced due to the drag of the surrounding ionic cloud.
However, for low surface electrostatic potential or low surface
charge, the reduction is only of a few percent [25,26], andis therefore neglected here. Note that the motion of ions, byvirtue of their small size, is much faster than the colloid,and is assumed to be in equilibrium throughout the colloid
motion. Their only effect on the colloid dynamics (apart fromthe reduction of its mobility, which we neglect) is through theequilibrium interaction potential, Eq. (6).
The presence of walls in the vicinity of the colloid, however, can modify the colloid mobility substantially due tohydrodynamic effects. As shown by Brenner [6], the Stokes’
law of a motion in the perpendicular direction to a single solid
wall, at distance d, is modified in the following way:F⊥ = 6π ηaλ(ζ )v⊥ ,

D∞ d0
(12)
2a
and the relations d = d0 + x and dr = d0 − x were used. In
Fig. 3, we compare between the approximated diffusion coefficient given by Eq. (11), and the “exact” numerical resultsof the diffusion coefficient of a colloid between two walls, for
motion in the perpendicular direction. The numerical data is
adapted from Fig. 6 of Ref. [8]. Here we show D⊥ /D∞ , whichis the inverse of λ, D⊥ /D∞ = 1/λ, whereas in Ref. [8], λ wascalculated [see its definition in Eq. (8)].
As seen in Fig. 3, for L/(2a) = 1.5, the approximatedexpression overestimates the “exact” one by ∼20% at the
midplane (x = 0), and becomes more accurate as we lookaway from the midplane. For L/(2a) = 1.25, the deviationin the midplane is of ∼10%, and decreasing away from the
midplane. For L/(2a) = 1.1, there is only one data point available at the midplane, and it deviates from the approximation
by only ∼2%. This strongly supports that in the d0
a limit(L ≈ 2a), the approximated Eq. (11) can by used.
The parallel motion is also described by a positiondependent diffusion coefficient D (x) although we do not
calculate it here [9].
D0 =A. Position-dependent diffusion coefficient (PDDC)(8)

B. Fokker-Planck equation

The dynamics of a colloid with diffusion coefficientsD⊥ (x) and D (x) in the perpendicular and parallel directions,042607-3

AVNI, KOMURA, AND ANDELMANPHYSICAL REVIEW E 103, 042607 (2021)

walls, where the current is defined by the continuity equation ∂ p/∂t = −∂ j/∂x. From Eq. (14), the current becomes
j = −D⊥ (x)[(Kx/kB T )p + ∂ p/∂x]. Notice that the conditions j = 0 at the walls are automatically satisfied by the factthat D⊥ vanishes there.
IV. MEAN-SQUARE DISPLACEMENT

Equation (15), together with the definitions of K and D0
in Eqs. (7) and (12), respectively, is the principal equationof this paper. By solving it, one can derive the MSD, which
can be measured in experiments. The MSD of an ensemble ofcolloids in equilibrium is
 d0
 d0
(x(t ) − x0 )2  =dx0dx (x(t ) − x0 )2 p(x, x0 ; t ),
−d0

−d0

(16)
FIG. 3. The diffusion coefficient D⊥ of a colloid between two flat
surfaces, normalized by its diffusion coefficient in a free medium,
D∞ = kB T /(6π ηa), as a function of x/d0 . Full black lines are
numerical results obtained in Ref. [8], and red dashed lines show
the approximated expression in Eq. (11), which relies on the linear
a limit. For L/(2a) =superposition approximation in the d , dr
1.1, a single data point was available in the numerical calculation:
D⊥ /D∞ = 0.049 at x = 0. The approximated expression in Eq. (11)yields D⊥ /D∞ (x = 0) = 0.05 at x = 0. The difference between thetwo values is smaller than the plot resolution.

respectively, and under potential field U (x), is governed by
the generalized Fokker-Planck equation [30],∂P(r, t )
∂
∂
1 ∂U
=
D⊥ (x)+
P(r, t )∂t∂x
kB T ∂x
∂x

 2
∂
∂2
P(r, t ),(13)+ D (x)+∂y2∂z2where P(r, t ) is the time-dependent probability distributionfunction (PDF) of the colloid position. The full threedimensional motion of the colloid is quite complex, because
the motion in the y-z plane is coupled to the motion along thex axis through the diffusion coefficient D (x). However, ourfocus is on the perpendicular motion in the x direction, whichby symmetry, does not depend on the position in the y-z plane.
Then one can integrate the equation both over y and z from−∞ to ∞, and obtain




∂
∂
1 ∂U
∂ p(x, t )=D⊥ (x)+
p(x, t ) ,
(14)
∂t
∂x
kB T ∂x∂x
where p(x, t ) is the reduced probability distribution function
p(x, t ) = dy dz P(r, t ).
Substituting the harmonic interaction potential obtained in
Eq. (6) and the PDDC of Eq. (11), we obtain



∂p
Kx
∂
x2
∂p
= D0
1− 2
p+
.
(15)
∂t
∂x
kB T
∂x
d0
All the electrostatic effects are captured by the term proportional to K. We note that, in general, rigid walls
impose zero current boundary conditions, i.e., j = 0 on thewhere p(x, x0 ; t ) is the probability of a colloid to beat x0 at time t = 0 and at x at time t. As the initial position, x0 , is drawn from an equilibrium distribution, we can write p(x, x0 ; t ) in terms of the conditional probability, p(x, x0 ; t ) = p(x, t|x0 , 0)peq (x0 ), wherepeq (x0 ) ∝ exp(−Kx02 /2kB T ). It then follows that the MSD is(x(t ) − x0 )2  =d0−Kx02 /2kB T(x(t ) − x0 )2 x0−d0 dx0 e.
d0−Kx02 /2kB T
dx
e
0
−d0

(17)

where Ox0 is the average over an ensemble of colloids with the same initial position, x0 , at t = 0, i.e.,Ox0 = dx p(x, t|x0 , 0)O, with p(x, t|x0 , 0) being a solution of Eq. (15) with the initial condition p(x, 0) = δ(x − x0 ).
The MSD resulting from Eq. (15) will be analyzed in depth inwhat follows.
A. Vanishing interaction potential limit

While Eq. (15) does not have an analytic solution, it can besolved for certain limits. In the vanishing interaction potential
limit, K is omitted altogether and we are left with

2x ∂ px2 ∂ 2 p∂p.
(18)= D0 − 2+ 1− 2∂t
d0 ∂x
d0 ∂x 2
The above equation, for a colloid at x = x0 at time t = 0, has
the solution [30]∞p(x, t|x0 , 0) =1
2
(2n + 1)e−D0 n(n+1)t/d0
2d0 n=0
× Pn

x0
x
Pn
,
d0
d0

(19)

where Pn (x) is the Legendre polynomial of order n. The PDF
of Eq. (19) is characterized by the following first and second
moments,x(t )x0 = x0 e−2D0 t/d0 ,
2

d02 
2
1 − e−6D0 t/d0 ,(20)
3
In the limit of t → ∞, the PDF reduces to a uniform distribution between the two walls that are positioned at x = ±d0 .

042607-4

x 2 (t )x0 = x02 e−6D0 t/d0 +2

BROWNIAN MOTION OF A CHARGED COLLOID IN …

PHYSICAL REVIEW E 103, 042607 (2021)

In the vanishing interaction-potential limit, K = 0, Eq. (17)simply becomes
 d0

1
2
dx0 x 2 (t )x0 − 2x0 x(t )x0 + x02 ,(x(t ) − x0 )  =2d0 −d0(21)
leading to the MSD(x(t ) − x0 )2  =2d02 
2
1 − e−2D0 t/d0 .
3(22)

The above MSD approaches 2d02 /3 in the limit of t → ∞.
B. Strong interaction-potential limit

In the opposite strong interaction-potential limit, the confinement due to electrostatic effects is strong, and the colloid
remains close to the midplane. Consequently, its diffusion
coefficient does not change, and equals to D0 . Then the termthat depends on d0 in Eq. (15) can be omitted, and we obtain




K
∂p
∂p
∂2 p
(23)
= D0
p+x
+ 2 .
∂t
kB T
∂x
∂x
Note that this is equivalent to the case of permeable walls,
where the walls interact electrostatically but not hydrodynamically. Additionally, note that a strong interaction potential U
does not contradict the assumption employed in Sec. II of asmall electrostatic potential ψ. This is evident, for example,
by the dependence of U on the colloid size [see Eq. (4)], which
is absent in ψ [see Eq. (A2)].
The above equation can be solved for a colloid at x = x0 attime t = 0, yielding [31]
1/2K
p(x, t|x0 , 0) =2π kB T (1 − e−2D0 Kt/kB T )K (x − x0 e−D0 Kt/kB T )2
. (24)
× exp −
2kB T (1 − e−2D0 Kt/kB T )FIG. 4. The apparent diffusion coefficient Dapp [see Eq. (30)],
scaled by D0 , as a function of α = Kd02 /kB T .
C. Short-time diffusion coefficient

For a known initial position, p(x, 0) is infinitely sharp(Dirac delta function). As a result, the term proportional to
the second derivative of p in Eq. (15) dominates the dynamics
for short times,∂p∂2 p= Deff (x0 ) 2 ,∂t∂xwith Deff (x0 ) ≡ D0 (1 − x02 /d02 ). The colloid experiences, fora short while, a free diffusion with diffusion coefficient
Deff (x), with the following conditional probability distributionfunction:
1/2



1
(x − x0 )2. (29)
p(x, t|x0 , 0) =exp −4π Deff (x0 )t
4Deff (x0 )tThis results in a linear time dependence of MSD,
(x − x0 )2  = 2Dappt, where Dapp is the apparent diffusion
coefficient
Dapp =For the averages, we can extend the integral range inOx0 = dx p(x, t|x0 , 0)O from ±d0 to ±∞. Then the firstand second moments are

kB T(1 − e−2D0 Kt/kB T ).
K

(25)

1
= D0 1 − +αK
2π kB Tα=1/2  ∞−∞

dx0 e−Kx0 /2kB T
2


× x (t )x0 − 2x0 x(t )x0 + x02 ,2(26)

yielding(x(t ) − x0 )2  =2kB T(1 − e−D0 Kt/kB T ),K

(27)which asymptotically approaches 2kB T /K for t → ∞. Equations (19) and (24) and their corresponding MSD functions,Eqs. (22) and (27), represent the two limiting dynamical
behaviors.

√
2/(π α)e−α/2,√erf( α/2)(30)where α is a dimensionless parameter defined by

For the MSD, Eq. (17) now reads(x(t ) − x0 )2  =d0−Kx02 /2kB T
Deff (x0 )−d0 dx0 e2
d0−Kx0 /2kB T−d0 dx0 e



x(t )x0 = x0 e−D0 Kt/kB T ,x 2 (t )x0 = x02 e−2D0 Kt/kB T +

(28)

Kd02
kB T

(31)

√
2
x
and erf(x) = (2/ π ) 0 e−z dz is the error function. Notice√is the ratio between two characteristic length scales,that α √d0 and kB T /K. While d0 is√ the length over which the
diffusion coefficient changes, kB T /K is the characteristic
interaction potential length scale, which is small for strong
potentials and large for weak ones.
In Fig. 4, Dapp /D0 is plotted as a function of α. Upon
increasing α, Dapp monotonically grows from Dapp /D0 = 2/3
obtained for α = 0 to Dapp /D0 = 1 valid for the α → ∞ limit.
We conclude that the colloid moves faster as the electrostatic
effects are stronger (large α). The change in Dapp occurs when
α is of order unity.

042607-5

AVNI, KOMURA, AND ANDELMANPHYSICAL REVIEW E 103, 042607 (2021)

FIG. 5. The mean-square displacement (MSD) of the colloid, (x − x0 )2  scaled by d02 , as a function of t scaled by d02 /D0 for (a) α = 0.1,(b) α = 10, and (c) and (d) α = 2. Red dots are the numerical solution of Eq. (15). Solid blue line is the vanishing interaction-potential limit[see Eq. (22)]. Dashed black line is the strong interaction-potential limit [see Eq. (27)]. In (d), the short-time limit of α = 2 is shown, together
with the prediction, (x − x0 )2  = 2Dappt with Dapp defined in Eq. (30) (thick green line).

the Euler method,

D. Numerical analysis

The full dynamical equation, Eq. (15), can be made dimensionless by using the dimensionless position x̃ = x/d0 andtime t˜ = (D0 /d02 )t aspi (t + t )

(32)





t 
pi+1 (t ) − pi (t )21−xi+1
αxi+1 pi+1 (t )+
x
x




pi (t ) − pi−1 (t )− 1 − xi2 αxi pi (t ) +
.
(33)
x

where α is given by Eq. (31). The vanishing interactionpotential limit, Eq. (18), and strong interaction-potential limit,
Eq. (23), correspond to α = 0 and to the limit α → ∞, respectively [in order to obtain Eq. (23) from Eq. (32), one needs
to redefine the normalized position and time].
We calculate the MSD of the particle for different α values
by solving Eq. (32) numerically. In the following we use thesame normalized variables as in Eq. (32), but drop the tilde
signs for brevity. We divide our space into N lattice points{x1 , . . . , xN }, where x1 = −1 and xN = 1, and define pi (t ) ≡p(xi , t ). For a given x0 , we start with a distribution that approximates a Dirac delta function, pi (0) = (2π δ)−1/2 exp[(xi −x0 )2 /2δ], with δ
1. We then iterate each time step using

The zero current boundary condition is guaranteed due tothe cancellation of the diffusion coefficient, proportional to(1 − xi2 ), at x1 and xN . We used a lattice spacing of x = 0.02,
and a time step t between 6.25 × 10−5 and 2 × 10−4 , depending on α. For each initial condition, x0 ∈ {−1, 1}, we
calculated (x(t ) − x0 )2 x0 and obtained the MSD throughEq. (17).
The numerically obtained full MSD is shown in Fig. 5 fordifferent α values, and is compared with the limits in Eqs. (22)and (27). Figure 5(a) shows that for α = 0.1, the MSD coincides, almost perfectly, with the vanishing interaction potential case (solid blue line), while in Fig. 5(b) for α = 10, it coincides with the strong interaction-potential case (dashed blackline). However, for an intermediate value, α = 2 in Fig. 5(c),the numerically obtained MSD deviates substantially from the



∂p
∂
∂p
=
(1 − x̃ 2 ) α x̃ p +
,
∂ t˜
∂ x̃
∂ x̃

= pi (t )+042607-6

BROWNIAN MOTION OF A CHARGED COLLOID IN …

PHYSICAL REVIEW E 103, 042607 (2021)

FIG. 6. The diffusion characteristic time τD , defined as the timeat which the MSD reaches half of its value in the t → ∞ limit, normalized by d02 /D0 , as a function of α. Solid blue line is the vanishinginteraction-potential limit, τD = d02 ln 2/(2D0 ), and dashed black lineis the strong interaction-potential walls limit, τD = d02 ln 2/(D0 α).
Red dots are the values calculated numerically from Eq. (15) andinterpolate between the two limits.

two limits. In the short-time limit, it lies between the two limiting cases, while in the long-time limit it has a significantly
lower value than both limits. In the short-time region shownin Fig. 5(d), the numerical calculation coincides with Eq. (30)(thick green line), and the MSD grows linearly with a slopethat increases with α. At saturation, the MSD for α = 0.1
equals ∼0.6d02 , for α = 2 equals ∼0.4d02 , and for α = 10, it
approaches ∼0.2d02 , indicating that the colloid motion is more
localized when the electrostatic force is strong.
A distinct difference between the MSD for different α
values is the characteristic time at which MSD saturates.
This is the time it takes the colloid to explore its available
space. Denoting τD as the time at which the MSD reaches
half of its maximal value (the diffusion characteristic time),
we obtain that in the α = 0 case, τD = d02 ln 2/(2D0 ) [seeEq. (22)], while in the α → ∞ limit, τD = kB T ln 2/(D0 K ) =d02 ln 2/(D0 α) [see Eq. (27)]. In Fig. 6, τD scaled by d02 /D0
is shown as a function of α. Similar to Dapp , the diffusion
characteristic time τD interpolates between the two limits
mentioned above, where the crossover occurs around α ∼ 1.
The decrease of τD as a function of α, seen in Fig. 6, isattributed to two effects. First, as α increases, the colloid
is more localized due to electrostatic forces and explores asmaller region, and second, its average diffusion coefficient,being further away from the walls, becomes larger.
E. Dependence of α on electrostatic properties

So far we have seen that the dynamics are determined by
the value of α. Substituting the value of K derived in Eq. (7),
we obtain
2π ad02 2σc σw cosh(κD d0 ) + σc2 + σw2
.
(34)
α=
ε0 εkB T
sinh2 (κD d0 )In an experimental setup, α can be controlled by changing
the ionic strength (salt concentration,n), and consequently the√
Debye length, λD ∼ 1/ n. The dependence of α on λD = κD−1
is plotted in Fig. 7, for reasonable colloid and solvent param-

FIG. 7. The natural logarithm of the dimensionless parameter
α = Kd02 /kB T , where K is given by Eq. (7), as a function of thescreening length, λD , on a logarithmic scale. The system parameters
are L = 5 μm, σc = 0.003e/nm2 , T = 300 K, and ε = 80. Blue lines
are for d0 = 500 nm (a = 2.0 μm), and red lines are for d0 = 100 nm(a = 2.4 μm). Solid lines are for uncharged walls, σw = 0, dashedlines are for walls with surface charge density, σw = σc , and dotdashed line are for σw = 10σc .

eters. As expected, α is large for large screening length anddecreases towards zero as λD decreases.
As we have shown, α = 1 signifies the crossover from a
PCCD dominated motion to an electrostatically dominatedone. For the range of parameters in Fig. 7, this crossover
occurs roughly when λD ∼ d0 /10. The screening length that
corresponds to the crossover increases with d0 , and for a fixedσc (σw ), it decreases as σw (σc ) increases.
V. SUMMARY AND DISCUSSION

We studied the dynamics of a charged colloid underrestricted confinement of two charged surfaces. The combination of electrostatic and hydrodynamic forces exerted by
the walls result in a unique behavior. This behavior can bequantified in terms of a dimensionless parameter α in Eq. (31)that determines the interplay between the electrostatic interaction and the position-dependent drag force. The parameter α
is also given by Eq. (34), showing how the colloid motion,
including its short-time behavior, Eq. (30), and long-time
behavior (Figs. 5 and 6), depends on the geometry and electrostatic properties. In particular, α can be tuned by changing
the screening length (Fig. 7) that is usually an easily controlled parameter in experiments. At small screening length(small α), the colloid moves slower than at larger screening
length (large α). In addition, as λD is increased the colloid
explores a smaller region in space throughout its motion. The
two effects lead to a decrease in the diffusion characteristic
time, τD , as λD is increased.
We note that the harmonic approximation of the interaction potential allows us to obtain analytical results. If we
use the full interaction potential, Eq. (4), in the dynamical
equation, we could not have expressed it in terms of a single
parameter, α, as in Eq. (32). This harmonic approximation
underestimates the strength of the interaction potential near
the walls. However, as long as d0 is of the order of the screen-

042607-7

AVNI, KOMURA, AND ANDELMANPHYSICAL REVIEW E 103, 042607 (2021)

ing length, or smaller, the approximation is valid except for the
regime very close to the wall where the interaction potential U
diverges. This small region near the walls does not affect the
MSD that integrates the motion throughout the entire space.
Moreover, the divergence of U close to the wall is unphysicalbecause at such close proximity the DH theory, Eq. (1), is nolonger valid. For d0  λD , the following limiting expression
should be used instead:
2Qσw −κD d0
U (x) ≈ecosh(κD x).
(35)
ε0 εκD2 a
However, as in this limit the interaction potential is very small
and can be ignored, there is no substantial difference in the
MSD between the exact and approximated potentials.
We also note that, in our analysis, the effect of the ioniccloud on the diffusion coefficient of the charged colloid wasneglected. In the free medium, the effect of the ionic cloud
was calculated and turned out to be very small (of only a few
percent) [25,26]. Since in a restricting confinement this effectwas not fully investigated [18,32], it is not considered here forsimplicity.
Although the above-mentioned calculations are approximated, we expect our results to qualitatively describe theBrownian motion of a charged colloid between confining
charged walls. This should be tested in future experiments thatwill focus on very restricted confinement, as considered here.
ACKNOWLEDGMENTS

We would like to thank R. Adar for useful discussions.
Y.A. is thankful for the support of the Clore Scholars Programme of the Clore Israel Foundation and the hospitality atTokyo Metropolitan University, where part of this work has
been conducted. S.K. acknowledges support by a Grant-inAid for Scientific Research (C) (Grants No. 18K03567 and
No. 19K03765) from the Japan Society for the Promotion of
Science, and support by a Grant-in-Aid for Scientific Researchon Innovative Areas “Information Physics of Living Matters” (Grant No. 20H05538) from the Ministry of Education,
Culture, Sports, Science and Technology of Japan. D.A. acknowledges support from the Israel Science Foundation (ISF)
under Grant No. 213/19.
APPENDIX A: DERIVATION OF EQ. (4)In the large sphere limit, the force acting on the colloid
from the right wall is independent of that acting from the left
wall, and vice versa. For each side, the force on the colloid isapproximated by the Derjaguin approximation [22] ∞dl (l ),(A1)
F (d ) ≈ 2π adwhere is the force per unit area exerted on the colloid, if thecolloid is considered as flat rather than spherical, and d is the
distance between the colloid and the wall (d = d for the left
wall and d = dr for the right wall).
The electrostatic potential in the Derjaguin approximation
can be calculated in a straightforward way. As it was performed multiple times for other geometries and boundary

conditions [22], we describe it only briefly here. Since thetwo sides of the colloid are decoupled, we start by solving the DH equation in Eq. (1) for a single wall, whileassuming that the colloid is flat (a → ∞). Denoting the distance from the left wall, without loss of generality, by ξ ,the boundary conditions are dψ/dξ |ξ =0 = −σw /(ε0 ε) anddψ/dξ |ξ =d = σc /(ε0 ε). Then the solution isψ (ξ ; d ) =σw cosh [κD (d − ξ )] + σc cosh(κD ξ ).
εε0 κD sinh (κD d )(A2)

The DH approximation is valid when the electrostatic
potential on the wall (ξ = 0) and on the colloid surface(ξ = d) is small compared to kB T /e [23]. This is satisfied
when σc , σwe/(B λD ) and σc , σwed/(lB λ2D ), where2B = e /(4π ε0 εkB T ) is the Bjerrum length. When the colloid gets very close to the wall, d → 0, the above Derjaguin
approximation fails. However, as discussed in Sec. V, thisregime does not affect our results.
One can show that the force per unit area, (d ), is independent of ξ and equals to [23,24]εε0(d ) =−ψ 2 (ξ ; d ) + κD2 ψ 2 (ξ ; d )22σc σw cosh (κD d ) + σc2 + σw2
=−
,
(A3)
2εε0 sinh2 (κD d )where ψ  = dψ/dξ . Substituting Eq. (A3) into Eq. (A1), we
obtain the force between the colloid and one surface. Integrating the force, we further get the interaction potential U1dwith a single wall, U1 (d ) = − δ dh F (h) where the lower
cutoff δ is an arbitrary distance. The total interaction potential of the colloid, when taking into account the two walls,
is U (x) = U1 (d0 + x) + U1 (d0 − x). After some algebra, we
obtainπacosh(κD d0 ) + cosh(κD x)
2σc σw ln
U (x) =
cosh(κD d0 ) − cosh(κD x)
εε0 κD2 2cosh(2κD d0 ).
+ σc + σw2 lncosh(2κD d0 ) − cosh(2κD x)(A4)
We note that since U (x) is constructed from integrating
twice, the effective spring constant K = U  |x=0 [see Eq. (6)]can be directly obtained from ,K = 4π a (d0 ).

(A5)

APPENDIX B: CHARGED WALLS WITH FIXED SURFACEELECTROSTATIC POTENTIAL

Assuming that the walls have a fixed surface electrostatic potential, V , we repeat the calculation of Appendix A
and Sec. II.
For the space between the left wall and the colloid,
the boundary conditions are ψ|ξ =0 = V and dψ/dξ |ξ =d =σc /(ε0 ε). The solution to Eq. (1) is042607-8

ψ (ξ ) =V ε0 εκD cosh [κD (d − ξ )] + σc sinh (κD ξ ),εε0 κD cosh (κD d )(B1)

BROWNIAN MOTION OF A CHARGED COLLOID IN …

PHYSICAL REVIEW E 103, 042607 (2021)

and following a similar calculation as done in Appendix A, the total interaction potential isπa(V ε0 εκD )2 − σc2 ln (cosh[κD (d0 + x)] cosh[κD (d0 − x)])
U (x) =
2εε0 κDκD (d0 + x)
κD (d0 − x)
−1
−1
tanh
tanh
+ tan
.
− 4σcV εε0 κD tan
2
2

Qualitatively, the behavior of U when varying V is similar toits behavior in Sec. II. The stability condition isδ > e−κD d0orδ < −eκD d0 ,(B3)where δ = V ε0 εκD /σc . Unlike Sec. II, here the stability condition does not coincide with the condition that the walls repel
the colloid at the midplane. The latter condition isδ>eκD d0 +1√
1 + e2κD d0orδ<eκD d0 −1.
√
1 + e2κD d0(B4)[1] P. S. Burada, P. Hänggi, F. Marchesoni, G. Schmid, and P.
Talkner, Chem. Phys. Chem. 10, 45 (2009).
[2] R. Phillips, J. Kondev, J. Theriot, and H. Garcia, Physical
Biology of the Cell (Garland Science, London, 2012).
[3] T. M. Squires and S. R. Quake, Rev. Mod. Phys. 77, 977 (2005).
[4] M. Krishnan, N. Mojarad, P. Kukura, and V. Sandoghdar,
Nature (London) 467, 692 (2010).
[5] F. Ruggeri and M. Krishnan, Phys. Rev. E 96, 062406 (2017).
[6] H. Brenner, Chem. Eng. Sci. 16, 242 (1961).
[7] A. J. Goldman, R. G. Cox, and H. Brenner, Chem. Eng. Sci. 22,
637 (1967).
[8] P. Ganatos, S. Weinbaum, and R. Pfeffer, J. Fluid Mech. 99, 739(1980).
[9] P. Ganatos, R. Pfeffer, and S. Weinbaum, J. Fluid Mech. 99, 755(1980).
[10] K. H. Lan, N. Ostrowsky, and D. Sornette, Phys. Rev. Lett. 57,
17 (1986).
[11] L. Lobry and N. Ostrowsky, Phys. Rev. B 53, 12050 (1996).
[12] M. Hosoda, K. Sakai, and K. Takagi, Phys. Rev. E 58, 6275
(1998).
[13] M. A. Bevan and D. C. Prieve, J. Chem. Phys. 113, 1228
(2000).
[14] P. Lançon, G. Batrouni, L. Lobry, and N. Ostrowsky, Physica A
304, 65 (2002).
[15] E. R. Dufresne, D. Altman, and D. G. Grier, Europhys. Lett. 53,
264 (2001).
[16] P. Huang and K. S. Breuer, Phys. Rev. E 76, 046307(2007).
[17] G. M. Kepler and S. Fraden, Langmuir 10, 2501 (1994).

(B2)

If the stability condition of Eq. (B3) is satisfied, we canapproximate the interaction potential around x = 0, to be U ≈Kx 2 /2, with the effective spring constant
2π a 2σcV εε0 κD sinh(κD d0 ) + (V εε0 κD )2 − σc2.
K=εε0
cosh2 (κD d0 )(B5)
As shown in Sec. IV D, the value of K determines the type of
colloid diffusion. For a fixed σc , K grows as V is increased,
and decays to zero for κD d0  1.

[18] S. L. Eichmann, S. G. Anekal, and M. A. Bevan, Langmuir 24,
714 (2008).
[19] S. Fringes, F. Holzner, and A. W. Knoll, Beilstein J.
Nanotechnol. 9, 301 (2018).
[20] M. I. M. Feitosa and O. N. Mesquita, Phys. Rev. A 44, 6677
(1991).
[21] N. A. Frej and D. C. Prieve, J. Chem. Phys. 98, 7552 (1993).
[22] J. Israelachvili, Intermolecular and Surface Forces, 3rd ed.
(Academic, New York, 2011).
[23] T. Markovich, D. Andelman, and R. Podgornik, in Handbookof Lipid Membranes: Molecular, Functional, and Materials
Aspects, edited by C. Safinya and J. Raedler (CRC Press, Taylor& Francis Group, 2021), Chap. 6.
[24] V. A. Parsegian and D. Gingell, Biophys. J. 12, 1192 (1972).
[25] H. Ohshima, T. W. Healy, L. R. White, and R. W. O’Brien, J.
Chem. Soc., Faraday Trans. 2 80, 1299 (1984).
[26] G. A. Schumacher and T. G. M. van de Ven, Faraday Discuss.
83, 75 (1987).
[27] L. P. Faucheux and A. J. Libchaber, Phys. Rev. E 49, 5158
(1994).
[28] T. Benesch, S. Yiacoumi, and C. Tsouris, Phys. Rev. E 68,
021401 (2003).
[29] B. Lin, J. Yu, and S. A. Rice, Phys. Rev. E 62, 3909 (2000).
[30] A. W. C. Lau and T. C. Lubensky, Phys. Rev. E 76, 011123(2007).
[31] H. Risken, The Fokker-Planck Equation (Springer, Berlin,
Heidelberg, 1996).
[32] B. Chun and A. J. C. Ladd, J. Colloid Interface Sci. 274, 687(2004).

042607-9