From diffusive mass transfer in Stokes flow to low Reynolds number Marangoni boats

We present a theory for the self-propulsion of symmetric, half-spherical Marangoni boats (soap or camphor boats) at low Reynolds numbers. Propulsion is generated by release (diffusive emission or dissolution) of water-soluble surfactant molecules, which modulate the air–water interfacial tension. Propulsion either requires asymmetric release or spontaneous symmetry breaking by coupling to advection for a perfectly symmetrical swimmer. We study the diffusion–advection problem for a sphere in Stokes flow analytically and numerically both for constant concentration and constant flux boundary conditions. We derive novel results for concentration profiles under constant flux boundary conditions and for the Nusselt number (the dimensionless ratio of total emitted flux and diffusive flux). Based on these results, we analyze the Marangoni boat for small Marangoni propulsion (low Peclet number) and show that two swimming regimes exist, a diffusive regime at low velocities and an advection-dominated regime at high swimmer velocities. We describe both the limit of large Marangoni propulsion (high Peclet number) and the effects from evaporation by approximative analytical theories. The swimming velocity is determined by force balance, and we obtain a general expression for the Marangoni forces, which comprises both direct Marangoni forces from the surface tension gradient along the air–water–swimmer contact line and Marangoni flow forces. We unravel whether the Marangoni flow contribution is exerting a forward or backward force during propulsion. Our main result is the relation between Peclet number and swimming velocity. Spontaneous symmetry breaking and, thus, swimming occur for a perfectly symmetrical swimmer above a critical Peclet number, which becomes small for large system sizes. We find a supercritical swimming bifurcation for a symmetric swimmer and an avoided bifurcation in the presence of an asymmetry.


Introduction
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. An important class of low Reynolds number swimming strategies generates interfacial fluid slipvelocities at the swimmer surface, which then lead to self-propulsion because the swimmer must be force-free. This class of swimming strategies comprises phoretic and Marangoni mechanisms. Phoretic mechanisms selfcreate gradients in concentration (self-diffusiophoresis) or temperature (self-thermophoresis) [1,2] which, in turn, give rise to interfacial fluid flow in a thin interaction layer [3].
There are two types of swimmers based on the Marangoni effect [4]: droplet swimmers with liquid interfaces, which can operate in the bulk and solid Marangoni boats or surfers operating at a liquid-air interface. The liquid droplet swimmer is fully immersed in a liquid that carries surfactant. Propulsion is generated by the Marangoni effect, which creates a slip a e-mail: jan.kierfeld@tu-dortmund.de (corresponding author) velocity from a surfactant concentration gradient along the entire liquid-liquid interface between swimmer and surrounding liquid. One typical mechanism to maintain such a surfactant gradient is that more surfactant is adsorbed at the front (in swimming direction) of the swimmer, which depresses the interfacial tension in the front [5][6][7]. In Ref. [8], an auto-diffusiophoretic mechanism coupled to advection [9,10] has been proposed to maintain the surfactant concentration gradient. This propulsion mechanism based on the Marangoni effect is utilized in different liquid Marangoni swimmers, for example, active liquid droplets or active emulsions [6], such as pure water droplets in an oil-surfactant medium (squalane and monoolein) [8] or liquid crystal droplets in surfactant solutions [6]. Many liquid Marangoni swimmers are spherically symmetric initially and swimming spontaneously breaks this symmetry. Beyond the instability, advection and/or preferred adsorption can produce sufficiently strong surfactant concentration gradients and swimming velocities to maintain advection and/or preferred adsorption [5,7,9]. Also, asymmetric shape changes can give rise to concentration gradients and sufficient swimming velocities to maintain asymmetric shapes [11,12]. Here, we consider Marangoni boats or surfers, which employ a different propulsion mechanism. Important examples are soap or camphor boats which have a long history [13]. The crucial difference to liquid Marangoni swimmers is that these boats or surfers operate at a liquid-air interface rather than in the bulk of the liquid. Propulsion is not caused by surfactants that are anisotropically distributed along the swimmer-liquid interface but by the anisotropic distribution of surfactant at the liquid-air interface along which the soap boat propels [14]. The surfactant molecules at the liquid-air interface are emitted or dissolved from the swimmer; this can be achieved by depositing them on the floating swimmer initially [15], by soaking the swimmer in surfactant [16][17][18][19][20][21], or by using a swimmer body made from dissolving surfactant [22]. There are many examples based on DMF (dimethylformamide) [23], alcohol [15,21], soap [21], camphor [16][17][18][19][20]24] or camphene [22] that have also been investigated quantitatively. In a companion paper [25], we discuss alginate capsules as versatile interfacial Marangoni swimmers working with many surface tensions reducing "fuels" in detail, in particular, polyethylene glycol (PEG)-loaded alginate capsules.
So far, Marangoni boats can be produced down to radii a ∼ 150 µm [25], and quantitative results are available down to a ∼ 1500 µm with Reynolds numbers Re ∼ 60, which is still above the low Reynolds number regime. Miniaturization is approaching the low Reynolds number regime, which is the regime we address in detail in the present paper. In a companion paper [25], we discussed low Reynolds number results more briefly and aimed to generalize to high Reynolds numbers using the concept of the Nusselt number in order to describe experiments on PEG-alginate capsule swimmers and camphor boats quantitatively. There is a related system of thermal Marangoni surfers propelled by the thermal Marangoni effect, which was successfully realized only recently [26]. Its theoretical description is equivalent to surfactant-driven Marangoni boats with thermal advection-diffusion replacing surfactant advection-diffusion. Because thermal diffusion coefficients are much higher and swimmer radii reach down to micrometers, this system operates at low Reynolds numbers.
The surfactant molecules are emitted or dissolved from the Marangoni boat, diffuse and advect to fluid flow in the water phase and adsorb to the air-water interface, eventually in interplay with evaporation for volatile surfactants. This creates surface tension gradients and Marangoni stresses on the fluid. Surface tension gradients give rise to a direct net propulsion force (direct Marangoni force in the following). Marangoni stresses on the fluid give rise to symmetrybroken Marangoni flows, which also contribute to (or impede) propulsion via hydrodynamic drag onto the swimmer (Marangoni flow forces in the following). Direct Marangoni forces propel into the direction of higher surface tension, i.e., lower surfactant concentration along the air-water-swimmer contact line. We note that this is opposite to the propulsion in the direc-  [5][6][7][8][9].
A full quantitative theory of Marangoni boats including hydrodynamics, surfactant advection, direct Marangoni forces and Marangoni flows is still elusive but numerical approaches exist [14,18,27,28]. Theoretical approaches ignore the advection of the surfactant concentration field [29][30][31] or even ignore hydrodynamic flow fields [16,17,24,32] or approximate it by uniform flow [20], which clearly oversimplifies the description of surfactant transport. Here, we focus on low Reynolds numbers as in Refs. [27,[29][30][31] and consider a halfspherical swimmer geometry (see Fig. 1), which can simplify the theoretical treatment because axial symmetry can be exploited in certain limits. Experimentally, half-spherical swimmers can be fabricated using the PEG-alginate system [25]. In the limit of weak Marangoni flows, for example, the fluid flow reduces to the well-known Stokes flow around a sphere. We fully include advection of the surfactant concentration field into our analysis as opposed to Refs. [29][30][31], where disks and spheres propelled by the soap boat mechanism have been considered previously. If advection is ignored, the formation of a concentration boundary layer at higher velocities, as it is well known from the related problem of mass transfer from a sphere in laminar Stokes flow [33][34][35][36], will be missed. This happens for velocities U a/D, where a is the sphere radius and D the surfactant diffusion constant, and is essential for the resulting swimming velocity, which is the quantity of main interest in this paper.
In order to obtain the Marangoni forces onto the swimmer, we have to calculate the surfactant concentration profile and have to revisit the problem of mass transfer from a sphere in laminar Stokes flow both for constant concentration and constant flux boundary conditions. In particular, analytical results on the relevant flux boundary conditions are missing in the literature. We fill this gap and derive results for concentration profiles and for the angular dependence of the Nusselt number both for isotropic and anisotropic emission. This allows us to calculate Marangoni propulsion forces both in the diffusive and in the advection-dominated regime and provide analytical results. There is a direct Marangoni force, which is propelling the swimmer into the direction of higher surface tension and a Marangoni flow force, as has been worked out by Masoud and Stone [37]. We will unravel, under which conditions the flow force contribution is propelling or dragging the swimmer. We can extend our analysis to large Marangoni propulsion (high Peclet number) and include effects from evaporation by approximative analytical theories. This part of the analysis is also the focus of a companion paper [25]. Therefore, this discussion is shortened in this paper. All analytical results are corroborated by numerical finite element calculations employing a novel iterative approach.
As opposed to previous work [14,18,[27][28][29][30][31], we consider here completely symmetric Marangoni boats with isotropic surfactant emission as motivated by the experiments in Ref. [25]. Swimming is established in a symmetry-breaking bifurcation. Similar swimming bifurcations have been analyzed by Michelin and coworkers [5,[7][8][9][10], but for a different type of swimmer, namely liquid droplet bulk Marangoni swimmers. For symmetric surface swimmers propelled by thermal Marangoni forces, a related symmetry-breaking effect has been observed in Ref. [38], where symmetric microbeads spontaneously circle around a heating laser beam. Our analysis allows us to obtain the swimming velocity of Marangoni boats as a function of Marangoni propulsion strength (Peclet number) and to analyze in detail the nature of the symmetry-breaking swimming bifurcation.

Model
We introduce coordinates such that the origin r = 0 is at the center of the planar surface of the half-sphere, 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, see Fig. 1. 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 parameterized 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 .
The general strategy to determine the swimming speed has been outlined in Ref. [25]. We first prescribe a stationary velocity U = U e z of the swimmer and analyze the following three coupled problems for its 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 evaporation.
(ii) Low Reynolds number fluid flow including both Stokes flow around the half-sphere and additional surfactant-induced fluid Marangoni flow. (iii) Diffusive surfactant release or surfactant dissolution from the swimmer and subsequent diffusion and advection.
We will start from the diffusion-advection problem (iii) in the presence of Stokes flow around a sphere and neglecting Marangoni flow. This will also give new results for the mass transfer from spheres in laminar Stokes flow. We later examine the additional effects of Marangoni flow and evaporation. The fully coupled problems can also be treated numerically. Solving these three coupled problems, we can obtain the Marangoni forces as a function of the prescribed velocity U from the surfactant concentration profile by employing the reciprocal theorem. Finally, the actual swimming velocity U = U swim is determined from the condition of a force-free swimmer, i.e., the force equilibrium between Stokes drag force, direct propelling Marangoni forces from the surface tension gradient along the air-water-swimmer contact line and Marangoni flow forces.
We begin with a short recapitulation of the governing equations [25].

Coupled adsorption, fluid flow and diffusion-advection problems
Regarding sub-problem (i), we use a local and linear relationship for the surface tension reduction (r is an interfacial vector with y = 0) by the local surfactant concentration difference c(r) with respect to a bulk concentration background value c 0 (the concentration at |r| → ∞); κ is a coefficient characterizing the propulsion strength. In formulating Eq. (1) locally, we assumed fast on and off kinetics of surfactant to the interface [39] such that the interfacial concentration Γ (r) is slaved to the bulk and only a passive "reporter" of the bulk subsurface concentration c(r)| y=0 . This is appropriate for water-soluble surfactants but, for example, the opposite limit of what has been considered in Refs. [29,31], where surfactant is strictly confined to the interface. Fast on and off kinetics also implies that an imbalance of flux to and from the interface into the bulk can only arise from an additional evaporating flux from the interface to the gas phase. In the general case including surfactant evaporation from the interface, the balance of fluxes to and from the interface gives where k is the rate constant for evaporation. This provides the boundary condition to the diffusion-advection sub-problem (iii) in the bulk. Regarding the low Reynolds number fluid flow subproblem (ii), 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 The flow field v(r) of an externally pulled half-sphere 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 v(r) =û(r, θ)e r +v(r, θ)e θ witĥ The total flow field v tot (r) also has no-slip boundary conditions on the surface of the sphere and 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. 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 generic Marangoni boats with 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, which act both on v M and v tot .
Surfactant diffusion and advection (iii) will play a central role. Surfactant molecules are emitted from the half-spherical surface S and diffuse in the liquid phase. At the same time, they are advected by the total fluid flow. In the stationary state, the bulk concentration field is governed by the diffusion-advection equation We consider two types of boundary conditions which seem most important for applications: (A) slow diffusional surfactant release on S leading to a constant flux boundary condition or (B) surfactant dissolution from the swimmer or surfactant production by some chemical reaction by the swimmer leading to a constant concentration boundary condition, (B) constant conc.: together with c(∞) = 0 and the no-flux boundary condition at the interface S Int . The surface flux α or the surface concentration c S is assumed to be only slowly changing on the time scales of the fluid flow and the surfactant diffusion and approximated as a constant for the calculation of quasi-stationary fluid flow and concentration fields. We non-dimensionalize sub-problems (i)-(iii) by measuring lengths in units of a, velocities in units of D/a, 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. 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. Therefore, we have to discuss both the diffusive caseŪ 1 and the advective caseŪ 1, even at low Reynolds numbers. 1 In many publications on the diffusion-advection problem, such as Refs. [33][34][35] but also in Refs. [20,28,29,31,40], U is called the Peclet number. Here, we define the Peclet number as Pe ≡Ūα, i.e., by the characteristic veloc-ityŪα = καa/Dμ for constant flux boundary conditions. We define it as Pe ≡Ūc s = κcSa/Dμ for constant concentration boundary conditions.Ūα andŪc s are characteristic velocities, where a typical direct Marangoni force FM ∼ κa 2 ∂rc(r = a) ∼ κa 2 α/D for constant flux boundary conditions or FM ∼ κac(r = a) ∼ κcSa/D for constant concentration boundary conditions is balanced by the typical Stokes drag force FD ∼ μaU . The Peclet number is a dimensionless measure of propulsion strength with these definitions.
ReM and Nu cannot be independently controlled but characterize the resulting solutions; the swimming velocityŪswim is determined by the force balance swimming condition Our dimensionless set of equations for problems (i)-(iii) becomes (B) const. conc.: with the dimensionless Peclet number whereṁ = 2πa 2 α is the mass loss per time of the swimmer. We also introduced the dimensionless Biot numberk governing possible evaporation. From Eq. (9d), 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. In the following, we will address the low Reynolds number regime implying that both Re 1 and Re M 1 such that both flow contributions fulfill the Stokes equation. Via the advection withv(ρ) +v M (ρ), 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 Table 1.

Marangoni forces, energy transduction and swimming condition
The half-spherical swimmer moving at velocity U must be force-free and is subject to three forces. First, there is the drag force, which is given by the standard Stokes drag for a half-sphere, F D = F D e z . In dimensionless form usingF ≡ F/Dμ, this is 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 L around the swimmer at y = 0,  [37]. In "Appendix" A, we discuss the reciprocal theorem in terms of energy transduction and find the result (A.10), which states that the mutual power input by Marangoni stresses via the Stokes flow field is completely transduced via the Marangoni flow force onto the sphere, while the power input by Marangoni stresses via the Marangoni flow field itself is completely dissipated. This energy transduction statement (A.10) is equivalent to the result derived by Masoud and Stone [37] for the Marangoni flow force directly from the reciprocal theorem. In the rest frame of the sphere, we obtain in dimensionless formF wherev(ρ)/Ū is the dimensionless Stokes flow field from (3a) and (3b) (in particular, this is independent ofŪ ) in the sphere frame. The total Marangoni forceF M,tot =F M +F M,fl is obtained by using Eqs. (13) and (14) and the Gauss theorem,F withc M (ρ) ≡ 2 π π 0 dθ cos θc(ρ, θ)| y=0 . The total Marangoni driving force has to be determined from the concentration fieldc(ρ, θ) of surfactant molecules (at the interface y = 0). Note that∇ S ·v(ρ) is the two-dimensional surface divergence of the 3D fluid velocity field; therefore,∇ S ·v(ρ) = 0 in general, although∇ ·v(ρ) = 0 for the 3D divergence of the stationary velocity field. The contribution from a constant velocityŪ e z of the whole fluid (if all the fluid would be dragged along by the particle) exactly cancels the direct Marangoni force in (15), and the velocityv(ρ) in the sphere frame determines the total force.
The sign of the Marangoni flow forceF M,fl determines whether it increases or decreases the direct Marangoni force into the direction of higher surface tension: , as in Refs. [29,37], we findF M,tot /Pe = −πA 1 = 1 2F M /Pe, i.e., the total Marangoni force is half the direct Marangoni force if only the first cos θ-component is relevant. Here, Marangoni flow forces drag and decrease the direct driving force (F M,fl < 0). This result will change as we (i) consider 3D diffusion and (ii) as symmetry breaking is only caused by advection, which can focus the concentration field and lead to higher Legendre components becoming relevant inc(ρ). As opposed to Ref. [29], we will find that both cases are possible. A backward force is found for steep radial gradients in the concentrationc(ρ), which is the case for high velocitiesŪ 1 in the advection-dominated regime, and a forward force is found at low velocitiesŪ 1 in the diffusive regime.
-Advection leads to a tangential e θ -component of ∇ Sc (ρ) pointing from the front to the rear corresponding to an increasing surfactant concentration toward the rear, which gives rise to a forward Marangoni flow ∼ −e θ . Accordingly, this increases the driving force . This effect dominates in the diffusive regime.
A radial e r -component of∇ Sc (ρ) pointing inward corresponding to a radially decaying surfactant concentration and, on the other hand, gives rise to to radially outward Marangoni flows. Because −e z · ∇ Sc (ρ) ∼ e z · e r ∝ cos θ in Eq. (14), this increases the direct force in the front (around θ = 0) but decreases it in the back (around θ = π). Advection gives rise to bigger surfactant concentrations in the back, which lead to bigger radial concentration gradients on the rear side (in some distance from the sphere because constant flux boundary conditions assure uniform radial gradients right at the surface of the sphere). Overall, the radial Marangoni flows in the back are stronger and decrease the direct force or increase the drag (F M,fl < 0). This effect dominates in the advective regime for constant flux boundary conditions (A) and is rather subtle, as can be seen from the fact that it is absent for the constant concentration boundary conditions (B), whereF M,tot =F M,fl > 0 always. Then, the constant concentration at the surface of the sphere leads to smaller radial concentration gradients on the rear side, because the concentration decay is stretched over a larger distance by advection. Then, radial Marangoni flows in the front are stronger and increase the direct force. -The last equality in (15) shows that the effect of including the Marangoni flow contribution is that the total Marangoni forces are dominated by the concentration profilec(ρ, θ) around ρ ∼ 2, where ρ −1 −ρ −3 assumes its maximal value. Concentration boundary layer profiles concentrated around ρ ≈ 1, as we will find for large swimmer velocitiesŪ > 1 in the advection-dominated regime, give a small total Marangoni force (because ρ −1 −ρ −3 ≈ 2(ρ−1) is small), i.e., Marangoni flows decrease the direct Marangoni driving force. -Long-range contributions as, for example, from a long advection tail can be important, even if they are limited to a small angular regime around θ ∼ π as for high velocities. The highest total force is obtained if a long-range − cos θ-component is present in the concentration profile, as we will find for small swimmer velocities; then, Marangoni flows increase the direct Marangoni driving force. This makes a Marangoni swimmer also susceptible to disturbances in its far-field as, for example, induced by other swimmers.
These results for Marangoni forces as a function ofŪ are inserted into the force balance or swimming condition in order to obtain an additional equation whose solution determines the actual swimmer velocityŪ =Ū swim as a function of the remaining control parameters Pe ("fuel" emission) and eventuallyk (evaporation).

Control parameters and parameter regimes
The non-dimensionalization reveals that the coupled problems (i)-(iii) and the Marangoni forces depend on three dimensionless control parameters (see also Table  1): first, the prescribed dimensionless velocity of the swimmerŪ ; second, the Peclet number Pe characterizing the strength α of the surfactant emission, and third, the Biot numberk characterizing the evaporation. A suitable Peclet number can be defined for both constant flux boundary conditions (A) and constant concentration boundary conditions (B). We also see that the Peclet number both controls the strength of the Marangoni flow via Eq. (9d) and the strength of all Marangoni forces. We note, however, thatF M /Pe and F M,tot /Pe still depend onŪ and Pe via the dependence ofc(ρ) on these parameters. Another important finding from non-dimensionalization is that the diffusion-advection problem (iii) with boundary conditions (i) decouples from the Marangoni flow problem (iib) for Pe Ū or Re M Re, where |v M | |v|, and we can neglect v M in the advection term. Then, the concentration profile is only determined by a classic diffusion-advection problem for mass transfer from a sphere in Stokes flow in the case of constant concentration boundary conditions (B) [33][34][35][36], but with unusual constant flux boundary conditions for case (A). It becomes axisymmetric and only depends on U . In this limit, the Marangoni flow field need not to be calculated in order to calculate the total Marangoni force for the swimming condition. This limit will be the starting point of several analytical calculations.
All in all, we have the following regimes for a symmetric Marangoni boat at low Reynolds numbers: -Ū < 1 and Pe < 1: The concentration profile is governed by diffusion, which is slightly perturbed by advection and described by a linear response in diffusion-advection (iii) with respect toŪ and Pe.
Only the linear response inŪ is relevant for symmetry breaking; therefore, the Marangoni flow can be neglected for the swimming problem. Only for Pe Ū , the Marangoni flow decouples from the advection problem and strict analytical analysis is possible. Swimming sets in (starting withŪ = 0) for a critical Peclet number Pe > Pe c ; if Marangoni flows forces are included, we find Pe c 1 and the symmetry-breaking bifurcation takes place within this regime. -Ū < 1 and 1 < Pe < Sc: All fluid flows are still at low Reynolds numbers, but Marangoni flows are relevant. The concentration profile is governed by symmetric Marangoni advection, which is slightly perturbed by a linear response in diffusion-advection (iii) with respect toŪ , which causes symmetry breaking and swimming. -1 <Ū < Sc and Pe < Sc: All fluid flows are still at low Reynolds numbers, but the concentration profile is governed by advection by the swimming flow forŪ > 1. Advection leads to the formation of a concentration boundary layer of width U −1/3 around the half-sphere forŪ > 1. Only for Pe Ū , the Marangoni flow decouples from the advection problem and strict analytical analysis is possible. For Pe >Ū , Marangoni flows are relevant to advection, in principle, but the surfactant is transported away by the swimming flow field via the concentration boundary layer before it can advect to the Marangoni flow field. There are, however, Marangoni flows in the advection tail, which will become relevant then. 1, the Marangoni flow is mostly radial at the interface because the radial concentration profile is only slightly perturbed by advection at the interface; it features a Marangoni roll (vortex) around the swimmer with an upward flow directly around the particle. For increasing Pe, the normalized Marangoni flow fieldv M (ρ)/Pe as plotted in Fig. 2 seems unchanged indicating a Marangoni flowv M (ρ) that is simply proportional to Pe in strength but otherwise independent of Pe.
At high velocitiesŪ 1, the Marangoni flow pattern changes because the concentration pattern develops the typical advection tail. As a result, there forms a vortex pair within the interface plane, which directs Marangoni flow from the tail to the front. In front of the particle, the flow reaches beneath the particle (around z = 5 in Fig. 2) and resurfaces behind the particle. This leads to a slightly distorted Marangoni vortex roll around the particle. Similar vortex patterns (with a vortex pair within the interfacial plane next to the swimmer) have been observed in Ref. [21], however, by particle image velocimetry (PIV) measurements at high Reynolds numbers. Again, for increasing Pe, the normalized Marangoni flow fieldv M (ρ)/Pe seems more or less unchanged in Fig. 2.
At smallŪ 1, the Marangoni flow forceF M,fl > 0 will increase the direct Marangoni force into positive zdirection because there is a net forward tangential component of∇ Sc (ρ) from the symmetry-breaking advection perturbation proportional toŪ ; the radial component of∇ Sc (ρ) increases the drag, but is slowly decaying at smallŪ and weaker.
At highŪ 1, on the other hand, there is a strong radial component in the concentration boundary layer around the swimmer, which increases the drag. This is created by the large radial component of∇ Sc (ρ) in the concentration boundary layer region of sizeŪ −1/3 and leads to a Marangoni flow forceF M,fl < 0 that decreases the direct Marangoni force. This effect is counter-intuitive as the large vortex pair suggests a strong forward Marangoni force on the large scale picture. The strong radial flows directly around the particle (which are stronger in the backward direction and, thus, dragging the particle) are not clearly visible on the larger scale in Fig. 2. The remaining total Marangoni force mainly comes from the net forward motion in the horizontal vortex pairs but will be weaker than the direct force.

Legendre decomposition for the decoupled limit Pe Ū
In the decoupled limit Pe Ū , the diffusion-advection problem becomes axisymmetric. Then,c =c(ρ, θ) only depends on the radial coordinate and one angular coordinate, and we can also employ a decomposition of the concentration field into Legendre polynomials with respect to the angle θ:c(ρ, θ) = ∞ n=0c n (ρ)P n (cos θ). As derived in "Appendix" B, the diffusion-advection equation (9e) only couples coefficientsc n (ρ) to coefficientsc n±1 (ρ) because the Stokes velocity field (3a) and (3b) can be written in terms of n = 1 polynomials only. We find the diffusion-advection equation (9e) in Legendre representation, for n = 0, 1, ..... For smallŪ 1, the Legendre coefficients will scale asc n (ρ) ∼Ū n and truncation of Legendre decomposition becomes an excellent approximation. This is one strategy for analytical progress in the linear response regime. In "Appendix" B, we also show how the Marangoni forces are expressed by the Legendre coefficients of the concentration field.
Both types of boundary conditions are completely isotropic and only n = 0 components are nonzero. We can include traditional soap boats into our description by introducing explicitly symmetry-breaking anisotropic flux components n > 0 into the boundary conditions, such as in the simplest generic case. Then, the soap boat emits preferentially on the lower half θ > π/2 in case (A) or produces surfactant preferentially on the lower half of its surface in case (B) as in a standard asymmetric soap boat. Such symmetry-breaking emission will give rise to an avoided swimming bifurcation.

Full iterative FEM solution
Numerically, we can consider the problems (i)-(iii) without further approximations at low Reynolds numbers, i.e., solve the coupled diffusion-advection problem and the Marangoni flow problem for a prescribed swimmer velocityŪ .
For the coupled problems of three-dimensional coupled diffusion-advection and Marangoni flow, we use an iterative scheme of three-dimensional FEM solutions to both problems, employing FEM-routines from Wolfram MATHEMATICA in a finite cylindrical or rectangular domain. We iteratively solve for the Marangoni flow field (iib) starting from an initial guess for the concentration profile; then, we solve the diffusion-advection equation (iii) with the resulting total flow field, which gives an improved approximation for the concentration profile. With this improved approximation, we go back into solving for the Marangoni flow field (iib) and start an iteration, which should converge 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. [41,42].
The FEM solution of the stationary equations (iib) and (iii) is obtained on a cylindrical or cubical irregular tetrahedral mesh. We use cubical (for example, with edge length 14 inxz-plane and height 7 inȳ-direction in Fig. 9) or cylindrical volumes (for example, with radius 8 inxz-plane and height 4 inȳ-direction in Fig. 2) for the FEM calculations. The maximal volume of mesh elements is 0.2, and the mean volume is 0.01. Mesh volumes are smaller (< 0.005) in the region −1 <ȳ < 0 below the interface to capture Marangoni advection. Because of the mirror symmetryx → −x, we only need to solve on half-cubes and half-cylindersx > 0 and apply Neumann boundary conditions ∂xc|x =0 = 0 and ∂xv M |x =0 = 0 to enforce the mirror symmetry. The boundary conditions at the outer boundaries are Dirichlet conditions for the concentrationc = 0 and the Marangoni flowv M = 0. For sufficiently large cubes or cylinders, these boundary conditions should not matter but we still have finite size effects. In particular, at large Peclet numbers this can trigger numerical instabilities if the Marangoni roll interferes with the system boundary.
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.

Two-dimensional FEM solution and Legendre representation for the decoupled limit Pe Ū
For Pe Ū , we obtain the decoupled limit, where Marangoni flow does not need to be calculated and the diffusion-advection problem becomes axisymmetric. Then,c =c(ρ, θ) only depends on the radial coordinate and one angular coordinate. We can solve the diffusion-advection problem in a two-dimensional angular representation using finite element methods (FEM), i.e., FEM-routines from Wolfram MATHEMATICA.
For a givenŪ , we can also employ the Legendre decomposition (17) of the diffusion-advection equation and calculate all functionsc n (ρ) by solving the resulting coupled ordinary differential equation boundary value problem. We use the MATLAB routine bvp4c for a domain 1 ≤ ρ ≤R = 300 [43,44] with Legendre components up to n = 61. In this way, we obtain all relevant coefficientsc n (ρ) to calculate all Marangoni forces for the force balance.

Diffusion-advection equation in the decoupled limit Pe Ū and mass transfer from a sphere in Stokes flow
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 becausev M can be neglected. We also neglect evaporation in the beginning. This problem is axisymmetric and equivalent to mass transfer from a full sphere in laminar Stokes flow [33][34][35][36], but with unusual constant flux boundary conditions for case (A). Therefore, we first derive new analytical results for concentration profiles and for the angular dependence of the Nusselt number for these boundary conditions, both for isotropic and anisotropic emission from the sphere. In the decoupled limit, the concentration profile only depends onŪ and is independent of Pe. Thus, the dimensionless Marangoni forces (13) and (15) only depend trivially linearly on Pe, butF M /Pe and F M,tot /Pe are independent of Pe as well. This will make analysis of the swimming condition (16) much easier.

Nusselt number
Diffusive release in an advecting 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 [36]. The average Nusselt number becomes Nu = 1 for a quiescent fluid (Ū = 0), where the flow is purely diffusivec 0 (ρ) ∝ 1/ρ; as soon as advection is present (Ū > 0), the current out of the sphere is increased resulting in Nu > 1. The Nusselt number thus measures how much the current out of the sphere is increased by advection over its purely diffusional value. It is an increasing function of the fluid velocityŪ .
The Nusselt number has been originally defined for constant concentration boundary conditions (B), for which the result is well known [33][34][35][36], with a prefactor that can be calculated analytically [33,36]. We address the Nusselt number also for constant flux boundary conditions (A) and find a very similar result (see Fig. 3) where the prefactor 0.65 is determined numerically from the data in Fig. 3. This result will be derived below. As opposed to the case of a constant concentration boundary condition, it is not possible to obtain an analytical result for the prefactor 0.65. Interestingly, the difference between both types of boundary conditions is small. We conclude that the Nusselt number characterizes the mass transport mechanism by the advecting fluid itself and is rather robust with respect to the emission mechanism (diffusive emission, dissolution or production by a chemical reaction at the surface) by which the transported molecules enter the advecting fluid. This is an

Main results for Marangoni forces
For constant flux boundary conditions (A), the main results for the Marangoni forces as a function of a prescribed velocityŪ arē where the numerical constant d M,B is obtained from the numerical results. Numerical results for these boundary conditions are also shown in Fig. 6. The numerical result in Fig. 6 clearly confirms the existence of just two regimes for both types of boundary conditions. At smallŪ 1, the Marangoni forces are linear inŪ for both types of boundary conditions and can be calculated as linear response in a perturbative approach. In this limit, diffusion dominates. ForŪ 1, on the other hand, advection dominates, and a concentration boundary layer forms around the half-sphere. There is a markedly different scaling for the total Marangoni force comparing both types of boundary conditions, which we will explain below. Figure 6 shows that direct and total Marangoni force reach maximal valuesF M ,F M,tot ∼ 0.15 πPe in the crossover regionŪ ∼ 1 between diffusive and advective transport.

Small velocityŪ , perturbation theory
At smallŪ 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 fieldc (0) (ρ) = 1/ρ atŪ = 0 in powers ofŪ . A first approach is a naive perturbation series Ansatz for each Legendre coefficient starting withc (0) 0 (ρ) = 1/ρ andc (0) n>0 (ρ) = 0. It turns out that this will work only in the "inner region" ρ < 1/Ū of a solution, because in the "outer region" ρ 1/Ū , the convection term can no longer be treated perturbatively, regardless how small U is [34]. The problem that arises in performing such a naive expansion is that alreadyc (1) 1 (ρ) does not vanish at infinity as required by the boundary conditions. What can be done, however, is to treat a finite system ρ <R and apply the boundary conditionsc 0 (R) = 0 andc n>0 (R) = 0 as in the numerical approach. The above results (24), (25) and (26) are obtained by this approach. We find excellent agreement between numerics and naive perturbation theory for such finite systems.
In an infinite system, the situation differs because in the "outer region" ρ 1/Ū the convection term can no longer be treated perturbatively [34]. These effects will only occur for system sizesR 1/Ū , which become extremely large in the perturbative regimeŪ → 0 of interest. To address this problem, in Ref. [34], a systematic expansion in inner and outer region and a matching procedure were performed for the constant concentration boundary condition (B), which is posed in typical heat and mass transport problems in laminar flow [33,34,36]. For the constant flux boundary condition (A), such calculations do not exist at the moment. We also adapt this more advanced matching procedure to the constant flux boundary condition (A). In "Appendix" C, we present the details of the perturbative approach, both the naive perturbation theory and the matching procedure. We find that in linear order inŪ , both approaches still agree in the inner region. For the total Marangoni force, there is a contribution ∝Ū lnR stemming from a ρ-integration of a ρ-independent contribution toc (1) 1 (ρ) in naive perturbation theory, see Eq. (25), and (26). In the framework of the matching procedure, this contribution becomes ∝ −Ū lnŪ as matching provides an upper cutoffR ∼ 1/Ū to the otherwise unchanged inner region.
Regardless of whether this contribution is regularized by system sizeR or by the boundary ρ ∼ 1/Ū of the inner region, the log-divergence of this linear contribution in the total Marangoni force is a remarkable result of these calculations. Because the linear term for the direct Marangoni force stays finite, its existence means that the Marangoni flow forces strongly increase the direct force forŪ 1.

Scaling arguments
For largeŪ 1, advection is strong and a concentration boundary layer of width Δr develops around the half-sphere. The width Δr is determined by the distance that a surfactant molecule can diffuse during the time Δt ∼ a/v(Δr/a) (see Eq. (3b)) that it takes to be transported along the sphere by advection: Δr 2 ∼ DΔt.
Because v(Δr/a) ∼ UΔr/a for Δr/a 1 because of the no-slip boundary condition (see Eq. (3b)), we find This is a classic result for the diffusion-advection problem for constant concentration boundary conditions [33,36], but also holds for constant flux boundary conditions.
We also stress that, for both types of boundary conditions, we find i.e., the local Nusselt number can be interpreted as the inverse local boundary layer width, which is also evident from its definition (23) as an inverse decay length if the concentration profile drops exponentially as a function of ρ.
We obtain in the rescaled variables forc =c(ξ, θ) i.e., a parameter-free equation confirming all boundary layer scaling results (28), (29) and (30). For the constant concentration boundary condition (B), the equations (32) can actually be solved analytically by a similarity transformation [33,36], i.e., with an Ansatzc(ξ, θ) = f (ξg(cos θ)) because this boundary condition is compatible to a boundary condition f (0) = 1 for the function f (x). Exact results can be obtained for the functions g(η) and f (x). An immediate consequence of the existence of such a solution is that the local Nusselt number is the inverse of the function g(cos θ), and that g(cos θ) is identical to the boundary layer width at angle θ because the function f (η) is exponentially decaying on a scale of order unity, This confirms the scaling (30) and (31). The exact results for the functions g(η) and f (x) also give the exact asymptotics of the Nusselt number in Eq. (22), Nu = −∂ ρc0 (ξ = 0) ∼ c 0Ū 1/3 with c 0 = 3 5/3 π 2/3 /8Γ (1/3) 0.624572 [34,36].
A similarity transformation is, however, not possible for the constant flux boundary conditions (A) ∂ ξc (0, θ) = −1, which is incompatible with the similarity Ansatzc(ξ, cos θ) = f (ξg(cos θ)). It turns out that we can reformulate the results for constant concentration boundary conditions in terms of a flux balance argument, which can also apply to the constant flux boundary conditions in order to obtain an approximative result for the local Nusselt number.

Flux balance argument for local Nusselt number
Here, we consider the balance of the diffusive flux out of the sphere at ρ = 1 with the advective flux assuming that a boundary layer Δρ 1 exists to which the advective flux is constrained. We also assume that by its definition (23), the local Nusselt number can be interpreted as an inverse decay length, which is to be identified with the boundary layer width Nu(θ) ∼ 1/Δρ(θ), see Eq. (31).
For the flux balance, we consider a volume from θ = 0 up to an angle θ around the sphere ρ = 1. The diffusive outflux from the sphere gives the particle influx into this volume. ForŪ 1, outflux from this volume is dominated by advection in θ-direction, which is limited to the boundary layer of thickness Δρ(θ) = Nu −1 (θ). Both influx and outflux have to balance in a stationary state. In order to show the flux balance explicitly, we integrate on both sides of equation (32) (for the unrescaledc rather thanc). The integrated diffusive term on the left hand side gives the diffusive influx The integrated advective term on the right-hand side gives the advective outflux where we used that outflux is confined to a boundary layer of size Δρ(θ) = Nu −1 (θ), see Eq. (31), in the last equality. The integrated Eq. (32) thus transforms into the flux balance I in (θ) = I out (θ). For constant concentration boundary conditions (B), we obtain after differentiating with respect to θ Apart from the undetermined constant, this is exactly the differential equation governing the scaling function g(cos θ) in the similarity solution [36], which confirms Nu(θ) = g(cos θ). The differential equation can be solved to give the well-known exact result [36] Nu(θ) = 2constŪ This can be directly integrated to give a new approximative result for the angular dependence of the Nusselt number, The normalized local Nusselt numbers Nu(θ)/Nu are plotted as black lines in Fig. 4. The agreement for largē U is excellent for constant concentration boundary conditions (B) and approximate for constant flux boundary conditions (A), as expected. The flux balance approach also confirms the scaling Nu(θ) ∼Ū 1/3 , see Eqs. (29) and (30).

Anisotropic emission
Finally, we want to discuss the effect of an anisotropic emission boundary condition using the example of an anisotropic diffusive flux as characterized by a parameterβ > 0 in Eq. (18). In general, we expect higher Marangoni forces, because these forces are caused by anisotropies in the concentration profile around the half-sphere. If anisotropies are present without the need to create them by advection, this increases the Marangoni forces as can also be seen in the numerical results in Fig. 7. These numerical results also show that increasing the anisotropic emission parameterβ beyondβ ∼ 1 erases the maximum in the Marangoni forces in the crossover regionŪ ∼ 1 between diffusive and advective transport. In the diffusive limitŪ 1, the anisotropy leads to an additional zeroth-order termc (0) 1 (ρ) = −β/2ρ 2 in the concentration field, which results in as derived in Appendix C (see Eqs. (C.5) and (C.6)). These perturbative results are in excellent agreement with numerical FEM results as can be seen in Fig. 7. For sufficiently smallŪ , the zeroth-order term dominates. If this term dominates, Marangoni flow forces decrease the direct force because 3β/32 <β/2; this is similar to the results of Ref. [29], where also an explicitly asymmetric situation was considered.
In the advective limitŪ 1, a boundary layer of width Δρ ∼Ū −1/3 determines the physics. On the scale of the boundary layer thickness, the concentration drops from its surface valuec(ρ, θ) to zero. For constant flux boundary conditions (A), this led to a  (36) and (37). The solid lines forŪ > 1 are the scaling results (38) and (39). Forβ = 0, we recover the results from Fig. 6 concentration levelc(ρ, θ) ∼ Δρ ∼Ū −1/3 (see Eq. (29)) at the sphere. In the presence of an explicitly symmetry-breaking emission ∂ ρc1 (ρ = 1) =β, this contribution will also decay on the scale of the boundary layer Δρ, and we expect a corresponding contribu-tionβŪ −1/3 to the concentration level at the sphere, c(ρ = 1, θ) ∼ (const +β)Ū −1/3 . Because the direct Marangoni force scales asF M /Pe ∼c(ρ = 1, θ), this leads toF which is in good agreement with numerical results as shown in Fig. 7. The total Marangoni force scaling is dominated by the advective tail, which led toF M,tot ∼ PeΔθc(ρ = 1, θ); we find This result is also in good agreement with numerical results as shown in Fig. 7.

Diffusion-advection with strong Marangoni flow Pe Ū
For a strong Marangoni flow, Pe Ū , the linear response regimeŪ 1 becomes modified. We first have to address the dominant Marangoni flow problem (iib), which determines the Marangoni flow v M . For Pe Ū , this is the dominant contribution to 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 as can be seen in Fig. 2. Because this solution lacks axisymmetry, a complete and analytical solution is no longer possible.
Applying mass conservationJ ∼ 2πcv M ρl c = const and the Marangoni boundary condition to concentration profile and Marangoni flow field in a concentration boundary layer of widthl c ∼ (ρ/v M ) 1/2 below the fluid interface S Int , we find a scaling [25] for strong Marangoni flows. Here, we will further test this result in numerical FEM solutions, We see that the advective currentj ∼cv M ∼ Pe 1/3 ρ −7/3 becomes smaller than the corresponding diffusive currentj D ∼ −∂ ρ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, 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 withc ∝ ρ −1 sets in. We can also introduce the dimensionless Marangoni number for the radial Marangoni flow, which exactly compares advective Marangoni current and diffusive current by definition, and see that ρ M is determined by the condition that the regime Ma > 1 determines the size of the Marangoni vortex.
We can test the predictions (40) and (41) in numerical FEM solutions, see Fig. 8. One problem is that, for large Peclet numbers, the finite size of the numerical system becomes too small to accommodate the Marangoni vortex of size ρ M ∼ Pe properly. This results in deviations of the interfacial Marangoni flow field from Eq. (41). The numerical results for the interfacial concentration field show excellent agreement with (40).
So far, we considered the leading order of our problem by settingŪ ≈ 0; going one order further, we get the  (41) and (40) linear response for smallŪ with the ansatzc =c (0) + Uc (1) withc (0) (ρ) given by (40). In the total flow v + v M , the Marangoni flow (41)  where the radial componentū of the Stokes flow is considered. Expanding up to first order inŪ , we find a scalinḡ which will give rise to a Marangoni force scalinḡ Numerical FEM results show that both prefactors are of order unity (but hard to quantify because of finite size effects), see Fig. 9. This shows that Marangoni flows depress the total driving force in the linear response regime by a factor Pe −2/3 because it is harder to break the symmetry in the presence of the strong Marangoni flow advection. Numerical results in Fig. 9 also show that the total Marangoni force is somewhat larger than the direct Marangoni force,F M,tot >F M . 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 strong

Marangoni flows Pe
Ū as the numerical results in Fig. 9 show. The flow field v will still give rise to a concentration boundary layer of thickness Δρ ∼Ū −1/3 around the sphere. On the scale of the boundary layer, the Marangoni flows v M are not yet developed; they develop only further away at 1 ρ < ρ M ∼ Pe because of the no-slip boundary condition for the Marangoni flow in (iib). Therefore, the results forŪ 1 are essentially unaffected by a strong Marangoni flow for Pe Ū .

Diffusion-advection in the presence of evaporation
In the presence of evaporation, we have a convective (Robin) boundary condition (9b), which is governed by the dimensionless Biot number (11), instead of the Neumann condition (9a), which is recovered for vanishing Biot numberk = 0. In general, evaporation of surfactant depletes the interface of surfactant and, thus, decreases the Marangoni driving forces (both direct and flow forces). For volatile camphor, we find a Biot numberk = ak/D ≈ 550 using results from Ref. [18], whereas other surfactants such as PEG are non-volatile and have a very small Biot number [25]. The Biot number can also be interpreted as an extrapolation length scale. 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 Ref. [25], we developed a qualitative scaling theory based on the assumption that the total evaporation flux balances the total emission flux of surfactant in a stationary state. In the diffusive regimeŪ 1, this leads tō In the advection-dominated limitŪ 1, we find In both limits, Marangoni forces are reduced by evaporation, because it reduces the surfactant concentration.

Swimming condition, symmetry breaking and speed
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 velocitȳ U at low Reynolds numbers. In particular, we know the Marangoni forces as a function of the prescribed velocityŪ .

Swimming condition
The swimming condition (16) gives an additional force balance relation between Marangoni forces andŪ , which has to be satisfied in the swimming state and determines the selected swimming speedŪ =Ū swim as a function of Peclet number Pe and Biot numberk. In general, the swimming velocity increases with Pe and decreases withk. The force balance condition can be interpreted such that intersections of the linear Stokes friction relation −F D = 3πŪ and the total Marangoni forceF M,tot = F M,tot (Ū ) relation give the swimming speedŪ =Ū swim . The resulting swimming state can only be stable if the Marangoni force curveF M,tot (Ū ) intersects the straight Stokes friction line 3πŪ from above. Then, a speed fluctuation δŪ > 0 will give rise toF M,tot < −F D such that friction dominates, and the swimming speed is decreased again.
All curves (F M,tot /πPe)(Ū ) in Figs. 6, 7 and 9 start linearly ∝Ū in the diffusive regimeŪ 1 and then cross over to sublinear growth and finally decrease in the advective regimeŪ > 1. Therefore, all intersection points with the linear Stokes friction function will represent stable swimming states, also if an anisotropic emission is included. These results remain unchanged if evaporation is included.
In the decoupled limit Pe Ū , the total Marangoni force is always trivially linear in Pe. Then, we can  Fig. 7 in the absence of evaporation and in the decoupled limit. For β = 0, the blue vertical line of data points ends in the critical Peclet number Pec at zero swimming speed. In the presence of an anisotropic emissionβ > 0, the swimming bifurcation in the diffusive regime becomes avoided resulting in an initial linear relationŪswim ∝ Pe, crossing over tō Uswim ∝ Pe 3/5 in the advective regime. The solid lines for U < 1 are derived from the analytical perturbative results (36) and (37). The solid lines forŪ > 1 are derived from the scaling results (38) and (39) directly obtain the swimming condition in the form Using the Marangoni forces from Fig. 7 in the decoupled limit, which include an asymmetric emissionβ and reduce to Eq. (25) forβ = 0 and inverting this relation, we obtain the swimming relation in Fig. 10.

Swimming bifurcation
Forβ = 0, i.e., a symmetrically emitting swimmer, we see a sharp spontaneous symmetry breaking above a critical Peclet number Pe c in the swimming relation in Fig. 10 (blue vertical line of data points). From Eq. (25), we obtain the existence of a symmetry-broken swimming state for Pe > Pe c ∼ 8/ lnR → 0, which approaches zero for large system sizes. Therefore, the symmetry is essentially always spontaneously broken in a large swimming vessel. The swimming bifurcation in the force balance is governed by the leading-order linear terms ∝Ū (from the drag force and the linear response regime of the Marangoni forces) and the next order correction ∝Ū 3 in the Marangoni force. Therefore, we expect a supercritical pitchfork bifurcation analogously to a φ 4 -theory for a second-order phase transition. In the presence of the additional symmetry-breaking emission rate,β > 0, which contributes a constantŪ 0 -term to the force balance. This corresponds to an additional symmetry-breaking field in the φ 4 -theory and gives rise to an avoided bifurcation. This bifurcation scenario is clearly reflected in Fig. 10. Figure 10 and Eq. (25) were, however, derived for the decoupled limit Pe Ū . At the swimming bifurcation, we have Pe = Pe c Ū ≈ 0, such that the feedback of Marangoni flows onto the diffusion-advection problem has to be taken into account, and the decoupling approximation should not be used. Then, Eq. (44) describes the Marangoni forces in the linear response regime, which further reduces the critical Peclet number to Pe c ∼ 1/(lnR) 3 → 0. In the presence of evaporation withk 1, as appropriate for surfactants such as camphor, the total Marangoni force is further depressed according to Eq. (45) resulting in an increased Pe c =k 3 /(lnR) 3 → 0, which is, however, still approaching zero for large swimming vessel sizes R.
For strong Marangoni flows and in the presence of evaporation, we still have a linear response of the Marangoni forces for smallŪ (see Fig. 9 and Eq. (45)) with higher-order correction terms competing with a linearlyŪ -dependent drag force. Therefore, the above supercritical bifurcation scenario should persist.

Swimming relation
For Pe > Pe c , a spontaneously symmetry-broken swimming state withŪ swim > 0 exists for a symmetrically emitting swimmer withβ = 0. Because the Marangoni force Eq. (25) remains approximately linear up toŪ ∼ O(1), as can also be seen in Fig. 6, the swimming velocity rises steeply for Pe Pe c and quickly enters the asymptotics for the advection-dominated regimē U swim 1 as can be clearly seen in Fig. 10. In the advective regime, we find the swimming rela-tionsŪ 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 this regime. Evaporation is significant fork Pe 1/5 and reduces the swimming speed because it reduces the driving Marangoni forces.
For an anisotropically emitting swimmer withβ > 0, the bifurcation is avoided, and we find a linear swimming relation for small Pe. In the vicinity of the bifurcation, the force balance can be written asŪ swim = (Pe/Pe c )Ū swim + Peβ/32, which results in the linear swimming relation This describes the linear relationsŪ swim ∝ Pe for small Pe in the swimming relation in Fig. 10. In the advective regime, we still have a crossover to the above swimming relations (48), but with a slightly increased prefactor, i.e.,Ū swim ∼ (const +β)Pe 3/5 , (50a)  [26]. Here, the thermal diffusion constant replaces the surfactant diffusion constant and is by a factor O(10 3 ) larger. Moreover, radii a ∼ 3 μm could be reached. At the same time, swimming velocities are still in the range above 10 3 −10 5 μm/s. These parameters correspond to dimensionless velocitiesŪ swim ∼ 2 × 10 −2 − 2, which is mostly in the diffusive regimeŪ 1 and at low Reynolds numbers Re ∼ 6 × 10 −6 . These swimmers were asymmetrically heated with a temperature difference ΔT across the swimmer corresponding to a constant concentration asymmetryc S,1 ∝ ΔT . We therefore expect to be in a constant concentration situation, which is analogous to the linear regime in the constant flux swimming relation in Fig. 10. This is in accordance with the theoretical results of Würger [30], because advection plays no role in this regime and agrees with the experimental observations in Ref. [26].
Our 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. Conceptually, there is no difference for a thermal Marangoni surfer as realized in Ref. [26]. In previous theoretical approaches to surfactant [29] or thermal [30] Marangoni boats, advection has been neglected. For surfactant driven Marangoni boats, this is typically a bad approximation as estimates in Ref. [25] show; for thermal Marangoni boats, this is typically justified as our above estimates show.
The three coupled problems of surfactant adsorption, low Reynolds number fluid flow and diffusion-advection of surfactant 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-water-swimmer contact line and the Marangoni flow force. We employ the reciprocal the-orem, which we could reinterpret in terms of energy transduction, to calculate the Marangoni forces.
Non-dimensionalization reveals that two dimensionless control parameters exist, the Peclet number (10), which is the dimensionless emission rate of surfactant, and the Biot number (11), which is the dimensionless evaporation rate. Evaporation is practically absent for PEG (Biot numberk 1), but strong for other frequently studied soap boat swimmers such as camphor boats (Biot numbersk ≈ 550 [18]). In Ref. [25], it is shown that evaporation is relevant to quantitatively understand the large differences in the swimming rela-tionŪ swim =Ū swim (Pe) between PEG-alginate swimmers and camphor boats from Ref. [20], but these Marangoni boats operate at moderate Reynolds numbers. For thermal Marangoni surfers [26], evaporation corresponds to a convective boundary condition for heat transfer from the water surface to the air; the corresponding convection coefficient will depend on the nature of the air flow that is applied to transfer heat, which is difficult to quantify. It also depends on the temperature difference to the surrounding air that can be established. Because the thermal Marangoni surfers from Ref. [26] mostly operate in the diffusive regimē U 1, we expect convection to reduce the Marangoni force according to Eq. (45) if the corresponding Biot numberk is sufficiently high.
Moreover, the dimensionless swimmer velocityŪ plays an important role as it controls the transition from a diffusive regimeŪ 1 to an advective regimē U 1. Non-dimensionalization of the coupled equations also shows a decoupling of the Marangoni flow problem for weak Marangoni flows Pe Ū . Then, the concentration field around the interfacial Marangoni swimmer with velocity U is essentially equivalent to the concentration field around a mass emitting sphere moving with velocity U through a bulk viscous fluid, which is a classical diffusion-advection problem. We developed solutions for this diffusion-advection problem for two types of boundary conditions which seem most important for applications: constant flux boundary conditions (A) for diffusive emission of surfactant from the swimmer and constant concentration boundary conditions (B) if the surfactant dissolves from the surface or is produced by a chemical reaction on the surface. We could obtain novel results for constant flux boundary conditions, which are unusual in the diffusion-advection literature. In particular, we could obtain qualitative results for the local Nusselt number by a novel flux balance argument. All theoretical results are supported by numerical FEM simulations.
Apart from extensive results for the decoupled limit Pe Ū , we also addressed strong Marangoni flow in the limit Pe Ū and evaporation on the basis of scaling arguments and numerical FEM simulations. This allowed us to obtain the Marangoni forces as a function of a prescribed swimmer speedŪ for all relevant situations, also including a possible anisotropic emission. For all cases, our theoretical results agree well with the numerical FEM calculations. Knowledge of the Marangoni forces is the basis to discuss the swimming bifurcation and swimming speed as a function of the Peclet number as main control parameter via the force balance condition.
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. The swimming bifurcation is a supercritical pitchfork bifurcation analogous to a second order symmetry-breaking phase transition, and the presence of an explicitly symmetry-breaking emission gives rise to an avoided bifurcation. Spontaneous symmetry breaking resulting in propulsion is possible by establishing an asymmetric surfactant concentration profile that is maintained by advection. The symmetry breaking mechanism is similar to what has been proposed for autophoretic swimmers [9,10] and liquid Marangoni swimmers [8] before.
In Eq. (48), 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). In Eq. (50), the result is extended in the presence of an explicitly symmetrybreaking emission. Then, a linear regime emerges in the diffusive limitŪ 1, which is caused by the avoided bifurcation. This regime is observed for the thermal Marangoni surfers in Ref. [26].

A Energy transduction and reciprocal theorem
In this Appendix, we discuss the reciprocal theorem in terms of energy transduction, in order to see how power input from Marangoni stresses via the flow fields is transduced to the sphere for propulsion. This will also provide an alternative derivation of the results obtained by Masoud and Stone [37] via the reciprocal theorem.
For this, we switch to the laboratory frame in the following. We first consider the dissipation rate Φ ≡ 2μ V dV eijeij of a solution of the Stokes equation in an arbitrary volume V . Here, eij = 1 2 (∂ivj + ∂jvi) is the strain tensor, and σij = −pδij + 2μeij is the stress tensor. The kinetic energy K = 1 2 ρ V dV v 2 of the fluid changes according to dK/dt = ∂V daivjσij − Φ, i.e., by the external power input P ∂V across the surface of the volume and by dissipation. In a stationary state, dK/dt = 0 and dissipation and power input are equal, (da is the outward normal to volume V ). Applying this equality to the Stokes flow field v and the liquid volume with the boundary ∂V = S + SInt and using σij = 0 at the liquid-air interface SInt, we find with the Stokes drag force FD ≡ − S daiσiz from Eq. This means that the power input by Marangoni stresses σM,ij on the interface SInt via the Marangoni flows vM,j (right hand side) is dissipated entirely within the fluid without transmitting any mechanical power onto the half-sphere because vM = 0 on S. Now, we consider the mutual dissipation rate for two solutions v (1) (r) and v (2) (r) to the Stokes equation. The mutual dissipation can be shown to be given by the mutual power input in a stationary state, The symmetry Φ (12) = Φ (21) leads directly to the reciprocal theorem ∂V daiv (2) j σ (1)

C Diffusion-advection equation, perturbation theory forŪ 1
For small fluid velocities,Ū 1, we can expand about the isotropic undisturbed diffusion solution atŪ = 0, which is given byc for each Legendre coefficient. It turns out that this will work only in the "inner region" ρ < 1/Ū of a solution, because in the "outer region" ρ 1/Ū , the convection term can no longer be treated perturbatively. In Ref. [34], a systematic expansion in inner and is used in the "inner region" ρ < 1/Ū . In the "outer region" ρ 1/Ū , the convection term can no longer be treated perturbatively, regardless how smallŪ is [34]. Here, we rescale σ ≡ ρŪ andC(σ, η) =c(σ/Ū, η) to obtain an equation in all orders inŪ and as a function of η or in all Legendre coefficients.