2. Governing equations of travelling wave electroosmosis

Consider the application of travelling wave charges

(2.1) \begin{equation} \sigma (x,t) = \sigma _0 e^{i (kx -\omega t)} \end{equation}

on an insulating wall, temporarily identified with the plane $z=0$ , cf. figure 1, with real frequency $\omega$ and real wave vector $\mathbf {k} = k \hat {\mathbf {x}}$ . We consider a $1:1$ electrolyte where each concentration species is $c_\pm =c_\infty + \delta c_\pm$ containing a perturbation $\delta c_\pm$ superposed on its uniform counterpart $c_\infty$ . Thus, the bulk charge distribution is $\rho = e(\delta c_+ - \delta c_-)$ while the salt distribution is $ s\sim 2 ec_\infty$ , to leading order, where $e$ is the proton charge. This implies that $s$ is also accompanied by a perturbation $e(\delta c_+ + \delta c_-)$ , which however is never invoked in this paper but it can be computed, need be. The validity of this approximation can only be examined a posteriori and in conjunction to the reductive form of the Nernst–Planck equation, as detailed below. We carry-out such an analysis in Appendix A.1.

Invoking the aforementioned approximation and the Poisson equation for the electric potential $\phi$ ,

(2.2) \begin{equation} \nabla ^2 \phi = - \frac {\rho }{\epsilon }, \end{equation}

where $\epsilon$ is the dielectric constant, the evolution of charge distribution,

(2.3) \begin{equation} \partial _t \rho = D\left [ \nabla ^2 \rho + \frac {e}{k_BT} \nabla \cdot \left ( s \nabla \phi \right ) \right ] - \mathbf {v} \cdot \nabla \rho , \end{equation}

reduces to the Debye–Falkenhagen equation (Bazant et al. Reference Bazant, Thornton and Ajdari2004)

(2.4) \begin{equation} \partial _t \rho = D\left [ \nabla ^2 \rho -\kappa ^2 \rho \right ] -\partial _x \rho \partial _z \psi + \partial _z\rho \partial _x \psi . \end{equation}

Here, $\mathbf {v}$ is the electrolyte velocity, $D$ is the charge diffusion coefficient (assumed to have the same value for both species), $k_BT$ the thermal energy and $\psi = \psi (x,z, t)$ is the streamfunction for a two-dimensional incompressible liquid, whereby

(2.5) \begin{equation} \mathbf {v} = (u, 0, w) = \left (\frac {\partial \psi }{\partial z}, 0, -\frac {\partial \psi }{\partial x}\right ), \end{equation}

and

(2.6) \begin{equation} \kappa = \left [ \frac {2e^2 c_\infty }{\epsilon k_BT} \right ]^{\frac {1}{2}} \end{equation}

is the inverse Debye length. How this approximation relates to its classical Debye–Hückel counterpart (where $\rho = -\epsilon \kappa ^2 \phi$ everywhere) is discussed in Appendix A.2.

Even in the absence of the last two (nonlinear advective) terms, (2.4) leads to a non-equilibrium charge distribution $\rho$ , where the charge changes both in time and in space, in contrast to the essentially equilibrium route taken in the literature (cf. Probstein Reference Probstein1994). Equation (2.4) is not new and has been derived earlier, for instance, by Cahill et al. (Reference Cahill, Heyderman, Gobrecht and Stemmer2004) and Mortensen et al. (Reference Mortensen, Olesen, Belmon and Bruus2005) following the older results of Ehrlich & Melcher (Reference Ehrlich and Melcher1982).

In the main body of this paper, we will adopt the low-Péclet-number limit, thus reducing (2.4) to

(2.7) \begin{equation} \partial _t \rho = D\left [ \nabla ^2 \rho -\kappa ^2 \rho \right ] . \end{equation}

We justify the small-Péclet-number approximation in Appendix B.1. We show in Appendix B.2 that retaining the advection terms in (2.4), based on the zero velocity mode only, does not lead to a significant change in the observables for moderate values of the Péclet number.

Finally, the streamfunction $\psi$ satisfies

(2.8) \begin{equation} \rho _l\partial _t \nabla ^2 \psi = \eta \nabla ^4 \psi + \partial _x \rho \partial _z \phi - \partial _z\rho \partial _x \phi , \end{equation}

where $\rho _l$ is the liquid electrolyte mass density and $\psi (x, 0, t) = 0 = \partial _z \psi (x,0, t)$ (the no-slip boundary condition) assuming temporarilly that the solid boundary is identified with the plane $z=0$ . If the flow extends to infinity, then its velocity is considered to have a finite value there.

All results of the present section depend on the (weakly) nonlinear character of the nonlinear torque (last two terms) in (2.8). Although, in general, the commensurate velocity field vanishes when averaged over the period of oscillation of the applied electric field, there are circumstances where this is not so. This is due to the presence of a zero mode that arises due to constructive interference of the charge and potential excitations arising in (2.8). The nonlinear effect described in (2.8) vanishes if the charge and the potential do not vary in the direction parallel to the wall (the $x$ -direction). Likewise, in the $\omega \equiv 0$ case, the nonlinear torque in (2.8) vanishes, on account of the connection $\rho = -\epsilon \kappa ^2 \phi$ (the Debye–Hückel approximation), or equivalently, the electric body force is the gradient of a scalar field and thus only affects the pressure distribution.

2.1. Boundary conditions

In this paper, we will predominantly consider Neumann boundary conditions for the Poisson equation (2.2), at a solid wall carrying a surface charge $\sigma$ as in (2.1), where

(2.9) \begin{equation} \hat {\mathbf {n}} \cdot \nabla \phi = -\frac {\sigma }{\epsilon }, \quad \textrm {at a charged wall}, \end{equation}

with the unit normal vector $\hat {\mathbf {n}}$ pointing into the liquid. We will consider also the case where the current of bulk charge $\rho$ vanishes at wall, that is,

(2.10) \begin{equation} \mathbf {J}\cdot \hat {\mathbf {n}} \equiv - D\left [ \nabla \rho + \frac {e}{k_BT} s \nabla \phi \right ] \cdot \hat {\mathbf {n}}=0 , \quad \textrm {at a charged wall}. \end{equation}

It is more convenient to replace (2.10) with a simpler expression. The same approximation that allowed the reduction of (2.3) into (2.4) leads the vanishing current condition (2.10) at the wall to become a boundary condition for the bulk charge $\rho$ ,

(2.11) \begin{equation} \hat {\mathbf {n}} \cdot \nabla \rho = \kappa ^2 \sigma , \quad \textrm {at a charged wall}. \end{equation}

Finally, appropriate conditions must be set at infinity if the domain is unbounded.

We should stress here that the system of equations (the Poisson–Nernst–Planck system) (2.2) and (2.4) with the Neumann boundary conditions does not have a solution when $k\equiv 0$ . Mathematically, the electric potential cannot satisfy all boundary conditions. However, the devices we have in mind are, say, at most $10$ cm long, with a wavenumber $k=10 \textrm { m}^{-1}$ . Past experiments on travelling wave electroosmosis employed electrode arrays with $k\sim 10^4 \textrm { m}^{-1}$ (Ramos et al. Reference Ramos, Morgan, Green, González and Castellanos2005) or even larger (Cahill et al. Reference Cahill, Heyderman, Gobrecht and Stemmer2004). We will thus operate within the interval of physically realistic wavenumbers $k\in (10^1, 10^5) \textrm { m}^{-1}$ , as we stated in table 1.

Table 1. Definitions of wavenumbers and parameter values. $\rho , k, K, \sigma$ and $\kappa$ as in Ajdari (Reference Ajdari1995). Our $k$ and $K$ correspond to the $q$ and $Q$ of Ajdari (Reference Ajdari1995). $p$ and ${\unicode{x1D6E5}}$ correspond to the $k$ and $\delta$ of Landau & Lifshitz (Reference Landau and Lifshitz1987, § 24) (replacing $D$ with $\nu$ ).

The notation used in (2.8) implies that $\rho$ and $\phi$ are real fields. In the subsequent discussion, we will employ the same notation to denote complex fields from now on.

2.2. Velocity scales

We will employ the following velocity scales:

(2.12) \begin{equation} u_0 = \frac {1}{2}\frac {\sigma _0^2}{\epsilon \eta \kappa }, \quad u_1 = \frac {\epsilon \kappa \phi _0^2}{2\eta }, \end{equation}

and, where appropriate, define $E \equiv {\sigma _0}/{\epsilon }$ to be the nominal electric field amplitude due to the interfacial charge distribution $\sigma _0$ (Neumann problem). Here, $\phi _0$ is a characteristic scale for the electric potential (Dirichlet problem). Dimensionless groups will be stated in raw form to avoid introducing new notation.

For brevity, we will use the phrases ‘Neumann velocities’ or ‘Dirichlet velocities’ to denote velocities obtained by solving the corresponding Poisson equation with Neumann or Dirichlet boundary conditions, respectively.

3. Boundary-driven flows in a semi-infinite space

We consider the configuration displayed in figure 1. Travelling wave surface charges $ \sigma (x,t) = \sigma _0 e^{i (kx -\omega t)}$ are applied on the wall $z=0$ bounding a $1:1$ electrolyte lying in the semi-infinite space $z\gt 0$ . The boundary conditions satisfied by the bulk charge $\rho$ and electric potential $\phi$ at a solid surface lying at $z=0$ with the unit vector $\hat {\mathbf {n}}=\hat {\mathbf {z}}$ are, from (2.9) and (2.11),

(3.1) \begin{equation} \partial _z \phi (x,z = 0, t) = - \frac {\sigma (x, t)}{\epsilon }, \quad \partial _z \rho (x,z = 0, t) = \kappa ^2 \sigma (x, t). \end{equation}

Assuming $\rho = \rho (z) e^{i (kx -\omega t)}$ , (2.7) reduces to $ \rho _{zz} + [i\omega /D -\kappa ^2 -k^2 ] \rho =0$ . To avoid clutter, we have introduced the notation $\rho _z\equiv \partial _z \rho$ etc. Thus, the charge distribution reads

(3.2) \begin{equation} \rho (x,z,t) = -\frac {i\sigma _0\kappa ^2}{P} e^{i(Pz +kx -\omega t)}, \end{equation}

where we assumed a vanishing bulk charge at infinity and introduced the notation

(3.3) \begin{equation} P \equiv P_1+iP_2 = \sqrt {p^2 -k^2 -\kappa ^2}, \quad i\omega = D p^2, \quad p = \frac {1+i}{{\unicode{x1D6E5}}}, \quad {\unicode{x1D6E5}} = \sqrt {\frac {2D}{\omega }}. \end{equation}

The notation employed in (3.3) for penetration depth ${\unicode{x1D6E5}}$ and complex wavenumber $p$ is analogous to that employed by Landau & Lifshitz (Reference Landau and Lifshitz1987, p. 84). These quantities here, however, refer to the wave-like form of the charge distribution away from a charged boundary, rather than the wavelike form of the velocity field away from a no-slip wall. The penetration depth is now determined by the inverse of the imaginary part of $P$ (taken here to be positive) and there is an oscillation whose wavevector is the real part of $P$ . In other words, if $P_1$ and $P_2$ denote the real and imaginary parts of $P$ in (3.3) ( $P = P_1+iP_2$ ), where $P_i$ are real and $P_2\gt 0$ , then

(3.4) \begin{equation} P_{1,2} = \frac {1}{\sqrt {2}} \sqrt {\sqrt {\left (\kappa ^{2}+k^2\right )^{2}+\frac {\omega ^{2}}{\textit {D}^{2}}}\mp (\kappa ^{2}+k^2)}, \end{equation}

where the minus/plus sign corresponds to the real/imaginary part of $P$ . Thus, the penetration depth is not determined by $2\pi /\kappa$ but by the length scale $2\pi /P_2\gt 0$ . With this notation, the charge distribution has the form

(3.5) \begin{equation} \rho (x,z,t) =- \frac {i\sigma _0\kappa ^2}{P_1+iP_2} e^{-P_2z}e^{i(P_1 z +kx -\omega t)}. \end{equation}

Thus, in addition to the exponential decay away from the wall, there is a plane wave whose wavevector is not parallel to the surface charge distribution wave vector (it is not parallel to the $x$ -axis) in (2.1) but lies in the direction $k \hat {\mathbf {x}} + P_1 \hat {\mathbf {z}}$ .

Similarly, assuming $\phi = \phi (z) e^{i (kx -\omega t)}$ , (2.2) reduces to $ \phi _{zz} -k^2 \phi = ({i\sigma _0\kappa ^2}/{\epsilon P}) e^{iPz},$ subject to the boundary condition (3.1) with surface charge (2.1) and $\phi =0$ at infinity (the alternative boundary condition of zero electric field at infinity leads to an identical result). Thus, the potential distribution reads

(3.6) \begin{equation} \phi (x,z,t) = \frac {\sigma _0}{\epsilon } \left ( 1 + \frac {\kappa ^2}{P^2 +k^2} \right ) \frac {e^{-|k|z}}{|k|} e^{i(kx -\omega t)} +\frac {1}{\epsilon (P^2 +k^2)} \rho , \end{equation}

with $\rho$ given by (3.5). Note that, as we mentioned in the previous section, setting $k=0$ , the boundary value problem does not have a solution since the electric potential $\phi$ cannot satisfy both boundary conditions. This is then reflected in the divergence of expression (3.6) as $k\rightarrow 0$ . For the purposes of this paper, however, where the physically realistic wavenumber $k$ lies in the interval $(10^1, 10^5) \textrm { m}^{-1}$ , this limit is never attained.

Also, see Appendix A.2 for a comparison between the result (3.6) and its Debye–Hückel counterpart.

It is now clear that the nonlinear torque (last two terms in (2.8)) involves the harmonics $e^{\pm i \theta }$ and $e^{\pm 2i\theta }$ , where $\theta = kx-\omega t$ , and the commensurate streamfunction $\psi$ satisfying (2.8) will also be composed of the same harmonics. Thus, the velocity vanishes when averaged over the period $2\pi /\omega$ of oscillations.

There is, however, a part of the nonlinear torque in (2.8) that is time-independent and does not vanish when averaged over the period of oscillations. This is caused by constructive interference of the various harmonics expressing the charge and potential distribution in (3.5) and (3.6), respectively, and gives rise to unidirectional pumping of liquid, its direction determined by the propagation direction of the charge distribution lying on the walls. We now investigate this zero mode.

Let $\psi = \psi (z)$ , and thus the velocity $ \mathbf {v} = u(z) \hat {\mathbf {x}}$ as in (2.5), and consider both $\psi$ and $u$ to be real fields, cf. figure 1. In Appendix H.1, we show that after averaging over the period of oscillation, the vorticity equation (2.8) reduces to

(3.7) \begin{equation} \partial _zu(z) = \frac {i\epsilon k}{4\eta } \left ( \phi ^*\phi _z - \phi \phi _z^*\right ), \end{equation}

where a star on $\phi$ denotes the complex conjugate of (3.6). It is clear that the complex exponentials $e^{i(kx -\omega t)}$ have cancelled out. This is the central formula of this paper. It states that the shear stress in a viscous electrolyte is given by the expression on the right-hand side (multiplied by viscosity) which resembles a probability current for a Schrödinger equation (Morse & Feshbach Reference Morse and Feshbach1953). We also verify the validity of (3.7) by starting directly from the time-dependent Stokes equations, see Appendix H.2.

Figure 2. Boundary layer of thickness (= penetration depth, see (3.3)) of $O(P_2^{-1})\sim \kappa ^{-1}$ for large $\kappa$ employing (3.9). Beyond the boundary layer, the liquid velocity rapidly acquires a constant value. $\kappa ^{-1}, D, \eta$ are the Debye length, charge diffusion coefficient and liquid viscosity, respectively. Here, we have taken $k =10^5 \textrm { m}^{-1}$ , $\kappa =10^7 \textrm { m}^{-1}$ $D = 10^{-9}\textrm { m}^{2}$ s⁻¹ and $\omega = D\kappa ^2$ . $u_0$ was defined in (2.12).

It is thus easy to determine the horizontal velocity component $u(z)$ . Employing the complex form of the field $\phi$ in (3.6), we solve (3.7) subject to the boundary conditions

(3.8) \begin{equation} u(z=0)=0, \quad u(z=\infty ) = \textrm {finite}. \end{equation}

Thus, the zero mode velocity $u(z)$ from (3.7) with (3.6) becomes

(3.9) \begin{equation} u(z) = \frac {i\epsilon k}{4\eta } \int _0^z \left ( \phi ^*\phi _z - \phi \phi _z^*\right ) {\textrm{d}}z. \end{equation}

The zero mode (3.9) inherits the boundary-layer structure of the field $\phi$ as displayed in figure 2. A boundary layer of thickness $O(P_2^{-1})\sim \kappa ^{-1}$ for large $\kappa$ separates the no-slip region at the wall to the constant value acquired by the velocity away from the wall.

In the semi-infinite space formulation of the present section, the observable of interest is the value taken by the electrolyte velocity zero mode (3.9) at $z=\infty$ . It reads

(3.10) \begin{equation} \frac {u(\infty )}{u_0} = \pm \frac {\kappa ^{3} \left (2 R^{3} \cos \left (3 \Theta \right )-4 \sin \! \left (2 \Theta \right ) R^{2} |k| +2 R (\kappa ^{2} - k^2) \cos \left (\Theta \right )+\kappa ^{2} |k| \cot \! \left (\Theta \right )\right )}{2 R^{2} \left (2 R^{2} k^{2} \cos \! \left (2 \Theta \right )-4R|k| (R^2 + k^2) \sin \! \left (\Theta \right ) -R^{4}-4 R^{2} k^{2}-k^{4}\right )}, \end{equation}

where the plus sign gives the velocity for positive $k$ and the negative sign for negative $k$ and, following Ehrlich & Melcher (Reference Ehrlich and Melcher1982), the horizontal velocity $u$ was scaled by $u_0$ , cf. (2.12). Here, $R$ and $\Theta$ are the amplitude and phase of the complex wavenumber $P = R e^{i\Theta }$ defined in (3.3) and given explicitly by

(3.11) \begin{equation} R = {\left [\left (\kappa ^{2}+k^2\right )^{2}+\frac {\omega ^{2}}{\textit {D}^{2}}\right ]}^{\frac {1}{4}}, \quad \Theta = -\frac {1}{2}\arctan \! \left [\frac {\omega }{\textit {D} \left (\kappa ^{2}+k^2\right )}\right ]+\frac {\pi }{2}. \end{equation}

In Appendix C, we discuss the relation of the velocity (3.10) with its counterpart that appeared in the past literature (Ehrlich & Melcher Reference Ehrlich and Melcher1982).

Figure 3. Horizontal zero-mode velocity evaluated at $z= \infty$ , from (3.9) (a) as a function of dimensionless frequencies of the travelling wave wall charge and (b) as a function of the travelling electric potential at the wall, the latter was derived in Appendix F. The curves are self-similar as described by (3.12) and (3.13), based however on different dimensionless groups. Introducing the scaling that appears in the label of the vertical axis of panel (b), the Dirichlet results centre about the (continuous) curve introduced by Ramos et al. (Reference Ramos, Morgan, Green, González and Castellanos2005), (10); (see (3.14)). The different $k$ curves in panel (a) centre about the Ehrlich & Melcher (Reference Ehrlich and Melcher1982) expression, cf. (C2). If the value of $k$ becomes sufficiently large, say $k\gt 10^5 \textrm { m}^{-1}$ , the amplitude of the curves in both panels decreases.

3.1. Self-similar behaviour of the slip velocity

The zero mode slip velocity (3.10) displays a self-similar behaviour for the dimensions of the parameters $\kappa , k, D, \eta , \omega$ define a matrix of rank three. Thus, there are two dimensionless combinations that can be formed (Bluman & Kumei Reference Bluman and Kumei1989; Panton Reference Panton1996). The average velocity (3.10) can then be expressed in the scaling form

(3.12) \begin{equation} u(z=\infty )= \frac {1}{2}\frac {\sigma _0^2}{\eta \epsilon \kappa } F\left (\frac {\omega }{D\kappa ^2}, \frac {k}{\kappa }\right ), \end{equation}

where $F$ is a nonlinear function of the two independent dimensionless parameters and the front factor is the velocity scale $u_0$ (2.12), as this is displayed in the vertical axis of figure 3.

In figure 3(a) we plot the ‘slip’ velocity (3.9), that is, the horizontal zero-mode velocity $u(z)$ evaluated at $z= \infty$ for various values of $k$ . It is thus seen that all curves with $k$ up to say $10^5 \textrm { m}^{-1}$ collapse into the single (continuous) curve introduced by Ehrlich & Melcher (Reference Ehrlich and Melcher1982), cf. (C2).

For comparison with results that have appeared in the more recent literature, infigure 3(b), we plot the velocity (3.9), but for Dirichlet boundary conditions, employing the potentials (F3), whose derivation was carried out in Appendix F. We thus display the slip velocity also in scaling form, although the dimensionless groups are now different,

(3.13) \begin{equation} u(z=\infty ) = \frac {1}{2}\frac {\epsilon k \phi _0^2}{\eta } G\left (\frac {\omega }{Dk\kappa }, \frac {k}{\kappa }\right ), \end{equation}

where $G$ is a nonlinear function of the two dimensionless variables and the front factor is the velocity scale $ku_1/\kappa$ (2.12), as this is displayed in the label of the vertical axis of figure 3(b). The curves appearing in figure 3(b) are thus the function $G$ where we only vary its first argument. Notice that, had we employed $\omega /(D\kappa ^2)$ instead of $\omega /(Dk\kappa )$ to non-dimensionalise $u(\infty )$ in (3.13), the curves in figure 3(b) for each different $k$ would be translated to the left and to the right of the centre curve.

All the curves in panel (b) collapse to the slip velocity (continuous curve) determined by Ramos et al. (Reference Ramos, Morgan, Green, González and Castellanos2005), figure 10, having the form

(3.14) \begin{equation} u_{\textrm {Ramos}} = \frac {\kappa u_1}{k} \frac {\Omega }{1+ \Omega ^2}, \end{equation}

where $\Omega = {\omega }/(D k \kappa)$ in our notation.

The self-similar behaviour of the velocity is important as it implies that results obtained with one experimental configuration can be employed to obtain estimates of results corresponding to a different experimental configuration without the need to repeat the experiment.

3.2. Estimates of $u(\infty )$

The preceding discussion, although informative, it does not provide an understanding of the magnitude of the effects, which we thus analyse in the present section.

Let us obtain an estimate of the velocity $u(\infty )$ in each case by considering the experimental values of Ramos et al. (Reference Ramos, Morgan, Green, González and Castellanos2005) and García-Sánchez et al. (Reference García-Sánchez, Ramos, Green and Morgan2006), where $k\sim 10^4 \textrm { m}^{-1}$ , voltage $\phi _0\sim 3$ V and employing standard parameter values for water $\epsilon = 7\times 10^{-10}$ F m⁻¹, $\eta = 0.001$ kg m⁻¹ s, $\kappa = 10^{7}\textrm { m}^{-1}$ .

A fair comparison between figures 3(a) and 3(b) can be furnished by considering the dimensionless value of the velocity at the peak of the graph in each panel, which occurs at different frequencies. Thus, in the Dirichlet case, we consider $\omega \sim Dk\kappa$ , thus giving the peak dimensionless velocity value $\sim 0.5$ in figure 3(b), leading its dimensional counterpart to be equal to

(3.15) \begin{equation} u(\infty ) = 1.57\; \textrm {cm} \, \textrm{s}^{-1} \quad \textrm {(Dirichlet)}, \end{equation}

which is rather high. However, for the same values as above but adopting $E =\sigma _0/ \epsilon = 10^5$ V m⁻¹, as done by Ehrlich & Melcher (Reference Ehrlich and Melcher1982), $\omega \sim D\kappa ^2$ , thus giving the peak dimensionless velocity value of ${\sim} 0.405$ in figure 3(a), leading its dimensional counterpart to be equal to

(3.16) \begin{equation} u(\infty ) = 141 \; \unicode{x03BC} \textrm {m}\,\textrm{s}^{-1} \quad \textrm {(Neumann)}, \end{equation}

which is within the order of magnitude of experimental velocities (cf. figure 7 of García-Sánchez et al. (Reference García-Sánchez, Ramos, Green and Morgan2006), but also figures 2 and 10 of Cahill et al. (Reference Cahill, Heyderman, Gobrecht and Stemmer2004)). Note that the electric field varies on the length scale $1/k$ , cf. (3.6), justifying the choice $E = 10^5$ V m⁻¹ for $k\sim 10^4 \textrm { m}^{-1}$ , in determining the estimate (3.16) for compatibility with the electric potential value $\phi _0\sim 3$ V employed in (3.15).

A recurrent theme in travelling wave electroosmosis studies is the high theoretical estimates (usually carried out by employing Dirichlet boundary conditions) in comparison to much lower, by a few orders of magnitude, experimental measurements (cf. Cahill et al. Reference Cahill, Heyderman, Gobrecht and Stemmer2004; Ramos et al. Reference Ramos, Morgan, Green, González and Castellanos2005). It is thus possible that the Neumann boundary conditions employed in this paper and their accompanying estimates, as in (3.16), may resolve this inconsistency.

3.3. Tails of $u(\infty )$

Both Neumann and Dirichlet tails in figure 3 are of $O(\omega )$ as $\omega \rightarrow 0$ :

(3.17) \begin{equation} u(\infty ) \sim \frac {\kappa \omega }{4D} \frac {k(7k^2 + 8\kappa ^2) + \sqrt {k^2 + \kappa ^2}(8k^2 + 4\kappa ^2)}{ \left ( \kappa ^2 + 2 k^2 + 2k \sqrt {k^2 + \kappa ^2} \right )^2 } \times \left \{ \begin{array}{ll} \frac {u_1}{k^2 + \kappa ^2} & \textrm {(Dirichlet)}, \\ \frac {u_0\kappa ^2}{(k^2 + \kappa ^2)^2 } & \textrm {(Neumann)}, \end{array} \right .\ \!\textrm {as} \!\!\quad \omega \rightarrow 0. \end{equation}

The high-frequency asymptotics of the exact expressions (3.12) and (3.13) are given by

(3.18) \begin{equation} u(\infty ) \sim \left \{ \begin{array}{ll} \frac {u_1}{\sqrt {2}} \left ( \frac {k}{\kappa }\right )^{3/2} \left (\frac {\omega }{Dk\kappa }\right )^{-3/2} & \textrm {(Dirichlet)}, \\[5pt] \frac {u_0}{\sqrt {2}} \left (\frac {\omega }{D\kappa ^2}\right )^{-3/2} & \textrm {(Neumann)}, \end{array} \right .\ \textrm {as} \quad \omega \rightarrow \infty . \end{equation}

Notice that the inherent frequency scalings in (3.18) agree with the scalings of the functions $G$ and $F$ in (3.13) and (3.12), respectively.

Figure 4. Travelling wave wall charges give rise to a nonlinear body force and torque (see (2.8)) in a $1:1$ electrolyte in a rectangular channel of width $2h$ leading to the appearance of a (zero mode) unidirectional fluid velocity in the $\hat {\mathbf {x}}$ direction, parallel to the channel walls that is quadratic with respect to the associated electric field and that does not vanish after averaging over the charge period of oscillation.

4. Boundary-driven flows in a channel of width $2h$

It is more physically realistic to consider the configuration displayed in figure 4. Travelling wave surface charges are applied on the channel walls at $z=\pm h$ enclosing a $1:1$ electrolyte. The boundary conditions satisfied by the potential $\phi$ and by the charge $\rho$ are

(4.1) \begin{equation} \partial _z \phi (x,z = \pm h, t) = \pm \frac {\sigma (x, t)}{\epsilon }, \quad \partial _z \rho (x,z = \pm h, t) =\mp \kappa ^2 \sigma (x, t), \end{equation}

with $\sigma = \sigma _0 e^{i (kx -\omega t)}$ . The choice of this set-up was dictated by results showing that the ideal and optimal configuration seems to require symmetrical electrode arrays with respect to the channel centre and with no phase lag between them (Yeh et al. Reference Yeh, Yang and Luo2011).

Assuming $\rho = \rho (z) e^{i (kx -\omega t)}$ and subject to the boundary condition (4.1) with surface charge (2.1), the equations to solve are identical to those of the semi-infinite space. Thus, the charge distribution reads

(4.2) \begin{equation} \rho (x,z,t) = \frac {\sigma _0\kappa ^2\cos Pz}{P\sin Ph} e^{i(kx -\omega t)} , \end{equation}

where $P$ is the complex wavenumber defined in (3.3).

Similarly, assuming $\phi = \phi (z) e^{i (kx -\omega t)}$ and subject to the boundary condition (4.1) with surface charge (2.1), the potential distribution reads

(4.3) \begin{equation} \phi (x,z,t) = \frac {\sigma _0}{\epsilon } \left ( 1 + \frac {\kappa ^2}{P^2 +k^2} \right ) \frac {\cosh kz}{k\sinh kh} e^{i(kx -\omega t)} +\frac {1}{\epsilon (P^2 +k^2)} \rho , \end{equation}

with $\rho$ given by (4.2). We note that, as is the case in the semi-infinite space, the above equations do not have a solution when $k=0$ . For the same reasons as before, we are not interested in this limit as the physically realistic wall charge excitation wavenumbers we have in mind are very large, lying in the interval $(10^1, 10^5) \textrm { m}^{-1}$ , as is stated in table 1.

The zero-mode velocity in the channel again satisfies the integrated momentum (3.7). Employing the complex form of the field $\phi$ in (4.3), we thus solve (3.7) subject to the no-slip boundary conditions

(4.4) \begin{equation} u(z=\pm h)=0. \end{equation}

The zero-mode velocity $u(z)$ from (3.7) with (4.3) becomes

(4.5) \begin{equation} u(z) = \frac {i\epsilon k}{4\eta } \int _h^z \left ( \phi ^*\phi _z - \phi \phi _z^*\right ) {\textrm{d}}z .\end{equation}

This is essentially a plug-like flow with very thin boundary layers of thickness ${\sim} \kappa ^{-1}$ located at $z =\pm h$ . The exact form of (4.5) is delegated to Appendix E (E3). The meaningful observable here is the average value of the zero mode over the width of the channel

(4.6) \begin{equation} \langle u \rangle = \frac {1}{2h} \int _{-h}^h {u(z)} {\textrm{d}}z. \end{equation}

As it stands, (4.6) is a double integral and could be awkward to evaluate, especially in cylindrical polars (see § 5). It can be significantly simplified by changing the order of integration and carrying out one of the integrations. It can thus be written as a single integral in the form

(4.7) \begin{equation} \langle u \rangle = -\frac {1}{2h} \frac {i\epsilon k}{4\eta }\int _{-h}^h \left ( \phi ^*\phi _z - \phi \phi _z^*\right ) z{\textrm{d}}z, \end{equation}

where the parity of $\phi$ with respect to the origin $z=0$ was taken into account. The exact form of (4.6) is delegated to Appendix E (E5).

4.1. Self-similar behaviour of the average velocity in the channel

The zero-mode velocity (4.7) displays a self-similar behaviour for the dimensions of the available parameters $h, \kappa , k, D, \eta , \omega$ define a matrix of rank three. From the available dimensionless combinations (Bluman & Kumei Reference Bluman and Kumei1989; Panton Reference Panton1996), we display the following in the average velocity (4.6):

(4.8) \begin{equation} \langle u \rangle = \frac {1}{2}\frac {\sigma _0^2}{\eta \epsilon k} f\!\left (\frac {\omega }{D\kappa ^2}, \kappa h\right ), \end{equation}

where $f$ is a nonlinear function of the two independent dimensionless parameters and the front factor is the velocity scale $\kappa u_0/k$ , cf. (2.12), as this is displayed in the label of the vertical axis of figure 5 (the scaling theory allows other combinations as well, for instance, see (4.14) below). The front factor in (4.8) implies that $\langle u \rangle$ scales as $k^{-1}$ . This behaviour is clearly displayed by the continuous curve in figure 14(b) of Appendix G. Such a behaviour for Neumann boundary conditions was met before (Liu et al. Reference Liu, Ren, Tao, Li and Wu2018, figure 6f). Therein, increasing the size of the system gave rise to a commensurate increase of the Coulomb force and the resulting slip velocity. The curve in figure 6(f) of this reference has exactly slope equal to $1$ , that is, the slip velocity scales linearly with respect to system size, and is thus inversely proportional to $k$ , as this is demonstrated by the front factor in (4.8) and the continuous curve in figure 14(b).

Figure 5. Self-similar behaviour of the zero mode velocity (4.6), averaged over the channel width $2h$ as a function of dimensionless frequencies (a) of the travelling wave wall charge or (b) of the travelling electric potential at the wall, the latter was derived in Appendix F. The curves are self-similar as described by relations (4.8) and (4.10). We emphasise that the dimensionless group $\omega /(Dk^2)$ employed in panel (b) differs from the group $\omega /(D\kappa ^2)$ employed in panel (a), and furthermore, it is different to the group $\omega /(Dk\kappa )$ employed for the Dirichlet case in the semiinfinite space (figure 3 b). In both panels, $h = 10^{-5}$ m, $D = 10^{-9}\textrm { m}^{2}$ s⁻¹ and in panel (a), we employed the fit (4.9) with the single parameter $\beta = 0.926$ .

A good single-parameter fit of these results that accounts for all dimensionless frequencies, including those present in the tails, is given by the continuous line in figure 5(a) and has the functional form

(4.9) \begin{equation} \langle u \rangle _{E\&M} = \frac {\beta }{ \kappa h} \frac {\frac {\omega }{D\kappa ^2} }{\left [\left (\frac {\omega }{D\kappa ^2}\right )^{2}+1\right ]^{\frac {5}{4}}}\cos \! \left (\frac {\arctan \left (\frac {\omega }{D\kappa ^2}\right )}{2}\right ), \end{equation}

i.e. it is the velocity (C2) scaled accordingly and the single-parameter value $\beta = 0.926$ .

Thus, experimental results obtained by employing a single configuration that agree with the above behaviour can be extrapolated theoretically to fit alternative experimental conditions and parameters without the need of repeating the experiments.

Employing the Dirchlet boundary conditions instead, we obtain the potential (F6) and thus the velocity (4.7) acquires the self-similar form

(4.10) \begin{equation} \langle u \rangle = \frac {1}{2}\frac {\epsilon k \phi _0^2}{\eta } g\!\left (\frac {\omega }{Dk^2}, \kappa h\right ), \end{equation}

where $g$ is a nonlinear function of the two independent dimensionless parameters and the front factor is the velocity scale $k u_1/\kappa$ cf. (2.12), as this is displayed in the vertical axis of figure 5(b). Notice however that the dimensionless group $\omega /(Dk^2)$ employed in (4.10) (figure 5 b) differs from the group $\omega /(D\kappa ^2)$ employed in (4.8) (figure 5 a) and furthermore, it is different to the group $\omega /(Dk\kappa )$ employed for the Dirichlet case in the semiinfinite space in (3.13) (figure 3 b).

4.2. Average velocity $\langle u \rangle$ estimates

We estimate the channel average velocity by considering the same material parameters as in § 3.2. Comparison between figures 5(a) and 5(b) can be furnished by considering the dimensionless value of the velocity at the peak of the graph in each panel, which occurs at different frequencies. In the Dirichlet case, we consider $\omega \sim Dk^2$ , thus giving the peak dimensionless velocity value of ${\sim} 0.45$ in figure 5(b), leading its dimensional counterpart to equal

(4.11) \begin{equation} \langle u \rangle = 1.4\; \textrm {cm}\ \textrm{s}^{-1} \quad \textrm {(Dirichlet)}. \end{equation}

Figure 6. (a) Self-similar behaviour of average velocity by varying only the scaled channel width (second argument of (4.8)). Continuous line denotes the fit (4.9). (b) Averaged channel width zero mode velocity (4.7) versus channel height $h$ , obtained with Neumann or Dirichlet boundary conditions (employing the potentials (4.3) and (F6), respectively) and scaled by either $u_0$ , $u_1$ or $ku_1/\kappa$ (cf. (2.12)). The Neumann velocity reaches a plateau for small channel heights, while its Dirichlet counterpart decays to zero. For larger values of $h$ , all curves reach plateaus. Curves obtained with a frequency corresponding to the peak velocity of figures 5(a) and 5(b), respectively.

However, for the same values as above but adopting $E = 10^5$ V m⁻¹, as done by Ehrlich & Melcher (Reference Ehrlich and Melcher1982), $\omega \sim D\kappa ^2$ , thus giving the peak dimensionless velocity value of ${\sim} 0.037$ in figure 5(a), leading its dimensional counterpart to equal

(4.12) \begin{equation} \langle u \rangle = 1.3 \textrm { cm}\ \mathrm{s}^{-1} \quad \textrm {(Neumann)}, \end{equation}

which are comparable to each other. Notice however that the velocity scaling of the Dirichlet problem $k u_1/\kappa$ decreases linearly with increasing $k$ and that the Neumann scaling $\kappa u_0/k$ , increases, cf. labels of the vertical axes of figures 5(a) and 5(b), respectively.

4.3. Average velocity $\langle u \rangle$ asymptotic behaviour

In figure 6, we display the characteristic behaviour of the average velocity (4.7) as the height $h$ of the channel varies. The Neumann velocity in general increases as $h$ decreases reaching a plateau at small channel heights (where the Debye layers from adjacent walls start to overlap). The Dirichlet velocity decays to zero in this limit. These behaviours are described by the limiting expressions

(4.13) \begin{align} \langle u \rangle \sim \left \{ \begin{array}{ll} \frac {k u_1}{\kappa } \frac {\frac {\omega }{Dk^2}}{3\left [ 1 + \left ( \frac {\omega }{Dk^2} \right )^2 \right ]} \left ( \kappa h\right )^2 & \textrm {(Dirichlet)}, \\\notag \frac {\kappa u_0}{k}\frac {2}{3} \frac {-\kappa ^2 P_1P_2 \left [ P_1^2 +(P_2+k)^2 \right ]^3\left [ P_1^2 +(P_2-k)^2 \right ]^3 }{|P|^4\left [ P^2+k^2 \right ]^3\left [ (P^*)^2+k^2 \right ]^3} & \textrm {(Neumann),} \end{array} \right .\ \textrm {as} \quad h\rightarrow 0, \end{align}

where $P$ is defined in (3.3), and $u_0$ and $u_1$ in (2.12). Combination of this result with the $k$ scaling behaviour, displayed in figure 14, leads to the conclusion that the Neumann conditions provide prominent average velocities in thin and long channels (where the wavelength of the charge excitation at the walls can become sufficiently long).

Both Neumann and Dirichlet tails in figure 5 are of $O(\omega )$ (whose coefficients are too long to include here), as $\omega \rightarrow 0$ . The high-frequency asymptotics of the exact expressions (4.8) and (4.10) are given by

(4.14) \begin{equation} \langle u \rangle \sim \left \{ \begin{array}{ll} u_1 \frac {1}{\sqrt {2}} \left (\frac {\omega }{Dk^2}\right )^{-3/2} \tanh kh & \textrm {(Dirichlet)}, \\ \frac {u_0}{\sqrt {2}\tanh kh} \left (\frac {\omega }{D\kappa ^2}\right )^{-3/2} & \textrm {(Neumann)}, \end{array} \right .\ \textrm {as} \quad \omega \rightarrow \infty . \end{equation}

It is clear that the expressions in (4.14) justify the frequency scalings of the functions $g$ and $f$ in (4.10) and (4.8), respectively.

Figure 7. Travelling wave wall charges give rise to a nonlinear body force and torque (see (2.8)) in a $1:1$ electrolyte in a cylindrical capillary of radius $a$ leading to the appearance of a unidirectional fluid velocity in the $\hat {\mathbf {z}}$ direction, along the capillary centre axis, that is quadratic with respect to the associated electric field and that does not vanish after averaging over the charge period of oscillation.

5. Boundary-driven flows in a cylindrical capillary

A still more physically realistic configuration is displayed in figure 4. Travelling wave surface charges $\sigma = \sigma _0 e^{i(kz - \omega t)}$ are applied to the wall of a capillary with circular cross-section of radius $a$ , enclosing a $1:1$ electrolyte. The boundary conditions satisfied by the potential $\phi$ and by the charge $\rho$ are

(5.1) \begin{equation} \partial _r \phi (z,r =a, t) = \frac {\sigma (z, t)}{\epsilon }, \quad \partial _r \rho (z,r=a, t) = - \kappa ^2 \sigma (z, t), \end{equation}

and we take the axis of the cylinder to lie in the $z$ -direction as displayed in figure 7. Assuming $\rho = \rho (r) e^{i (kz -\omega t)}$ and $\phi = \phi (r) e^{i (kz -\omega t)}$ , the evolution equation for the charge and Gauss law reduce to

(5.2) \begin{equation} \rho _{rr} + \frac {1}{r} \rho _r + \left [ \frac {i\omega }{D} - k^2 - \kappa ^2 \right ] \rho = 0 \quad \textrm {and} \quad \phi _{rr} + \frac {1}{r} \phi _r - k^2 \phi = - \frac {\rho }{\epsilon }, \end{equation}

respectively. With boundary conditions (5.1), we obtain

(5.3) \begin{equation} \rho (z,r,t) = \frac {\sigma _0\kappa ^2J_0( Pr)}{PJ_1( Pa)} e^{i(kz -\omega t)} , \end{equation}

where $P$ is (again) the complex wavenumber (3.3) and $J_0$ , $J_1$ are Bessel functions of the first kind. Likewise, subject to the boundary condition (5.1), the potential distribution reads

(5.4) \begin{equation} \phi (z,r,t) = \frac {\sigma _0}{\epsilon } \left ( 1 + \frac {\kappa ^2}{P^2 +k^2} \right ) \frac {I_0(kr)}{k I_1( ka)} e^{i(kz -\omega t)} +\frac {1}{\epsilon (P^2 +k^2)} \rho , \end{equation}

with $\rho$ given by (5.3), and $I_0$ and $I_1$ are modified Bessel functions of the first kind (Bessel functions of imaginary argument).

As before, the velocity field consists of terms that vanish after averaging over the period of oscillation of the fields. The exception is the zero-mode velocity, satisfying an equation analogous to (3.7). Starting from the time-dependent Stokes equations and considering the velocity field to have the form $\mathbf {v} = u(r) \hat {\mathbf {z}}$ (as depicted in figure 7), we obtain

(5.5) \begin{equation} u_{rr} + \frac {1}{r} u_r + \frac {ik}{4\eta } (\rho \phi ^* - \rho ^* \phi ) =0. \end{equation}

Using the second of (5.2) and performing one integration, we arrive at

(5.6) \begin{equation} \partial _ru(r) = \frac {i\epsilon k}{4\eta } \left ( \phi ^*\phi _r - \phi \phi _r^*\right ) ,\end{equation}

which is the cylindrical counterpart of (3.7). Employing the complex form of the field $\phi$ in (5.4), we thus solve (5.6) subject to the no-slip boundary condition

(5.7) \begin{equation} u(r=a)=0. \end{equation}

The zero-mode velocity $u(r)$ from (5.6) with (5.4) becomes symbolically

(5.8) \begin{equation} u(r) = \frac {i\epsilon k}{4\eta } \int _a^r \left ( \phi ^*\phi _r - \phi \phi _r^*\right ){\textrm{d}}r, \end{equation}

which is the cylindrical counterpart of (3.9) and (4.5) (the constant of integration is zero to maintain $\partial _r u |_{r=0} =0$ . Otherwise, the flow would have a cusp (discontinuous derivative of $u$ ) at $r=0$ ). The same steps followed in the semi-infinite space problem of § 3 lead to a long expression for $u(r)$ which is a plug-like flow with very thin boundary layers of thickness of ${\sim} \kappa ^{-1}$ located at the cylinder wall $r=a$ . The meaningful observable here is the average value of the zero mode over the capillary cross-sectional area, defined as

(5.9) \begin{equation} \langle u \rangle = \frac {1}{\pi a^2} \int _{0}^a {u(r)} 2\pi r {\textrm{d}}r, \end{equation}

where $u(r)$ is given by (5.8). Equation (5.9) is a double integral and involves integrals of Bessel function pairs which, in general, cannot be calculated in closed form (for exceptions, see for instance Luke (Reference Luke1962)). We circumvent this complication by exchanging the order of integration in the double integral appearing in (5.9). Thus, (5.9) can be written as a single integral

(5.10) \begin{equation} \langle u \rangle = \frac {-1}{a^2} \frac {i\epsilon k}{4\eta } \int _{0}^a \left ( \phi ^*\phi _r - \phi \phi _r^*\right ) r^2 {\textrm{d}}r, \end{equation}

and is calculated numerically by simple quadrature. In figure 8, we display the zero mode velocity averaged over the channel width (5.10) versus the frequency of charge oscillation $\omega$ scaled by the frequency $D\kappa ^2$ , where $D$ is the diffusion coefficient of the charge distribution.

Figure 8. Average velocity in a cylinder, cf. figure 7, with travelling wave charge distribution (that is, Neumann boundary conditions for the Poisson equation). Self-similar behaviour of the zero mode velocity (5.9), averaged over the cross-sectional area of the cylinder of radius $a$ as a function of the dimensionless frequency $\omega /(D\kappa ^2)$ of the travelling wave wall charge. The curves are self-similar as described by relation (5.11). Continuous line denotes the fit (5.12) by adopting the single-parameter value $\beta = 1.9$ . $a = 10^{-5}$ m, $D = 10^{-9}\textrm { m}^{2}$ s⁻¹.

5.1. Self-similar behaviour of the average velocity in the cylindrical capillary

The zero mode velocity (5.10) displays a self-similar behaviour for the same reason discussed in the channel case. From the available dimensionless combinations (Bluman & Kumei Reference Bluman and Kumei1989; Panton Reference Panton1996), we display the following in the average velocity (5.10):

(5.11) \begin{equation} \langle u \rangle = \frac {1}{2}\frac {\sigma _0^2}{\eta \epsilon k} \mathcal {F}\!\left (\frac {\omega }{D\kappa ^2}, \kappa a\right ), \end{equation}

where $\mathcal {F}$ is a nonlinear function of the two independent dimensionless parameters and the front factor is the velocity scale $\kappa u_0/k$ cf. (2.12), as this is displayed in the label of the vertical axis of figure 8. A good single-parameter fit of these results that accounts for all parameter values, including the tails in figure 8, has the functional form

(5.12) \begin{equation} \langle u \rangle _{E\&M} = \frac {\beta }{ \kappa a} \frac {\frac {\omega }{D\kappa ^2} }{\left [\left (\frac {\omega }{D\kappa ^2}\right )^{2}+1\right ]^{\frac {5}{4}}}\cos \! \left (\frac {\arctan \left (\frac {\omega }{D\kappa ^2}\right )}{2}\right ), \end{equation}

i.e. it is the velocity (C2) scaled accordingly and we have adopted the single-parameter value $\beta = 1.9$ in figure 8. In figure 8, we plot the function $\mathcal {F}$ from (5.11) versus the dimensionless (Debye) frequency for a number of values of the excitation wavenumber  $k$ .

We estimate the cylinder average velocity by considering the same material parameters as in § 3.2. Adopting $E = 10^5$ V m⁻¹, as done by Ehrlich & Melcher (Reference Ehrlich and Melcher1982), $\omega \sim D\kappa ^2$ , thus giving the peak dimensionless velocity value of ${\sim} 0.0715$ in figure 8, leading its dimensional counterpart to equal

(5.13) \begin{equation} \langle u \rangle = 2.5 \textrm { cm}\, \mathrm{s}^{-1}. \end{equation}

Notice, however, that this estimate decreases with increasing $k$ .

6. Numerical results

In this section, we solve the Poisson–Nernst–Planck–Navier–Stokes differential equations numerically with a finite-element package (Comsol). We are interested in establishing the existence of the zero mode velocity and how does it compares in magnitude and its general trend with the exact solutions obtained in the previous section.

We employ the cylindrical capillary geometry of § 5 amended so as to incorporate effects of boundaries in the longitudinal direction of the cylinder as displayed in figure 9 and explained in more detail in the supplementary materials addendum. We employ the low-Péclet-number approximation to the Poisson–Nernst–Planck–Navier–Stokes differential equations which, for the material parameters employed in this paper, is a valid approximation as shown in Appendices B.1 and B.2.

Numerically, we solve the full Poisson–Nernst–Planck–Navier–Stokes system and obtain the species distributions, electric field and velocity, without other approximations. For instance, the time-dependent Stokes equations read

(6.1) \begin{equation} \rho _l\frac {\partial \mathbf {v}}{\partial t} = -\nabla p + \eta \nabla ^2 \mathbf {v} + \rho \nabla \phi , \quad \nabla \cdot \mathbf {v} = 0. \end{equation}

Figure 9. Configuration of the domain for the numerical solution of the Poisson–Nernst–Planck–Navier–Stokes system in § 6.

The set-up of the system is shown in figure 9. The cylinder has length $L$ and is connected to two reservoirs of fixed charge concentration and zero electric potential (left-most and right-most edges of the figure). The boundary conditions in the lateral surface of the cylinder are the same as those employed in § 5 of the manuscript and we employ the real part of the wall charge density (thus, $\sigma (z,t) = \sigma _0 \cos (kz - \omega t)$ ). The intermediate regions connecting the two reservoirs to the cylinder have boundary conditions of zero normal electric field and zero normal current $\hat {\mathbf {n}}\cdot \mathbf {J}_\pm$ , where

(6.2) \begin{equation} \mathbf {J}_\pm = -D\left [ \nabla c_\pm \pm \frac {e}{k_BT} c_\pm \nabla \phi \right ] \end{equation}

and $\hat {\mathbf {n}}$ is the normal vector to the interface. We employ the parameters

(6.3) \begin{align} a & = 10^{-5} \textrm { m}, \quad L =2\times 10^{-4} \textrm { m}\ \textrm { or }\ L = 4\times 10^{-4} \textrm { m}, \quad D = 10^{-9} \textrm { m}^2 \ \textrm {s}^{-1}, \notag\\ \sigma _0 & = 7\times 10^{-5} \textrm { C m}^{-2},\quad \kappa = 2.3 \times 10^6 \textrm { m}^{-1},\quad k = 0.8 \times 10^5 \textrm { m}^{-1},\notag\\ 2c_\infty & = 10^{-3} \textrm { mM}. \end{align}

Figure 10. (a) Comparison of the exact zero mode velocity (5.9) (continuous curve, calculated with a cylinder of infinite length) with two of its numerical counterparts (circles and triangles) obtained with the finite-element package Comsol in the geometry displayed in figure 9 by employing the parameter values (6.3) (with cylinder aspect ratio $L/a=20$ and $L/a=40$ ). In general, the agreement between the exact and the numerical results is very good near the tails. The difference between numerical and exact solutions is attributed to hydrodynamic entrance and edge effects when the aspect ratio $L/a$ is relatively small. As the aspect ratio $L/a$ increases, however, the numerical solution converges to the exact as seen by comparing the circle with the triangle data points. (b) Following figure 10(b) of Cahill et al. (Reference Cahill, Heyderman, Gobrecht and Stemmer2004), we display the same data points/continuous curve as in panel (a) but now dividing the velocities by their maximum value. In general, the agreement between exact and numerical results is excellent.

To obtain the zero mode, we average the time-dependent velocity field over the cylinder cross-section, the cylinder length $L$ (cf. figure 9) and over $N$ periods of oscillation of the applied field, where $N$ is an integer. The details are described in the supplementary materials addendum.

The circle and triangle symbols in figure 10(a) denote the velocity obtained numerically with the configuration of figure 9 and parameters (6.3). In general, there is excellent agreement with the exact solution (5.9) (continuous curve) at the tails. Away from the latter, the numerical and exact results agree well in an order-of-magnitude basis and share the same general trend.

The difference seen between the averaged zero-mode velocity obtained numerically here and the exact model, as seen in figure 10(a), can mostly be attributed to the requirement of finite length geometry for the numerical configuration. Indeed, we have verified that by increasing $L$ , the numerical results approach their analytical counterparts as can be seen by comparison of the circle and triangle data points in figure 10(a) with its exact counterparts. One explanation for this behaviour is that small length $L$ of the capillary tube is mainly associated with a hydrodynamic entrance flow and edge effects, while the analytical result assumes that the flow is maintained in a domain of infinite length. Figure 10 is generated by employing the parameter values (6.3), where the tube aspect ratio (length-to-radius) is $L/a = 20$ and $L/a=40$ . In the supplementary materials addendum, we provide the corresponding data tables. A full parametric analysis of travelling wave electroosmosis in a finite-length capillary tube is beyond the scope of this paper.

Finally, following figure 10(b) of Cahill et al. (Reference Cahill, Heyderman, Gobrecht and Stemmer2004), we compare the exact and numerically obtained zero-mode velocities by scaling each one of the curves in figure 10(b) by their respective maximum. Now, there is an excellent agreement between the exact and the numerically obtained zero mode velocities as displayed infigure 10(b).

In the supplementary materials addendum, we briefly discuss the effect a finite Péclet number has on the velocity magnitude.