Shape-directed rotation of homogeneous micromotors via catalytic self-electrophoresis

The pursuit of chemically-powered colloidal machines requires individual components that perform different motions within a common environment. Such motions can be tailored by controlling the shape and/or composition of catalytic microparticles; however, the ability to design particle motions remains limited by incomplete understanding of the relevant propulsion mechanism(s). Here, we demonstrate that platinum microparticles move spontaneously in solutions of hydrogen peroxide and that their motions can be rationally designed by controlling particle shape. Nanofabricated particles with n-fold rotational symmetry rotate steadily with speed and direction specified by the type and extent of shape asymmetry. The observed relationships between particle shape and motion provide evidence for a self-electrophoretic propulsion mechanism, whereby anodic oxidation and cathodic reduction occur at different rates at different locations on the particle surface. We develop a mathematical model that explains how particle shape impacts the relevant electrocatalytic reactions and the resulting electrokinetic flows that drive particle motion.

Supplementary Figure 2: Particle tracking and data analysis. Example of particle tracking methodology for a single particle with number of fins n = 3, fin asymmetry c = 2, size a = 5.77 µm, and fin length b = 0.3 in 10 wt % H 2 O 2 . (a) Measured position of a single particle over time. The coordinates of two points on the particle were tracked for 5 seconds at 6 frames per second. The inset shows an image of the particle taken by optical microscopy with the position of each tracked area highlighted. The plot shows the x-coordinate of both tracked areas. (b) Angular orientation vs. time for the same particle. Orientation was determined from the x-and y-coordinates of the two tracked areas on the particle. Data was fit to a line by the least squares method in order to determine the angular velocity of the particle.

Supplementary Note 1: Non-dimensionalization of the Model
The governing equations and boundary conditions outlined in the Methods can be non-dimensionalized using the following characteristic scales: [x] = λ reaction-diffusion length (2) [φ] = k B T /e thermal potential (3) [u] = U 0 characteristic velocity (4) Here, the bulk ion concentration is related to the peroxide concentration as The reaction-diffusion length λ is defined as The characteristic diffusivity D ± is defined as the average diffusivity of the dominant ions H + and Finally, the characteristic velocity U 0 is defined as In the subsequent sections, we reformulate the governing equations and boundary conditions using these scales and discuss the various dimensionless parameters that emerge. For simplicity of notation, we use the same symbols for dimensional and dimensionless variables alike; the distinction should be clear from the context.
Species Transport. The dimensionless concentration of species i is governed by the conservation equation where Pe = U 0 λ/D ± is the Peclet number. We anticipate that the Péclet number will be very small such that convective species transport can be neglected. The boundary conditions far from the particle become where the dimensionless peroxide concentration is γ = (C H 2 O 2 /K p ) 1/2 1. The dimensionless reaction rates are The dimensionless current densities at the particle surface are given by where the parameter µ characterizes the relative rates of reaction and diffusion at the particle surface This equation defines the reference potential Φ 0 from which the Stern layer voltage is measured, This potential difference ∆Φ must be determined self-consistently to ensure the condition of zero net current The boundary conditions at the particle surface become Electrostatics. The dimensionless potential φ is governed by where β = λ D /λ is the ratio between the Debye screening length, λ 2 D = εk B T /2e 2 C ∞ ± , and the reaction-diffusion length λ. Far from the particle, the potential decays to zero; at the particle surface, we assume a constant surface potential Φ s , Hydrodynamics. The dimensionless velocity and pressure are governed by the Stokes equations including an electric body force f e ∇ · u = 0 (22) The boundary conditions at the particle surface S and far from the particle are u(x) = 0 for |x| → ∞.
Finally, the particle velocities U and Ω are determined by the conditions of no net force and torque on the particle, accounting for contributions due both to hydrodynamic and electric stresses, Here, σ = −pδ + (∇u + ∇u T ) is the hydrodynamic stress tensor and σ e = EE − 1 2 E 2 δ is the Maxwell stress tensor. This condition follows from the application of the divergence theorem to the Stokes equation (23), noting that f e = ∇ · σ e .

Symbol
Value We hypothesize that the reaction parameter µ is sufficiently large that the electrochemical reactions are limited by the transport of HO − 2 to the platinum surface (µ 1). Even in this limit, the reaction does little to reduce the local concentration of peroxide, which is present in great excess (γ 1). We will therefore assume that the concentration of peroxide is everywhere constant and equal to its bulk value C H 2 O 2 = γ. Under these conditions, the model has a well defined solution for the 1D concentration profiles in the vicinity of a planar platinum surface (Fig. 4b). The rate-limiting species HO − 2 is continuously generated in the bulk and diffuses to the surface at a rate Concentration variations extend over the reaction-diffusion length λ before approaching the equilibrium bulk values. This prediction compares favorably with electrochemical measurements of the partial current densities of anodic oxidation and cathodic reduction of hydrogen peroxide on platinum for C 0 supporting electrolyte at pH = 5.9. 6 At the mixed potential, the measured currents are ca. 1 A/m 2 , which corresponds to a peroxide consumption rate of 1 × 10 −5 mol/m 2 ·s. The above expression predicts a comparable rate of 0.93 × 10 −5 mol/m 2 ·s.
The rate of H 2 O 2 decomposition by this electrochemical mechanism represents a small fraction of the overall decomposition rate, which includes both electrochemical and non-electrochemical processes. The overall rate of peroxide decomposition by platinum-polystyrene Janus particles (2 µm diameter) in 10 wt% H 2 O 2 (3.0 M) was measured previously to be (8 ± 4) × 10 10 molecules per second per particle (0.02 mol/m 2 ·s). 7 The model above predicts a rate of 5.1 × 10 −5 mol/m 2 ·s (ca. 400 times smaller).
Using the above results, we can evaluate the approximation that the peroxide concentration is everywhere constant. The overall rate of peroxide decomposition is of order 0.02 mol/m 2 ·s, which must be balanced by the flux of peroxide to the surface. The diffusive flux of peroxide to a particle of size a is approximately is the difference in the peroxide concentration between the surface and the bulk. Equating these two estimates implies a concentration difference of ∆C H 2 O 2 ∼ 10 mM, which is much smaller than the bulk concentration of C H 2 O 2 = 3 M.

Supplementary Note 3: Equilibrium Dissociation of H 2 O and H 2 O 2 .
Far from the particle, the concentrations of H + , HO − 2 , and OHare determined by the equilibrium dissociation of H 2 O and H 2 O 2 . The conditions for dissociation equilibrium are given by where K w = 1.82 × 10 −16 M and K p = 2.40 × 10 −12 M are the dissociation coefficients for water and peroxide, respectively. Additionally, we have the following constraints due to species conservation and electroneutrality, denote the analytical concentrations of water and peroxide, respectively. For a given peroxide concentration, Supplementary Equations (30)-(33) can be solved to determine the concentrations of the respective dissociation products.
Under the relevant conditions (C 0 , the concentrations of water and peroxide are nearly identical to their analytical values-that is, With this approximation, the ion concentrations can be expressed as where is dimensionless parameter measuring the amount of ions due to water dissociation relative to that of peroxide dissociation. In experiment, this parameter is small ω ∼ 10 −3 1; below, we assume that ω → 0. Equilibrium ion concentrations for the self-ionization of hydrogen peroxide in water at 25 • C. Here, the peroxide concentration is computed from the peroxide weight fraction using the reported densities of aqueous peroxide solutions. 1

Supplementary Note 4: Rotation of Twisted Star Particles.
To estimate the angular velocity Ω = Ωe z of a twisted star particle, we decompose the Stokes flow u around the active particle into two simpler flows denoted u and u . In the first, the particle rotates with velocity Ω = Ωe z in the absence of any body force f = 0. In the second, a body force f = f e drives fluid flow around a stationary particle with Ω = 0. Owing to the linearity of the Stokes equations, these auxiliary flows can be superimposed to obtain the desired flow, u = u +u . Each of these flows is associated with a net torque on the particle, such that the net torque on the active particle can be written as L = L + L . Owing the symmetry of the particle, all torques are oriented in the z-direction-for example, L = Le z . According to equation (27), the net torque on the active particle is zero, L = 0. Moreover, the torque on the rotating particle is linearly related to its velocity as L = −RΩ, where R > 0 is a resistance coefficient. By estimating the torque on the stationary particle L , we can therefore determine the velocity Ω as Below we discuss our estimates of the resistance R and the torque L .
Resistance, R. In an unbounded fluid, the rotation of a disk-like particle of radius a can be approximated by that of a oblate ellipsoid of equal radius and vanishing thickness, R = 32 3 ηa 3 . In experiment, however, the particles rotate at some height h above a planar substrate. For small separations h a, the resistance to rotation is due largely to shear forces acting on the particle's bottom surface where b is the dimensionless fin length. For b 1, we obtain the approximation used in the main text-namely, R = πηa 4 /2h.
Torque, L . For thin disk-like particles rotating at a small height h above a planar substrate, the torque (37) can be approximated by neglecting the particle edges and integrating over the particle's top surface S t In this approximation, the contribution due to the electric stress is zero for particles with uniform surface potentials, for which n · σ e = 1 2 E 2 n. As noted in the main text, we do not compute the tangent components of the hydrodynamic stress directly. Instead, we compute the radial component of the surface stress σ rz on a thin circular disk and then map them onto the surface of the twisted star particles. This heuristic approach is discussed in more detail in the following section.  Table 1. Here, lengths are scaled by the reaction diffusion length λ = 140 nm, current by eD ± C ∞ ± /λ = 10 A/m 2 , and velocity by U 0 = C ∞ ± k B T λ/η = 1.0 mm/s. These parameters correspond to the experimental conditions of Howse et al. 8 for a polystyrene sphere of radius a = 0.81 µm coated on one hemisphere with a 5.5 nm thick layer of platinum. The particle is immersed in an aqueous solution of 10 wt% H 2 O 2 (C H 2 O 2 = 3.0 M or 7.4 vol%). The polystyrene hemisphere is assumed to be chemically inert with a zeta potential equal to that of the platinum hemisphere. The particle moves toward its polystyrene hemisphere with a computed speed of 2.2 µm/s. For comparison, Figure 4a of Howse et al. shows an experimental propulsion velocity of ca. 3 µm/s. 8 We emphasize that the model has no free parameters; however, several parameters are uncertain. Notably, an increase in the particle zeta potential could help to reconcile the 40% difference between experiment and theory. Further work is required to investigate the impact of different model parameters on the predicted propulsion velocity.

Supplementary Note 5: Conformal Mapping and the Role of Particle Shape.
Here, we consider how both the particle shape and the stress distribution influence the magnitude and direction of the torque L that drives particle rotation. For thin particles (δ a), the torque L z is approximated as where r = xe x + ye y is the position vector from the center of the particle face, and S t is the planar surface on the top side of the particle. To compute the integral, we parameterize the position vector as r(ρ, s) where ρ and s are orthogonal curvilinear coordinates (Supplementary Figure 10). The curve ρ = 1 corresponds to the perimeter of the particle, and ρ = 0 to the particle origin. We will assume that the stress depends only on ρ and is oriented along the negative ρ-directionthat is, e z · σ (r) = −σ (ρ)e ρ with σ (ρ) > 0. This heuristic approximation is exact for the circular disk, for which the coordinates (ρ, s) are identical to the polar coordinates (r, θ). For a weakly twisted star (small b and c), we hypothesize that this approach correctly predicts the direction and the characteristic magnitude of the torque on the particle. We note, however, that the magnitude and/or direction of the force density could vary with position along the perimeter s. Such variations are likely to be small when the reaction-diffusion length λ is much smaller than the radius of curvature of particle perimeter.
As detailed below, we find that the torque may be positive or negative depending on the stress distribution σ (ρ). When the stress is concentrated within a thin region near the edge of the disk, the torque on an R particle is positive, resulting in CCW rotation. When the stress is distributed across a thicker region around the particle parameter, the torque is negative, resulting in CW rotation (as in experiment). These qualitative trends agree with previous studies of similar particles rotating via asymmetric acoustic streaming (see below). 9 Orthogonal Curvilinear Coordinates (ρ, s) The orthogonal coordinates ρ(r) and s(r) satisfy the Cauchy-Riemann equations ∂ρ ∂x = ∂s ∂y and ∂s ∂x = − ∂ρ ∂y , which ensures that the unit vectors e ρ and e s in the ρ and s directions are everywhere orthogonal. In particular, the unit vectors e and scale factors h are defined in the usual manner as By construction, the unit vectors satisfy the following conditions to create a right-handed orthogonal coordinate system e ρ · e s = 0 and e ρ × e s = e z .
To construct the coordinate ρ(r), we first solve for an auxiliary function v(r) that satisfies the Poisson equation on a two-dimensional domain bounded by the particle contour C, Supplementary Figure 10: Orthogonal curvilinear coordinate system (ρ, s) for a twisted star particle with c = 1, n = 3, and b = 0.3. The particle face corresponds to the region 0 < ρ < 1 and 0 < s < 2π.
The function v(r) is set to zero along the perimeter of the particle, and a point sink is positioned at its origin. For a circular domain of unit radius, this problem can be solved analytically to obtain v(r) = ln r for 0 < r ≤ 1.
Note that the radial coordinate is given by r = exp(v). For the twisted star domain, the solution v(r) can be computed numerically to determine the generalized radial coordinate ρ(r) = exp(v(r)), which is zero at the particle origin and one along its perimeter. Finally, the coordinate s(r) can be computed from ρ(r) by numerical integration of the Cauchy-Riemann equations (42).

Thick Regions
When the stress is distributed across the whole particle face, the predicted direction of rotation agrees with the experimental observations. In particular we consider the scenario in which the stress is directed radially inward in the negative e ρ direction with a constant magnitude across the particle face-that is, e z · σ = −σ o e ρ for 0 < ρ < 1. The coordinate ρ(r) is computed numerically on a twisted star domain using Mathematica's finite element solver, and the unit vector e ρ is computed using Supplementary Equation (44). The torque on the particle is then computed by numerical evaluation of the integral (41) in Mathematica. Supplementary Figure 11 shows the computed torque as a function of the parameters c and n characterizing the shape of the twisted star particles. The qualitative trends are consistent with the experimental observations. The torque (which is proportional to the angular velocity) is negative for positive values of the asymmetry parameter c, increases in magnitude with increasing asymmetry, and is largely insensitive to the number of fins n.

Thin Regions
When the stress is concentrated near the particle perimeter, the direction of rotation reverses. Here, we consider the scenario in which the stress is given by e z · σ = −σ o e ρ for ρ o < ρ < 1, where Supplementary Figure 11: Computed torque on a twisted star particle subject to a uniform surface stress in the negative e ρ direction, e z · σ = −σ o e ρ . The torque is plotted as a function of the asymmetry parameter c for different numbers of fins n with b = 0.3. The inset shows the torque density (scaled by σ o a) on a particle with c = 1 and n = 3.
the position ρ o (s) corresponds to a small distance λ from the particle perimeter. In this case, the integral (41) for the torque can be expressed as where the limit of integration ρ o (s) is defined implicitly by the condition We focus on the limit where the length λ is much smaller than the particle radius a and develop a perturbation theory to compute the torque at successive orders in λ.
Perturbation Theory. The torque is expanded in a Taylor series in λ about λ = 0 As detailed below, the first order contribution is identically zero; we aim to compute the second order contribution to determine the direction of rotation. To develop this approximation, we will require additional series expansions in the coordinate ρ about ρ = 1 where ρ = ρ − 1. Here, the zeroth order contribution is simply the particle perimeter introduced in the main text. The first order contribution is the unit vector normal to the particle perimeter. These and other higher order contributions are detailed below. From the expansion (51), we can derive analogous expansions for the coordinate vectors a, scale factors h, and unit vectors e. Substituting these expansions into Supplementary Equation (48) for the torque and integrating with respect to ρ, we obtain where ρ o = ρ o − 1. To convert this expansion in ρ o into the desired expansion in λ, we make use of Supplementary Equation (49) to obtain This expansion can be inverted to obtain a series expansion for ρ o in powers of λ Curvilinear Coordinates. We develop a local approximation for the coordinate system (ρ, s) using a series expansion about a point (x o , y o ) on the particle contour at which ρ = 1 and s = s o . We start from the general Taylor expansions .
Results. For the parametric shapes used here, the first order contribution is identically zero Physically, this result implies that the application of a uniform force density along the particle perimeter directed normal to the perimeter results in zero net torque on the particle. At second order in λ, the torque has four possible contributions; however, numerical analysis suggests that only one is nonzero Supplementary Figure 12 plots this leading order contribution to the torque as a function of the asymmetry parameter c for different numbers of fins n. Notably, the direction of rotation is opposite of that for the case of thick regions (Supplementary Figure 11) with R particles (c > 0) rotating CCW (Ω > 0) and S particles (c < 0) rotating CW (Ω < 0). λ Supplementary Figure 12: Computed torque on a twisted star particle subject to a surface stress in the negative e ρ direction within a narrow region of thickness λ along the particle perimeter. The torque is plotted as a function of the asymmetry parameter c for different numbers of fins n with b = 0.3. The inset shows the torque density (scaled by σ o a) on a particle with c = 1, n = 3, and λ = 0.05a.

Discussion
To drive particle rotation, a surface force density of the form e z · σ = σ (ρ)e ρ must be distributed over a region of finite thickness. When this region is located along the particle perimeter, the twisted star particles with c > 0 are predicted to rotate in the CCW direction (Supplementary Figure 12). This behavior agrees with previous observations of particle rotation driven by acoustic streaming at high frequencies. 9 In that system, the time-averaged force density on the particle was concentrated within a thin region near the particle perimeter and directed inward to the particle center. Interestingly, hydrodynamic simulations at lower frequencies predict that the particle rotation should reverse direction as the thickness of the viscous boundary layer increases. This prediction is consistent with the above result for thick regions (Supplementary Figure 11). For the Pt-spinners, the stress distribution lies somewhere between these limiting scenarios. Supplementary Figure 13 shows the stress computed numerically for the circular disk. Although the stress is largest near the perimeter of the disk, there are significant contributions across the entire face. Consequently, when the torque is computed using this distribution (as in Figure 4 of the main text), the direction of rotation agrees with that of the uniform stress.  Figure 13: Computed stress e z ·σ = σ zr (r)e r on a circular disk of radius a = 41λ and thickness δ = 0.72λ as a function of the radial position r (scaled by the reaction diffusion length λ). The inset shows the torque density (scaled by C ∞ ± k B T a) on a twisted star with b = 0.3, c = 1, and n = 3 obtained by mapping the stress from the circular disk as σ zρ (ρ)e ρ = σ zr (r)e r . The rounded edge of the disk of width δ/2 is omitted in the mapping.