Microfluidic pumping by micromolar salt concentrations

Abstract

An ion-exchange-resin-based microfluidic pump is introduced that utilizes trace amounts of ions to generate fluid flows. We show experimentally that our pump operates in almost deionized water for periods exceeding 24 h and induces fluid flows of μm s⁻¹ over hundreds of μm. This flow displays a far-field, power-law decay which is characteristic of two-dimensional (2D) flow when the system is strongly confined and of three-dimensional (3D) flow when it is not. Using theory and numerical calculations we demonstrate that our observations are consistent with electroosmotic pumping driven by μmol L⁻¹ ion concentrations in the sample cell that serve as ‘fuel’ to the pump. Our study thus reveals that trace amounts of charge carriers can produce surprisingly strong fluid flows; an insight that should benefit the design of a new class of microfluidic pumps that operate at very low fuel concentrations.

---

Fluid, solute, and colloid transport on the microscale pose a significant challenge, due to external pressure-driven pumping requiring the pump itself to withstand large forces. To circumvent this issue, a range of microfluidic pumps has recently been developed,^1–21 most of which exploit self-generated solute gradients. Typically, a solute-gradient-based (osmotic) pump consists of a source/sink of solute molecules, close to the surface of a sample cell. The solutes are produced/consumed either by chemical reactions on the surface of the pump,^{3–5,9,11–15,17–24} the pump slowly dissolving,^8,22,24 or exchange reactions taking place within the pump.^16,23 This sets up a concentration gradient in the fluid and along the surface of the sample cell. The interaction between the solutes and the surface causes a force on the fluid, which—coupled with the spatial heterogeneity of the solutes—leads to fluid flow, in a process that is referred to as osmosis.²⁵

In osmotic pumps, the solute thus acts as ‘fuel’, which enables the pump to move fluid around. Such pumps generate relatively small forces applied over a much larger range of the fluid through long-ranged concentration gradients, thus overcoming the issues that face external pressure-driven pumps. Depending on the nature of the surface–solute interactions, neutral or Coulombic, the pump is either diffusioosmotic or electroosmotic. However, there is strong evidence that solute-density²⁰ and thermal-convection¹⁸ effects can also play a role for large pumps.

A fundamental problem for microfluidic pumping based on osmosis is the need for a solute (gradient) in the fluid medium that also contains the material to be transported, since solutes can interact with the transported material. For instance, pumps that utilize catalytic decomposition of hydrogen peroxide^3,4,13 or hydrazine⁵ are not biocompatible and these solutes will also react (unfavorably) with other materials. The working of pumps that instead employ enzymatic reactions to convert biomolecules, such as urease,^18,20 will be inhibited when material that is transported reacts with the solute itself. This also limits their use in transporting biological material, which typically interacts with such biomolecules. To solve this issue, it is desirable to design pumps that are driven by unreactive solutes, preferably in minimal amounts, and that have a limited impact on their environment.

In this paper, we introduce a microfluidic pump that accomplishes this goal. Our pump is experimentally shown to function in almost completely deionized water for periods of over 24 h. We study the fluid flow by means of tracer velocimetry (close to the bottom of the sample cell) and show that the pumping speed is in the μm s⁻¹ range over hundreds of μm. The dependence of this flow on the size of the pump and the added salt concentration in the system is also characterized. It is further experimentally demonstrated that solute-density and thermal-convection effects do not play a role in our system. We therefore hypothesize that our pump operates on trace amounts of ions present in the bulk fluid, by exchanging one species of ion for another, thereby generating a diffusion potential which drives electroosmotic flow. This sets it apart from other microfluidic pumps that generate flow by slow dissolution of the pump itself, see, e.g., ref. 8 and 22. Specifically, our pump only modifies the identity of the ionic species in the bulk, whereas dissolving pumps increase the bulk ion concentration.

Furthermore, we show experimentally that the decay of the flow velocity can be modified by changing the geometry of the sample cell on the mm length scale. The far-field, power-law decay of the speed with the radial distance r is either quasi-2D (∝r⁻¹) for small cell heights (≤2 mm), or 3D (∝r⁻²) for tall cells (≥10 mm). Even in the quasi-2D regime, our system displays almost time-independent (steady-state) fluid flow profiles. This is surprising, as 2D diffusive systems are not expected to exhibit steady-state solutions.

We interpret our experimental findings using a combination of the numerical finite-element method (FEM) and analytic calculations. It is shown that the experimental observations can indeed be understood by the resin exchanging trace amounts of cations from its surroundings with protons from its interior. We estimate the relevant trace cation concentration to be in the low micromolar range. The experimental observations are further shown to be consistent with an electroosmotic pumping mechanism: the difference in ion mobility between the protons and the exchanged cations sets up a diffusion potential that causes flow toward the exchange resin in the absence of a net electrical current. The mechanism is the same as previously found for similar ion-exchange pumps¹⁶ as well as dissolving pumps.^8,22 However, our results indicate that ion-exchange-resin-based microfluidic pumps have a surprisingly small lower bound to the ion concentration under which they can operate, which we chart in this paper.

In our numerical work, we directly model the electroosmotic flow generated by ion exchange in the geometry of the experiment. We employ steady-state solutions for the concentration fields, electrostatic potential, and fluid velocity using the FEM. These computations go far beyond the thin electrostatic screening limit that is typically considered for such systems and give insight into the flow throughout the cell. Using analytic theory, we investigate the time dependence of the flow in the quasi-2D far field. We use our analytic theory to prove that the experimentally observed steady-state flow can be explained by the fact that the flow is driven by concentration gradients. While the relevant solute concentrations evolve over time and have no steady state, the concentration gradients become time-independent beyond a characteristic, system-dependent diffusion time that we identify.

Finally, we can explain the scaling of our results with cell height in terms of interaction between the out-of-equilibrium ion fluxes and the confining geometry. Here, we observe qualitative, but not quantitative, agreement between the experiments and the numerical calculations. In the experiment, the power-law decay of the flow sets in unexpectedly close to the ion-exchange resin. We argue that this is due to the neglect of solute transport by advection in our calculations, which is necessary to make progress in both numerical and analytic theory. Accurately modeling the near-field effect of advection will be important to understanding the formation and performance of swimmers comprised of mobile ion-exchange resins and inert particles²⁶ and therefore presents challenges for future study.

Our results on ion-exchange-resin-based microfluidic pumps lead to the startling finding that trace amounts of ions are sufficient to generate significant fluid flow, which is driven by diffusion-potential electroosmosis. This insight should prove instrumental for the design of new microfluidic pumps operating in close-to-deionized water, which is the natural and often desirable environment in which to perform experiments. It furthermore provides compelling evidence that the effect of small amounts of charge and minute ionic fluxes may have significant consequences in other systems, such as chemically self-propelled colloids.

1 Experiments

In this section, we describe the experimental setup for a single ion-exchange-resin pump and characterization of the tracer properties used in our velocimetry measurements. We also provide quantification of a wide range of resin pumps and tracers to show the generality of our findings. Finally, we study the impact of added salt on the pumping.

1.1 Tracer characterization

Polystyrene (PS) tracers were used for the velocimetry (PIV) measurements of our ion-exchange-resin pump. Stock PS particle suspensions (Microparticle GmbH, Germany) were diluted with distilled water and thoroughly deionized using ion-exchange resin (Amberlite K306, Carl Roth GmbH + Co. KG, Karlsruhe, Germany). The electrophoretic mobilities of the PS tracers were determined using micro-electrophoresis in a custom-built, disposable setup. For this setup, we employed a Perspex cell 45 mm in height and with a square cross section (10 mm edge length). Based on the geometry of Uzgiris^27,28 two platinum electrodes of width 1 mm were mounted into the center of the Teflon® cap sealing the cell. This ensures sufficient electrode-wall distances to effectively reduce stray-field-driven electroosmosis at the cell walls. The electrode spacing was set to 1 mm to obtain homogeneous electric fields and square-wave alternating voltages of ±1 V were applied by a function generator (PeakTeck 4060 by PeakTeck GmbH, Germany). The cell was mounted on the stage of a micro-electrophoresis instrument (Mark II, Rank Bros. Bottisham, Cambridge, UK) supplying ultramicroscopic illumination and particle tracks were imaged using exposure times of 3 s on a consumer digital single-lens reflex camera (DSLR; D800, Nikon, Japan).

Electrophoresis of the tracers was performed in the horizontal direction, while the particles sedimented in the vertical direction due to gravity. Thus the trace of a single particle has a saw-tooth shape. The mobility of the particles μ_E was calculated from the averaged velocity in the horizontal direction v_E and is given by μ_E = v_E/E, with E the amplitude of the electric field. The obtained values of μ_E are listed in Table 1. The values of μ_E are relatively low due to the non-monotonic scaling of the electrophoretic mobility at low salt concentrations.²⁹

Table 1 The electrophoretic mobility μ_E of the PS tracer particles used in this work. The first column gives the label, the second the diameter, and the third the value of μ_E

Type (—)	Diameter (μm)	μ E (10^−8 m^2 V^−1 s^−1)
PS1	1.7 ± 0.1	−2.0 ± 0.2
PS7	7.6 ± 1.0	−2.6 ± 0.3
PS10	10.4 ± 0.9	−2.5 ± 0.3
PS15	15.2 ± 0.9	−2.5 ± 0.2
PS15COOH	15.5 ± 0.2	−2.1 ± 0.2

---

1.2 Velocimetry for the ion-exchange-resin pump

For the characterization of the ion-exchange-resin pumps via tracer velocimetry, we constructed custom sample cells of radius R = 10 mm and several heights H out of poly(methyl methacrylate) (PMMA) rings attached to a microscopy glass slide and covered with another glass slide (soda lime glass of hydrolytic class 3 by VWR International), see the sketch in Fig. 1a. The glass slides were washed with alkaline solution (Hellmanex® III, Hellma Analytics) by sonication for 30 min, then rinsed with tap water, and finally washed several times with doubly distilled water (distilled using a Quartz Hareaus Destamat; the conductivity was measured to be 55 nS cm⁻¹). Spherical cationic resin beads (CGC50×8, Purolite Ltd, UK; exchange capacity 1.7 eq. L⁻¹§) with radii ranging from 10 to 50 μm were carefully glued to the bottom glass slide with a tiny amount of two-component glue (UHU plus sofortfest, UHU GmbH, Germany), which was then set aside for 24 h to allow the glue to completely solidify. One resin bead was glued in each sample cell.

Fig. 1 The ion-exchange resin and sample cell. (a) Sketches of the geometry, showing top and side view of the cylindrical sample cell, with radius R and height H. A zoom-in shows the exchange resin, with radius r_R in red, a polystyrene (PS) tracer in green flowing along a blue flow line. In this paper, radial distance r is measured from the center of the resin. (b) Experimentally measured trajectories for the PS tracers (PS7) toward the resin (center of image), as shown in a top-view microscopy image. Blue arrows indicate the paths of the tracers, which are obtained from our image analysis.

The sample cell for the ion-exchange-resin pump experiments was loaded with a dilute PS-tracer suspension, prepared according to the above deionization procedure. It was subsequently mounted on the stage of an inverted scientific microscope (DMIRBE, Leica, Germany), and observed in bright field, typically at 5× magnification. Images were shot with a DSLR and videos recorded with standard video equipment at frame sizes of 5.2 Mpix and frame rates of 30 fps. We imaged an area with cross-section of (typically) larger than 1000 μm, slightly above the bottom glass plate, focusing on the average hovering height of gravitationally settled PS tracers, see Fig. 2.

Fig. 2 Typical microscopy image of the tracers (red outlined circles) around the ion-exchange resin (black circle). The red outlines are obtained using the fit algorithm explained in the text.

The resin bead glued to the glass slide displayed significant fluid pumping, as evidenced by the PS tracers moving toward the individual resin beads; Supplemental Movie “Exchange_Resin_Pump.avi” gives an example of this for tracers close to the resin (ESI†). These tracers come in from far away along the substrate, move up vertically from the substrate close to the resin, then move radially away from the resin, subsequently sediment to the substrate away from the resin, and finally move back toward the resin along the substrate, leading to a recirculation of the tracer particles. Along their path the tracer speed varies as a function of r. The radial dependence of U_PS was determined from the tracer positions in successive frames of the recorded movies. These positions were extracted using an in-house Python code. In brief: the circular perimeter of each particle was extracted using standard edge-detection methods, and then fitted to a circle using the Hough transform,³⁰ implemented in the OpenCV function HoughCircle, see Fig. 2. Tracer positions in consecutive frames were compared to determine radial velocity. The velocity of a given tracer species for a specific ion-exchange-resin bead size was measured for 80–100 individual PS particles for each bead and the results averaged over some 40–50 beads.

We made use of the following expression to determine the tracer speed U_PS

where r̂ is the 2D unit vector pointing from the resin to the tracer, Δs⃑ is the displacement of the tracer between frames (time between frames Δt), “·” is the inner product, and 〈⋯〉 indicates averaging over all tracers that are a distance r from the resin's center.

The results of our velocimetry are shown in Fig. 3, which provides U_PS as a function of the radial distance. Two regimes can be distinguished. For r ≲ 75 μm, there is a slight increase in the tracer speed, followed by a maximum and subsequent decrease (this is more evident in Fig. 4). For r ≳ 75 μm the speed decreases with a power law and is appreciable over at least 300 μm. For sample cells with a height of H = 1 mm, we find that U_PS ∝ r^−0.9±0.1 in the far field (H = 0.5 mm, U_PS ∝ r^−1.2±0.1; H = 2 mm, U_PS ∝ r^−1.1±0.1), while for the sample with height H = 10 mm, the fitted decay is U_PS ∝ r^−2.2±0.3.

Fig. 3 The speed of the tracer U_PS as a function of the distance r for several values of the sample cell height H, an ion-exchange resin with radius r_R = 22.5 μm, and PS7 tracers. The symbols show the experimentally measured values; the standard error is given for each data point. The gray dashed lines serve as a guides to the eye for the power-law decay.

Fig. 4 Velocity of tracer particles U_PS as a function of radial distance r for a cell height of H = 1.0 mm. (a) Three different sized PS particles are pumped by resin beads with r_R = 22.5 μm. (b) Three different sized resin beads are used to form a pump, the tracer particle is PS7 for each. (c) Two PS tracers of the same size, but with different μ_E, are pumped by resin beads with r_R = 22.5 μm. In all panels, the gray dashed line serves as a guide to the eye for the power-law decay and the standard error is given for each data point.

We concentrate on understanding the far-field regime throughout this paper, as in the near field there are several competing effects, including electrophoresis, local flow, and interaction with the substrate, which complicate understanding of the physics. For instance, it is difficult to assess on the strength of our experiments and the theory what the cause of the apparent near-field maximum in the tracer speed is. Fluid incompressibility could explain the decrease in speed close to the resin, i.e., an increasingly upward-directed component of the near-field flow requires a decrease in the horizontal component. However, other possibilities cannot be excluded at this time.

1.3 Resin size, tracer properties, and salt concentration

In this section we demonstrate the generality of the fluid pumping by ion-exchange resins. The systematic quantification of the tracer speed U_PS as a function of r is shown in Fig. 4.

• In Fig. 4a we vary the size of the PS tracers for tracers which have similar electrophoretic mobility. In the far-field region there is a power-law decay of the tracer speed, which is insensitive to the type of tracer used within the error bar. This shows that in the far-field the size of the tracer does not play a role.

• Fig. 4b shows results for three resin sizes (radius r_R). A larger resin induces a stronger electroosmotic flow over a larger range. We analyzed the far-field tracer speed by fitting the three curves using power-law decays. Then, we established the speed at an arbitrary far-field distance (r = 150 μm) as a function of the size. For these three data points, we found a linear dependence through the origin U_PS(r = 150 μm) ≈ r_R × (5.5 ± 0.5) × 10⁻² s⁻¹; the offset ≈0.15 μm s⁻¹ is negligible—similar scaling was observed for other far-field distances. This strongly indicates that the process is diffusion limited. By diffusion limited, we mean the upper speed limit imposed by the rate at which ions can diffuse towards the resin bead from the bulk reservoir. In this limit the flux through the particle surface j_s,dl (per unit area) is determined by the diffusivity of the ions D* and the concentration far away ρ*, with the familiar diffusion-limit scaling j_s,dl ∝ D*ρ*/r_R (ref. 31). The speed is proportional to the total flux through the resin, i.e., U_PS ∝ 4πr_R²j_s,dl ∝ D*ρ*r_R, giving the linear dependence with r_R observed in the experiment.

• In Fig. 4c, we vary the electrophoretic mobility of the tracers, but not their size. It is evident that these tracer particles have the same velocity within the error bar in the power-law regime. This shows that the results are reproducible with nominally similar (μ_E is comparable within the error bar), but possibly slightly different particles.

Finally, we added KCl solution (Merck KGaA, Germany) to the sample cell for H = 1 mm and the r_R = 22.5 μm resin beads. Fig. 5 shows the change in tracer speed: adding 5 μmol L⁻¹ KCl increases U_PS, adding 10 μmol L⁻¹ instead, increases the speed further. That is, a higher concentration of exchangeable ions induces stronger flow. However, at a KCl concentration of 80 μmol L⁻¹, the velocity of tracer particles is effectively zero (therefore not shown here). For the 80 μmol L⁻¹ sample, we also do not observe any Brownian motion of the tracer beads. This indicates that the beads have become firmly stuck to the sample cell wall, probably because of the increased electrostatic screening at this higher salt concentration. Therefore, we cannot use this 80 μmol L⁻¹ data to infer a drop in pumping speed at higher salt concentration. This is in line with similar findings for chemically-propelled swimmers in ref. 32.

Fig. 5 Tracer speed U_PS as a function of radial distance r for different added KCl concentrations: no added salt (green plusses), 5 μmol L⁻¹ (purple crosses), and 10 μmol L⁻¹ (yellow triangles), for a cell with H = 1.0 mm and a resin with r_R = 22.5 μm. The gray dashed line serves as a guide to the eye for the power-law decay and the standard error is given for each data point.

1.4 An inverted pump

We inverted our setup to check whether solute density variations or thermal convection effects played a role in our system, as is the case in ref. 18 and 20. That is, we glued the resin to the top glass slide and examined the movement of the tracers. In order to ensure that the tracers were at the top cover slide, we modified the overall density of the solution by adding glycerol (water : glycerol mass ratio of 1 : 0.3) to slightly exceed the density of our PS particles. We used tracers with a diameter of 3.3 μm here and we increased the size of the ion-exchange resin to r_R = 250 μm, in order to increase the speed of the tracers in this mixture of higher viscosity (approximately double that of water).

Supplemental movie “Inverted_Resin_Pump.avi” shows the result of this experiment (ESI†). It is clear that inverting the pump did not change the direction in which the tracers move toward the resin. While we increased the overall density of the mixture, this should not affect the possible density variations induced by ion exchange. Our experiment therefore rules out density variation effects.

2 General considerations

We performed a theoretical/numerical analysis of the pump to gain understanding of the fluid flow observed in our experiments and to show that the observed fluid pumping is indeed caused by trace amounts of cations in the sample cell. In this section, we start with several general considerations to provide a background for our calculations.

• The resin is designed to exchange H⁺ for other cations in the bulk, with a capacity of 1.7 eq. L⁻¹. There are no decomposition-type surface reactions. Nor does our pump itself dissolve, as is the case in ref. 8 and 22.

• The surface of the sample cell is charged and we measured the zeta potential of the bottom glass slide to be ζ ≈ −(105 ± 5) mV.¶ The negative surface charge is due to the dissociation of surface groups, which release cations into the bulk. However, because the sample cell is filled with deionized water and carefully rinsed before preparation, there will be very few non-protonic cations present.

• Dissolved CO₂ forms carbonic acid and thus creates cations in the form of protons that screen the wall charge. We measured a pH ∼ 5.4 at the onset of the experiment, consistent with typical pH values for water in equilibrium with atmospheric CO₂.³⁵ The cations associated with CO₂ dissociation (protons) are the same as the ions inside the ion-exchange resin, so they cannot contribute to electroosmosis via ion exchange.

One might assume that the exchange resin cannot exchange protons for other cations, as any non-protonic cations released from the cell walls will have been washed away during preparation of the sample cell, according to our second point. This would then prevent ion-exchange-based electroosmosis.

However, we will argue that a very low concentration of cations remaining in the bulk after cleaning would be sufficient to fuel pumping. We estimate here the residual concentration of ions that would be required in this case. During a period of 24 h,|| a spherical pump of r_R = 22.5 μm in radius exchanges at most a part of its ionic content via an (assumed) constant surface flux density j_s. Let us further assume 25% of the original content to account for an unmodified pumping speed over the course of the experiment. Then it follows from the resin size and exchange capacity that j_s ≲ 5 × 10⁻⁸ mol m⁻² s⁻¹. This value is reasonable, as similar numbers are found for self-electrophoretic Janus swimmers that move at speeds comparable to our maximum U_PS (ref. 31). To make this level of exchange possible, the non-protonic cation concentration in the sample cell has to be at least 1 × 10⁻⁷ mol L⁻¹—the total number of exchanged ions over 24 h and the volume of the sample cell were used to arrive at this number. Again, we can assume only a fraction of the total ions present are exchanged. This would lead to an estimate for the cation concentration of ρ* = 1 μmol L⁻¹ (if 10% is exchanged). Such low ion concentrations could be contributed by “impurity” cations released from the glass slides into the bulk fluid following rinsing. Taking ρ* and a typical value for cation diffusivity of D* ≈ 2 × 10⁻⁹ m² s⁻¹, the surface flux in the diffusion-limited regime is j_s,dl ≈ D*ρ*/r_R = 4 × 10⁻⁸ mol m⁻² s⁻¹ (using K⁺ for the cationic contaminant, which has a typical cation diffusivity). Our estimate for j_s is thus in the physically reasonable regime, close to the diffusion limit, in accordance with our experimental result.

This proposed mechanism of generating fluid flow by ion exchange is as follows. Exchanged protons moving away from the resin have a higher mobility than the to-be-exchanged cations moving toward the resin—H⁺ has the highest mobility of any ion—and thus the protons have a greater diffusive flux. To prevent bulk charge separation due to the difference in diffusivity, an electric field (E-field) is formed to compensate for the difference in diffusive flux with a migrative flux (via the E-field). The E-field points toward the resin and prevents charge separation in the bulk, by slowing down the H⁺ and accelerating the cations, such that the total fluxes remain equal and opposite throughout. Since the associated electrostatic potential stems from a difference in ion diffusivity (equivalently mobility) this mechanism is “diffusion potential” based.^16,25 The E-field acts on all ions in the system. However, the E-field only exerts a significant force on the fluid in a screening layer close to the chamber boundaries, where there is an excess of cations. In the electrically neutral bulk, the E-field has a vanishingly small effect. The result is that the E-field drives fluid flow along the glass slides towards the resin, which then, through incompressibility, generates a backflow outwards along the horizontal center-plane of the chamber.

We verified that this simple picture and our estimates can indeed give rise to the observed flow speed and direction using finite-element-method (FEM) simulations representative of the experimental geometry, see Sections 4 and 5. We also used linearized, analytic theory to study the regime where the geometry can be considered as quasi-2D, i.e., for r > H, see Section 6. In both cases, we solved the associated time-dependent electrokinetic equations, which we discuss next.

3 Electrokinetic equations

To model the electroosmotic flow around the ion-exchange resin, we require three coupled equations, collectively known as the electrokinetic equations: Nernst–Planck for the solutes, Poisson for the electrostatics, and Stokes for the fluid flow, together with boundary conditions for the respective problems. We explain the three equations in detail below and discuss the boundary conditions for the simulations in Section 4 and the approximations made for the theory in Section 6, respectively.

The Nernst–Planck equation describes the diffusion and migration of the solute species. Here, we consider three ionic solute species in the fluid, protons H⁺, potassium K⁺, and chloride Cl⁻. The protons are loaded in the exchange resin and the choice of the two other ions is arbitrary. Here, we selected two ions with almost equal diffusivities D_K⁺ ≈ D_Cl⁻ (ref. 36) to avoid the complication of additional (but relatively small) diffusion potentials. These three species are indexed by i ∈ {H⁺,K⁺,Cl⁻}, respectively. We write ρ_i for the time- and space-dependent concentration fields and D_i for the molecular diffusivities. Then the flux of each species is given by

where u⃑ is the fluid velocity (accounting for advection), k_B is Boltzmann's constant, T the temperature, e the elementary charge, z_i the valency, Φ the electrostatic potential, and the gradient operator. The continuity equation is given bywhere ∂_t denotes the time derivative. For the steady-state problem ∂_tρ_i = 0.

Before we move on to the other equations, we should comment on two simplifying assumptions typically made in the above description.

• We ignored the advective contribution to the flow in eqn (2) in all our calculations. We consider the Péclet number, which give the ratio of diffusion to advection, to examine whether this is reasonable. A simple estimate is as follows: using a typical length scale of H for the development of fluid flow in this problem, a typical non-protonic ion diffusivity of D ≈ 2.0 × 10⁻⁹ m² s⁻¹ (ref. 36), and a typical speed of Ū_PS ≈ 1.0 μm s⁻¹, we arrive at Pe = Ū_PSH/D ≈ 5. This indicates that the value of Pe is probably high,** so the advective term should not be ignored in eqn (2). However, due to the computational complexity of our FEM calculations,†† as well as the need to linearize our analytic theory, this approximation must be made in order to make progress. As we will see, the understanding of the physics of the resin pump is not strongly affected by this reduction.

• We have ignored bulk ionic association–dissociation reactions, as described in ref. 31, which would have entered on the right-hand side of eqn (3) as coupled chemical source and sink terms. In the physical system, bulk exchange will lead to coupling of the H⁺ flux coming from the ion-exchange resin and H₂O and OH⁻ present in solution via H₂O ⇌ H⁺ + OH⁻.‡‡ The main effect of these bulk reactions would be to replace the relevant diffusion rate D_H⁺ with an effective rate

with ρ^∞_i the concentration “very far” away from the resin.³¹ For the experimental pH ≲ 5.4, [OH⁻] ≪ [H⁺], so D_av ≈ D_H⁺, to within 2%. Hence, it is justified to ignore the bulk reactions here.

The electrostatic potential fulfills the Poisson equation, which is given by

where ε₀ is the vacuum permittivity and ε_r the (spatially constant) relative permittivity. It should be noted that the ρ_i and Φ in the Poisson equation are time dependent and these two quantities provide the coupling between the Nernst–Planck and Poisson equations. For completeness, we introduce the electrostatic screening (Debye) length κ⁻¹ here via

Finally, we have the incompressible Stokes equations to describe the fluid flow. These read

with η the viscosity of the fluid, p the hydrostatic pressure, and f⃑ the body-force density. Here, u⃑ and f⃑ are time-dependent quantities. Given the electrostatic potential and the densities of all (charged) species, we can specify the body-force density on the fluid to close the problem asThis expression was obtained by first-order expansion of the chemical potential around thermodynamic equilibrium, which gives the gradient of the chemical potential as a driving force.³⁸ The specific choice of this driving force is to eliminate the spurious flow due to inexact cancellation of pressure and electrostatic interactions in FEM calculations. It is, however, completely equivalent³⁸ to the more commonly used expressionThe only difference between eqn (9) and (10) is the interpretation of the hydrostatic pressure: eqn (9) does not, while eqn (10) does include the ideal-gas contribution from the dissolved solutes.³⁸

4 Finite-element model of the pump

In this section we describe the boundary conditions for the above equation system and the specific choices made for the FEM modeling. Throughout, we used COMSOL Multiphysics® Solver 5.2a to numerically solve the electrokinetic equations for a model setup of the experimental geometry.

We considered a 3D cylindrical portion of the microscopy cell, with the resin located on the symmetry axis of the cylinder. Due to the rotational symmetry of our setup, which corresponds closely to the experiment, the simulations could be performed on a quasi-2D axisymmetric domain, see Fig. 6. We considered two domains in order to simulate small sample heights H ≲ 1 mm (a) and large sample heights H ≳ 2 mm (b) for the steady-state electrokinetic equations; we will come back to this in Section 5. The latter domain is a half-open domain, which we will refer to as “unbound”.

Fig. 6 Schematic of the rotationally symmetric (a) top-bound and (b) “unbound” simulation domain, the red dashed line shows the symmetry axis. The resin is modeled as a hemisphere of radius r_R (lower-left corner). Cation exchange is modeled by an inward/outward directed flux j_s of cations (K⁺) and protons (H⁺), respectively, see eqn (12). A constant surface charge density σ_wall is imposed on the bottom (and top wall) and on the resin. All solid surfaces form no-slip boundaries for the hydrodynamics. The right-most boundary (orange line (a) or a circular arc (b)) is an “open boundary” for the hydrodynamic problem and we impose a pre-computed electrostatic profile and ion distributions on it, as explained in the text. Cut lines are used to emphasize that the domain is much larger than the resin, see Fig. 7.

Let us first describe the simulation domain that most accurately represents the experiment, see Fig. 6a. The bottom and top of the simulation domain correspond to the glass slides of the sample cell, the height of the sample cell H is fully resolved. The radius of the simulated geometry is R_sim. The spherical resin (experiment) is modeled as a hemisphere of radius r_R attached to the lower boundary (substrate). We chose a hemi-spherical resin, rather than a fully spherical one—as in the experiment—for simulation convenience. Specifically, the choice of a hemispherical resin facilitates the use of quadrilateral elements for the mesh, see inset to Fig. 7. This meshing would not be possible for a resin sphere in contact with the substrate, as is likely the case in the experiment, due to the cusp-like feature present in that geometry.

Fig. 7 Example of the fine mesh used for our FEM calculations (H = 1 mm). Close to the electrically charged surfaces (within 6 Debye lengths) quadrilateral elements are used, see insets. The rest of the domain is composed of triangular mesh elements.

Quadrilateral elements are necessary, since we use the spurious-flow reducing method of ref. 39–41. This method consists of finely meshing several Debye lengths (in our case 6) using such elements around the walls and the rest of the domain using triangular elements, see Fig. 7. Quadrilaterals have the advantage that larger aspect ratios are permitted than for triangles, before numerical instabilities become important in FEM. It should be noted that our choice of a hemispherical resin will only affect the near-field flow around the resin. Beyond a certain distance, the resin can be considered a point source for protons and sink for potassium ions, respectively, and the details of its shape thus become irrelevant. This far field is the regime of interest to us.

At the edge of the domain, there is an “open boundary” for the hydrodynamic problem. This implies that there is no fluid momentum flux through the boundary. Since there is no convective momentum transport in the Stokes equations, there can be flow, but no stress normal to the boundary. This is a standard technique to model a piece of a domain that is embedded in a larger physical region, without modeling the full geometry, but allowing for the flow lines not to be closed within the domain. The unbound simulation domain (Fig. 6b) is the same as the top-bound domain, but replaces the top glass slide with a hemispherical (open-boundary) domain.

We now provide the expressions for the boundary conditions used in the FEM model.

• For all solute species, no-penetration conditions are imposed in the Nernst–Planck equation on the bottom/top of the cell

where n̂ is the unit normal to the boundary pointing into the fluid.

For the resin, we only impose no-penetration conditions for Cl⁻. The exchange of H⁺ and K⁺ is modeled via out- and influx on the resin, respectively. To be precise, we impose the following flux boundary condition

where k_ex is the ion-exchange rate coefficient, which we need to determine by fitting to the experimental data. Note that we have assumed that the exchange is determined entirely by the cation concentration close to the resin. This is probably valid as long as the H⁺ concentration inside the resin is much larger than the cation concentration outside, which is the case for a fresh resin.

At the outer edge of the domain, the orange line in Fig. 6a, we impose concentration profiles for the ions that are based on the Poisson–Boltzmann solution for a two plate geometry with height H and surface charge σ_wall. For the geometry of Fig. 6b, the solution to the Poisson–Boltzmann equation for a single plate was used.

• For the Poisson equation, we impose a constant surface charge density σ_wall on all solid surfaces via

where the electrostatic potential is assumed to be evaluated at the boundary. The surface charge density is obtained from the experimental zeta-potential measurement through the Grahame equation⁴²with ζ the zeta potential. Note that it is not clear what the most appropriate boundary conditions are for the resin surface, hence we chose the same boundary condition (eqn (13)) as on the other surfaces for computational convenience.

• The boundary conditions for the incompressible Stokes equation are no-slip boundary conditions on all solid surfaces (u⃑ = 0), that is the bottom and top of the cell, plus the resin. At the outer edge of the domain, a no-normal stress boundary condition is applied, which reads

with T denoting transposition and [Doublestruck I] the identity matrix.

A final detail for the FEM solver is that we used polynomial ansatz functions of order 2 for the electrostatic, order 2 for the diffusion/migration, and order 3 + 2 for hydrodynamic equations. This is necessary in order to further reduce spurious flows. Despite these measures, as well as decoupling the solute and solvent problems by our low Pe assumption, extremely fine meshes are required, see Fig. 7, that push the boundaries of modern computational platforms in order to obtain convergence and sufficiently smooth results within a reasonable time.

4.1 Tracer speed

We determined the tracer speed from the solution of the above system of electrokinetic equations with boundary conditions as follows. The speed

U_PS = u⃑·r̂ + μ_EE⃑·r̂, (16)

is comprised of an advective term, which is captured by u⃑, and a component deriving from the electrophoretic mobility μ_EE⃑. We evaluate the velocity and electric field at a constant “equilibrium height”, h*, where gravity balances electrostatic repulsion from the wall. Throughout, we used a constant height of h* = 5 μm. The exact height in the experiment is difficult to measure, presumably varies locally, and changes with the environment. We therefore varied h* between 4 μm and 10 μm to check how our specific choice affected the result. The resulting speed profiles turned out to be virtually the same in this range. This is because the fluid flow velocity, which is the major component in the tracer speed, varies over a typical length scale of order H ≫ h*. Note that eqn (16) treats the tracer particle as if it were a point-like object, i.e., it does not perturb the flow and electric fields by its presence. In general we found that including the second term in eqn (16) does not significantly modify the U_PS, leading us to conclude that advection indeed dominates over electrophoretic effects for the tracer motion.

4.2 Parameter choices

We made the following parameter choices to simulate the experimental system. For the geometry of the simulation setup we typically chose: r_R = 25 μm for the radius of the resin and R_sim = 3 mm for the radius of the cylindrically symmetric domain. This choice was a compromise between the size of the sample cell, which is too large to numerically simulate in its entirety, and a domain size on which the power-law decay in the fluid velocity was observable in the steady-state FEM calculations. The height of the domain was chosen to match the relevant experimental setup, e.g., H = 1 mm, with the open simulation domain representing the H = 10 mm domain, as we explain in Section 5.

The fluid represents water at room temperature (T = 298.15 K), which has a mass density of ρ_f = 1.0 × 10³ kg m⁻³, viscosity η = 8.9 × 10⁻⁴ Pa s, and relative permittivity ε_r = 78.4. The diffusivities of the ionic species are D_H⁺ = 9.3 × 10⁻⁹ m² s⁻¹ (ref. 43) and D_K⁺ = D_Cl⁻ = 2.0 × 10⁻⁹ m² s⁻¹ (ref. 36). The bulk concentration of impurities was chosen to be , in line with our estimates from the experiment, and or , with . The ion-exchange rate coefficient k_ex = 3.08 × 10⁻⁶ m s⁻¹ was obtained by fitting the near-field velocity to the experiment for H = 1 mm in Fig. 3. The surface charge density σ_wall = −4.03 × 10⁻⁴ C m⁻² was computed from the experimentally measured zeta potential ζ ≈ −0.1 V using the Grahame eqn (14), see ref. 42.

5 Finite-element results

The FEM-computed fluid flow for the steady-state problem is shown in Fig. 8 for a large portion of the sample cell; we used H = 1 mm. Both on the top and bottom wall the fluid flow is radially inward, due to the electroosmotic driving near the walls, with swirl-like patterns forming in the middle of the cell, due to the incompressibility of the fluid.

Fig. 8 Visualization of the FEM-calculated fluid flow in a H = 1 mm sample cell over a radial range of 2 mm. The resin is shown in gray in the bottom-left corner, the direction of the fluid flow is shown with arrows, and color indicates the magnitude of the local velocity.

Fitting the near-field tracer speed U_PS for a cation concentration of ρ = 1 μmol L⁻¹ to the H = 1 mm profile in Fig. 3, we found that k_ex = 3.08 × 10⁻⁶ m s⁻¹ is sufficient to match the experimentally observed near-field speed—compare Fig. 3 and 9a. We used this parameter throughout our simulations. This gives rise to an average surface flux of j_s ≈ 7 × 10⁻⁸ mol m⁻² s⁻¹, which corresponds closely to our back-of-the-envelope estimate in Section 2. This shows that the experimentally observed tracer speeds can indeed be explained by ion exchange of trace amounts of cationic impurities in the μmol L⁻¹ range.

First, we verified that our steady-state solution for the quasi-2D domain gives a reasonable result, when compared to the time-dependent simulations. We considered a cell height of H = 0.2 mm for this problem. This choice allowed us to reduce the number of mesh elements required compared to the typical experimental height H = 1.0 mm and thereby improve the computational time sufficiently to access second time scales. Fig. 9a shows several time-dependent tracer speed curves, up to the maximum time of t = 10 s that we could access with our FEM calculations (several days of computer run time). Note that for these times, the long-time, far-field power-law decay has not yet set in. We estimate the time for this decay to set in, using the time it takes H⁺ ions to diffuse a distance H: t_H⁺ = H²/D_H⁺ ≈ 400 s. This time is short on the time scale of the experiment, but too long to access via FEM calculations, which is why we consider analytic theory in Section 6. Nevertheless, the near-field solution has begun to converge to the steady-state after 10 s. Considering the relatively short time scales compared to the length of the experiment, on which convergence should take place, we are justified in neglecting the time-dependence in the FEM calculations.

Fig. 9 The tracer speed U_PS obtained using our numerical calculations for several configurations of the system. (a) The time-dependent solution for U_PS as a function of the radial distance r for several times t, for PS7 tracers in a sample cell with height H = 0.2 mm. The steady-state solution is given by the purple curve (t = ∞). The power-law decay sets in only for long times, see the far-right of the image. (b) The steady-state solution for several cell heights. The gray dashed lines serve as guides to the eye for the power-law decay.

Second, we considered the far-field r⁻¹ scaling in our steady-state simulations in Fig. 9b. For increasing H there is an increasingly large intermediate range of 3D decay with r⁻². This is to be expected, because the minimum length (equivalently the time, in the time-dependent problem) that the ions travel before the quasi-2D decay sets in will increase with the height. This also explains the experimental observation of r⁻² scaling for the H = 10 mm sample cell in Fig. 3, since for this height the transition time can be estimated to be t_H⁺ = H²/D_H⁺ ≈ 1.0 × 10⁴ s, which is long compared to the experimental time scale (typical velocity measurements take place within an hour after sample preparation).

When solving the steady-state system of equations on the geometry of Fig. 6a, the flow field always decays with r⁻¹ in the far field. This is because in the steady-state problem a sufficient amount of time has passed for the ions to “become aware” of their confinement, i.e., stationarity is analogous to t → ∞ in the time-dependent problem. Therefore, the experimentally obtained transition from quasi-2D to 3D decay by increasing the height H, cannot be observed in such a simulation. To observe a r⁻² decay in the steady-state far field, an unbound domain must be simulated instead. Hence the need for the geometry of Fig. 6b.

Finally, comparing Fig. 3 and 9b, we find that the regime in which the power-law decay sets in is much closer to the resin in the experiment, i.e., around 75 μm. This is counterintuitive, since on the basis of simple geometry arguments one would expect the ions to become aware of their 2D confinement when the distance they have diffused becomes comparable to the confining height. This suggests that there are mechanisms by which the ions are transported faster than through diffusion alone. A clear candidate is advection via the fluid flow, as the flow field around the resin, see Fig. 8, causes significant vertical displacement of the near-resin ions, provided the Péclet number is sufficiently large. We already estimated in Section 3 that this is likely the case. Unfortunately, the effect of advection cannot be incorporated in this work due to current limitations in computational performance for our FEM modelling, so that a quantitative match between theory and experiment is left for future study. Nevertheless, the qualitative behavior that we do capture together with our geometric arguments already provide important insights into the pumping mechanism.

6 Analytical calculations

In this section, we present an approximate, analytical solution to the electrokinetic equations on the domain of the sample cell. This allows us to obtain the time dependence and radial scaling of the flow in the far field. Our solution makes use of the equations provided in Section 3, where we already made the following simplifying assumptions: (i) Advection can be neglected. (ii) The only ions present are H⁺, K⁺, and Cl⁻. (iii) The diffusivities of K⁺ and Cl⁻ ions are equal. In order to make progress analytically, we also require: (iv) The perturbation of the ionic concentrations from their equilibrium distributions due to ionic fluxes generated by the resin bead are small compared to the background ionic concentrations—this is likely to be strictly valid only at short times or far from the central bead. (v) Far enough from the resin bead, the solution can be treated as 2D, and the bead as a point (δ-function) source. This means that our geometry is essentially a 2D disk, rather than the 3D cylinder segment of the FEM setup, with the fluxes of the species independent of the vertical position in the sample cell. (vi) The electrostatic Debye screening length κ is much shorter than the relevant length scales of the problem, which are of order H. (vii) The resin bead produces constant, equal and opposite fluxes of H⁺ outwards and K⁺ inward.

Under the above assumptions, we can combine eqn (2) and (3) to obtain

∂_tx_Cl⁻ = D_Cl⁻∇²_2D(x_Cl⁻ − ψ), (19)

where we define dimensionless concentrations x_i ≡ (ρ_i − ρ^∞_i)/ρ^∞_i, the dimensionless electrostatic potential ψ = Φe/(k_BT), and ∇²_2D is the 2D Laplacian. From the linear approximation (iv), we have kept only terms up to linear order in x_i and ψ, and from the no-advection approximation (i), we have neglected the term in u⃑. The final term in eqn (17) and (18) represents the steady production of H⁺ and consumption of K⁺ at the origin. Here, the 2D radial vector r⃑ = x⃑ + y⃑ (r = |r⃑|) and δ_2D is the 2D δ-function, which is normalized so that , with the integral running over the whole plane. Γ is the total production rate of H⁺ in molecules s⁻¹. Note that we do not make Γ dependent on ρ_K⁺, as in the FEM model, to avoid complicating our calculation.

Linearizing the Poisson equation, eqn (5) and defining the dimensionless background concentrations gives

where the inverse Debye length κ is given by eqn (6).

We now apply the thin-Debye-layer approximation (vi). Since we are interested in distances from the origin r ≫ H, the capillary height, this approximation can be quantified as κr ≫ 1, and the approximation involves making an expansion to lowest order in the small parameter 1/(κr). Now, from eqn (17)–(19), we must have that ψ and x_i are of similar order. However, if we multiply eqn (20) by r², we see that the right-hand side is of order unity since ∇²_2D = [scr O, script letter O](r⁻²), but the left-hand side is of order κ²r². This means that, for consistency, the sum on the left-hand side of eqn (20) must be zero to lowest order in 1/(κr), i.e.,

That is, the charge density is approximately zero everywhere outside a thin Debye layer close to the capillary surface. Note that this does not mean that ∇²_2Dψ = 0, as from eqn (20) it follows that the leading order term in ∇²_2Dψ is equal to the finite next-to-leading-order term in the charge density. For simplicity, we write α_H⁺ = α and, from the condition of charge balance in the background concentrations, α_K⁺ = 1/2 − α and α_Cl⁻ = 1/2. Then eqn (21) can be rewritten as

2x_H⁺α + x_K⁺(1 − 2α) = x_Cl⁻. (22)

Eqn (17)–(19) together with eqn (22) represent a linear system of equations that we will now solve. From linear combinations of eqn (17)–(19), we can eliminate ψ, and eqn (22) allows us to also eliminate x_Cl⁻. This leaves us with two equations for x_H⁺ and x_K⁺

We solve these equations using the ansatz functions

where A, B, C, and E are constants, and f₁(r,t) and f₂(r,t) are given by the time-integral of the Green's function for the 2D diffusion equationn,^20,44 which accounts for the constant point-source at the originwhere the D̃_m, m ∈ {1,2} are effective diffusivities to be determined. Directly solving for these effective diffusivities in general yields extremely unwieldy expressions, so instead we first make the simplifying assumption (iii) that D_K⁺ = D_Cl⁻. This gives

D̃₁ = D_K⁺; (28)

We then solve for the constant terms in eqn (26), yielding

A = 0; (30)

B = 1; (31)

Plugging these constants into eqn (25) and (26) and using x_i ≡ (ρ_i − ρ^∞_i)/ρ^∞_i allows us to obtain the time-dependent ion density profiles. To calculate ψ, we take the time derivative of eqn (22) and use this to eliminate the left-hand side from eqn (17)–(19) to obtain, after some algebraUsing Φ = k_BTψ/e this yields

The potential in eqn (35) will generate an equal slip-velocity on both the upper and lower surfaces of the channel; we are sufficiently far away from the resin that the asymmetry caused by it being glued to the bottom wall should not strongly affect the flow field. From eqn (35) we then obtain

Strictly speaking, u⃑_slip is the velocity at the outer edge of the Debye layer. However, in the thin Debye limit, we can take u⃑_slip to be the fluid velocity on the wall itself. In the bulk of the channel, the fluid flow will vary with a typical length scale H, the capillary height. Therefore, we can use u⃑_slip as a good estimate for the velocity u⃑_PS(z) of a tracer particle located a small distance z ≪ H above the capillary surface, with a fractional error [scr O, script letter O](z/H) < 10⁻² for the channels used here.

Fig. 10 shows the tracer speed U_PS as a function of the distance as calculated using eqn (36) for a sample cell with H = 0.2 mm and other parameters the same as in Fig. 9a. Since the theory does not use the same expression for the surface fluxes (vii), we selected a molecular exchange rate of Γ = 2 × 10⁸ s⁻¹ to give reasonable speed correspondence when compared to the results in Fig. 9a. At long times, i.e., for t ≫ r²/D̃₂, the exponential term approaches unity, so we obtain the 1/r scaling observed experimentally. At short times, the exponential term dominates, so we obtain the rapid radial decays seen in the FEM calculations in Fig. 9a.

Fig. 10 The tracer speed U_PS obtained using our analytic theory as a function of the radial distance r for several times t, for PS7 tracers in a sample cell with height H = 0.2 mm, pH = 7, and Γ = 2 × 10⁸ s⁻¹; this models the setup of Fig. 9a. The steady-state solution is given by the purple curve (t = ∞). The gray dashed line serves as a guide to the eye for the power-law decay.

Note that, as we would expect for a 2D system, the electrostatic field in eqn (35) (as well as the ion density profiles) is not in steady state—the integral approaches log(t) in the limit of large t. Nevertheless, the flow field in eqn (36) does approach a steady-state solution, see Fig. 10, because it scales with the gradient , which is in steady state. This is in line with our experimental observations, where we observed the same tracer-velocity trends after 24 h of pumping, albeit with decrease in speed by a factor of two. The latter can be attributed to depletion of the trace amounts of cations in the cell or a reduction in the effectiveness with which the resin exchanges ions. These results further underpin our conclusion that the microfluidic pumping is driven by ion exchange of trace amounts of cationic species in the sample cell.

The correspondence to the FEM calculations (Fig. 9a) is semi-quantitative, despite the additional simplifying assumptions. However, it should be noted that we cannot make a prediction for the near field using the theory, due to the quasi-2D assumption (v), nor do we account for advection of the solutes (i). This means that we cannot expect, nor do we find, quantitative agreement with the experiment.

7 Discussion and outlook

In summary, we have introduced and characterized an ion-exchange-resin-based microfluidic pump. The striking feature of this pump is that it operates in μmol L⁻¹ ionic concentrations for periods exceeding 24 h and yet manages to produce fluid flows with speeds of several μm s⁻¹ over hundreds of μm, without strongly modifying its environment.

We demonstrated that our pump uses trace amounts of cations to generate fluid flow using a combination of tracer velocimetry experiments, analytic electrokinetic theory, and finite-element-method simulations. These together show that fluid flow is achieved via electroosmosis, by the exchange of cations for protons in its interior. The difference in ionic mobility between the cations and the protons, for which they are exchanged, sets up a diffusion potential that points towards the resin and causes fluid flow in this direction. The speed of pumping can be modified by varying the bulk cation concentration in the μmol L⁻¹ range.

Our pump has several advantages over other microfluidic pumps that also exploit diffusion-potential-based electroosmosis, e.g., small pieces of salt that slowly dissolve.^8,22 Firstly, ion-exchange pumps only modify the nature of the bulk ions, not their net concentration. Secondly, as a pump dissolves, it might change shape, thereby inducing an undesirable directionality to the pumping. Our system does not have this disadvantage, as our spherical ion-exchange resins retain their shape throughout. Thirdly, ion exchange using protons as the exchangeable cation has the advantage of setting up significant diffusion potentials (and hence flow), due to the strong difference in diffusivity between the proton and any exchanged cationic species. Finally, the ion-exchange-resin pump functions for very long times in a low-ionicity medium—over 24 h—compared to the much shorter operating times of dissolving micropumps, which were indicated to be around 20 min in ref. 22.

The range of our pump can be tuned via the height of the sample cell to give rise to either 3D or quasi-2D decay of the far-field flow velocity, i.e., power-law decay with the relevant exponent. We have thus demonstrated that significant microfluidic pumping can be achieved at very low fuel (cation) concentrations and can be sensitively tuned via the geometry. This tunability can be exploited to self-assemble single colloidal crystals.²³

In our modeling of the experiment it proved necessary to ignore the advective contributions to the ionic fluxes. This simplifying assumption is the likely cause of the quantitative (but not qualitative) differences between the theory and experiment. We argue in favor of including advective (fluid flow) contributions to the ionic fluxes in any future modeling of these systems, as this is important in understanding the near-field fluid flow around the ion-exchange resin, and to extract the kinetics of ion exchange from far-field flow and concentration profile measurements.

Presently, we are only able to indicate that the ion-exchange process is likely diffusion limited in our system. Future experimental work will focus on pH measurements to quantify the exchange process, while the nature of the decay in these concentration profiles will be further examined using numerical approaches. For the latter, the use of a boundary-layer approach to rescale the high Péclet number regime and make these problems computationally tractable will be explored. Furthermore, capturing the near-field flow accurately will be relevant to understanding the formation of self-assembled cooperative swimmers based on mobile ion-exchange resins and tracer particles.²⁶

In conclusion, our system showcases the significance of very small ionic concentrations and fluxes in microfluidic settings. This suggests that such fluxes may be responsible for flow and motion in a much wider range of out-of-equilibrium systems, such as for chemical swimmers and in biological processes, and should be considered in future modeling thereof.

Acknowledgements

JdG gratefully acknowledges financial support from a Marie Skłodowska-Curie Intra European Fellowship (G.A. No. 654916) within Horizon 2020. JdG and CH, as well as RN, DB, and TP, further thank the DFG for funding through the SPP 1726 “Microswimmers—From Single Particle Motion to Collective Behavior” (HO1108/24-1 & Pa459/18-1).

References

1. H. Andersson, W. Van Der Wijngaart, P. Nilsson, P. Enoksson and G. Stemme, Sens. Actuators, B, 2001, 72, 259.
2. J. Santiago and D. Laser, J. Micromech. Microeng., 2004, 14, R35.
3. T. Kline, W. Paxton, Y. Wang, D. Velegol, T. Mallouk and A. Sen, J. Am. Chem. Soc., 2005, 127, 17150.
4. W. Paxton, P. Baker, T. Kline, Y. Wang, T. Mallouk and A. Sen, J. Am. Chem. Soc., 2006, 128, 14881.
5. M. Ibele, Y. Wang, T. Kline, T. Mallouk and A. Sen, J. Am. Chem. Soc., 2007, 129, 7762.
6. S. Chang, V. Paunov, D. Petsev and O. Velev, Nat. Mater., 2007, 6, 235.
7. A. Nisar, N. Afzulpurkar, B. Mahaisavariya and A. Tuantranont, Sens. Actuators, B, 2008, 130, 917.
8. M. Ibele, T. Mallouk and A. Sen, Angew. Chem., Int. Ed., 2009, 48, 3308.
9. I.-K. Jun and H. Hess, Adv. Mater., 2010, 22, 4823.
10. A. Shields, B. Fiser, B. Evans, M. Falvo, S. Washburn and R. Superfine, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 15670.
11. Y. Hong, M. Diaz, U. Córdova-Figueroa and A. Sen, Adv. Funct. Mater., 2010, 20, 1568.
12. T. Hogg and R. Freitas, Nanomedicine, 2010, 6, 298.
13. A. Solovev, S. Sanchez, Y. Mei and O. Schmidt, Phys. Chem. Chem. Phys., 2011, 13, 10131.
14. H. Zhang, K. Yeung, J. Robbins, R. Pavlick, M. Wu, R. Liu, A. Sen and S. Phillips, Angew. Chem., Int. Ed., 2012, 51, 2400.
15. V. Yadav, H. Zhang, R. Pavlick and A. Sen, J. Am. Chem. Soc., 2012, 134, 15688.
16. A. Reinmüller, E. Oğuz, R. Messina, H. Löwen, H. Schöpe and T. Palberg, J. Chem. Phys., 2012, 136, 164505.
17. A. Farniya, M. Esplandiu, D. Reguera and A. Bachtold, Phys. Rev. Lett., 2013, 111, 168301.
18. S. Sengupta, D. Patra, I. Ortiz-Rivera, A. Agrawal, S. Shklyaev, K. Dey, U. Córdova-Figueroa, T. Mallouk and A. Sen, Nat. Chem., 2014, 6, 415.
19. M. Esplandiu, A. Afshar Farniya and A. Bachtold, ACS Nano, 2015, 9, 11234.
20. I. Ortiz-Rivera, H. Shum, A. Agrawal, A. Sen and A. Balazs, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 2585.
21. C. Zhou, H. Zhang, Z. Li and W. Wang, Lab Chip, 2016, 16, 1797.
22. J. McDermott, A. Kar, M. Daher, S. Klara, G. Wang, A. Sen and D. Velegol, Langmuir, 2012, 28, 15491.
23. R. Niu, E. C. Oğuz, H. Müller, A. Reinmüller, D. Botin, H. Löwen and T. Palberg, Phys. Chem. Chem. Phys., 2017, 1 10.1039/C6CP07231C.
24. D. Velegol, A. Garg, R. Guha, A. Kar and M. Kumar, Soft Matter, 2016, 12, 4686.
25. J. L. Anderson, Annu. Rev. Fluid Mech., 1989, 21, 61.
26. R. Niu, D. Botin, A. Reinmüller and T. Palberg, 2016, arXiv 1610.04723, 1.
27. E. Uzgiris, Rev. Sci. Instrum., 1974, 45, 74.
28. E. Uzgiris, Prog. Surf. Sci., 1981, 10, 53.
29. B. Midmore, G. Pratt and T. Herrington, J. Colloid Interface Sci., 1996, 184, 170.
30. P. Hough, International conference on high energy accelerators and instrumentation, 1959.
31. A. Brown, W. Poon, C. Holm and J. de Graaf, Soft Matter, 2017, 1 10.1039/C6SM01867J.
32. A. Brown and W. Poon, Soft Matter, 2014, 10, 4016.
33. T. Palberg, T. Köler, B. Sieber, H. Schweinfurth, H. Reiber and G. Nägele, J. Phys.: Condens. Matter, 2012, 24, 464109.
34. Á. Delgado, F. González-Caballero, R. Hunter, L. Koopal and J. Lyklema, J. Colloid Interface Sci., 2007, 309, 194.
35. F. Millero, Geochim. Cosmochim. Acta, 1995, 59, 661.
36. H. Harned and R. Nuttall, J. Am. Chem. Soc., 1949, 71, 1460.
37. U. Córdova-Figueroa and J. Brady, Phys. Rev. Lett., 2008, 100, 158303.
38. G. Rempfer, G. Davies, C. Holm and J. de Graaf, J. Chem. Phys., 2016, 145, 044901.
39. G. Rempfer, S. Ehrhardt, C. Holm and J. de Graaf, Macromol. Theory Simul., 2016, 1 DOI:10.1063/1.4968596.
40. G. Rempfer, S. Ehrhardt, N. Laohakunakorn, G. Davies, U. Keyser, C. Holm and J. de Graaf, Langmuir, 2016, 32, 8525.
41. P. Kreissl, C. Holm and J. de Graaf, J. Chem. Phys., 2016, 144, 204902.
42. D. Grahame, Chem. Rev., 1947, 41, 441.
43. CRC handbook of chemistry and physics, ed. W. M. Haynes, CRC Press, Boca Raton, USA, 93rd edn, 2013.
44. S. Shin, E. Um, B. Sabass, J. Ault, M. Rahimi, P. Warren and H. Stone, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 257.

---

Footnotes

† Electronic supplementary information (ESI) available: Movies for the ion-exchange-resin pumps. See DOI: 10.1039/c6sm02240e

‡ These authors contributed equally to this work.

§ The eq. stands for amount of charge that the resin can exchange (here per liter of resin): 1 eq. = 1 mol of monovalent ions, 0.5 mol of divalent ions, etc.

¶ We used Doppler velocimetry with PS tracers,³³ originally designed to measure bulk electrokinetics in colloidal suspensions, to determine the electroosmotic mobility μ_wall of cleaned glass slides using a custom made cell with exchangeable sides for the top and bottom. Standard electrokinetic theory was used to calculate the zeta potential from the mobility.³⁴

|| We experimentally measured a tracer speed decrease of only a factor 2 over a 24 h period, justifying the assumption of almost constant pumping. The shape of the speed profile remained unchanged.

** We believe that in practice the Pe number is likely to be self-limiting to a value ≈1, as found for chemically propelled swimmers.³⁷ This is because the high-concentration-gradient region centered around the colloid would be expelled into the bulk of the channel by a strong advective current—see the direction of the flow lines in Fig. 8—where it would no longer contribute strongly to electrophoretic flow generation.

†† Making the Pe = 0 assumption allows us to split the solute and solvent problems and solve them in series, rather than in parallel, leading to a strong reduction in the required mesh resolution and therefore of computer memory.

‡‡ There will be similar association–dissociation reactions involving dissolved CO₂ in water, as well as other species. The arguments provided here apply equally to these reactions.

---