Abstract

Colloidal spheres with a partial platinum surface coating perform auto-phoretic motion when suspended in hydrogen peroxide solution. We present a theoretical analysis of the self-propulsion velocity of these particles using a continuum multi-component, self-diffusiophoretic model. With this model as a basis, we show how the slip-layer approximation can be derived and in which limits it holds. First, we consider the differences between the full multi-component model and the slip-layer approximation. Then the slip model is used to demonstrate and explore the sensitive nature of the particle’s velocity on the details of the molecule-surface interaction. We find a strong asymmetry in the dependence of the colloid’s velocity as a function of the level of catalytic coating, when there is a different interaction between the solute and solvent molecules and the inert and catalytic part of the colloid, respectively. The direction of motion can even be reversed by varying the level of the catalytic coating. Finally, we investigate the robustness of these results with respect to variations in the reaction rate near the edge between the catalytic and inert parts of the particle. Our results are of significant interest to the interpretation of experimental results on the motion of self-propelled particles.

Over the past decade, research into self-propelled particles (SPPs) has undergone rapid development and attracted strong interest from the scientific community, due to the inherent out-of-equilibrium nature of their behavior. [1, 2] Examples of naturally occurring self-propelled ‘particles’ include humans, [3, 4, 5] birds, [6] fish, [7] insects, [8, 9] spermatozoa, [10, 11, 12] bacteria, [13, 14, 15] and algae. [16, 17] The principle of autonomous movement is thus relevant across many orders of magnitude in size and speed. Some of the most interesting behavior of SPPs is observed when the particles are colloidal in size (∼1 nm – 1 μm) and suspended in a liquid medium, since there will be a competition between Brownian (thermal) motion and the self-propulsion mechanism. [18] Moreover, the viscous nature of the suspending medium strongly impacts the strategies by which autonomous movement may be achieved mechanically. For low-Reynolds-number swimmers this movement must be nonreciprocal, as put forward by Purcell in his scallop theorem. [19]

There exists another set of colloidal SPPs that differs strongly from the aforementioned microorganisms, namely, artificial particles that move by self-generated solute gradients. Artificial SPPs were pioneered by the work of Ismagilov et al.[20] and Paxton et al.[21] These early studies have led to the development of a wide range of artificial SPPs that move by utilizing a variety of gradient-based propulsion mechanisms, including: bimetallic nanorods that move by self-electrophoresis, [21, 22] Pt-coated colloidal spheres moving by self-diffusio- or self-electrophoresis, [23, 24, 25, 26, 27, 28] Au-coated colloidal spheres that move by thermophoresis [29, 30] or by thermally induced local mixing-demixing transitions, [31] hollow cones that expel oxygen bubbles, [32, 33] and many others. A common theme for most gradient-based SPPs is the decomposition of hydrogen peroxide at a platinum surface, which acts as a catalyst to the reaction. These particles are considered ideally suited as model systems for out-of-equilibrium phenomena [2] and may be utilized as functional components in micro-fluidic devices, [34] or be employed for medical purposes. [35] However, there are many open questions with regards to the way these particles achieve their motion [23, 24] and how this motion may be most effectively controlled, [36, 37, 38] that must be answered to allow for the successful implementation of these particles in real-world applications.

From the modeling perspective, a lot of attention has gone into the description of the self-propulsion mechanism of these particles. In particular, it is accepted that self-electrophoresis plays a dominant role for the self-propulsion of bi-metallic nanorods. [39, 40] The mechanism by which Pt-coated colloidal spheres achieve motion in hydrogen-peroxide solution is, however, not as well understood. This has led to heated debate on the type of phoretic motion by which these particles self-propel. [23, 24] The traditional view is that self-diffusiophoresis can be used to explain the experimentally observed motion in these systems, [26, 25] which has resulted in a large number of studies into the specifics of this mechanism.

The majority of the theoretical models are based on a coarse-grained, continuum-level approach that employs the dilute-limit slip-layer formalism, see, e.g., Refs. [41, 42]. From a simulation/computation perspective, particle-based models are more commonly used. [43, 44, 27, 45, 46, 47] Investigations into the properties of the self-diffusiophoretic model have considered the influence of a large number of parameters. For instance, the influence of the reaction mechanism on the size-velocity dependence of the SPPs was matched to experiment. [26, 25] In addition, there have been more fundamental studies which unified the continuum- and particle-based descriptions [42, 48] and considered the influence of advective and reactive processes on the self-propulsion. [49] Moreover, the influence of particle shape [50, 45, 51] and catalytic coating [52, 53, 49] have been considered, as well as the nature of the interaction potential. [54, 53]

A common theme in continuum modeling is that typically only the ‘effective’ interaction between (one of) the solutes and the surface of the particle is taken into account for the slip-based self-propulsion. In addition, the interaction between the solute and the inert, as well as the solute and the catalytic surface area is often assumed the same. There are theoretical studies that have touched upon such differences, see, e.g., Refs. [55, 56, 53], but a full analysis is still lacking. The use of an identical interaction potential for the inert and catalytic surface in theory stands in stark contrast with the more common particle-based simulation models, wherein anisotropy of the solvent-surface interaction is often a crucial ingredient to achieve movement. [43, 44, 27, 45, 47] Furthermore, from a physical perspective, the assumption of a homogeneous molecule-surface interaction is questionable.

In this manuscript, we solve a set of multi-component diffusion-advection equations [57] with appropriate boundary conditions [41, 42] to describe the self-diffusiophoretic movement of a self-propelled particle (SPP). We depart from the traditional single-component dilute-limit approximation in order to maintain momentum conservation and incompressibility of the total fluid at high fuel (hydrogen-peroxide) and product (water and oxygen) concentrations. We derive the continuity equations governing the dynamics of the individual species from a chemical potential and show the approximations that lead to a self-consistent model. The role of the solvent-surface potential is elucidated and it is shown how the slip-layer approximation follows from these equations. In addition, we consider the differences with previously established multi-component models. Our advection-diffusion-reaction model is subsequently used to analyze the motion of a catalytically coated colloid suspended in hydrogen-peroxide solution. In particular, we investigate the sensitivity of the self-propulsion velocity on the surface-molecule interaction potential and the level of catalytic coating for physically reasonable parameters. This part of the investigation is in a similar vein as the work of Sabass et al.[54], but takes a different approach to quantifying the effect of the surface-molecule interaction by coupling it to the level of catalytic coating.

We find a strong sensitivity of the SPP’s velocity on the molecule-surface interaction by taking into account the gradients of all the species involved in the reaction and by assigning different interaction potentials between the molecules and the inert and catalytic side of the particle, respectively. In contrast to the result of Refs. [50, 53], we show that the typical dependence of the swimming velocity on the level of catalytic surface coating is asymmetric in this parameter. Moreover, we find a specific set of parameters for which the SPP does not move when it is half coated, and moves in the direction of the catalytic cap or away from it, when it is less than or more than half coated, respectively. We also introduce a method to quickly determine the velocity dependence of a SPP for small perturbations of the molecule-surface interaction potential for low Péclet numbers. Finally, we briefly consider the sensitivity of the self-propulsion velocity on the details of the local reactivity of the catalytic cap and we find that for reasonable parameters our results are robust with respect to sizable changes. Our findings have significant implications for the realization of SPP-based applications and interpretation of experimental results.

The remainder of the document is structured as follows. In Section 2 we introduce our multi-component model for bulk fluids. The boundary conditions required to achieve self-propulsion for arbitrary shapes are discussed in Section 3. This general discussion is exemplified on the basis of a spherical Janus SPP suspended in hydrogen peroxide solution in Section 4. In Section 5 the dilute limit of the multi-component self-propulsion model is analyzed and the slip model is recovered. This analysis is built upon in Section 6, where for low Péclet numbers, we introduce a method to determine the influence of small changes in the molecule-surface interaction potential on the self-propulsion velocity. We subsequently present our results in Section 7, which is partitioned into several subsections. We introduce our system parameters and numerical solution approach in Section 7.1. This is followed by a comparison between the full multi-component result and the dilute-limit slip model result for the motion of a colloidal SPP in Section 7.2. We then analyze the nature of the flow field around the SPP and the dependence of its speed on the catalytic surface coverage in Section 7.3. This analysis is extended upon in Section 7.4, where we consider the effect of heterogeneities in the surface-molecule interaction on the SPP’s speed. Our last result is presented in Section 7.5, wherein we briefly consider the influence of variations in the reactivity of the catalytic coating. Finally, we summarize and discuss our results in Section 8 and present an outlook.

In the following, we consider a three component mixture, consisting of hydrogen peroxide H2O2 (fuel), oxygen O2 (reactant), and water H2O (medium + reactant). However, in this paragraph we will start with a more general formulation of the multi-component model.

2.1 Remarks on Multi-Component Modeling

The choice of a multi-component description for the fluid is related to the composition of the ‘total fluid’ typically used in experiments, which contains at least three species in the case of a hydrogen-peroxide solution. In addition, most current theoretical descriptions do not take properly into account the differences in the molecular masses of the species. Finally, at the high concentrations of fuel used in experiment, typically up to 10% b.v. of H2O2, [23, 24, 25, 26, 27, 28] make it unclear how the incompressibility constraint on the total fluid is preserved within the diffusion-advection formulation of the behavior of the individual species.

An accurate description of the physics in a multi-component system uses incompressibility as a constraint that alters the behavior of the individual components. The Maxwell-Stefan multi-component formalism is such a model, that describes a total fluid that satisfies the Navier-Stokes equations, comprised of separate species which satisfy complex diffusion-advection relations that ensure conservation of local fluid density. [58, 57, 59] This model is, however, exceedingly complex and leaves many free parameters. For instance, the values of the species’ cross-diffusion coefficients are typically not known and not measurable in experiments. In this section, we therefore present a reduced multi-component model that takes into account some of the complexities in this type of system, whilst leaving the number of free parameters to a minimum.

We assume a bulk liquid of infinite extent; boundary conditions will be imposed at a later stage. All components that comprise the fluid are fully miscible and no demixing takes place, for example, dissolved oxygen cannot cross the solvation threshold to form bubbles. The latter would require, e.g., the addition of a Cahn-Hilliard-type description of the mixing-demixing transition, which needlessly complicates the system for our purposes.

2.2 Incompressibility and the Fluid Density Profile

In the following we will consider a multi-component fluid for which the mass and momentum transport is described using the Navier-Stokes equations. For any flow characterized by the Navier-Stokes formalism (both compressible and incompressible), the continuity equation of the fluid is given by

∂∂tρ(r)+∇⋅(ρ(r)u(r))=0,

(1)

where ρ(r) is the mass density, ∂/∂t denotes partial differentiation with respect to time, ∇ is the spatial differential operator, ⋅ denotes the dot product (such that ∇⋅ is the divergence), u(r) is the advective velocity, and the time dependence of the quantities is implicit. The incompressibility condition is given by

This equivalence leads to a subtlety in the definition of incompressibility. Note that (∂/∂t)ρ(r)=0 and ∇ρ(r)=0 do not have to be satisfied simultaneously, only Eq. (3) has to hold for a flow field to be incompressible. That is, while the flow of a material can be incompressible, the material itself does not need to have a homogeneous density distribution.

In the following we will choose to impose a homogeneous density distribution of the material, i.e., ρ≡ρ(r). This is a reasonable assumption since we can prepare our system with a homogeneous density distribution and our boundary conditions all conserve local density (reactions are mass conserving), as we will see. Therefore, in the initial state we find ∇ρ(r)=0; consequently incompressibility guarantees (through (∂/∂t)ρ(r)=0) that the density remains homogeneous. In the derivation below this condition may be relaxed, however, the final result will not be substantially affected.

2.3 Chemical Potential and Diffusive Flux

We start our derivation from the expression for the chemical potential [60] of the species that comprise a non-ideal, multi-component mixture

μk(r)=μ0k+kBTlogγk(r,{xl})+kBTlogxk(r)+Φk(r),

(4)

where the subscript k is used to identify the k-th species out of N total species; μ0k is the reference chemical potential of that species; kB is the Boltzmann constant and T is the temperature; γk(r,{xl}) is the activity coefficient of species k with r the position coordinate; Φk(r) is an external potential acting on molecules of species k; and xk(r) is the molar fraction of species k, which is defined as

xk(r)≡nk(r)∑Nj=1nj(r),

(5)

with nk(r) the particle density of species k. The activity coefficient accounts for non-ideality and is therefore dependent on the concentrations of the other species as well, which is denoted here by the use of the set {xl}. In the following we specify γk(r,{xl}) to be constant in r; this will necessitate further approximations to be made in order to recover incompressibility of the total fluid in the following. The first strong reduction is therefore made here, as the fluid is considered an ideal mixture; effectively we set γk(r,{xl})=1.

In order to derive the continuity equations for the k-th species, we need to introduce several quantities. The mass density of the total fluid is given by ρ and is the sum of the component’s mass densities ρk(r). We may now introduce the mass fractions of the species as

ωk(r)≡ρk(r)ρ.

(6)

Species k has a molecular mass of mk, which gives us the equivalence ρk(r)=mknk(r) and allows us to write

m(r)≡(N∑k=1ωk(r)mk)−1,

(7)

if we assume that there is a field m(r) such that ρ=m(r)n(r), with n(r) the total number density. Here, m(r) is a quantity that has the properties of a number-density averaged molecular mass. That is, an average molecular mass of the fluid, based on the composition of the fluid into species with different molecular masses. Note that if all molecular masses are equal, then m(r) is position (and time) independent, since the sum over ωi(r) is unity by definition. Because we are interested in an incompressible total fluid, it is convenient to recast the expression for the molar fractions in terms of the above quantities

xk(r)=ωk(r)m(r)mk.

(8)

This results in the following expression for the gradient (denoted by ∇) of the molar fraction

∇xk(r)=m(r)mk∇ωk(r)+ωk(r)mk∇m(r).

(9)

Using the above quantities and Eq. (9), it follows that the expression for the number density flux J∗k(r) of species k, in the frame of reference that is co-moving with the total fluid, can be written as

J∗k(r)

=

νknk(r)fk(r)

(10)

=

−νknk(r)∇μk(r)

(12)

=

−Dknk(r)kBT∇μk(r);

=

Extra open brace or missing close brace

−Dkρmk(∇ωk(r)+ωk(r)∇~Φk(r)),

where νk is the mobility of the k-th species, which is related to that species’ diffusion constant Dk via the Einstein-Smoluchowski relation νk=Dk/kBT; fk(r) is the force per particle acting on the component; and ~Φk(r)≡Φk(r)/kBT. Here, we used the first order expansion for small deviations of μk(r) from μ0k, in order to write fk(r)=−∇μk(r) in Eq. (12). [59]

2.4 Mass Flux and Continuity

The expression in Eq. (12) allows us to write the kinematic mass flux (in the co-moving frame) as

j∗k(r)

=

mkJ∗kρ

=

−Dk(∇ωk(r)+ωk(r)∇m(r)m(r)+ωk(r)∇~Φk(r)).

Note that the kinematic mass flux has the physical dimension of velocity. Further note that when all species have the same molecular mass, the mass-correction (second) term drops out of the above equations, which allows us to recover the traditional Fickian diffusion.

The kinematic mass flux in the co-moving frame j∗k(r) and the kinematic mass flux in the laboratory frame jk(r) are related according to the following transformation

jk(r)=j∗k(r)+ωk(r)u(r),

(14)

where u(r) is the velocity of the total fluid, which accounts for the advective contribution. The time-dependent continuity equations for a system, in which there are no bulk reactions, read

∂∂tωk(r)+∇⋅jk(r)=0.

(15)

Note that bulk reaction terms may be straightforwardly accounted for in our formalism by introducing source and sink terms into Eq. (15). The time-dependent continuity equation for the k-th species in the co-moving frame now follows from the relation in Eq. (14)

∂∂tωk(r)+∇⋅j∗k(r)+u(r)⋅∇ωk(r)+ωk(r)∇⋅u(r)=0,

(16)

The above set of equations corresponds to the simplified diffusion-advection model put forward in Ref. [58].

2.5 Momentum Transport

For momentum transport in the total fluid, we are interested in the low Reynolds number limit. That is, Re≪1, where Re≡au/η, with a the relevant length scale of the problem, u the speed of the object or fluid, and η the dynamic viscosity of the fluid. For colloidal SPPs typically used in experiments, Re≪1 is satisfied. [23, 24, 25, 26, 27, 28] This implies that we can work with the linearized form of the Navier-Stokes equations for the total fluid

ρ∂∂tu(r)

=

−∇p(r)+ηΔu(r)+f(r);

(17)

∇⋅u(r)

=

0,

(18)

where p(r) is the pressure, Δ the tensorial Laplacian, and f(r) the force per volume acting on the fluid. Equation (18) gives the incompressibility constraint for the total fluid. The force term f(r) derives from the particle fluxes of the individual species [59]

f(r)=N∑k=1nk(r)fk(r)=N∑k=1J∗k(r)νk=N∑k=1ρj∗k(r)mkνk.

(19)

Using Eq. (LABEL:eq:fluxsp) it now follows that

f(r)

=

−kBTρ(∇1m(r)+1m(r)∇m(r)m(r))

(20)

−kBTN∑k=1ρωk(r)mk∇~Φk(r);

=

−kBTN∑k=1ρωk(r)mk∇~Φk(r).

The first two terms exactly cancel each other for an incompressible fluid. This is a nice feature of the multi-component model that takes into account the molecular mass of the species and utilizes incompressibility, which is missing in the traditional dilute-limit approximation.

2.6 Incompressibility and Flux

The continuity equation of the total fluid, see Eq. (18), should also follow by summing the continuity equations of the individual species, see Eq. (16), in order for the model to be self-consistent in the incompressibility assumption. From the properties of ωk(r) it follows that a sufficient condition to obtain incompressibility of the total fluid from the individual continuity equations, see Eq. (16), is given by

N∑k=1j∗k(r)=0.

(21)

This is also physically reasonable as otherwise there would be mass flow with respect to the frame co-moving with the fluid (local center of mass).

To ensure that Eq. (21) is satisfied, constraints on the species’ fluxes must be imposed, as the current expressions do not necessarily result in incompressibility of the total fluid. In particular, there is insufficient coupling between the individual mass fractions ωk(r) to ensure that this constraint is satisfied. This is typically the point where cross-diffusion terms are introduced in order to recover incompressibility. [59, 58, 57] However, the introduction of cross-diffusion terms leads to complicated relations between the components and diffusion coefficients that cannot be (easily) extracted from experimental data.

Here, we propose a different route to salvage the ‘simple’ Fickian diffusion model. We first consider the force-free part of the flux. We may write this flux as

j∗,0k(r)=−Dk(∇ωk(r)+ωk(r)∇m(r)m(r)).

(22)

If we impose

j∗,01(r)≡−N∑k=2j∗,0k(r),

(23)

then we obtain

N∑k=1j∗,0k(r)=0.

(24)

Here, we assumed that one of the components has a mass fraction that is substantially greater than the other components, namely the solvent. This component, say it is labeled k=1, is singled out and employed to cancel momentum, flux, and density related inconsistencies, when there are no forces acting. We must still take care of contribution coming from the potentials.

To ensure zero net flux in the co-moving frame, some correction force fcorr(r) must be added to the gradients of the potentials. We make the following ansatz for the total force ftot,k(r) acting on particles of the k-th species in the co-moving frame

ftot,k(r)=−∇~Φk(r)−fcorr(r).

(25)

By plugging in ftot,k(r) into Eq. (LABEL:eq:fluxsp) and utilizing Eq. (23), the following relation can be derived for the correction force such that Eq. (21) is satisfied and that incompressibility of the total fluid is recovered

fcorr(r)=−∑Nk=1Dkωk(r)∇~Φk(r)∑Nk=1Dkωk(r),

(26)

We may thus write for the co-moving kinematic mass flux

j∗k>1(r)

=

−Dk(∇ωk(r)+ωk(r)∇m(r)m(r)

(27)

−Dk(−ωk(r)ftot,k(r));

j∗1(r)

=

−N∑k=2j∗k(r);

(28)

ftot,k(r)

=

−(∇~Φk(r)−∑Nk=1Dkωk(r)∇~Φk(r)∑Nk=1Dkωk(r)).

Note that Eq. (28) now follows from the previously imposed condition of Eq. (23) and the properties of ftot,k(r).

2.7 Further Properties and Reductions

The above flux model has the desirable property that throughout space, only the continuity equations for N−1 of the species have to be solved simultaneously, as is expected for an incompressible system. Moreover, the potential acting on the solvent is not ignored, as would be the case if the condition of Eq. (28) had been imposed without modification of the force acting on the fluid. Note that this modification is only appropriate for the flux in the co-moving frame. The force acting on the fluid, see Eq. (20), should not be modified. This becomes clear by considering the case wherein all Dk and mk are the same, for which imposing the correction would lead to a zero net force in the momentum transport equation.

Since we have shown that the continuity equation for the total fluid can be derived from the individual continuity equations of the species for the above flux expression, we can now split these two parts, writing

∇⋅u(r)=0

(30)

for the fluid and

∂∂tωk(r)+∇⋅j∗k(r)+u(r)⋅∇ωk(r)=0,

(31)

for the diffusive species, respectively.

2.8 Overview of the Model

Summarizing, the multi-component model for the bulk fluid satisfies the following equations under the assumption of a homogeneous incompressible medium

ρ∂∂tu(r)

=

−∇p(r)+ηΔu(r)

(32)

−kBTρN∑k=1ωk(r)mk∇~Φk(r);

∇⋅u(r)

=

0;

(33)

∂∂tωk(r)

=

−∇⋅j∗k(r)−u(r)⋅∇ωk(r);

(34)

j∗k>1(r)

=

−Dk(∇ωk(r)+ωk(r)∇m(r)m(r)

(35)

−Dk(−ωk(r)ftot,k(r));

j∗1(r)

=

−N∑k=2j∗k(r);

(36)

ftot,k(r)

=

−(∇~Φk(r)−∑Nk=1Dkωk(r)∇~Φk(r)∑Nk=1Dkωk(r));

(37)

m(r)

=

(N∑k=1ωk(r)mk)−1.

(38)

This set of equations is fully consistent with the incompressibility assumption for the total fluid and dependent on a minimal number of fields: u(r), p(r), ωk(r), and ~Φk(r); and parameters: T, η, ρ, Dk, and mk.

Note that this model presents a significant departure from the typical dilute-limit approximation, since the difference in mass of the molecular species is taken into account and incompressibility is incorporated in a self-consistent manner. However, our modifications and underlying assumptions require the solute species to remain dilute in the solvent. It is therefore not of the same class as the more realistic, albeit far more complex Maxwell-Stefan formulation for higher solute concentrations.

Thus far, we have only concentrated on the description of the ‘bulk’ fluid, not the mechanism by which self-propulsion is achieved. To obtain self-diffusiophoretic motion of a particle in our multi-component description we require only a set of boundary conditions. We are interested in the stationary state of the system and we therefore solve the time-independent variants of Eqs. (32-38). By solving these equations simultaneously, the self-propulsion velocity and induced fluid flow profile in the stationary state may be obtained.

At the surface of the particle we require a no-slip boundary condition. That is,

u(s)=0,

(39)

where s is a point on the surface. We also require a flux boundary condition on (part of) the surface to cause the system to go out of equilibrium. For a multi-component system there can be a large number of (non-linear) reactions taking place simultaneously between the surface and the various components. The only requirement on the reactions is that there is no net mass flux into the surface

N∑k=1jk(s)⋅^n(s)=0,

(40)

where ^n(s) is the normal vector to the surface at the point s. If there were a net mass flux into the surface, the no-slip boundary condition of Eq. (39) would be violated, due to the relation between jk(r) and u(r). In physical terms, the condition of Eq. (40) imposes that the reactions are mass conserving. The choice for the reaction scheme is otherwise completely free.

There are a few caveats to the formulation of the system we presented above.

There is no fluid motion in the stationary state when all ~Φk(r)=0. That is, in the absence of an interaction between the surface of the SPP and the solvent/solutes, there is no back-coupling to the total fluid of the heterogeneous solute distribution caused by the surface reaction. This leads to zero fluid flow, in accordance with the results of Ref. [42].

Solving the time-independent variants of Eqs. (32-38) requires the existence of a stationary state, for which there is an inertial transformation between the frame in which the particle moves and the frame in which the particle is stationary and the fluid moves.

The existence of a inertial transformation to a co-moving frame implies that the particle performs rectilinear motion that is unaccelerated. This imposes restrictions on the shape of the particle and surface coating heterogeneity. Namely, the particle must be such that it does not experience a torque. A sufficient condition is that the particle is rotationally symmetric in the shape and surface coating. Moreover, there are restrictions on the reaction mechanisms that are permitted, e.g., the condition of stationarity forbids oscillating reactions.

Finally, some care needs to be taken in setting up the boundary conditions far away from the colloid, to ensure that a stationary state with non-zero self-propulsion velocity may be achieved, i.e., fuel must not be depleted.

In this section, we concentrate exclusively on a single spherical SPP of radius a, that is suspended in a three-component mixture of water, hydrogen peroxide, and oxygen (N=3), to better illustrate our auto-phoretic model. This colloid is partially coated with platinum to allow for a catalytic decomposition reaction of the hydrogen peroxide to take place; the uncoated part of the particle is assumed to be nonreactive. The sphere is centered in and fixed to the origin of our coordinate system. We use spherical coordinates with r the distance measured to the origin, θ the azimuthal angle, and ϕ the polar angle. The system is rotationally symmetric around the z-axis and ϕ is measured with respect to this axis, whereas θ describes rotations around it. Let the normal vector to the sphere’s surface (radial unit vector) be denoted by ^n and the polar-angle tangent vector by ^t; the Cartesian unit vectors are given by ^x, ^y, and ^z.

Using these notations the boundary conditions at the surface of the sphere are as follows. We include a flux boundary condition to describe a simple linear reaction 2 H2O2→ 2 H2O + O2 for hydrogen peroxide decomposition at the platinum surface. The reaction kinetics satisfy

∂∂t[H2O2]

\lx@stackrelPt=

−k[H2O2];

(41)

∂∂t[H2O]

\lx@stackrelPt=

−∂∂t[H2O2];

(42)

∂∂t[O2]

\lx@stackrelPt=

−12∂∂t[H2O2],

(43)

where k is the reaction rate. The set of equations, Eqs. (41-43), is used to treat bulk catalytic decomposition and needs to be converted to a surface model, considering the specific geometrical properties of the thin Pt-coating of the SPP, which we will do shortly. Other choices for the reaction mechanism by which hydrogen peroxide is decomposed at a platinum surface are possible, as set out in Ref. [40]. A multiple-step reaction mechanism through a variety of intermediate Pt-complexes is a more probable route for the decomposition. [61, 62] However, there is some indication that the rate-limiting step in the reaction process is linear in the hydrogen-peroxide concentration. [40] Taking the above considerations into account, we may write for the flux boundary conditions

ji(r)⋅^n|r=a=3∑k=1R(ϕ)aikωk(r)|r=a,

(44)

where the aik are stoichiometric/mass coefficients and R(ϕ) is the surface reaction rate. The relation to the bulk reaction rate k is given by R(ϕ)=k(ϕ)/Sc, where k(ϕ) is the bulk reaction rate and Sc is the catalytic surface area for which this rate is achieved. The ϕ-dependence of the reaction rate allows us to treat one part of the SPP’s surface as inert and one as catalytic. The stoichiometric/mass coefficients are conveniently summarized using the following matrix notation

A=[aij]=⎛⎜
⎜
⎜
⎜
⎜
⎜⎝mH2OmH2O200−100mO22mH2O200⎞⎟
⎟
⎟
⎟
⎟
⎟⎠,

(45)

where we used the assignment i=1 (H2O), i=2 (H2O2), and i=3 (O2). Note that the two zero columns imply that only a single decomposition reaction takes place.

A no-slip boundary condition is imposed at the surface of the particle

u(r)|r=a=0,

(46)

Finally, we set the following boundary conditions at ‘infinity’, the edge of our simulation/computation domain. The mass fractions of the various species are set to their ‘reservoir’ values

ωi(r)|r↑∞=ωi,res,

(47)

where ωi,res is the reservoir mass fraction. This ensures the formation of a stationary state, i.e., there is no depletion of species as there is a continuous influx of fuel and outflow of products at the boundary. For the total fluid, the pressure boundary condition is

p(r)|r↑∞=0,

(48)

and the velocity boundary condition reads

u(r)|r↑∞=v,

(49)

where v is the self-propulsion velocity of the colloid, for which we solve. Note that Eq. (49) is a direct consequence of being in the frame co-moving with the SPP.

Before we move on to the results let us return to an N component model. Solving Eqs. (32-38) simultaneously is a non-trivial task, even for the simple system of a partially catalytic sphere. In literature the slip model is often used to allow for analytic treatment of the self-diffusiophoretic properties of (spherical) SPPs. [41, 55, 63, 50, 42, 25, 53, 64, 51, 49] Here, we show the additional assumptions that are required to arrive at an analytically tractable slip model from our set of multi-component equations.

The slip model is based on an expansion of Eqs. (32-38) in the small layer around the colloid where the molecule-surface interaction potential is non-negligible. Let δ be the ‘range’ of the longest-ranged surface-molecule interaction potential. That is, for |r−s|>δ we have Φk(r)≪kBT, where s is the closest point on the surface to r. Let κ−1(s) be the local curvature, then we can write λ(s)=δκ(s). Note that κ−1(s)=a for a sphere. We further introduce the concept of the local Damköhler number, which gives the ratio of the reaction rate and the diffusive mass transfer rate

Dak(s)=R(s)Dkκ(s),

(50)

where R(s) is the reaction rate. In order to apply the slip-layer approximation, an expansion of Eqs. (32-38) in terms of λ(s) and Da(s), must be accurate to first order. This implies that λ(s)≪1 and Dak(s)≪1, for all points on the surface and for all k. Sharifi-Mood et al.[53] investigated the validity of the first-order approximation and concluded that its range of applicability depends strongly on the nature of the interaction potential and the size of the colloid, however, λ(s)<10−3 would be in the right regime.

When the first-order approximation is valid, we obtain local, instantaneous equilibrium in the small layer around the particle in the direction perpendicular to the surface. For all k the following holds

j∗k(r)⋅^n(s)=0,

(51)

where s is the closest point to the surface of the particle. [42] For convenience, we introduce the notation z=|r−s| for the distance perpendicular to the surface. Then Eq. (51) can be written as

j∗k(z)≡j∗k(s+z^n(s))⋅^n(s)=0.

(52)

Using the assumption of Eq. (52) the ωk(z)≡ωk(s+z^n(s)) are solved within the layer to obtain the species profiles perpendicular to the surface. However, since there are coupling terms in our expressions for j∗k(z), this is a non-trivial problem.

If we assume that ∇m(r)≈0, the equations for the kinematic mass fluxes decouple and the above condition reduces to

j∗k(z)=−Dk(∂∂zωk(z)−ωk(z)Ftot,k(z))=0,

(53)

which can be solved analytically for ωk(z), with Ftot,k(z)≡Ftot,k(s+z^n(s))⋅^n(s). Note that this implies that the molecular masses are roughly equal. The solution can be written as

ωk(z)=ωk(δ)exp(∫zδFtot,k(z′)dz′),

(54)

where it is assumed that Ftot,k(r>δ)=0. We still require

ω1(z)=1−N∑k=2ω(z).

(55)

This solution can be inserted into the momentum balance of the Stokes’ equation normal to the surface, which (up to first order, see Ref. [42]) reads

∂∂zp(z)=−kBTρN∑k=1ωk(z)mk∂∂z~Φk(z),

(56)

where p(z)≡p(s+z^n(s)). We may solve this equation for the pressure decay perpendicular to the surface. Note that we can pull the ~m≡mk out of the sum, since we had already assumed the difference between the molecular masses to be small in establishing Eq. (53). However, by plugging the expression of Eq. (54) into Eq. (56), the system is still not analytically tractable, due to the complex form of Ftot,k(z).

We therefore make the additional assumptions that the diffusion coefficients are similar ~D≡Dk and that ω1(z)≫ωk>1(z), which implies that ω1(z)≈1. This is exactly the dilute-limit approximation. Then we may write

Ftot,k(z)

=

−(∂∂z~Φk(z)−∂∂z~Φ1(z));

(57)

~Φ∗k(z)

≡

~Φk(z)−~Φ1(z);

(58)

⇒ωk(z)

=

ωk(δ)exp(−~Φ∗k(z));

(59)

∂∂zp(z)

=

−kBTρ~m(∂∂z~Φ1(z)+N∑k=2ωk(z)∂∂z~Φ∗k(z));

(60)

p(z)

=

−kBTρ~m(~Φ1(z)

(61)

−kBTρ~m(+N∑k=2ωk(δ)(e−~Φ∗k(z)−1)),

where we introduced the ‘effective’ potential ~Φ∗k(z), which is the potential of species k with respect to that of the solvent. In addition, we assumed that ~Φ∗k(δ)=0 in writing Eq. (59). Equation (61) follows from Eq. (58) by making use of the boundary conditions at infinity. Our expression differs from the ones given by Sharifi-Mood et al.[53], as the surface-solvent potential is not weighted by the ratio of the molecular masses. We consider the expression derived above to be the correct form for the effective potential in the dilute limit.

The solution for the pressure can be plugged into the Stokes equation for motion parallel to the surface. [42] Say that y is the coordinate parallel to the surface and v(y) is the velocity parallel to the surface, then

η∂2∂z2v(y)−∂∂yp(y)=0.

(62)

Now we make use of the fact that only ωk(δ) has a y dependence (to first order) in Eq. (61). The velocity of the fluid at a distance δ away from the surface is therefore given by

vδ(y)=−kBTρη~mN∑k=2(∂∂yωk(δ))∫δ0s(e−~Φ∗k(s)−1)ds,

(63)

where terms of order O(δ2) have been ignored in evaluating the partial integration required to arrive at the above form. Note that the above integral can be extended to infinity, since ~Φ∗k(z)=0 for z>δ.

Combining the above results and translating them back to the three-dimensional (3D) result, we may now write for the boundary conditions at the edge of the slip layer

u(s+δ^n(s))⋅^n(s)

=

0;

(64)

u(s+δ^n(s))⋅^t(s)

=

−N∑k=2ξk(s)

(65)

−N∑k=2×^t(s)⋅∇ωk(s+δ^n(s)),

where the interaction of the fluid molecules and the surface is taken into account by the coupling parameter

ξk(s)=ρkBTη~m∫∞0t(exp(−~Φ∗k(s+t^n(s)))−1)dt.

(66)

In Eq. (65), ^t(s) represents the tangent vectors to the surface at s. N.B. There are two orthogonal tangent vectors at any two-dimensional (2D) surface in a 3D space. It is understood here that ^t(s) is shorthand notation for both of these vectors. This makes Eq. (65) a set of two conditions, rather than a single condition. Further note that the above conditions are dependent on the position on the surface s and the length of the slip layer δ. When δ≪κ−1(s), we may effectively take the limit δ↓0 of Eqs. (64,65) and assume that the slip velocity is imposed at the surface of the colloid.

Figure 1: Sketch of a z-axisymmetric spherical particle, for which the surface is partitioned into stripes. On these stripes the effective molecule-surface interaction potential ~Φ∗k(s,ϕ) experienced by species k is constant in ϕ, the polar angle, but may have a distance dependence r. The azimuthal angle, in which the figure is symmetric, and the x- and y-axis are shown for completion. The stripes are identified by the label ϕi,k, where i refers to the index of the stripe and the subscript k is used since the striped pattern need not coincide for different species. In this example two solute species are assumed; the stripes for species k=2 are delimited using dashed lines and for species k=3 using solid lines.

Finally, let us return to the spherical particles that we consider in this manuscript. N.B. Our spheres are assumed rotationally symmetric in the z-axis. The magnitude of the velocity of a spherical SPP that follows from the slip-model can be written in terms of the following surface integral

v=−14πa2N∑k=2∬Sξk(ϕ)(^t(ϕ)⋅∇ωk(ϕ))(^t(ϕ)⋅^z)ds,

(67)

where S denotes the sphere’s surface and ds is the measure for the surface integration. [42] The parameters ^t(ϕ), ωk(ϕ), and ξk(ϕ) are the z-axisymmetric forms of the respective 3D expressions, which are now functions only of the polar angle ϕ. If we assume that the surface of the particle can be partitioned into i stripes, for which the effective molecule-surface interaction ~Φ∗k(r,ϕ) is piecewise constant, see Fig. 1, then Eq. (67) can be rewritten further. For species k, the i-th integration surface is denoted using Si,k and the relevant slip factor using ξk(ϕi,k) which is given by Eq. (66) evaluated in a point in the region labeled by ϕi,k. The velocity is then given by

v

=

−N∑k=2∑iξk(ϕi,k)

Missing or unrecognized delimiter for \left

where the ξk(ϕi,k) act as multiplicative constants on the respective domains of integration. This is a particularly useful form for the analysis of the effect of the molecule-surface interaction.

We conclude our discussion of the theory of self-diffusiophoretic drive by considering perturbations of the surface-molecule interaction potential for the slip model. A perturbation formalism can be used to quantify the effect of modifying the surface-molecule interaction in an experiment. For instance, by changing the material of the inert surface or adding surfactants to the system that adsorb to part of the SPP’s surface the interaction potential is modified. If the impact of such a modification is small, the effect on the velocity of the colloid can be straightforwardly determined by utilizing the slip-layer solution for the unperturbed potential. Moreover, the perturbation theory can be used to quickly assess the change in speed and direction of an SPP by varying the surface-molecule interaction over a range of parameters. The relevance of this approach will become more clear in the results section.

Let us assume that the interaction ~Φ1(ϕ) for the solvent remains unperturbed, i.e., only the interaction for the solutes may be modified. We further assume that ξk(ϕ) is piecewise constant in ϕ and that a solution to Eqs. (32-38) with boundary conditions, Eqs. (64,65), is given for a specific set of non-zero effective molecule-surface interaction potentials ~Φ∗k(r,ϕi).

For small perturbations of the effective interaction potential ~Φ′k(r,ϕi), where we keep the position of the piecewise constant domains fixed, the change with respect to the unperturbed effective potential ~Φ∗k(r,ϕi) can be expressed using the effective coupling constant

Veff,i,k=∫∞0r(exp(−~Φ′k(r,ϕi))−1)dr∫∞0r(exp(−~Φ∗k(r,ϕi))−1)dr.

(69)

For small Péclet numbers one may safely assume that the effect of advective back coupling in Eq. (34) is negligible, whereby small Péclet number we mean Pe=av/~D≪1 with ~D as before. N.B. Colloidal SPPs are typically in the low Péclet number regime, as can be readily determined by examining the experimental parameters. [23, 24, 25, 26, 27, 28] We refer to Ref. [49] for a detailed discussion of the effects of a finite Péclet number on the self-propulsion. For Pe≪1, the distribution of species around the colloid does not change substantially by the introduction of the perturbation for the speed at which the particle is moving. This allows us to treat the

Ci,k≡∬Si,k(^t(ϕ)⋅∇ωk(ϕ))(^t(ϕ)⋅^z)4πa2ds

(70)

as constants with respect to the perturbation and write the self-propulsion velocity as a simple weighted sum over the interaction prefactors ξk(ϕi) and species distribution terms

v=−∑iN∑k=1ξk(ϕi)Ci,k.

(71)

This leads to the following expression for the approximate self-propulsion velocity for the perturbed potentials in terms of the original solution

v′=−∑iN∑k=1Veff,i,kξk(ϕi)Ci,k,

(72)

where for the Ci,k we can use the values precomputed for the original potentials. Equation (72) is a powerful tool to evaluate the influence of perturbations in the molecule-surface interaction potential for small-Péclet-number SPPs, as we will see in the following.

7.1 System Parameters and Methods

In order to study the effect of the molecule-surface interaction on the self-diffusiophoretic swimming speed of a colloidal SPP, we use the following quantities as a base set. We consider N=3 species for the multi-component model, with i=1 (H2O), i=2 (H2O2), and i=3 (O2). The molecular mass of the components is given by m1=18 u, m2=34 u, and m3=32 u, respectively, with u=1.66⋅10−27 kg the atomic mass unit. The diffusion coefficients are given by D1=2.3⋅10−9 m2 s−1[65], D2=1.2⋅10−9 m2 s−1[66, 67], D3=1.9⋅10−9 m2 s−1[68, 69, 70], respectively, where we used the values for (self-)diffusion in water. The colloid is assumed to be a=0.5⋅10−6 m in radius. We set the density of hydrogen peroxide to 10% b.v. A simple linear catalytic surface reaction is used as in Eq. (44), with the stoichiometric and mass coefficients as in Eq. (45), i.e., 2H2O2→ 2 H2O + O2, where the reaction rate is given by

R(ϕ)=2.0⋅10−5ms−1×{[]ccl0ϕ≤π−α1ϕ>π−α,

(73)

with α the area of the colloid covered by the catalyst, see Fig. 2. The choice for this specific reaction rate is such that it reproduces the hydrogen-peroxide-depletion rates given in Ref. [23] for this H2O2 concentration; we obtain a consumption of 1.0⋅1011 molecules/second per SPP. From the mass density of water 1.0⋅103 kg m−3 and hydrogen peroxide 1.45⋅103 kg m−3 at room temperature, it then follows that for this volume fraction of H2O2, the total density of the fluid can be approximated by ρ=1.045⋅103 kg m−3. Similarly, we can use the formalism of Ref. [71] to compute the dynamic viscosity of the mixture. Using the viscosity of water 1.0⋅10−3 kg m−1 s−1 and hydrogen peroxide 1.245⋅10−3 kg m−1 s−1 at room temperature, we arrive at η=1.017⋅10−3 kg m−1 s−1 for the mixture’s dynamic viscosity. For these parameters we find that ω1,res=0.861, ω2,res=0.139, and ω3,res=0.000.

Figure 2: (color online) Sketch of the level of catalytic surface coating. The figure is rotationally symmetric in the z-axis, as indicated by the arrow, and the catalytic surface (thick red line) coverage is given by the angle α. The rest of the surface ϕ∈[0,π/2−α] is assumed to be inert.

The above choices leave only four free parameters. The catalytic coverage of the colloid α and the three surface-molecule interaction potentials ~Φk(r). Both the full multi-component and the slip-layer model can be solved using, e.g., finite-element methods or kinetic lattice-Boltzmann schemes. In this manuscript, we solve our system of equations using the COMSOL Multiphysics 4.4 solver suite. For the slip model, we combine the slip boundary condition of Eqs. (64-66) with the full multi-component solution in the region where the surface-molecule interaction is zero. We refer to this as the hybrid slip-layer model. Finally, we consider the dilute-limit approximation, where we combine the slip model with Eqs. (64-66), but ignore the m(r)-related term in Eq. (35).

7.2 Comparison between the Full Multi-Component Model and the Dilute-Limit Slip Model

Let us start by assuming that the interaction between the three species of molecules is uniform over the surface. We choose the following form for the surface interaction potential ~Φ3(r) between the SPP’s surface and the solvated oxygen

f(r)

=

1r3−1(a+b)3;

(74)

g(r)

=

f(r)−∂∂rf(a+b)(r−(a+b));

(75)

~Φ3(r)

=

{0r<a∨r≥a+bg(r)/g(a)a≤r<a+b,

(76)

where a is the colloid radius, set to 0.5 μm here, b is 1 nm, and f(r) and g(r) are auxilliary functions. This corresponds to a repulsion that decays as ~Φ3(z)∝1/z3, with z the distance between the surface and the molecule, as before. The potential is set to 1 kBT at contact (z=0) and decays to 0 over a length of 1 nm, in such a way that (∂/∂z)~Φ3(z)=0 at this distance; this decay is accomplished by construction via the auxilliary functions. We assume that the other interaction potentials are zero, i.e., the ones for the water and hydrogen peroxide. The particles thus only interact via soft short-ranged potentials, i.e., their interaction with the wall is that of point particles, which is modified by a slight repulsion in the case of the oxygen. The total force acting on the species, see Eq. (37), can now be straightforwardly calculated.

Note that the above choice constitutes a rather simplified set of interaction potentials. Employing hard or strongly divergent potentials at the surface is not admissible in our model, since this would interfere with the way our incompressibility condition is constructed and the ideal gas behavior of the solutes. Any potential that attracts species sufficiently to the surface of the colloid would cause the solute to be strongly accumulated near the surface, since our theory does not incorporate an excluded-volume term. N.B., a strongly repulsive term would also cause problems, due to the back-coupling in the flux equations by the correction term, see Eq. (37). Moreover, employing such potentials would require a re-examination of the reaction scheme in the flux boundary conditions, since it would no longer be obvious at which distance the molecules can react with the surface. Therefore, it is ill-advised to incorporate anything but the soft-potential parts of the surface-molecule interactions, which is the reason behind the above choice.

It should also be emphasized that solving the full multi-component model for a system, where the interaction length is a factor of 10−3 smaller than the colloid diameter, is technically challenging. However, it is necessary in order to compare the validity of the slip-model results, as the slip-model model is only applicable when there is such a strong separation of length scales. [53] We therefore made a judicious choice in specifying the range of the interaction potential to be 1 nm, rather than ∼2 Å, to keep the model computationally manageable.

Figure 3: (color online) The dependence of the self-propulsion velocity v on the imposed hydrogen peroxide concentration, expressed in a by-volume percentage. The red line shows the self-propulsion velocity that follows from solving the full multi-component model, the blue circles show the result of the hybrid slip-layer approximation, and the green squares show the result of the dilute slip-layer approximation.

We can assess the quality of the slip-model approximation, by solving the full set of multi-component equations and comparing this to the slip-model result for the same interaction parameters. In order to make the comparison, the requirements for the validity of the slip-layer approximation were checked and were found to be satisfied to within reasonable tolerance. In Fig. 3, we consider the dependence of the terminal velocity on the concentration of hydrogen peroxide, while keeping all parameters the same for the three models (multi-component, hybrid, and dilute-limit). Note that the full multi-component model and our hybrid multi-component slip-layer approximation model match over the entire range of densities considered here. This is a surprising result, as it implies that use of the hybrid model is permissible for very high levels of solute concentrations. However, one should be careful not to generalize this result too far, as such correspondence may break down for other parameter sets or less well-behaved interaction potentials. Moreover, the use of too high concentrations will violate the approximations that went into our multi-component model.

Note that the dilute-limit approximation gives a qualitatively similar result, however, the velocity of the SPP is systematically underestimated by a factor of 30%. This can be attributed to the absence of the cross coupling based on m(r) in Eq. (35). Interestingly, the dilute limit result does not approach the solution of the full multi-component model in the limit of zero solute concentration. This implies that the cross-coupling term in Eq. (35) is not sub-dominant even in the dilute limit and could significantly impact the self-propulsion velocity, when there are substantial differences in the mass of the molecular species. The level of correspondence is, however, good enough to justify the use of the equal-mass dilute-limit approximation for a hydrogen-peroxide solution.

For further discussion on the topic of multi-component modeling in the dilute limit, we refer the interested reader to the work by Sharifi-Mood et al.[53] In Ref. [53] the effect of the ratio between the interaction length and the colloid radius on the self-propulsion speed is examined and found to be significant when the radius of the colloid is sub-micron in size.

7.3 Slip-Model Self-Propulsion for Hard Surface-Molecule Interactions

From here on, we consider only the hybrid slip model. The reason for this is, that is it far less computationally expensive to derive results using this approximation and it therefore allows us to cover a larger region of parameter space. Moreover, we can employ less well-behaved interaction potentials, since the slip prefactors ξk(ϕ) can be precomputed analytically and numerical stability in the slip-layer is not a prerequisite. In fact, it is commonplace to relax the conditions on the potentials significantly, in order to permit hard-core interactions. [42] From a technical point, this is questionable, but we decided to follow suit here, in order to relate our result to previous studies. In this section, we still assume that the surface-molecule interaction is homogeneous. The molecule-surface interaction is of the form

~Φk(r,ϕ)={0r−a>ak∞r−a≤ak,

(77)

where the ak is the molecular radius of the k-th component. It should be noted that Eq. (58) is only well-defined for hard interaction potentials, when we specify the way these can be subtracted, again being indicative of the inherent problems of using such potentials. In particular, to avoid infinities in ξk(ϕ), the solvent must have a smaller hard radius than all solutes, i.e., a1<ak>1. We then define the difference between the two infinities to be zero, such that the contribution of Eq. (58) to ξk(ϕ) is zero over the range a<r<a+a1. Note that hard potentials can only be introduced in this manner into the slip-layer approximation, as the more complex expression for the total force of Eq. (37) would result in an ill-posed problem.

In the following we typically use a1=1.4 Å, a2=1.9 Å, and a3=1.8 Å. [72] Figure 4a shows concentration profile of oxygen around the colloid for α=π/2 and the above parameter choices. The maximum concentration of oxygen around the platinum cap is approximately 9 mM, in agreement with the oxygen-production-rate measurements of Ref. [23]. The corresponding fluid velocity profile is shown in Fig. 4b, which gives