Diffusiophoretic SelfPropulsion for Partially Catalytic Spherical Colloids
Abstract
Colloidal spheres with a partial platinum surface coating perform autophoretic motion when suspended in hydrogen peroxide solution. We present a theoretical analysis of the selfpropulsion velocity of these particles using a continuum multicomponent, selfdiffusiophoretic model. With this model as a basis, we show how the sliplayer approximation can be derived and in which limits it holds. First, we consider the differences between the full multicomponent model and the sliplayer approximation. Then the slip model is used to demonstrate and explore the sensitive nature of the particle’s velocity on the details of the moleculesurface 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 selfpropelled particles.
5pt
0.84(0.08,0.93) ©2015 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.
I Introduction
Over the past decade, research into selfpropelled particles (SPPs) has undergone rapid development and attracted strong interest from the scientific community, due to the inherent outofequilibrium nature of their behavior. [1, 2] Examples of naturally occurring selfpropelled ‘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 ( nm – 1 m) and suspended in a liquid medium, since there will be a competition between Brownian (thermal) motion and the selfpropulsion mechanism. [18] Moreover, the viscous nature of the suspending medium strongly impacts the strategies by which autonomous movement may be achieved mechanically. For lowReynoldsnumber 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 selfgenerated 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 gradientbased propulsion mechanisms, including: bimetallic nanorods that move by selfelectrophoresis, [21, 22] Ptcoated colloidal spheres moving by selfdiffusio or selfelectrophoresis, [23, 24, 25, 26, 27, 28] Aucoated colloidal spheres that move by thermophoresis [29, 30] or by thermally induced local mixingdemixing transitions, [31] hollow cones that expel oxygen bubbles, [32, 33] and many others. A common theme for most gradientbased 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 outofequilibrium phenomena [2] and may be utilized as functional components in microfluidic 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 realworld applications.
From the modeling perspective, a lot of attention has gone into the description of the selfpropulsion mechanism of these particles. In particular, it is accepted that selfelectrophoresis plays a dominant role for the selfpropulsion of bimetallic nanorods. [39, 40] The mechanism by which Ptcoated colloidal spheres achieve motion in hydrogenperoxide solution is, however, not as well understood. This has led to heated debate on the type of phoretic motion by which these particles selfpropel. [23, 24] The traditional view is that selfdiffusiophoresis 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 coarsegrained, continuumlevel approach that employs the dilutelimit sliplayer formalism, see, e.g., Refs. [41, 42]. From a simulation/computation perspective, particlebased models are more commonly used. [43, 44, 27, 45, 46, 47] Investigations into the properties of the selfdiffusiophoretic model have considered the influence of a large number of parameters. For instance, the influence of the reaction mechanism on the sizevelocity dependence of the SPPs was matched to experiment. [26, 25] In addition, there have been more fundamental studies which unified the continuum and particlebased descriptions [42, 48] and considered the influence of advective and reactive processes on the selfpropulsion. [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 slipbased selfpropulsion. 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 particlebased simulation models, wherein anisotropy of the solventsurface interaction is often a crucial ingredient to achieve movement. [43, 44, 27, 45, 47] Furthermore, from a physical perspective, the assumption of a homogeneous moleculesurface interaction is questionable.
In this manuscript, we solve a set of multicomponent diffusionadvection equations [57] with appropriate boundary conditions [41, 42] to describe the selfdiffusiophoretic movement of a selfpropelled particle (SPP). We depart from the traditional singlecomponent dilutelimit approximation in order to maintain momentum conservation and incompressibility of the total fluid at high fuel (hydrogenperoxide) 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 selfconsistent model. The role of the solventsurface potential is elucidated and it is shown how the sliplayer approximation follows from these equations. In addition, we consider the differences with previously established multicomponent models. Our advectiondiffusionreaction model is subsequently used to analyze the motion of a catalytically coated colloid suspended in hydrogenperoxide solution. In particular, we investigate the sensitivity of the selfpropulsion velocity on the surfacemolecule 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 surfacemolecule interaction by coupling it to the level of catalytic coating.
We find a strong sensitivity of the SPP’s velocity on the moleculesurface 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 moleculesurface interaction potential for low Péclet numbers. Finally, we briefly consider the sensitivity of the selfpropulsion 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 SPPbased applications and interpretation of experimental results.
The remainder of the document is structured as follows. In Section II we introduce our multicomponent model for bulk fluids. The boundary conditions required to achieve selfpropulsion for arbitrary shapes are discussed in Section III. This general discussion is exemplified on the basis of a spherical Janus SPP suspended in hydrogen peroxide solution in Section IV. In Section V the dilute limit of the multicomponent selfpropulsion model is analyzed and the slip model is recovered. This analysis is built upon in Section VI, where for low Péclet numbers, we introduce a method to determine the influence of small changes in the moleculesurface interaction potential on the selfpropulsion velocity. We subsequently present our results in Section VII, which is partitioned into several subsections. We introduce our system parameters and numerical solution approach in Section VIIA. This is followed by a comparison between the full multicomponent result and the dilutelimit slip model result for the motion of a colloidal SPP in Section VIIB. 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 VIIC. This analysis is extended upon in Section VIID, where we consider the effect of heterogeneities in the surfacemolecule interaction on the SPP’s speed. Our last result is presented in Section VIIE, wherein we briefly consider the influence of variations in the reactivity of the catalytic coating. Finally, we summarize and discuss our results in Section VIII and present an outlook.
Ii The MultiComponent Model for Bulk Fluid
In the following, we consider a three component mixture, consisting of hydrogen peroxide HO (fuel), oxygen O (reactant), and water HO (medium + reactant). However, in this paragraph we will start with a more general formulation of the multicomponent model.
Iia Remarks on MultiComponent Modeling
The choice of a multicomponent 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 hydrogenperoxide 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 HO, [23, 24, 25, 26, 27, 28] make it unclear how the incompressibility constraint on the total fluid is preserved within the diffusionadvection formulation of the behavior of the individual species.
An accurate description of the physics in a multicomponent system uses incompressibility as a constraint that alters the behavior of the individual components. The MaxwellStefan multicomponent formalism is such a model, that describes a total fluid that satisfies the NavierStokes equations, comprised of separate species which satisfy complex diffusionadvection 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’ crossdiffusion coefficients are typically not known and not measurable in experiments. In this section, we therefore present a reduced multicomponent 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 CahnHilliardtype description of the mixingdemixing transition, which needlessly complicates the system for our purposes.
IiB Incompressibility and the Fluid Density Profile
In the following we will consider a multicomponent fluid for which the mass and momentum transport is described using the NavierStokes equations. For any flow characterized by the NavierStokes formalism (both compressible and incompressible), the continuity equation of the fluid is given by
(1) 
where is the mass density, denotes partial differentiation with respect to time, is the spatial differential operator, denotes the dot product (such that is the divergence), is the advective velocity, and the time dependence of the quantities is implicit. The incompressibility condition is given by
(2) 
which in combination with Eq. (1) implies that equivalently
(3) 
This equivalence leads to a subtlety in the definition of incompressibility. Note that and 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., . 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 ; consequently incompressibility guarantees (through ) that the density remains homogeneous. In the derivation below this condition may be relaxed, however, the final result will not be substantially affected.
IiC Chemical Potential and Diffusive Flux
We start our derivation from the expression for the chemical potential [60] of the species that comprise a nonideal, multicomponent mixture
(4) 
where the subscript is used to identify the th species out of total species; is the reference chemical potential of that species; is the Boltzmann constant and is the temperature; is the activity coefficient of species with the position coordinate; is an external potential acting on molecules of species ; and is the molar fraction of species , which is defined as
(5) 
with the particle density of species . The activity coefficient accounts for nonideality and is therefore dependent on the concentrations of the other species as well, which is denoted here by the use of the set . In the following we specify to be constant in ; 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 .
In order to derive the continuity equations for the 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 . We may now introduce the mass fractions of the species as
(6) 
Species has a molecular mass of , which gives us the equivalence and allows us to write
(7) 
if we assume that there is a field such that , with the total number density. Here, is a quantity that has the properties of a numberdensity 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 is position (and time) independent, since the sum over 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
(8) 
This results in the following expression for the gradient (denoted by ) of the molar fraction
(9) 
Using the above quantities and Eq. (9), it follows that the expression for the number density flux of species , in the frame of reference that is comoving with the total fluid, can be written as
(10)  
(12)  
where is the mobility of the th species, which is related to that species’ diffusion constant via the EinsteinSmoluchowski relation ; is the force per particle acting on the component; and . Here, we used the first order expansion for small deviations of from , in order to write in Eq. (12). [59]
IiD Mass Flux and Continuity
The expression in Eq. (12) allows us to write the kinematic mass flux (in the comoving frame) as
Note that the kinematic mass flux has the physical dimension of velocity. Further note that when all species have the same molecular mass, the masscorrection (second) term drops out of the above equations, which allows us to recover the traditional Fickian diffusion.
The kinematic mass flux in the comoving frame and the kinematic mass flux in the laboratory frame are related according to the following transformation
(14) 
where is the velocity of the total fluid, which accounts for the advective contribution. The timedependent continuity equations for a system, in which there are no bulk reactions, read
(15) 
Note that bulk reaction terms may be straightforwardly accounted for in our formalism by introducing source and sink terms into Eq. (15). The timedependent continuity equation for the th species in the comoving frame now follows from the relation in Eq. (14)
(16) 
The above set of equations corresponds to the simplified diffusionadvection model put forward in Ref. [58].
IiE Momentum Transport
For momentum transport in the total fluid, we are interested in the low Reynolds number limit. That is, , where , with the relevant length scale of the problem, the speed of the object or fluid, and the dynamic viscosity of the fluid. For colloidal SPPs typically used in experiments, is satisfied. [23, 24, 25, 26, 27, 28] This implies that we can work with the linearized form of the NavierStokes equations for the total fluid
(17)  
(18) 
where is the pressure, the tensorial Laplacian, and the force per volume acting on the fluid. Equation (18) gives the incompressibility constraint for the total fluid. The force term derives from the particle fluxes of the individual species [59]
(19) 
Using Eq. (LABEL:eq:fluxsp) it now follows that
(20)  
The first two terms exactly cancel each other for an incompressible fluid. This is a nice feature of the multicomponent model that takes into account the molecular mass of the species and utilizes incompressibility, which is missing in the traditional dilutelimit approximation.
IiF 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 selfconsistent in the incompressibility assumption. From the properties of it follows that a sufficient condition to obtain incompressibility of the total fluid from the individual continuity equations, see Eq. (16), is given by
(21) 
This is also physically reasonable as otherwise there would be mass flow with respect to the frame comoving 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 to ensure that this constraint is satisfied. This is typically the point where crossdiffusion terms are introduced in order to recover incompressibility. [59, 58, 57] However, the introduction of crossdiffusion 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 forcefree part of the flux. We may write this flux as
(22) 
If we impose
(23) 
then we obtain
(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 , 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 comoving frame, some correction force must be added to the gradients of the potentials. We make the following ansatz for the total force acting on particles of the th species in the comoving frame
(25) 
By plugging in 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
(26) 
We may thus write for the comoving kinematic mass flux
(27)  
(28)  
Note that Eq. (28) now follows from the previously imposed condition of Eq. (23) and the properties of .
IiG Further Properties and Reductions
The above flux model has the desirable property that throughout space, only the continuity equations for 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 comoving frame. The force acting on the fluid, see Eq. (20), should not be modified. This becomes clear by considering the case wherein all and 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
(30) 
for the fluid and
(31) 
for the diffusive species, respectively.
IiH Overview of the Model
Summarizing, the multicomponent model for the bulk fluid satisfies the following equations under the assumption of a homogeneous incompressible medium
(32)  
(33)  
(34)  
(35)  
(36)  
(37)  
(38) 
This set of equations is fully consistent with the incompressibility assumption for the total fluid and dependent on a minimal number of fields: , , , and ; and parameters: , , , , and .
Note that this model presents a significant departure from the typical dilutelimit approximation, since the difference in mass of the molecular species is taken into account and incompressibility is incorporated in a selfconsistent 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 MaxwellStefan formulation for higher solute concentrations.
Iii SelfDiffusiophoretic Motion of Objects of Arbitrary Shape
Thus far, we have only concentrated on the description of the ‘bulk’ fluid, not the mechanism by which selfpropulsion is achieved. To obtain selfdiffusiophoretic motion of a particle in our multicomponent description we require only a set of boundary conditions. We are interested in the stationary state of the system and we therefore solve the timeindependent variants of Eqs. (3238). By solving these equations simultaneously, the selfpropulsion velocity and induced fluid flow profile in the stationary state may be obtained.
At the surface of the particle we require a noslip boundary condition. That is,
(39) 
where 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 multicomponent system there can be a large number of (nonlinear) 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
(40) 
where is the normal vector to the surface at the point . If there were a net mass flux into the surface, the noslip boundary condition of Eq. (39) would be violated, due to the relation between and . 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 . That is, in the absence of an interaction between the surface of the SPP and the solvent/solutes, there is no backcoupling 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].

The existence of a inertial transformation to a comoving 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 nonzero selfpropulsion velocity may be achieved, i.e., fuel must not be depleted.
Iv A SelfDiffusiophoresing Sphere in Hydrogen Peroxide
In this section, we concentrate exclusively on a single spherical SPP of radius , that is suspended in a threecomponent mixture of water, hydrogen peroxide, and oxygen (), to better illustrate our autophoretic 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 the distance measured to the origin, the azimuthal angle, and the polar angle. The system is rotationally symmetric around the 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 and the polarangle tangent vector by ; the Cartesian unit vectors are given by , , and .
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 HO 2 HO O for hydrogen peroxide decomposition at the platinum surface. The reaction kinetics satisfy
(41)  
(42)  
(43) 
where is the reaction rate. The set of equations, Eqs. (4143), is used to treat bulk catalytic decomposition and needs to be converted to a surface model, considering the specific geometrical properties of the thin Ptcoating 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 multiplestep reaction mechanism through a variety of intermediate Ptcomplexes is a more probable route for the decomposition. [61, 62] However, there is some indication that the ratelimiting step in the reaction process is linear in the hydrogenperoxide concentration. [40] Taking the above considerations into account, we may write for the flux boundary conditions
(44) 
where the are stoichiometric/mass coefficients and is the surface reaction rate. The relation to the bulk reaction rate is given by , where is the bulk reaction rate and 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
(45) 
where we used the assignment (HO), (HO), and (O). Note that the two zero columns imply that only a single decomposition reaction takes place.
A noslip boundary condition is imposed at the surface of the particle
(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
(47) 
where 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
(48) 
and the velocity boundary condition reads
(49) 
where is the selfpropulsion velocity of the colloid, for which we solve. Note that Eq. (49) is a direct consequence of being in the frame comoving with the SPP.
V The Dilute Limit Slip Model for SelfDiffusiophoresis
Before we move on to the results let us return to an component model. Solving Eqs. (3238) simultaneously is a nontrivial 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 selfdiffusiophoretic 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 multicomponent equations.
The slip model is based on an expansion of Eqs. (3238) in the small layer around the colloid where the moleculesurface interaction potential is nonnegligible. Let be the ‘range’ of the longestranged surfacemolecule interaction potential. That is, for we have , where is the closest point on the surface to . Let be the local curvature, then we can write . Note that 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
(50) 
where is the reaction rate. In order to apply the sliplayer approximation, an expansion of Eqs. (3238) in terms of and , must be accurate to first order. This implies that and , for all points on the surface and for all . SharifiMood et al. [53] investigated the validity of the firstorder approximation and concluded that its range of applicability depends strongly on the nature of the interaction potential and the size of the colloid, however, would be in the right regime.
When the firstorder approximation is valid, we obtain local, instantaneous equilibrium in the small layer around the particle in the direction perpendicular to the surface. For all the following holds
(51) 
where is the closest point to the surface of the particle. [42] For convenience, we introduce the notation for the distance perpendicular to the surface. Then Eq. (51) can be written as
(52) 
Using the assumption of Eq. (52) the are solved within the layer to obtain the species profiles perpendicular to the surface. However, since there are coupling terms in our expressions for , this is a nontrivial problem.
If we assume that , the equations for the kinematic mass fluxes decouple and the above condition reduces to
(53) 
which can be solved analytically for , with . Note that this implies that the molecular masses are roughly equal. The solution can be written as
(54) 
where it is assumed that . We still require
(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
(56) 
where . We may solve this equation for the pressure decay perpendicular to the surface. Note that we can pull the 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 .
We therefore make the additional assumptions that the diffusion coefficients are similar and that , which implies that . This is exactly the dilutelimit approximation. Then we may write
(57)  
(58)  
(59)  
(60)  
(61)  
where we introduced the ‘effective’ potential , which is the potential of species with respect to that of the solvent. In addition, we assumed that 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 SharifiMood et al. [53], as the surfacesolvent 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 is the coordinate parallel to the surface and is the velocity parallel to the surface, then
(62) 
Now we make use of the fact that only has a dependence (to first order) in Eq. (61). The velocity of the fluid at a distance away from the surface is therefore given by
(63) 
where terms of order 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 for .
Combining the above results and translating them back to the threedimensional (3D) result, we may now write for the boundary conditions at the edge of the slip layer
(64)  
(65)  
where the interaction of the fluid molecules and the surface is taken into account by the coupling parameter
(66) 
In Eq. (65), represents the tangent vectors to the surface at . N.B. There are two orthogonal tangent vectors at any twodimensional (2D) surface in a 3D space. It is understood here that 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 and the length of the slip layer . When , we may effectively take the limit of Eqs. (64,65) and assume that the slip velocity is imposed at the surface of the colloid.
Finally, let us return to the spherical particles that we consider in this manuscript. N.B. Our spheres are assumed rotationally symmetric in the axis. The magnitude of the velocity of a spherical SPP that follows from the slipmodel can be written in terms of the following surface integral
(67) 
where denotes the sphere’s surface and is the measure for the surface integration. [42] The parameters , , and are the 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 stripes, for which the effective moleculesurface interaction is piecewise constant, see Fig. 1, then Eq. (67) can be rewritten further. For species , the th integration surface is denoted using and the relevant slip factor using which is given by Eq. (66) evaluated in a point in the region labeled by . The velocity is then given by
where the act as multiplicative constants on the respective domains of integration. This is a particularly useful form for the analysis of the effect of the moleculesurface interaction.
Vi The Slip Model and Perturbations of the SurfaceMolecule Potential
We conclude our discussion of the theory of selfdiffusiophoretic drive by considering perturbations of the surfacemolecule interaction potential for the slip model. A perturbation formalism can be used to quantify the effect of modifying the surfacemolecule 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 sliplayer 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 surfacemolecule 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 for the solvent remains unperturbed, i.e., only the interaction for the solutes may be modified. We further assume that is piecewise constant in and that a solution to Eqs. (3238) with boundary conditions, Eqs. (64,65), is given for a specific set of nonzero effective moleculesurface interaction potentials .
For small perturbations of the effective interaction potential , where we keep the position of the piecewise constant domains fixed, the change with respect to the unperturbed effective potential can be expressed using the effective coupling constant
(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 with 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 selfpropulsion. For , 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
(70) 
as constants with respect to the perturbation and write the selfpropulsion velocity as a simple weighted sum over the interaction prefactors and species distribution terms
(71) 
This leads to the following expression for the approximate selfpropulsion velocity for the perturbed potentials in terms of the original solution
(72) 
where for the we can use the values precomputed for the original potentials. Equation (72) is a powerful tool to evaluate the influence of perturbations in the moleculesurface interaction potential for smallPécletnumber SPPs, as we will see in the following.
Vii Results
Viia System Parameters and Methods
In order to study the effect of the moleculesurface interaction on the selfdiffusiophoretic swimming speed of a colloidal SPP, we use the following quantities as a base set. We consider species for the multicomponent model, with (HO), (HO), and (O). The molecular mass of the components is given by u, u, and u, respectively, with kg the atomic mass unit. The diffusion coefficients are given by m s [65], m s [66, 67], m s [68, 69, 70], respectively, where we used the values for (self)diffusion in water. The colloid is assumed to be 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., 2HO 2 HO O, where the reaction rate is given by
(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 hydrogenperoxidedepletion rates given in Ref. [23] for this HO concentration; we obtain a consumption of molecules/second per SPP. From the mass density of water kg m and hydrogen peroxide kg m at room temperature, it then follows that for this volume fraction of HO, the total density of the fluid can be approximated by kg m. Similarly, we can use the formalism of Ref. [71] to compute the dynamic viscosity of the mixture. Using the viscosity of water kg m s and hydrogen peroxide kg m s at room temperature, we arrive at kg m s for the mixture’s dynamic viscosity. For these parameters we find that , , and .
The above choices leave only four free parameters. The catalytic coverage of the colloid and the three surfacemolecule interaction potentials . Both the full multicomponent and the sliplayer model can be solved using, e.g., finiteelement methods or kinetic latticeBoltzmann 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. (6466) with the full multicomponent solution in the region where the surfacemolecule interaction is zero. We refer to this as the hybrid sliplayer model. Finally, we consider the dilutelimit approximation, where we combine the slip model with Eqs. (6466), but ignore the related term in Eq. (35).
ViiB Comparison between the Full MultiComponent Model and the DiluteLimit 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 between the SPP’s surface and the solvated oxygen
(74)  
(75)  
(76) 
where is the colloid radius, set to 0.5 m here, is 1 nm, and and are auxilliary functions. This corresponds to a repulsion that decays as , with the distance between the surface and the molecule, as before. The potential is set to 1 at contact () and decays to over a length of 1 nm, in such a way that 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 shortranged 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 excludedvolume term. N.B., a strongly repulsive term would also cause problems, due to the backcoupling in the flux equations by the correction term, see Eq. (37). Moreover, employing such potentials would require a reexamination 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 illadvised to incorporate anything but the softpotential parts of the surfacemolecule interactions, which is the reason behind the above choice.
It should also be emphasized that solving the full multicomponent model for a system, where the interaction length is a factor of smaller than the colloid diameter, is technically challenging. However, it is necessary in order to compare the validity of the slipmodel results, as the slipmodel 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 nm, rather than Å, to keep the model computationally manageable.
We can assess the quality of the slipmodel approximation, by solving the full set of multicomponent equations and comparing this to the slipmodel result for the same interaction parameters. In order to make the comparison, the requirements for the validity of the sliplayer 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 (multicomponent, hybrid, and dilutelimit). Note that the full multicomponent model and our hybrid multicomponent sliplayer 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 wellbehaved interaction potentials. Moreover, the use of too high concentrations will violate the approximations that went into our multicomponent model.
Note that the dilutelimit 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 in Eq. (35). Interestingly, the dilute limit result does not approach the solution of the full multicomponent model in the limit of zero solute concentration. This implies that the crosscoupling term in Eq. (35) is not subdominant even in the dilute limit and could significantly impact the selfpropulsion 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 equalmass dilutelimit approximation for a hydrogenperoxide solution.
For further discussion on the topic of multicomponent modeling in the dilute limit, we refer the interested reader to the work by SharifiMood et al. [53] In Ref. [53] the effect of the ratio between the interaction length and the colloid radius on the selfpropulsion speed is examined and found to be significant when the radius of the colloid is submicron in size.
ViiC SlipModel SelfPropulsion for Hard SurfaceMolecule 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 wellbehaved interaction potentials, since the slip prefactors can be precomputed analytically and numerical stability in the sliplayer is not a prerequisite. In fact, it is commonplace to relax the conditions on the potentials significantly, in order to permit hardcore 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 surfacemolecule interaction is homogeneous. The moleculesurface interaction is of the form
(77) 
where the is the molecular radius of the th component. It should be noted that Eq. (58) is only welldefined 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 , the solvent must have a smaller hard radius than all solutes, i.e., . We then define the difference between the two infinities to be zero, such that the contribution of Eq. (58) to is zero over the range . Note that hard potentials can only be introduced in this manner into the sliplayer approximation, as the more complex expression for the total force of Eq. (37) would result in an illposed problem.
In the following we typically use Å, Å, and Å. [72] Figure 4a shows concentration profile of oxygen around the colloid for and the above parameter choices. The maximum concentration of oxygen around the platinum cap is approximately mM, in agreement with the oxygenproductionrate measurements of Ref. [23]. The corresponding fluid velocity profile is shown in Fig. 4b, which gives in the laboratory frame, wherein the fluid is stationary and the particle moves. The velocity achieved for this system is m s. That is, the particle is moving in the direction of the platinum cap with a speed of m s, in agreement with the estimate given in Ref. [23]. The calculated speed is substantially smaller than is typically observed in experiment for these parameters [23, 24, 25, 26, 27, 28], in particular Ref. [23] establishes the velocity to be around m s. Moreover, in the experiment the particle is moving in the direction of the uncoated surface. Unfortunately, it is difficult to extract information about the nature of the selfpropulsion mechanism from the mismatch between this diffusiophoretic result and the experimental observations, as will become clear in the following.
The flow field in Fig. 4b displays the expected potential farfield contribution for forcefree SPPs with finite extent. [73] Here, the term potential signifies a dipole with finite extent. However, a combination of a dipolar and higherorder modes is found close to the surface. This result is similar to the thincap limit for selfthermophoresing particles studied in Refs. [74, 29]. We found that the potentialswimmer approximation only holds for distances . Figure 4c shows the flow field in the comoving frame. For the reasonable parameter choices made here, the dominant contribution to the selfpropulsion velocity comes from a small region near the divide between the inert and catalytic material, which is only nm in size, see Fig. 4c. This can be attributed to the fact that the concentration gradient along the surface is the largest in this region, as is illustrated in Fig. 4d,e, which show the mass fraction and its gradient along the surface, respectively.
Note that is extremely spiked in this region, as is further evidenced by the inset to Fig. 4e where a closeup of is shown. From our numerical results it is clear that this feature remains present at small length scales. However, it is not immediately evident whether this feature is a cusp or not, due to the complexity of the coupled system of equations that is solved simultaneously. Finally, it should be remarked that in applying the sliplayer approximation, the gradient is evaluated at a distance away from the surface, whereas Fig. 4e shows the result at the surface (). Without advection, crosscoupling, and forcecorrection, the differential equation for the solute is of the Harmonic form (). For this equation the solution can be expressed in terms of a Legendrepolynomial expansion with radial dependence according to , with the summation index. Therefore, the convergence of the expansion at is improved by the geometric factors