Surfactant-loaded capsules as Marangoni microswimmers at the air–water interface: Symmetry breaking and spontaneous propulsion by surfactant diffusion and advection

Abstract We present a realization of a fast interfacial Marangoni microswimmer by a half-spherical alginate capsule at the air–water interface, which diffusively releases water-soluble spreading molecules (weak surfactants such as polyethylene glycol (PEG)), which act as “fuel” by modulating the air–water interfacial tension. For a number of different fuels, we can observe symmetry breaking and spontaneous propulsion although the alginate particle and emission are isotropic. The propulsion mechanism is similar to soap or camphor boats, which are, however, typically asymmetric in shape or emission to select a swimming direction. We develop a theory of Marangoni boat propulsion starting from low Reynolds numbers by analyzing the coupled problems of surfactant diffusion and advection and fluid flow, which includes surfactant-induced fluid Marangoni flow, and surfactant adsorption at the air–water interface; we also include a possible evaporation of surfactant. The swimming velocity is determined by the balance of drag and Marangoni forces. We show that spontaneous symmetry breaking resulting in propulsion is possible above a critical dimensionless surfactant emission rate (Peclet number). We derive the relation between Peclet number and swimming speed and generalize to higher Reynolds numbers utilizing the concept of the Nusselt number. The theory explains the observed swimming speeds for PEG–alginate capsules, and we unravel the differences to other Marangoni boat systems based on camphor, which are mainly caused by surfactant evaporation from the liquid–air interface. The capsule Marangoni microswimmers also exhibit surfactant-mediated repulsive interactions with walls, which can be qualitatively explained by surfactant accumulation at the wall. Graphic Abstract


Introduction
Designing and understanding self-propelling biological or artificial microswimmers is the basis for the physics of active systems.Swimming on the microscale is governed by low Reynolds numbers and requires special propulsion mechanisms which are effective in the presence of dominating viscous forces.The first propulsion principle that comes to mind are shape-changing swimmers, which deform their body in a cyclic way in order to propel.At low Reynolds numbers, the cyclic deformation pattern of a swimmer must not be invariant under time-reversal according to the scallop theorem [1].In nature, many different examples of deformation swimmers can be found such as bacteria, algae and spermatozoa [2].Realizing this concept in synthetic microswimmers is often difficult as the scallop theorem requires control of at least two parameters.
Shape changing swimmers force the surrounding fluid via no-slip boundary conditions on the surface of their moving parts.Another successful class of synthetic mi-croswimmers are phoretic swimmers, which actively create slip velocities at their surface.Self-propelling phoretic swimmers autonomously create gradients in external fields such as concentration of a "fuel" or temperature, which in turn give rise to symmetry-breaking interfacial fluid flow in a thin interaction layer [3].This fluid flow constitutes an effective slip velocity leading to propulsion [4,5].Examples of such autophoretic swimmers are thermophoretic or diffusiophoretic swimmers, which generate gradients in temperature or concentration of interacting particles along their body.Self-diffusiophoretic swimmers generate a nonvanishing interfacial slip velocity on the particle surface via asymmetries in the solute concentration field and a short range interaction between solute and swimmer [3].Diffusiophoretic models typically neglect advection of the fuel concentration [6][7][8], but this has been included in Refs.[9,10].A lot of different aspects of swimmer behavior have been studied for self-diffusiophoretic swimmers such as efficiency [11], confinement effects [6,8] cargo trans-port [7,12], or the rich behavior during collisions with walls [13,14].
While diffusiophoresis creates concentration gradients within the liquid surrounding the swimmer, concentration gradients or surface active molecules (surfactants) within the interface of liquid swimmers can also generate symmetry-breaking interfacial forces based on the Marangoni effect [15].These propulsion mechanism based on the Marangoni effects are utilized in different liquid Marangoni swimmers, such as active liquid droplets or active emulsions [16].Examples are pure water in an oil-surfactant medium (squalane and monoolein) [17] or liquid crystal droplets in surfactant solutions [16] but many other systems can be generated making this a versatile route to microswimmer production.This type of Marangoni swimmer is a liquid drop fully immersed in a liquid carrying surfactant, and propulsion is generated by the Marangoni effect along the liquid-liquid interface between swimmer and surrounding liquid, where a surfactant concentration gradient is maintained.In Ref. [17], an auto-diffusiophoretic mechanism [9,18] has been proposed to maintain the surfactant concentration gradient.Another mechanism that has been proposed is increased adsorption of surfactant at the front (in swimming direction) of the swimmer, which depresses the interfacial tension in the front [16,19,20].This gives rise to a Marangoni stress towards the back (where the interfacial tension is higher).The Marangoni stress forces the surrounding fluid towards the back of the swimmer resulting in a swimmer motion towards the front of the swimmer.For all proposed mechanisms, the liquid swimmer autonomously maintains an increased surfactant concentration in the front of its interface with the surrounding liquid, and it propels in the direction of higher surfactant concentration at its own interface.
The self-phoretic and Marangoni swimming mechanisms discussed so far do not generate net forces on the swimmer but non-vanishing slip velocities on the particle surface via asymmetries in a temperature or solute concentration field.There is another class of self-propelling swimmers partly based on the Marangoni effect and with a long history [21], which are so-called soap or camphor boats (or surfers), which we call Marangoni boats in the following.Marangoni boats are moving at the liquid-air interface [22]; typically, they are solid swimmers and operate at the centimeter scale.They are often used as a popular demonstration experiment for the Marangoni effect [23].As "fuel" serve surface active molecules, which are deposited on the floating swimmer [23] or in which the swimmer is soaked [24][25][26][27][28][29], or the swimmer body itself is made from dissolving surfactant [30].There are many examples based on DMF (dimethylformamide) [31], alcohol [23,29], soap [29], camphor [24][25][26][27][28]32] or camphene [30] that have also been characterized quantitatively.The surfactant molecules are emitted or dissolved from the swimmer and a radial concentration gradient is established at the air-water interface by diffusion, eventually aided by evaporation for volatile surfactants.The radial concentration gradient creates (i) surface tension gradients and (ii) Marangoni stresses on the fluid.This leads, however, not necessarily to swimming as long as the surface tension is symmetric and uniform around the swimmer.The surface tension is pulling in normal direction on the closed air-water-swimmer contact line.A uniform surface tension cancels along any arbitrarily shaped closed three-phase contact line, but a gradient in surfactant concentration along the contact line can generate a net propulsion force.We call this net force generated by surface tension gradients direct Marangoni force in the following.Also symmetry-broken Marangoni flows created by the Marangoni effect can contribute to (or impede) the propulsion via hydrodynamic drag onto the swimmer surface.We denote the resulting forces that Marangoni flows exert by Marangoni flow forces in the following.The Marangoni boat mechanism is relying on both types of forces.If surfactant emission is anisotropic the boat is, in general, propelled into the direction of higher surface tension, i.e., lower surfactant concentration along the airwater-swimmer contact line.We note that this is opposite to the propulsion in the direction of higher surfactant concentration for the active liquid swimmers discussed before.The Marangoni boat mechanism is also employed by some insects (rove beetle and Velia) [33] to propel on the water surface.There are also recent experiments [34] and theoretical work [35] on a closely related system of thermally driven Marangoni boats or surfers.
A full quantitative theory of Marangoni boats including hydrodynamics, surfactant advection, direct Marangoni forces and Marangoni flows is still elusive despite previous progress [22,26,36,37].Some theoretical approaches ignore the advection [35,38,39], several ignore the hydrodynamic flow fields [24,25,32,[40][41][42] or approximate it by uniform flow [28], which clearly oversimplifies the description of surfactant transport.In particular, on the numerical side, a recent paper of Kang et el.[37] provides progress by including advection fully into the numerical solution for an anisotropic Marangoni boat.A theoretical description is complicated by the fact that most of the Marangoni boats operate at higher Reynolds numbers, and fluid flow generated during Marangoni propulsion is typically featuring vortices [29,37].Miniaturization to the microscale leads to low Reynolds numbers.Therefore, miniaturization is not only attractive for possible applications but also provides a starting point for the development of hydrodynamic theories, as the simpler linear Stokes equation holds for fluid flow at low Reynolds numbers.This has been initiated in Refs.[35,36,38,39].
Another question is regarding the role of intrinsic anisotropy, namely, whether isotropic swimmers with no intrinsically defined motion direction are also capable of a spontaneous motion which then spontaneously breaks the symmetry of the system.This question has been answered positively for autophoretic swimmers [9,18], where it has been shown that advection by the surrounding fluid can maintain the necessary gradients in fields and/or concentrations above a critical strength of the advection (characterized by a dimensionless Peclet number).Liquid Marangoni swimmers are always symmetric by construction and have to maintain an increased surfactant concentration in the front of their interface by adsorption of surfactant (or micelles) or by autophoretic effects [17,19,20].For the Marangoni boats the question regarding spontaneous symmetry breaking has been addressed experimentally in Ref. [28], where symmetric camphor disks have been shown to propel and swimming velocities have been shown to be largely independent of intrinsic swimmer anisotropy.So far, a theoretical answer is missing for Marangoni boats.
Here, we present a combined experimental and theoretical approach.We try to further approach the microscale by synthesizing alginate capsules as swimmer bodies, which provide a porous matrix that can accept surface active molecules.Several weakly surface active fuels are tested, among which polyethylene glycol (PEG) turns out to be the most effective.The PEG-alginate swimmers exhibit fast and prolonged propulsion.In general, we find prolonged propulsion only if spreading molecules are water-soluble as for PEG; then the air-water interface can regenerate by the fuel being dissolved in water.The PEG swimmers are approximately half-spherical, i.e., symmetric; therefore, we can address the question of spontaneous motion for a symmetric swimmer design.Moreover, a half-spherical shape turns out to be very convenient for theoretical modelling, and has also been employed in Ref. [37].For small Reynolds numbers, this geometry allows for a complete theoretical description of Marangoni boat propulsion by analyzing the coupled problems of surfactant diffusion and advection, fluid flow, which includes surfactant-induced fluid Marangoni flow, and surfactant adsorption at the air-water interface; we also include a possible evaporation of surfactant.The swimming speed is determined from the balance of Marangoni forces (both direct forces from surface tension gradients and from Marangoni flow forces) and drag forces.We can address the problem of spontaneous symmetry breaking and predict the swimmer's speed in a stationary state.This solution gives also hints how to generalize to higher Reynolds numbers using the concept of the Nusselt number, for which many results are known phenomenologically.
On the experimental side, we find further effects, such as the repulsive interaction of PEG-alginate swimmers with walls and the tendency to move in curved trajectories, which can be explained in the framework of the Marangoni boat mechanism.2 Alginate based capsule swimmers

Swimmer synthesis and characterization
The synthesized capsules show typical propelling mechanisms similar to phenomena observed for the insect class of Microvelia.Our artificial microswimmers consist of PEG droplets, which were surrounded by thin alginate shells.For the preparation of these particles, we first form an aqueous PEG-alginate composite solution (standard: w PEG300 = 0.5%, w alginate = 0.5%).A droplet of this mixture is then deposed on the surface of an aqueous CaCl 2 solution (standard: w CaCl2•2 H2O = 0.5%).The Ca 2+ ions serve as cross-linker and induce, within several microseconds, the gelation of the alginate membranes according to the box-egg model [43][44][45][46].Immediately after the formation of these particles, the capsules start to swim along the water surface.
Dripping microliter amounts of alginate into a crosslinker salt solution containing counterions starts an ionotropic gelation and produces approximately half-spherical alginate gel capsules of millimeter radius (see Fig. 2) [43].We report results for a ∼ 1500 µm; radii a ∼ 150 µm can be reached.For alginate gelation, different salt solutions can be used containing divalent cations such as CaCl 2 , CuCl 2 , or BaCl 2 solutions.
Adding surfactant to the alginate solution before dripping automatically loads the porous gel capsule with surfactant molecules.For suitable surfactants, capsules start to propel spontaneously on the air-water interface directly after dripping.The simple dripping technique allowed us to test many different "fuels": successfully propelling fuels are polyethylene glycols (PEGs) with molar weights 200−35000 g/mol, alcohols, acetic acid (stronger acids lead to protonization of alginate and subsequent coagulation), and organic solvents.A complete list of successfully tested fuels substances is given in  PEG (or polypropylene glycol (PPG)), in particular PEG 300, exhibit the best results regarding propulsion speed and propulsion duration; the reason is a suitable combination of diffusion constant, solubility, but also gelation properties of the alginate-PEG mixture.Corresponding monomers and dimers (ethylene glycol, propylene glycol, diethylene glycol) also exhibit good swimming properties but with lower speed and duration.It is particularly important for a prolonged propulsion that the fuel substance lowers the air-water surface tension but is also water-soluble such that it dissolves in the water reservoir after spreading in order to regenerate the air-water interface.Evaporation from the air-water interface is another mechanism to achieve such a regeneration, which is at work in camphor boats [26][27][28]32].Strong surfactants and detergents, such as sodium dodecyl sulfate, generate spreading pressures that can rupture the alginate capsule.Moreover, they quickly saturate the air-water interface such that concentration gradients and, thus, swimming can not be established.In the following, we report results for solutions of alginate and PEG 300 (w alginate = 0.5% and w PEG300 = 0.5%) dripped into a CaCl 2 crosslinker solution (w CaCl2•2 H2O = 0.5%).
Alginate gels have a porous structure [43][44][45].Scanning electron microscopy (SEM) of the alginate capsules reveals their porosity and also a certain roughness on the microscale with asperities on the capsule surface (see Fig. 3).The pores are essential for the slow diffusive emission of surfactant from the capsule [44,45].PEG diffusion through the porous alginate matrix is much slower than PEG diffusion in water; therefore, PEG should be released with a slowly varying controlled diffusive current that is limited by its slow diffusion in the alginate.The shape of the capsule and the spatial distribution of pores on the surface can break the overall spherical symmetry and give rise to small anisotropies in the emission, in principle.x-position in cm x-position in cm y-position in cm y-position in cm Fig. 5. Swimming trajectory of two PEG-alginate swimmers in a cylindrical container.Color-coded is the sign of the curvature, blue/red trajectories curve clockwise/counter-clockwise.Swimmers prepared according to same protocol exhibit individually different curving behavior (left: mostly counter-clockwise, right: mostly clockwise).Reflections at walls are of different duration.

Swimming motion
The alginate-PEG swimmers exhibit a fast and sustained motion.The swimming motion was observed in a cylindrical dish (diameter 24 cm) for up to 20 min.The swimmers exhibit typical speeds U swim ∼ 2 − 3 cm/s corresponding to 10 − 20 swimmer sizes per second (see Fig. 4); after 20 min, velocities U swim ∼ 1 cm/s can still be measured.This swimming performance is comparable to camphor boats [27,28] and active liquid droplets [16,17].
A typical swimming trajectory (lasting 84 s) far from a wall is shown in Fig. 4. We obtained this trajectory from a single particle tracking analysis (using ImageJ); typical swimming velocities are U swim ∼ 2 − 3 cm/s corresponding to 10 − 20 swimmer sizes per second.This corresponds to moderate Reynolds numbers Re = ρU swim 2a/µ ∼ 60 (with the swimmer diameter 2a 3000 µm as length scale and the viscosity and density of water, µ 10 −3 Pas and ρ = 10 3 kg/m 3 ).
Swimming trajectories such as in Fig. 4 and in confinement in Fig. 5 exhibit phases with a characteristic curvature.Marking the swimmer with elongated plastic fragments shows that the elongated fragment is not turning with respect to the direction of motion, i.e., the curving of the trajectory is correlated with a reorientation of the swimmer.This is a hint that during curved swimming the swimming direction is linked to the orientation of the particle and, therefore, that the spherical symmetry is slightly broken by irregularities in the porous structure of the alginate particle (see SEM pictures in Fig. 3).The swimming direction is selected by dominating pores which determine a preferred direction of emission and, thus, propulsion by the resulting surfactant gradients.Curving itself can be caused by additional torques from asperities of the alginate capsule where surfactant is emitted preferentially in the tangential direction.A similar mechanism is at work at camphor-driven rotors [42,47].This is supported by the finding that the curving behavior of swimmers prepared by the same protocol (such as the swimmers in Fig. 5) is individually different and seems to depend on small differences between irregularities acquired in the preparation process.Recently, also vortex shedding at Reynolds numbers Re ∼ 100 − 200 have been proposed to cause curving of trajectories [29].Swimmers are also repelled by walls and reverse their direction of motion normal to the wall.In a course of a collision in normal direction, the swimmer keeps, however, its orientation while the direction of motion is reversed, i.e., during normal wall collisions the swimming direction also reverses with respect to the particle orientation.Swimming direction reversal has also been observed for camphor boats [22,24,25].Swimming trajectories in Fig. 5 also show collisions with walls that last longer; these collisions can also feature a reorientation of swimmer, similar to what has been predicted for self-diffusiophoretic swimmers [13].

Swimming mechanism
The order of magnitude of swimming speeds can only be explained as a result of a modulation of the large liquid-air surface tension.Marangoni mechanisms based on surface tension variations within the gel-liquid interface between alginate capsule and surrounding water are unlikely because the interfacial tensions and, thus, also Marangoni stresses, are too small for solid-liquid or gel-liquid interfaces.This hints at a Marangoni boat propulsion mechanism for the alginate-PEG swimmers.
There is further experimental evidence supporting the Marangoni boat mechanism: (a) Sinking capsules stop swimming which excludes a phoretic or Marangoni mechanism based only on the swimmer-liquid interface such as the active liquid droplet mechanisms [16,17].(b) Only watersoluble spreading molecules lead to prolonged propulsion because they allow regeneration of the air-water interface by re-dissolving after spreading, which is crucial to establish concentration gradients at the air-water interface.(c) Local Wilhelmy plate surface tension measurements demonstrate surface tension modulations depending on the distance to a swimmer; this demonstrates that additional surfactant is emitted close to the swimmer.(d) Particle image velocimetry (PIV) measurements and selective staining with anilin show fluid motion consistent with surfactant spreading by surface tension reduction.(e) Swimming speed depends on the diffusive mass outflux.We will develop a theory for this dependence in the theoretical part of the paper, which describes the data (without free fitting parameters).(f) Repulsion and direction reversal without reorientation of the swimmer can be explained by an accumulation of surfactant emitted by the swimmer in front of the wall.This points to a motion opposite to the surfactant concentration gradient, while the direction of motion is not completely fixed relative to swimmer orientation.
More details regarding points (d) and (e) are given below.All these results suggest that the swimmer diffusively emits surfactant which reduces the surface tension.The swimmers are spherically symmetric to a good approximation and this symmetry is strongly broken by the concentration profile in the fast moving state.The only available mechanism for symmetry breaking in the moving state is by advection to the surrounding moving liquid, which selects a swimming direction spontaneously.
The results regarding curved trajectories and wall collision suggest that spherical symmetry is not perfect and large pores in the capsule shell can select a weakly preferred propulsion direction and link capsule orientation to swimming direction.This weak link can be deleted during a normal collision with a wall, when the swimmer reverses direction without changing orientation.

PIV measurements
PIV measurements were performed with Polymethyl methacrylate (PMMA) tracer particles with sizes between 30 and 50 nm and visualize the fluid flow close to the air-water interface.Figure 6 shows the results directly after swimmer synthesis by dripping (A), i.e., in the initial starting phase of the swimmer and (B) shortly after the swimming started.
In Fig. 6(A) we observe strong radial spreading of surfactant by initial Marangoni flows.Then, the symmetry is spontaneously broken when swimming is initiated and Fig. 6(B) shows the fluid surface flow in the initial swimming stage.During swimming we still observe radial Marangoni flows (blue) but the fluid flow around the swimming object creates a tangential backwards component (green) because two counter-rotating vortices form; moreover, Karman-like vortices appear on the rear side (red).Vortex formation demonstrates that the fluid motion happens at moderate Reynolds numbers Re ∼ 60.Similar vortex structures have also been observed in Ref. [29] for disks propelled by alcohol.Nevertheless, Reynolds number are moderate (Re 200) such that we can expect a steady fluid flow (eventually with boundary layer separation from the sphere and stationary Marangoni vortices).Only at higher Reynolds numbers Re > 200, we expect unsteady or even turbulent flow around a sphere [48].

Mass outflux and velocity measurements
We propose that the swimming motion is caused by surfactant that is diffusively emitted by the PEG-alginate capsule.Therefore, a relative slow reduction of the total mass of the PEG-alginate capsules should be measured, which also correlates with the swimming speed.Overall spherical symmetry of the capsule implies that the emission current density α is uniform on the capsule surface.
Quantitative measurements of the mass outflux are difficult.In Ref. [28] this has been achieved only indirectly by measuring the increase in surfactant in the surrounding solution.Here, we measure the mass outflux directly by removing swimmers (prepared according to the same protocol) after times t = 1, 2, 3, ...min from the swimming solution, dry freezing the swimmers to completely remove water from the alginate hydrogel, and determine their weight, which gives the mass m(t) of the swimmer at times t = 1, 2, 3, ...min.Figure 7(left) shows the results for the mass averaged over 10 swimmers, error bars are the standard deviation.
Diffusive outflux through a porous shell of thickness h (h < a) approximately follows an exponential decay.The emission current density is α ≈ −D s (c i − c 0 )/h, where D s is the diffusion constant in the gel, c i the interior PEG concentration and c 0 the exterior PEG concentration in solution.We assume that PEG has to diffuse over a fixed distance h for release; more refined release models use a time-dependent diffusion distance [28,49,50].For a halfsphere, this results in ċi = 3α/a = −3D s (c i −c 0 )/ha and an exponential decay of c i (t) and, thus, m(t).This motivates an error-weighted least-square fit with an exponential which describes the data well with an "empty" mass m ∞ ≈ 180 µg, a mass loss m 0 ≈ 182 mug, and a time constant τ m ≈ 9.95 min (see Fig. 7(left)).This confirms a slow diffusional surfactant release, i.e., α is changing on a large time scale τ m ; this time scale is large compared to any microscopic time scale of the fluid flow and the surfactant diffusion.Therefore, we expect that fluid flow and surfactant concentration are always in a quasi-stationary state during swimming, i.e., adiabatically following a slowly changing α.Differentiating with respect to time gives the mass outflux ṁ as a function of time (see Fig. 7(middle)).
The corresponding swimming velocity U swim of the PEG-alginate swimmers is measured via the single particle tracking analysis.Figure 7(right) shows the results for the velocity averaged over 10 swimmers prepared according to the same protocol as for the mass measurements, error bars are the standard deviation.The data clearly shows a fast initial drop of the velocity in a first phase, followed by a slower decay in a second phase.During the first phase, slow diffusion of surfactant through the porous alginate matrix might not be necessary yet and gelation of the capsule alginate shell might also be incomplete.The existence of several swimming phases has also been observed for camphor disks in Ref. [27].The second phase should be characteristic for propulsion triggered by slow diffusional release as described above.This motivates a fit with two exponentials.The resulting error-weighted leastsquare fit describes the data well as shown in Fig. 7(right).The first phase has a time constant τ u,0 0.37 min (and u 0 0.74 cm/s), whereas the second phase has τ u,1 15.89 min (and u 1 2.27 cm/s), which is comparable with τ m .This further supports that diffusive release of surfactant causes the propulsion during the second phase.
3 Theoretical results

Model
In the theoretical part of the paper we focus on the dependence of swimming speed on diffusive surfactant release.So far, this important question has not received attention in the literature.The strategy to calculate the swimming speed is as follows.We first prescribe a stationary velocity U = U e z of the swimmer and analyze the following three coupled problems for their stationary state: (i) Surface tension reduction by surfactant adsorption at the air-water interface; depending on the volatility of the surfactant we also need to include a possible evaporation of surfactant.PEG is not volatile but the theory will apply to the physics of the Marangoni boat mechanism in general and should also explain quantitative  results on camphor boats from Ref. [28].As opposed to PEG, camphor is a volatile surfactant which quickly evaporates from the air-water interface.(ii) Fluid flow, which includes both the fluid flow induced by motion of the half-spherical capsule and additional surfactant-induced Marangoni flow inside the fluid.(iii) Diffusive surfactant release from the swimmer and subsequent diffusion and advection.
Solving these three coupled problems we can obtain the Marangoni forces as a function of the prescribed velocity U from the surfactant concentration profile.Finally, the actual swimming velocity U = U swim is determined by the force equilibrium between drag force, direct propelling Marangoni forces from the surface tension gradient along the air-water-swimmer contact line, and Marangoni flow forces, which can increase either the drag or the direct Marangoni propulsion force.The fluid flow part (ii), the drag force, and also the Marangoni forces in the force balance strongly depend on the Reynolds number.Although the Reynolds number for the PEG-alginate swimmers is moderate (Re ∼ 60), we will first develop a low Reynolds number theory, and try to address higher Reynolds numbers afterwards using phenomenological results for the Nusselt number.
We introduce coordinates such that the origin r = 0 is at the center of the circular planar solid surface of the half-sphere, and the liquid-air interface is at y = 0 (with y < 0 being the liquid phase) and e z will coincide with the spontaneously selected swimming direction.We also use spherical coordinates such that θ = 0 is the swimming direction and the interfacial plane is located at φ = 0, π (y = 0).The half-sphere has radius a such that the contact line is at r = a and φ = 0, π (and parametrized by θ).We denote the half-spherical surface of the swimmer by S, the circular air-water-swimmer contact line by L, and the liquid-air interface outside the swimmer as S Int , see Fig. 8.

Surface tension reduction by surfactant adsorption and evaporation
We start with problem (i), which is independent of the Reynolds number.In equilibrium, the surfactant concentration Γ (r) at the interface y = 0 for a given bulk subsurface concentration c(r) is given by Langmuir adsorption (with the adsorption equilibrium constant K L and the maximal surfactant surface concentration Γ max ).In Lang- muir adsorption, we assume ideal behaviour of the surfactant molecules.According to the Gibbs adsorption isotherm, the interfacial surface tension γ is related by dγ = −RT Γ d ln c to surface concentration and bulk concentration [51].Together with the Langmuir equation, this leads to the Szyszkowski equation (with the gas constant R = N A k B and Γ max in mol per area), formulated for small local variations around a background, c(r) = c 0 + ∆c(r).Small surfactant concentration variations thus linearly reduce the surface tension, where r is an interfacial vector with y = 0. We choose the background c 0 as the bulk value c 0 = c(∞) for |r| → ∞.
In formulating eq. ( 5) locally, we already assumed that the on and off kinetics of surfactant to the interface is fast such that equilibrium can be assumed to be established instantaneously at every point r on the interface.Then also the surface concentration Γ (r) is slaved to the bulk and only a passive "reporter" of the bulk subsurface concentration c(r)| y=0 , and we do not have to solve a separate dynamics for Γ (r) in the interface.This assumption is typically valid for small surfactant molecules [52], in particular for water-soluble spreading molecules such as PEG.The assumption also implies that there is no flux imbalance within the interface, and also the bulk diffusive flux to the interface S Int has to vanish, which provides the corresponding boundary condition to the diffusion-advection problem (iv) in the bulk.Here, D is the surfactant diffusion constant in the bulk liquid.The surface concentration Γ (r) should also be small enough to avoid saturation of the air-water interface, which also requires water-soluble spreading molecules such as PEG.So far, we did not consider the possibility of surfactant evaporation from the interface.This enters the balance of fluxes to and from the interface (see Fig. 8), and we have to replace eq. ( 6) by where k is the rate constant for evaporation.

Fluid flow at low Reynolds numbers
We consider the rest frame of the swimmer and linearly decompose the total fluid flow field into a field v(r), which is the flow field of a half-sphere pulled with velocity U e z through the liquid and a correction v M (r) from Marangoni flows, v tot (r) = v(r) + v M (r).For low Reynolds numbers, both v(r) and v M (r) (and the associated pressure fields) fulfill the incompressibility condition ∇ • v = 0 and the linear Stokes equation µ∇ 2 v = ∇p, where µ is the fluid viscosity.The Stokes equations for v(r) and v M (r) are decoupled because of linearity; this will be different at high Reynolds numbers.The flow field v(r) of an externally pulled half-sphere has no-slip boundary conditions on the surface of the sphere, stress-free boundary conditions at the liquid-air interface, and v(∞) = −U e z at infinity.The total flow field v tot (r) also has no-slip boundary conditions on the surface of the sphere, assumes v tot (∞) = −U e z at infinity, but is subject to Marangoni stresses at the liquid-air interface.Consequently, the difference v M (r) = v tot (r) − v(r) from Marangoni flows has no-slip boundary conditions on the surface of the sphere, has vanishing velocity v M (∞) = 0 at infinity, and is subject to Marangoni stresses at the liquid-air interface.Moreover, for all three flow fields, there is no normal flow across the liquid-air interface v tot,y (r)| y=0 = v y (r)| y=0 = v M,y (r)| y=0 = 0. We will assume that the liquid-air interface remains flat, even if the sphere moves.This requires that typical viscous forces remain small compared to interfacial stress, µU γ, which is fulfilled with µU ∼ 10 −5 N/m for U ∼ 1 cm/s and γ ∼ 0.07 N/m for the air-water interface.We also neglect a possible curvature of the interface from wetting effects.
The Marangoni flow is caused by tangential Marangoni stresses at the liquid-air interface y = 0, (8) which act both on v M and v tot .
At low Reynolds numbers, the flow field v(r) is given by "half" (y < 0) the Stokes flow field around a sphere, which automatically fulfills the boundary condition v y (r)| y=0 = 0 for symmetry reasons.In spherical coordinates, the axisymmetric Stokes flow field is

Surfactant diffusion and advection
Surfactant molecules are emitted from the half-spherical surface S and diffuse into the liquid phase.At the same time, they are advected by the total fluid flow.In the stationary limit, the bulk concentration field is governed by the diffusion-advection equation Because of the slow diffusional surfactant release the appropriate boundary condition on S is a constant flux boundary condition, together with c(∞) = 0 and the no-flux boundary condition ( 6) at the interface S Int .The flux α is only slowly changing (on the time scale τ m ) and approximated as a constant for the calculation of quasi-stationary fluid flow and concentration fields.

Drag and Marangoni forces at low Reynolds numbers
The half-spherical swimmer moving at velocity U is subject to three forces.First, there is the drag force, which is, at low Reynolds numbers, given by the Stokes drag for a half-sphere, Second, there is the direct Marangoni propulsion force F M = F M e z from integrating the surface stress ∆γ(r) = −κc(r) along the air-water-swimmer contact line, Third, there is the Marangoni flow force F M,fl = F M,fl e z , which is by definition the force transmitted by fluid stresses of the Marangoni flow onto the sphere, For low Reynolds numbers, we can apply the reciprocal theorem to the flow fields v and v M and their associated stress tensors to calculate the Marangoni flow force without explicitly calculating the Marangoni flow v M , as has been shown in detail in Ref. [53].This gives the identity 0 = S+SInt da i v j σ M,ij , which leads to a Marangoni flow force The total Marangoni force is obtained by using ( 13) and the Gauss theorem, Because ρ −1 − ρ −3 > 0 for ρ > 1, the total Marangoni force is always positive for concentration profiles, which are increasing towards the rear side.This implies that the total Marangoni force is always propulsive.i.e., points in the same direction as the imposed velocity U regardless of its absolute value.This is a necessary condition for selfpropulsion.When the particle is pulled by an external force, this also implies that the total Marangoni force will always support the pulling force instead of increasing the drag.The Marangoni flow contribution F M,fl , however, can have both signs.For F M,fl > 0, the flow force increases the direct Marangoni force resulting in F M,tot > F M ; for F M,fl < 0, the flow force is directed backwards and increases the drag force resulting in F M,tot < F M .As opposed to Ref. [38], we will find that both cases are possible.
In the stationary swimming state, drag and total Marangoni force have to be balanced, such that the swimmer is force-free.This is the swimming condition that finally determines the actual swimmer velocity U = U swim .

Non-dimensionalization
To proceed, we make the system of coupled equations governing our sub-problems (i)-(iii) and the Marangoni forces dimensionless by measuring lengths in units of a, velocities in units of D/a, concentrations in units of αa/D, and forces in units of Dµ.We introduce The prescribed dimensionless velocity Ū of the swimmer is the first control parameter of the problem, 1 which is related to the Reynolds number, Re = 2 Ū /Sc, via the Schmidt number Sc ≡ µ/ρD.

Our dimensionless set of equations for problems (i)-(iii) becomes (i)
− ∇c(ρ) • n out ȳ=0 ≈ 0 without evaporation, (20a) with the dimensionless Peclet-number where ṁ = 2πa 2 α is the mass loss per time of the swimmer (see Fig. 7).The Peclet number is a dimensionless measure of propulsion strength.Typical values for the PEG-alginate swimmer are very high, Pe ∼ 10 7 (see Tab. 2).We also introduced the dimensionless Biot number governing possible evaporation, which is practically absent for PEG.From eq. (20d), we see that the Peclet number Pe determines the velocity scale of the Marangoni flow field.Therefore, we can also assign a Reynolds number Re M = 2Pe/Sc = RePe/ Ū to the Marangoni flow.Typical values for the PEG-alginate swimmer are Re M ∼ 10 4 (see Tab. 2), which agrees with the experimentally observed turbulent features of Marangoni flows (see Fig. 6).Via the advection with v(ρ) + vM (ρ), the concentration field c(ρ) depends both on the dimensionless velocity scale Ū of the Stokes field and the dimensionless velocity scale Pe of the Marangoni flow field, in general.All dimensionless parameters are summarized in Tab. 2, along with typical values for the PEG-alginate swimmers and in comparison to camphor boats according to Ref. [28].
In the following, we will solve the problems (i)-(iii), in order to obtain the dimensionless direct and total Marangoni forces (see eqs. ( 13) and ( 17)) from the concentration profiles by for a prescribed swimmer velocity Ū .Finally, force balance gives the dimensionless version of the swimming condition (18), − FD = 3π Ūswim = FM (Pe, Ūswim ) + FM,fl (Pe, Ūswim ), (25) which then selects the actual swimmer velocity Ū = Ūswim as a function of the remaining control parameters Pe ("fuel" emission) and eventually k (evaporation).
The non-dimensionalization reveals that the coupled problems (i)-(iii) and the Marangoni forces depend on three dimensionless control parameters (see also table 2): First, the prescribed dimensionless velocity of the swimmer Ū ; second, the Peclet number Pe characterizing the strength α of the surfactant emission, and third, the Biot-number k characterizing the evaporation.We also see that the Peclet number both controls the strength of the Marangoni flow via eq.( 20d) and the strength of all Marangoni forces.We note, however, that FM /Pe and FM,tot /Pe still depend on Ū and Pe via the dependence of c(ρ) on these parameters.
Another important finding from non-dimensionalization is that the diffusion-advection problem decouples from the Marangoni flow problem for Pe Ū , where we can neglect v M in the advection term.Then, the concentration profile is only determined by Stokes flow, becomes axisymmetric, and only depends on Ū .In this limit, the Marangoni flow field need not to be calculated in order to calculate the total Marangoni force via eq.( 24).

Numerical methods
Numerically, we only address the low Reynolds number regime.In general, we consider the problems (i)-(iii), i.e., solve the coupled diffusion-advection problem and the Marangoni flow problem for a prescribed swimmer velocity Ū .From the solution for the concentration field, we then calculate the Marangoni forces as a function of Ū and Pe in order to finally solve the force balance swimming condition.We use an iterative finite element scheme to solve the full coupled problem; this approach is explained in detail in Appendix A.

Low Reynolds numbers Re
1 are realized for Ū Sc/2, which can still be much larger than unity as typical Schmidt numbers for surfactants in aqueous solutions are of the order of 1000 (see Tab. 2).Therefore, we have to discuss both the diffusive limit Ū 1 and the advective limit Ū 1. 3.4.1 Decoupled limit Pe Ū First, we will consider the limit Pe Ū , where the diffusion-advection problem for a half-sphere with prescribed velocity U decouples from the Marangoni flow problem.We also focus on the case in the absence of evaporation first, as it is appropriate for the PEG-alginate swimmer.
Diffusive release from a moving emitter or from a resting emitter in a fluid flow can be characterized by the average Nusselt number (or Sherwood number Sh), which is the dimensionless ratio of the total emitted flux and the typical diffusive flux [57].A Nusselt number of one is realized for purely diffusive transport, Nusselt numbers much larger than one indicate strong advective transport.
In the decoupled limit Pe Ū , we find for the Nusselt number where c 0 (ρ) ≡ 1 2 π 0 dθ sin θc(ρ, θ) is the zeroth Legendre coefficient.There are two regimes, a diffusive regime for Ū 1 characterized by a Nusselt number close to one and an advective regime for Ū 1 where the Nusselt number becomes large, which is also clearly supported by the numerical results in Fig. 9.The result ( 27) can be derived analytically [58], apart from the value of the prefactor 0.65, which we determined numerically from the data in Fig. 9.A short derivation based on scaling arguments for the advective regime is presented below.The numerical results in Fig. 9 show perfect agreement and clearly confirm the existence of just two regimes.
The Nusselt number has been originally defined for constant concentration boundary conditions c(1, θ) = 1, For constant concentration boundary conditions, the result is well-known [54][55][56][57] and very similar (see Fig. 9), with a prefactor that can be calculated analytically.This indicates that literature results for the Nusselt number for constant concentration boundary conditions also apply to our situation of constant flux boundary conditions, which is an important insight that we will assume to also hold for higher Reynolds numbers below.
In the decoupled limit Pe Ū , we also find find a diffusive and an advective regime for the Marangoni forces FM,tot where numerical constants d M and d M,tot are obtained from the numerical results, see Fig. 10, and R is the radial system size.Again, the numerical results (Fig. 10) show perfect agreement and clearly confirm the existence of just two regimes, a diffusive and an advective regime.Direct and Fig. 10.Marangoni forces FM/πPe and FM,tot/πPe as a function of imposed velocity Ū together with corresponding concentration profiles (in the z x-plane) in the decoupled limit Pe Ū .All results are from numerical FEM solutions of the axisymmetric diffusion-advection equation in two-dimensional angular representation with ρ < R = 30.In the advective regime Ū 1 a concentration boundary layer develops (see eq. ( 31)).
total Marangoni force reach maximal values FM , FM,tot ∼ 0.15πPe in the crossover region Ū ∼ 1 between diffusive and advective transport.
Also the results ( 29) and ( 30) can be derived analytically from a calculation of the concentration field [58], apart from the value of the numerical constants.Here we present a short derivation based on scaling arguments.In the diffusive limit Ū 1, there is a linear response of the concentration field, which leads to a linear response of the Nusselt number and Marangoni forces.The coefficients can be calculated by perturbation theory about the concentration field c (0) (ρ) = 1/ρ at Ū = 0 in powers of Ū .A remarkable result of this calculation is that the linear term for the total Marangoni force diverges logarithmically with the system size R, see eq. ( 30), while the linear term for the direct Marangoni force stays finite.This means that the Marangoni flow forces strongly increase the direct force for Ū 1; such a behaviour could not be found in Ref. [38].For very large system sizes R 1/ Ū , the large scale cutoff R in (30) will be replaced by 1/ Ū because the convection term can no longer be treated perturbatively in the region ρ 1/ Ū , regardless how small Ū is [55].In the limit of strong advection Ū 1, a concentration boundary layer develops around the half-sphere, as can be clearly seen in the concentration profiles in Fig. 10).Its width ∆r is determined by the distance that a surfactant molecule can diffuse during the time ∆t ∼ a/v(∆r/a) (see eq. (9b)) it takes to be transported along the sphere by advection: ∆r 2 ∼ D∆t.Because v(∆r/a) ∼ U ∆r/a (see eq. ( 9b)), we find This is a classic result for the diffusion-advection problem for constant concentration boundary conditions [54,57], but also holds for constant flux boundary conditions.Because the concentration will drop in radial direction from its value at the surface S of the half-sphere to zero within the concentration boundary layer of width ∆ρ, we also have 1 = −∂ ρ c(ρ = 1, θ) ∼ c(ρ = 1)/∆ρ, which leads to a scaling of the symmetry-breaking concentration level at the surface S of the half-sphere for constant flux boundary conditions.For strong advection, the Marangoni forces decrease as a function of Ū because the concentration boundary layer width ∆ρ, in which forces are generated, and the concentration level around the sphere, to which forces are proportional, both decay with velocity as Ū −1/3 .The scaling property (32) for c(ρ = 1, θ) directly explains the results (27), Nu ∼ 1/c(ρ = 1, θ) ∼ Ū 1/3 , for the Nusselt number and (29), FM /Pe ∼ c(ρ = 1, θ) ∼ Ū −1/3 , for the direct Marangoni force in the advective limit Ū 1.They are clearly confirmed by all numerical results in Figs. 9 and 10.The result for the total Marangoni force (30) seems to deviate from this advective scaling Here, the expected scaling from the radial concentration boundary layer of width ∆ρ is FM,tot /Pe ∼ ∆ρ 2 c(ρ = 1, θ) ∼ Ū −1 (see eq. ( 24)), which is clearly not found in the numerics (yellow line in Fig. 10).The reason is that this contribution is actually only sub-dominant.The leading contribution comes from the advective tail of angular width ∆θ ∼ Ū −1/3 ; the width of the tail follows from the scaling of the stream function ψ ∝ r 2 ∆θ 2 in the tail and ψ ∝ 3∆r 2 sin 2 θ/2 in the boundary layer and the fact that a fluid particle should follow a Stokes flow stream line ψ = const in the advective limit.Therefore, the dominant contributions in eq. ( 24) are FM,tot ∼ Pe∆θc(ρ = 1, θ) ∼ Ū −2/3 in agreement with the numerical results in Fig. 10.
Comparing the curves for direct and total Marangoni force in Fig. 10, we observe a crossing such that FM < FM,tot in the diffusive regime Ū 1, while FM > FM,tot in the advective regime Ū 1.This means that the Marangoni flow force FM,fl = FM,tot − FM increases the propulsion force in the diffusive regime Ū 1 but decreases the propulsion force (or increases the drag) for Ū 1.This subtle result is related to the structure of the Marangoni flows, which are generated by the surfactant concentration gradients ∇S c(ρ) ȳ=0 within the liquid-air interface (see eq. ( 20d)) and can be qualitatively rationalized with the help of eq. ( 15) for the Marangoni flow force.We can decompose the surfactant concentration gradient into tangential and radial components Because of advection the tangential component points from the front to the rear corresponding to an increasing surfactant concentration towards the rear (∂ θ c > 0).It gives rise to a forward Marangoni flow and F M,fl > 0 in eq. ( 15) because −e z • ρ −1 ∂ θ ce θ ∝ sin θ > 0. This is the dominating effect in the diffusive regime Ū 1, where the perturbation theory gives to leading linear order in Ū a contribution of the form Ū c(1) ∝ − Ū c1 (ρ) cos θ resulting in ρ −1 ∂ θ c = Ū ρ −1 c1 (ρ) sin θ > 0. The front-directed tangential Marangoni flow components give rise to the blue flow directions in the schematic in Fig. 11 (top).
The radial component points inward (∂ ρ c < 0) because of the radially decaying surfactant concentration.This gives rise to to radially outward Marangoni flows.Because −e z • ∂ ρ ce r ∝ cos θ in eq. ( 15), this increases the direct force in the front (around θ = 0) but decreases it in the back (around θ = π).Advection leads to bigger surfactant concentrations in the rear, which also result in bigger radial concentration gradients on the rear side and lead to an overall decrease of the direct force FM,fl < 0 and, thus, an increased drag.This is the dominating effect in the advective regime Ū 1, where a concentration boundary layer of width ∆ρ ∼ Ū −1/3 forms around the half-sphere, which results in steep radial concentration gradients that are stronger on the rear side.The stronger radial Marangoni flow components in the rear are indicated by the blue arrows in Fig. 11 (bottom); this phenomenon is also visible in the experimental PIV measurements in Fig. 6(B) during motion.There are also slightly bigger radial concentration gradients in the diffusive regime Ū 1, but they are sub-dominant for the slow radial decay of the function c1 (ρ) in the absence of a radial concentration boundary layer.

Strong Marangoni flow Pe Ū
For weak Marangoni flow, Pe Ū , we could decouple the problem and to regimes, a diffusive or linear response regime for Ū 1 and an advective regime for Ū 1.Now we increase the Peclet number Pe and, thus, the Marangoni flow.For a strong Marangoni flow, Pe Ū , the linear regime Ū 1 becomes modified.We first have to consider the dominant Marangoni flow problem (iib), which determines the main component of the fluid flow in the diffusion-advection problem (iii).The Marangoni flow pattern is a stationary Marangoni vortex ring around the spherical swimmer below and parallel to the fluid interface S Int .Because this solution lacks axisymmetry an analytical solution is no longer possible.Nevertheless, we can obtain novel scaling results for concentration profile and Marangoni flow field in a concentration boundary layer of width l c below the fluid interface S Int along similar lines as Refs.[59,60].
The surfactant is emitted from the sphere with the fixed current density j| S = 1.Advection by the Marangoni flow vM concentrates the surfactant in the boundary layer of width lc below S Int .It takes a time t ∼ r/v M to reach a radial distance r.During this time, the surfactant diffuses over a distance l c ∼ (Dt) 1/2 ∼ (Dr/v M ) 1/2 or lc ∼ (ρ/v M ) 1/2 in vertical y-direction, which sets the boundary layer width lc .Because we are at low Reynolds numbers, the laminar boundary layer below the the fluid interface S Int is of the size δ ∼ a ( δ ∼ 1) set by the sphere.The laminar boundary layer governs the decay of the Marangoni flow field v M in y-direction.
Moreover, we have mass conservation of the emitted surfactant.The total advective flow J ∼ 2πcv M ρ lc below the interface at distance ρ will always equal the original flow J = 2π that is emitted at the half-sphere, In addition, the Marangoni boundary condition (see eq. (20d) gives a second equation for concentration and velocity, which follows from the definition of the laminar boundary layer width δ.Combining both eqs.( 34) and ( 35), we find a differential equation for c(ρ), which we solve with the boundary condition c(∞) = 0 resulting in We see that the advective current jM ∼ cv M ∼ Pe 1/3 ρ −7/3 becomes smaller than the corresponding diffusive current jD ∼ −∂ ρ c ∼ Pe −1/3 ρ −5/3 for ρ > Pe.Then our assumption of advective transport breaks down, and this should mark the boundary of the Marangoni advection dominated region.Therefore ρ M ∼ Pe should be the scaling of the size of the Marangoni vortex around the sphere for low Reynolds numbers.At larger distances, a crossover to diffusive transport with c ∝ ρ −1 sets in.So far, we considered the leading order of our problem by setting Ū ≈ 0; going one order further, we get the linear response for small Ū with the ansatz c = c(0) + Ū c(1) with c(0) (ρ) given by (36).In the total flow v + v M , the Marangoni flow (37) is the zeroth order result, v M = v (0) M , while the Stokes swimming flow v = v (1) is linear in Ū .In an advection dominated situation, the constant flux relation (34) still holds in the presence of Stokes flow, where ū(ρ) is the radial component of the Stokes flow.
Expanding up to first order in Ū and using (34) for the leading order, we find This contribution is symmetry-breaking; inserting this scaling of the concentrations into the relations ( 23) and ( 24for the Marangoni forces, we obtain We checked these predictions numerically Fig. 12 by using our iterative FEM approach (see Appendix A), which is possible up to Pe ∼ 50 and find good agreement, in particular, for the predictions FM /Pe ∝ Pe −2/3 and FM,tot /Pe ∝ Pe −2/3 , which will be most important for the swimming speed relation (see insets in Fig. 12).Moreover, these numerical FEM results show that both prefactors in eq. ( 40) are of order unity but hard to quantify because of finite size effects.This shows that Marangoni flows depress the total driving force in the linear response regime by a factor Pe −2/3 reflecting the fact that it is harder to break the symmetry in the presence of the strong Marangoni flow advection.Numerical results in Fig. 12 also show that the total Marangoni force is somewhat larger than the direct Marangoni force, FM,tot > FM .In this respect, our previous results for linear response regime for Pe Ū remain unchanged: the Marangoni flow force increases the direct force.
In the advection-dominated regime Ū 1, on the other hand, results are essentially not affected by increasing the Peclet number Pe into the regime of strong Marangoni flows Pe Ū , as the numerical results in Fig. 12 show: all curves for the Marangoni forces converge to the previous diffusion-advection results for Ū 1.The reason for this behavior is that the flow field v will still give rise to a concentration boundary layer of thickness ∆ρ ∼ Ū −1/3 around the sphere, which determines the concentration

Evaporation
In the presence of evaporation, the boundary condition for the diffusion-advection problem changes at the air-water interface S Int .We then have the convective (Robin) boundary condition (20b), which is governed by the dimensionless Biot number (22), instead of the Neumann condition (20a), which is recovered for vanishing Biot number k = 0.In general, evaporation of surfactant depletes the interface of surfactant and, thus, decreases the Marangoni driving forces (both direct and flow force).
For volatile camphor, k ≈∼ 10 −4 m/s has been suggested [26], which corresponds to a high Biot number k = ak/D ≈ 550 for the camphor disks from [28].PEG, on the contrary, has a negligible Biot number as it is not volatile.As a consequence of the new convective boundary condition, the concentration profile will fall off exponentially perpendicular to the interface in the outward direction on a dimensionless extrapolation length scale ∆ȳ ∼ 1/ k given by the inverse of the Biot number.
In the presence of evaporation, we can develop a qualitative scaling theory based on the assumption that the total evaporation flux balances the total emission flux of surfactant in a stationary state, which gives in dimensionless quantities and determines the derivatives ∂ ȳ c(ρ)| ȳ=0 at the scaling level.Via the convective boundary condition (20b), this also determines the surface concentration c(ρ)| ȳ=0 .Moreover, the convective boundary condition should reduce to our previous results for the Neumann condition for small Biot-numbers k, where the evaporation flux jev = k c(ρ)| ȳ=0 is smaller than the dominating transport flux, which is the diffusive or Marangoni flux for Ū 1 and the convective flux for Ū 1.For the diffusion or Marangoni dominated situation for Ū 1, the concentration and, thus, evaporation is distributed over the whole interface S Int , i.e., there is no concentration boundary layer around the half-sphere.Therefore, flux balance (41) leads to ∂ ȳ c(ρ)| ȳ=0 ∼ O(1) resulting in c(ρ)| ȳ=0 ∼ 1/ k because of the convective boundary condition (20b).Then, also FM ∼ 1/ k and FM,tot ∼ 1/ k.Moreover, the evaporation flux dominates over the diffusive or Marangoni fluxes (which are O(1)) only for k > k0 with a crossover value k0 = O(1), which determines the crossover to the non-evaporative case.Our numerical results in Fig. 13 suggest k0 1.Therefore, we expect We checked these predictions numerically in Fig. 13 by using our iterative FEM approach (see Appendix A) and find good agreement.The plots in the bottom row (yellow symbols) show that the dependence on k for Ū 1 is described very well by FM = FM k=1 / k and FM,tot = FM,tot k=1 / k suggesting k0 1.For the advection-dominated situation for Ū 1 the concentration and, thus, evaporation is present only in the concentration boundary layer of radial width ∆ρ ∼ Ū −1/3 around the half-sphere and in the advection tail.Therefore, the flux balance (41) leads to ∂ ȳ c(ρ)| ȳ=0 ∼ Ū 1/3 and c(ρ)| ȳ=0 ∼ Ū 1/3 / k because of the convective boundary condition (20b).Then, also FM ∼ Ū 1/3 / k and FM,tot ∼ Ū 1/3 / k if the evaporation flux dominates over the convective fluxes.The convective flux at the interface and at the boundary layer (ρ ≈ 1 + ∆ρ) is in radial direction ju = Ū ū(ρ) cos θc(ρ) ∼ Ū ∆ρ 2 c(ρ) ∼ Ū 1/3 c(ρ) and jv = Ū v(ρ) sin θc(ρ) ∼ Ū ∆ρ sin θc(ρ) ∼ Ū 2/3 sin θc(ρ) in θ-direction.In the advection tail (of angular width ∆θ ∼ Ū −1/3 ), this also leads to jv ∼ Ū 1/3 c(ρ).Therefore, the evaporation flux starts to dominate over the convective radial flux and the flux in θ-direction in the advection tail for k > Ū 1/3 ; only then we see the effects of evaporation.Therefore, we expect for Ū 1.We checked these predictions numerically in Fig. 13 and find good agreement.The plots in the bottom row (blue symbols) show that the dependence on k for Ū 1 agrees very well with eq. ( 43).In summary, we see a reduction of all Marangoni forces by evaporation both in the linear response regime Ū 1 but also in the regime Ū 1 of strong advection.In both regimes, evaporation reduces the surfactant concentration, which decreases the Marangoni forces.

Force balance and swimming condition
Now we have a rather complete picture of the solution of problems (i)-(iii), i.e., diffusion-advection coupled to hydrodynamics for a prescribed swimmer velocity Ū at low Reynolds numbers.The main result are the Marangoni forces as a function of the prescribed velocity Ū .The swimming condition (18) or (25) gives an additional force balance relation between Marangoni forces and Ū , which has to be satisfied in the swimming state and determines the swimming speed Ū = Ūswim as a function of Peclet number Pe and Biot number k.In general, the swimming velocity increases with Pe and decreases with k.
First, we consider small Pe, i.e., small surfactant emission rates and see whether a swimming state with spontaneously broken symmetry can exist.For k ≈ 0, as appropriate for the PEG-alginate swimmer, we find from eq. ( 30) that a solution for the swimming condition exists above a critical Peclet number Pe > Pe c ∼ 8/ ln R → 0, which approaches zero for large system sizes.Therefore, the symmetry is essentially always spontaneously broken in a large swimming vessel.Spontaneous symmetry breaking resulting in propulsion is possible by establishing an asymmetric surfactant concentration profile that is maintained by advection and can produce the necessary Marangoni forces.Equation ( 30) is valid only in the decoupled limit Pe Ū .At the swimming bifurcation, we have Pe = Pe c Ū ≈ 0, however, such that the feedback of Marangoni flows onto the diffusion-advection problem has to be taken into account, and the decoupling approximation cannot be used.Then, eq. ( 40) describes the Marangoni forces in the linear response regime, which further reduces the critical Peclet number to Pe c ∼ 1/(ln R) 3 → 0. In the presence of relevant evaporation k 1, as appropriate for camphor, the total Marangoni force is depressed according to eq. ( 42) Fig. 13.Iterative three-dimensional FEM results for FM/πPe (top) and FM,tot/πPe (middle) as a function of Ū for Pe = 0 and Biot numbers k = 0 − 400 for a half-cylindrical system with ρ < 15, x > 0, −4 < ȳ < 0, Blue dashed lines are results for k = 0 from FEM solutions to the axisymmetric diffusionadvection equation in two-dimensional angular representation with R = 30.Forces in the diffusive regime Ū 1 are reduced according to eq. (42).Forces in the advective regime Ū 1 are reduced according to eq. ( 43).Bottom row: FEM results for FM/πPe and FM,tot/πPe as a function of Biot number k for Ū = 10 −4 , 100 in comparison to scaling results in eqs.( 42) and (43).
resulting in an increased Pe c ∼ k3 /(ln R) 3 → 0, which is, however, still approaching zero for large swimming vessel sizes R. Immediate propulsion in all experiments is in accordance with a bifurcation with a small Pe c .Moreover, we observe an intermittent stopping of the swimming motion only in the very end (after 20 min or more) before the the swimming motion stops completely (because the fuel has been consumed).This confirms a small critical value Pe c below which Pe drops only for very small emission current densities α.Small irregularities can already break the symmetry and give rise to an avoided bifurcation and select a fixed swimming axis with respect to the particle orientation, which is also observed in the experiments.
For Pe > Pe c , a spontaneously symmetry-broken swimming state with Ūswim > 0 exists.Because the Marangoni force eq. ( 30) remains approximately linear up to Ū ∼ O(1), as can also be seen in Fig. 10, the swimming velocity rises steeply for Pe Pe c and quickly enters the asymptotics for the advection-dominated regime Ūswim 1.Here, we find the swimming relations Ūswim ∼ Pe 3/5  for k Pe 1/5 , (44) Ūswim ∼ k−3/4 Pe 3/4  for k Pe 1/5 .
Also in this regime, we have Pe Ūswim such that Marangoni flows are strong, but this has little influence on the swimming speed because of the concentration boundary layer that forms in the advective regime.Evaporation is significant for k Pe 1/5 and reduces the swimming speed because it reduces the driving Marangoni forces.The swimming relations (44) and (45) are shown in Fig. 14 as dotted yellow and dotted blue lines, respectively, together with the experimental results for our PEG-alginate swimmers and camphor boats from Boniface et al. [28].We see clearly, that the experimentally observed swimming speed differs, because these swimmers operate at higher Reynolds numbers.

High Reynolds numbers
We have developed a complete picture of the solution of problems (i)-(iii), i.e., diffusion-advection coupled to hydrodynamics for a prescribed swimmer velocity Ū at low Reynolds numbers, including evaporation.Low Reynolds numbers Re = 2 Ū /Sc 1 are realized for Ū Sc/2, which can still be much larger than unity for typical Schmidt numbers for surfactants in aqueous solutions (see Tab. 2).For the relevant Marangoni propulsion forces, the following picture has emerged from our analysis at low Reynolds numbers.There is a diffusive regime for Ū 1, which becomes modified by strong Marangoni flows for Peclet numbers Pe Ū , and there is an advective regime for Ū 1, which is essentially unchanged in the presence of strong Marangoni flows for Pe Ū (see Fig. 12).Both regimes are modified in the presence of evaporation if the Biot number is k ≥ 1 in the diffusive regime and of k Ū 1/3 in the advective regime (see Fig. 13).High Reynolds numbers occur for large velocities Ū Sc/2 and, therefore, always deep in the advective regime Ū 1.At low Reynolds numbers, the concentration boundary layer of dimensionless width ∆ρ ∼ Ū −1/3 (see eq. ( 31)) determines the results for the Marangoni forces in this advective limit (see eqs. ( 27), ( 28), ( 29), ( 30), (40), and ( 43)).
In order to generalize to higher Reynolds numbers, we realize that the concentration boundary layer width is closely related to the Nusselt number.By definition (26), Nu = −∂ ρ c0 (1)/c 0 (1), the Nusselt number is an inverse extrapolation length, which we expect to be the inverse concentration boundary layer width, The result Nu ∼ Ū 1/3 from eqs. ( 27) and ( 28) confirms this result both for constant flux and constant concentration boundary conditions at low Reynolds numbers, and we conjecture it to hold also at higher Reynolds numbers.Phenomenologically, the Nusselt number is well-studied also for high Reynolds number [61], both for heat (Nu T in the following) and for mass transfer (Nu in the following, also Sherwood number Sh in the literature), and we can draw on these results in order to develop a theory for the concentration boundary layer and the Marangoni forces.Up to moderate Reynolds numbers Re 200, the physics is governed by additional laminar (viscous) boundary layers that appear around a sphere in fluid flow, which typically have a width δ ∝ Re −1/2 [62,63].The viscous boundary layer scaling can be rationalized by generalizing our above scaling argument for the concentration boundary layer leading to eq. (31).The important difference is that the velocity field close to the sphere changes from v(∆r) ∼ U ∆r/a for Stokes flow to v(∆r) ∼ U ∆r/δ for laminar boundary layer flow with a no slip boundary condition.This leads to This scaling result is in accordance with phenomenological results for the Nusselt number by Ranz and Marshall [64] (Re = 2 Ū /Sc) (Re = 2 Ū /Sc and with Sc replacing the Prandtl number Pr for the mass transfer Nusselt number).
Because the concentration will drop in radial direction from its value at the surface S of the half-sphere to zero within the concentration boundary layer of width ∆ρ, and, thus, 1 = −∂ ρ c(ρ = 1) ∼ c(ρ = 1)/∆ρ for constant flux boundary conditions, the scaling of the concentration boundary layer width (46) also gives rise to i.e., the symmetry-breaking concentration level at the sphere is inversely proportional to the Nusselt number (generalizing eq. ( 32)).Therefore, the direct Marangoni force (23) at higher Reynolds numbers.The total Marangoni force does no longer follow from a reciprocal theorem.In terms of an energy balance, the reciprocal theorem can be interpreted as the absence of mutual dissipation between swimming flow and Marangoni flow [58].Therefore, the power input by surface Marangoni stresses into the swimming flow is, transmitted without loss as power input by the Marangoni flow force onto the swimmer.For higher Reynolds numbers, the mutual dissipation is no longer zero but there are additional viscous terms appearing, which are connected to the vorticity of the flow.This suggests that the Marangoni stresses at the interface become less effective in generating a Marangoni flow force because of this additional dissipation.Therefore, we simply neglect the Marangoni flow force (or assume that the Marangoni flow force is sub-dominant) and only consider the direct Marangoni force (50) at high Reynolds numbers, in the following.
Likewise, the existence of viscous boundary layers around the half-sphere modifies the drag force.On phenomenological grounds, it has been suggested that F D = D c π 2 µaU with D c 6Nu T [65], where Nu T is the Nusselt number for heat transport, resulting in Using the Ranz and Marshall correlation (48), we find from the force balance FD + FM = 0 Pe ≈ 3NuNu T Ūswim , Ūswim ∼ Sc 1/3 Pr −1/6 Pe 1/2 (52) in the absence of evaporation.In the presence of evaporation, we use FM ∼ FM k=0 Nu/( k + Nu) (cf.eq. ( 43)) to find Both results ( 52) and ( 53) are also shown in Fig. 14 together with the experimental data on PEG-alginate and camphor Marangoni boats.

Comparison with experiment
Force balance for low and high Reynolds numbers results in a characteristic Pe-Ūswim -relation for the swimmer with characteristic power laws.For low Reynolds numbers, these are relations (44), Ūswim ∝ Pe 3/5 , without evaporation and (45), Ūswim ∝ Pe 3/4 , in the presence of strong evaporation.For higher Reynolds numbers, we find relations (52), Ūswim ∝ Pe 1/2 without evaporation and ( 53), Ūswim ∝ Pe 2/3 for strong evaporation.Also the experiment on PEG-alginate swimmers and the camphor boats from Boniface et al. [28] take place at higher Reynolds numbers; parameter value estimates for these experiments are summarized in Tab. 3, the resulting dimensionless parameters in Tab. 2. Our experimental results for the mass release mass ṁ(t) as a function of time (see Fig. 7(middle)) and the corresponding swimming velocity U swim (t) (see Fig. 7(right)) of the PEG-alginate swimmers give the red line in Fig. 14 in the Pe-Ūswim parameter plane.We also show experimental results for camphor boats from Boniface et al. [28] (black data points from experiments varying the radius and black line from time-dependent swimming data).Figure 14 compares these experimental results with our theoretical results for the appropriate parameter values for the PEG-alginate swimmers and the camphor boats from Tabs. 2 and 3.The swimming relations in Fig. 14 are the main result of the paper.For the PEG-alginate swimmer (red line) we see good agreement between the high Reynolds number theory in the absence of evaporation, i.e., with Biot number of k = 0 (yellow line).The corresponding low Reynolds number theory (dotted yellow line) gives significantly lower swimming velocities at the same Peclet number.This outcome is what we expected based on the above estimate of moderate Reynolds numbers Re ∼ 60 for the PEG-alginate swimmers and based on the non-volatility of PEG.We also see that the slower second phase of the swimming motion of the PEG-alginate swimmers is described slightly better by our theory (left part of the red line in Fig. 14), which is in accordance with our initial observation that the time constants for mass release and swimming velocity agree only in the second phase.In the first phase, the Peclet number Pe necessary to achieve the measured swimming velocity is slightly lower than predicted by our theory, i.e., a more efficient propulsion.This is a hint that some of our theoretical assumptions could be violated during the first phase, for example, regarding the adsorption equilibrium, which might not yet be established starting from an initially "empty" air-water interface, which could give rise to steeper concentration gradients and more efficient propulsion.
For volatile camphor, a Biot number of k ≈ 550 has been suggested in Ref. [26], which we use in Fig. 14 to compare with the experimental data of Boniface et al. on camphor disks [28] (black data points and black line).Again, we obtain good agreement with the high Reynolds number theory (blue line).The disk geometry differs from the half-spherical geometry we discussed in detail, but we expect that the swimming relation will only differ by numerical factors of order unity.The corresponding low Reynolds number theory (dotted blue line) significantly underestimates swimming velocities, and a theory without evaporation (dashed blue line) overestimates swimming velocities.

Discussion and Conclusion
We presented an experimental realization of alginate capsule self-propulsion at the air-water interface by loading the alginate capsule with surfactant molecules during synthesis.Self-propulsion of these capsule swimmers is based on a Marangoni boat mechanism.Alginate is bio-compatible and widely used for capsule production, which are interesting aspects for further applications.The versatile and for PEG-300, we find a fast and sustained motion with swimming speeds U swim ∼ 2 − 3 cm/s over 20 min and more.The swimming speed corresponds to several swimmer diameters per second and is comparable or superior to other self-phoretic or microswimmers [4] or active liquid droplets [16].In general, we find prolonged propulsion only if spreading molecules are water-soluble as the PEG molecules are; then the air-water interface can regenerate by the fuel being dissolved in water.Evaporation from the air-water interface is another mechanism to achieve regeneration, which is utilized in camphor boats [26][27][28]32].We conclude that a mechanism that regenerates the airwater interface, such as water-solubility or evaporation of surfactants, is crucial for prolonged propulsion.We could produce alginate swimmers down to radii of several hundreds of micrometers, which is slightly above the realm of low Reynolds numbers.Future work could address further miniaturization of capsules.
Starting from low Reynolds numbers, we developed a theory for Marangoni boat propulsion of a completely symmetric, half-spherical, surfactant emitting swimmer.The theoretical description comprises the coupled problems of surface tension reduction by surfactant adsorption at the air-water interface including the possibility of surfactant evaporation, fluid flow (both Marangoni flow and flow induced by swimmer motion), diffusion and advection of the surfactant.In particular, advection is systematically included in our approach and turns out to be essential for all swimmer velocities U D/a ( Ū 1).These three problems are first solved for prescribed swimmer velocity U ; the actual swimming velocity U swim is determined by force balance between the drag force, the direct Marangoni force from the surface tension contribution at the air-waterswimmer contact line, and the Marangoni flow force.We find that Marangoni flows can either act to increase the direct Marangoni force (at low velocities) or to increase the drag (at higher velocities).For low Reynolds numbers, all theoretical results are supported by numerical FEM simulations.Non-dimensionalization shows that the swimmer is controlled by two dimensionless control parameters, the Peclet number (21), which is the dimensionless emission rate of surfactant, and the Biot number (22), which is the dimensionless evaporation rate.Evaporation is practically absent for PEG, but strong for other frequently studied Marangoni boat swimmers, such as camphor boats [26].
We showed that a spontaneous symmetry breaking, i.e., a spontaneous transition into a swimming state is possible also for a completely symmetric swimmer above a critical Peclet number.Spontaneous symmetry breaking resulting in propulsion is possible by establishing an asymmetric surfactant concentration profile that is maintained by advection.We find that the critical Peclet number for this transition approaches zero logarithmically for large system sizes, Pe c ∝ 1/(ln R) 3 .The possibility of such a spontaneous symmetry breaking has been pointed out for autophoretic swimmers [9,18] and liquid Marangoni swimmers [17] before.Also in these systems, advection by the surrounding fluid can maintain the necessary concentration gradients in fields and/or concentrations.
In eqs.( 44) and ( 45), we obtain the power laws governing the swimming velocity as a function of Peclet and Biot number, which are Ūswim ∝ Pe 3/5 , without evaporation (PEG) and Ūswim ∝ k−3/4 Pe 3/4 , in the presence of strong evaporation (camphor).This demonstrates that additional evaporation reduces swimming speed.
Experimentally realizable PEG-alginate or camphor swimmers are operating at moderate Reynolds numbers around 60 or more.Accordingly, we generalized the theoretical approach to higher Reynolds numbers by using the concept of the Nusselt number, for which many results at higher Reynolds numbers are known phenomenologically.This might also account for some effects related to the formation of vortices around the swimmer during propulsion at higher Reynolds numbers (see PIV-results in Fig. 6 and Ref. [29]).Finally, we obtained the swimming relations ( 52) and (53), which give Ūswim ∝ Pe 3/4 , without evaporation (PEG) and Ūswim ∝ k−2/3 Pe 1/2 , in the presence of strong evaporation (camphor).We find a good quantitative fit (without any free fitting parameters) with our own experimental results on PEG-alginate swimmers and the results of Ref. [28] on camphor swimmers in Fig. 14.This is the main result of this paper.Future work should extend the numerical approach to higher Reynolds numbers in order to verify our scaling results for the swimming relation using, for example, the methods introduced in Ref. [37].There are several aspects of the self-propulsion of PEG-alginate capsules, where we presented first experimental results but which deserve a much more detailed investigation in future work: curved trajectories, interactions with container walls, and swimmer-swimmer interaction.
Curved trajectories as observed in Figs. 4 and 5 with a swimming direction of the swimmer that is, at least, weakly linked to its orientation, while the orientation of the swimmer is slowly turning can only be explained by small asymmetries of capsules induced by irregularities in the pore structure.This view is supported by the individual character of the turning characteristics of different swimmers (see Fig. 5).Future work should explore the relation between capsule irregularities and turning statistics in more detail.Experimentally, purely rotary systems could be constructed [42].
PEG-alginate swimmers are repelled by walls.In normal collisions we observe direction reversal without reorientation of the swimmer.In the framework of the Marangoni boat mechanism, this can be explained by an accumulation of surfactant emitted by the swimmer in front of the wall because of the zero flux boundary condition at the wall.Surfactant accumulation creates a gradient in surfactant concentration towards the wall, and the swimmer reverses direction if advection and accumulation balance without changing its orientation.This behavior is similar to what has been observed for asymmetric [24] and symmetric [25] camphor boats [22].During the collision the orientation of the swimmer particle does not change, while the swimming direction reverses; therefore, the swimming direction also reverses relatively to the particle orientation.This is consistent with a weak symmetry breaking by small irregularities in the pore distribution, which give rise to many possible metastable propulsion directions.A perturbation as during surfactant accumulation and direction reversal at the wall can easily cause a change between these propulsion directions.There are more oblique collisions (see Fig. 5), which take longer and can feature a reorientation of the swimmer.The underlying mechanisms could be similar to the reorientation mechanisms of self-diffusiophoretic swimmers [13,14] but this issue also requires future work.
Finally, we have experimental evidence that PEG-alginate swimmers interact with each other via their surfactant concentration fields.Similar observations have been made already in Refs.[22,40,66], mostly in channel geometries.We monitored different collisions between swimmers, where we could find both attraction and repulsion.In the framework of the Marangoni boat mechanisms, where swimmers prefer to move opposite to surfactant concentration gradients, we expect that swimmers are repelled by their surfactant tails, which represent traces of high concentration.This predicts a kind of "chemo-repellent" behavior with respect to the tails.For the interaction between swimmers, the concentration-dependence of the surface-activity κ (the c 0dependence in eq. ( 5)) can also play an important role [40].The topic of swimmer interactions appears to be very rich and important for applications regarding the swarming of PEG-alginate swimmers; it deserves a much more detailed investigation in future work.
6 Authors contributions H.R. and A.-K.F.conceived the experiments.A.-K.F.performed the experiments.J.K. and H.E. developed the theoretical model, performed the analytic calculations and numerical simulations.J.K. wrote the manuscript with support from H.E. and H.R. All authors discussed the results and commented on the manuscript.

A Numerical methods
For Pe Ū , the Marangoni flow problem decouples and the advection problem becomes axisymmetric.Then, c = c(ρ, θ) only depends on the radial coordinate and one angular coordinate and we solve the diffusion-advection problem on a two-dimensional rectangular domain in the ρ-cos θ plane.Typically we use ρ < 30 and FEM-routines from Wolfram MATHEMATICA; we use irregular triangular meshes with a mean area of mesh elements of 0.00015 in the ρ-cos θ plane (the maximal area is 0.005).We checked that discretization effects are negligible.
For Pe Ū , we have to solve the coupled problems of Marangoni flow field (iib) and diffusion-advection equation (iii) with the appropriate boundary conditions from the adsorption problem (i).Numerically, we only address low Reynolds numbers for the fluid flow such that the Stokes equation applies.For the coupled problem we use an iterative scheme of three-dimensional FEM solutions both for Marangoni flow and for diffusion-advection, again employing FEM-routines from Wolfram MATHEMATICA.For the Stokes flow we use the analytical result (eqs.(9a) and (9b)).In the iterative scheme we solve for the Marangoni flow field (iib) starting from an initial guess for the concentration profile (typically c(0) (ρ) = 1/ρ).Then we use this Marangoni flow field in the advecting fluid flow field in the FEM solution v(1) M (ρ) of the diffusion-advection equation (iii), which gives in turn an improved approximation c(1) (ρ) for the concentration profile.With this improved approximation we go back into solving for the Marangoni flow field (iib) to obtain the next iteration v(2) M (ρ), and so on.This results in improving approximations c(n) (ρ) and v(n) M (ρ), and the iteration typically converges within n = 5 − 10 iterations to the final Marangoni flow field and surfactant concentration field.The iterative approach has the advantage that the Marangoni boundary condition in the fluid flow problem (iib) is a fixed one at each iterative step and only adjusts over the iteration; the coupling of the two problems is correctly established over the iteration.Similar iterative numerical schemes for coupled problems have been applied successfully in Refs.[67,68].The iterative scheme works well up to Pe ∼ 50 for the system sizes we use.As outlined in the section on strong Marangoni flow in the limit Pe Ū , Marangoni flows are fully established on the scale ρ M ∼ Pe, which is growing with the Peclet number.Therefore, numerical problems typically arise for higher Peclet numbers when the Marangoni roll becomes strongly distorted by the system boundaries, which gives rise to numerical instabilities and a failure of the iteration procedure.This is why the data presented in Fig. 12 is limited to the range Pe = 0 − 50.
The FEM solution of the stationary equations (iib) and (iii) is obtained on a cylindrical or cubical irregular tetrahedral mesh.We use cubical volumes (for example, with edge length 14 in xz-plane and height 7 in ȳ-direction in Fig. 12) for the FEM calculations.Cylindrical volumes can also be easily implemented.The maximal volume of mesh elements is 0.2, the mean volume is 0.01.Mesh volumes are smaller (< 0.005) in the region −1 < ȳ < 0 below the interface to capture Marangoni advection.Again, we checked that discretization effects are small.Because of the mirror symmetry x → −x, we only need to solve on half-cubes x > 0 and apply Neumann boundary conditions ∂ x c| x=0 = 0 and ∂ x vM | x=0 = 0 to enforce the mirror symmetry.The boundary conditions at the outer boundaries are Dirichlet conditions for the concentration c = 0 and the Marangoni flow vM = 0.For sufficiently large cubes or cylinders these boundary conditions should not matter but we still have finite size effects, in particular at larger Peclet numbers Pe > 50 as explained above.In particular, our systems are always several particle radii long in every direction such that strong confinement effects as observed in Ref. [37] for system sizes below two particle radii should be absent.
We are interested in the resulting symmetry-breaking Marangoni forces caused by a symmetry-breaking swimming motion as a function of the velocity Ū .At small Ū there is the problem that artificial symmetry breaking from lattice irregularities/defects is often larger than symmetry breaking by swimming.Therefore we average all measured quantities over two simulations with Ū and − Ū to cancel artificial symmetry-breaking effects.This step is crucial to obtain accurate results.

Fig. 1 .
Fig. 1.Schematic of the PEG-alginate capsule.The watersoluble "fuel" or spreading molecule PEG is incorporated during alginate capsule synthesis in the core and diffusively emitted during swimming.

Fig. 2 .
Fig. 2. Synthesis of PEG-alginate swimmers by pipetting microliter amounts of PEG-alginate solution into crosslinker solution.Side view of a PEG-alginate swimmer showing its half-spherical shape.

Fig. 4 .
Fig. 4. Typical swimming trajectory of a PEG-alginate swimmer in the crosslinker solution.Color-coded is the swimming velocity Uswim, the trajectory shows the first 84 s of swimming.

Fig. 6 .
Fig. 6.PIV measurements of a PEG-alginate swimmer (A) after dripping and (B) in motion.The velocity scale bar refers to the fluid velocities, which are also color-coded for velocity.The white dashed arrows indicate the direction of the particle motion.On the right we indicate the fluid motion in (B) featuring radial Marangoni flow (blue), combined with flows corresponding to two counter-rotating vortices created by particle motion (green), and downstream Karman-like vortices (red).

Fig. 7 .
Fig. 7. Left: Mass of the PEG-alginate swimmers as a function of time with an exponential fit m(t) = m∞ + m0 exp(−t/τm) (see text); error bars (light blue) denote the standard deviation.Middle: Resulting mass loss − ṁ as a function of time.Right: Corresponding velocity of the swimmer averaged over 10 swimmers; error bars (light blue) denote the standard deviation.The orange line is a fit Uswim(t) = u0 exp(−t/τu,0) + u1 exp(−t/τu,1) (see text) motivated by the existence of two swimming phases.

Fig. 8 .
Fig. 8. Side view (top) and top view (bottom) of the halfspherical Marangoni swimmer geometry with surfactant concentration field c(r) and coordinates.

Fig. 9 .
Fig. 9. Average Nusselt number as a function of Ū for constant flux and constant concentration boundary conditions in the decoupled limit Pe Ū .All results are from numerical FEM solutions of the axisymmetric diffusion-advection equation in two-dimensional angular representation with ρ < R = 30.

Fig. 11 .
Fig. 11.Schematic of concentration field (as concentration contour lines) and Marangoni flow field (arrows) and resulting direct force FM and Marangoni flow force F M,fl in the diffusive ( Ū 1) and advective ( Ū 1) regime.The dominant Marangoni flow contributions are marked in blue.In the diffusive regime Marangoni flows increase the direct force, in the advective regime they decrease the direct force.

Fig. 12 .
Fig.12.Iterative three-dimensional FEM results for FM/πPe (top) and FM,tot/πPe (bottom) as a function of Ū for Pe = 0 − 50 for a cubic system with −7 < ȳ < 0, 0 < x < 7, −7 < z < 7. Insets show the slopes FM/ Ū πPe and FM,tot/ Ū πPe as a function of Pe calculated from the results for Ū = 0.1.Artificial symmetry breaking from lattice irregularities/defects is prevented by averaging all measured quantities over two simulations with U and −U .Blue open circles are results for Pe = 0 from FEM solutions to the axisymmetric diffusionadvection equation in two-dimensional angular representation with R = 30.The slope in the linear response regime for Ū 1 is reduced according to eq. (40).Results for Ū 1 are essentially not affected by strong Marangoni flows Pe Ū .

Table 1 .
Fuel substances leading to successful alginate capsule propulsion

Table 2 .
Dimensionless parameters.Re or Ū , Sc, Pe, and k are control parameters of the problem.ReM and Nu cannot be independently controlled but characterize the resulting solutions; the swimming velocity Ūswim is determined by the force balance swimming condition.

Table 3 .
Estimates of experimental parameters.