Power-law electrokinetic behavior as a direct probe ofeffective surface viscosity
Yuki Uematsua,b,∗, Roland R. Netzb , Douwe Jan BonthuisbaDepartment of Chemistry, Kyushu University, Fukuoka 819-0395, Japanb
Fachbereich Physik, Freie Universität Berlin, 14195 Berlin, Germany

Abstract
An exact solution to the Poisson-Boltzmann and Stokes equations is derived
to describe the electric double layer with inhomogeneous dielectric and viscosity profiles in a lateral electric field. In the limit of strongly charged surfaces
and low salinity, the electrokinetic flow magnitude follows a power law as afunction of surface charge density. Remarkably, its exponent is determinedby the interfacial dielectric constant and viscosity, the latter of which has
eluded experimental determination. This provides a novel method to extractthe effective interfacial viscosity from standard electrokinetic experiments.
We find good agreement between our theory and experimental data.
Keywords: Electrokinetics, interfacial slip, liquid-solid interface, electric
double layer
Introduction. — Electro-osmosis is the motion of an electrolyte solution induced by an electric field along a surface [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],
whereas electrophoresis is the electric-field-induced motion of a colloid suspended in an electrolyte [1, 2, 3, 13, 14, 15]. These electrokinetic effects play
important roles in industry and bio-systems [16, 17, 18]. Originating in thenanometer-wide electric double layer, both effects depend sensitively on the
structure and the dynamic properties of this interfacial layer. Simulation
studies demonstrate that the surface layer exhibits an increased viscosity
at hydrophilic surfaces [19], in agreement with early theories used to explain electrokinetic experiments [3]. At hydrophobic surfaces, however, a∗Corresponding author
Email address: uematsu@chem.kyushu-univ.jp (Yuki Uematsu)

Preprint submitted to Chemical Physics Letters

December 19, 2016

depletion layer and a finite slip length are observed, indicating a decreased
effective viscosity [4, 5, 6, 7, 8, 9, 19, 20, 21]. The interfacial dielectric profile of water has been studied using molecular dynamics (MD) simulations[4, 5, 6, 7, 8, 19, 20, 21, 22, 23, 24, 25, 26], revealing a highly inhomogeneous and oscillating profile near the interface [22, 23]. Whereas an effective
interfacial dielectric constant can be readily obtained from capacitance measurements, experimental determination of the interfacial viscosity has been
controversial [27, 28, 29].
The traditional model describing electro-osmosis and electrophoresis isbased on the Poisson-Boltzmann (PB) and Stokes equations [1, 2, 3]. However, this model fails to capture many experimental observations, such as
the saturation of the electrokinetic flow with rising surface charge density[3, 30, 31, 32]. The discrepancies have largely been resolved by computing
numerical solutions to the PB and Stokes equations, including the dielectricand viscosity profiles extracted from simulations [24, 25]. Yet the numericalsolutions cannot be used to unequivocally extract the interfacial profiles fromexperimental data without prior knowledge.
In this Letter, we construct an exact solution and derive analytical asymptotic expressions for the electrokinetic flow beyond the linear theory [33].
With these expressions, information about the subnanometer-scale interfacial
properties can be directly acquired from standard electrokinetic experiments.
In particular, in the limit of strongly charged surfaces and low salinity, we obtain an analytical power-law relation for the magnitude of the electrokinetic
flow.
Central result. — We model the interfacial structure by a layer of width z ∗ ,which has a dielectric constant εint and viscosity ηint different from the bulk
values ε and η. In addition, we include specific adsorption of monovalent ions
into the interfacial layer, expressed by a constant chemical potential kB T αwith kB T being the thermal energy. Using this box model, we derive fromthe PB and Stokes equations (explained below) that in the limit of strongly
charged surfaces and low salinity, the electrokinetic surface charge density
σek scales with the bare surface charge density σ0 as
!1−γ
p
r α
e ε 2εint ε0 kB T arccosh eα ε/εintp|σ0 |γ ,(1)|σek | ≈∗αεintez e ε/εint − 1with ε0 being the electric permittivity of vacuum and e the elementary charge.
Notably, Eq. (1) exhibits a power-law behavior with exponent γ ≡ εint η/εηint .
2

-̀ek [e/nm2]10.2

̀=0.0740.1

0.1
̀=0.015

0.01
0.01

(a)

0.1
-̀0 [e/nm2]

full (hydrophilic, MD)
full (hydrophobic, MD)
full (hydrophilic, Fit)
full (hydrophobic, Fit)

1

0.05
0.08

(b)

0.8
2

-̀0 [e/nm ]experiment (hydrophilic)experiment (hydrophobic)
power law (hydrophilic, Fit)
power law (hydrophobic, Fit)

Figure 1: The electrokinetic surface charge density σek for a c0 = 1 mM electrolyte solution
as a function of the bare surface charge density σ0 (double-logarithmic). The coloredlines show the full analytical solutions [Eqs. (7), (8), (11), and (12)] with parameters
obtained in MD simulations (solid lines, denoted MD) and with fitted parameters (brokenlines, denoted Fit). The colored dotted lines denote Eq. (1) with fitted parameters. The
black broken line denotes σek = σ0 . We set ε = 78 [35] and T = 298 K. The interfacial
parameters are summarized in Table 1. Because the experimental data have been published
as continuous lines [31], we have extracted data points every 0.5 µC/cm2 (≈ 0.03 e/nm2 )for σ0 < 0 (symbols). In order to cancel the surface conduction correction [31], thecorrhydrophilic data have been reconstructed according to σek = sinh(arcsinh(σek)/1.55),corr
with σek corresponding to the value reported in reference [31]. (b) is a magnification of(a).

Fitting Eq. (1) to experimental data uniquely determines the equation’s prefactor, as well as the exponent γ. Because εint can be derived from capacitancemeasurements [22, 34], the procedure yields the interfacial viscosity ηint .
Fig. 1 shows a double-logarithmic plot of σek . The experimental data
(symbols) have been obtained with TiO2 (hydrophilic) and AgI (hydrophobic) in KNO3 aqueous solutions (explained in the caption of Fig. 1) [31].
Without interfacial effects σek = σ0 (black broken line). The deviation of
the data from the black broken line reveals the importance of the interfacial
effects. The colored solid and broken lines denote the solutions defined by
Eqs. (7) and (8) (see below), and the colored dotted lines denote Eq. (1). Theinterfacial parameters of the model are obtained in two different ways. First,
3

we take the parameters obtained from MD simulations of generic hydrophilic
and hydrophobic solid planar surfaces, which have been independently verified by comparison with multiple sets of experimental data on surface capacitance [22, 23], surface conduction and electrokinetics [24, 25]. Even though
the simulations were performed for diamond surfaces instead of the surface
types used in the experiments that we compare with, the advantage of thisapproach is the total absence of fit parameters. The specific adsorption parameter α is expected to depend sensitively on the combination of ions andsurfaces used, and is set to α = 0 due to a lack of literature data. For
hydrophilic surfaces, the slip length b and ηint are measurable in MD simulations [19]. Therefore, we use z ∗ = b(η/ηint − 1)−1 , and εint is determined by
the value of the surface capacitance εint /z ∗ [22, 23, 24, 25]. For hydrophobic surfaces, slip is caused by a low-density depletion layer, which can bemodeled by a slip boundary condition or by a low-viscosity layer. Either
description constitutes an approximation that captures the large-scale hydrodynamic properties well, but at the molecular scale, both approximations
are imperfect. We choose εint = 1, determining z ∗ through the capacitance,and determine ηint from b by ηint = ηz ∗ /(b + z ∗ ). For both surface types weobserve a fair agreement between the experimental data and our analyticaltheory without fit parameters (Fig. 1, solid lines). However, this procedure
without fit parameters predicts a higher slope than the one exhibited by theexperimental data for hydrophobic surfaces (blue squares).
As a second approach, we extract the interfacial properties from a fit to
the experimental data. At the hydrophilic surface, the exponent γ extracted
from the MD simulations provides a good fit to the experimental data (Fig. 1,
dotted lines). Excellent agreement is achieved by fitting only α (Fig. 1,
broken lines), resulting in the moderate value of α = 1.1. At the hydrophobic
surface, we fit Eq. 1 to the high-σ0 part of the curve shown in Fig. 1(b). As
the quality of the decades-old data does not allow for a 4-parameter fit, we
fix εint and fit z ∗ , ηint /η and α. Good agreement with the experimental
data is achieved using a 3.2 Å-wide low-viscosity layer (resulting slip lengthb = 1.6 nm). The differences between the MD and fit parameters are notunexpected, given that the MD simulations where performed for diamond andthe experiments for AgI. The stronger repulsion of ions from the interfacial
layer at hydrophobic surfaces (α = 10) than at hydrophilic surfaces (α = 1.1)
is consistent with ionic potentials of mean force calculated by MD simulationsof hydrophobic and hydrophilic self-assembled monolayers [25]. The fitted
and MD parameters used in Fig. 1 are summarized in Table 1.
4

Table 1: Interfacial parameters from molecular dynamics simulations (MD) [19, 22] andfrom fitting (Fit). The interfacial capacitance εint /z ∗ and the slip length b are shown forcomparing the difference between MD and Fit.

εinthydrophilichydrophobic

MD
Fit
MD
Fit

4.4
4.4
1.0
1.0

z∗
ηint /ηα(nm)0.44 3.7[19] 0.0
0.44
3.7
1.1
0.12 0.054 0.0
0.32
0.17
10

εint /z ∗b(nm−1 )
(nm)
10[22] -0.32[19]10-0.32
8.3[22]
2.1[19]3.1
1.6

Our theory and the comparison with experimental data reveal that theelectrokinetic surface charge σek exhibits a power law behavior with a small
but non-zero exponent (γ = 0.074 for the hydrophobic and γ = 0.015 for thehydrophilic surface), determined by the subnanometer-wide interfacial structure. The fitting procedure works particularly well for hydrophilic surfaces,
confirming the value of ηint derived from MD simulations. Note that charged
hydrophilic surfaces are also the most relevant in biological and industrial applications. Therefore, the newly-derived Eq. (1) provides a method to measure the parameters characterizing the interfacial structure, and in particular
the elusive interfacial viscosity, by standard electrokinetic experiments.
Model. — To derive Eq. (1) we consider a planar interface, so that theelectrostatic potential and all density distributions only depend on the perpendicular coordinate z. In the presence of an inhomogeneous local dielectric
profile, the Poisson equation reads [22, 23]


d
d
(2)
ε⊥ (z)ε0 ψ(z) = −ρ(z),dzdzwhere ε0 is the vacuum permittivity, ψ is the electrostatic potential, and ρ
is the ionic charge density. The dielectric profile ε⊥ (z), calculated using MD
simulations, shows strong oscillations [22, 23], but can be cast into a box
profile [22, 23, 24, 25, 33] to make the model analytically tractable,εint for 0 < z < z ∗ ,ε⊥ (z) =(3)
εfor z > z ∗ ,where εint is the interfacial dielectric constant, ε is the dielectric constantof bulk water, and z ∗ is the width of the interfacial layer [see Fig. 2(a)].
5The charge density equals ρ(z) = e[c+ (z) − c− (z)], where e is the elementary charge, and c+ (c− ) is the concentration of cations (anions). Neglecting
steric ion effect [36, 37, 38, 39], we model monovalent ions whose concentrations obey a Boltzmann distribution, c± (z) = c0 exp [∓Ψ(z) − µ± (z)], wherec0 is the bulk salt concentration and Ψ = eψ/kB T . Furthermore, µ+ (z)and µ− (z) are non-electrostatic contributions to the potential, for which weuse a box profile, µ± (z) = αθ(z ∗ − z), where α is an adjustable parameterand θ(z) is the Heaviside function [see Fig. 2(a)]. The boundary conditionon the surface reads dψ/dz|z=0 = −σ0 /εint ε0 , where σ0 is the bare surface
charge density, and far away from the surface ψ|z→∞ = 0. At the boundary
between the interfacial and the diffuse layer, we use continuity of the potential ψ(z)|z=z∗ −0 = ψ(z)|z=z∗ +0 , and continuity of the displacement field
εint dψ/dz|z=z∗ −0 = εdψ/dz|z=z∗ +0 . For an inhomogeneous viscosity profile,the hydrodynamics equation reads [24, 25]dudη⊥ (z)+ ρ(z)Ek = 0,
(4)
dz
dz
where η⊥ (z) is the viscosity profile, u(z) is the tangential velocity profile
and Ek is the applied tangential electric field. The zeta potential ζ is defined as the solvent flow infinitely far away from the surface [24, 25, 33]
ζ = −ηu|z→∞ /εε0 Ek . Combining Eqs. (2) and (4), and using the boundary
conditions du/dz|z→∞ = 0 and u|z=0 = 0, the zeta potential follows asZη ∞ ε⊥ (z) dψ
dz.
(5)
ζ=−
ε 0 η⊥ (z) dz
We approximate the viscosity profile by a box function,
ηint for 0 < z < z ∗η⊥ (z) =η
for z > z ∗ ,(6)where ηint is the viscosity in the interfacial layer [see Fig. 2(a)]. The width ofthe viscous interfacial layer generally differs from the dielectric one [24, 25],
depending on the assumed values of the interfacial dielectric constant εint
and viscosity ηint . For analytical tractability, we use the same width for
both, allowing εint and ηint to vary.
Clearly, the ion distributions in the interfacial layer are only well reproduced when using the full dielectric profile. Our approach is valid, nevertheless, because we calculate the potential difference and the velocity far away
6

Figure 2: (a) Illustration of the model of the interfacial layer with width z ∗ , which has adielectric constant εint , additional chemical potential kB T α, and viscosity ηint . The surface
has a bare surface charge density of σ0 . The region of z > z ∗ has a dielectric constantε and viscosity η. (b) The nondimensional electrostatic potential Ψ(z) for different bulk
salt concentrations (monovalent) and σ0 = 0.1 e/nm2 . We set ε = 78 [35], T = 298 K,z ∗ = 0.1 nm (broken line), εint = 1, and α = 0.0.

from the surface. In other words, we define the interfacial dielectric andviscous profile such that those quantities can be reproduced by continuum
modeling. It has been extensively verified that this box-based continuum
model effectively reflects the complex subnanometer interfacial structure well[22, 23, 24, 25]. It is this proven correspondence that allows us to use theanalytically tractable model presented here to reverse the calculation, and to
extract the effective microscopic parameters from macroscopic experiments
without relying on MD simulations and numerical calculations.
From Eqs. (3), (5) and (6) the zeta potential follows as [40]
ζ = γψ0 + (1 − γ) ψ ∗ ,(7)using the surface potential ψ0 = ψ(0) and the interfacial-boundary potentialψ ∗ = ψ(z ∗ ), and defining γ = εint η/εηint . Based on the Grahame equation,
the electrokinetic surface charge density is defined as [24, 25]


p
eζ
.
(8)
σek = 8εε0 kB T c0 sinh
2kB T
7For homogeneous systems, in which ηint = η and εint = ε, we obtain σek = σ0 .
In other words, deviations from σek = σ0 are a measure for interfacial effects.
Exact solution. — Eq. (2) leads to(ε/εint eα ) κ2 sinh Ψ for 0 < z < z ∗ ,
2
dΨ 
(9)
=
 2
dz 2
∗
κ sinh Ψ
for z > z ,where κ−1 is the Debye length defined via κ2 = 2e2 c0 /εε0 kB T . For z > z ∗ , thegeneral solution of Eq. (9) which satisfies the boundary condition at z → ∞,is given in terms of Ψ∗ = Ψ(z ∗ ) as∗

Ψ(z) = 2 ln

1 + e−(z−z )κ tanh(Ψ∗ /4)
for z > z ∗ .
1 − e−(z−z∗ )κ tanh(Ψ∗ /4)(10)

Next, we solve Eq. (9) in the range of 0 < z < z ∗ . Eq. (9) is symmetricfor cations and anions, so we restrict our discussion to a negatively charged
surface. The analytical solution of Eq. (9) for 0 < z < z ∗ is [41, 42]Ψ(z) =

i


h  p

−2arctanh cn −κ ε/εint eα (z − z ∗ ) − F arccos [− tanh(Ψ∗ /2)] 1 − p 1 − p ,





for 0 ≤ p ≤ 1,h  √ p


i



−2arctanh sn − pκ ε/εint eα (z − z ∗ ) + F arcsin [− tanh(Ψ∗ /2)] 1 − p−1 1 − p−1 ,for 1 < p,(11)

where cn(x|m) and sn(x|m) are Jacobian elliptic functions [42], F (ϕ|m) is
the elliptic integral of the first kind [42], and p is an integration constantwhich satisfies the relation p = [(eα ε/εint ) − 1](cosh Ψ∗ − 1)/2. The boundary
condition at z = 0 leads to

−

 p



p

−2κ ε/εint eα ds κ ε/εint eα z ∗ − F arccos [− tanh(Ψ∗ /2)] 1 − p 1 − p ,





for 0 ≤ p ≤ 1,


eσ0
=
√ p




εint ε0 kB T
√ p2 pκ ε/εint eα dcpκ ε/εint eα z ∗ + F arcsin [− tanh(Ψ∗ /2)] 1 − p−1 1 − p−1 ,for 1 < p,(12)

where ds(x|m) and dc(x|m) are Jacobian elliptic functions [42]. Fig. 2(b)shows the electrostatic potential profile, eq. (11), for different bulk salt
concentrations and σ0 = 0.1 e/nm2 . We use the parameters ε = 78 [35],T = 298 K, εint = 1, z ∗ = 0.1 nm, and α = 0. In the interfacial layer, thepotential decays faster than the region of z > z ∗ due to the low dielectric
constant. Both |Ψ0 | and |Ψ∗ | decrease with increasing salinity.
8

1

5

10

4

0

10

-̀*

3

10-1

linear behavior

2
10-2

1

-3

10

0

10-4 -4
10 10-3 10-2 10-1 100 101(b)-̀0 [e/nm2]

-1
0

(a)

0.1

1mM
10mM
100mM
1M

0.2

-̀0 [e/nm2]

Figure 3: Eq. (12) in terms of the interfacial-boundary potential Ψ∗ (solid lines). We use
the parameters ε = 78 [35], T = 298 K, εint = 1, z ∗ = 0.1 nm, and α = 0. (a) is a linearscale graph, whereas (b) is a double-logarithmic one. Dotted lines in (a) are calculated bylinear analysis [33], and broken lines in (b) denote Eq. (13) for 1 mM and 10 mM.

From asymptotic analysis of Eqs. (11) and (12) in the limit of strongly
charged surfaces and low salinity, the interfacial-boundary potential saturates
at the value
#
"
p
2
α ε/ε
4arccosh
e
int
,
(13)
Ψ∗ ≈ − ln(ε/εint − e−α )(ε/εint )κ2 z ∗2and the surface potential Ψ0 = Ψ(0) becomes"2 #eα εinteσ0Ψ0 ≈ − ln.
κ2 εεint ε0 kB T

(14)

Substituting Eqs. (13) and (14) into Eq. (7), and substituting Eq. (7) into
Eq. (8) leads to Eq. (1).
Fig. 3 shows the interfacial-boundary poten∗
tial Ψ for different salt concentrations as a function of the bare surface
charge density σ0 . In Fig. 3(a), the solid lines denote Eq. (12), and thedotted lines denote the linearized solution [33] for parameters ε = 78 [35],T = 298 K, εint = 1, z ∗ = 0.1 nm, and α = 0. For weakly charged surfaces, Ψ∗ increases linearly with the surface charge density. However, abovea specific σ0 depending on the salinity, Ψ∗ shows a nonlinear dependence.
Fig. 3(b) (double-logarithmic) clearly shows the saturation of Ψ∗ in the limitof strongly charged surfaces. Horizontal broken lines indicate Eq. (13). The
saturation of Ψ∗ is not unique for this modified PB equation, but also occurs
9

20

20

-̀0

15

15

1mM
10mM
100mM
1M

10
ln 10

10
5
5

0
-5
0

(a)

0.1

0
10-3

0.2
2

(b)

-̀0 [e/nm ]

10-2

10-1

100
2

101

-̀0 [e/nm ]Figure 4: The surface potential Ψ0 as a function of the bare surface charge density σ0(solid lines). We use the parameters ε = 78 [35], T = 298 K, εint = 1, z ∗ = 0.1 nm, and
α = 0. (a) is a linear-scale graph, whereas (b) is a semi-logarithmic one. Dotted lines in(a) are calculated by linear analysis [33], and broken lines in (b) denote Eq. (14).

for a homogeneous dielectric constant [2]. Fig. 4 shows the surface potential
Ψ0 at a given salt concentration. In Fig. 4(a) (linear scale), the solid lines
denote the numerical result obtained by solving Eqs. (11) and (12), and thedotted lines denote the linearized solution [33]. Broken lines in Fig. 4(b)(semi-logarithmic scale) indicate the dependence |Ψ0 | ∼ ln(σ02 /c0 ) for large|σ0 |, which follows from Eq. (14).
Summary and conclusions. — We derive an analytical expression for theelectrokinetic surface charge density σek using a universal model of the structure of the interfacial layer. The model consists of box profiles for the viscosity, the dielectric constant and the specific ion adsorption. In the limitof strongly charged surfaces and low salinity, |σek | increases as a functionof the bare surface charge density as a power law with nonuniversal exponent γ = εint η/εηint . The exponent is larger for hydrophobic surfaces (whereηint < η) than for hydrophilic surfaces (where ηint > η), in agreement withexperimental data. Fitting the analytical power law to experimental data ofσek as a function of the bare surface charge density σ0 gives direct access tothe interfacial quantities z ∗ , εint , α and in particular the elusive effective interfacial viscosity ηint . Therefore, our analytical expression provides a novel
method to extract the subnanometer-scale structure of the interfacial layer
directly from conventional macroscopic experiments. As the available experimental data is limited and dates back to mid last century, we hope that our10present results will stimulate new experimental studies.
Acknowledgements
YU was supported by Grant-in-Aid for JSPS Fellows 13J00063 and 16J00042.
[1] S. S. Dukhin and B. V. Derjaguin, Surface and Colloid Science, ed. E.
Matijrvic, vol. 7, Wiley, New York, (1974).
[2] W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal Dispersions, Cambridge University Press, Cambridge, (1989).
[3] J. Lyklema, Fundamentals of Interface and Colloid Science, vol.3, Academic Press, London (1995).
[4] L. Joly, C. Ybert, E. Trizac, and L. Bocquet, Phys. Rev. Lett. 93 (2004)
257805.
[5] L. Joly, C. Ybert, E. Trizac, and L. Bocquet, J. Chem. Phys. 125 (2006)
204716.
[6] A. Ajdari and L. Bocquet, Phys. Rev. Lett. 96 (2006) 186102.
[7] D. M. Huang, C. Cottin-Bizonne, C. Ybert, and L. Bocquet, Phys.
Rev. Lett. 98 (2007) 177801.
[8] D. M. Huang, C. Cottin-Bizonne, C. Ybert, and L. Bocquet, Langmuir
24 (2008) 1442–1450.
[9] C. I. Bouzigues, P. Tabeling, and L. Bocquet, Phys. Rev. Lett. 101(2008) 114503.
[10] O. Bonhomme, O. Liot, A.-L. Biance, and L. Bocquet, Phys. Rev. Lett.
110 (2013) 054502.
[11] L. Joly, F. Detcheverry, and A.-L. Biance, Phys. Rev. Lett. 113 (2014)088301.
[12] S. R. Maduar, A. V. Belyaev, V. Lobaskin, and O. I. Vinogradova,
Phys. Rev. Lett. 114 (2015) 118301.

11[13] K. Kim, Y. Nakayama, and R. Yamamoto, Phys. Rev. Lett. 96 (2006)
208302.
[14] V. Lobaskin, B. Dünweg, M. Medebach, T. Palberg, and C. Holm,
Phys. Rev. Lett. 98 (2007) 176105.
[15] G. Giupponi and I. Pagonabarraga, Phys. Rev. Lett. 106 (2011) 248304.
[16] M. S. Kilic, M. Z. Bazant, and A. Ajdari, Phys. Rev. E 75 (2007)
021502.
[17] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Phys. Rev. Lett. 106(2011) 046102.
[18] A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell, and L. Bocquet, Nature 494 (2013) 455–458.
[19] C. Sendner, D. Horinek, L. Bocquet, and R. R. Netz, Langmuir 25(2009) 10768–10781.
[20] J.-L. Barrat and L. Bocquet, Phys. Rev Lett. 82 (1999) 4671–4674.
[21] D. M. Huang, C. Sendner, D. Horinek, R. R. Netz, and L. Bocquet,
Phys. Rev. Lett. 101 (2008) 226101.
[22] D. J. Bonthuis, S. Gekle and R. R. Netz, Phys. Rev. Lett. 107 (2011)
166102.
[23] D. J. Bonthuis, S. Gekle and R. R. Netz, Langmuir 28 (2012) 7679–
7694.
[24] D. J. Bonthuis, R. R. Netz, Langmuir 28 (2012) 16049–16059.
[25] D. J. Bonthuis and R. R. Netz, J. Phys. Chem. B 117 (2013) 11397–
11413.
[26] K. Takae and A. Onuki, J. Chem. Phys. 143 (2015) 154503.
[27] Y. Zhu and S. Granick, Phys. Rev. Lett. 87 (2001) 096104.
[28] O. I. Vinogradova and G. E. Yakubov, Langmuir 19 (2003) 1227–1234.

12

[29] C. Cottin-Bizonne, B. Cross, A. Steinberger and E. Charlaix, Phys.
Rev. Lett. 94 (2005) 056102.
[30] J. Lyklema, J. T. G. Overbeek, J. Colloid Sci. 16 (1961) 501–512.
[31] J. Lyklema, Colloid Surf. A 92 (1994) 41–49.
[32] J. Lyklema, Curr. Opin. Colloid Interface Sci. 15 (2010) 125–130.
[33] D. J. Bonthuis, Y. Uematsu, and R. R. Netz, Phil. Trans. R. Soc. A
374 (2016) 20150033.
[34] O. Stern, Zeitschrift für Elektrochemie, 30 (1924) 508.
[35] B. B. Owen, R. C. Miller, C. E. Milner, and H. L. Cogan, J. Phys.
Chem. 65 (1961) 2065.
[36] J. J. Bikerman, Philos. Mag. 33 (1942) 384–397.
[37] M. Eigen and E. Wicke, J. Phys. Chem. 58 (1954) 702–714.
[38] I. Borukhov, D. Andelman, and H. Orland, Phys. Rev. Lett. 79 (1997)
435–438.
[39] I. Borukhov, D. Andelman, and H. Orland, Electrochim. Acta. 46(2000) 221–229.
[40] Y. Uematsu and T. Araki, J. Chem. Phys. 139 (2013) 094901.
[41] S. H. Behrens and M. Borkovec, Phys. Rev. E 60 (1999) 7040–7048.
[42] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, 9th ed. Dover Publications, New York, (1972).

13