Ionic screening and dissociation are crucial for understanding chemical self-propulsion in polar solvents

Polar solvents like water support the bulk dissociation of themselves and their solutes into ions, and the re-association of these ions into neutral molecules in a dynamic equilibrium, e.g., H2O2 " H + + HO2 . Using continuum theory, we study the influence of these association–dissociation reactions on the self-propulsion of colloids driven by surface chemical reactions (chemical swimmers). We find that association–dissociation reactions should have a strong influence on swimmers’ behaviour, and therefore should be included in future modelling. In particular, such bulk reactions should permit charged swimmers to propel electrophoretically even if all species involved in the surface reactions are neutral. The bulk reactions also significantly modify the predicted speed of chemical swimmers propelled by ionic currents, by up to an order of magnitude. For swimmers whose surface reactions produce both anions and cations (ionic self-diffusiophoresis), the bulk reactions produce an additional reactive screening length, analogous to the Debye length in electrostatics. This in turn leads to an inverse relationship between swimmer radius and swimming speed, which could provide an alternative explanation for recent experimental observations on Pt-polystyrene Janus swimmers [S. Ebbens et al., Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 020401]. We also use our continuum theory to investigate the effect of the Debye screening length itself, going beyond the infinitely-thin-screening-length approximation used by previous analytical theories. We identify significant departures from this limiting behaviour for micron-sized swimmers under typical experimental conditions and find that the approximation fails entirely for nanoscale swimmers.


Introduction
The 20th century witnessed a revolution in condensed matter physics, due to the ready availability of well-characterised colloidal particles (1 nm to 10 mm in size). These particles are often viewed as 'large atoms': they are small enough to be subject to Brownian motion, and thus to all the machinery of equilibrium statistical physics, but large enough that their microscopic dynamics and interactions can be observed and tuned. Studying colloidal particles has led to fundamental breakthroughs. Most notably, the observation and subsequent understanding of Brownian motion in colloidal systems 1 led to acceptance of the molecular picture of matter.
Moving into the 21st century, physicists have started to recruit colloids to tackle systems that are intrinsically out-of-equilibrium, specifically where the components are themselves self-propelled. This is the field of 'active matter'. A wide range of novel, selfpropelled colloids [2][3][4][5][6][7][8][9] have been synthesised-see example sketches in Fig. 1 for two designs relevant to this work. Such self-propelled Fig. 1 Cartoon of the two paradigmatic chemical swimmers discussed in the text. Both swimmers move at a few mm s À1 in 10% H 2 O 2 solution, powered by the decomposition of H 2 O 2 on their surfaces. (a) Bimetallic (typically goldplatinum) rod, 2 of typical length 2 mm and width 300 nm. The accepted propulsion mechanism for these swimmers is via a H + current, as shown. (b) Platinum-polymer (usually polystyrene) Janus sphere, 3 of typical radius 1 mm. colloids are intrinsically out of equilibrium-they continuously transform chemical, thermal, or electromagnetic energy into directed motion-and recent work has focussed on using these systems to experimentally explore exciting non-equilibrium phenomena such as phase separation and collective motion. [10][11][12][13] In parallel with this research, much work has gone into understanding the experimental propulsion mechanisms at the level of surface chemical reactions. [2][3][4][14][15][16][17] Working out how these tiny motors function is a fundamental problem in its own right. However, understanding the propulsion mechanism is also an essential first step in understanding the experimental collective behaviour. This is because unlike biological swimmers such as E. coli, where the biochemical reactions responsible for propulsion take place internally, synthetic swimmers are usually propelled by surface chemical reactions that directly modify the chemical, electrostatic, or temperature fields of their surroundings. These fields modify the propulsion speed of other swimmers, generating so-called 'phoretic' swimmer-swimmer interactions in addition to the hydrodynamic and contact interactions experienced by all swimmers. [18][19][20][21] As these phoretic interactions are directly coupled to the chemical reactions responsible for propulsion, knowing how artificial swimmers self-propel is essential for understanding their collective behaviour.
This bottom-up approach contrasts with the tactic employed by most theoretical modellers of active matter, which is to explore the phenomenology arising from minimal or effective models of swimmer-swimmer interactions, by which we mean models that have no explicit connection with microscopic parameters. Minimal models include the Vicsek model 22 or the active Brownian model 23,24 which considers only contact forces, while more complex, but still effective, hydrodynamic [25][26][27] and phoretic models 18,[28][29][30][31] are now also being explored. This research is valuable in unlocking generic non-equilibrium physics principles and is often able to reproduce experimental behaviour (phase separation, etc.) surprisingly well. However, it is often unclear whether such agreement is due to judicious tuning of parameters, or whether there is a real correspondence in generic physics so that details do not matter. In lieu of a general theory for nonequilibrium physics, only in-depth knowledge of the microscopic physics for specific swimmers is likely to allow the resolution of this issue. If a microscopically justified model predicts the same phenomenon as generic models, then details may indeed not matter. On the other hand, if a microscopic model fails, then new physics is indicated.
That said, current understanding of self-propulsion mechanisms is often very incomplete, and hence the necessary foundations for models which can reproduce multiparticle behaviour from a bottom-up perspective are lacking. In particular, most research so far has focussed on unravelling the surface chemistry of the swimmer itself. 2,[14][15][16][17][32][33][34] This is understandable and necessary, but it has meant that other aspects, such as the chemistry of the bulk solvent, have been neglected.
In this article, we show that taking into account the chemistry of common polar solvents such as water and hydrogen peroxide, significantly, and often qualitatively, modifies the predicted propulsion behaviour of almost all self-propelled synthetic swimmers. In this first study, we limit our discussion to the propulsion of single particles, because this is a necessary first step in understanding more complex behaviour. For most of the paper, we also do not go into the details of the surface chemistry, which are not well understood. This allows us to highlight the bulk solvent effects, at the expense of explicit predictions for the propulsion speed. However, in Section 7, in order to compare our results with experiments, we do examine the predicted propulsion speed obtained with suitably simple assumptions for the surface chemistry. One of the main conclusions of our paper is that more detailed experimental studies of both the surface and bulk chemistry are crucial for a detailed understanding of self-propulsion.

Self-electrophoresis
We discuss here the most experimentally typical self-propelled colloids, which we term 'chemical swimmers'. They are most easily defined by example. Fig. 1 shows two chemical swimmers, both powered by the catalytic decomposition of hydrogen peroxide on their surfaces. Because the colloid surface is anisotropic, this reaction produces chemical gradients which, via interaction with the particle surface, eventually lead to self-propulsion.
We say 'eventually' because the propulsion mechanism of these swimmers is somewhat involved. For the example given in Fig. 1a, H 2 O 2 decomposition does not occur just by the simple chemical reaction but also occurs partially electrochemically, with two half reactions taking place preferentially on the Au or Pt surfaces 2,4 for example These half reactions produce a proton gradient outside the colloid, which generates a local electric field. The colloid surface, like most surfaces in water, is charged, so this electric field causes electroosmotic flow over the colloid surface, leading to self-propulsion. This propulsion mechanism is called 'selfelectrophoresis'. 4 The electric field also generates a proton current outside, which is balanced by an electric current inside the conductive swimmer. A large body of experimental evidence confirms that selfelectrophoresis is the appropriate propulsion mechanism for these bimetallic swimmers. For example, their propulsion speed scales inversely with salt concentration, 11,16 which is expected from a simple application of Ohm's law. Recent results 16,17 indicate that self-electrophoresis is also the appropriate propulsion mechanism for the type of colloid shown in Fig. 1b, which has a single metallic coating. This is at first surprising because there is no obvious mechanism for producing the ionic gradient needed for self-electrophoresis. However, geometrical differences between the equator and pole of the catalytic coating, such as thickness variation, may couple to the half-reaction rates in (R2) and so provide the necessary asymmetry. 17 In this paper, we go further and show that these effects are not limited to swimmers that can support ionic currents. All swimmers in aqueous solution are likely to be self-electrophoretic to a major degree, whatever their surface reaction mechanism.

Surface chemistry
Before we discuss the effects of the bulk solution, we point out one general difficulty with self-electrophoresis that will also apply to other complex propulsion mechanisms. This is that the relevant surface reaction rates are extremely hard to measure. The overall reaction rate (R1) can be easily obtained by measuring reactant or product concentrations, 16,17 but for self-electrophoresis the important rate is the proton production rate (R2), and this is likely to make up only a tiny proportion of the overall reaction, with the remainder proceeding via (R1). 2 Measuring the rate of an individual reaction pathway like (R2) is challenging, and has not yet been done, to our knowledge, for any selfpropelled particle.
This would not be a problem if we could predict these rates. However, surface catalysis is a notoriously sensitive phenomenon, and these surface reaction rates are likely to vary unpredictably with almost every parameter, e.g., pH, ionic strength, and surface roughness. 35 Reversing the argument, the only currently available method of estimating these reaction rates is from the particle propulsion velocity itself. That is, with a sufficiently accurate microscopic model, the surface reaction rates can be inferred from the propulsion speed. 32 The catalytic chemistry of micro-and nanoparticles is of huge industrial importance, so this provides another major motivation for obtaining a detailed theoretical understanding of self-electrophoresis.

Bulk chemistry
At first glance, the chemistry of the bulk solution is much simpler than that of the surface. However, a polar solvent such as water presents two complications which have not yet been fully taken into account. The first of these is electrostatic screening. In so-called phoretic mechanisms, such as self-electrophoresis, fluid flow is generated in a layer around the particle surface. In self-electrophoresis, the thickness of this interaction layer is given by the electrostatic screening or Debye length k À1 : outside this screening layer, the free charge density, which is responsible for fluid flow, decays rapidly to zero.
Analytical studies typically make use of a thin-screening approximation ka c 1, where a is the swimmer radius, because this dramatically simplifies the calculation of propulsion speed, flow fields, etc. 16,17,33,36 However, this assumption is not generally valid: k À1 is of order 100 nm for an experimentally typical 3 M H 2 O 2 solution, 16 and active colloids typically range in size from 10 nm 37 to 10 mm. 15 In this paper, we show that the thin-screening approximation can be dropped from analytical calculations, and that this dramatically reduces the predicted propulsion speeds, by up to several orders of magnitude for nanoscale swimmers. To the best of our knowledge, this result has not been shown even with numerics, such as the finite element method (FEM), for which the thin screening approximation is not employed. 14,34 Here, we use FEM calculations to verify the analytical results.
The second complication of aqueous and similar polar environments is that the solvent is not chemically inert. It is an 'active fluid' that can be driven out of chemical equilibrium by the reactions on the particle surface. This consideration has been appreciated for biological fluids such as the cytoplasm, 20 where biomolecules are continuously synthesized and broken down, but it is also true for simple fluids like water which permit the ionic dissociation of both themselves and any polar solutes, e.g., H 2 O " H + + OH À and H 2 O 2 " H + + HO 2 À . The implications of these reactions for self-electrophoresis are the main focus of this paper. The most striking implication is that a gradient of a neutral molecule like H 2 O 2 will result in ionic gradients, here of H + and HO 2 À ions, which will themselves produce electric fields. This means that a surface reaction with only uncharged species like H 2 O 2 or H 2 O can itself generate self-electrophoretic propulsion. It is worth highlighting that these bulk reactions should also qualitatively modify interparticle interactions. To demonstrate this, we describe an effect called 'reactive screening', 20 which will underlie much of our later discussion. We illustrate this effect with a simple 1D model, see Fig. 2. Let an uncharged molecule of diffusivity D be produced uniformly at a plane surface z = 0 and consumed in the bulk (z 4 0) with rate g. The steady-state concentration profile c(z) then obeys D @ 2 c @z 2 ¼ gc; where diffusion, on the left, is balanced by consumption, on the right. The solution to eqn (1) is an exponentially screened concentration profile, c = c 0 exp(Àqz), with c 0 the concentration at the surface, and the reactive screening length q À1 = (D/g) 1/2 . This uncharged model has been applied to the diffusiophoresis of small particles inside a biological cell, where the relevant bulk reactions are the breakdown of biomolecules in the cytoplasm. 20 With self-electrophoresis, as we shall see, reactive screening Fig. 2 Schematic representation of the simple 1D molecular screening model introduced in the text. A wall releases molecules (red triangles) which diffuse with D (black arrow) and are consumed in the bulk with a rate g (blue symbol). This leads to an exponential decay of the concentration, as indicated using the dashed red line and continuum-level red gradient.
can also exponentially screen the electrostatic potential, effectively turning off the long-range electrostatic interactions, which would otherwise be inevitable. Hence, this reactive screening is a qualitatively new effect of bulk reactions, which cannot be ignored a priori even at the level of phenomenological models of collective motion.

Overview of main results
The theory of self-electrophoretic propulsion is mathematically involved, even without the introduction of additional bulk reactions, so we will use this section to sketch out our main results in advance. This will necessarily skim over or simplify many relevant details that will be addressed fully in later sections.

Overall framework
Our main mathematical result is that the self-electrophoretic propulsion speed of an arbitrary, uniformly charged, spherical swimmer can be written, if a suitable linearization is applied, in the form where U SM is the 'standard model' propulsion speed assuming the thin screening limit without bulk reactions, explored for example in ref. 33 and 38. U SM depends on, among other parameters, the surface reaction rates j s , the salt concentration c salt , and the surface charge density s. Here we introduce two new factors, F and B, to account for realistic electrostatic screening and bulk reactions, respectively. These factors depend on the dimensionless parameters ka and qa, with k and q the inverse electrostatic and inverse reactive screening lengths, and a the swimmer radius.

Electrostatic screening
The new dimensionless factor F(ka) is exactly analogous to the well-known function f (ka) (Henry's function 39 ) that controls the speed of a particle undergoing electrophoresis in an external field via U ext = m E E N , with m E = zef(ka)/Z the electrophoretic mobility, e the dielectric constant, z the particle's surface potential, Z the solution viscosity, and E N the external electric field. 39,40 Both f and F are plotted in Fig. 3. For ka { 1, i.e., for small particles or low salt concentration, F decreases rapidly, scaling as (ka) 3 . This is different from external electrophoresis because of the different geometries of the driving fields-a uniform field for external electrophoresis compared with a dipole for self-electrophoresis. The implication of this a 3 scaling is that, other things being equal, nanoswimmers should swim much slower than microswimmers. Experimentally, however, nanoswimmers are found to swim faster than equivalent microswimmers. 37 From this we conclude that other things are not equal: either the surface reaction rates are much larger for nanoswimmers, or the standard self-electrophoresis theory does not apply for these small swimmers. If this issue can be resolved, which we do not attempt here, it will likely also give insight into the related phenomenon of directed motion in nanoscale biological enzymes. 41,42

Bulk reactions
The effects of ionic dissociation depend upon the nature of the surface reaction responsible for propulsion. A common feature is the importance of the reactive screening length q À1 which controls the propulsion behaviour through the parameter qa. We can understand why qa is the relevant parameter as follows: for qa { 1, the swimmer is smaller than the reactive screening length, so any molecules produced at the swimmer surface will diffuse away or return to the swimmer surface before they have time to react. In this 'reactionless limit', the swimmer will behave as though there are no bulk reactions, which is the usual, tacit assumption. For swimmers larger than the reactive screening length, qa c 1, we are in a 'reactive limit' where the bulk ionic reactions dominate the behaviour. For typical experimental conditions, e.g., 3 mol L À1 H 2 O 2 , we find a reactive screening length q À1 E 70 nm, which is in the centre of the experimental range of swimmer radii. 15,37 Both the reactionless and reactive limits, and the intermediate regime (qa E 1), should therefore be experimentally relevant.
We now explain the effect of bulk reactions on specific types of swimmer. The overall surface reaction we focus on is the H 2 O 2 decomposition reaction (R1). As we have seen, this overall reaction can occur through several different pathways. We therefore define three model swimmers, shown in Fig. 4 with surface reactions that are representative of these different pathways. A real swimmer might exhibit any or all of these.
The upper panels (a-c) of Fig. 4   Second, we assume that the chemical fluxes are dipolar rather than monopolar, i.e., they both leave and enter the particle surface. For a uniformly charged particle, which we assume for simplicity, only the dipolar component of the flux contributes to the propulsion speed, 33 so our choice of a dipolar surface-flux profile does not affect our results, and simplifies the argument. Third, we consider only self-electrophoresis here, and ignore 'neutral self-diffusiophoresis', which is propulsion generated by a direct non-electrostatic interaction between a neutral species, such as O 2 and the swimmer surface, 3,43 and which is typically much weaker than self-electrophoresis. 16,34,44 Hence the model swimmer in (a) does not move, because, without bulk reactions, a surface reaction involving only uncharged species cannot generate electric fields, and therefore cannot produce self-electrophoresis.
In (b) a surface proton flux generates an electric field via reaction (R2). We assume the particle is positively charged, and it then swims in the direction indicated by the white arrow. This corresponds to the standard self-electrophoresis model, e.g., for Au-Pt swimmers. 2 The electric field has a dipolar form, like the proton flux.
In (c), we have a third mechanism, with equal fluxes of H + and HO 2 À ions. There is no net electrical current for this swimmer since there are equal positive and negative fluxes. However, self-propulsion still occurs. This is because the two ions diffuse at different rates (H + faster than HO 2 À ), and this creates a so-called diffusion potential, which acts to prevent net charge separation. The diffusion potential leads to an associated (self-generated) electric field, which then produces motion via electrophoresis in the usual way (white arrow). This propulsion mechanism is called 'ionic diffusiophoresis', 45 and is typically used to model swimmers composed of solid salts, which generate propulsion through dissolution of the swimmer itself, e.g., AgCl(s) -Ag + (aq) + Cl À (aq). 46 We include this model here because ionic diffusiophoresis may contribute to the propulsion of Pt-Janus swimmers, for example via with subsequent recombination of H + and OH À in the bulk. However, note that reaction (R3) is not the reaction shown in Fig. 4c, where HO 2 À is used instead of OH À for modelling simplicity.
The lower half of Fig. 4(d-f) shows the effect on each of these swimmers of a single ionic reaction occurring in the bulk, aqueous phase (R4) As we mentioned before, this reaction will only begin to have a significant effect when we are in the reactive, qa 4 1 regime. In Fig. 4, the reactive screening length q À1 is indicated by the dashed line: for these particular swimmers, qa E 3. The white arrows show the qualitative effect of this reaction on the propulsion speed, which is different for each of the swimmers. For (a-d), the reaction generates propulsion, for (b-e) the reaction increases the propulsion speed, and for (c-f) the propulsion speed falls. We now briefly explain the reason for these effects.
In the absence of a swimmer, reaction (R4) is in a state of dynamic equilibrium. If a swimmer consumes or produces molecules on either side of this equilibrium, then this will push the reaction out of equilibrium, and the system will respond so as to reduce the effect of that perturbation: this is Le Chatelier's principle. Thus, in Fig. 4d, the H 2 O 2 flux injected from the particle surface is to the left of the equilibrium, so the H 2 O 2 partially dissociates into ions, producing ionic fluxes. In Fig. 4e the proton flux is to the right of equilibrium, so there is net ionic recombination in the bulk to give an H 2 O 2 flux (this small H 2 O 2 flux is not shown because it does not significantly contribute to self-electrophoresis) and an HO 2 À flux in the opposite direction to the original proton flux. In Fig. 4f, the proton and HO 2 À fluxes are both to the right of equilibrium, so these both recombine with counterions in the bulk to give an H 2 O 2 flux instead of the ionic fluxes. For (a-d), the new ionic fluxes produce a diffusion potential, which generates motion. Hence a swimmer without any electrochemical reactions on its surface can still exhibit self-electrophoretic propulsion. Crucially, it will also display the experimental behaviour that would be expected of a selfelectrophoretic swimmer, e.g., propulsion speed scaling inversely with salt concentration (via the U SM factor in eqn (2)). This means that the kind of ionic behaviour observed in ref. 16, 17 and 32 does not a priori require an electrochemical surface reaction. In practice, however, we find that, because of the weak dissociation of H 2 O 2 , the simple non-electrochemical surface reaction mechanism in Fig. 4a-d cannot account for the magnitude of the experimentally observed propulsion in, e.g., ref. 16: genuine self-electrophoretic propulsion is still required.
In (b-e), the important point is that chemical reactions conserve charge, and this also implies the conservation of electrical current. There is a net electrical current in (b), and because this current is conserved it will have the same magnitude with or without bulk reactions. It is only the identity of the current-carrying ions which changes: in this case, the current becomes partially carried by HO 2 À ions travelling in the opposite direction, see Fig. 5. As we discuss later, the propulsion speed scales inversely with the diffusivity of the current-carrying ion. In this case HO 2 À diffuses approximately 10 times slower than H + , and this is why the speed increases. In fact, in the appropriate environment of high pH (= high HO 2 À concentration), the predicted speed increases ten-fold because the current becomes entirely carried by HO 2 À ions.
In (c-f), on the other hand, there is no net electric current to be conserved and both anions and cations react freely with their counterions in the bulk. Hence, far from the swimmer, the ionic gradients become vanishingly small, with a resultant drop in propulsion speed compared to the case without bulk reactions. In detail, the presence of ions in the bulk, due to the surface reactions, generates a diffusion potential (similar to the situation in Fig. 4d). However, since the ions in Fig. 4f can recombine through bulk reactions, the further one is from the swimmer surface, the fewer ions generated by the surface reaction remain to induce the diffusion potential. This shows up as an exponentially screened potential, with screening length q À1 . This also affects the swimming speed, because we find that the magnitude of the diffusion potential scales with the thickness of the screening layer, leading to a scaling of U p 1/(qa) for large swimmers. This lowers the predicted propulsion speed by up to a factor of approximately 100 for the largest E10 mm radius swimmers.

Comparison with experiments
As previously discussed, measurements of relevant surface reaction rates are not currently available. This makes it difficult for us (or others) to predict propulsion speeds that can immediately be compared with experiments. Nevertheless, we will present some speed calculations with the simplest assumption of fixed surface reaction rates in Section 7. These comparisons indicate that the U p 1/a scaling observed with Pt-polymer Janus particles 15 might come from an ionic diffusiophoretic mechanism like that discussed above, see Fig. 4c-f. In addition, our results indicate that the speed of nanoscale swimmers is too high to be explained by self-electrophoresis with fixed surface reaction rates, see Section 3.2.

Summary
In brief, we find that ionic reactions and electrostatic screening should have very significant and system-dependent effects on the propulsion of a wide range of chemical swimmers. These effects include increasing or decreasing the predicted speed by several orders of magnitude, as well as the qualitatively new behaviour of reactive screening. In addition, we find that even swimmers with no ionic surface reaction can behave as though they are self-electrophoretic. The remainder of this paper provides a detailed account of the theory that gives rise to these results, and compares them to experiments as far as current data allow.

Theoretical model
In this and the following two sections, we present a quantitative model of self-electrophoresis. Here in Section 4, we will lay out the general theoretical model and detail how this will be applied to the specific H 2 O 2 reaction system described above. In order to obtain analytical results we also linearize our theory. We will then apply this model to obtain explicit results, first for a system with only surface reactions (Section 5), and then with bulk reactions (Section 6).

General model
The standard theoretical approach to self-electrophoresis involves coupling the chemical fluxes arising from reactions on the particle surface to bulk differential equations (Nernst-Planck, Poisson, and Navier-Stokes). 14,34 This treatment generally ignores bulk chemical reactions by assuming that each chemical species is conserved. We adopt the standard approach, but include bulk reactions by coupling chemical fluxes to local reaction rates. We solve this model numerically using COMSOL. Separately, and in common with previous work, 33 we also linearize the model to obtain an analytical approximation. Unlike in previous work, the analytical solution does not require the assumption of a thin electrostatic screening layer.
We consider a spherical swimmer of radius a and uniform surface charge density s, see Fig. 6. The electrostatic boundary condition of such a particle iŝ with f the electrostatic potential field and e the dielectric constant of the fluid (the dielectric constant of the particle is assumed to be zero). Here (s) and n indicate evaluation at, and the normal out of, the particle surface. We have chosen a uniform, dielectric boundary condition for simplicity. However, in Appendix A.5 we show formally that, with an appropriate choice of surface potential, an equipotential (conducting) surface gives the same swimming speed as a dielectric. We do not deal with mixed dielectric/conducting particles here, but this should not qualitatively affect the basic physics of the self-propulsion. Propulsion is generated by reactions on the swimmer surface. These reactions produce and consume N different chemical species, labelled l = 1. . .N. The surface production (or, if negative, consumption) rate per unit area of each species is j s l (y), and is a function only of y, the polar angle with respect to the symmetry axis ẑ, see Fig. 6. The surface reaction rates can be equated to the bulk flux j l of each species out of the particle surface, giving the boundary condition These bulk fluxes obey the classical Nernst-Planck equation 34 j l ¼ c l u À D l rc l À D l z l e k B T c l rf; with z l , D l , and c l , respectively the valence, diffusivity, and concentration field of each chemical species, u the fluid flow field, e the fundamental charge, k B Boltzmann's constant, and T temperature. Physically, eqn (5) expresses the bulk fluxes as linear sums of advective, diffusive, and conductive terms respectively. Eqn (5) is the standard flux expression used in studies of self-electrophoresis. Its main simplification is the neglect of cross-coupling terms between the molecular fluxes, and this is valid as long as we are in the dilute limit with relatively small ionic gradients, 44 which is true here. Without bulk reactions, conservation of chemical species would require that the bulk fluxes are incompressible vector fields, i.e., rÁj l = 0. This is the standard approach, see ref. 14 and 34, and Section 5 here. Bulk chemical reactions can be incorporated by writing instead where R l is the local rate at which each chemical is produced (if negative, consumed), in chemical reactions. In general, R l depends on the local concentration of all chemical species involved in reactions with species l. Note that chemical reactions are charge-conserving, i.e., X l z l R l ¼ 0; everywhere, and combining this condition with eqn (6) implies the conservation of electrical current where the electrical current i ¼ e P l z l j l .
Infinitely far from the particle, the chemical concentrations are labelled c N l , and are determined by equilibrium equations and charge neutrality. The other boundary conditions at infinity are j N l = 0, f N = 0, u N = 0 (in the lab frame), and p N = p atm , where p is the hydrostatic pressure field and p atm is the atmospheric pressure, whose absolute value does not affect the calculations.
The electrostatic potential f is determined by the Poisson equation with charge density r e ¼ e P l z l c l . The interaction of the electrostatic potential and the unbalanced charge density (r e a 0) generates a force density and this drives fluid flow via the Stokes equations for lowinertia, incompressible flow Finally, the swimmer is not held in place, so fluid flow around it will cause it to move with some propulsion velocity U. This propulsion velocity is determined by the condition that there is no net force acting on the total system of swimmer plus fluid out to infinity. 47 This force-free condition is simply a reflection of the fact that all the forces are internal to this total system-there are no long-range, external forces like gravity. The force-free condition can be translated into an expression for U by using the Lorentz reciprocal theorem, which is a restatement of the Stokes equations in integral form. This gives a closed-form expression for the propulsion velocity 48 where r is the distance from the particle centre, r andŷ are unit vectors in the r and y directions, and the scalar speed U is defined by U = Uẑ. The volume integral is over the region outside the particle.

Numerical solution
We solve the full non-linear model numerically using FEM implemented in COMSOL. To do this, we make several modifications to the above equation system. In particular, we define a new force density f FEM to replace f in eqn (10). The two quantities are related by Fig. 6 Diagram of a model swimmer, highlighting the distinction between bulk and surface parameters. In the bulk, we have an electrostatic potential field f, chemical concentration fields c l , a pressure field p, and fluid velocity u.
Far from the particle, these fields approach uniform values (superscript N). On the particle surface, the uniform surface charge density s and the nonuniform molecular fluxes out of the surface j s l set boundary conditions for the bulk potential and concentration fields, respectively. The particle, of radius a, is axisymmetric around the z axis, so the surface fluxes are parameterized by the polar angle y. The swimming velocity, which we calculate, is U = Uẑ. This redefinition does not influence the result of the calculation, but limits spurious flows related to numerical artefacts in the electrostatics. 49 Our calculations are performed on an axisymmetric spherical domain of radius L = 10a + 25k À1 , which we verified to be sufficient to eliminate most finite-size effects in our speed calculations. This frame co-moves with the colloid. We impose no-slip at the colloid surface, and on the edge of the domain we employ the same boundary conditions as the theory has at infinity. For the fluid velocity we impose a no-stress condition on the edge of the domain with n the normal to the domain, * denoting transposition, and I the 3D identity matrix. This is equivalent to imposing a forcefree condition on the swimmer-fluid system. 44,50 Our technique is to first obtain approximate numerical solutions for the electrostatic and concentration fields in the absence of advection, so neglecting the first term in eqn (5). This approach is justified because experimental swimmers generally have low Péclet numbers, i.e., molecular diffusion D dominates over advection. The Péclet number is defined as Pe = Ua/D, with U E 10 mm s À1 , a E 1 mm, and D E 10 À9 m 2 s À1 typical for experiments on microswimmers, leading to Pe E 0.01.
The flow field is then computed self-consistently on the domain by employing the force density, eqn (14), following from the concentration and potential fields. The speed of the swimmer is determined by taking the average of the fluid velocity on the edge of the domain: U = Àhui r=L , where U is in the lab frame and u in the co-moving frame. We subsequently verified the low-Pe approximation by solving the fully coupled equations (with advection) directly in a limited number of cases, which gave agreement to within a few per cent. See Appendix C for full details of the numerical calculations.

Analytical solution
We also linearize the model to provide an analytical solution.
To do this we assume that the fields f and c l have only small deviations from their values in the uncharged, unreactive state where f = 0 and c l = c N l everywhere (for f, this assumption corresponds to the usual Debye-Hückel approximation, We then expand the model to linear order in the small dimensionless parameters c = fe/(k B T) and x l = (c l À c N l )/c N l . Applying this linearization to eqn (5) gives where the advection term has been dropped entirely because u scales quadratically with the small parameters (eqn (10) contains a product of r e and f, which are both small). We must also Taylor expand the production rates R l to linear order, i.e., where O (Á) means 'of order Á', and the elements are components of a matrix k which we can call the linear reaction matrix. Its meaning will become clearer when we consider specific reactions. From eqn (6) and (16), we have, to linear order X m k lm x m ¼ Àc 1 This set of N equations, together with the Poisson equation, which we rewrite as makes up a system of N + 1 linear differential equations in N + 1 fields (x l , l A 1. . .N and c).
This system of equations is soluble in a spherical geometry by standard spectral methods, and the electrostatic potential field so obtained can then be used to calculate the propulsion speed by evaluating the integral in eqn (13). In doing this, we make the further usual assumption of a relatively small driving field. 39 That is, if we define f = f eq + f sr where f eq is the electrostatic potential in the absence of surface reactions and f sr is the additional potential generated by these reactions then f sr { f eq . As a result, the surface reaction rates j s l , which only come into f sr , contribute linearly to the final velocity. The algebra required to solve eqn (19) and (20) is significant, so we go through this explicitly in Appendix A.

Specific H 2 O 2 reaction model
The chemical reaction system we consider is the simplified version of the H 2 O 2 reaction system described in Section 3. On the particle surface, H 2 O 2 decomposes into O 2 and H 2 O. For simplicity, however, we ignore both products of this reaction: O 2 because it is electrically neutral and does not dissociate, and H 2 O because it dissociates much less than H 2 O 2 -the respective equilibrium constants 51 are K eq,H 2 O = 1.0 Â 10 À14 mol L À1 (pH = 7) and K eq,H 2 O 2 = 2.5 Â 10 À12 mol L À1 . In the bulk, we ignore any slow decomposition of H 2 O 2 via reaction (R1), and the only bulk reaction we consider is the ionic dissociation reaction (R4), which we rewrite here We therefore have only three chemically active species, with associated subscripts in brackets: H 2 O 2 (), H + (+), and HO 2 À (À).
Protonation reactions like (R4 0 ) are normally extremely rapid, with kinetics controlled by the diffusion and collision of the ions, 52 and with simple first order rate expressions where we estimate the association rate constant to be k as = 4.9 Â 10 10 mol À1 L s À1 using the Smoluchowski-Debye theory for diffusion-limited reactions, see Appendix B. This theoretical value agrees closely with experimentally measured rates for similar reactions, 52 e.g., H + + HCO 3 À " H 2 CO 3 in water has k as = 5 Â 10 10 mol À1 L s À1 (ref. 53).
The dissociation rate constant k dis = 0.12 s À1 is then determined from the equilibrium constant K eq = k dis /k as = 2.5 Â 10 À12 mol L À1 (ref. 51). Far from the particle, the system is in equilibrium, so we have The production rates are R + = R À = k forward À k reverse and R = k reverse À k forward , and linearizing using eqn (18) gives the linear reaction matrix where the order of rows and columns is , +, À.
The three reactive species have diffusivities D = 1.7 Â 10 À9 m 2 s À1 (ref. 54), D + = 9.3 Â 10 À9 m 2 s À1 (ref. 55), and D À = 0.9 Â 10 À9 m 2 s À1 (ref. 56). We also have two unreactive ions, which we take to be Na + and Cl À with diffusivities D Na + = 1.3 Â 10 À9 m 2 s À1 and D Cl À = 2.0 Â 10 À9 m 2 s À1 (ref. 57). Because these ions are not involved in chemical reactions at the surface or in the bulk, their concentration fields are in equilibrium with the electrostatic potential. The implication is that the diffusivity of these ions does not contribute to the propulsion speed in the linear regime. We show this mathematically in Appendix A.1.
The chemical concentrations at infinity are determined by the chemical equilibrium, eqn (22) and by charge balance These two equations (eqn (22) and (24)) connect five concentrations, so we can set three concentrations freely. In practice, we choose instead to set the H 2 O 2 concentration, the total ionic strength, and the pH. The reaction scheme presented here is the simplest possible that gives the necessary freedom: bulk ionic dissociation reactions require at least three reactive species, and the two non-reactive ions are necessary to allow the ionic strength to be modified independent of other parameters. For the variable parameters, our base set, used unless specified otherwise, is 1 mmol L À1 salt, i.e., c 1 Na þ ¼ c 1 Cl À ¼ 1 mmol L À1 , a = 500 nm, and c N = 3 mol L À1 . For these parameters, k À1 = 10 nm, and c N + = c N À = 3 Â 10 À6 mol L À1 . These values were chosen because micron-sized particles and H 2 O 2 concentrations of order 3 mol L À1 are experimentally typical, 16 while the 1 mmol L À1 baseline salt concentration allows us to scan a wide range of the important parameter ka for realistically sized particles.
Meanwhile, the surface reactions are specified by surface fluxes j s l , l A {, +, À}, of the three active species. We consider the three model swimmers shown in Fig. 4, referred to as: S , the nominally neutral swimmer; S + , powered by a proton current; and S = , powered by ionic diffusiophoresis. As mentioned above, only the dipolar part of the fluxes, that is the 1st Legendre component, contributes to the propulsion speed of uniformly charged swimmers, 33 so we include only this term by setting j s l j s l,1 cos y for each surface flux where j s l,1 is a constant coefficient.
For S , only j s ,1 is finite; for S + , only j s +,1 is finite; and for S = , j s +,1 and j s À,1 are equal and finite, with j s ,1 = 0. Our model makes a number of simplifications. This includes those chemical simplifications already discussed, as well as the neglect of potential contaminants such as CO 2 , which also undergo ionic dissociation. The main purpose of this paper is to illustrate the physical principles behind the effect of bulk reactions on self-propulsion. These physical principles will also apply to a more complex and realistic H 2 O 2 reaction system, as well as to other chemical systems. 58,59 We also neglect any dependence of the surface parameters on environmental conditions, so, for example, we take s = constant, independent of pH, salt concentration etc. This does not imply that surface parameters are independent of the environment; it is just that detailed knowledge of this dependence is currently lacking. The 'pure' effect of bulk reactions which we capture will occur in addition to any such interdependence.

Electrophoresis without bulk reactions
Before discussing the effect of bulk ionic reactions, it is important to set out the basic theory for propulsion by self-electrophoresis without such reactions. This theory has been set out multiple times before for the limit of vanishing electrostatic screening length, ka c 1, ref. 14, 33, 34, 36, and 60. Here we extend the theory to include the effect of a finite ka.
For comparison, we first write down the standard results for electrophoresis in an external, linear field, 40 Fig. 3a (top). We align the external field along the z axis and consider the velocity U ext ẑ of a uniformly charged spherical particle with radius a and small surface (z) potential 40 in an externally imposed electric field E = Àrf. Far from the particle, E is a constant linear field E = E N ẑ. The particle is suspended in an aqueous solution of a monovalent salt, e.g., NaCl. In a weak field, particle velocity is proportional to electric field strength, U ext = m E E N , with m E called the electrophoretic mobility. The standard expression for m E for small z is 39,40 where f (ka) is Henry's function, 39 which accounts for electrostatic screening and depends only on ka, the ratio between particle radius a and the electrostatic screening length k À1 . The function f is plotted in Fig. 3b: it has constant limits of f (N) = 1, (the Debye or Smoluchowski limit) corresponding to high salt concentration or large particles, and f (0) = 2/3, (the Hückel limit) corresponding to small particles or non-polar solvents. Eqn (26), typically in either the high or low ka limit, is the expression commonly used to compute colloidal z potentials from mobility measurements in, e.g., commercial Zetasizers. In self-electrophoresis, the independent parameters are the surface reaction rates, and therefore the ionic fluxes, rather than the electric field. To facilitate understanding, we translate the expression for external electrophoresis into these terms. We write down expressions for the inverse electrostatic screening length the ionic conductivity and the concentration-averaged diffusivity together with Ohm's law which relates the electric field to the ionic current density at infinity i N ẑ and which we can rewrite as Combining eqn (31) with eqn (26) we then have Note that for electrophoresis in an external field, the particle speed is inversely proportional to the concentration-averaged diffusivity À D. We now compare eqn (32) with the analogous expression for the most well-studied self-electrophoretic swimmer, a protonpowered bimetallic swimmer. 14,36 Consider a spherical swimmer of radius a, with surface charge density s, a surface proton flux j s + (y) = j s +,1 cos y and no bulk reactions: we call this model S NBR + . In this case eqn (19) and (20) can be easily solved to yield, after some algebra where the first and second terms are f eq and f sr , the potentials generated, respectively, by the surface charge and the surface reactions. The Á Á Á indicate additional, electrostatically screened terms that are necessary to match the electrostatic boundary conditions, but which make no contribution to the propulsion, see Appendix A.5. The propulsion speed is obtained by evaluating eqn (13) with eqn (33) to give where F is, like f, a function of ka only. The full form of F is given in Appendix A.3. Eqn (34) corresponds closely to eqn (32), the particle velocity with external electrophoresis, and we compare these expressions factor by factor: (I) The relevant current density i N becomes Àj s +,1 e/3 because of the exclusive dependence of the propulsion speed on the first Legendre component of the flux 33 discussed in Section 4.4.
(II) The relevant diffusivity À D becomes D + because, for selfelectrophoresis in steady state, the ionic current can only be carried by the active ion involved in reactions at the particle surface, in this case H + . There can be no net flux of the other ions, or they would build up at the particle surface. In fact, the other ions are in local equilibrium with the electrostatic potential f, i.e., c Na þ ¼ c 1 Na þ exp Àf=k B T ð Þetc., from standard Debye-Hückel theory, 40 and the swimmer behaviour therefore cannot depend on their dynamic properties at all. The appropriate version of Ohm's law for the self-electrophoretic swimmer is thus not eqn (31), but instead 61 which depends only on the mobility of the active ion, and the electrical current. Eqn (33) and the propulsion speed have the same dependencies. This difference between external-and selfelectrophoresis has been confirmed in numerical calculations 60,61 and the inverse scaling of self-propulsion speed with the diffusivity of the active ion will be crucial for understanding the effects of bulk ionic reactions in Section 6.
(III) We have chosen to parameterize our model in terms of s rather than z because this is the most natural choice from a microscopic point of view. Much of the charge on the surface, both of conducting and dielectric particles, is due to surfaceabsorbed groups, leaving s fixed as other parameters, such as k, vary. This has been demonstrated experimentally for dielectric particles. 62 Nevertheless, the experimental evidence indicates that self-electrophoretic propulsion speeds scale with k À2 (ref. 11, 16 and 32) which is consistent with a fixed z, not a fixed s, though to our knowledge there is no microscopic justification for this. Since we are most interested in bulk effects, we do not insist upon a particular surface parameterization, and eqn (25) can be used to translate our results into a parameterization where z is fixed, which gives at small radius, a speed scaling as a 2 /k rather than a 3 as in eqn (36) below. Note that this point is distinct from the choice between conducting and dielectric boundary conditions on the particle surface, which is discussed in Section 4.1.
(IV) We have replaced f (ka) with an equivalent expression for self-electrophoresis, F(ka), shown in Fig. 3b. In the thinscreening limit, F(ka -N) = 1, and eqn (34) then agrees with previous self-electrophoresis results in the thin-screening limit, 33 except that ref. 33 incorrectly assumes that the propulsion is controlled by the total ionic diffusivity % D, (see point II above). In the opposite limit, F(ka -0) = (ka/2) 3 , so that for small ka the propulsion speed scales with a 3 The reason that F(0) -0, while f (0) is finite, is the different geometry of the driving currents, as illustrated in Fig. 3a. For self-electrophoresis, the driving potential is a local, dipolar field which decays over a length of order a, see eqn (33), whereas in external electrophoresis the driving potential is infinite in extent. Therefore, in self-electrophoresis, additional factors of a in the propulsion speed are to be expected. ‡ Several of the features of eqn (34) have been verified experimentally for bimetallic swimmers 2,4 like the Au-Pt swimmers in Fig. 1a, which explains the wide-acceptance of the selfelectrophoretic model for this system. As discussed above, this equation has also been found to be applicable to singlecatalyst swimmers such as Pt-polystyrene Janus particles, 11,16,17 suggesting that these swimmers are also powered by proton currents. 16,17 The additional screening parameter F(ka), is more problematic. It predicts that the speed of a swimmer will drop off sharply as ka decreases, Fig. 3. This drop-off is significant for surprisingly large ka: F(10) E 0.5, while F(1) o 0.1. This shows that the thin-screening limit, which is commonly employed, 16,17,21,33,63 is not justified even for the common situation of a 1 mm radius swimmer in 3 mol L À1 H 2 O 2 , where ka E 10 (ref. 16). However, to the best of our knowledge, there is no experimental evidence for this drop off. In fact, a small number of experiments show a larger speed for nanoswimmers 37,64 than is typical for microswimmers. 3,65 We discuss this experimental comparison in more detail in Section 7.

Electrophoresis with bulk reactions
We now examine the effect of bulk reactions on the propulsion of model swimmers, in particular of reaction (R4 0 ), H 2 O 2 " H + + HO 2 À , on the three model swimmers depicted again for convenience in Fig. 7a. In Section 6.1 we write down the general form of expressions for the swimming speed when bulk reactions are included, before focussing on the effect of two experimentally relevant parameters-swimmer radius and H + concentration-on the swimming speed, in Sections 6.2-6.5. In Section 7 we will compare our predictions with experimental observations.

General form of the solutions
Mathematically, the bulk reactions make it impractical to solve even the linearized problem by hand. Instead, we solve the system of equations, eqn (19) and (20) in the form of eqn (2). Here, † indicates a particular swimmer type, i.e., † A {, +, =} and j s † is the appropriate surface flux density for that swimmer. We define j s † = j s ,1 for the S swimmers, and j s † = j s +,1 for the S + and S = swimmers. The use of D + in the denominator of eqn (37)  Note that the expression in square brackets in eqn (37) is identical to eqn (34). This emphasizes that all the propulsion mechanisms, S , S + , and S = , are really forms of self-electrophoresis, and display all the responses to, e.g., salt-concentration, particle radius, and surface charge, which standard self-electrophoresis models predict. The inclusion of bulk reactions just adds a new layer of phenomena on top of this behaviour.

Influence of bulk reactions on swimming
We study the bulk reactions by varying two common experimental parameters: particle radius, and proton concentration, i.e., pH, Fig. 7b and c, respectively. Including bulk reactions (solid curves = analytic; solid symbols = FEM numerics) introduces a range of effects compared to the case with no bulk reactions (broken, horizontal lines), with qualitatively different behaviour for the three swimmer models.
Examining Fig. 7b first, the bulk reactions permit propulsion of the neutral swimmer S (inset), and B increases with radius, saturating for large radii. However, the magnitude of B always remains smaller than that of the other swimmers by a factor of order 10 À6 . In practice, this is typically partially compensated for by the much larger flux of the neutral species.
For the proton-current-driven swimmer S + , B + shows plateaux at both large and small radius, with the large-radius plateau approximately twice as high. For the ionic-diffusiophoretic swimmer S = , B = scales inversely with radius for large radius, but has a plateau at small radius. Meanwhile, varying the proton concentration c N + , as in Fig. 7c, produces a peak in B and B = , and decreases the overall value of B = by at least a factor of 5 compared to without reactions. For S + there are again two plateaux, at high and low c N + , with the low c N + plateau now a factor of approximately 10 higher than the other. The main control parameter for all these effects is qa, and there is a qualitative change of behaviour for all three swimmers at qa E 1: the vertical lines on Fig. 7b are for qa = 1. In Table 1, we write down the bulk parameters for each of the model swimmers in the limits qa { 1 and qa c 1. The full analytical expressions, which are lengthy, are provided in Appendix A.4, but the basic physics can be understood from the limiting behaviour. For the table, we have also assumed weak ionic dissociation, i.e., K eq { c N which is valid here, and thin electrostatic screening, ka c 1. These assumptions also apply to the analytical expressions given in the rest of this section. Table 1 matches Fig. 7 in all but one respect, which is the scaling of B at qa { 1, and this difference occurs because the assumption ka c 1 does not hold for small a in Fig. 7b. The parameters a and d* will be defined below. For  Fig. 7c, qa E 7 or larger, so we will assume that this figure is always in the qa c 1 limit.
To understand the results shown in Fig. 7 and Table 1, we will examine the bulk reactions in terms of three physical principles: reactive screening, the composition of the electrical current, and the dissociation of the neutral flux. Though we focus on these underlying principles, which are crucial for understanding the effect of bulk chemical reactions on any swimmer, this structure also allows us to discuss the three model swimmers in a logical order: S = , S + and S .

Reactive screening (model S = )
If an ion is released from the particle surface, it will react and come into local equilibrium with the surrounding solution. The characteristic distance over which this approach to equilibrium occurs can be called a 'reactive screening length' q À1 . As for the simple model discussed in Section 2.3, the reactive screening length is a balance between molecular diffusion and the reaction rate. However, the expression for q is more complex than in the simple model. We find Mathematically, q corresponds to one of the eigenvalues of the linear system of equations, eqn (19) and (20), see Appendix A.2. For our base parameter set, we obtain a screening length q À1 = 74 nm. Just as for the simple model, reactive screening gives an exponential decay of chemical concentrations with distance from the particle surface. The inclusion of charged species means that the electrostatic potential can now also become screened. However, we observe this reactive electrostatic screening only for S = swimmers, where the two ions released from the surface both react with oppositely charged ions in the bulk solution, causing an exponential decay in the resulting diffusion potential, see Fig. 8a.
For S + swimmers, no such screening is observed, see Fig. 8b. This is because the electrical current is conserved, so cannot be screened, and hence the associated electrical field also retains its unscreened dipolar form. Similarly, the H 2 O 2 concentration field c remains unscreened because this field is approximately conserved in the weak-dissociation limit. This results in an unscreened electrostatic potential field for S swimmers (not shown).
For S = swimmers, reactive screening also explains the 1/(qa) scaling of B = at high qa, as we show with a simple scaling argument: from eqn (4) and (5), we expect a fixed ratio between the surface reaction rates and the concentration gradients normal to the surface. For example, at r = a @c þ @r independent of other parameters. For qa c 1, the concentration decays exponentially away from the surface c + p exp(Àqr). (40) Differentiating this equation with respect to r gives qc + /qr E Àqc + and comparing this with eqn (39) yields Since the diffusion potential is proportional to the ionic concentrations (c + , c À ), and the propulsion speed U = is proportional to the diffusion potential, we have U = p q À1 . On the other hand, without bulk reactions the only relevant length scale is a, so a similar argument gives U NBR = p a. Therefore, one obtains B = p U = /U NBR = = 1/(qa). Physically, for qa c 1, the concentration flux only has the small screening length q À1 over which to set up a diffusion potential, whereas without bulk reactions a length of order a is available. This 1/(qa) scaling immediately explains the B = p 1/a scaling in Fig. 7b. We can also understand the peak in B = in Fig. 7c by noting that the screening length q À1 vanishes both for high c N + and for high c N À in eqn (38). Since c N À scales inversely Table 1 The bulk mobility factors predicted in the thin-screening ka c 1, ka c qa, and low dissociation c N + ,c N À { c N limits, for low, qa { 1 and high qa c 1 reaction rates. In both limits the prefactor should be multiplied by the relevant expression in the right-hand columns. The full expressions are given in Appendix A.4 with c N + due to the ionic equilibrium of eqn (22), this means that q À1 vanishes at either end of the c N + scale. As B = p q À1 , it too vanishes at either extreme and is peaked for intermediate c N + . Physically, at either end of the c N + scale, the high concentration of ions screens electric fields, preventing the formation of a diffusion potential.

Composition of the electrical current (model S + )
The total electrical current in the bulk is a conserved quantity, and is therefore not screened. However, the individual ionic fluxes making up that current are not conserved, and the bulk reactions modify the identity of the current-carrying ions. As discussed in Section 5, this is important because the swimming speed scales inversely with the diffusivity of the current-carrying ion or ions. In our system, an initially pure proton current will be partially replaced by HO 2 À ions travelling in the opposite direction, as illustrated in Fig. 5. In the reactive, qa c 1 limit we can calculate the composition of this electrical current relatively simply. From this, we will obtain the propulsion speed of the S + swimmer. In the qa c 1 limit, at any point outside the thin reactive screening layer, the ions released from the surface will have had time to come into equilibrium with each other. This is equivalent to requiring that the chemical production rates vanish, i.e., R l = 0. In the linear approximation, see eqn (17), this means X m k lm x m ¼ 0: In other words, the deviations in concentration x l of each of the reactive chemical species are coupled by the reaction matrix k given in eqn (23). This concentration coupling also implies a coupling of the chemical fluxes, in the same way that charge conservation implies the conservation of electrical current. Consider the linearized flux equation, eqn (16). Multiplying both sides by k ml /(D l c N l ) and summing over l yields À X l k ml j l D l c 1 l ¼ r X l k ml x l þ rc X l k ml z l : Then, from eqn (42), with l and m exchanged, the first term on the right vanishes, while charge conservation in reactions implies P l k ml z l ¼ 0 (see eqn (7)), so the second term vanishes too.
Hence, the general flux coupling equation is X l k ml j l D l c 1 For our specific system, substituting the expression for k from eqn (23) into eqn (44) then gives The physical meaning of eqn (45) is that each of the molecular fluxes has a characteristic scale set by D l c N l , and that, with this scaling, the relationship between the currents is set by the stoichiometry of the bulk reactions. Eqn (45) can be rearranged to give each of the ionic fluxes j AE in terms of the conserved quantities i and j where D* is the concentration-averaged diffusivity of the active ions § and the dimensionless factor a, which specifies the equilibrium decomposition of a neutral current into ionic currents, is The meaning of the first term in eqn (46), which is relevant for S + swimmers, is that the electric current is carried by a fixed proportion of H + ions travelling in one (positive) direction, and a counter current of HO 2 À ions in the opposite (negative) direction. In the second term, which is relevant for S swimmers, the neutral flux j continuously dissociates into H + and HO 2 À ions, producing small, equal fluxes of these ions, which travel with the neutral flux. If we are also outside the electrostatic screening length, which is the case in the ka c 1 limit, we also have a zero-charge-density condition, which reads X l c 1 l x l z l ¼ 0:

View Article Online
Just as above, but now multiplying eqn (16) by z l /D l , we can derive a direct relationship between the electric field and the chemical fluxes 61 rf ¼ À e k 2 e X l j l z l D l ; (50) which, combined with eqn (46), yields a version of Ohm's law for the reactive limit Comparing the first term of this equation with eqn (35) for self-electrophoresis without bulk reactions, we see that they are identical apart from the switch from D + to D*.
For an S + swimmer, j = 0, and for high qa all of the electric field outside a thin screening layer will be determined by the first term of eqn (51). From this we can understand the 1/D* factor which appears in B + , see Table 1. Just as without bulk reactions (see eqn (34)) the propulsion speed is inversely proportional to the diffusivity of the current-carrying ion. However, the current is now made up of two ions, with a total effective diffusivity D*. This explains the 2Â speed increase in Fig. 7b: D + at low a is replaced by D* at high a, and for c N + = c N À , which is the case in Fig. 7b, D* = (D + + D À )/2 E D + /2.
To understand the effect of varying c N + , Fig. 7c, we examine the form of D* in eqn (47). At high c N + , D* E D + , while at low c N + (=high c N À ), D* E D À . Physically, this is again simply a result of the relative number of each ionic species: if there is an overwhelming number of protons in solution, then the ionic current must be carried predominately by protons. This explains the factor of D + /D À E 10 speed difference between the two plateaux for B + in Fig. 7c.

Dissociation of the neutral flux (model S )
To understand the dissociation of a purely neutral flux, we examine the parameter a in eqn (48). The form of a can be explained by the fact that in the absence of a net electrical current, e.g., for S swimmers, the ionic currents are constrained by j + = j À . This means that the total ionic flux will be limited by whichever ion has the lower value of D l c N l , as this ion will contribute most to the flux balance in eqn (45). Hence the parameter a, like q À1 , vanishes at the extreme ends of the c N + scale: at low c N + it is limited by the low proton concentration, and at high c N + , by the low HO 2 À concentration.
The dissociation of the neutral flux generates a diffusion potential. Hence, the prefactor for B in Table 1 is made up of two factors: a and (d + À d À )/(d + d À ), the latter of which controls the diffusion potential just as for B = . The peak in B as a function of c N + then follows directly from the behaviour of a. Interestingly, both S and S = show peaks in speed at intermediate c N + , but for two different reasons. For S = , the reason is that the reassociation of ions is slowest at intermediate concentrations. For S , the reason is that the least conductive fraction of the solution limits the total carrying capacity, and this effect is strongest at either extreme in pH.

Comparison with experiments
We now compare our theoretical predictions with experimental results, in so far as this is possible at present. We stress here again the lack of understanding of the surface chemistry, and in particular of the effect of experimental parameters on the ionic surface reaction rates. We have not attempted to predict these reaction rates, so we cannot immediately test our theoretical predictions. In this section alone, we will make the simple and typical 3,15,16 assumption that the surface properties, i.e., surface reaction rates and surface charge densities, vary only with fuel concentration and are otherwise constant, so that U p k À3 F(ka)B(qa) for all swimmers (eqn (37)). This allows us to make some suggestive comparisons with experiments. We note that more complex reaction rate dependencies based on electrochemical modelling of the surface have been proposed previously. 14,17,36 However, to the best of our knowledge, these models do not have independent experimental justification, or experimental validation from speed experiments beyond that achieved by the assumption of uniform surface properties.
Independent of bulk reactions and swimmer type, we predict a speed scaling with a 3 for small particles. In particular, for a proton powered swimmer U NBR + , the predicted speed dependence is as shown in Fig. 9a (solid curve), due mostly to the new electrostatic screening parameter F(ka)-the bulk reactions do not significantly modify the form of this curve. Here, we have matched the solution parameters to those of the typical microswimmer experiments of ref. 65, choosing the (constant) surface parameters to give a speed of 7 mm s À1 (blue circle) for a 1 mm radius bimetallic sphere, as in that paper (see Appendix B for the surface parameters used). We also plot nanoswimmer data (black triangle) from ref. 37, where a swimmer of radius 10 nm had U = 650 mm s À1 . These two experiments used similar concentrations of H 2 O 2 , but differed in salt concentration. Thus, to predict the nanoswimmer results, we also plot a theoretical curve (dashed), with the salt concentration modified to match those of the nanoswimmer experiment, 37 but keeping Fig. 9 Comparison between theory and experiments. (a) The predicted speed of a swimmer powered by a proton current, in the presence of chemical reactions, with parameters chosen to match typical measurements on microparticles (red solid, theory; blue circle, experiment 65 ); and, with the same surface parameters, but bulk solution parameters chosen to match experiments on a nanoswimmer (red dashed, theory; black triangle, experiment 37 ). (b) The speed of an S = swimmer plotted against experimental data on Janus-Pt microswimmers. 15 The experimental error bars are smaller than the data points.
the surface parameters the same as in the microswimmer experiments. 65 If the assumption of constant surface properties holds true, then this dashed curve should agree with the experimental value for nanoswimmers. Instead, there is a clear disparity amounting to several orders of magnitude. This is not the result of our linear approximation: we find a good match between analytics and numerics up to values of s and j s + higher than those used in plotting Fig. 9a (Appendix C). The disparity is also not significantly reduced if we assume uniform z rather than uniform s, see Section 5. The discrepancy could be explained in at least two ways. It may be that self-electrophoresis is not the correct propulsion mechanism for bimetallic nano-swimmers. It has recently been found that nanometre scale biological enzymes also exhibit self-propulsion, 41,42 and a range of mechanisms has been proposed for this propulsion, 66 some of which might also apply to bimetallic nano-swimmers. Alternatively, it may be that the assumption of constant surface properties is inappropriate. That is, the proton current density could be much higher for these nano-swimmers than for micro-swimmers. Whatever the explanation, this discrepancy highlights the need for more systematic studies of identical or comparable swimmers over wide parameter ranges, as in ref. [15][16][17], and for independent measurements of the relevant ionic reaction rates.
For Pt-polystyrene Janus particles, such systematic studies do exist. 15 These show a U p a À1 scaling for 0.2 mm o a o 5 mm (this scaling has also been observed over a narrower range for some bimetallic swimmers 65 ). Self-electrophoresis S + , reaction (R2), is currently the preferred mechanism for Pt-polystyrene Janus swimmers, 16,17 but comparison of this 1/a scaling with Fig. 7b suggests self-ionic diffusiophoresis S = as an alternative mechanism, corresponding to reaction (R3). This is plausible: ion release without net electrical currents, which would correspond to reaction (R3), has previously been observed for H 2 O 2 decomposition on Pt. 67 This mechanism would also avoid the conceptual difficulty of producing a net ionic current in singlecatalyst systems. 16,17 However, when we plot the experimental data from ref. 15 against our theoretical predictions for S = propulsion, Fig. 9b-which is again scaled to match the experimental data for 1 mm radius swimmers, see Appendix B-we see that the fit fails at small a, again due to the F(ka) parameter. It is possible that evaluation of the complete H 2 O 2 -H 2 O reaction system would provide a better fit, but this goes beyond the scope of this work.
Note that the 1/a scaling has previously been explained by postulating that the overall surface reaction rate j s is limited by diffusion, 15 and therefore scales as 1/a just from geometrical arguments. However, the diffusion-limit implies a large flux density j s E D c N /a, which for a 1 mm radius swimmer in 3 M H 2 O 2 , as in ref. 15 requires j s E 3 Â 10 24 m À2 s À1 . So far, only much smaller rates, j s E 10 22 m À2 s À1 , have been measured, both by us 16 and by the authors 17 of ref. 15. Therefore, these swimmers do not appear to be in the diffusion limited regime, so this explanation for the 1/a scaling cannot hold.
Next, we have previously calculated the values of the uncharged flux j s and the charge density s for Pt-coated Janus swimmers. 16 We estimated that the propulsion speed of such swimmers was too high to be explained by a purely uncharged reaction like (R1). 16 This estimate did not allow for bulk ionic reactions. However, including these reactions, we calculate in Appendix C that such a mechanism could still only account for E5% of the observed speed of these swimmers. Hence, a model with just surface reaction (R1), even with bulk dissociation, cannot explain the propulsion of H 2 O 2 -powered swimmers, so that such swimmers probably still require more complex ionic surface reaction schemes like (R2) and (R3). Nevertheless, purely neutral-surfacereaction mechanisms could still be relevant for swimmers powered by more dissociative fuels, such as hydrazine. 59 Turning to the effect of pH, there have been two suggestive studies, 16,59 but no systematic investigation. First, we found that NaOH reduced the swimming speed of Pt-polystyrene Janus swimmers, but that this effect was much weaker than the speed reduction due to NaCl. 16 This is consistent with our prediction that increasing pH at fixed Debye length should raise the swimming speed for any of the 3 swimming mechanisms discussed. Raising the pH corresponds to moving left from the Â symbol in Fig. 7c; for all swimming mechanisms the value of B increases in this direction.
Second, the silica-iridium swimmers of ref. 59 show a clear spike in speed as a function of fuel (hydrazine) concentration similar in form to the peaks in Fig. 7c. This spike could be due to modulation of pH by the reaction product ammonia, which would imply that either neutral self-diffusiophoresis or self-ionic diffusiophoresis dominates this swimmer's self-propulsion. In both these experiments, however, variation of the reaction rates with pH could also explain the results, 68 so further systematic study is necessary.
Finally, we note that the bulk ionic dissociation rates are also not well known, see Section 4.4. Our estimate of the association constant k as is a diffusion limited rate, and is therefore an upper limit, though for a proton exchange reaction, this is likely an accurate estimate. 52 Nevertheless, it might be argued that the effects of bulk reactions will not be observed in practice for lower reaction rates, as we require qa ] 1. However, since q À1 p k as À1/2 , see eqn (38), for a much lower value of k as , say 1000 times lower, we still find q À1 B 2 mm, well within the experimental microswimmer size range. Further, the value of q À1 becomes arbitrarily small in both high and low pH solutions, see eqn (38).

Conclusion
In this article, we have theoretically explored the influence of common polar solvents such as water and hydrogen peroxide on the propulsion behaviour of chemically-propelled, synthetic microswimmers. We have focussed on two unavoidable properties of such polar solvents-electrostatic screening and ionic dissociation-and calculated their effect on the swimming speed of a wide range of microswimmers propelled by chemical reactions on their surfaces. These effects have not been studied systematically before; nevertheless, they are highly significant, and including these effects can modify predicted swimming speeds by several orders of magnitude.
By ionic dissociation, we mean the breaking up of neutral molecules, including water, into charged species, for example, H 2 O 2 " H + + HO 2 À . One of our main predictions is that this kind of ionic dissociation reaction allows even microswimmers whose surface chemistry does not involve any ions, e.g., swimmers propelled by the simple decomposition of hydrogen peroxide, 2H 2 O 2 -2H 2 O + O 2 , to generate ionic gradients and thereby electric fields. The implication of this is that all microswimmers in water should experience some degree of self-electrophoresis, i.e., propulsion via self-generated electric fields. This is significant because self-electrophoresis is much more efficient than other putative propulsion mechanisms and is likely to dominate over them. Put simply: our results imply that all swimmers in aqueous solution are likely to be selfelectrophoretic to a large degree.
The second major prediction of our work is that for some chemically-propelled swimmers, ionic dissociation reactions may result in a kind of exponential 'reactive-screening'. 20 Electrical and chemical concentration fields generated by surface reactions on microswimmers are usually taken to decay slowly into the bulk solution, that is, as a power-law with distance. Ionic dissociation can instead produce a short-ranged exponential decay of these fields, just as in electrostatic screening. This is significant since these chemical and electrical fields are implicated in inter-swimmer interactions and collective behaviour, and the interaction range will play a crucial role in this behaviour.
Our third prediction relates to electrostatic screening itself. Most analytical treatment of microswimmers has focussed on the thin-screening limit, where the electrostatic screening length is much smaller than the swimmer size. For very small swimmers, this limit does not apply, and we find that this massively reduces the predicted swimming speed. This is important because experiments on nanoscale swimmers 37 show that these in fact swim faster than microswimmers, in apparent contradiction to our predictions. This opens up the exciting possibility that nanoscale swimmers move by entirely novel mechanisms compared to their microscopic counterparts.
Finally, the general conclusion that we draw from our results is that more experimental work is required to understand selfpropulsion mechanisms. The effect of ionic dissociation in particular depends crucially on the type of surface reactions which are responsible for propulsion-and the details of these reactions remain almost universally unknown. What is most urgently required in this regard is the independent measurement of surface reaction rates, which is challenging, and has so far only been achieved in the simplest of cases. However, recent results with electroosmotic pumps 69 suggest that such measurements will not long remain beyond our reach. We particularly hope that our theoretical results will lead to renewed efforts in this direction.
Looking ahead, our results suggest that a deeper understanding of self-propulsion will lead to greater insights into swimmer-swimmer interactions and collective effects. This is particularly relevant to synthetic swimmers, as their propulsion is closely coupled to their interactions through self-generated electrostatic, chemical, and hydrodynamic flow fields. We have shown here that reactive screening can qualitatively change the electrostatic interactions between swimmers. A detailed follow-up study will look explicitly at such interactions. Further theoretical work will focus on applying our calculations to fully realized experimental systems, e.g., mixed metal-dielectric swimmers.

Appendix A Calculation of the analytical solution
In this appendix, we explicitly calculate the propulsion speed for a general swimmer. This calculation is based on the linearized model described in Section 4, but now with a more concise notation, A.1. From this linearized model, we determine first the electrostatic potential fields, A.2, then the propulsion speed, A.3, as set out in Section 4. In A.4, we apply this general calculation to determine the speed of the model swimmers presented in the main text. Finally, in A.5, we demonstrate the equivalence of uniform charge and uniform potential boundary conditions for the calculation of the propulsion speed.

A.1 The linearized model
We begin with the linearized model described in Section 4.3. For notational convenience, we define a composite dimensionless parameter y l by combining the linearized potential c and concentration x l fields With this notation, the linear system of equations, eqn (19) and (20) is given by Eqn (53) represents a system of N + 1 linear equations. However, several of the species, typically inactive ions such as Na + or Cl À , may not be involved in any bulk or surface reactions, and we will now show that these inactive species can be eliminated. We specify that the first N 0 indices (N 0 o N) correspond to the reactive species. For the remaining, unreactive species, all the bulk reaction coefficients, k lm are zero, and there is no surface flux, so eqn (53) can only be satisfied if y l = Àz l y 0 , l 4 N 0 .
This is the linear approximation to the Boltzmann distribution, which one expects, since these unreactive species should be in equilibrium with the electric field. Using eqn (54), these ions can be eliminated from the remaining N 0 + 1 parts of eqn (53) to yield k À2 r 2 y l ¼ where k is the inverse Debye screening length with the Bjerrum length, l B = e 2 /(4pek B T), and where w l is a dimensionless ionic concentration Eliminating the inactive ions makes it clear that the motion of the swimmer cannot depend on the diffusivity of these ions, and is only affected by them through the value of k and through charge balance. Finally, linearizing the boundary conditions in eqn (3) and (4) givesn A.2 Calculation of the electrostatic potential Eqn (55) has the form of a matrix equation with components corresponding to the chemical concentrations and the electrostatic potential, so it is convenient to introduce some additional matrix notation. The bold font is reserved for real-space vectors, such as the fluid velocity u, while vectors in this concentrationpotential space will be underlined. A general vector À t will have N 0 + 1 components labelled t l , while a matrix T will have (N 0 + 1) Â (N 0 + 1) components labelled T lp . A point in the concentration-potential space is specified by the vector À y, with components y l , as defined in eqn (52). Using this notation, we can rewrite eqn (55) as which can be solved by finding the N 0 + 1 eigenvectors of the matrix M, with eigenvalues m p . These eigenvectors define a new basis, in which M is diagonal. Defining À w as the representation of À y in this basis, we have where the matrix G is diagonal, with components G lp = d lp g p , where d lp is the Kronecker delta and g p ¼ ffiffiffiffiffi m p p is an inverse screening length. For the model system described in the main text, the unique values of g p are k, 0 and q, as given in eqn (38). Here, for clarity, we use the index p to refer to the screening lengths, and the indices l or m to refer to the concentrations and potentials, even where these are dummy indices.
Eqn (60) is a series of N 0 + 1 independent Helmholtz equations, and the full solution to this equation is just a vector of individual solutions to the Helmholtz equation. In spherical polar coordinates, these solutions have the form 20,70 w p ¼ X n w p;n P n ðcos yÞ a r nþ1 T n g p r À Á T n g p a À Áe ÀgpðrÀaÞ ; with w p,n an as yet undetermined surface coefficient, P n the Legendre polynomial of order n, 71 y the polar angle, and T n ðxÞ ¼ X n s¼0 2 s n!ð2n À sÞ! s!ð2nÞ!ðn À sÞ! x s : We refer to the Legendre components by the subscript n throughout, and where we have multiple subscripts, the Legendre subscript shall be preceded by a comma. Transforming back into the original coordinate frame linearly combines the solutions in eqn (61), so that the final form for the electrostatic potential is f ¼ X p;n f p;n P n ðcos yÞ a r nþ1 T n g p r À Á T n g p a À Áe ÀgpðrÀaÞ ; with analogous expressions for each concentration field. Here, f p,n are surface coefficients which we will now determine. Transformation back into the original coordinate system is achieved with a transformation matrix K where each element K lp of K is equal to the lth component (in the original coordinate system) of the pth eigenvector. Applying this transformation to eqn (61) gives y l ¼ X p;n K lp w p;n P n ðcos yÞ a r nþ1 T n g p r À Á T n g p a À Áe ÀgpðrÀaÞ : The boundary conditions specified in eqn (58) can also be rearranged into a matrix equation where À b is a vector specifying each of the boundary fluxes or charge density. We define the harmonic components À b n of À b by b ¼ P n P n ðcos yÞb n , with analogous expressions defining B n . The solution to the boundary conditions is found by inverting eqn (66) to yield where where Y(ka,qa) depends on the relationship between the 3 length scales a, k À1 and q À1 and is Fðka; qaÞ FðkaÞ : With the limits k c q, and either qa c 1 or qa { 1, we obtain the expressions given in Table 1.

A.5 A note on electrostatic boundary conditions
In this section, we show that a particle with fixed, uniform surface charge s has the same propulsion velocity as an equivalent particle with fixed, uniform surface potential z, as long as To do this, we first need to show that modifying the electrostatic boundary conditions of the particle has only a limited effect on the fields of concentration and potential; namely, that such modifications can only introduce electrostatic fields corresponding to the equilibrium Debye-Hückel solutions around passive colloids.
We take a swimmer, in a given chemical environment, and apply to it three sets of boundary conditions. Boundary conditions (1) and (2), with corresponding solutions À y (1) and À y (2) , have equal chemical flux boundary conditions (equal surface reaction rates), but have arbitrary, different electrostatic boundary conditions. Boundary condition (3) consists of a no flux condition on all species (no surface reactions), and the electrostatic boundary condition y (3) 0 (s) = y (2) 0 (s) À y (1) 0 (s).
Since there are no fluxes through this particle's surface, each chemical species is in equilibrium, and the solution to this boundary condition is just the equilibrium, Debye-Hückel solution where the equilibrium potential field c must satisfy both the electrostatic boundary condition, eqn (85), and the Debye-Hückel equation One can then show by direct substitution of eqn (86) into eqn (55), that the solutions to the three boundary problems are related by À y (2) À À y (1) = À y (3) . In particular, f (2) À f (1) = c, which implies, from eqn (87) Hence the difference f (2) À f (1) between the electric potential fields of particles (1) and (2) corresponds to an equilibrium Debye-Hückel solution around a passive colloid. As before, we make the assumption of a small driving field, f sr { f eq , where f = f eq + f sr . Now, consider two particles (1 0 ) and (2 0 ), with equal surface reactions, but where (1 0 ) has uniform surface charge density s, and (2 0 ) has uniform surface potential z, with s and z satisfying eqn (84). In this case, the two equilibrium fields are equal, i.e., f (1 0 )eq = f (2 0 )eq , and are given by eqn (78). Subtracting this equality from eqn (88) yields, for the remaining, non-equilibrium part of the potential In other words, the difference in the reaction-generated electrostatic potential field between (1 0 ) and (2 0 ) is an equilibrium, Debye-Hückel type field, which has an inverse screening length g p = k. Since there is a (k À g p ) factor in eqn (79), we see that such a field can have no effect on the propulsion speed. This proves the assertion that, to linear order, a particle with fixed, uniform surface charge s will have the same propulsion velocity as an equivalent particle with fixed, uniform surface potential z, as long as eqn (84) is satisfied. In fact, one can make a more general statement, which we will not prove. For any two particles (1 0 ) and (2 0 ), with equal arbitrary shape, surface reactions, and equilibrium (possibly non-uniform) fields f eq , not only the propulsion speed but the entire flow field will be the same. A physical justification for this conclusion is that if the interaction between one equilibrium field (f eq ) and another (the difference field between (1 0 ) and (2 0 )) could generate fluid flow, then this would constitute a perpetualmotion machine. Analogous conclusions have also been drawn for electrophoresis in an external field. 72

B.1 The ionic association constant
An important parameter in our calculations is the ionic reaction association constant k as in (R4 0 ), H + + HO 2 À " H 2 O 2 in water, see Section 4.4. We were unable to find a value for this constant in the literature. However, reactions involving the transfer of a proton or a hydroxyl ion are normally sufficiently fast to be diffusion limited. 52 It has been shown, 73 that the diffusion-limited rate constant between two species, A and B, with diffusivities D A , D B , and valences z A , z B , which react at a short distance r AB is 74 Here, f (z A z B ,r AB ) is a modifier for charged species For reactions between oppositely charged species, over a typical reaction distance in water of r AB = 0.2 nm, 74 f (À1,r AB ) = 3.59. For the reaction between H + (species A) and HO 2 À (species B), this yields, using the diffusivities quoted in the main text, k as = 4.9 Â 10 10 mol À1 L s À1 .

B.2 Comparison with experiments
For the comparison of self-electrophoretic micro-and nanoswimmers in Fig. 9, we take experimental parameters from ref. 65 (microswimmers) and ref. 37 (nanoswimmers). For microswimmers, we use a = 1 mm, c N = 1.5 mol L À1 and no added salt. For nanoswimmers, we use a = 15 nm, c N = 1.5 mol L À1 , pH 7 and 5 mmol L À1 NaCl (this has the same Debye length as 1 mmol L À1 trisodium citrate, which was used in practice 37 ). For both swimmers, we take s = 10 À2 e nm À2 (ref. 21) and j s +,1 = 1.66 Â 10 À7 mol m À2 s À1 , which is chosen to match the microswimmer speed in ref. 65.
For comparison between S = and polystyrene-Pt Janus particles, we take experimental parameters from ref. 15, which are c N = 3 mol L À1 with no added salt. We used s = 10 À2 e nm À2 again and j s +,1 = 6.42 Â 10 À6 mol m À2 s À1 , which is chosen to match the experimental propulsion speed at a = 1 mm s À1 .

C Finite element method calculations
In this section we give additional details for the numerical FEM calculations discussed in Section 4.2. FEM calculations are performed using the COMSOL 5.1 Multiphysics Modelling package.
We employed the following strategies to accelerate the calculations and obtain high quality results. (i) The solutions were obtained in a 2D cylindrically symmetric geometry. (ii) We ignored the advective coupling term in eqn (5). This allowed us to split the problem into electrostatic plus hydrodynamic parts, as for the linear theory, and thus solve the uncoupled equations more efficiently. This approach is justified, since the Péclet number (Pe) t 10 À2 for typical experimental systems. We also verified this directly, by including the advective coupling term in a subset of the data points, finding good agreement. (iii) We created a physics-specific mesh, see Fig. 10, on which we solved the system. Quadrilateral elements were used out to a distance of 3k À1 from the colloid surface. These elements grow exponentially in size with increasing distance, whilst maintaining a constant number along the tangential direction. The remainder of the domain was meshed with triangular elements which grow larger with distance from the colloid. This approach is necessary to ensure convergence of the model. (iv) The following polynomial orders were employed for the test functions: electrostatics (3), diffusion (5) and hydrodynamics (2 + 3). These higher orders proved necessary to reduce spurious flow (see also ref. 44).
(v) Finite-size scaling was employed to check for artefacts arising from the finite extent of the simulation domain, we found that for L = 10a + 25k À1 the effects on the speed of the particle were negligible. (vi) Mesh refinement was used for several simulations to determine the dependence of our result on the element size. (vii) We also varied the tolerance on the residual for a few cases to verify that our solutions had sufficiently converged.
To verify the analytic results, we first performed calculations with sufficiently low values of the surface charge density and flux to remain in the linear regime. These j s and s are given in Table 2. Fig. 7 in the main text shows that there is excellent correspondence between the theory and FEM calculations in this regime. Different fluxes were used for the different propulsion models because the low efficiency of type S and S = propulsion mean that numerical errors become significant more quickly as the flux density is reduced for these models.
In addition, the FEM calculations and the linearized theory produce essentially identical electrostatic potential fields. Fig. 11a illustrates this for type S + electrophoresis. Note that we had to use a much smaller computational domain than we typically use (L = 6a here, rather than L = 10a + 25k À1 ), in order to show details in Fig. 11a. This means that the deviation from the theory, which stems from the f = 0 boundary condition on the edge of the domain (Fig. 11b), occurs closer to the particle than in our regular calculations. However, the potential and Fig. 10 The mesh on which the FEM calculations are performed. This particular mesh was generated for radius a = 0.5 mm and a salt concentration of 10 À5 mol L À1 , but illustrates the generic features of all the meshes. The rotational symmetry of the simulation domain is exploited to calculate on a quasi-2D domain: the symmetry axis is indicated by the dashed red line. The domain typically has a radius L = 10a + 25k À1 in size. This domain is subdivided into two pieces on which triangular and quadrilateral elements are used. In a range of 3k À1 around the colloid the domain consists of quadrilaterals, which grow in size geometrically, see the zoom-in (blue box). Beyond this range the elements are triangular and are allowed to grow out linearly to best fit the domain boundary and reduce the overall number of elements. Table 2 The charge densities and the first Legendre components of the surface flux densities used in Fig. 7 in the main text and Fig. 12 here. The flux densities have units mol m À2 s À1 , and the charge densities have units e nm À2 . The final column gives the product of s and the relevant non-zero flux density, with units e mol m À2 nm À2 s À1 s j s 7 S 3 Â 10 À1 0 0 1 0 À4 3 Â 10 À5 S + 0 3 Â 10 À7 0 1 0 À4 3 Â 10 À11 S = 0 3 Â 10 À5 3 Â 10 À5 10 À4 3 Â 10 À9 12 S 1.5 Â 10 À2 0 0 1 0 À2 1.5 Â 10 À4 S + 0 1.5 Â 10 À5 0 1 0 À2 1.5 Â 10 À7 S = 0 1.5 Â 10 À5 1.5 Â 10 À5 10 À2 1.5 Â 10 À7 flow-fields decay sufficiently rapidly that this does not affect the potential near the particle, or the propulsion speed beyond a few per cent. We can also use the FEM to go beyond the linear approximation. We defer to future work a systematic investigation of the non-linear behaviour, and here focus on the propulsion speed for selected experimentally relevant values of the surface charge density and chemical fluxes. These values are taken from measurements on the Pt-polystyrene Janus swimmers in ref. 21, and are listed in Table 2. The neutral flux density j s ,1 is that which would be produced by a Janus particle which uniformly consumes H 2 O 2 on one hemisphere at a rate G = 8 Â 10 10 molecules per second per particle. This rate was measured for a = 1 mm radius particles in 3 mol L À1 H 2 O 2 . 16 The surface charge density is taken from the electrophoretic mobility measurements made on the same particles in ref. 16. The ionic fluxes are unknown, but we arbitrarily set j s AE,1 = 10 À3 j s ,1 , so that S + electrophoresis gives a speed of order 100 mm s À1 , which is somewhat larger than typical experimental values, E10 mm s À1 for Au-Pt spherical microswimmers. 10,65 Hence, our results should overestimate the non-linear behaviour of the propulsion speed. Note that though the ionic flux densities for the experimentally realistic case are sometimes lower than those for the linear case, the product of charge density and surface flux is always greater in the experimentally realistic case, Table 2. Fig. 12a and b, both with 1 mmol L À1 NaCl, correspond to Fig. 7 in the main text. We see that the analytical theory continues to match the FEM calculations well even for these realistic values of the flux and charge densities. However, many experiments are performed with no added salt, and as shown in Fig. 12c, the agreement worsens as the salt concentration falls. This is to be expected, since it is low salt that generates a high-z, large-screening-length regime where linear approximations break down. 75 In fact, with 0 mmol L À1 NaCl, the dimensionless zeta-potential ze/(k B T) = 5.6 for these particles, well beyond the Debye-Hückel regime of ze/(k B T) { 1. Nevertheless, for all propulsion types, the agreement remains semi-quantitative between simulations and theory over the whole radius range for 0 mmol L À1 NaCl, Fig. 12d.
From Fig. 12d, we obtain a speed of 0.5 mm s À1 for type S electrophoresis with particles of radius a = 1 mm, no salt, and Fig. 11 (a) Comparison of normalized surface-reaction-generated potential fields f sr for (left) FEM and (right) linearized theory, for S + self-electrophoresis, with j s +,1 = 3 Â 10 À7 mol m À2 s À1 , and with other conditions as in the base parameter set (see main text). The radius of the simulation domain L = 3 mm = 6a here. (b) Normalized radial decay of f sr for linearized theory ( ) and FEM calculations ( ) along X-X 0 in (a). Fig. 12 Propulsion speed for realistic parameters (given in this text) for type S 0 ( , inset), S + ( ), and S = (K) propulsion, from analytical theory (solid curves) and FEM simulations (symbols). For (a) [H + ] at fixed k equivalent to 1 mmol L À1 NaCl, (b) particle radius with 1 mmol L À1 NaCl (c) NaCl concentration, (d) particle radius with 0 mmol L À1 NaCl. In (d), the black arrow indicates the experimental point from ref. 16 referred to in the text.