PHYSICAL REVIEW E 98, 032604 (2018)

Linear response functions of an electrolyte solution in a uniform flow
Ram M. Adar
Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Ramat Aviv, Tel Aviv 69978, Israel
and Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University, Tokyo 192-0397, Japan

Yuki Uematsu
Department of Chemistry, Kyushu University, Fukuoka 819-0395, Japan

Shigeyuki Komura*
Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University, Tokyo 192-0397, JapanDavid Andelman†
Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Ramat Aviv, Tel Aviv 69978, Israel(Received 19 June 2018; published 13 September 2018)
We study the steady-state response of a dilute monovalent electrolyte solution to an external source with aconstant relative velocity with respect to the fluid. The source is taken as a combination of three perturbations:an external force acting on the fluid, an externally imposed ionic chemical potential, and an external charge
density. The linear response functions are obtained analytically and can be decoupled into three independent
terms, corresponding to (i) fluid flow and pressure, (ii) total ionic number density and current, and (iii) charge
density, electrostatic potential, and electric current. It is shown how the uniform flow breaks the equilibrium
radial symmetry of the response functions, leading to a distortion of the ionic cloud and electrostatic potential,
which deviates from the standard Debye-Hückel result. The potential of a moving charge is underscreened in its
direction of motion and overscreened in the opposite direction and normal plane. As a result, an unscreened
dipolar electric field and electric currents are induced far from the charged source. We relate our generalformalism to several experimental setups, such as colloidal sedimentation.
DOI: 10.1103/PhysRevE.98.032604
I. INTRODUCTION

Ionic solutions are found in a wide range of biologicalsystems and are used for a plethora of industrial applications
and processes. For ionic solutions out of thermodynamic equilibrium, many well-known physical phenomena rely on theinterplay between Coulombic interactions, hydrodynamics,
and thermal diffusion [1–4]. For example, an applied electricfield can induce an electrolyte flow in a capillary (electroosmosis) [5–11]. Similarly, a salt concentration gradient in acapillary causes a water flow (diffusio-osmosis) [12,13] andelectric currents (osmotic current) [12,14].
For colloidal suspensions, transport of charged colloids
can be achieved by applying an external electric field (electrophoresis) [15–19]. The distortion of the electric double
layer near the colloidal surface, in the presence of the appliedfield, significantly decreases the electrophoretic mobility forhigh surface potentials [15]. Moreover, when charged colloids
sediment under gravity in an electrolyte solution, a potential
difference builds along the direction of gravity, called sedimentation potential [20,21].
Electrokinetic processes can be described by the linear
response theory of the system to external sources, when the*†

komura@tmu.ac.jp
andelman@post.tau.ac.il

2470-0045/2018/98(3)/032604(10)sources are sufficiently weak. The response in the quiescentstate is well known when the hydrodynamics and electrostatics are completely decoupled. In such cases, the linear
response of the velocity field with respect to a point force
is described by the Oseen tensor [22,23] and that of theelectrostatic potential to a point charge is described by the
Debye-Hückel (DH) theory [24].
In this work, we investigate a different scenario and consider the linear response of a bulk electrolyte in a uniform
flow. Such a flow can correspond to an electrolyte flowing
past a stationary object, e.g., an optically trapped colloid orparticle. Alternatively, we can think of a source moving inan otherwise stationary electrolyte. This source, for example,
can be an active biomolecule or the tip of an atomic force
microscope (AFM).
We present several generalized response functions of asystem in a uniform and stationary flow, taking into accountthe combination of three different external sources: (i) aforce density acting on the fluid, (ii) an externally imposedionic chemical potential, and (iii) an external charge density.
Although the response of the velocity field is still independentof the electrostatic interactions, we show that the response of
charge density, ionic number density, and their currents arecoupled to the uniform flow. In particular, the response of
the charge density to an external point charge is described bya nonlinear function of the fluid velocity, extending the DH
screening theory for a uniform flow.

032604-1©2018 American Physical Society

ADAR, UEMATSU, KOMURA, AND ANDELMANPHYSICAL REVIEW E 98, 032604 (2018)

(b)

(a)

+- ++

F
u0

(c)

+

+

P,u

u0

+
+
+
-

-

- +

n,J

-

+

q
-

u0

+

-

+

,c,j

FIG. 1. A schematic illustration of three possible external sources (in red) and their corresponding induced fields (in gray box). All thesources are at rest, while the electrolyte flows with a constant velocity u0 = u0 ẑ. (a) A thin rod exerting a force density, F, in the negative z
direction, leading to a pressure, P , and velocity field, u, in the electrolyte. (b) A tip of a pipette containing the electrolyte with a different ionicchemical potential, μ, causing changes in the total ionic number density, n, and number density current, J. (c) A colloid with charge density,q, producing an electrostatic potential, ψ, as well as a charge density, c, and electric current, j . The partition into sources and consequent
fields is described in detail in Sec. II C. Any surface effects stemming from the boundaries of the sources are neglected.

The outline of this paper is as follows. In Sec. II, we present
our model for a dilute electrolyte and a general external perturbation with a constant relative velocity. The basic electrohydrodynamic equations of our model are derived (Sec. II B),
and the linearized scheme is described (Sec. II C). Next,
in Sec. III, the general response functions are derived andanalyzed in detail. They include the hydrodynamic response(Sec. III A), the number density response (Sec. III B), and theelectric response (Sec. III C). In Sec. IV, we summarize ourresults and relate our findings to possible experiments.

II. MODEL
A. Three sources and their responses

Consider a dilute ionic solution, consisting of 1:1 monovalent cations and anions of total bulk concentration n0 ,immersed in a continuum solvent of dielectric permittivity
ε. The solvent is modeled as an incompressible fluid withviscosity η. The ions are assumed to be pointlike, and thefriction coefficient of both cations and anions with the solvent
is ζ . The system is held at constant temperature T . Under
these conditions, the homogeneous solution is in thermal
equilibrium and satisfies local electroneutrality.
The homogeneous electrolyte can be perturbed by an externally controlled source. When the object is at rest, the systemreaches a new equilibrium state. However, when it is mobile
at a constant velocity, a steady state can be established, where
the system is out of equilibrium, but all physical quantities
are time independent. We explore this latter scenario of amobile external perturbation with a relative velocity withrespect to the ionic solution. For convenience, the frame ofreference is chosen such that the perturbation is at rest, while
the electrolyte flows with a velocity u0 .
As a general perturbation, we consider the combination of
three possible sources depicted in Fig. 1: (a) an external force
density, F, e.g., a force exerted by a thin rod immersed in thesolution; (b) an externally imposed ionic chemical potential,
μ, e.g., a tip of a pipette containing the electrolyte solution

with a different ionic concentration; and (c) an external charge
density (per unit volume), q, e.g., a small charged colloid.
We assume that surface effects from the sources are
negligible, such that no boundary conditions are imposed.
This assumption is appropriate for small sources or forthe linear response in the far-field away from the sources.
For example, a colloidal probe AFM with a fixed scanningvelocity can be modeled as a combination of a point charge,q, and a point force density, F, originating from the drag on
the colloidal probe.
Each of the three sources affects the electrolyte differently
and can be identified with different fields that characterize theelectrolyte, as is clarified in detail in Sec. II C and indicatedin Fig. 1. The force density in Fig. 1(a), F, perturbs theelectrolyte velocity, u, and pressure, P . In Fig. 1(b), theexternally imposed chemical potential, μ, modifies the total
ionic number density, n, and current, J, defined asn = n+ + n− ,J = n+ v + + n− v − − nu.

(1)

Here n± are the densities of cations and anions (n± = n0 in
the homogeneous bulk), and v ± are their velocities. Note thatthe current is evaluated in the moving frame of reference.
The external charge density of Fig. 1(c), q, results in anelectrostatic potential, ψ, a charge density (per unit volume),c, and an electric current, j . The latter two are defined asc = e(n+ − n− ),j = e(n+ v + − n− v − ) − cu,(2)where e is the elementary charge. Similarly to the number
density current J, the electric current j is defined according
to the ionic relative velocities, v ± − u.
B. Electrohydrodynamic equations

The response of the seven fields, P , u, n, J, ψ, c, andj , to the three sources, F, μ, and q, is obtained by solving
seven coupled differential equations; the Stokes equation and032604-2

LINEAR RESPONSE FUNCTIONS OF AN ELECTROLYTE …

incompressibility condition for the fluid, Poisson’s equation
for the electrostatic potential, force balance equations forthe two ionic densities, and the corresponding continuity
equations for the two currents. All of these equations canbe derived consistently within a single framework, using
Onsager’s variational principle [1,25]. Below, the equations
are presented and discussed in detail.
First, the fluid is incompressible:∇ · u = 0.

(3)

For low Reynolds numbers (no inertia), the electrolyte satisfies the Stokes equation,
η∇ 2 u − ∇P − c∇ψ = −F.

(4)

The first term in Eq. (4) originates from the solvent viscosity,
η. The third term stems from the solute charge, c, and couples
the hydrodynamics with the electric variables.
It is possible to decompose the pressure P and the force
density F into hydrodynamic and thermal terms. The pressure
is given by the sum of the hydrodynamic solvent pressure andsolute pressure as explained below. Similarly, the force density is the combined external force density on the electrolyte,
f , and diffusive thermal force, i.e., F = f − n∇μ. As mentioned above, μ is the externally imposed chemical potentialof the ions, taken to be the same for the cations and anions.
The electrostatic potential, ψ, in Eq. (4) satisfies Poisson’sequation (in SI units),ε∇ 2 ψ + c = −q,(5)where ε is the dielectric permittivity of the medium and q is
the external charge density. Including both the contributionsof the electric field, E = −∇ψ, and the hydrodynamic drag,we obtain the force balance equations for the cations andanions as

PHYSICAL REVIEW E 98, 032604 (2018)

In principle, additional source terms can be included in theright-hand side of Eq. (8) [26,27]. However, in this study,
we neglect such extra ionic source terms, as well as any
possible chemical reactions that lead to additional fluxes in
the continuity equation.
In the absence of electrolyte flow and an externally imposed chemical potential, i.e., u = 0 and μ = 0, our setof equations reduces to the steady-state Poisson-NernstPlanck equations [5]. In thermodynamic equilibrium, theseequations reduce further to the Poisson-Boltzmann equation[5,28]. Finally, for small electrostatic potentials, the PoissonBoltzmann equation can be approximated by its linear form,
which is the DH equation [28].
It is convenient to rescale all variables and define thefollowing dimensionless quantities:
ζ λD
4π lB
P
,
u≡
u,

η≡
η,
n0 kB T
kB T
ζ
μ
λD
ζ λD
n


J, 
μ≡
,
F≡
F,
J≡

n≡ ,
n0
n0 kB TkB T
n0 kB Teψcqj ≡ ζ λD j , 
≡
ψ
, 
c≡
q≡
,
,
kB T
en0
en0 kB T
en0
(9)
≡
P

where lB = e2 /(4π εkB T ) is the Bjerrum length and λD =(4π lB n0 )−1/2 is the Debye screening length (n0 is defined as
the combined concentration of cations and anions together).
Similarly, the position vector is rescaled as r ≡ r/λD so that ≡ λD ∇.
∇Hereafter, the tilde notation is omitted and all variables
are treated as dimensionless quantities. Hence, the set of
equations, Eqs. (3)–(8), can be written in a compact form asη∇ 2 u − ∇P − c∇ψ = −F,∇ 2 ψ + c = −q,−en+ ∇ψ − ζ n+ (v + − u) − kB T ∇n+ = n+ ∇μ,(6)

J + c∇ψ + ∇n = −n∇μ,where kB T is the thermal energy. We combine the force
balance equations above into two new equations in terms ofthe number density current, J, and the electric current, j :j + n∇ψ + ∇c = −c∇μ,en− ∇ψ − ζ n− (v − − u) − kB T ∇n− = n− ∇μ,∇ · u = 0,
∇ · J + u · ∇n = 0,

ζ J + c∇ψ + kB T ∇n = −n∇μ,ζ j + e2 n∇ψ + kB T ∇c = −c∇μ.

(7)These equations couple the number density, n, and current, J,
with the electric variables, ψ, c, and j (electrostatic potential,
charge density, and electric current, respectively).
Note that the term kB T ∇n± in Eq. (6) that is carried over to
Eq. (7) can be regarded as the gradient of van ’t Hoff ideal gas
pressure. Here we assume that the system is well describedby such a term, as it is not far from thermal equilibrium. It
follows from the ideal-gas pressure form that steric effects andfluctuations of the electrostatic potential beyond the meanfield treatment are not included within our framework.
Finally, we combine the continuity equation for each of thetwo ionic flows, ∇ · (n± v ± ) = 0, and arrive at the following
equations in terms of the two currents:∇ · J + u · ∇n = 0,
∇ · j + u · ∇c = 0.

(8)

∇ · j + u · ∇c = 0.

(10)

A general solution of these equations requires to specify theboundary conditions in terms of fixed (and charged) objects
and interfaces. However, as we are interested only in bulk
properties far from any boundaries, such conditions will notenter into our study.
C. Linearized equations

In the absence of sources, F = q = μ = 0, all the equations in Eq. (10) are homogeneous and describe a uniform
bulk electrolyte in a constant flow. The ions are at rest in
the moving fluid reference frame. They are distributed homogeneously and satisfy local electro-neutrality. This special
solution is given by c = ψ = J = j = 0 while u = u0 and
n = P = 1 (for the dimensionless variables).
We now analyze the system response to a small perturbation, as is described above, via a linearization of Eq. (10)032604-3

ADAR, UEMATSU, KOMURA, AND ANDELMANPHYSICAL REVIEW E 98, 032604 (2018)

around the homogeneous solution. Each field, (r ), canbe expanded as (r ) = 0 + 1 (r ) + · · · , where 0 is thehomogeneous field and 1 (r ) is the linear correction in
the presence of sources. All the linear corrections as well as
the sources are considered to be of a similar small magnitude.
We keep only terms that are linear in this small magnitudeand neglect quadratic or higher-order terms. The linearization
procedure is described in detail in Appendix A.
The linearized equations can be written as three decoupledsets of equations. The first set includes the hydrodynamic
variables [see Fig. 1(a)] and reads

response has a convolution form:(r ) = d 3 r  K (r − r  )S (r  ) .

Below we present separately the response functions foreach of the three sets: (i) the response of the hydrodynamic
variables, u1 and P1 , to the source F; (ii) the response of thenumber density and current variables, n1 and J 1 , to the source
μ; and (iii) the response of electric variables, ψ1 , c1 , and j 1 , tothe source q. For a full derivation of these response functions,
see Appendix A.

η∇ 2 u1 (r ) − ∇P1 (r ) = −F(r ),∇ · u1 (r ) = 0.

A. Response of hydrodynamic variables(11)

These are the Stokes equation and incompressibility condition
for the electrolyte. They describe how the force density, F,
induces a pressure gradient that results in a fluid flow. Theseequations are independent of electrostatics, while the ionicnumber density (and not charge density) enters the definitionof the pressure, P , and force density, F.
The second set of equations corresponds to the number
density, n1 , and current, J 1 [see Fig. 1(b)]:J 1 (r ) + ∇n1 (r ) = −∇μ(r ),∇ · J 1 (r ) + u0 · ∇n1 (r ) = 0.

(12)

The first equation describes how the externally imposed chemical potential, μ, induces a number density gradient, ∇n1 ,and, consequently, yields a number density current, J 1 . The
second equation in Eq. (12) implies that the current J 1 isa compressible field, since it is defined relative to the fluid
velocity. Combining the two equations yields a steady-state
convection-diffusion equation [3] for the number density, n1 .
Finally, we obtain for the electric variables [see Fig. 1(c)]∇ 2 ψ1 (r ) + c1 (r ) = −q(r ),j 1 (r ) + ∇c1 (r ) = −∇ψ1 (r ),∇ · j 1 (r ) + u0 · ∇c1 (r ) = 0.

(14)

(13)

The three equations above describe how an external charge, q,
produces an electric field and a charge density gradient, whose
combined effect leads to an electric current. The first equationin Eq. (13) is Poisson’s equation that remains unchanged by
the linearization procedure. The bottom two equations havethe same structure as in Eq. (12), where ψ1 plays the roleof μ. The reason is that the chemical potential and number
density are conjugate variables, just as the charge density andelectrostatic potential are.

Consider an external force density imposed on the fluid,
F, as in Fig. 1(a). Solving Eq. (11) for the hydrodynamic
response (Appendix A), we arrive ati
P1 (k) = − 2 k · F(k),
k


1
kk
u1 (k) = 2 I − 2 · F(k),
ηk
k

(15)

where I is the identity tensor, k = |k|, and kk is the dyadic
product (second-rank tensor). Note that Eq. (15) is independent of ionic properties.
Transforming the hydrodynamic responses of Eq. (15) into
real space yields the response functionsr,
KP (r ) =4π r 3rr 
1 
Ku (r ) =I+ 2 ,(16)8π ηr
r
where r = |r|. We emphasize that KP is a vector and Ku is asecond-rank tensor. As the linearized hydrodynamic equations
are independent of electrostatics, Eq. (16) restores the Oseen’s
result with Ku being the Oseen tensor [22,23]. The response
to a point force source F = −δ(r )ẑ (ẑ is a unit vector along
the z direction) of Eq. (16) is illustrated in Fig. 2. As the Oseen’s result is well known, we present it only for completenessand do not discuss it further.
B. Response of number density and current

Consider an externally imposed chemical potential, μ, of
the ionic species, as in Fig. 1(b). The response to such a sourceis given by solving Eq. (12) (see also Appendix A),

III. LINEAR RESPONSE FUNCTIONS

k2
n1 (k) = − 2μ(k),k + iu0 · k(u0 · k)kμ(k).
J 1 (k) = 2
k + iu0 · k

The solutions to each of the three sets of linear equations in Eqs. (11)–(13) can be written in terms of linear
response functions. Such solutions are conveniently found via
the Fourier transform of the fields and sources, where theFourier
 3 transform of any function (r ) is defined as (k) =d r (r ) exp(−ik·r ). In Fourier space, the response of thesystem in terms of a variable  to an external source S
is given by the product, (k) = K (k)S (k), where K (k)is the response function (kernel). In real space, the same

We see here that the two fields are related by J 1 (k) =−(u0 · k̂) k̂ n1 (k), where k̂ = k/k is a unit vector.
In the absence of flow, u0 = 0, the response reduces to
n1 = −μ and J 1 = 0. Such a response can be obtained fromthe equilibrium Boltzmann distribution for small μ values because n1 = e−μ − 1 ≈ −μ. Otherwise, for a nonzero velocity,
the direction of u0 (taken in the z direction) defines an axisof symmetry, and the usual spherical symmetry becomes anazimuthal one (body of revolution around u0 ).

032604-4(17)

LINEAR RESPONSE FUNCTIONS OF AN ELECTROLYTE …

PHYSICAL REVIEW E 98, 032604 (2018)

P

u1

( )

10 (b)

z

z

5
0
5
10
10

5

x

0

x

5

10

FIG. 2. Response of the hydrodynamic fields, pressure P1 and fluid velocity u1 , to a point force source at the origin, F = −δ(r )ẑ [seeEq. (16)]. All quantities are dimensionless, according to the definitions of Eq. (9). The response is plotted in the xz plane, and can be extended
to the entire space due to the azimuthal symmetry around the z axis. (a) Contour plots of equal pressure. The red regions in the upper half-plane
correspond to negative values, P1 < 0, while the blue regions in the bottom half-plane correspond to positive ones, P1 > 0. Note that the innerarea corresponds to larger absolute values of P1 and is left white for clarity sake. (b) Stream lines of the velocity u1 . The arrows indicate theflow direction.

The real-space response of the number density, n1 , andcurrent, J 1 , are given by the two response functions

obtained by solving Eq. (13) are given by

Kn (r ) = −δ(r ) + g(r ),c1 (k) = −K J (r ) = −∇g(r ),(18)

where Kn is a scalar, K J is a vector, and the function g is
given byu0 r 2 − z(2 + u0 r ) u0 (z−r )/2g(r ) = u0e.
(19)
8π r 3
We note that the Green’s function for the corresponding
two-dimensional convection-diffusion equation was solved inRef. [29].
The expressions of Eqs. (18) and (19) for n and J are
plotted in Fig. 3. It is convenient to interpret this case byconsidering a source moving in an otherwise stationary electrolyte with a velocity −u0 . Due to friction, ions are pushed
by the source outside of its trajectory into an outer region. The
boundary of this region can be obtained from Eq. (19) as thecontour of n1 = 0 and is given byr=

2 cos θ.
u0 1 − cos θ(20)

In the above, θ is the azimuthal angle between r and ẑ, and thiscontour is drawn as a black line in Fig. 3(a). Consequently, awake is formed, as is seen in Fig. 3(b). Ions within the wake
experience a rotational flow, while ions in the outer region
flow away from the moving source, with an average velocity
in the negative z direction.
C. Response of electric variables

We now consider the third case where the electrolyte is
perturbed by an external charge source, q, as in Fig. 1(c).
In response, mobile ions surround the source and an electric
current is formed. The charge density and electric current

j 1 (k) =1q(k),1 + iu0 · k + k 2(u0 · k̂) k̂
q(k).
1 + iu0 · k + k 2(21)

For k = 0, Eq. (21) yields c1 (k = 0) = −q(k = 0), satisfying
electroneutrality. In the limit of u0 = 0, on the other hand, thestandard DH result is obtained, i.e., c1 (k) = −q(k)/(1 + k 2 ).
The charge density and electric current of Eq. (21) are related
to each other by j 1 (k) = −(u0 · k̂) k̂ c1 (k), just as the number
density and current in Eq. (17) are related.
Performing the inverse Fourier transform of Eq. (21) results in the real-space response functions
⎛
⎞
2
u
z
1
u
0
Kc (r ) = −exp ⎝− 1 + 0 r ⎠,4π r
2
4K j (r ) = ∇[u0 · ∇Kc (r )],(22)

where Kc is a scalar and K j is a vector. The first lineof Eq. (22) describes the “relaxation effect” of the ioniccloud. Ions are dragged in the direction of the flow, and theotherwise-spherical charge density of mobile ions is stretchedin the flow direction, as is plotted in Fig. 4(a). This effect isknown as the Dorn effect [21] in the case of sedimentation,
where counterions accumulate behind the sedimenting particle. A similar distortion of the ionic cloud has been previously
reported in the presence of a shear flow [30].
As a consequence of the relaxation effect, the charged
source together with the ionic cloud have a nonzero dipole
moment. Assuming a pointlike source q(r ) = Qδ(r ), we
integrate Kc (r )r and find the dipole moment −Qu0 . Notethat all the physical quantities here are dimensionless, as is
explained above. For example, a negatively charged source

032604-5

ADAR, UEMATSU, KOMURA, AND ANDELMANPHYSICAL REVIEW E 98, 032604 (2018)n1

J1
15 (b)

z

z

10
5
0
5
10

x

5

0

x

5

10

FIG. 3. Response of the number density, n1 , and current field, J 1 , to a point source at the origin, μ = δ(r ), for the fluid velocity u0 = 2ẑ[see Eqs. (18) and (19)]. (a) Contour plots of equal number densities n1 . The contour n1 = 0 is drawn by the black line [see Eq. (20)]. Notethat the inner (uncolored) areas correspond to even larger absolute values of n1 . (b) Stream lines of the number density current, J 1 .

moving in the negative z direction results in a dipolar momentin the positive z direction, as can be inferred from Fig. 4(a).
At large distances, this dipole dominates the electric field, as
is further discussed below.
While the ionic cloud is well approximated by a dipole
far from the sources, the exact form of the charge density
becomes important closer to it. By examining the argumentof the exponent in Eq. (22), it is possible to approximate theshape of the ionic cloud by a prolate ellipsoid, satisfying theequationρ2
u0 a 2
1


z
−
+
= 1.
a2
2
a 2 1 + u20 /4(23)√minor semiaxis and a 1 + u20 /4 is the major one. Both the
major axis and one of the foci of the ellipsoid increase with
u0 . Such a description fits well with the ionic cloud illustrated
in Fig. 4(a).
Another important feature of the nonequilibrium steadystate is the formation of an electric current, as is given by
the second line of Eq. (22). An example for such a current isplotted in Fig. 4(b). The symmetric closed loops around the zaxis resemble those of a dipolar field. This observation can bebetter understood by examining Eq. (21) for large distances.
Keeping only the lowest order of k leads to the electric current
response, (u0 · k̂) k̂. The inverse Fourier transform of thisfunction is obtained as a second derivative of the kernelIn the above equation, ρ = x 2 + y 2 is the absolute value ofthe in-plane two-dimensional vector ρ = (x, y), while a is theK j (r ) ≈11[u0 − 3(u0 · r̂ )r̂],
4π r 3(24)

115 (b)

z

z

10
5
0
5
10

x

5

0

x

5

10

FIG. 4. Response of the charge density, c1 , and electric current, j 1 , to a negative point charge source at the origin, q = −δ(r ) for the fluid
velocity u0 = 2ẑ [see Eq. (22)]. (a) Contour plots of equal charge densities c1 . Note that the inner (uncolored) area corresponds to even larger
absolute values of c1 . (b) Stream lines of the electric current j 1 .
032604-6

PHYSICAL REVIEW E 98, 032604 (2018)

10

10

5

5

0

0

z

z

LINEAR RESPONSE FUNCTIONS OF AN ELECTROLYTE …

5

5

10

10
10

5

0

x

5

10

10

5

0

x

5

10

FIG. 5. Response of of the electrostatic potential, ψ1 , and electric field, E 1 , to a negative point charge source at the origin q = −δ(r )for the fluid velocity u0 = 2ẑ [see Eqs. (27) and (28)]. (a) Contour plots of equipotential values of ψ1 . Note that the inner (uncolored) area
corresponds to even larger absolute values of ψ1 . (b) Stream lines of the electric field E 1 .

which indeed corresponds to a field of a dipolar moment −u0 .
The charge density and electric current of Eqs. (21)and (22) are accompanied by an electrostatic potential, ψ1 .
By solving Eq. (13), we find that


1
1
q(k).
(25)
ψ1 (k) = 2 1 −k
1 + iu0 · k + k 2
The prefactor 1/k 2 is the Fourier transform of the Coulomb
kernel, 1/(4π r ). The first term inside the parenthesis corresponds to the direct effect that a charged source has on thepotential. As is evident from Eq. (21), the second term corresponds to the ionic cloud. Furthermore, substituting ψ1 (k) in
Eq. (25) into Eq. (21) yields
k2
ψ1 (k),c1 (k) = − 2
k + iu0 · k(u0 · k)k
ψ1 (k),
j 1 (k) = 2
k + iu0 · k

14π ρ 2 + z2 ∞1 − exp[h(ρ, z − z )] −z /u0
1
−
dz 
e
,
4π u0 0
ρ 2 + (z − z )2(27)

where we have denotedh(ρ, z) =

u0
z−
2



u2
1 + 0 (ρ 2 + z2 ).
4

IV. SUMMARY AND DISCUSSION(26)

demonstrating the same dependence on ψ1 as that of n1 andJ 1 on μ [see Eq. (17)], as mentioned before.
Equation (25) generalizes the DH result for sources moving with a relative velocity with respect to the solvent. We
transform it back to real space to obtainKψ (ρ, z) =plane, leading to an over-screened potential in the positive z
direction and an underscreened one in the negative direction.
The potential is also overscreened in the normal plane, as
indicated above by the factor (1 + u20 /4)1/2 .
The equipotential contours of ψ1 are plotted in Fig. 5(a),while the field lines of the resulting electric field, E 1 =−∇ψ1 , are plotted in Fig. 5(b). The figure highlights thedistinct features of the steady-state solution both in the nearand far fields. In the near field, the equipotential contoursof ψ1 resemble the ellipsoidal charge density of Fig. 4(a),and describe an anisotropic screened interaction. The far-field
behavior, on the other hand, can be described by the dipole
moment of the charge density, as discussed before. Therefore,
the electrostatic interaction is long ranged and not screened asin the DH equilibrium case.

(28)

The relaxation effect, therefore, carries over from the charge
density to the electrostatic potential. The steady flow breaks
the spherical symmetry between the z direction and the xy

In this paper, we have presented a general formalism toanalyze the steady-state linear response of bulk electrolyte solutions to an externally imposed source. Our focus was on theeffect of constant relative velocity between the electrolyte andsource. In the absence of such a velocity, the derived response
functions reduce to the Oseen’s result for the hydrodynamic
variables and the DH one for the electrostatic variables.
The steady flow is shown to break the natural sphericalsymmetry of the bulk and to define an azimuthal symmetryaround the direction of motion. Namely, the hydrodynamic
drag in presence of flow distorts the ionic cloud from a sphereto an ellipsoidal shape (relaxation effect) and, consequently,
modifies the electric field close to the source and far from it.
Close to the source, the elongated ionic cloud exerts a net
electric force in the direction of u0 . If we consider a sourceof total charge Q moving with velocity −u0 , the electric force
thus opposes the direction of motion. For small velocities, thisforce scales as ∼Q2 u0 . This result can be obtained by expanding the electrostatic potential of Eq. (27) to linear order in u0 ,and integrating the electric field over the surface of a small
sphere around the origin. Such a force can be interpreted as032604-7

ADAR, UEMATSU, KOMURA, AND ANDELMANPHYSICAL REVIEW E 98, 032604 (2018)

an increase in the source friction coefficient. More generally,
the change in viscosity of charged suspensions is known as
the electroviscous effect [20].
Far from the source, the anisotropic charge density together
with the charged source can be described by an electric dipole
moment −Qu0 , where Q is the total external charge. Therefore far-field behavior of the electric field E 1 and electric
current j 1 are well approximated by dipolar fields. Such longranged fields qualitatively differ from the equilibrium ones,
where the electric field decays exponentially and no current
arises.
The above effect is expected to be pronounced especiallyfor large u0 values. This dimensionless parameter quantifies
the role of convection over diffusion, and can be interpreted
as the Péclet number [3,29], with the Debye length, λD ,playing the role of the characteristic length scale in the system.
The interplay between convection and diffusion is discussed
further below. For the sake of clarity, we return to physicalquantities in their natural units. In particular, we refer to u0
below in units of velocity.
According to Eq. (9), high velocities correspond to u0 >kB T /(ζ λD ) = D/λD , where D is the diffusion coefficient.
When this condition is satisfied, the thermal energy of ionsis not sufficient to overcome the hydrodynamic drag acrossthe Debye length, λD , and ions cannot form a screening ioniccloud. As an alternative interpretation, high velocity values
refer to u0 > λD /τ , where we have introduced the relaxation
time (also called “ambipolar time”), τ = λ2D /D, as the typicaltime for ions to diffuse throughout the Debye length [5].
Therefore, at high velocities, ions move too rapidly to screena charged object properly.
We estimate the velocity scale as mentioned above, byconsidering a dilute aqueous solution at room temperature
with n0 = 2 mM. The corresponding Debye length is λD
10 nm. The diffusion coefficient of simple ions in such asolution is D ≈ 10−9 m2 /s [31], leading to τ 100 ns. Thenthe ratio λD /τ yields a velocity of u0 √0.1 m/s, which isrelatively large. However, as λD /τ ∼ 1/( εζ ), the velocity
becomes much smaller for solvents that are more viscous thanwater even if they are less polar.
As an example for such a viscous solvent, we mention
glycerol with permittivity ε 40ε0 (where ε0 is the vacuum
permittivity) and viscosity η 1.4 Pa s at T = 293 K [32].
For the same ionic concentration of the aqueous solution
above, the velocity u0 100 μm/s is obtained. Velocities oforder 100 μm/s are accessible in experiments and have beenused in colloidal-probe AFM measurements [33–35] as therelative velocity between probe and object.
We note that the relaxation effect within our framework is
continuous in u0 , without any critical behavior. Furthermore,
it is not expected to play a role in standard experiments
designed to measure equilibrium forces. In these setups, it iscustomary to measure forces over a range of velocities, inorder to ensure that the force is velocity independent [36].
In this manner, possible electrohydrodynamic effects are
surpassed.
As a possible experimental setup to capture the dipolar electric interaction, we propose to measure the relativetranslation between two identical, spherical and charged colloids, undergoing sedimentation in a large container. The

hydrodynamic forces between two identical spheres in an
unbound fluid induce no relative translation [23,37], leaving
the mutual electric force as the main possible origin for sucha motion. For micron-sized colloids sedimenting in a diluteelectrolyte (see, e.g., Ref. [38]), the electric force between
colloids can be of the same order of magnitude as the external
gravitational force, making this experimental setup feasible inrelation to our theory.
At this point, we would like to clarify some of the limitations of our framework. First, beyond the linear approximation, quadratic terms such as c∇ψ appear in the Stokes
and force balance equations, mixing these three decoupledequations. Second, our framework was derived for a diluteelectrolyte. Beyond the dilute limit, steric effects between
the ions become important, and can be described by a latticegas model. Including the corresponding pressure term in theformalism yields the modified Poisson-Nernst-Planck framework [39].
In addition, ionic correlations also become important forhigh ionic concentrations. They affect not only the pressure,but also the solvent viscosity [40–42] and the dielectric
constant [43–45]. However, we note that within our linear
framework, such a dependence on the ionic concentration
amounts to a mere shift in the values of the dielectric constantand viscosity, corresponding to the homogeneous value n0 .
Consequently, all the results derived in this work remainunchanged.
As a possible future extension of our model, we recall thatthe independent response to the externally imposed chemical
potential, μ, and charge density, q, holds only for anions andcations with equal friction coefficients, as considered here.
In principle, cations and anions can have different friction
coefficients, ζ+ = ζ− , due to their different sizes and chemical
properties [31]. Such a friction asymmetry couples between
the number density and electric variables. For example, anexternal charge density, q, imposes oppositely directed forceson the cations and anions. Due to their different mobilities,
a number density current is generated with its corresponding number density. The response functions for asymmetricfriction coefficients are derived in Appendix B, and will be
investigated further in a future work.
Another natural extension of our linear-response theory
would be to include boundaries, such as flat surfaces. The
generalization of the Oseen tensor for a flat surface with ano-slip boundary condition is known as the Blake tensor [46].
Our result for the electrostatic potential in Eq. (27) can besimilarly modified for a flat boundary condition by the use of
the “image charge” method. This scenario will be exploredelsewhere.
As was mentioned before, the linear response to moving
sources is relevant for several colloidal systems. Examples
vary from sedimenting colloids, colloidal probe AFM or
manipulation of optically trapped colloids in solution. In
addition, this work can be of use in biological systems, wherecharged molecules interact in aqueous environments underflow. For example, it was shown that charged biomolecules are
only partially screened due to the presence of electrodiffusion
current flow in aqueous pores [47]. We hope that our generalframework can further applied in a wide range of physical andbiological systems.

032604-8

LINEAR RESPONSE FUNCTIONS OF AN ELECTROLYTE …
ACKNOWLEDGMENTS

We thank H. Diamant, M. Urbakh, and M. Bazant for
fruitful discussions and helpful suggestions. R.A. thanks the
hospitality of Tokyo Metropolitan University, where part ofthis research was conducted under the TAU-TMU cotutorial
program. Y.U. was supported by a Grant-in-Aid for the Japan
Society for the Promotion of Science (JSPS) Research Fellow
No. 16J00042. He also thanks the hospitality of Tel Aviv
University, where this work has been completed. S.K. acknowledges support by a Grant-in-Aid for Scientific Research(C) (Grant No. 18K03567) from the JSPS. D.A. thanks the
ISF-NSFC (Israel-China) joint research program under Grant
No. 885/15 for partial support.

PHYSICAL REVIEW E 98, 032604 (2018)

This is a set of linear equations to be solved in Fourier spaceasi
P1 (k) = − 2 k · F(k),
k


kk
1
u1 (k) = 2 I − 2 · F(k),ηkkk2
n1 (k) = − 2μ(k),k + iu0 · k(u0 · k)kμ(k),J 1 (k) = 2
k + iu0 · k
ψ1 (k) =APPENDIX A: LINEARIZATION SCHEME

Consider a small perturbation to the homogenous solution
due to three weak sources: an external force density, F, acharge density, q, and an externally imposed chemical potential, μ. The resulting fields, which have been rescaled in
Eq. (9), can be expanded around the homogeneous solution asu ≈ u0 + u1 ,P ≈ 1 + P1 ,n ≈ 1 + n1 ,J ≈ J 1,c ≈ c1 ,j ≈ j 1,k 2 + iu0 · kq(k),2
0·k+k )k 2 (1 + iu

k2
ψ1 (k),c1 (k) = − 2
k + iu0 · k(u0 · k)k
ψ1 (k).
j 1 (k) = 2
k + iu0 · k(A4)The above solution is written in Eqs. (15), (17), (25), and (26),and the real-space solutions can be obtained by taking theappropriate convolutions.

(A1)

APPENDIX B: ASYMMETRIC FRICTION COEFFICIENTSThe linear corrections, denoted by the subscript “1,” areconsidered to be of a small magnitude, comparable to that of
the sources.
Expanding Eq. (10) to first order leads to the following set
of seven equations corresponding to Eqs. (11)–(13):Generally, cations and anions may have different friction
coefficients, ζ± . Such a friction asymmetry results in a slightly
modified set of equations. The force balance equations for theions in Eq. (6) is now replaced byη∇ u1 (r ) − ∇P1 (r ) = −F(r ),en− ∇ψ − ζ− n− (v − − u) − kB T ∇n− = n− ∇μ.

ψ ≈ ψ1 .

−en+ ∇ψ − ζ+ n+ (v + − u) − kB T ∇n+ = n+ ∇μ,2∇ ψ1 (r ) + c1 (r ) = −q(r ),2Note that the original physical quantities (and not dimensionless ones) are used in the above equation. Following thesame scheme as in Sec. II B, it is possible to define thedimensionless variables for the number density n and charge
density c. In the scaling of the velocities [see Eq. (9)], we use
the average friction coefficient, ζ = (ζ+ + ζ− )/2. In addition,
we define the dimensionless friction asymmetry parameter:

J 1 (r ) + ∇n1 (r ) = −∇μ(r ),j 1 (r ) + ∇c1 (r ) = −∇ψ1 (r ),∇ · u1 (r ) = 0,
∇ · J 1 (r ) + u0 · ∇n1 (r ) = 0,
∇ · j 1 (r ) + u0 · ∇c1 (r ) = 0.

(A2)

Because all the coupling terms are of second order in ourcorrection fields, the above equations are decoupled. The fluid
velocity u1 and pressure P1 depend solely on F. The number
density n1 and corresponding current J 1 depend solely on μ.
The charge density c1 and electric current j 1 depend solely on
q. The potential ψ is related to the charge density source q via
Poisson’s equation.
In Fourier space, Eq. (A2) becomesξ=ζ+ − ζ−.
ζ+ + ζ−(B2)

The above parameter ξ couples between the number density and charge density fields. Explicitly, Eq. (10) is replacedby the following set of seven equations:η∇ 2 u − ∇P − c∇ψ = −F,∇ 2 ψ + c = −q,
J + ξ j + c∇ψ + ∇n = −n∇μ,−ηk 2 u1 (k) − ikP1 (k) = −F(k),j + ξ J + n∇ψ + ∇c = −c∇μ,k 2 ψ1 (k) − c1 (k) = −q(k),∇ · u = 0,J 1 (k) + ikn1 (k) = −ikμ(k),∇ · J + u · ∇n = 0,

j 1 (k) + ikc1 (k) = −ikψ1 (k),∇ · j + u · ∇c = 0.

ik · u1 (k) = 0,
ik · J 1 (k) + iu0 · kn1 (k) = 0,
ik · j 1 (k) + iu0 · kc1 (k) = 0.

(B1)

(A3)

(B3)

While most of the equations remain unchanged, the third andfourth equations contain a new term that is linear in ξ . It isevident that Eq. (10) is restored for ξ = 0.

032604-9

ADAR, UEMATSU, KOMURA, AND ANDELMANPHYSICAL REVIEW E 98, 032604 (2018)

Employing the same linearization scheme, the linear correction terms can be obtained in Fourier space. The hydrodynamic fields, P1 and u1 , are given by the Oseen’s result in
Eq. (15). The number density and charge density fields are
now given by

B(k) =

ψ1 (k) =c1 (k) = ξ B(k) μ(k) + A(k) ψ1 (k),j 1 (k) = −(u0 · k̂) k̂ c1 (k),(B4)

where the functions A(k) and B(k) are defined as
k 2 (k 2 + iu0 · k),A(k) = − 2
k + (iu0 · k)2 + (ξ u0 · k)2[1] M. Doi, Soft Matter Physics (Oxford University Press, Oxford,2013).
[2] S. S. Dukhin and B. V. Derjaguin, in Surface and Colloid
Science, edited by E. Matijrvic (Wiley, New York, 1974),
Vol. 7.
[3] J. H. Masliyah and S. Bhattacharjee, Electrokinetic and Colloid
Transport Phenomena (Wiley, New Jersey, 2006).
[4] D. J. Shaw, Introduction to Colloid and Surface Chemistry(Butterworths, London, 1980).
[5] M. Z. Bazant, K. Thornton, and A. Ajdari, Phys. Rev. E. 70,
021506 (2004).
[6] T. M. Squires and M. Bazant, J. Fluid Mech. 509, 217 (2004).
[7] M. Z. Bazant and T. M. Squires, Phys. Rev. Lett. 92, 066101(2004).
[8] L. Joly, C. Ybert, E. Trizac, and L. Bocquet, Phys. Rev. Lett.
93, 257805 (2004).
[9] L. Joly, C. Ybert, E. Trizac, and L. Bocquet, J. Chem. Phys.
125, 204716 (2006).
[10] A. Ajdari and L. Bocquet, Phys. Rev. Lett. 96, 186102(2006).
[11] Y. Uematsu and T. Araki, J. Chem. Phys. 139, 094901 (2013).
[12] J. L. Anderson, Annu. Rev. Fluid Mech. 21, 61 (1989).
[13] H. J. Keh and H. C. Ma, Langmuir 23, 2879 (2007).
[14] A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T.
Purcel, and L. Bocquet, Nature 494, 455 (2013).
[15] R. W. O’Brien and L. R. White, J. Chem. Soc. Faraday Trans.
II 74, 1607 (1978).
[16] S. R. Maduar, A. V. Belyaev, V. Lobaskin, and O. I.
Vinogradova, Phys. Rev. Lett. 114, 118301 (2015).
[17] K. Kim, Y. Nakayama, and R. Yamamoto, Phys. Rev. Lett. 96,
208302 (2006).
[18] V. Lobaskin, B. Dünweg, M. Medebach, T. Palberg, and C.
Holm, Phys. Rev. Lett. 98, 176105 (2007).
[19] G. Giupponi and I. Pagonabarraga, Phys. Rev. Lett. 106, 248304(2011).
[20] H. Ohshima, Biophysical Chemistry of Biointerfaces (John
Wiley & Sons, New York, 2010).
[21] F. Booth, J. Chem. Phys. 22, 1956 (1954).
[22] C. W. Oseen, Neuere Methoden Und Ergebnisse in Der Hydrodynamik (Akademische Verlagsgesellschaft, Leibzig, 1927).

(B5)

Finally, the electrostatic potential ψ1 is related to q and μ
by the relation

n1 (k) = A(k) μ(k) + ξ B(k) ψ1 (k),J 1 (k) = −(u0 · k̂) k̂ n1 (k),k 2 (iu0 · k).
2
2
0 · k) + (ξ u0 · k)

(k 2 + iu

q(k) + ξ B(k)μ(k).
k 2 [1 − A(k)](B6)

Equations (B4) and (B5) are the generalizations ofEqs. (17) and (26), respectively, for asymmetric friction coefficients, ζ+ = ζ− . In both cases, the currents J and j are
related to n and c, respectively, by a factor of −(u0 · k̂) k̂,
as given by the corresponding continuity equations. The response of the n and J fields to their direct source μ andindirect one ψ1 is found to be the same as the response ofc and j to their direct source ψ1 and indirect one μ.

[23] J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Martinus Nijhoff, The Hague, 1983).
[24] P. Debye and E. Hückel, Phys. Z. 24, 185 (1923).
[25] L. Onsager, Phys. Rev. 38, 2265 (1931).
[26] D. E. Elrick, Australian J. Phys. 15, 283 (1962).
[27] A. Okubo and M. J. Karweit, Limnol. Oceanogr. 14, 514 (1969).
[28] D. Andelman, in Handbook of Physics of Biological Systems,
edited by R. Lipowsky and E. Sackman (Elsevier Science,
Amsterdam, 1995), Vol. I, Chap. 12.
[29] J. Choi, D. Margetis, T. M. Squires, and M. Z. Bazant, J. Fluid
Mech. 536, 155 (2005).
[30] D. A. Lever, J. Fluid Mech. 92, 421 (1979).
[31] P. Vanýsek, in CRC Handbook of Chemistry and Physics, edited
by D. R. Lide, 83rd ed. (CRC Press, Boca Raton, FL, 2002).
[32] G. Åkerlöf, J. Am. Chem. Soc. 54, 4125 (1932); J. B. Segur and
H. E. Oberstar, Indus. Eng. Chem. 43, 2117 (1951).
[33] T. Lai, Y. Meng, W. Zhan, and P. Huang, J. Adhes. 94, 334(2018).
[34] A. Darwiche, F. Ingremeau, Y. Amarouchene, A. Maali, I.
Dufour, and H. Kellay, Phys. Rev. E 87, 062601 (2013).
[35] C. D. F. Honig and W. A. Ducker, Phys. Rev. Lett. 98, 028305(2007).
[36] R. Lhermerout and S. Perkin, Phys. Rev. Fluids 3, 014201(2018).
[37] T. Goldfriend, H. Diamant, and T. A. Witten, Phys. Rev. E 93,
042609 (2016).
[38] J. C. Crocker and D. G. Grier, Phys. Rev. Lett. 77, 1897
(1996).
[39] M. S. Kilic, M. Z. Bazant, and A. Ajdari, Phys. Rev. E 75,
021503 (2007).
[40] G. Jones and M. Dole, J. Am. Chem. Soc. 51, 2950 (1929).
[41] H. Falkenhagen and E. L. Vernon, Z. Phys. 33, 140 (1932).
[42] L. Onsager and R. M. Fuoss, J. Phys. Chem. 36, 2689 (1932).
[43] L. Onsager, J. Am. Chem. Soc. 58, 1486 (1936).
[44] J. G. Kirkwood, J. Chem. Phys. 4, 592 (1936).
[45] J. B. Hasted, D. M. Ritson, and C. H. Collie, J. Chem. Phys. 16,
1 (1948).
[46] J. R. Blake, Proc. Gamb. Phil. Soc. 70, 303 (1971).
[47] Y. Liu, J. Sauer, and R. W. Duton, J. Appl. Phys. 103, 084701(2008).

032604-10