Colloidal bubble propulsion mediated through viscous flows

Abstract

Bubble-propelled catalytic colloids stand out as a uniquely efficient design for artificial controllable micromachines, but so far lack a general theoretical framework that explains the physics of their propulsion. Here we develop a combined diffusive and hydrodynamic theory of bubble growth near a spherical catalytic colloid, that allows us to explain the underlying mechanism and the influence of environmental and material parameters. We identify two dimensionless groups, related to colloidal activity and the background fluid, that govern a saddle-node bifurcation of the bubble growth dynamics, and calculate the generated flows analytically for both slip and no slip boundary conditions on the bubble. We finish with a discussion of the assumptions and predictions of our model in the context of existing experimental results, and conclude that some of the observed behaviour, notably the ratchet-like gait, may stem from peculiarities of the experimental setup rather than fundamental physics of the propulsive mechanism.

---

I. Introduction

Since their introduction in 2004,¹ catalytic micromotors have been subject of intense theoretical and experimental investigation.^2–4 This interest is driven in part by their potential in a wide range of applications, including biomedical uses such as drug delivery⁵ or stem cell transplantation.⁶ Initially, investigations focused mostly on theoretical^7–9 and material-science aspects.^10,11 Recently, practical uses of these active motors have been achieved,¹² facilitated in part by the introduction of bio-compatible designs.¹³ In many cases, these microrobots are immersed in a fluid medium, and exploit (or are subject to) the particular physics of hydrodynamics at such small scales. To this end, several physical mechanisms to generate propulsion have been proposed, including external actuation with magnetic fields¹⁴ and ultrasound,¹⁵ as well as autonomous swimming through the exploitation of phoretic effects,^16,17 and bubble propulsion.^18–21 While phoretic propulsion relies on a fluid flow that is driven by an osmotic imbalance of ions produced on the swimmer's surface in an asymmetric fashion, bubble propulsion is a variant of this mechanism in which the produced fluxes are so large that locally a supersaturation is achieved, which makes the nucleation of a gas bubble energetically favourable. Rather than the weak coupling between chemical gradients and fluid mechanics,¹⁶ it is then a direct mechanical coupling with the bubble growth dynamics that propels the swimmer. As a result, the achieved propulsion speeds are much higher, on the order of mm s⁻¹, rather than μm s⁻¹.⁴

Perhaps due to these empirical observation of fast motion, theoretical work investigating bubble-propelled swimmers has focused for the most part on inertial models, such as momentum transfer from the detachment of nanobubbles,¹⁸ and inertial bubble collapse.¹⁹ Attempts have been made to rationalise bubble propulsion dynamics by comparing with the phenomenology of bubble growth as simulated with lattice-Boltzmann methods and investigated experimentally,^22,23 but since these had been performed in a regime where the Reynolds number is at least of order unity, it remains doubtful whether they capture the correct physics. Indeed, in order to compensate for the discrepancy of scales, any quantitative comparisons have required the adjustment of empirical coefficients for the forces that growing bubbles exert on boundaries, normally of order unity, to several orders of magnitude larger.^19,20 Other models for bubble-propulsion have given significant attention to surface tension, mass and momentum conservation, but did not incorporate fluid mechanical modelling of the interaction between bubble and colloid.^24,25

While flow physics were incorporated numerically in the specific case of conical microrockets,^21,26,27 no theoretical work has been devoted towards gaining a fundamental understanding of the propulsion mechanism for the classical design of a spherical colloidal catalyst powered by a collection of nucleated bubbles. Modelling needs to include (1) the physics of bubble growth itself, which due to the small scales exhibit peculiarities that are different from macroscopic dynamics^28,29 and (2) the transport within the realm of inertia-less hydrodynamics, in which propulsion is not achieved by momentum exchange with the surrounding fluid, but through non-reciprocal changes in the geometry,^30,31 here due to bubble growth.

The canonical experimental setup is illustrated Fig. 1, adapted from ref. 19, including (a) the production of the microswimmer, (b) an image of the swimmer featuring the smooth Pt-covered reactive surface, (c) a time-lapse of the bubble growth dynamics, and (d) the resulting trajectory. Bubbles are nucleated on the surface of the spherical colloid and propulsion occurs in a two-step ratchet-like periodic motion where slow pushing due to diffusive bubble growth is followed by fast retraction due to sudden bubble disappearance. In this work, we present a theoretical study to address this two-step motion. We first derive a coupled model of diffusive bubble growth and inertialess hydrodynamics in terms of experimentally relevant parameters in Section II. Here we also include a discussion of the model assumptions and their experimental relevance. We then solve this model, first for the instantaneous diffusive growth dynamics and instantaneous coupling with fluid dynamics, followed by their combined time-dependent evolution in Section III. We comment briefly on the conservation of momentum during inertial bubble collapse in Section III D and make a prediction for a new type of bubble-propelled swimmer that could be experimentally realised. We conclude with a summary in Section IV.

Fig. 1 Schematic of the experimental setup by ref. 19. (a) Manufacture of beads, (b) scanning electron microscope image of the swimmer. (c) Snapshots of a microbead showing the bubble growth and burst processes and the bead motion behavior. (d) The trajectory of the bead. The red arrows denote the direction of the trajectories of the bead after bubble burst, and the green arrows represent the direction of the trajectories of the bead during bubble growth. Reprinted with permission from Manjare, Yang and Zhao, Phys. Rev. Lett., 2012, 109, 128305.¹⁹ Copyright 2012 by the American Physical Society.

II. Mathematical modelling of bubble growth and propulsion

A. Motivation

We begin by building some intuition for the physics of the growth process based on the typical length and time scales involved. The colloidal particle and the bubble are small, with a relevant length scale , while the experimentally observed time scale of the bubble growth cycle is about t ∼ 0.1 s.¹⁹ Assuming a surface tension γ ≈ 7.2 × 10⁻² N m⁻¹ as for a clear air–water interface, a mass density for the fluid of ρ_l ≈ 1 × 10³ kg m⁻³ and a dynamic fluid viscosity of μ ≈ 1 × 10⁻³ Pa s, we can compute the dimensionless Reynolds and Capillary numbers asFrom this we deduce that inertial effects and the contribution of the dynamic pressure to the total pressure in the bubble are both negligible. This implies that the bubble does remain spherical throughout its growth (as is observed in experiments). We can thus write at all time that the bubble pressure p_bubble is given by the Laplace law²⁹where p^∞ is the hydrostatic background pressure and R_B is the radius of the bubble. This is the inertia-less version of the general Rayleigh–Plesset equation that is canonically used for the description of bubble growth. In this limit the bubble growth rate Ṙ_B does not enter the equation. As a result, the bubble size varies quasi-steadily with the bubble pressure.

The pressure in the bubble p_bubble is in turn related to the (molar) gas density ρ_gvia the ideal gas law as

where is the universal gas constant and T is the absolute temperature, which we assume to be uniform and constant in time. To determine the density ρ_g(t) as a function of time, we thus need to solve for the molar concentration field c of dissolved gas surrounding the system, and then compute the mass flux into the bubble.

The two fundamental physical processes at work in the system are therefore (i) the production and transport of gas into the bubble, leading to its growth, and (ii) hydrodynamic interactions between the bubble and the colloidal particle that convert bubble growth into propulsion. In general these two problems are coupled mathematically, since the flux-dependent bubble growth rate determines the fluid flows, and these in turn can affect the transport of dissolved gas through advection. However, in the limit where gas transport is dominated by diffusion the latter effect is negligible, so there is only a one-way coupling: a quasi-steady problem may be solved first to find the instantaneous bubble growth rate, and its solution can then be fed into the equations of fluid mechanics to determine the propulsion velocity. In our discussion of parameter estimates in Section II E, we show that this condition is adequately satisfied for the canonical decomposition of hydrogen peroxide in which the surface of the colloid acts as a catalyst and the produced gas is oxygen (2H₂O₂ → O₂ + 2H₂O) as used in many experimental studies.^2,4,19

B. Assumptions

In the course of developing and analysing this model we make a number of assumptions, which we briefly summarise here.

Since the most commonly used experimental swimmers feature a Pt-coated surface acting as a catalyst for the chemical reaction,^4,19 we assume that the reactive surface is smooth and not porous. While examples of nanoporous colloids exist,^32,33 these have not to our knowledge been shown to self-propel, or designed with the goal of achieving self-propulsion. While the raw silica beads used in the production of phoretic swimmers do appear to exhibit significant surface roughness in SEM images, this roughness is not visible after Pt-coating (see Fig. 1(b)). It is conceivable that the nature of the coating process tends to smooth out such irregularities, and we hence ignore them in this work.

Furthermore, previous theoretical models of (phoretic) swimming have occasionally included explicitly a thin ionic interaction layer of width ∼1 nm above the reactive surface.^16,34,35 Since this layer is negligibly thin compared to the colloid and bubble radii, we also ignore it here. While it may be important locally at the contact point, this effect would amount only to a higher-order correction to the global fluxes computed in this study. For the reaction itself, we assume in the main text that a uniform constant flux of solute is produced on the colloid, and additionally discuss quantitatively the higher-order correction due to considering only a constant reaction rate in Appendix B.3.

For mathematical simplicity, we further assume that there is only a single bubble present in the problem, and we consider only the dynamics post-nucleation. We note indeed that experimentally, while there are occasionally multiple bubbles forming, there is usually one that emerges to dominate in size.¹⁹ At early stages, this process is likely driven by Ostwald ripening.³⁶ At later stages, local perturbations of the concentration field due to the presence of small bubbles may lead to short-term deviations on the time scale of dissolution. Modelling these additional bubbles mathematically is challenging, and it is unclear whether they enhance or inhibit propulsion. Their contribution would thus deserve an analysis that is beyond the scope of the present paper. Additionally, we ignore any variation in surface tension on the bubble that may arise, for instance, from a spatial heterogeneities in the chemical reaction.

We also assume that the colloidal particle is spherical, and that it remains in contact with the bubble throughout the growth process. In reality, there may be a partial or complete wetting layer. In Appendix C.4 we demonstrate that a thin gap only affects the hydrodynamics at higher order than the global flows. The same is true for the concentration of reaction product in the fluid, since there is no singularity at the contact point and the solution varies continuously when a gap is introduced. At length scales of less than ∼1 nm, the Rayleigh–Plesset eqn (2) may need to be modified to include effects of adsorption,³⁷ but these effects would again only be relevant very locally at the contact point and thus a higher-order correction to the findings of this study. A discussion of partial wetting and its consequences for bubble detachment is offered in Appendix D.

Finally we assume, also for mathematical simplicity, that the colloid is radially symmetric and all parts of its surface equally participate in the production of solute. As we will show, the gradients of the concentration field surrounding the bubble determine the bubble growth rate, and so the concentration boundary condition on the far side of the colloid can reasonably be expected to be of limited importance. However, since the net flux of solute produced by the colloid will appear as an important parameter, , this parameter will need to be rescaled when a particle with lower active surface coverage is considered.

C. Derivation of the model equations

We use tangent sphere coordinates {μ, ν, θ} as defined in Appendix A to conveniently describe the geometry, with surfaces of constant ν defining spheres that touch at the origin. The bubble is considered to be the sphere defined by ν = ν_B := R_B⁻¹, and the colloidal particle of radius R_c to be the sphere defined by ν = −ν_c := −R_C⁻¹. A sketch of the geometrical setup is provided in Fig. 2. Note that in the figure we scale lengths by the colloid radius R_c, in line with the dimensionless formulation of the problem that will be introduced in Section II D.

Fig. 2 Growth of a spherical bubble attached to a spherical colloidal particle: sketch of the geometry in coordinates scaled by the colloid radius R_c. The bubble (light sphere) has radius R_B(t) and grows at a rate Ṙ_B while the colloidal particle (dark sphere) has unit radius. Both bubble and colloid are assumed spherical and touch at the origin of a cylindrical coordinate system {ρ, z}, which is equivalent to the point at infinity in tangent-sphere coordinates {μ, ν}. The whole system translates with a velocity U along the axis of symmetry (μ = 0), defined positive when the bubble pushes the colloid.

The colloidal particle acts as a catalyst for a chemical reaction. The concentration of the reaction product, c, satisfies the quasi-steady diffusion equation as

D∇²c = 0, (4)

where D is the molecular diffusivity. The problem is axisymmetric and the solution c(μ,ν) is therefore independent of the azimuthal angle θ (Fig. 2). The concentration on the bubble surface is set by Henry's law,³⁸ which states that concentration is proportional to the partial pressure of the corresponding gas inside the bubble, i.e. p_bubble. Assuming that the bubble contains only the reaction product, eqn (2) leads to

c = K_H⁻¹ (p^∞ + 2γν_B), at ν = ν_B, (5)

where K_H is the constant of Henry's law, which depends on the dissolved gas.

As discussed in Section II B, we assume that the entire surface participates equally in the reaction. We model this with zeroth-order kinetics (0ok), which assumes a constant catalytic flux per unit area that is independent of other factors. This results in a Neumann boundary condition for the concentration product of the form

A more realistic description is first-order kinetics (1ok), which involves a second chemical solute (the ‘fuel’) of concentration c, that is converted into a product (concentration c) on the colloid surface at a constant rate. We find that this does not significantly change the c field and only leads to an effective rescaling of global catalytic activity, so a discussion of this boundary condition is relegated to Appendix B.3.

Finally, the boundary condition far from the colloidal particle is

where c^∞ is the background concentration of reaction product. Normally, the value of c^∞ is determined by applying Henry's law at the interface with the surrounding atmosphere.

The instantaneous concentration field for a given bubble radius is determined by the solution to eqn (4) along with the boundary conditions in eqn (5)–(7). The instantaneous growth rate of the bubble is then determined by mass conservation,

which states that the rate of change of the bubble mass (left) is equal to diffusive mass flux into the bubble (right).

Meanwhile, since the Reynolds number is very small, the fluid dynamics surrounding the system are governed by the linear and quasi-steady incompressible Stokes equations,³⁰

∇p = η∇²u, ∇·u = 0, (9)

where η is the dynamic fluid viscosity. The Stokes problem described above is linear with respect to its driving kinematics, the bubble growth velocity Ṙ_B; the resulting instantaneous translation velocity of the bubble-colloid pair, U(t), is linear Ṙ_B, with a constant of proportionality that depends only on the geometry, i.e. the ratio of bubble and colloid radii:The function f can be calculated by imposing that the colloid-bubble system remains force-free.³⁰

As boundary conditions for the fluid flow, we assume an undisturbed fluid in the bulk (i.e.u → 0 as ) and a no slip condition on the colloidal particle. On the bubble, we solve for both cases of a no stress (i.e. free surface) and a no slip condition, the latter to account for a possible rigidification of the bubble due to surfactants. Combining eqn (8) and (10) with the solution to eqn (4) for the instantaneous concentration field, this determines the velocity U(t) of the bubble-colloid pair at each point in time.

D. Dimensionless formulation and sketch of the solution method

We now reformulate the problem in dimensionless form by introducing the following scalings for lengths, concentration and time respectively,that is we scale lengths by the colloid radius, concentrations by the typical value on a bubble the same size as the colloid, and time by a natural scale for the diffusive growth described by eqn (8). As we show in Appendix B, the general solution to the diffusion equation, eqn (4), that decays far from the colloid, eqn (7), may be written exactly aswhere c̃^∞ = c^∞K_HR_c/2γ and the functions A(s) and B(s) are such that the integral converges in the range μ ≥ 0, −1 ≤ ν ≤ ν_B. Here J₀ denotes the zeroth order Bessel function of the first kind. For zeroth-order kinetics, the two remaining boundary conditions eqn (5) and (6) simplify to

c − c̃^∞ = ξ + ν_B at ν = ν_B, (13)

where we define the dimensionless parameters ξ and asPhysically, ξ is a dimensionless measure of the saturation of the background fluid, scaled by capillary effects (=0 when saturated, >0 otherwise). We may therefore interpret ξ as a measure of how ‘soluble’ gas is in the background fluid relative to the bubble. Intuitively, a larger value of ξ is therefore conducive to bubble dissolution, and hence limits bubble growth. Meanwhile is interpreted as a scaled catalytic flux out of the colloid, and it is expected that a large value of will be conducive to bubble nucleation and growth.

Applying the boundary conditions eqn (13) and (14) to eqn (12) leads to a regular second-order ODE boundary value problem that can be solved numerically for given values of the parameters and ξ. The solution for A may then be substituted into the mass conservation equation, eqn (8), which simplifies to

where we have introduced the dimensionless group ζ = ξ + c̃^∞ = p^∞R_c/2γ. As a corollary we have the condition for a steady state,which is satisfied when there is no net flux of gas in or out of the bubble. A detailed derivation of the equations for this diffusive problem and their solution is provided in Appendix B.

For the fluid mechanics we use a Stokes streamfunction Ψ with the ansatz

where Ψ_s ∝ Ṙ_B is a singular contribution that accounts for expansion of the bubble, and where the coefficients Ã(s) to D̃(s) are found by applying the hydrodynamic boundary conditions on bubble and colloid (they are thus implicitly a function of the instantaneous bubble radius R_B).

We additionally require that the bubble-colloid system remains force-free at all times, which may be shown to reduce to the condition (see Appendix C). Since the boundary conditions are linear in both U and Ṙ_B, the same is true for Ψ, and so we may split the full solution for the flow as the linear superposition of a ‘motility’ problem ∝U and a ‘growth’ problem ∝Ṙ_B. The function f in eqn (10) is then given by

Full details, including entirely analytical expressions for the functions B̃_mot(s) and B̃_gro(s), are provided in Appendix C.

E. Parameter estimates

In order to provide an experimentally relevant model, it is essential to identify the typical ranges and variability of the relevant dimensionless parameters in practical conditions. We take γ ≈ 7.2 × 10⁻² N m⁻¹ as for an air–water interface, and p^∞ ≈ 1.0 × 10⁵ Pa for atmospheric pressure. The gas constant is fixed as and at room temperature T ≈ 300 K approximately. For oxygen gas dissolved in water at 25 °C, D ≈ 2.4 × 10⁻⁹ m² s⁻¹ and K_H ≈ 7.7 × 10⁴ J mol⁻¹ are standard values.³⁹ From this, Henry's law calculates a background molarity c^∞ ≈ 2.7 × 10⁻¹ mol m⁻³, which is relatively large due to the percentage of oxygen in the atmosphere, which is well known to be 21% (= K_Hc^∞/p^∞).³⁹

For the activity , it is harder to find an estimate due to strong variability of conditions and reactions driving the propulsion in the relevant experiments. In the limit of zero Damköhler number (i.e. when first-order kinetics reduce to zeroth-order), is given by , which is the product of reaction rate k and bulk reaction fuel molarity c^∞_f. In the classical platinum-catalytic H₂O₂ decomposition we can assume a molarity c^∞_f ≈ 2.9 × 10² mol m⁻³ (corresponding to 1 wt% H₂O₂, a low estimate) and k ≈ 4.1 × 10⁻⁵ m s⁻¹.^40,41 Finally, we consider colloid radii R_c between 1 μm and 50 μm, as relevant to all published experiments to do (the colloidal size is the experimentally most variable parameter).

Combining all of these, we find the typical range for ξ is , while is in the range . Furthermore, the growth time scales as , which is consistent with experimental observations.^19,42 Since ξ scales as ∼R_c, and and t̃ scale as ∼R_c², larger values apply to larger colloids, and due to the quadratic dependence on R_c, the value of depends quite sensitively on the colloid size. As we shall see in Section III A, this is the fundamental reason why bubble propulsion is limited to sufficiently large colloids.

In the development of our model we assumed that the bubble grows in a quasi-steady fashion, i.e. we neglected advective terms in eqn (4). This is appropriate if the Péclet number Pe is small, i.e. if the time scale of diffusive transport is much shorter than the time scale of fluid flow. As we see from eqn (10), the time scale of fluid flow is linked to the time scale of bubble growth, which allows us to justify this assumption a posteriori by demonstrating that it is self-consistent with its prediction. Quantitatively this is the requirement that

where we have assumed diffusion on the scale of the colloid. With our parameter estimates, we obtain Pe ≈ 0.097, which is indeed small compared to unity, and even more so when diffusion on length scales smaller than the colloid is concerned. Thus the quasi-steady assumption is approximately valid and advective transport only leads to a small correction. It is important to note however that for gases other than oxygen, K_H can be substantially different. While for other common gases such as N₂ and H₂ we also have Pe ≪ 1, this is not the case e.g. for CO₂ (Pe ≈ 2.4) or NH₃ (Pe ≈ 750),²⁸ in which case it would be essential to solve the fully coupled nonlinear model for the transport of the reaction product.

III. Analysis of bubble growth and propulsion dynamics

A. Solute concentration fields and stationary points of the growth dynamics

We first solve eqn (4) for the instantaneous concentration field as a function of the dimensionless parameters , ξ and the bubble radius R_B, and identify the stationary points of the bubble size by means of the condition in eqn (17). The results are shown in Fig. 3. For small values of , the bubble always shrinks to zero. Above a critical activity, a saddle-node bifurcation leads to the existence of three regimes with two fixed points such that

bubble shrinks < R_min < bubble grows < R_max < bubble shrinks. (21)

At the critical point we have approximately that and R_B ∼ ξ^−0.5.

Fig. 3 Left: Stationary points of the bubble growth dynamics (solid lines) as a function of the parameters and ξ. Dashed and dash-dotted lines indicate theoretical estimates for R_max (stable) and R_min (unstable) respectively. Right: Concentration fields with c-contours in red, the colloid in black, and the bubble in white, for 4 points in space with ξ = 10.

The behaviour of this system can be understood from a theoretical point of view by considering the competing effects that lead to bubble growth and shrinkage. Small bubbles have high capillary pressure, which enforces a high concentration on their surface through Henry's law. Only if an even larger concentration around the bubble can be sustained through catalytic flux from the colloid can there be bubble growth. Hence, the lower fixed point is achieved when there is a balance between these effects, i.e.

For more details of this calculation, we refer to Appendix B.2. For a small bubble to grow and propel the colloid, R_min needs to be sufficiently small. Since and ξ ∼ R_c, this explains quantitatively why bubble propulsion is only observed for strongly catalytic colloids.¹⁹

In the opposite case, the bubble reaches its final size when influx from the colloid balances with outflux due to diffusive dissolution. This corresponds to the condition

where the origin of these scalings is also discussed in Appendix B.2.

In Fig. 3 we plot the numerically determined stationary bubble sizes (solid lines) against the theoretical predictions from eqn (22) and (23) (dashed lines) as a function of the dimensionless groups and ξ, and find that there is excellent agreement when . The disagreement becomes significant near the bifurcation point, when bubble and colloid are of roughly the same size and geometric details become important.

B. Relationship between bubble growth rate and propulsion

Focusing next on propulsion, we now calculate the function f in eqn (10) for a no slip condition on the colloid, a flow decaying far away on the system, and two different boundary conditions on the bubble: a ‘rigid’ bubble with no slip as on the colloid, and a ‘shear-free’ bubble that supports no tangential shear stress. While the second condition is canonical for an air–water interface in most contexts, it is a well-known fact that surfactants and impurities can lead to a rigidification of these interfaces under experimental conditions.^28,29 Since it is not entirely clear which boundary condition applies for this system, and they differ slightly in their physical predictions, we include both in our discussion below.

The results are illustrated in Fig. 4, showing f = U/Ṙ_B as a function of R_B for no slip (blue solid line) and no shear stress (red dash-dotted line) boundary conditions. In both cases, the ratio of propulsion speed to bubble growth rate U/Ṙ_B increases monotonically with the bubble size, so for a fixed growth rate, large bubbles are more efficient than small ones at pushing the colloid. While the mathematical expressions for the drag and propulsion coefficients are very complicated, the asymptotic behaviour for large and small bubbles can be understood intuitively. As R_B → ∞ (ν_B → 0), both the drag and propulsion scale as R_B and we recover the drag coefficients for a single sphere in Stokes flow, i.e. 6πηR_B or 4πηR_B depending on the boundary condition (no slip or no shear). In this limit, the colloid's contribution to drag is negligible and it is just pushed outwards with velocity U = Ṙ_B by the growing bubble that remains stationary itself. Conversely, as R_B → 0 the drag asymptotes to that of an isolated colloid, 6πη, while the propulsion force decays as R_B⁻². This is line with expectations of the force exerted by a point source on a rigid sphere, which scales as Q/d with Q the source strength and d the distance to the colloid centre.⁴³ Interestingly, it makes no difference in this limit whether the boundary condition on the bubble is no slip or shear-free. The difference is most pronounced when bubble and colloid are of equal size, and because a rigid surface has more drag the propulsion is up to 50% more efficient when the action of surfactants is strong.

Fig. 4 Ratio of propulsion speed, U, to bubble growth rate, Ṙ_B, for no slip (blue solid line) and no shear boundary conditions (red dash-dotted line) as a function of bubble size.

C. Time-dependent propulsion dynamics

1. Evolution of a single growth cycle and comparison with experiments. We are now in a position to combine all of our results to find the colloid propulsion speed, U, and distance travelled, X, over one bubble growth cycle as a function of the dimensionless parameters and ξ, which we recall represent catalytic flux magnitude and solubility of gas in the background fluid relative to the bubble.

To learn about the generic behaviour of the system, we initialise the bubble at a size slightly above the unstable stationary point R_min, calculate the growth rate according to eqn (16) and the rate of displacement Ẋ = U from eqn (10), and integrate the system forward in time using first-order explicit Euler time-stepping. A typical example of the resulting dynamics is shown in Fig. 5, where we display the time-variation of R_B, Ṙ_B, X and Ẋ. We note that it follows from time-integration of eqn (10), that the displacement X is slaved to the instantaneous bubble size. Since small bubbles propel the colloid insignificantly, the bubble nucleation size has little influence over the final displacement. Hence, the final value of X is determined, to within a good approximation, only by the final bubble size and the hydrodynamic boundary condition on the bubble. Consequentially, the only variation in the dynamics shown in Fig. 5 for different values of the parameters are in the variation of the final bubble size (discussed previously in Section III A), and the time scales of growth.

Fig. 5 Example of propulsion over one bubble growth cycle for intermediate values of the parameters, and ξ = 10^0.5. Shown are the bubble radius R_B (purple) and travelled distance X (green), as well as their time-derivatives (dashed lines) against logarithmic time. Growth and propulsion rates increase rapidly on very short time scales, peak at intermediate ones, and decay as the bubble asymptotically reaches its stable equilibrium size at long times.

The logarithmic representation of the growth rate reveals three distinct regimes at short, intermediate and long times. At very short times, when the bubble size is close to R_min, velocities and displacements are small. However, since the deviation from an unstable fixed point generically follows exponential growth,⁴⁴ the dynamics transition rapidly to an intermediate regime in which the growth rate R_B peaks. In this regime the bubble typically crosses from being smaller than the colloid to being larger, which also maximises the colloid velocity U. At long times, the approach to the stable fixed point R_max is again exponential, and no significant propulsion occurs.

The striking prediction of our diffusive growth model is thus that parameter values that allow for the growth of very small bubbles simultaneously allow these bubbles to grow orders of magnitude larger than the colloid, unless the experiment is performed in a narrow region of the parameter space in which ξ is large and not much larger than its critical value. This was illustrated in Fig. 3 and rationalised with scaling arguments for the stationary points in eqn (22) and (23). Due to the scalings of these dimensionless parameters with the colloid radius, eqn (15), this is difficult to achieve in practice unless the experiment is carried out in a high-pressure environment. Indeed, both bubbles that are much larger than the colloid and a slow asymptotic approach to a stable equilibrium are typically not observed in experiments, even though the catalytic activity is clearly strong enough to allow for the nucleation of small bubbles.¹⁹

Aiming to address this apparent paradox, we considered both more realistic models for the catalytic reaction, and more detailed descriptions of the contact physics between bubble and colloid. Neither have brought to light any credible intrinsic mechanism that leads to these observations. Accounting for the depletion of reaction fuel with a non-zero Damköhler number leads to an effective decrease of the activity parameter , but not to a different scaling law for the maximal bubble radius (see Appendix B.3). Likewise, an analytical asymptotic solution for the pressure at the contact reveals that there is no divergent hydrodynamic force pulling the system apart, and both a scaling analysis and a simplified model with a non-singular contact demonstrate that bubble is too small for entrainment or buoyancy to be significant (see Appendix C.4 and D). We therefore conjecture that the premature disappearance of bubbles that was reported in ref. 19 is in fact due to an external perturbation, most likely the coalescence of the bubble with the air–water interface of the experimental setup. This idea is in fact supported by analogous observations in later work.⁴²

While a detailed model of bubble growth in confinement is beyond the scope of this study, in the context of micro-rockets it has been shown numerically and experimentally that confining bubble growth is beneficial in achieving efficient propulsion.^21,45 Interpreting the bubble disappearance at the air–water interface as a confinement effect, we note that this is another example of this principle, albeit indirectly by limiting ‘natural’ bubble growth to the point in which it is efficient in propelling the colloid forwards.

In order to effectively make a prediction for experimental systems we thus need introduce a phenomenological cut-off for the bubble size into our model. Since the colloid sets a natural scale in terms of confinement, we choose R_max = 1 for this; note that it could also be larger if the colloid is in a deep layer.¹⁹

2. Dependence of time-averaged propulsion on catalytic activity and background saturation. From a practical point of view, an interesting question concerns the propulsive efficiency of the system as a function of its parameters. As established in Section III B, the total distance travelled, X, depends only on the final bubble size. Fixing this at the cut-off R_B = 1 leads to a prediction of the colloid displacement over one growth cycle as

X_{rigid bubble}/R_C ≈ 0.25, X_{shear-free bubble}/R_C ≈ 0.18, (24)

independently of the parameters and ξ.

Another quantity of interest is the time scale on which this displacement is achieved, or equivalently, the average speed of the colloid. We define this average 〈U〉 by the distance X travelled divided by the time required to reach the cut off size from a small perturbation away from the unstable fixed point R_min. The dependence of this average velocity on the two relevant parameters ( and ξ) is displayed in Fig. 6. Due to the rapidness of the initial exponential growth phase, these results are robust towards the exact size of the initial perturbation. For reference, we recall from our parameter estimates in Section II E that U = 1 in dimensionless units corresponds to the physical velocity

which is on the order of 1 mm s⁻¹. The dependence on the catalytic activity of the colloid is approximately linear above a cut-off , below which the bubble does not grow (Fig. 6, left). Close to this cut-off there is a small deviation from linear behaviour as the slow asymptotic approach to R_max becomes significant. Unsurprisingly, a larger activity leads to faster propulsion. Within the scope of our quasi-static model this growth is unbounded – in practice fast velocities will eventually lead to non-negligible Péclet and Reynolds number that regularise this behaviour. Additionally, a stronger solubility of gas in the fluid ξ has a negative effect on propulsion (Fig. 6, right).

Fig. 6 Simulated average propulsion speed 〈U〉 over one bubble growth cycle for rigid boundary conditions on the bubble surface, as a function of flux (left) and background solubility ξ (right). Above a certain ξ-dependent threshold, 〈U〉 increases linearly with . Conversely, the propulsion speed decreases as the solubility in the bulk increases.

D. Conservation of momentum during inertial bubble collapse

Finally, we briefly address the second phase of the propulsion cycle, which is bubble collapse. As established in Section III C1, the bubble size is limited in practice by the depth of the fluid layer in which the colloid is immersed, and so it is likely that the coalescence with the top interface leads to its disappearance. This process is very fast, and so for a brief moment the dynamics become inertial, and the colloid retracts a short distance.¹⁹ While a numerical simulation of this process would require faithful modelling of the geometry and non-linear flows, and is beyond the scope of this paper, we would like to present a simple thought experiment based on conservation of linear momentum that suggests interesting dynamics might happen for a colloid that is lighter than water.

Specifically, if we call the retraction distance d, we can consider conservation of linear momentum before and after the bubble collapse to find d as a function of the bubble radius before collapse, and of the densities of the colloid and the surrounding fluid. Since there is no net momentum change, the centre of mass of the whole system remains approximately stationary. The colloid is displaced a distance d backwards, while fluid is displaced a distance 2R_B − d forwards (assuming no inflow from the sides). Denoting by ρ_c and ρ_l the densities of the colloid and the surrounding fluid respectively we then have that

This is illustrated by a sketch in Fig. 7 (top row). Of course, this simplified argument ignores many important factors such as directionality of flows, time-dependence and confinement effects, and should therefore only be understood as a rough estimate for d as a function of the density ratio ρ_c/ρ_l. Within these limitations, we see then that for a heavy colloid (i.e. one with ρ_c ≫ ρ_l) we have d ≪ R_B, and so the retraction distance is much smaller than the final bubble radius, which is itself on the same order as the propulsion distance. As a result, the net displacement is such that the bubble pushes the colloid over one growth cycle.

Fig. 7 Thought experiment on bubble retraction due to inertial collapse. Based on momentum considerations, the retraction distance d is shorter than the propulsion distance if the colloid is heavy. But if the colloid is lighter than water, a net pulling by the bubble might be achieved if inertial retraction outweighs pushing due to bubble growth.

Conversely, one could imagine a scenario where ρ_c < ρ_l, where the simple argument predicts that d > R_B. This would correspond to a floating colloid, and it is illustrated in Fig. 7 (bottom row). This suggests that a floating colloid may actually be net pulled by the bubble, since the inertial retraction due to collapse could outweigh non-inertial propulsion due to growth. Whether this is possible in practice remains to be seen, but we hope to see this hypothesis tested experimentally.

IV. Discussion

In this work we have developed a combined diffusive-hydrodynamic model for the propulsion of a catalytic colloid by a growing gas bubble. We have identified key non-dimensional parameters of the system and their relation to experimentally accessible quantities, in particular the colloid size. With this, we have gained new understanding of the conditions which allows the nucleation of small bubbles that mark the departure from purely phoretic modes of propulsion into a regime where bubble growth allows colloids to reach much higher speeds. In particular, we identified that conditions that allow small bubbles to grow simultaneously allows them to grow much larger than the colloid, unless an external event or forcing leads to their detachment or disappearance. Finally, we combine all of our results to find a prediction for the propulsion speed of the system as a function of the catalytic activity of the colloid, and the solubility of gas in the surrounding medium. We conclude by predictions for an as yet unrealised setup in which a light catalytic colloid can surf on water being pulled by the collapse of bubbles that it generates itself.

As discussed in Section II B, our results are obtained under a number of simplifying assumptions in an effort to analyse the fundamental mechanisms on the simplest possible system. We purposefully do not address the question of bubble nucleation, nor do we consider multiple bubbles due to the mathematical complexity of the fluid mechanics. Instead, we demonstrate how fundamental physical ingredients can be combined to give novel predictions to be tested quantitatively in experiments. More advanced modelling of the dynamics e.g. at the contact point will be of interest once the leading-order predictions obtained here have been validated. In particular, nanoporous colloids with potentially sophisticated bubble nucleation and contact line dynamics present an interesting future research direction.

Our results pave the way for further experimental and theoretical investigations of bubble propulsion. Since the focus has been on gaining an analytical understanding of the fundamental physics, our study has been mostly limited to the geometry of a bubble-colloid pair in isolation. Yet as discussed confinement plays a crucial role by clearing the bubble from the system. Further investigation may reveal the extent to which this process can be optimised, and how the substrate on which heavy colloids are typically located affects propulsive efficiency.

Appendix A: properties of tangent sphere coordinates

We define tangent sphere coordinates {μ, ν, θ} in terms of usual cylindrical coordinates {ρ, θ, z} asTheir ranges are μ ≥ 0, , and θ ∈ (0, 2π). They are essentially an ‘inversion’ of cylindrical coordinates. As such, μ and ν have units of inverse length, and the coordinate origin and point at infinity have been exchanged. Surfaces of constant ν define spheres with radius 1/|ν| and their centre located on the z-axis at z = 1/ν. Thus all surfaces of constant ν touch the plane z = 0 at the spatial origin, and the sign of ν determines whether they are located in the half space z > 0 (for ν positive), or z < 0. The surface ν = 0 corresponds to the plane z = 0. Similarly, surfaces of constant μ define toroids around the z-axis with circular cross-sectional radius 1/μ. The surface defined by μ = 0 is degenerate and corresponds to the z-axis. Finally, θ is the azimuthal angle as in cylindrical coordinates.

The scale factors are

and the basis vectors areHence and in particularregardless of the sign of ν.

Appendix B: detailed solution of the steady diffusive growth model

B.1 Zeroth-order kinetics

We assume non-dimensional scalings. Recalling from the main text, we aim to solve the problem

∇²c = 0, (B1)

subject to the boundary conditions

c − c̃^∞ = ξ + ν_B at ν = ν_B, (B2)

Defining G(μ, ν) asthe diffusion equation eqn (B1) reduces toThus a solution that satisfies the boundary condition at spatial infinity, eqn (B4), is determined by a solution for G that remains bounded as . Using separation of variables this yields G ∼ exp(±sν)J₀(sμ), where J₀ denotes the zeroth order Bessel function of the first kind (a second family of solutions, proportional to the divergent Y₀(sμ), is discarded). Hencefor some functions A(s) and B(s) that are such that the integral converges in the range μ ≥ 0, −1 ≤ ν ≤ ν_B.

To solve the steady diffusive problem, we need to find the functions A(s) and B(s). The boundary condition on the bubble, eqn (B2), gives

while the one on the colloid, eqn (B3), yieldswhere we have defined the functions

α(s) = e^−sA(s) − e^sB(s), β(s) = e^−sA(s) + e^sB(s). (B10)

The right-hand side of these relations is reminiscent of the Hankel transform which is defined asMoreover, andHowever, in order to apply these relations to our boundary condition, we need to eliminate μ² on the RHS of eqn (B9). To this end we define the operatorwhich satisfiesIn particular, this allows us to writewhere the second equality follows from two integrations by parts if we place the condition that α(0) is bounded and α → 0 as s → ∞. We may hence rewrite the boundary conditions to the diffusive problem asand apply a Hankel transform to find

(ν_B + ξ)e^−sν_B = α(s)sinh(s(1 + ν_B)) + β(s)cosh(s(1 + ν_B)), (B18)

which gives

β = (ν_B + ξ)e^−sν_B sech(s(1 + ν_B)) − tanh(s(1 + ν_B))α, (B21)

where primes denote derivatives with respect to s and the boundary conditions are α(0) bounded and α → 0 as s → ∞. The ODE is singular at s = 0, but this singularity is balanced by the forcing. Substituting a series solution of the formreadily shows that the correct lower boundary condition is .

Finally, we note that we can use eqn (A5) and (B12) to find the flux integral

which demonstrates that the flux out of the colloid (with ν < 0) is proportional to , while the flux into the bubble (with ν > 0) is proportional to . Any residual between these is emitted to spatial infinity.

B.2 Example solutions

a. No bubble. In the limit ν_B = ∞, i.e. a vanishing bubble, the system is radial symmetric about the centre of the colloid, and so we expect to recover the fundamental solution to Laplace's equation. The ODE for α, eqn (B20) reduces toand β = −α. We note that this limit cannot be recovered smoothly, since there is always a boundary layer in the forcing for s ≲ 1/(1 + ν_B) (corresponding to the contact region between bubble and colloid). This requires us to drop the boundary condition for α′(0) and just demand that the solution is bounded at s = 0 and ∞. In that case the ODE has the particular integralThe two complementary functions for α are e^s and e^sEi(−2s) which diverge at s = ∞ and s = 0 respectively and are therefore discarded. Henceas expected. From this we can conclude that the concentration scale on the colloid surface set by the flux condition is . If Henry's law requires the concentration on a hypothetical bubble surface to be larger than this value, then the bubble will tend to dissipate. Furthermore, the flux integrals (see eqn (B23)),inform us that there is no net flux into a surface of constant ν > 0 (as expected, since there is no bubble to absorb gas), and also that the flux out of the colloid in non-dimensional units is , or per unit area. This again confirms expectations.

b. A giant bubble. A giant bubble corresponds to the limit ν_B → 0. We note that without making any assumptions about ν_B, the boundary condition on the colloid, eqn (B21), implies that

β ≈ 2(ν_B + ξ)e^{−s(2ν_B+1)} − α for s ≫ (1 + ν_B)⁻¹, (B29)

⇒ A(s) ≈ (ν_B + ξ)e^−2sν_B for s ≫ (1 + ν_B)⁻¹. (B30)

Here, the range restriction on s expresses the validity of the approximation only on the ‘far’ side of the bubble. Due to our choice of length scaling, the range s ≲ (1 + ν_B)⁻¹ actually describes a boundary layer wherein the colloid touches the bubble and perturbs the flux of solute locally. It is therefore more instructive to reinsert dimensions and instead take the limit of a vanishing colloid, R_c → 0. It may then be seen that the boundary condition on the colloid, eqn (B3), reduces to the condition B(s) = 0, whereupon the asymptotic for A(s) becomes exact and extends to the whole range. Introducing the scaled variable w = (R_c/R_B)s we have (in dimensional units)where ξ* = ξ/R_c remains finite when the limit R_c → 0 is taken. In analogy to the case of no bubble, we again obtain the fundamental solution to Laplace's equation for the concentration field,We note that the net catalytic flux out of the bubble scales as 2γD(1 + ξ*R_B)/K_H. It therefore grows without bound with the bubble size unless the background fluid is completely saturated (ξ* = 0). Hence, unless this condition is met, there is an upper limit to bubble growth for any catalytic flux .

B.3 First order kinetics

a. Model and method of solution. In this section we examine first-order kinetics, which are a more realistic model of the chemical reaction on the colloid. Instead of assuming production of a constant flux of reaction product (oxygen), , we instead assume the presence of a second, ‘fuel’ field (hydrogen peroxide), c_f, that is converted on the colloid surface at a constant rate, k. This leads to two boundary conditions of the formIt is additionally necessary to solve diffusive dynamics for the fuel field, D_f∇²c_f = 0, and specify another boundary condition on the bubble. Since hydrogen peroxide is well soluble in water, we can assume that the bubble consists exclusively of reaction product, i.e. that there is no fuel flux in or out of the bubble. ThenFinally, we assume that far away the fuel concentration asymptotes to some bulk value c^∞_f that supplies the reaction.

Using the same scalings as before, the four boundary conditions on bubble and colloid then become

and feature two new non-dimensional parameters,The Damköhler number Da expresses the ratio of reaction speed to the diffusive recruitment of reaction fuel. A large Damköhler number indicates that fuel is depleted faster than it is replenished, which may be expected to inhibit bubble growth. The parameter Δ is just the ratio of diffusivities of product and fuel and not very important for the dynamics.

We note that if we define another parameter analogously to in zeroth order kinetics, with replaced by kc^∞_f, we can write the product boundary condition on the colloid as

which is more reminiscent of its form for zeroth order kinetics, eqn (B3). This limit then corresponds formally to Da → 0 with constant (and ), i.e. a slow reaction with a large fuel supply. In summary, first order kinetics are determined by four dimensionless parameters, , that generalise the two dimensionless parameters of zeroth order kinetics.

To solve the problem, we use the solution ansatz

which satisfies Laplace's equation and the boundary conditions at spatial infinity by construction. We note that here A_f and B_f have been defined with a prefactor that leads to more elegant dimensionless expressions. Substituting, rearranging and Hankel-transforming as before leads to four coupled second order differential equations,

β = (ν_B + ξ)e^−sν_B sech(s(1 + ν_B)) − tanh(s(1 + ν_B))α, (B44)

where α_f and β_f are defined as

α_f(s) = e^ν_BsA_f(s) − e^−ν_BsB_f(s) (B45)

β_f(s) = −e^−sA_f(s) + e^sB_f(s). (B46)

The first two of these equations uncouple from the rest and we can solve them numerically to find α_f and β_f. The boundary conditions are such that both α_f and β_f decay to zero at infinity (from the integration by parts), and that and (from requiring regularity of the solution at s = 0). Then with β_f in hand we can solve the second pair of equations to find α and β, with and α → 0 as s → ∞. Both of these steps are straightforward using a standard boundary value problem solver. This completes the solution. The case of zeroth-order kinetics is smoothly recovered as Da → 0, in which case α_f, β_f → 0, and the equations and boundary conditions for the product field reduce appropriately.

b. Interpretation. For typical values of the parameters (see Section II E), we can estimate the new parameters as

Da = 2.9R_c × 10⁻² μm⁻¹, Δ = 0.59, (B47)

and so the Damköhler number typically ranges from . We illustrate the effect of Da in Fig. 8. For small values, the concentration fields are essentially indistinguishable from zeroth-order kinetics, as one might expected. However, even for Da of order unity, the effect of including first-order kinetics is mostly an overall reduction of the flux out of the colloid, resulting in a shift of the line of stationary points in the phase diagram. In particular, there is no excessive perturbation to the c field in the gap between bubble and colloid, despite the depletion of fuel there. As a result, the time-dependent dynamics described in the main text remain unchanged, save for a quantitative shift in the colloid activity that is required to achieve a certain bubble size.

Fig. 8 Effect of the Damköhler number Da on the bubble growth dynamics. Left: Dashed and dash-dotted lines indicate theoretical estimates for R_max (stable) and R_min (unstable) respectively. The effect of a non-zero Damköhler number is approximately a shift of the nullclines towards larger values of flux . Right: Illustration of the fuel (top) and solute (bottom) concentration fields for small and large Da.

Appendix C: detailed solution of the quasi-steady fluid mechanical model

C.1 Setup

The classical work on Stokes flows in this geometry was carried out in ref. 46 for the motion of two rigid spheres in contact, and generalised in ref. 47 to fluid spheres. Here we broadly follow ref. 46 in developing the solution ansatz. Since the flow is axisymmetric, we employ a stream function to solve the problem. In cylindrical coordinates {ρ, θ, z}, such a function Ψ is defined by

u_ρ = ρ⁻¹∂_zΨ, u_z= −ρ⁻¹∂_ρΨ. (C1)

It may then be shown that this flow satisfiesin tangent sphere coordinates. Then, it can be shown that the Stokes equations reduce to

Λ⁴Ψ = 0, (C3)

where the Stokes operator Λ² is given by

Λ² = ρ∂_ρ(ρ⁻¹∂_ρ) + ∂_zz (C4)

A solution to the Stokes equations Ψ may be constructed from solutions to Λ²ϕ = 0 by adding linearly independent resonant terms. By considering the form of the Λ² operator in spherical and cylindrical coordinates respectively it may be shown that one way to achieve this is by writing

Ψ = (ρ² + z²)ϕ + zχ, (C6)

where Λ²ϕ = Λ²χ = 0. In this case we have irrotational flow, Λ²Ψ = 0, if and only if

ϕ + 2ρϕ_ρ + 2zϕ_z + χ_z = 0. (C7)

The reason we choose this representation among several equivalent ones (such as Ψ = ϕ + zχ) is due to it leading to a convenient form of the stream function later on.

In order to solve Λ²ϕ = 0 we write

It may then be shown⁴⁶ thatand with the ansatz F = exp(±sν)f(sμ) we see that f satisfies Bessel's equation of order one. The solution bounded for μ → 0 (i.e. no radial flow on the z-axis) may hence be written aswhere the functions f₊ and f₋ are determined by the boundary conditions. Using the definition of the coordinates μ and ν in eqn (C6) we then have a general solution of the formwhere the coefficients A to D are functions of s and to be determined from the boundary conditions. We stress that these are different from the coefficients of the diffusive problem, but in this section we drop the notation with a ∼ used in the main text to avoid clutter. We note however, that since J₁(sμ) → 0 as μ → 0, this ansatz is only capable of describing flows that have Ψ ≡ 0 on the axis of symmetry. This is the case for pure translation, but not for bubble growth since volume is not conserved. For this reason it is necessary to add a singular correction Ψ_s that satisfies Ψ_s = 2Ṙ_B/ν_B² for z ≥ 2R_B and Ψ_s = 0 for z ≤ 0 on the axis of symmetry. This is achieved by adding the stream function of a point source at z = R_B.⁴⁸ In terms of {μ, ν} coordinates this isThe complete stream function is henceFrom ref. 30 we have the formula 4–14.18 for the force on the body in the z direction,where n is the outward normal, η is the fluid viscosity and dι is the tangential length increment with a negative z-component (see Fig. 4–5.1 in ref. 30). Thus,both on the top and on the bottom sphere. We then find the forces on bubble, colloid, and the dimer to beThe final result can alternatively be obtained using the limit formula (4–14.19 in ref. 30),where .

C.2 Solution for no slip conditions on the bubble

Suppose the bubble-colloid system is moving with velocity U in the negative z-direction (i.e. with the colloid at the front), and the radius R_B of the bubble grows at a rate Ṙ_B. Then in the lab frame, the flow has velocity zero at infinity, −Uẑ on the colloid, and on the bubble it is Ṙ_B(ẑ − ) − Uẑ. Using eqn (A3), the velocity on the bubble ν = ν_B is found to bewhile on the colloid at ν = −1 it isBy comparing with the expression for the stream function in these coordinates, eqn (C2), we find the boundary conditionswhere a constant of integration has been chosen so that Ψ → 0 as μ → ∞ and μ → 0 with ν < 0, and Ψ → 2Ṙ_B/ν_B² as μ → 0 with 0 < ν < ν_B, thus matching the behaviour of Ψ_s. We solve for the motility and growth problems separately.

Writing , we note that the ν-derivative of the streamfunction is related to by

Hence,At the same time, it follows from a direct Hankel transform of Ψ thatEvaluated at ν = {ν_B, −1} these two relations give an algebraic system of four equations that determines the four coefficients in the solution. For the motility problem,

−(A − C)s⁻¹ sinh s + (B − D)s⁻¹ cosh s = 2U(s⁻² + s⁻¹)e^−s, (C29)

−[Cs⁻¹ + (B − D)]sinh s + [Ds⁻¹ + (A − C)]cosh s = 2Ue^−s, (C30)

(A + Cν_B)s⁻¹ sinh sν_B + (B + Dν_B)s⁻¹ cosh sν_B = 2U(s⁻² + ν_Bs⁻¹)e^−sν_B, (C31)

[Cs⁻¹ + (B + Dν_B)]sinh sν_B + [Ds⁻¹ + (A + Cν_B)]cosh sν_B = −2Uν_Be^−sν_B, (C32)

while for the growth problem the right hand side of these equations readsThe solution to this is quite messy, and since we only need the coefficient B for the force calculation, we will only report that one here. For the motility problem,while for the growth problem,where w = e^−2s. As s → 0, the coefficient A_gro ∼ −s⁻² → −∞ diverges, which according to eqn (C16) and (C17) leads to an infinite attractive force between the two spheres. The regularisation of this force is discussed in Appendix C.4. Unlike for the diffusive problem, there is no need to solve an ODE since Ψ and ∂Ψ/∂ν are specified on the same surfaces. This is no longer the case if we replace the boundary condition on the bubble with no shear stress.

C.3 Solution for no stress condition on the bubble

In this section we consider an alternative boundary condition on the bubble ν = ν_B. As in the rigid case, we have a no penetration condition for the normal -component of the velocity. In terms of the stream function, this reads (cf.eqn (C24)).The condition for the tangential [small mu, Greek, circumflex]-component is replaced by a no stress condition, . With some algebra it may be shown that this is equivalent toThe boundary conditions on the colloid remain unchanged. Using the no penetration condition eqn (C36) to simplify the right-hand side of eqn (C37), we arrive at

With the boundary conditions posed in this form, it turns out that there is an elegant path to a fully analytical solution. First, note that by performing two Hankel transforms of different orders we have the following identity,

Furthermore, from differentiating the general form of the streamfunction eqn (C13) we have thatwhere we define Greek coefficients as(α is defined only for completeness). After substituting the no stress condition eqn (C41) and applying the identity eqn (C42) this simplifies toIn analogy with the diffusive case (cf.eqn (B13)), we define the self-adjoint operatorwhich satisfiesIn particular, this allows us to writewhere the second equality follows from two integrations by parts if δ is bounded at 0 and ∞. This allows us to rewrite eqn (C45) asor, after applying another Hankel transform,The other three boundary conditions eqn (C38)–(C40) (no slip on the colloid and no penetration on the bubble) yield three linear equations which can be used to express α, β and δ in terms of γ. Then, we arrive at the modified Bessel equation,Considering the behaviour for small and large s, we find that the appropriate boundary conditions for a regular solution are γ(0) = − 2(U − Ṙ_B) and γ → 0 as s → ∞, which is consistent with the integration by parts performed earlier. Since γ is clearly linear in both U and Ṙ_B, a complete solution to the hydrodynamic problem for any bubble size ν_B may be found by solving eqn (C51). In fact, it is possible to do so by inspection. The result is

γ(s) = −2(U − Ṙ_B)e^−ν_Bs. (C52)

Bootstrapping our way up again to find the coefficient B we find the valueswhere w = e^−2s. For efficient numerical evaluation, we recommend that asymptotic expressions should be used in form of a Taylor expansion in s for s ≲ min(0.1, 0.1/ν_B) and the leading order in w, [scr O, script letter O](w^min(1,ν_B)), for s ≳ 10. As s → 0, the coefficient A_gro ∼ s⁻² → ∞ diverges, which according to eqn (C16) and (C17) leads to an infinite repulsive force between the two spheres. The origin and regularisation of this force is discussed in the next section.

C.4 Discussion of the pressure in the gap between bubble and colloid

In this paper we consider a setup in which the bubble and the colloid touch at a single point. As discussed in the main text, solving the hydrodynamic problem leads to infinite forces acting on the individual spheres, even though the net force on the system is finite. Since this is unrealistic, it is necessary to justify the physical validity of our results. In this section we investigate this divergence by solving for the flow near the contact point using a local asymptotic expansion of the Stokes equations. We demonstrate that the pressure diverges quadratically with a sign depending on the boundary condition, and hence that the divergence of the force due to having a singular contact point is logarithmic. As such, it is easily regularised even for small contact areas or separation between the spheres.

For clarity, we stick to dimensional units in this section. Working in cylindrical coordinates {ρ, z}, the surface of the bubble is given by z = h⁺ = R_B(1 − cos θ⁺) where tan θ⁺ = ρ/R_B. Expanding the cosine for small ρ we arrive at the leading order expression h⁺ = ρ²/2R_B. Similarly, we can expand the bubble surface as z = −h⁻ = −ρ²/2R_c. Since the correct asymptotic expansion of the Stokes equation depends on the shortest length scale, we define r₀ = min(R_b(t), R_c) and ε = ρ/r₀, which is our asymptotic expansion parameter. Geometrically we then have that at leading order

The length scales are given by ρ ∼ r₀ε and z ∼ r₀ε². Since ε ≪ 1 by assumption, we have ∂_z ≫ ∂_ρ as in lubrication theory. By inspection, the leading-order solution for velocity and pressure is given byIntroducing the scaled variablewe find that the leading order terms in the incompressible Stokes equations become

u₁ + ∂_z̄w₂ = 0. (C61)

Additionally we have a continuity condition of the deforming bubble interface due to growth,

For a rigid bubble, the boundary conditions to this system are u₁ = w₂ = 0 at z̄ = −r₀/R_c and u = −Ṙ_B(ν − z) at z = h⁺, where ν = cos θ⁺z − sin θ⁺ρ. At leading order this second condition becomes

It may be checked that this is compatible with the continuity condition, and we note that Dh/Dt is positive, so material points on the bubble move away from the colloid as the bubble expands. In particular, ∂_th and u_z|_z=h have opposite signs. Solving the Stokes equations with these boundary conditions leads to messy expressions for u₁(z̄) and w₂(z̄), as well as the dominant pressure term,Hence p ∼ ρ⁻² as ρ → 0, so the area integrated stress exerted by the spheres diverges logarithmically. Since the pressure is strictly negative, it corresponds to an attraction of the two spheres. Intuitively, this may be explained by visualising a stretching of the bubble surface as that of a balloon. While for fixed ρ the gap width h reduces as the bubble expands, individual points on the bubble surface move radially outwards at , and upwards at . Because of incompressibility fluid is sucked into the gap to compensate for the entrainment by the expanding bubble surface, creating an attraction.

This is not the case if the bubble surface exerts no shear stress. In this case the upper boundary conditions are replaced by ∂_z̄u₁ = 0 at z̄ = r₀/R_B and the continuity condition. This time the solution for the pressure is

which is strictly positive. Therefore, the bubble and colloid repel each other with this boundary condition, as one might expect intuitively. Differently from the rigid case, we haveso material points on the bubble move in to close the gap as the bubble expands. As before, the divergence of the pressure is quadratic, which means the repulsive force is easily regularised.

Appendix D: discussion of possible physics leading to premature bubble disappearance

As outlined in Section III C1, our model predicts the slow diffusive growth of the bubble to a final size that, in large regions of the parameter space, is significantly larger than the colloid. In experiments, these bubbles are generically not observed,^19,20 which led us to question possible mechanisms of premature bubble detachment. To this end it is necessary to revisit an assumption we have made in developing this model, namely that the bubble and colloid are both spherical and touch at an infinitesimally small point. If there is partial wetting, the bubble may form a spherical cap on the colloid with a small but non-zero contact radius that is determined by the surface properties of the colloid and the coefficients of surface tension between the gas, liquid, and solid materials that meet there. This is illustrated in Fig. 9. Here we introduce two angles, the contact angle θ, measured in the liquid phase, and the coverage angle ϕ which is related to the radius of the contact line r_c = sin ϕ in scaled units.

Fig. 9 Sketch of a more realistic geometry of bubble growth in which the bubble forms a spherical cap on the colloid. The contact angle θ and coverage angle ϕ are exaggerated for illustrative purposes. An internal pressure detachment force F_p and a contact line adhesion force F_γ compete in attaching the bubble to the colloid.

As a first check, we asked whether buoyancy would be sufficient to detach the bubble from the surface. For this, we compared the scaling of adhesive force along the contact line against buoyancy on the bubble volume as a function of the contact radius r_c. The point at which they balance is known as the Fritz radius r_f,²⁹ and the bubble detaches only if the true contact is smaller. From the scales of our problem we can estimate the Fritz radius as

where we assumed a length . Despite the L³-scaling rapidly increasing this radius for larger bubbles, this is so small that the bubble needs to cover only a vanishing fraction of the colloid surface in order to be tethered. In conclusion it seems unlikely that buoyancy is the cause of premature detachment.

A second check concerns the internal pressure in the bubble pushing against the colloid, F_p, compared to the contact line adhesive force from surface tension F_γ. Since the pressure difference between the bubble and the surrounding fluid is dominated by capillary pressure, we can approximate these as

where n is the normal to the colloid and t is tangential to the bubble. We ask whether the relative strength of these two forces depends on the angles involved, and whether this leads to a condition on the bubble radius R_B to which the angles θ and ϕ are geometrically related. It turns out however the reverse is true, and these forces always balance. Specifically, by considering the triangle formed by the sphere centres and the contact point it follows from the law of sines thatand the result F_p = F_γ follows from an easy substitution. While this may seem surprising, it just means that the contact angle θ is not controlled dynamically by the state of the system, but by the surface and material properties, which are unknown. They would be even more important when a moving contact line, and possibly pinning or surface roughness played a role. We also made one implicit assumption in this argument, which is that the fluid pressure is uniform. As shown in Appendix C.4, this is not the case in the gap, which could tip the balance in this argument. However, since the sign of this gap pressure is not dynamic and the bubble does not detach immediately, we suspect that the true boundary condition is approximately no slip, and the gap pressure indeed acts to further glue the bubble to the colloid.

In conclusion, we could not find an argument based on the intrinsic physics of the growth process that supports premature detachment. However, in separate experiments involving bubble formation in a thin film with immersed catalytic colloids, it was observed that coalescence with an external air–water interface can lead to the sudden disappearance of the bubble.⁴² While the depth of this film was not reported in previous work,¹⁹ we strongly suspect that the same external limitation of the system was causing the bubble disappearance there, rather than sudden inertial collapse that is hard to justify in an otherwise entirely non-inertial system. Arguing along similar lines, hydrodynamic interactions with the rigid interface on which the colloids are located may also hinder bubble growth. A detailed numerical investigation of this confinement however lies beyond the scope of this study.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We would like to thank Javier Rodriguez-Rodriguez for useful discussions. This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 714027 to S. M.). We thank Julien Husson for creating the graphical abstract.

References

1. W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St Angelo and Y. Cao, et al., Catalytic nanomotors: autonomous movement of striped nanorods, J. Am. Chem. Soc., 2004, 126(41), 13424–13431.
2. J. Elgeti, R. G. Winkler and G. Gompper, Physics of microswimmers – single particle motion and collective behavior: a review, Rep. Prog. Phys., 2015, 78(5), 056601.
3. P. Illien, R. Golestanian and A. Sen, ‘Fuelled’ motion: phoretic motility and collective behaviour of active colloids, Chem. Soc. Rev., 2017, 46(18), 5508–5518.
4. J. L. Moran and J. D. Posner, Phoretic self-propulsion, Annu. Rev. Fluid Mech., 2017, 49, 511–540.
5. B. J. Nelson, I. K. Kaliakatsos and J. J. Abbott, Microrobots for minimally invasive medicine, Annu. Rev. Biomed. Eng., 2010, 12, 55–85.
6. S. Jeon, S. Kim, S. Ha, S. Lee, E. Kim and S. Y. Kim, et al., Magnetically actuated microrobots as a platform for stem cell transplantation, Sci. Robot., 2019, 4(30), eaav4317.
7. R. Golestanian, T. Liverpool and A. Ajdari, Designing phoretic micro-and nano-swimmers, New J. Phys., 2007, 9(5), 126.
8. A. Chamolly and E. Lauga, Stochastic dynamics of dissolving active particles, Eur. Phys. J. E: Soft Matter Biol. Phys., 2019, 42, 1–15.
9. A. Chamolly, E. Lauga and S. Tottori, Irreversible hydrodynamic trapping by surface rollers, Soft Matter, 2020, 16(10), 2611–2620.
10. S. Tottori, L. Zhang, F. Qiu, K. K. Krawczyk, A. Franco-Obregón and B. J. Nelson, Magnetic helical micromachines: fabrication, controlled swimming, and cargo transport, Adv. Mater., 2012, 24(6), 811–816.
11. A. Walther and A. H. Müller, Soft, Nanoscale Janus Particles by Macromolecular Engineering and Molecular Self-assembly, Janus Particle Synthesis, Self-Assembly and Applications, Royal Society of Chemistry, 2012, pp. 1–28.
12. J. Li, X. Li, T. Luo, R. Wang, C. Liu and S. Chen, et al., Development of a magnetic microrobot for carrying and delivering targeted cells, Sci. Robot., 2018, 3(19), eaat8829.
13. F. Soto, E. Karshalev, F. Zhang, B. Esteban Fernandez de Avila, A. Nourhani and J. Wang, Smart materials for microrobots, Chem. Rev., 2021, 122(5), 5365–5403.
14. A. Ghosh and P. Fischer, Controlled propulsion of artificial magnetic nanostructured propellers, Nano Lett., 2009, 9(6), 2243–2245.
15. F. Mou, Y. Li, C. Chen, W. Li, Y. Yin and H. Ma, et al., Single-Component TiO2 Tubular Microengines with Motion Controlled by Light-Induced Bubbles, Small, 2015, 11(21), 2564–2570.
16. S. Michelin and E. Lauga, Phoretic self-propulsion at finite Péclet numbers, J. Fluid Mech., 2014, 747, 572–604.
17. J. L. Moran and J. D. Posner, Electrokinetic locomotion due to reaction-induced charge auto-electrophoresis, J. Fluid Mech., 2011, 680, 31–66.
18. J. G. Gibbs and Y. P. Zhao, Autonomously motile catalytic nanomotors by bubble propulsion, Appl. Phys. Lett., 2009, 94(16), 163104.
19. M. Manjare, B. Yang and Y. P. Zhao, Bubble driven quasioscillatory translational motion of catalytic micromotors, Phys. Rev. Lett., 2012, 109(12), 128305.
20. M. Manjare, B. Yang and Y. P. Zhao, Bubble-propelled microjets: Model and experiment. The, J. Phys. Chem. C, 2013, 117(9), 4657–4665.
21. G. Gallino, F. Gallaire, E. Lauga and S. Michelin, Physics of Bubble-Propelled Microrockets, Adv. Funct. Mater., 2018, 28(25), 1800686.
22. Z. Yang, T. N. Dinh, R. Nourgaliev and B. Sehgal, Numerical investigation of bubble growth and detachment by the lattice-Boltzmann method, Int. J. Heat Mass Transfer, 2001, 44(1), 195–206.
23. L. Zeng, J. Klausner and R. Mei, A unified model for the prediction of bubble detachment diameters in boiling systems-I. Pool boiling, Int. J. Heat Mass Transfer, 1993, 36(9), 2261–2270.
24. V. M. Fomin, M. Hippler, V. Magdanz, L. Soler, S. Sanchez and O. G. Schmidt, Propulsion mechanism of catalytic microjet engines, IEEE Trans. Robot., 2013, 30(1), 40–48.
25. L. Li, J. Wang, T. Li, W. Song and G. Zhang, Hydrodynamics and propulsion mechanism of self-propelled catalytic micromotors: Model and experiment, Soft Matter, 2014, 10(38), 7511–7518.
26. Z. Wang, Q. Chi, L. Liu, Q. Liu, T. Bai and Q. Wang, A viscosity-based model for bubble-propelled catalytic micromotors, Micromachines, 2017, 8(7), 198.
27. Z. Wang, Q. Chi, T. Bai, Q. Wang and L. Liu, A dynamic model of drag force for catalytic micromotors based on Navier-Stokes equations, Micromachines, 2018, 9(9), 459.
28. D. Lohse and X. Zhang, et al., Surface nanobubbles and nanodroplets, Rev. Mod. Phys., 2015, 87(3), 981.
29. J. Rodrguez-Rodrguez, A. Sevilla, C. Martnez-Bazán and J. M. Gordillo, Generation of microbubbles with applications to industry and medicine, Annu. Rev. Fluid Mech., 2015, 47, 405–429.
30. J. Happel and H. Brenner, Low Reynolds number hydrodynamics: with special applications to particulate media, Springer Science & Business Media, 2012, vol. 1.
31. E. Lauga and T. R. Powers, The hydrodynamics of swimming microorganisms, Rep. Prog. Phys., 2009, 72(9), 096601.
32. W. Zhai, Y. Song, Z. Gao, J. B. Fan and S. Wang, Precise synthesis of polymer particles spanning from anisotropic Janus particles to heterogeneous nanoporous particles, Macromolecules, 2019, 52(9), 3237–3243.
33. C. Kaewsaneha, P. Tangboriboonrat, D. Polpanich, M. Eissa and A. Elaissari, Janus colloidal particles: preparation, properties, and biomedical applications, ACS Appl. Mater. Interfaces, 2013, 5(6), 1857–1869.
34. B. Sabass and U. Seifert, Dynamics and efficiency of a self-propelled, diffusiophoretic swimmer, J. Chem. Phys., 2012, 136(6), 064508.
35. N. Sharifi-Mood, J. Koplik and C. Maldarelli, Diffusiophoretic self-propulsion of colloids driven by a surface reaction: the sub-micron particle regime for exponential and van der Waals interactions, Phys. Fluids, 2013, 25(1), 012001.
36. S. Inoue, Y. Kimura and Y. Uematsu, Ostwald ripening of aqueous microbubble solutions, J. Chem. Phys., 2022, 157(24), 244704.
37. J. Zhou, C. Wei and Y. Dong, How solid-liquid adsorption affects the liquid-vapor interface in pores, Vadose Zone J., 2022, 21(5), e20214.
38. M. S. Plesset and A. Prosperetti, Bubble dynamics and cavitation, Annu. Rev. Fluid Mech., 1977, 9(1), 145–185.
39. D. R. Lide, CRC handbook of chemistry and physics, CRC Press, 2004, vol. 85.
40. W. Yu, C. Batchelor-McAuley, X. Chang, N. P. Young and R. G. Compton, Porosity controls the catalytic activity of platinum nanoparticles, Phys. Chem. Chem. Phys., 2019, 21(36), 20415–20421.
41. X. Chang, C. Batchelor-McAuley and R. G. Compton, Hydrogen peroxide reduction on single platinum nanoparticles, Chem. Sci., 2020, 11(17), 4416–4421.
42. F. Yang, M. Manjare, Y. Zhao and R. Qiao, On the peculiar bubble formation, growth, and collapse behaviors in catalytic micro-motor systems, Microfluid. Nanofluid., 2017, 21(1), 6.
43. A. Chamolly and E. Lauga, Stokes flow due to point torques and sources in a spherical geometry, Phys. Rev. Fluids, 2020, 5(7), 074202.
44. S. H. Strogatz, Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering, CRC Press, 2018.
45. W. Gao, A. Uygun and J. Wang, Hydrogen-bubble-propelled zinc-based microrockets in strongly acidic media, J. Am. Chem. Soc., 2012, 134(2), 897–900.
46. M. Cooley and M. O'Neill, On the slow motion of two spheres in contact along their line of centres through a viscous fluid, in Mathematical Proceedings of the Cambridge Philosophical Society, Cambridge University Press, 1969, vol. 66, pp. 407–415.
47. L. Reed and F. Morrison Jr, The slow motion of two touching fluid spheres along their line of centers, Int. J. Multiphase Flow, 1974, 1(4), 573–584.
48. S. Michelin, E. Guérin and E. Lauga, Collective dissolution of microbubbles, Phys. Rev. Fluids, 2018, 3(4), 043601.

---