2.1. General framework of the model

Our model describes applications such as D-SFA or CP-AFM. We consider two solid spheres of the same radius $2R$ , physically identical to each other, separated by a small distance $D\ll R$ . In the framework of the Derjaguin approximation (Israelachvili Reference Israelachvili2011) and the related approximations developed below, this geometry is equivalent to a sphere of radius $R$ and a plane. We assume a full symmetry around the axis $Oz$ defined by the sphere centres, with the origin $O$ (a fixed point) located midway between the two spheres, and a plane symmetry across $z=0$ . The cylindrical coordinates are denoted as $(O,r,z)$ , the gap between the surfaces at the radial distance $r$ is $h(r,t)$ . The system defined in this manner is depicted in figure 1. A list of the main variables used in the paper is provided in Appendix A.

The solids are immersed in a solution of monovalent electrolyte with bulk number density $n_0$ (number of cations and anions per unit volume). The solids are assumed to be dielectric and their surface bears a surface charge $\sigma$ . The relative velocity of the spheres is denoted as $\dot D$ . In §§ 4 and 6 we will focus on the case of a harmonic driving at angular frequency $\omega$ and small amplitude $h_0$ around a fixed distance $D$ , corresponding to a relative sphere velocity $\dot D=h_0\omega \text{Re} [e^{i\omega t}]$ .

The quantity of interest is the force $F$ acting on the upper sphere along the $Oz$ axis of symmetry. We consider only contributions to the force arising from electrostatics in the mean-field approximation, hydrodynamics and electrokinetics. We refer to § 7 for a discussion of other effects and range of validity of our model.

Figure 1.

Schematics of the system geometry: an electrolytic solution is confined between two identical spheres of curvature radius $2R$ and whose separation $D(t)$ varies in time.

Our model uses the standard formulation of the EK theory in the regime of slowly varying fields, i.e. it neglects the magnetic field and wave propagation. The transport of ions is governed by the Nernst–Planck equation. The number density of cations and anions is denoted as $n^+$ and $n^-$ , their diffusion coefficient $D^+$ and $D^-$ , their electrochemical potential $\mu ^+=k_{\textit{B}}T \ln (n^+/n_0)+eV$ and $\mu ^-=k_{\textit{B}}T \ln (n^-/n_0)-eV$ , where $T$ is the temperature, $k_{\textit{B}}$ the Boltzmann constant, $V$ the electric potential and $e$ the elementary charge.

In D-SFA ( $R\sim 3\,\textrm {mm}$ ) and CP-AFM ( $R\sim 10$ ${\unicode[Arial]{x03BC}}$ m) experiments, the Reynolds number is usually very low. An upper limit is $\textrm{Re} \leq \rho \textit{v}_r D/\eta$ with $\rho$ the liquid density, $\eta$ the liquid viscosity and $\textit{v}_r= \dot D\sqrt {R/D}$ a maximum radial velocity. The driving rate in local probe measurements is usually at most $\dot D/D \lt 0.1\,$ s $^{-1}$ , and the working distance $D \lt R/100$ . With the density and viscosity of water, the Reynolds number does not exceed $10^{-3}$ . The flow will thus be described by the incompressible Stokes equation. We introduce the hydrodynamic stress tensor and the electrostatic Maxwell stress tensor:

(2.1) \begin{eqnarray} \unicode{x1D64E}^{\mathcal H} = -P\unicode{x1D644} + \frac {\eta }{2} ({\boldsymbol{\nabla }}\boldsymbol{v}+^t{\boldsymbol{\nabla }} \boldsymbol{v}) \quad \vert z \vert \lt h/2 , \end{eqnarray}

(2.2) \begin{eqnarray} \unicode{x1D64E}^{\mathcal M} = \varepsilon _i\left ({-}\frac {\Vert \boldsymbol{\nabla } V \Vert ^2}{2}\unicode{x1D644}+\boldsymbol{\nabla } V \otimes \boldsymbol{\nabla } V \right ) \quad \begin{cases} \varepsilon _i=\varepsilon _{\textit{l}} \ \ \textrm {if} \ \vert z \vert \lt h/2,\\ \varepsilon _i=\varepsilon _{\textit{s}} \ \ \textrm {if} \ \vert z \vert \gt h/2, \end{cases} \end{eqnarray}

where $\boldsymbol{v}=\textit{v}_r \boldsymbol{e}_r+\textit{v}_z\boldsymbol{e}_z$ is the flow velocity, $\unicode{x1D644}$ the identity matrix, $P$ the pressure in excess of the atmospheric pressure, $\eta$ the liquid viscosity which will be assumed independent of the ion concentrations and $\varepsilon _{\textit{l}}$ , $\varepsilon _{\textit{s}}$ the dielectric permittivities of the liquid and of the solids, respectively. With these notations, the Poisson, Nernst–Planck and Stokes equations, as well as the conservation laws for the ions and the incompressible fluid, are

(2.3a) \begin{align} \boldsymbol{\nabla} ^2 V = - \frac {e (n^+-n^-)}{\varepsilon _{\textit{l}}} \quad \vert z \vert \lt h/2 \quad \boldsymbol{\nabla} ^2 V =0 \quad \vert z \vert \gt h/2 , \end{align}

(2.3b) \begin{align} \boldsymbol{j}^{\pm }=-\frac {D^{\pm }}{k_{\textit{B}}T} n^{\pm } \boldsymbol{\nabla } \mu ^{\pm } + n^{\pm } \boldsymbol{v} , \end{align}

(2.3c) \begin{align} \textrm {div} (\unicode{x1D64E} ^{\mathcal M}+\unicode{x1D64E} ^{\mathcal H}) = \boldsymbol{0} = -\boldsymbol{\nabla }\! P + \eta \boldsymbol{\nabla} ^2 \boldsymbol{v} - e (n^+-n^-) \boldsymbol{\nabla } V , \end{align}

(2.3d) \begin{align} \frac {\partial n^{\pm }}{\partial t} =-\textrm {div} \boldsymbol{j}^{\pm } \quad \textrm {div} \boldsymbol{v} = 0 . \end{align}

The boundary conditions at the symmetry plane $z=0$ (quantities at $z=0$ are hereafter denoted with the subscript m), correspond to invariance by the transformation $z \rightarrow -z$ :

(2.4) \begin{equation} z=0, \quad \textit{v}_{z,\textit{m}}=0 , \quad \frac {\partial V}{\partial z}\bigg |_{\textit{m}}=0 , \quad \frac {\partial \mu ^{\pm }}{\partial z}\bigg |_{\textit{m}}=0 , \quad \frac {\partial \textit{v}_r}{\partial z}\bigg |_{\textit{m}} = 0 . \end{equation}

At the liquid–solid interface $z=h/2$ we will use the continuity of the electric potential and the Maxwell–Gauss equation for the normal component of the electrical displacement. We will also assume a no-slip boundary condition. Quantities at $z=h/2$ are denoted with the subscript $w$ and satisfy

(2.5a) \begin{align} V_{{w}^-}=V_{{w}^+}, \end{align}

(2.5b) \begin{align} - \varepsilon _{\textit{s}} \boldsymbol{\nabla } V \vert _{{w}^+} \boldsymbol{\cdot }\boldsymbol{n} + \varepsilon _{\textit{l}} \boldsymbol{\nabla } V \vert _{{w}^-}\boldsymbol{\cdot }\boldsymbol{n}= \sigma ,\end{align}

(2.5c) \begin{align} \boldsymbol{v}_{{w}}= \frac {\dot D}{2} \boldsymbol{e}_z,\end{align}

(2.5d) \begin{align} (\boldsymbol{j}^{\pm }_{{w}}-n^{\pm } \boldsymbol{v}_{{w}})\boldsymbol{\cdot }\boldsymbol{n}=j^{\pm ,\textit{source}}. \end{align}

Here $\boldsymbol{n}$ denotes the unit normal to the upper sphere surface, pointing from the liquid into the upper sphere. The superscripts ${w}^-$ and ${w}^+$ refer, respectively, to the liquid and solid sides of the interface. The source terms $j^{\pm ,\textit{source}}$ account for the charge regulation phenomena at a dielectric surface when the surface charge is not fixed, and describe the associated adsorption/emission of ions.

Without loss of generality, we take $\sigma \gt 0$ in the following. We also consider that $V(r=\infty )=0$ , i.e. no electric field applied in the reservoir.

The total stress is continuous through the interface. We denote $\boldsymbol{F}$ the force that an external operator has to apply in order to maintain the prescribed driving velocity $\dot {D}$ of the upper sphere, and $F=\boldsymbol{F} \boldsymbol{\cdot }\boldsymbol{e}_z$ its component on the $z$ -axis. The force balance on the upper sphere projected on the $z$ -axis gives

(2.6) \begin{equation} F = \boldsymbol{e}_z \boldsymbol{\cdot }\left ( \oint _{S} \left (\boldsymbol{S}^\mathcal{H} + \boldsymbol{S}^\mathcal{M}\right )\big |_{{w}^-} \boldsymbol{\cdot }({-}\boldsymbol{n}) \,\mathrm{d}S -\oint _{S} \boldsymbol{S}^\mathcal{M}\big |_{{w}^+} \boldsymbol{\cdot }({-}\boldsymbol{n}) \,\mathrm{d}S \right ) \end{equation}

where $\boldsymbol{n}$ is the unit normal defined above (pointing from the liquid into the upper sphere). The first integral is the traction exerted by the liquid on the sphere and the second the traction exerted by the Maxwell stress in the dielectric solid.

Our model calculates the force response within two further approximations:

1. (i) the Derjaguin approximation (Israelachvili Reference Israelachvili2011) and its consequences discussed in the next subsection;

2. (ii) the approximation of local thermodynamic equilibrium (LTE) across the thickness $h$ of the liquid film, developed in § 3.

2.2. The Derjaguin approximation

The Derjaguin approximation assumes that the interacting surfaces are nearly parallel. Therefore, the unit normal vector to the surfaces in the interaction region is well approximated by the axial direction: $\boldsymbol{n} \approx \boldsymbol{e}_z$ . While the condition for the approximation is well established and widely used for equilibrium interactions, its application to electrohydrodynamic interactions requires more careful consideration.

Equilibrium interactions. For electrostatic interactions in an electrolyte, the characteristic decay length is the Debye length $\ell _{\textit{D}}$ . In the sphere–plane geometry with minimal gap $D$ , electrostatic interactions become negligible when the local gap exceeds a few Debye lengths. Using the parabolic approximation for the gap,

(2.7) \begin{equation} h (r, t) = D(t) + \frac {r^2}{2R} \end{equation}

the lateral extent of the electrostatic interaction region is $r \sim \sqrt {R\ell _{\textit{D}}}$ . The maximum slope of the surfaces over this lateral distance is

(2.8) \begin{equation} \left |\frac {\textrm{d}h}{\textrm{d}r}\right |_{r \sim \sqrt {R\ell _{\textit{D}}}} \sim \frac {\sqrt {R\ell _{\textit{D}}}}{R} = \sqrt {\frac {\ell _{\textit{D}}}{R}} \ll 1 .\end{equation}

Hydrodynamic interactions. The hypothesis of nearly parallel surfaces corresponds to the condition for the lubrication approximation of the Stokes equation. There is no intrinsic characteristic decay length for the hydrodynamic interaction as the pressure field in a lubrication flow between curved surfaces decays algebraically with the separation distance $D$ . However, Brenner (Reference Brenner1961) showed that the hydrodynamic pressure becomes negligible beyond a lateral distance of order $\sqrt {RD}$ from the apex. This defines the lateral extent of the hydrodynamic interaction region. The maximal slope of the surfaces at this lateral distance is

(2.9) \begin{equation} \left |\frac {\textrm{d}h}{\textrm{d}r}\right |_{r \sim \sqrt {RD}} \sim \frac {\sqrt {RD}}{R} = \sqrt {\frac {D}{R}} \ll 1 .\end{equation}

Physical consequences. Defining $\lambda = \max (\ell _{\textit{D}}, D)$ as the maximum interaction range in the normal direction, the condition $\lambda \ll R$ ensures a separation of length scales between normal variations (over scale $\sim \lambda$ ) and lateral variations (over scale $\sim \sqrt {R\lambda }$ ). For any quantity $a$ that varies significantly across the film thickness the following inequality holds:

(2.10) \begin{equation} \left \vert \frac {\partial a}{\partial r} \right \vert \ll \left \vert \frac {\partial a}{\partial z} \right \vert \quad \vert z \vert \le \frac {h}{2} .\end{equation}

An important consequence is that the Poisson equation in the liquid ((2.3a ), $\vert z \vert \lt h/2$ ) and its boundary condition (2.5b ) reduce to their one-dimensional form:

(2.11) \begin{equation} \frac {\partial ^2 V}{\partial z^2} = -\frac {e (n^+-n^-)}{\varepsilon _{\textit{l}}} \quad \left .\frac {\partial V}{\partial z}\right \vert _{{w}^-} = \frac {\sigma }{\varepsilon _{\textit{l}}} .\end{equation}

It is shown in Appendix B that the normal electric field in the dielectric solids is of the same order as the radial field $(\partial V/\partial r) (h^-/2)$ , which is negligible according to inequality (2.10). By integrating equation (2.11) along $z$ the local electroneutrality condition is obtained:

(2.12) \begin{equation} \int _0^{h (r,t) /2}e (n^+-n^-) (r,z,t) \textrm {d}z=-\sigma .\end{equation}

With $\boldsymbol{n} \approx \boldsymbol{e}_z$ and inequality (2.10), the force per unit surface due to the Maxwell stress tensor at the interface is

(2.13) \begin{align} \unicode{x1D64E}^{\mathcal M}_{{w}} \boldsymbol{n} \approx S^{\mathcal M}_{zz} \boldsymbol{e}_z \approx \frac {\varepsilon _i}{2}\left (\frac {\partial V}{\partial z} \right ) ^2 \boldsymbol{e}_z .\end{align}

This contribution is negligible in the dielectric solids (see Appendix B). With $\boldsymbol{n}\approx \boldsymbol{e}_z$ and $\mathrm{d}S \approx 2\pi r\,\mathrm{d}r$ , the force (2.6) reduces to

(2.14) \begin{align} F & = - \int _0^{\infty }2\pi r (S ^{\mathcal M}_{zz}+S ^{\mathcal H}_{zz}) (r,h^-/2) \textrm {d}r \nonumber \\ &= - \int _0^{\infty }2\pi r \left [-P+\eta \frac {\partial \textit{v}_z}{\partial z} +\frac {\varepsilon _{\textit{l}}}{2}\left (\frac {\partial V}{\partial z} \right ) ^2 \right ] \left (r,\frac {h^-}{2}\right ) \textrm {d}r. \end{align}

We show in Appendix B that within our approximations the quantity in brackets is uniform across the film thickness and that the force reduces to

(2.15) \begin{equation} F=\int _0^\infty 2\pi r P (r,0) \textrm {d}r .\end{equation}

3.1. Conditions for time scale separation

In addition to the separation of length scales $h_0\ll D \ll R$ , we assume that the diffusion time $\tau _z$ across the thickness of the electrolyte film is much shorter than the time scale of the driving:

(3.1) \begin{equation} \tau _z = \frac {D^2}{D^\pm } \ll \tau =\frac {D}{\dot D} \quad \textrm {or} \quad \tau =\frac {2\pi }{\omega } .\end{equation}

For a continuous-like driving at velocity $\dot D$ , the time scale separation is equivalent to the condition of small normal Péclet number. For an oscillatory driving at amplitude $h_0\ll D$ , the small-normal-Péclet condition is not sufficient and the time scale separation requires a condition on the driving frequency.

In these conditions of spatial and time scales separations the transport equation (2.3d ) reduces to $\partial _z\! (n^{\pm }\partial _z \mu ^{\pm } )=0$ to the leading order in $\tau _z/\tau$ (see also Appendix B.3). Together with the boundary condition at $z=0$ this gives $\partial \mu ^{\pm }/\partial z=0$ to the leading order in $\tau _z/\tau$ , which is the approximation of LTE: $\mu ^{\pm }$ are uniform across the film thickness. This relies only on the midplane symmetry (2.4) and not on the condition at the walls: the possible ion source in (2.5d ) feeds only the thickness-integrated balances, so the result holds whether or not surface sources are present.

The range of driving rate, or driving frequency, for the time scale separation depends on the working range $D$ and on the diffusion coefficient of ions. Taking the case of a KCl electrolyte of average diffusion coefficient 2.0 $\times 10^{-9}\textrm {m}^2$ s⁻¹ in water at 300 K (Haynes Reference Haynes2016), typical values are as follows.

The D-SFA. With a typical working distance $D\lesssim 1\,$ µm (Garcia et al. Reference Garcia, Barraud, Picard, Giraud, Charlaix and Cross2016), the characteristic time is $\tau _z\lesssim 10^{-3}\,\textrm {s}$ . For working frequencies typically spanning between $5\,\textrm {Hz}$ and $500\,\textrm {Hz}$ , the condition $\tau _z \ll 2\pi /\omega$ is well satisfied.

The CP-AFM. With a working distance $D\lesssim 100\,\textrm {nm}$ (Liu et al. Reference Liu, Klaassen, Zhao, Mugele and van den Ende2018; Rodríguez Matus et al. Reference Rodríguez Matus, Zhang, Benrahla, Majee, Maali and Würger2022) the characteristic time is $\tau _z\lesssim 10^{-5}\,\textrm {s}$ . The hypothesis of local equilibrium across the film is appropriate for driving rates or working frequencies up to ${\sim} 100\,\textrm {kHz}$ .

3.2. Thermodynamic potentials

We introduce the nominal Debye length $\ell _{\textit{D}}$ associated with the electrolyte concentration $n_0$ in the reservoir, the nominal Gouy length $\ell _{\textit{G}}$ associated with the surface charge $\sigma _0$ of non-interacting surfaces at rest, the Bjerrum length $\ell _{\textit{B}}$ , and the dimensionless parameter $\alpha$ (Saville Reference Saville1977; Israelachvili Reference Israelachvili2011; Adar & Andelman Reference Adar and Andelman2018):

(3.2) \begin{equation} \ell _{\textit{B}}=\frac {e^2}{4\pi \varepsilon _{\textit{l}} k_{\textit{B}}T}, \quad \ell _{\textit{D}}=\sqrt {\frac {\varepsilon _{\textit{l}} k_{\textit{B}} T}{2e^2 n_0}}=\sqrt {\frac {1}{8\pi \ell _{\textit{B}} n_0}}, \quad \frac {2}{\ell _{\textit{G}}}=\frac {\sigma _0 e}{\varepsilon _{\textit{l}} k_{\textit{B}}T}, \quad \alpha =\frac {2\ell _{\textit{D}}}{\ell _{\textit{G}}} .\end{equation}

In order to decouple electrical from chemical effects in the transport equations, sum and difference potentials are introduced as

(3.3) \begin{equation} \mu _{\textit{s}}=\frac {\mu ^++\mu ^-}{2}=k_{\textit{B}}T\ln \frac {n (r,t)}{n_0}, \quad W=\frac {\mu ^+-\mu ^-}{2e}, \quad \chi (r,z,t) =\frac {e (V-W)}{k_{\textit{B}}T} .\end{equation}

The average chemical potential $\mu _{\textit{s}}$ and electrical potential $W$ are uniform along $z$ . The quantity $n (r,t) =n_0 e^{\mu _{\textit{s}}/k_{\textit{B}}T}$ is the equivalent bulk concentration at position $r$ and time $t$ . Using $n^{\pm }=n_0\exp [(\mu ^{\pm }\mp eV)/k_{\textit{B}}T]$ the sum and difference ion densities are derived:

(3.4) \begin{equation} n_{\textit{s}}=n^++n^-=2n (r,t) \cosh \chi, \quad n_{\textit{e}}=n^+-n^-=-2n (r,t) \sinh \chi. \end{equation}

The Poisson equation (2.3a ) and its boundary condition (2.5a ) become

(3.5a) \begin{align} \frac {\partial ^2\chi }{\partial z^2}=\frac {\sinh \chi }{\lambda _{\textit{D}}^2}, \quad \lambda _{\textit{D}} = \sqrt {\frac {\varepsilon _{\textit{l}} k_{\textit{B}} T}{2e^2 n(r,t)}}=\sqrt {\frac {n_0}{n(r,t)}}\ell _{\textit{D}}, \end{align}

(3.5b) \begin{align} \left .\frac {\partial \chi }{\partial z}\right \vert _{h/2(r,t)}= \frac {2}{\lambda _{\textit{G}}}, \quad \frac {2}{\lambda _{\textit{G}}}=\frac {\sigma (r,t) e}{\varepsilon _{\textit{l}} k_{\textit{B}}T} .\end{align}

Thus, $\chi (z)$ satisfies a Poisson–Boltzmann equation along $z$ with local parameters. The analytical solution $\chi$ of (3.5a ) for a positive surface charge is (Ninham & Parsegian Reference Ninham and Parsegian1971)

(3.6a) \begin{align} \chi (r,z,t) =\chi _{\textit{m}} (r,t) -\ln [\textrm {cd}(z/2 \lambda _{\textit{D}} \sqrt {k} \vert k^2) ]^2 \quad k (r,t) =e^{-\chi _{\textit{m}} (r,t)}, \end{align}

(3.6b) \begin{align} \cosh \chi _{{w}} (r,t) -\cosh \chi _{\textit{m}} (r,t)= \frac {2 \lambda _{\textit{D}}^2}{\lambda _{\textit{G}}^2}= \frac {\alpha ^2}{2}\frac {n_0}{n(r,t)} \frac {\sigma (r,t)^2}{\sigma _0^2}, \end{align}

where $\textrm {cd} (u \vert k^2)$ is the Jacobi elliptical function of argument $u$ and parameter $k^2$ , and $k(h,n,\sigma )$ is determined by the boundary condition (3.5b ) (see Appendix D).

At a large distance from the sphere–plane axis the EDLs do not overlap, the midplane potential vanishes, and $k=1$ . In the following we denote $\psi =\chi _{\textit{eq}}$ the equilibrium potential $\chi$ when the surfaces are not moving ( $n\equiv n_0$ everywhere). In our hydrodynamic condition of no-slip at the solid surfaces, the equilibrium surface potential of a non-interacting surface reduces to its nominal zeta potential, and satisfies the Grahame (Reference Grahame1953) equation:

(3.7) \begin{equation} \tilde \zeta = \frac {e\zeta }{k_{\textit{B}}T}=\psi _{w,\textit{0}}, \quad \cosh \tilde \zeta = \frac {\alpha ^2}{2}+1, \quad \sinh \frac {\tilde \zeta }{2}=\frac {\alpha }{2}. \end{equation}

3.3. Transport equations

The transport equations are written in terms of the sum and difference variables by introducing the average diffusion coefficient and the diffusion coefficient contrast:

(3.8) \begin{equation} D_{\textit{s}}=\frac {D^++D^-}{2}, \quad \delta = \frac {D^+-D^-}{D^++D^-}, \quad D^{\pm }=D_{\textit{s}} (1\pm \delta ). \end{equation}

Similarly we define the sum and difference number flow densities as

(3.9) \begin{equation} \boldsymbol{j}_{\!\textit{s}}=\boldsymbol{j}^++\boldsymbol{j}^- ,\quad \boldsymbol{j}_{\!\textit{e}}=\boldsymbol{j}^+-\boldsymbol{j}^- .\end{equation}

Using (3.3) and (3.4) we get the identity

(3.10) \begin{equation} e (n^+-n^-) \boldsymbol{\nabla } V=en_{\textit{e}} \boldsymbol{\nabla } W +n_{\textit{s}}\boldsymbol{\nabla } \mu _{\textit{s}}-k_{\textit{B}}T \boldsymbol{\nabla } n_{\textit{s}} \end{equation}

and write (2.3d ), (2.3b ) and (2.3c ) as

(3.11a) \begin{align} \frac {\partial n_{\textit{s}}}{\partial t}= - \textrm {div} \boldsymbol{j}_{\!\textit{s}}, \quad \boldsymbol{j}_{\!\textit{s}}=n_{\textit{s}} \boldsymbol{v} - \frac {D_{\textit{s}}}{k_{\textit{B}}T}\left \{(n_{\textit{s}}+\delta n_{\textit{e}}) \frac {\textrm{d}\mu _{\textit{s}}}{\textrm{d}r} + e (n_{\textit{e}}+\delta n_{\textit{s}}) \frac {\textrm{d}W}{\textrm{d}r} \right \}\boldsymbol{e}_r, \end{align}

(3.11b) \begin{align} \frac {\partial n_{\textit{e}}}{\partial t}= - \textrm {div} \boldsymbol{j}_{\!\textit{e}}, \quad \boldsymbol{j}_{\!\textit{e}}=n_{\textit{e}} \boldsymbol{v} - \frac {D_{\textit{s}}}{k_{\textit{B}}T}\left \{(n_{\textit{e}}+\delta n_{\textit{s}}) \frac {\textrm{d}\mu _{\textit{s}}}{\textrm{d}r} + e (n_{\textit{s}}+\delta n_{\textit{e}}) \frac {\textrm{d}W}{\textrm{d}r} \right \}\boldsymbol{e}_r, \end{align}

(3.11c) \begin{align} \textrm {div} \boldsymbol{v} = 0, \quad \eta \frac {\partial ^2 \textit{v}_r}{\partial z^2} = \frac {\textrm{d} \varPi }{\textrm{d} r} + en_{\textit{e}} \frac {\textrm{d} W}{\textrm{d} r}+ (n_{\textit{s}}-2n) \frac {\textrm{d} \mu _{\textit{s}}}{\textrm{d} r}, \end{align}

(3.11d) \begin{align} \varPi = P- (n_{\textit{s}}-2n) k_{\textit{B}}T. \end{align}

In the Stokes equation (2.3c ) the radial component $\varDelta _r \boldsymbol{v}$ has been neglected according to the Derjaguin approximation, to obtain (3.11c ). We show in Appendix B that the osmotically compensated pressure $\varPi$ defined in (3.11d ) is uniform across the film thickness.

The three terms on the right-hand side of the Stokes equation (3.11c ) represent three distinct driving forces for the flow. The first term, $\textrm{d}\varPi /\textrm{d}r$ , is uniform across the film cross-section and induces a Poiseuille-like velocity profile: this is the pressure-driven flow. The second term, $en_{\textit{e}} \textrm{d}W/\textrm{d}r$ , represents the electrical driving force, where $-\textrm{d}W/\textrm{d}r$ is the radial electric field that develops in the liquid film. This term vanishes outside the EDLs (where $n_{\textit{e}}=0$ ) and generates an electro–osmotic flow. The third term, $(n_{\textit{s}}-2n) \textrm{d}\mu _{\textit{s}}/\textrm{d}r$ , is the driving force due to the chemical potential gradient. This term also vanishes outside the EDLs and induces a diffusio–osmotic flow (Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984).

3.4. Conservation laws integrated over the thickness

The above equations are integrated across the thickness by defining

(3.12) \begin{eqnarray} N_{\textit{s}} (r,t) =\int _0^{h (r,t) /2} (n_{\textit{s}}-2n_0) \textrm {d}z, \quad N_{\textit{e}} (r,t) =\int _0^{h (r,t) /2} n_{\textit{e}} \textrm {d}z = -\frac {\sigma }{e}, \end{eqnarray}

(3.13) \begin{eqnarray} J_{\textit{s,e}}= \int _0^{h (r,t) /2} \boldsymbol{j}_{\textit{s,e}}\boldsymbol{\cdot }\boldsymbol{e}_r \textrm {d}z, \quad J_{\textit{v}}=\int _0^{h (r,t) /2}\textit{v}_r (z) \textrm {d}z. \end{eqnarray}

The value of $N_{\textit{e}}$ reflects the local electroneutrality (2.12).

In view of (3.11) it is clear that the fluxes $(J_{\textit{v}},J_{\textit{s}},J_{\textit{e}})$ are linear functions of the gradients of the thermodynamic potentials $(\varPi ,\mu _{\textit{s}},W)$ . We show in Appendix C that these fluxes are related to the above gradients by the $3\times3$ symmetric matrix $\unicode{x1D64F}$ as

(3.14a) \begin{align} \begin{pmatrix} J_{\textit{v}}\\[3pt] J_{\textit{s}}-2n (r,t) J_{\textit{v}} \\[3pt] eJ_{\textit{e}}\end{pmatrix} = -\unicode{x1D64F} \frac {\partial }{\partial r} \begin{pmatrix} \varPi \\[3pt] \mu _{\textit{s}} \\[3pt] W \end{pmatrix}, \quad \unicode{x1D64F} =\frac {\boldsymbol{\varGamma }}{\eta } +\frac {D_{\textit{s}}\unicode{x1D647}}{k_{\textit{B}}T}, \end{align}

(3.14b) \begin{align} \varGamma _{ij}= \int _0^{h (r,t) /2} \varGamma _i\varGamma _j \textrm {d}z \quad (i,j) \in (\mathrm{v},\mathrm{s},\mathrm{e}) ,\end{align}

(3.14c) \begin{align} \varGamma _{\textit{v}} (z) = z ,\quad \varGamma _{\textit{s}} (z) = \int _0^z \left (n_{\textit{s}} -2n\right ) \textrm {d}z' \quad \varGamma _{\textit{e}} (z) = \int _0^z e n_{\textit{e}} \textrm {d}z' , \end{align}

(3.14d) \begin{align} L_{\textit{v}\textit{i}}=L_{i\textit{v}}=0 , \end{align}

(3.14e) \begin{align} L_{\textit{ss}}=\frac {L_{\textit{ee}}}{e^2}=N_{\textit{s}}+n_0 h+ N_{\textit{e}} \delta, \quad L_{\textit{se}}=L_{\textit{es}}=e\left [N_{\textit{e}}+ (N_{\textit{s}}+n_0 h) \delta \right ]. \end{align}

The two components $\boldsymbol{\varGamma }$ and $\unicode{x1D647}$ of the transport matrix $\unicode{x1D64F}$ reflect the convective and the diffusive contributions in the ion flux expressions (3.11a ) and (3.11b ). Due to the Onsager reciprocity principle and to the second principle of thermodynamics, the matrix $\unicode{x1D64F}$ is symmetric and positive. The diagonal coefficients are always positive, as well as the determinant of the matrix and its partial diagonal determinants. This can be expressed by inequalities involving the three factors of merit $\xi _{ij}$ ,

(3.15) \begin{eqnarray} -1 \le \xi _{ij}=\frac {T_{ij}}{\sqrt {T_{ii}T_{jj}}} \le 1 , \end{eqnarray}

(3.16) \begin{eqnarray} \frac {\textrm {Det} \unicode{x1D64F}}{T_{\textit{vv}}T_{\textit{ss}}T_{\textit{ee}}}=1-\xi _{\textit{ve}}^2-\xi _{\textit{es}}^2-\xi _{\textit{vs}}^2+2\xi _{\textit{ve}}\xi _{\textit{vs}}\xi _{\textit{se}} \ge 0 . \end{eqnarray}

The three conservation laws (3.11) integrated over the thickness can be written as

(3.17a) \begin{align} \frac {\partial (h/2)}{\partial t}=\frac {\dot D}{2}= - \frac {1}{r} \frac {\partial r J_{\textit{v}}}{\partial r} , \end{align}

(3.17b) \begin{align} \frac {\partial N_{\textit{e}}}{\partial t}=\int _0^{h (t) /2}\left (\frac {\partial n_{\textit{e}}}{\partial t} +\frac {\partial \boldsymbol{j}_{\!\textit{e}}.\boldsymbol{e}_z}{\partial z}\right ) \textrm {d}z = - \frac {1}{r} \frac {\partial r J_{\textit{e}}}{\partial r} +\frac {\partial N_{\textit{e}}^{\textit{source}}}{\partial t} , \end{align}

(3.17c) \begin{align} \frac {\partial N_{\textit{s}}}{\partial t}= - \frac {1}{r} \frac {\partial r (J_{\textit{s}}-2nJ_{\textit{v}})}{\partial r} -\frac {2}{r} \frac {\partial r (n-n_0) J_{\textit{v}}}{\partial r} +\frac {\partial N_{\textit{s}}^{\textit{source}}}{\partial t} , \end{align}

where the source terms depend on the electrical properties of the solid–electrolyte surfaces. The source terms in (3.17b ) and (3.17c ) account for the ions received from the wall per unit area and per unit time, respectively, in total number (cations plus anions) and excess charge (cations minus anions). They are derived from the surface fluxes defined in (2.5d ) by the sum-and-difference combinations:

(3.18) \begin{equation} \frac {\partial N_{\textit{s}}^{\textit{source}}}{\partial t} = -(j^{\textit{+,source}} + j^{\textit{-,source}}) \quad \frac {\partial N_{\textit{e}}^{\textit{source}}}{\partial t}= -(j^{\textit{+,source}} - j^{\textit{-,source}}) .\end{equation}

Typical electrical conditions at the solid/electrolyte interface are as follows.

1. (i) Metallic surfaces, for which the surface charge can change without ion exchange with the liquid, thanks to metal conduction. Perfect metallic surfaces involve no source terms, a constant potential at the surface, and a vanishing electric field $\textrm{d}W/\textrm{d}r=0$ . We will not consider further the case of metallic surface in this work.

2. (ii) Dielectric surface of fixed surface charge, for which the source terms also vanish. Due to the global electroneutrality (3.12), (3.17b ) integrates in $J_{\textit{e}}=0$ .

3. (iii) Dielectric surfaces with charge regulation phenomena. Any change $\textrm{d}\sigma$ in the surface charge results from exchange of ions with the electrolyte, thus $\textrm{d}\sigma =e(\textrm{d}N^{-source}-\textrm{d}N^{+source})=-e\textrm{d}N_{\textit{e}}^{\textit{source}}$ . As a result (3.17b ) also integrates in $J_{\textit{e}}=0$ . The source term $ \partial N_{\textit{s}}^{\textit{source}}/\partial t$ depends on the particular charge regulation mechanism at play, its dependency on the local potential and chemical potentials, as well as its dynamics.

Equations (3.17) are further reduced by noting that

1. (i) $\partial h/\partial t=\dot D$ does not depend on $r$ therefore (3.17a ) can be integrated once, leading to $J_{\textit{v}}=- r\dot D/4$ ;

2. (ii) in the parabolic approximation $r^2=2R(h-D)$ and at all time $r$ d $r=R$ d $h$ . At each time $t$ , the variable $r\in [0,\infty [$ can be changed in the variable $h\in [D(t),\infty [$ .

Using the relations (3.14a ) and eliminating $\textrm{d}W/\textrm{d}r$ , (3.17a ) and (3.17c ) become a set of two second-order nonlinear differential equations,

(3.19a) \begin{align} &\left (T_{\textit{vv}}-\frac {T_{\textit{ve}}^2}{T_{\textit{ee}}} \right ) \frac {\partial \varPi }{\partial h} +\left (T_{\textit{vs}}-\frac {T_{\textit{ve}}T_{\textit{se}}}{T_{\textit{ee}}} \right ) \frac {\partial \mu _{\textit{s}}}{\partial h} =\frac {R\dot D}{4} , \end{align}

(3.19b) \begin{align} &\frac {\partial }{\partial h} \left \{(h-D)\left [ \left (T_{\textit{vs}}-\frac {T_{\textit{ve}}T_{\textit{se}}}{T_{\textit{ee}}} \right ) \frac {\partial \varPi }{\partial h} +\left (T_{\textit{ss}}-\frac {T_{\textit{se}}^2}{T_{\textit{ee}}} \right ) \frac {\partial \mu _{\textit{s}}}{\partial h} +n_0 (e^{\mu _{\textit{s}}/k_{\textit{B}}T}-1)\frac {R\dot D}{2} \right ] \right \} \nonumber \\ & \quad =\frac {R}{2} \left ( \frac {\partial N_{\textit{s}}}{\partial t} -\frac {\partial N_{\textit{s}}^{\textit{source}}}{\partial t} \right )\!, \end{align}

to be solved at each time for $h\in [D(t),\infty [$ , with boundary conditions of vanishing $\varPi$ and $\mu _{\textit{s}}$ at $h=\infty$ and finite second derivatives at $h=D$ . The source term has to be specified by the charge regulation conditions and dynamics. The transport coefficients $T_{ij}$ and the ion excess $N_{\textit{s}}$ are functions of $n=n_0 e^{\mu _{\textit{s}}/k_{\textit{B}}T}$ , $h$ and $\sigma (h,t)$ , given by (3.12) and (3.14). Note that the coefficient $T_{\textit{ee}}$ indeed represents the conductance of the electrolyte film and cannot vanish.

Once the system (3.19) solved, the pressure $\varPi$ is determined from (3.11d ), and the electrohydrodynamic force (3.20) is

(3.20) \begin{equation} F (t) =2\pi R\int _{D(t)}^\infty \left [\varPi +2n k_{\textit{B}}T(\cosh \chi _{\textit{m}} -1)\right ]\textrm {d}h .\end{equation}

3.5. Contributions to the electrohydrodynamic force

The force $F(t)$ can be split in different contributions according to their physical origin.

1. (i) Equilibrium DLVO force. According to (3.19), at equilibrium $\varPi \equiv 0$ , $n \equiv n_0$ and $\chi =\chi _{\textit{eq}}=\psi$ . The equilibrium pressure in the midplane (3.11d ) is the electrostatic disjoining pressure: $P_{\textit{eq}}(r,0)=\pi _{\textit{d,e}}(h(r))=2n_0k_{\textit{B}}T(\cosh \psi _{\textit{m}}-1)$ . In the DLVO theory (Derjaguin & Landau Reference Derjaguin and Landau1941; Verwey & Overbeek Reference Verwey and Overbeek1948) the van der Waals force and the electrostatic force are considered to act independently and to add to each other. The equilibrium DLVO force is the sum of the van der Waals force and of the electrostatic force (2.15),

   (3.21) \begin{equation} \quad F_{\textit{DLVO}} = -\frac {A_{\textit{SLS}}R}{6D^2}+2\pi R \int _D^\infty \pi _{\textit{d,e}}(h) \textrm{d}h, \quad \pi _{\textit{d,e}}(h)=2n_0k_{\textit{B}}T(\cosh \psi _{\textit{m}}-1) ,\end{equation}

   where $A_{\textit{SLS}}$ is the solid–liquid–solid Hamaker constant. The equilibrium midspace plane potential $\psi _{\textit{m}}$ depends only on $h$ , $n_0$ and $\sigma$ . In case of charge regulation, $\sigma =\sigma _{\textit{eq}}(h)$ is a function of the thickness $h$ of the film.

2. (ii) Reynolds lubrication force. In the absence of surface charge the transport coefficients and the ion excess vanish. The dynamic pressure is given by the usual Reynolds lubrication equation:

   (3.22) \begin{equation} \frac {\textrm{d}\varPi _{\textit{Rey}}}{\textrm{d}h} =\frac {R\dot D}{4T_{\textit{vv}}}=\frac {6 \eta R\dot D}{h^3} .\end{equation}
   The Reynolds lubrication equation contributes to the so-called Reynolds force acting during the drainage of a neutral liquid between two spheres of radius $2R$ or a sphere of radius $R$ and a plane:
   (3.23) \begin{equation} F_{\textit{Rey}}=-\frac {6\pi \eta R^2 \dot D}{D} .\end{equation}
   At each time the Reynolds force depends only on the gap $D$ and the velocity $\dot D$ , and is a pure damping force (proportional to $\dot D$ ).

3. (iii) The EK force. Equations (3.19) are linear in terms of $\partial \varPi /\partial h$ , therefore the pressure can be split in three terms, $\varPi =\varPi _{\textit{Rey}}+\varPi _{\textit{EK}}+\varPi _{\textit{DK}}$ defined from

   (3.24a) \begin{align} \frac {\partial \varPi _{\textit{EK}}}{\partial h} =\frac {\partial \varPi _{\textit{Rey}}}{\partial h}\frac {\xi _{\textit{ev}}^2}{(1-\xi _{\textit{ev}}^2)} , \end{align}
   (3.24b) \begin{align} \frac {\partial \varPi _{\textit{DK}}}{\partial h} = -T_{\textit{a}}(h,t)\frac {\partial \mu _{\textit{s}}}{\partial h},\quad T_{\textit{a}}=\frac {T_{\textit{vs}}}{T_{\textit{vv}}}\frac {1-\xi _{\textit{ev}}\xi _{\textit{es}}/\xi _{\textit{vs}}}{1-\xi _{\textit{ev}}^2}, \end{align}

   with zero-value at $h=\infty$ . The EK pressure $\varPi _{\textit{EK}}$ depends exclusively on the EK transport coefficients. It defines an EK contribution $F_{\textit{EK}}$ to the force:

   (3.25) \begin{equation} F_{\textit{EK}}=-12\pi \eta R^2 \dot D \int _D^\infty \textrm {d}h\int _h^{\infty } \frac {\xi _{\textit{ev}}^2}{h^{\prime 3}(1-\xi _{\textit{ev}}^2)} \textrm {d}h' .\end{equation}
   The EK force adds to the DLVO and the Reynolds forces. The figure of merit $\xi _{\textit{ev}}^2(h,n(h,t),\sigma (h,t))$ is fully defined from the relations given in Appendix C and the subsequent appendices. As it is always positive and less than unity, the EK force always increases the energy dissipation with respect to the Reynolds force acting alone.

   Note that if one neglects the changes in the average ion concentration $n_{\textit{s}}$ and the associated transport, then $n(r,t) \equiv n_0$ everywhere, the ion excess $N_{\textit{s}}$ is constant, its time-derivative is zero and the midplane potential $\chi _{\textit{m}} = \psi _{\textit{m}}$ keeps its equilibrium value. Within this approximation, the electroviscous force is reduced to the sole EK force, which in addition is a pure damping as the transport coefficient $\xi _{\textit{ev}}$ in (3.25) keeps its equilibrium value. We discuss in § 5 the physical origin of $F_{\textit{EK}}$ , and the variation of the EK damping with the gap between the surfaces in the case of a constant and uniform surface charge.

4. (iv) The DK force. The remaining contribution to the pressure $\varPi _{\textit{DK}}$ defines the last contribution to the force,

   (3.26) \begin{equation}\quad F_{\textit{DK}}=2\pi R\int _D^\infty \textrm {d}h\left [2nk_{\textit{B}}T(\cosh \chi _{\textit{m}}-1)-\pi _{\textit{d,e}}(h)+\int _h^\infty T_{\textit{a}}(h^{\prime })\frac {\partial \mu _{\textit{s}}}{\partial h^{\prime }} \textrm {d}h^{\prime } \right ] \quad \end{equation}

   where $\chi _{\textit{m}}$ is the function of $(h,n,\sigma )$ given by (3.6) and $T_{\textit{a}}(h,n,\sigma )$ is given by (3.24b ). The chemical potential $\mu _{\textit{s}}$ is solution of the second-order differential equation obtained by substituting (3.19a ) into (3.19b ):

   (3.27a) \begin{align} &\frac {\partial }{\partial h}\left \{(h-D)\left (T_{\textit{b}}\frac {\partial \mu _{\textit{s}}}{\partial h}+\frac {R\dot D}{4}\left [T_{\textit{a}}+2n_0(e^{\mu _{\textit{s}}/k_{\textit{B}}T}-1) \right ] \right )\right \}\nonumber\\ &\quad = \frac {R}{2} \left ( \frac {\partial N_{\textit{s}}}{\partial t} -\frac {\partial N_{\textit{s}}^{\textit{source}}}{\partial t} \right )\!, \end{align}
   (3.27b) \begin{align} T_{\textit{b}}(h,n,\sigma )=T_{\textit{ss}}\frac {1-\xi _{\textit{ev}}^2-\xi _{\textit{es}}^2-\xi _{\textit{vs}}^2+2\xi _{\textit{ev}}\xi _{\textit{es}}\xi _{\textit{vs}}}{1-\xi _{\textit{ev}}^2}= \frac {\textrm {Det} \unicode{x1D64F}}{T_{\textit{vv}}T_{\textit{ee}}(1-\xi _{\textit{ev}}^2)}. \end{align}

   We call $F_{\textit{DK}}$ the DK force as it takes into account all contributions associated with the changes in the local ions concentration: the additional dynamic pressure required to balance the diffusio–osmotic flow and the modified electro–osmotic flow, the changes in osmotic pressure in the midplane, and the eventual release of ions in solution due to dynamic surface charge regulation processes. The latter have to be specified by a chemical model for the charge regulation. The (3.27a ) governing the chemical potential is not only highly nonlinear but also history-dependent, due to the effect of ions accumulation in the right-hand side. Therefore, the DK force depends on the whole drainage history, and can be calculated only for an entire drainage driving function $D(t)$ starting from $t=0$ .

In the remainder of this article, we focus on small-amplitude oscillatory drainage around a nominal intersphere separation. Such conditions are readily accessible with most local probe instruments such as CP-AFM and D-SFA set-ups. In addition, small-amplitude oscillatory forcing is an appropriate driving to study electroviscous force as it limits the impact of their nonlinear and history-dependent character.

We will not consider the case of dynamic charge regulation processes, and we will assume that the surface charge is frozen to its equilibrium value during the harmonic oscillation.

Finally for simplicity, we will treat in detail the case of a uniform surface charge, and we will indicate how to take into account eventual effects of static (not time-dependent transport) charge regulation effects if they appear in the equilibrium DLVO force.

4.1. Electrohydrodynamic impedance

In this section the spheres oscillate around a constant gap value $D$ with a small amplitude $h_0\ll D$ and frequency $\omega /2\pi$ . The quantity of interest is the component of the force at the exciting frequency in the permanent regime, which is proportional to $h_0$ . We use the complex notation for all the quantities proportional to $h_0$ which oscillate at the frequency $\omega /2\pi$ . More precisely we will write

(4.1) \begin{align} h(r,t)=h(r)+h_0\text{Re} [e^{i\omega t}], \quad h(r)=D+\frac {r^2}{2R}, \end{align}

(4.2) \begin{align} F (D,t) =F_{\textit{eq}} (D) - h_0\text{Re} [Ze^{i\omega t}]+O (h_0^2,t), \quad Z(D,\omega ) \in \mathbb C ,\end{align}

where $F_{\textit{eq}} (D)$ refers to the equilibrium situation ( $h_0=0$ ) and $Z(D,\omega )$ is the electrohydrodynamic impedance of the system. The real part of $Z$ represents the conservative contribution to the force, while the imaginary part accounts for the dissipative contribution (the negative sign in the definition of $Z$ ensures positive real and imaginary components for usual stiffness and damping). In the following, $h$ will designate the time-averaged profile of the gap $h(r) = D + r^2/2R$ , and no longer the total profile $h(r,t)$ .

The mechanical impedances associated with the DLVO force and the Reynolds lubrication force are, respectively, a pure stiffness (Israelachvili & Tabor Reference Israelachvili and Tabor1972) and a pure damping:

(4.3a) \begin{align} & Z_{\textit{DLVO}}=-\frac {\partial F_{\textit{DLVO}}}{\partial D}=-\frac {A_{\textit{SLS}}R}{3D^3}+2\pi R \pi _{\textit{d,e}}(D), \quad Z_{\textit{Rey}}=\frac {6i\pi \eta \omega R^2}{D}, \end{align}

(4.3b) \begin{align} & \pi _{\textit{d,e}}(D)=2n_0k_{\textit{B}}T (\cosh \psi _{\textit{m}}(D)-1) . \end{align}

4.2. The EK impedance $Z_{\textit{EK}}$

The linearisation of (3.25) with $\dot D=\text{Re} (h_0i\omega e^{i\omega t})$ shows that the EK impedance is also a pure damping:

(4.4) \begin{equation} Z_{\textit{EK}}=12i\pi \eta \omega R^2\int _D^\infty \textrm {d}h\int _{h}^\infty \frac {\xi ^2_{\textit{ev,eq}}}{h^{\prime 3}(1-\xi ^2_{\textit{ev,eq}})}\textrm {d}h^{\prime }. \end{equation}

Here the EK coefficients of the transport matrix are the equilibrium coefficients, i.e. the figures of merit $\xi _{\textit{ev,eq}}$ are taken at their equilibrium value $\xi _{\textit{ev,eq}}=\xi _{\textit{ev}}(h,n_0,\sigma _{\textit{eq}}(h))$ . It is of interest to explicit the damping by introducing lengths scaled by $2 \ell _{\textit{D}}$ such as $\tilde D=D/2 \ell _{\textit{D}}$ , $\tilde h=h/2 \ell _{\textit{D}}$ , $\tilde h_0=h_0/2 \ell _{\textit{D}}$ , as well as a dimensionless transport matrix $\unicode{x1D669}$ at equilibrium (see also (C14) in Appendix C),

(4.5a) \begin{align} t_{ij}=2\tilde \varGamma _{ij}+2K \tilde L_{ij}, \quad K=\frac {D_{\textit{s}} \pi \eta \ell _{\textit{B}}}{k_{\textit{B}}T}, \quad \tilde \varGamma _{ij}=\int _0^{\tilde h/2}\tilde \varGamma _{i} (\tilde z) \tilde \varGamma _{j} (\tilde z) \textrm {d}\tilde z, \end{align}

(4.5b) \begin{align} \tilde \varGamma _{\textit{v}} (\tilde z) =\tilde z ,\quad \tilde \varGamma _{\textit{s}} (\tilde z) =\int _0^{\tilde z} (\cosh \psi -1) \textrm {d}\tilde z' ,\quad \tilde \varGamma _{\textit{e}} (\tilde z) =\int _0^{\tilde z} ({-}\sinh \psi ) \textrm {d}\tilde z' ,\end{align}

(4.5c) \begin{align} \tilde L_{\textit{ss}}=\tilde L_{\textit{ee}}=\tilde \varGamma _{\textit{s}}\left (\frac {\tilde h}{2}\right ) +\frac {\tilde h}{2}+\delta \tilde \varGamma _{\textit{e}}\left (\frac {\tilde h}{2}\right )\!, \quad \tilde L_{\textit{se,es}} =\tilde \varGamma _{\textit{e}}\left (\frac {\tilde h}{2}\right ) +\delta \left (\tilde \varGamma _{\textit{s}}\left (\frac {\tilde h}{2}\right ) +\frac {\tilde h}{2}\right )\!,\end{align}

where all coefficients are calculated with the equilibrium surface charge and potential. The non-dimensional diffusivity $K$ was introduced in the study of the conductance of EDLs (Bikerman Reference Bikerman1933, Reference Bikerman1935; Squires & Bazant Reference Squires and Bazant2004), and quantifies the ratio of the electromigration velocity and the electro–osmotic velocity for an ion of average mobility. The analytic expressions of the $\varGamma _{\textit{s,eq}}$ and $\varGamma _{\textit{e,eq}}$ functions as well as the relations between the $t_{ij}$ and the $T_{ij}$ are given in Appendix C. The EK impedance $Z_{\textit{EK}}$ can be expressed as

(4.6) \begin{equation} Z_{\textit{EK}}=i\omega \frac {\pi \eta R^2}{2\ell _{\textit{D}}}\mathcal Z_{\textit{EK}}, \quad \mathcal Z_{\textit{EK}}= \int _{\tilde D}^\infty \textrm {d}\tilde h\int _{\tilde h}^\infty \frac {12\xi _{\textit{ev}}^2}{\tilde h^{\prime 3}(1-\xi _{\textit{ev}}^2)}\textrm {d}\tilde h'. \end{equation}

The amplitude of the EK damping can be compared with the Reynolds damping by defining the EK damping factor $\tilde f_{\textit{EK}}$ as

(4.7) \begin{equation} \tilde f_{\textit{EK}}=\frac {Z_{\textit{EK}}}{Z_{\textit{Rey}}}=2\tilde D\int _{\tilde D}^\infty \textrm {d}\tilde h\int _{\tilde h}^{\infty } \frac {\xi _{\textit{ev,eq}}^2}{\tilde h^{\prime 3}(1-\xi _{\textit{ev,eq}}^2)} \textrm {d}\tilde h', \quad \xi _{\textit{ev,eq}}^2=\frac {12\,t_{\textit{ev}}^2}{\tilde h^3 t_{\textit{ee}}}. \end{equation}

The factor $\tilde f_{\textit{EK}}$ is a dimensionless quantity which depends only on dimensionless combinations of variables: the reduced distance $D/2\ell _{\textit{D}}$ , and three dimensionless parameters which are further discussed in the paragraph (4.4) below. The number $K$ and the diffusion contrast $\delta$ gather all the transport coefficients in the electrolyte. They characterise the amplitude of the diffusion-migration part $\tilde {\unicode{x1D647}}$ of the transport matrix with respect to its convective part $\tilde {\boldsymbol{\varGamma }}$ . The electrical properties are fully characterised by the number $\alpha$ defined in (3.2). Note that in the case where the equilibrium DLVO force $F_{\textit{DLVO}}(D)$ exhibits charge regulation, then $\alpha$ should be a function of $h$ determined consistently from the equilibrium DLVO force.

Finally, as both Reynolds and EK impedances are dampings, their ratio does not depend on frequency.

4.3. The DK impedance $Z_{\textit{DK}}$

The DK impedance is obtained by linearising (3.27a ) and (3.26). We do not consider dynamic charge regulation and ignore the ion source term in (3.26). The chemical potential $\mu _{\textit{s}}$ (3.19a ) vanishes at equilibrium and is to the leading order proportional to $h_0$ , thus all other terms in the left-hand side of (3.27a ) reduce to their equilibrium value.

The variation of the number of ions per unit surface $\textrm{d}N_{\textit{s}}$ is expressed by introducing the derivatives

(4.8) \begin{equation} c_{h}(\tilde h)=\frac {1}{n_0}\left .\frac {\partial N_{\textit{s}}(h,n)}{\partial h}\right \vert _{h,n_0} \quad c_{n}(\tilde h)=\frac {1}{\ell _{\textit{D}}}\left .\frac {\partial N_{\textit{s}}(h,n)}{\partial n}\right \vert _{h,n_0} \end{equation}

whose expressions are given in (E12). One also introduces a dimensionless chemical potential response $m_{\textit{s}}$ , a characteristic frequency $\omega _{\textit{c}}$ and a frequency ratio as

(4.9) \begin{equation} \frac {\mu _{\textit{s}}(\tilde h,t)}{2k_{\textit{B}}T \tilde h_0}=\text{Re} [m_{\textit{s}} (\tilde h)e^{i\omega t}] \quad m_{\textit{s}} \in \mathbb C\quad \omega _{\textit{c}}=\frac {D_{\textit{s}}}{2 R \ell _{\textit{D}}} \quad \beta =\frac {K\omega }{4\omega _{\textit{c}}}. \end{equation}

With these dimensionless quantities the linearised equation (3.27a ) becomes

(4.10) \begin{equation} \frac {\partial }{\partial \tilde h}\left \{(\tilde h-\tilde D)\left (t_{\textit{b}}(\tilde h)\frac {\partial m_{\textit{s}}}{\partial \tilde h}+i\beta t_{\textit{a}}(\tilde h) \right )\right \} = i\beta (c_{h}(\tilde h)+c_{n}(\tilde h)m_{\textit{s}}) ,\end{equation}

where the dimensionless $t_{\textit{a}}$ and $t_{\textit{b}}$ are obtained from the equilibrium values of the functions $T_{\textit{a}}$ , $T_{\textit{b}}$ defined in (3.24b ), (3.27b ) by replacing the $T_{ij}$ with the $t_{ij}$ . The angular frequency $\omega _{\textit{c}}$ is built on the average diffusion time of the ions over the lateral extent $\sqrt {2R\ell _{\textit{D}}}$ of the EDLs overlap when the surfaces are at contact. The linear differential equation (4.10) is to be solved with the boundary conditions of vanishing $m_{\textit{s}}$ at infinity and finite second derivative at $\tilde h=\tilde D$ .

The linearised equation (3.27a ) for the force $F_{\textit{DK}}$ is

(4.11) \begin{equation} F_{\textit{DK}}=2\pi R \int _D^\infty \textrm {d}h \left [2n_0\mu _{\textit{s}}\left (\cosh \psi _{\textit{m}}-1+k_{\textit{B}}T\frac {\partial (\cosh \chi _{\textit{m}})}{\partial \mu _{\textit{s}}} \right ) + \int _{h}^\infty T_{\textit{a}}(h^{\prime }) \frac {\textrm {d} \mu _{\textit{s}}}{\textrm {d} h^{\prime }}\textrm {d}h' \right ]\!. \quad \quad \end{equation}

We further show in Appendix E that the quantity $\cosh \psi _{\textit{m}}-1 + k_{\textit{B}}T\partial (\cosh \psi _{\textit{m}}) / \partial \mu _{\textit{s}}$ is equal to the function $c_{h} (\tilde h)$ , so that the DK impedance $Z_{\textit{DK}}$ writes as

(4.12) \begin{equation} Z_{\textit{DK}}=8\pi R n_0k_{\textit{B}}T \mathcal Z_{\textit{DK}}, \quad \mathcal Z_{\textit{DK}}=-\int _{\tilde D}^\infty \textrm {d}\tilde h \left [ m_{\textit{s}}c_{h}(\tilde h)+\int _{\tilde h}^\infty t_{\textit{a}}(\tilde h^{\prime })\frac {\textrm{d}m_{\textit{s}}(\tilde h^{\prime })}{\textrm{d}\tilde h^{\prime }}\textrm {d}\tilde h^{\prime } \right ]\!. \end{equation}

Clearly the DK impedance is not a pure damping and has also a real component, i.e. a stiffness. In order to compare the DK damping to the Reynolds damping, we define the DK damping factor $\tilde f_{\textit{DK}}$ in analogy with (4.7):

(4.13) \begin{equation} \tilde f_{\textit{DK}}=\frac {\text{Im} (Z_{\textit{DK}})}{\text{Im} (Z_{\textit{Rey}})}=\beta ^{-1}\text{Im} (\mathcal Z_{\textit{DK}})\frac {\tilde D}{6}. \end{equation}

4.4. Summary

In summary, our theory describes the mechanical impedance $Z$ of the electroviscous forces in the oscillatory drainage flow of a monovalent electrolyte between dielectric surfaces under uniform and constant (not time-dependent) surface charge and no-slip boundary conditions as

(4.14a) \begin{align}\text{Re} (Z) =8\pi R n_0k_{\textit{B}}T\left (\frac {1}{2}[\cosh \psi _{\textit{m}} (\tilde D,\alpha ) -1]+\text{Re} (\mathcal Z_{\textit{DK}}) \right )\!, \end{align}

(4.14b) \begin{align} \text{Im} (Z) = \frac {6 \pi \eta \omega R^2}{D} (1+\tilde f_{\textit{EK}}+\tilde f_{\textit{DK}}), \end{align}

where $\mathcal Z_{\textit{DK}}$ , $\tilde f_{\textit{EK}}$ and $\tilde f_{\textit{DK}}$ are given by (4.12), (4.7) and (4.13). These non-dimensional quantities depend on the normalised distance $\tilde D=D/2\ell _{\textit{D}}$ , and on four non-dimensional parameters as follows.

1. (i) The parameter $\alpha =2\ell _{\textit{D}}/\ell _{\textit{G}}$ characterises the nominal surface charge $\sigma _0$ of non-interacting surfaces in the solution of concentration $n_0$ . More specifically $\alpha$ is the dimensionless electric field at the non-interacting surfaces, denoted as $\mathcal E$ in Adar & Andelman (Reference Adar and Andelman2018). It is related by the Grahame (Reference Grahame1953) equation to the dimensionless nominal zeta potential: $\displaystyle \alpha = 2 \sinh (\tilde \zeta /2)$ . At 300 K, the value $\alpha =1$ corresponds to a nominal surface potential of $24\,\textrm {mV}$ , $\alpha =10$ to $\zeta =115\,\textrm {mV}$ , $\alpha =100$ to $\zeta =230\,\textrm {mV}$ .

   If the interacting surfaces exhibit static charge regulation, with a surface charge that does not depend on time, then $\alpha$ should be replaced by a local thickness-dependent parameter $\alpha _{\textit{l}}=\alpha \sigma _{\textit{eq}}(h)/\sigma _0$ determined from the equilibrium DLVO force.

2. (ii) The parameter $K=D_{\textit{s}}\pi \eta \ell _{\textit{B}}/ (k_{\textit{B}}T)$ characterises the average ionic mobility; $K^{-1}$ is also referred to as an intrinsic Péclet number (Yariv et al. Reference Yariv, Schnitzer and Frankel2011). If anions and cations have a similar diffusion coefficient, $K$ is related to the hydrodynamic radius of the ions defined by the Stokes–Einstein relation $D_{\textit{s}}=k_{\textit{B}}T/6\pi \eta r_{\textit{H}}$ by $K= \ell _{\textit{B}}/6 r_{\textit{H}}$ . For a KCl electrolyte of average diffusion coefficient 2.0 $\times 10^{-9}\textrm {m}^2\,\textrm{s}^{-1}$ in water at 300 K (Haynes Reference Haynes2016), $K=1.0$ .

3. (iii) The parameter $\delta = (D^+-D^-) / (D^++D^-)$ quantifies the contrast between the diffusion coefficients of coions and counterions.

4. (iv) In addition the DK impedance depends on the frequency ratio $\beta =K\omega /4\omega _{\textit{c}}$ related to the average ionic diffusion time over the lateral extent $\sqrt {2R\ell _{\textit{D}}}$ of the EDLs overlap at spheres contact. For D-SFA with $R=3\,\textrm {mm}$ and a KCl electrolyte of mM concentration, the system frequency $f_{\textit{c}}=\omega _{\textit{c}}/2\pi$ is $f_{\textit{c}} \simeq 5\,\textrm {Hz}$ . For CP-AFM with $R=10$ ${\unicode{x03BC}}$ m in the same conditions, $f_{\textit{c}}\simeq 1.6\,\textrm {kHz}$ .

In § 5 and 6 we study the EK and DK contributions to the electroviscous impedance over a very extensive range of reduced distance $D/2\ell _{\textit{D}}$ in order to highlight the main tendencies. It should be emphasised that in a given experiment the portion of the distance range which has a physical meaning is reduced, and depends on the experimental configuration. We discuss in §§ 7.1 the limitation of our approach in terms of real distances.

Furthermore, our result as formulated in (4.14a ) express the electroviscous impedances in terms of relative increase with respect to non-electroviscous forces, i.e. the DLVO stiffness and the Reynolds damping. In order to have a rough idea of the relative impact of the electroviscous effects in a given experimental configuration, the experimentalist should calculate the dimensionless numbers corresponding to the experimental parameters and then refer to figures below whose caption starts by ‘Relative EK damping’, ‘High frequency DK stiffness’ or ‘Low frequency DK damping’. Note that the probe radius enters only in frequency considerations, and for determining the higher distance limit compatible with the LTE approximation. As an example, an experiment performed with NaCl at 1 mM in water on glass surfaces (with surface charge 8 mC m $^{-2}$ ) corresponds to $\alpha =4.5$ and $K=0.8$ , the DK damping is negligible in these conditions and the maximum effect of EK damping should not exceed 10 % of the Reynolds damping. For more accurate calculation Python routines are given in the Supplementary material available at https://doi.org/10.1017/jfm.2026.12093.

5.1. Non-overlapping EDLs

When the EDLs are sufficiently far from each other the transport coefficients split as the sum of a bulk contribution across the half-thickness of the electrolyte film, and an excess or surface contribution, due to transport in the EDLs. The transport coefficients are affine functions of $h$ of the form $t_{ij}=a_{ij}\tilde h+b_{ij}$ , with the $a_{ij}$ and $b_{ij}$ given in Appendix F. This holds for all transport coefficients except for the permeability which is always equal to $t_{\textit{vv}}=\tilde h^3/12$ (i.e. $T_{\textit{vv}}=h^3/24 \eta$ ) whatever the thickness $h$ of the film.

• (i) Regime I: non-overlap-low-Dukhin.

At very large distance, the bulk contribution to the coefficients $T_{\textit{ev}}$ and $T_{\textit{ee}}$ dominates over the EDLs one. The conductance of the (half-)film is $T_{\textit{ee}}=(2n_0 e^2D_{\textit{s}}/k_{\textit{B}}T)h/2$ , the electro–osmotic flow coefficient is $T_{\textit{ev}}=-(\varepsilon _{\textit{l}} \zeta /\eta )h/2$ , giving $\xi _{\textit{ev}}^2=3\tilde \zeta ^2/(4K\tilde h^2)$ . The leading term of $\mathcal Z_{\textit{EK}}$ decays in $D^{-3}$ , as first calculated by Bike & Prieve (Reference Bike and Prieve1990)

(5.1) \begin{equation} \frac {\xi _{\textit{ve}}^2}{1-\xi _{\textit{ve}}^2} \simeq \frac {3\tilde \zeta ^2}{4K\tilde h^2}, \quad \mathcal Z_{\textit{EK,I}}=\frac {3 \tilde \zeta ^2}{4K \tilde D^3}, \quad \frac {Z_{\textit{EK,I}}}{i\omega }=\frac {6\pi R^2k_{\textit{B}}T}{D^3n_0 D_{\textit{s}}}\left ( \frac {\varepsilon _{\textit{l}}\zeta }{e}\right )^2. \end{equation}

In region I the EK damping $Z_{\textit{EK,I}}$ increases monotonically with the surface charge. It does not depend on the fluid viscosity $\eta$ , and it is inversely proportional to the bulk conductivity of the electrolyte $\gamma _{\textit{bulk}}=2n_0 e^2D_{\textit{s}}/k_{\textit{B}}T$ .

Regime I requires both coefficients $T_{\textit{ev}}$ and $T_{\textit{ee}}$ to be dominated by the bulk. This is always the case for the electro–osmotic coefficient in the condition of non-overlapping EDLs ( $\vert b_{\textit{ev}}/a_{\textit{ev}}\tilde D \vert \le \vert b_{\textit{ev}}/a_{\textit{ev}} \vert \le 1$ ). This property stems from the fact that the electro–osmotic velocity in the bulk has always a larger amplitude than in the EDLs, due to the no-slip boundary condition. As for the conductance $T_{\textit{ee}}$ , the relative importance of bulk and EDLs contributions has been discussed by several authors (Bikerman Reference Bikerman1933; Lyklema Reference Lyklema1991; Squires & Bazant Reference Squires and Bazant2004), and is generally expressed in terms of the Dukhin number, defined as the ratio of these two contributions (Squires & Bazant Reference Squires and Bazant2004),

(5.2) \begin{equation} \textit{Du}(\tilde h)=\frac {b_{\textit{ee}}}{a_{\textit{ee}} \tilde h}=\frac {4}{\tilde h} \left (1+\frac {1}{2K}-\frac {\delta }{\gamma }\right )\sinh ^2 \left (\frac {e \zeta }{4k_{\textit{B}}T}\right ) \end{equation}

with $\gamma =\tanh (\tilde \zeta /4)$ . Regime I is thus a low Dukhin number regime. Its frontier $\displaystyle \tilde D_{\textit{I}-\textit{II}}=4 \sinh ^2 ({\tilde \zeta }/{4} ) (1+ ({1}/({2K}))-({\delta }/{\gamma }) )$ is plotted in the phase diagram figure 2(a) in the $(\tilde D,\alpha )$ plane, for the case of ions of identical diffusivity $\delta =0$ and various values of $K$ . The region I extends at larger $\tilde D$ and lower $\alpha$ from this frontier.

• (ii) Regime II: non-overlap-high-Dukhin.

The non-overlap-high-Dukhin region I develops at high $\alpha$ and intermediate distance range $1\lt \tilde D \lt \tilde D_{\textit{I-II}}$ . In this region the figure of merit $\xi _{\textit{ev}}^2 \simeq 3\tilde \zeta ^2/(4b_{\textit{ee}}\tilde h)$ scales as $1/\tilde h$ and the EK damping is

(5.3a) \begin{align} \frac {\xi _{\textit{ve}}^2}{1-\xi _{\textit{ve}}^2} \simeq \frac {3\tilde \zeta ^2}{4b_{\textit{ee}}\tilde h}, \quad \mathcal Z_{\textit{EK,II}}=\frac {3 \tilde \zeta ^2}{2b_{\textit{ee}} \tilde D^2} ,\quad \frac {Z_{\textit{EK,II}}}{i\omega }=\frac {12\pi R^2( \varepsilon _{\textit{l}}\zeta )^2}{D^2\ T_{\textit{ee,surf}}}, \end{align}

(5.3b) \begin{align} \quad T_{\textit{ee,surf}}=\gamma _{\textit{bulk}} \frac {b_{\textit{ee}}\ell _{\textit{D}}}{K}=\gamma _{\textit{bulk}}4\ell _{\textit{D}}\left ( 1-\frac {\delta }{\gamma }+\frac {1}{2K}\right )\sinh ^2 \frac {\tilde \zeta }{4} ,\quad \gamma _{\textit{bulk}}=\frac {2n_0 e^2D_{\textit{s}}}{k_{\textit{B}}T}. \end{align}

Regime II is observed on figure 2(b) where $\mathcal Z_{\textit{EK}}(\tilde D)$ is plotted at the fixed value $\alpha =100$ for various value of $K$ . It can be seen that the EK damping in regime II is not always proportional to the inverse of the electrolyte conductivity: this is due to the role of the electro–osmotic convection of ions which dominates the surface conductance at low diffusivity (low $K$ ).

Experimentally, the measurement of the EK damping $Z_{\textit{EK,II}}$ in regime II gives access to the surface conductance if the $\zeta$ -potential of the surfaces is known (or measured in regime I). Thus, (5.3a ) opens a way to test the Bikerman relation (5.3b ) without adjustable parameter, as the equilibrium value of the surface charge could be measured on the very same surfaces through the DLVO force.

Figure 3.

(a) Normalised hydrodynamic distance $\tilde D_{\textit{H}}= \tilde D/ (1+\mathcal Z_{\textit{EK}}/\mathcal Z_{\textit{Rey}})$ as a function of $\tilde D$ , for $K=0.01$ , $\delta = 0$ , and various values of $\alpha$ . The dashed line is the affine function $\tilde D_{\textit{H}}=\tilde D-2\tilde b$ , expected in the case of a stagnant fluid layer of thickness $b$ on each surface. (b) The dimensionless thickness of the apparent stagnant layer at the walls $\tilde b=\tilde \zeta ^2/8b_{\textit{ee}}$ , characterising the effect of the EK damping in the regime II, plotted in colour scale as a function of $K$ on the y-axis and $\alpha$ on the x-axis.

Figure 3. Long description

Panel A: A line graph shows the normalized hydrodynamic distance as a function of D with K equal to 0.01. The x-axis represents D, and the y-axis represents the normalized hydrodynamic distance. Multiple lines are plotted for different values of alpha (100, 30, 3, 1, 0.3), showing how the normalized hydrodynamic distance changes with D. The dashed line represents the affine function D/2, expected in the case of a stagnant fluid layer of thickness B on each surface. Panel B: A heatmap illustrates the dimensionless thickness of the apparent stagnant layer at the walls, characterizing the effect of the EK damping in regime II. The x-axis represents alpha, and the y-axis represents K. The color scale on the right indicates the dimensionless thickness values, ranging from 0 to 0.9.

More specifically, it can be of interest to characterise the excess damping due to the EK effect by the apparent hydrodynamic distance $D_{\textit{H}}$ between the surfaces. The apparent hydrodynamic distance is defined as the virtual position of the surfaces that would provide the same value of the overall damping under a pure lubrication process: $-6\pi \eta R^2 /D_{\textit{H}}=-6\pi \eta R^2 /D-Z_{\textit{EK}}/i\omega$ . Figure 3(a) plots $\tilde D_{\textit{H}}$ as a function of $\tilde D\in [0,8]$ for various values of $\alpha$ , at the fixed value $K=0.01$ chosen so as to emphasise the effect. One observes that when the regime II is well established (here $\alpha \ge 3$ ) $\tilde D_{\textit{H}}$ is essentially an affine function of $\tilde D$ of the form $\tilde D_{\textit{H}} \simeq \tilde D-2 \tilde b$ . This behaviour stems from the $D^{-2}$ scaling of $Z_{\textit{EK,II}}$ , leading to $D_{\textit{H}}=(1/D+Z_{\textit{EK,II}}/(6i\pi \eta R^2 \omega ) )^{-1}\simeq D-2(\varepsilon _{\textit{l}}\zeta )^2/\eta T_{\textit{ee,surf}}$ . Thus, regime II is conveniently described by apparent stagnant layers of thickness $b=(\varepsilon _{\textit{l}}\zeta )^2/\eta T_{\textit{ee,surf}}$ on each solid surface, and the experimental measurement of the thickness of these apparent stagnant layers, or `negative slip effect’, could provide a direct determination of the surface conductance.

Figure 3(b) plots in colour scale the non-dimensional thickness $\tilde b=\tilde \zeta ^2 /8b_{\textit{ee}}$ of the pseudostagnant layer as a function of $\alpha$ and $K$ . It culminates at a value close to unity, corresponding to a maximum thickness $2\ell _{\textit{D}}$ of apparently non-flowing fluid on each surface. Counterintuitively, $\tilde b$ decays monotonically with $\alpha$ and $K$ , due to the fact that the surface conductance increases more rapidly with the surface charge than does the $\zeta$ -potential.

Figure 3(a) also shows the striking non-monotonic behaviour of the EK damping with the surface charge. For a given value of $\tilde D$ , the hydrodynamic distance $\tilde D_{\textit{H}}$ decreases with increasing $\alpha$ when $\alpha \lt 3$ (i.e. regime II is not yet reached) and then increases for larger values of $\alpha$ . This charge inversion in the behaviour of the EK effect was also mentioned in Zhao et al. (Reference Zhao, Zhang, van den Ende and Mugele2020) and takes its origin in the electrical conductance of the film. At low Dukhin number (regime I) the film conductance is set by the electrolyte bulk conductivity; the EK damping steadily grows as its motor mechanisms which are the streaming current and the electro–osmotic flow. But at high Dukhin number the surface conductance takes over and grows more rapidly with the surface charge than does the $\zeta$ -potential, restoring very efficiently the charge equilibrium and leading to the decrease of $\mathcal Z_{\textit{EK}}$ with increasing surface charge.

5.2. Overlapping EDLs

The EDLs start to overlap significantly when $\tilde D \lesssim 1$ . Following Adar & Andelman (Reference Adar and Andelman2018) the equilibrium ions distribution in the confined EDLs falls into three domains (see also Appendix G).

1. (i) The Gouy–Chapman domain.

When $\alpha \gtrsim \pi /2$ and $\pi ^2/2\alpha \lesssim \tilde D \lesssim \pi$ , the midpotential at equilibrium is given by the Gouy–Chapman approximation $\psi _{\textit{m}}=2\ln (\pi /\tilde D)$ . The electrostatic disjoining pressure $\pi _{\textit{d,e}}(D)$ scales as $D^{-2}$ (Adar & Andelman Reference Adar and Andelman2018): $\pi _{\textit{d,e}}(D) \simeq \pi k_{\textit{B}}T/(2\ell _{\textit{B}} D^2)$ . In this case $\mathcal Z_{\textit{EK}}$ obeys a scaling law very close to regime II, except that the term $\tilde \zeta ^2$ of the prefactor is modified by $(\tilde \zeta \ln \tilde D)^2$ (see Appendix G). The variation of $\mathcal Z_{\textit{EK}}$ in the Gouy–Chapman regime is hardly distinguished from the non-overlap-high-Dukhin regime II in figure 2(b).

1. (ii) Overlap-Limited-Debye–Hückel domain (regime IV).

When $2\alpha \lesssim \tilde D \lesssim 1$ the equilibrium potential in the confined EDLs is always less than unity and the ion distributions are described by the Debye–Hückel linearisation of the Poisson–Boltzmann equation. In this domain the disjoining pressure $\pi _{\textit{d,e}}(D)$ also scales as $D^{-2}$ : $\pi _{\textit{d,e}}(D) \simeq 2 k_{\textit{B}}T \ell _{\textit{D}}^2/(\pi \ell _{\textit{B}} \ell _{\textit{G}}^2 D^2)$ . The EK impedance varies as

(5.4) \begin{equation} \frac {\xi _{\textit{ve}}^2}{1-\xi _{\textit{ve}}^2}\simeq \frac {\alpha ^2}{12 K} ,\quad \mathcal Z_{\textit{EK,IV}}=\frac {\alpha ^2}{2K\tilde D}, \quad \frac {Z_{\textit{EK,IV}}}{i\omega }=\frac {\pi R^2 k_{\textit{B}}T\sigma ^2}{D_{\textit{s}}n_0 e^2D}. \end{equation}

This scaling law is plotted as thin line in the region IV of figure 2(c) for the values $\alpha = K =0.01$ . We find that it describes well the region IV in the phase diagram (figure 2 a). It should be mentioned that $\mathcal Z_{\textit{EK}}$ is very low as it is proportional to $\alpha ^2$ , and represents only a very low fraction of the Reynolds damping, as discussed in the next paragraph.

1. (iii) Ideal-gas domain (regime III).

When $\alpha \gtrsim \pi /2$ and $\tilde D \lesssim \pi ^2/2\alpha$ or when $\alpha \lesssim \pi /2$ and $\tilde D \lesssim 2\alpha$ , the confined electrolyte is well described by the ideal-gas model (Adar & Andelman Reference Adar and Andelman2018) or Donnan equilibrium (Donnan Reference Donnan1924), which assumes a uniform ion concentration across the film thickness. In this domain the disjoining pressure $\pi _{\textit{d,e}}(D)$ scales as $D^{-1}$ : $\pi _{\textit{d,e}}(D) \simeq k_{\textit{B}}T /(\pi \ell _{\textit{B}} \ell _{\textit{G}} D)$ . The electro–osmotic velocity profile is parabolic, the average electro–osmotic velocity grows as $h^2$ (see Appendix G), the figure of merit $\xi _{\textit{ev}}^2$ is proportional to $\tilde h$ , and the EK damping varies to the leading order as $-\ln {\tilde D}$ ,

(5.5) \begin{equation} \quad \ \frac {\xi _{\textit{ve}}^2}{1-\xi _{\textit{ve}}^2}\simeq \frac {\alpha \tilde h}{12 K (1-\delta )}, \quad \mathcal Z_{\textit{EK,III}}=-\frac {\alpha \ln \tilde D}{K (1-\delta )}, \quad \frac {Z_{\textit{EK,III}}}{i\omega }=\frac {2\pi R^2 k_{\textit{B}}T\sigma }{D^-e} \ln \frac {2\ell _{\textit{D}}}{D} ,\end{equation}

where $D^-$ is the diffusion coefficient of the counterions. The law (5.5) is plotted in figure 2(b) for $(\alpha ,K) = (100,0.01)$ and in figure 2(c) for $(\alpha ,K) = (0.01,0.01)$ . We find that it adequately describes the variations of $\mathcal Z_{\textit{EK}}$ at very small distance for all values of $(\alpha , K, \delta )$ . It should be emphasised that the law (5.5) involves the mobility $K (1-\delta )$ of the counterions only, not the one of the coions.

5.3. Comparison with the Reynolds damping

We use the relative amplitude factor $\tilde f_{\textit{EK}}$ defined in (4.7), which is also expressed as $\tilde f_{\textit{EK}}= (\tilde D/6) \mathcal Z_{\textit{EK}}$ , to compare the EK damping with the Reynolds damping. Three main points emerge from this comparison illustrated in figure 4.

(i) At intermediate distances the relative amplitude of the EK damping varies in a non-monotonic way with the dimensionless surface charge. This is due to the surface conductance and was discussed in § 5.1.

Figure 4.

(a) Relative EK damping $\tilde f_{\textit{EK}}=\mathcal Z_{\textit{EK}}\tilde D/6$ as a function of the normalised distance $\tilde D=D/2\ell _{\textit{D}}$ for $\alpha =100$ , $\delta =0$ and various values of $K$ . The regions I, II and III of the plot correspond, respectively, to the large, intermediate and short distance variations of $\tilde f_{\textit{EK}}$ for the value $K=0.01$ . (b) Plot of $\tilde f_{\textit{EK}}$ as a function of $\tilde D$ for $K=0.1$ , $\delta =0$ and various values of $\alpha$ . (c) Maximum value of $\tilde f_{\textit{EK}}$ obtained for $\alpha \ge 10$ , plotted as a function of the normalised diffusion coefficient of the counterions $K (1-\delta )$ for $\delta \in ({-}0.7,-0.3,0,0.3,0.7)$ . Values of $K$ are reported in the legend, for each $K$ the values of $\delta$ can be read on the x-axis. The dashed line in $[K (1-\delta ) ]^{-1}$ is a visual guide.

(ii) The EK damping factor increases in general as $K^{-1}$ . However, at large values of $\alpha$ this growth saturates, as the surface conductance becomes dominated by the EK flow and no longer depends on the ion electromigration contribution.

(iii) The damping ratio $\tilde f_{\textit{EK}}$ always reaches a maximum value, which is obtained when EDLs overlap. The maximum occurs at the entrance in the ideal-gas region upon decreasing the distance, and is due to the fact that the electro–osmotic flow becomes Poiseuille-like when the electric charge is uniform in the thickness, which represents a loss in flow rate. This causes the EK damping to quasisaturate (in $\ln \tilde D$ ) in the ideal-gas regime, whereas the Reynolds damping diverges as $\tilde D^{-1}$ . When $\alpha \lt 1$ the value of the maximum $\textrm {max} (\tilde f_{\textit{EK}})$ is proportional to $\alpha ^2$ , so that $\mathcal Z_{\textit{EK}}$ represents only a very small fraction of the overall damping. However, for $\alpha \ge 5$ a remarkable property is that the maximum value $\textrm {max} (\tilde f_{\textit{EK}})$ is independent of $\alpha$ , and depends only on the mobility of the counterions $K(1-\delta )$ (see figure 4 c). The surface gap at which the maximum of $\tilde f_{\textit{EK}}$ is reached, however, decays as $\alpha ^{-1}$ , and might fall in a distance range lacking physical significance (as discussed in § 7).

Overall the EK damping may reach 100 % of the Reynolds damping for $K (1-\delta ) \le 0.3$ , which corresponds to counterions of hydrodynamic radius equal or larger than $0.4\,\textrm {nm}$ in water at 300 K. Furthermore, if the cations and anions have a significant mobility contrast, measuring the maximum value $\textrm {max}(\tilde f_{\textit{EK}})$ may allow one to determine the sign of the surface charge, a property that cannot be distinguished from the DLVO force measurement alone.

6.1. Variations with frequency: the spring-and-dashpot damper

At high frequency, the ions follow adiabatically the position of the sphere independently of its velocity. The modification of the disjoining pressure, as well as the diffusio–osmotic flow induced by the spatial gradient of the modified local concentration, are in phase with the surface position. The additional pressure $\varPi _{\textit{DK}}$ needed to counterbalance the diffusio–osmotic flow is thus also in phase with the displacement. In this high frequency limit, one expects the stiffness $\text{Re} (Z_{\textit{DK}})$ to dominate the DK impedance. At lower frequency, this picture is smeared out by the diffusion of ions.

These characteristics are typical of a damper. Indeed, the variations of the real and imaginary part of $\mathcal Z_{\textit{DK}}$ as a function of $\beta =K\omega /4\omega _{\textit{c}}$ , shown in figure 5(a) for a particular value of ( $\tilde D, \alpha , K)$ , correspond very well to a damper made of spring-and-dashpot in series (with no inertia), of the form $Z= (k+ik^2/\lambda \omega ) / (1+k^2/\lambda ^2\omega ^2)$ (denoting $k$ the spring stiffness and $\lambda$ the friction coefficient). At high frequency the impedance is dominated by a real part independent of frequency, whereas at low frequency the leading contribution is an imaginary part proportional to $\omega$ . It is important to remind that the quantity $\beta ^{-1} \text{Im} (\mathcal Z_{\textit{DK}}) =6\tilde f_{\textit{DK}}/\tilde D$ determines the relative DK damping coefficient $\tilde f_{\textit{DK}}$ in (4.13). Thus, $\tilde f_{\textit{DK}}$ becomes frequency independent at low frequency, which facilitates the comparison of the DK damping to the Reynolds damping in § 6.3.

Figure 5.

(a) Real and imaginary parts of $\mathcal Z_{\textit{DK}}$ as a function of $\beta =K\omega /4\omega _{\textit{c}}$ for $\tilde D=1$ , $\alpha =10$ and $K=1$ . (b) The ratio $\omega _0/\omega _{\textit{c}}=4\beta _0/K$ for which $\text{Re} (\mathcal Z_{\textit{DK}}) =\text{Im} (\mathcal Z_{\textit{DK}})$ at the reduced distance $\tilde D=1$ , plotted as a function of $K$ for various values of $\alpha$ . (c) Real part of $\mathcal Z_{\textit{DK}}$ plotted as a function of the reduced distance $\tilde D$ for $K=1$ , $\alpha =10$ and various values of $\omega /\omega _{\textit{c}}$ . (d) Imaginary part of $\beta ^{-1}\mathcal Z_{\textit{DK}}$ plotted as a function of $\tilde D$ for $K=1$ , $\alpha =10$ and various values of $\omega /\omega _{\textit{c}}$ .

Figure 5. Long description

Panel A: A line graph shows the real and imaginary parts of the DK impedance as a function of beta for specific values of alpha, K, and D-tilde. The x-axis represents beta on a logarithmic scale, and the y-axis represents the real and imaginary parts of the DK impedance on a logarithmic scale. The graph includes two data series: one for the real part and one for the imaginary part, with different markers and colors. Panel B: A line graph plots the ratio of frequencies for which a certain condition holds at a reduced distance as a function of K for various values of alpha. The x-axis represents K on a logarithmic scale, and the y-axis represents the ratio on a logarithmic scale. Different data series are shown for different values of alpha. Panel C: A line graph shows the real part of the DK impedance as a function of the reduced distance for specific values of K, alpha, and various values of omega over omega-c. The x-axis represents the reduced distance on a logarithmic scale, and the y-axis represents the real part of the DK impedance on a logarithmic scale. Different data series are shown for different values of omega over omega-c. Panel D: A line graph plots the imaginary part of the DK impedance as a function of the reduced distance for specific values of alpha, K, and various values of omega over omega-c. The x-axis represents the reduced distance on a logarithmic scale, and the y-axis represents the imaginary part of the DK impedance on a logarithmic scale. Different data series are shown for different values of omega over omega-c.

The asymptotic limits of the dominant components of $\mathcal Z_{\textit{DK}}$ are obtained by solving the differential equation (4.10) in the limits $\beta =\infty$ (superscript hf) and $\beta =0$ (superscript lf). One obtains

(6.1a) \begin{align} \text{Re} [m_{\textit{s}}^{\textit{hf}}(\tilde h)] =\frac {t_{\textit{a}} (\tilde h) -c_{h} (\tilde h) +t_{\textit{a}}'(\tilde h) (\tilde h-\tilde D)}{c_n(\tilde h)},\end{align}

(6.1b) \begin{align} \frac {d \text{Im} (\beta ^{-1}m_{\textit{s}}^{\mathrm{\ell} {f}})}{\textrm{d}\tilde h} = \text{Im} (\beta ^{-1}m_{\textit{s}}^{\mathrm{\ell} {f}}) ' =-\frac {t_{\textit{a}} (\tilde h)}{t_{\textit{b}} (\tilde h)}+\frac {G (\tilde D) -G (\tilde h)}{t_{\textit{b}} (\tilde h) (\tilde h - \tilde D) },\nonumber\\ G (\tilde h) =\int _{\tilde h}^\infty c_{h} (\tilde h^{\prime }) \textrm {d}\tilde h'. \end{align}

The asymptotic limits of $\mathcal Z_{\textit{DK}}$ are then obtained from (4.12) and computed by a numerical integration (details of implementation in the Supplementary material):

(6.2a) \begin{align} \text{Re} (\mathcal Z_{\textit{DK}}^{hf}) =\int _{\tilde D}^\infty \left ( (t_{\textit{a}}-c_{h}) \text{Re} (m_{\textit{s}}^{hf}) + \int _{\tilde h}^\infty t_{\textit{a}}' (\tilde h^{\prime })\text{Re} (m_{\textit{s}}^{hf}) \textrm {d}\tilde h'\right ) \textrm {d}\tilde h, \end{align}

(6.2b) \begin{align} \beta ^{-1}\text{Im} (\mathcal Z_{\textit{DK}}^{\ell \textit{f}}) =\int _{\tilde D}^\infty \left (c_{h} (\tilde h) \int _{\tilde h}^\infty \text{Im} (\beta ^{-1}m_{\textit{s}}^{\mathrm{\ell} {f}}) ' \textrm {d}\tilde h'- \int _{\tilde h}^\infty t_{\textit{a}} (\tilde h^{\prime }) \text{Im} (\beta ^{-1}m_{\textit{s}}^{\mathrm{\ell} {f}}) ' \textrm {d}\tilde h'\right ) \textrm {d}\tilde h. \end{align}

Figures 5(c) and 5(d) show the convergence of $\text{Re} (\mathcal Z_{\textit{DK}})$ and $\beta ^{-1}\text{Im} (\mathcal Z_{\textit{DK}})$ to the above asymptotic limits for the particular value $(\alpha ,K) = (10,0.1)$ . From this frequency study, one can draw the following important conclusions.

1. (i) The asymptotic limits are always upper bounds of the corresponding quantities: $\text{Re} (\mathcal Z_{\textit{DK}}) \le \text{Re} (\mathcal Z_{\textit{DK}}^{hf})$ and $\tilde f_{\textit{DK}} \le \tilde f_{\textit{DK}}^{\ell \textit{f}}$ . This property is expected for a spring-and-dashpot in series, and is verified for all distances and parameters $\alpha$ and $K$ . Therefore, in § 6.2 we use $\text{Re} (\mathcal Z_{\textit{DK}}^{hf})$ to compare the DK stiffness to the DLVO stiffness, keeping in mind that it is the highest possible DK effect. Indeed, at low frequency $\text{Re} (\mathcal Z_{\textit{DK}})$ decays as $\omega ^2$ and becomes negligible. Symmetrically we use the low-frequency limit $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ to compare the DK damping to the Reynolds damping in § 6.3, remembering that $\tilde f_{\textit{DK}}$ decays in $1/\omega ^2$ and becomes negligible at high frequency.

2. (ii) The crossover frequency $\omega _0$ separating the low and high-frequency domain of $\mathcal Z_{\textit{DK}}$ is essentially of the order of $\omega _{\textit{c}}$ when $\tilde D \le 10$ . The value of $\omega _0/\omega _{\textit{c}}=4\beta _0/K$ is determined from the value of $\beta _0$ at which $\text{Re} (\mathcal Z_{\textit{DK}})$ and $\text{Im} (\mathcal Z_{\textit{DK}})$ become equal (see figure 5 a). Figure 5(b) indeed shows that the ratio $\omega _0/\omega _{\textit{c}}$ obtained at $\tilde D=1$ is essentially close to unity for a wide range of values of $K$ and $\alpha$ . Figure 5(c) shows that in the distance range $\tilde D \le 10$ the high-frequency stiffness is essentially reached when $\omega /\omega _{\textit{c}}\ge 5$ while the asymptotic low-frequency damping is essentially reached for $\omega /\omega _{\textit{c}}\le 0.2$ (figure 5 d).

3. (iii) At larger distances, the asymptotic limits are reached for lower frequency. More specifically, the convergence at large distance $\tilde D$ is essentially obtained for $\omega \simeq D_{\textit{s}}/(4R \ell _{\textit{D}} \tilde D)$ , which is the inverse of the diffusion time over the lateral extension of the hydrodynamic contact $\sqrt {2RD}$ .

Thus, at moderate distance in terms of $\ell _{\textit{D}}$ the system time governing the DK response is the ionic diffusion time over the lateral extension $\sqrt {2R\ell _{\textit{D}}}$ of the overlapping EDLs. From a practical perspective, the relevant frequency domain for a distance range $D\lesssim 100$ nm and a KCl electrolyte with concentration of 1 mM ( $\ell _{\textit{D}} \sim 10$ nm) is as follows.

The CP-AFM with a probe radius $R=1$ μm: $3\,\textrm {kHz} \leq f \leq 75\,\textrm {kHz}$ is the intermediate frequency range, fully developed DK damping is reached at $ f\leq 3\,\textrm {kHz}$ , and fully developed DK stiffness would require $f \geq 75\,\textrm {kHz}$ . Note that this latter frequency is close to the upper limit frequency of 100 kHz for the LTE to remain valid at the distance range considered, thus the fully developed DK stiffness is hardly reached with a 1 μm probe radius. For CP-AFM with a probe radius $R=100$ μm, the above frequencies for the DK electroviscous effect are divided by the ratio 100 and fully fall in the range of the LTE condition.

The D-SFA with a probe radius $R=3$ mm: the intermediate frequency range is $1\,\textrm {Hz} \leq f \leq 25\,\textrm {Hz}$ , the DK damping is fully established at $f \leq 1\,\textrm {Hz}$ and the DK stiffness fully established at $f \geq 25\,\textrm {Hz}$ . The above frequencies scale as $\sqrt {n_0}$ while the corresponding distance range scales as $1/\sqrt {n_0}$ .

6.2. High-frequency DK stiffness

Figure 6.

(a) High frequency DK stiffness $\text{Re} (\mathcal Z_{\textit{DK}}^{hf})$ as a function of $\tilde D$ for $\alpha =100$ and various values of $K$ . (b) Solid lines: high-frequency DK stiffness $\text{Re} (\mathcal Z_{\textit{DK}}^{hf})$ for various values of $\alpha$ and $K=0.1$ ; dashed lines: non-dimensional DLVO stiffness $\mathcal Z_{\textit{DLVO}}=Z_{\textit{DLVO}}/8\pi R n_0k_{\textit{B}}T$ for $\alpha =1$ , 10 and 100. For higher values of $\alpha$ the DLVO stiffness is essentially the same as for $\alpha =100$ . For clarity the van der Waals contribution to the DLVO stiffness has not been considered.

Figure 6 represents the high-frequency stiffness $\text{Re} (\mathcal Z_{\textit{DK}}^{hf})$ as a function of $\tilde D$ for various values of $K$ .

When EDLs overlap, the DK stiffness arises from the accumulation of ions in the midplane and the associated excess osmotic pressure. This direct interaction, however, fades away with the distance, and for non-overlapping EDLs the stiffness is due to the additional pressure gradient needed to compensate the flow induced by the gradient of chemical potential $\textrm{d}\mu _{\textit{s}}/\textrm{d}r$ . A primary diffusio–osmotic flow is $-T_{\textit{vs}}\textrm{d}\mu _{\textit{s}}/\textrm{d}r$ and always descends along the concentration gradient, because $T_{\textit{vs}}$ is always positive. But the concentration gradient also induces a diffusio–osmotic electrical current $I_{\textit{DO}}=-T_{\textit{es}} \textrm{d}\mu _{\textit{s}}/\textrm{d}r$ . The electric field generated to cancel this diffusio–osmotic current creates in turn a secondary electro–osmotic flow $-T_{\textit{ev}}/T_{\textit{ee}}I_{\textit{DO}}=T_{\textit{ev}}T_{\textit{es}}/T_{\textit{ee}}\textrm{d}\mu _{\textit{s}}/\textrm{d}r$ . As the product $T_{\textit{es}}T_{\textit{ev}}$ is always positive, this secondary electro–osmotic flow opposes the primary diffusio–osmotic flow and climbs up the concentration gradient. The relative importance of these two opposite contributions is responsible for the shoulders appearing in the variations of the stiffness with the distance in figure 6.

The main feature of the DK stiffness is its algebraic decay with the distance when EDLs do not overlap, qualitatively different from the DLVO stiffness which decays exponentially. This property arises from the affine variation of the transport coefficients with the film thickness (except for the film permeability, which scales as $h^3$ ). The accumulation of excess ions at large distance has also an affine variation with $h$ : $\textrm {d}N_{\textit{s}} \simeq h (n-n_0)$ , corresponding to $c_{h} (\tilde h) \simeq 0$ and $c_{n} (\tilde h) \simeq 2 \tilde h$ . At very large distance the $T_{\textit{a}}$ function (3.24b ) essentially reduces to $T_{\textit{a}} \simeq T_{\textit{vs}}/T_{\textit{vv}}$ and the DK stiffness decreases as $\tilde D^{-4}$ :

(6.3) \begin{equation} \text{Re} (\mathcal Z_{\textit{DK}}^{hf}) \simeq \frac {42 a_{\textit{vs}}^2}{5\tilde D^4}, \quad a_{\textit{vs}}=-\frac {1}{2} \ln (1-\gamma ^2), \quad \text{Re} (Z_{\textit{DK}}^{hf}) \simeq \frac {336 \pi R\ln ^2 (1-\gamma ^2)\varepsilon _{\textit{l}}^2(k_{\textit{B}}T)^3}{5 e^4 n_0 D^4}. \end{equation}

This large distance behaviour depends only on the electrical characteristic of the surface and does not involve any diffusion process. At short distance, the stiffness saturates to a value which is also independent of diffusion (see figure 6 a). Overall, the influence of diffusion in the high-frequency stiffness is very limited and appears only in an intermediate distance range.

Figure 7.

(a) Low frequency DK damping factor. Solid lines: $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ (4.14a ) as a function of $\tilde D$ , for $\alpha =100$ and various values of $K$ ; dashed lines: EK damping factor $\tilde f_{\textit{EK}}$ in the same conditions. (b) Solid lines: factor $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ for $K=0.1$ and various values of $\alpha$ . Dashed lines: factor $\tilde f_{\textit{EK}}$ in the same conditions. (c) Maximum value of $K\boldsymbol{\cdot }\tilde f_{\textit{DK}}^{\ell \textit{f}}$ over the distance, as a function of $\alpha$ , for various values of $K$ .

Figure 7. Long description

Panel A: A line graph shows the low-frequency limit of the DK damping factor and the EK damping factor as functions of the reduced distance for a fixed alpha value and different values of K. The x-axis represents the reduced distance and the y-axis represents the damping factors. Solid lines indicate the DK damping factor, while dashed lines indicate the EK damping factor, the two in log scales. Panel B: Another line graph displays the damping factors for a fixed K value and various values of alpha. The x-axis represents the reduced distance and the y-axis represents the factor. Solid lines show the factor, and dashed lines show the damping factors in the same conditions, the two in log scales. Panel C: A line graph illustrates the maximum value of the DK damping factor over the distance as a function of alpha for various values of K. The x-axis represents alpha, and the y-axis represents the maximum value of the factor, the two in log scales.

It is of interest to compare the DK stiffness with the DLVO stiffness. As quoted in (4.14a ), the dimensional prefactors are the same for both stiffness and the comparison involves essentially the distance $\tilde D$ and the parameter $\alpha$ , independently of the sphere radius. One sees on figure 6(b) that at distances sufficiently large the DK stiffness always exceeds the DLVO stiffness. Furthermore, the inversion of magnitude occurs essentially at a distance where the high-frequency DK stiffness does not depend on $K$ . For $\alpha =1$ the crossing distance is approximately $\tilde D=8$ , at a very low value of the stiffness. However, for higher values of $\alpha$ the DK stiffness exceeds the DLVO stiffness at distance $\tilde D \le 3.5 - 4$ , for stiffness value $\sim 0.05\times 8\pi Rn_0k_{\textit{B}}T$ at $\alpha \gt 20$ .

6.3. Low-frequency DK damping

Figure 7 represents the low-frequency DK damping factor $\displaystyle \tilde f_{\textit{DK}}^{lf}= ({\tilde D}/{6}) \beta ^{-1} \text{Im} (\mathcal Z_{\textit{DK}}^{\ell \textit{f}})$ as a function of $\tilde D$ for various values of $K$ and $\alpha$ . The EK damping factor $\tilde f_{\textit{EK}}$ is plotted in dashed lines for comparison.

For all values of $\alpha$ and $K$ , $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ exhibits a similar variation with $\tilde D$ : a linear increase at small $\tilde D$ , a maximum (or two maxima) around $\tilde D \sim 1$ , and a severe decay at $\tilde D \ge 5$ . The DK damping is in general significantly lower than the EK damping $\tilde f_{\textit{EK}}$ . However, in the small and intermediate distance range $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ increases proportionally to the inverse of the ions mobility $1/K$ (see figure 7 a) without saturation at small $K$ , unlike the EK damping. Therefore, at small $K$ values, $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ exceeds significantly $\tilde f_{\textit{EK}}$ in the distance range $1 \le \tilde D \le 10$ . For very small values $K \le 0.01$ , $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ may even exceed $\tilde f_{\textit{EK}}$ by one order of magnitude and become comparable to the Reynolds damping at $\tilde D \approx 3$ .

As the EK damping, the DK damping shows a non-monotonic behaviour with the surface charge parameter $\alpha$ (see figure 7 b), but this non-monotonic variation is restricted to the small distance range $\tilde D \le 2$ . In this small distance range, we find that $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ increases as $\alpha ^2$ at small $\alpha$ , goes through a maximum around $\alpha =10$ , and then decreases when $\alpha$ is further increased. For the distance range $\tilde D \gt 1$ the trend is reversed, and $\tilde f_{\textit{DK}}^{\ell \textit{f}}$ resumes a monotonic increase with $\alpha$ . It is remarkable that for $\alpha \ge 3$ , all $\tilde f_{\textit{DK}}^{\ell \textit{f}} (\tilde D)$ curves have a common point around $\tilde D \approx 1.5$ , whatever the value of $\alpha$ . We find that the amplitude of the DK damping factor at this common point is around $10^{-2}/K$ .

Finally, the maximum value of $K \times \tilde f_{\textit{DK}}^{\ \ell f}$ over all distances is shown in figure 7(c) as a function of $\alpha$ . For moderate and high values of $\alpha \in [3,1000]$ , the maximum DK damping factor lies in a narrow interval $[10^{-2}K^{-1}, 3. 10^{-2}K^{-1}]$ . The peak value is found around $\tilde D \approx 0.5$ for $\alpha \in [3,30]$ and around $\tilde D \approx 3$ for $\alpha \ge 30$ .