Hydrodynamic efficiency limit on a Marangoni surfer

Abstract

A Marangoni surfer is an object embedded in a gas–liquid interface, propelled by gradients in surface tension. We derive an analytical theorem for the lower bound on the viscous dissipation by a Marangoni surfer in the limit of small Reynolds and capillary numbers. The minimum dissipation can be expressed with the reciprocal difference between drag coefficients of two passive bodies of the same shape as the Marangoni surfer, one in a force-free interface and the other in an interface with surface incompressibility. The distribution of surface tension that gives the optimal propulsion is given by the surface tension of the solution for the incompressible surface and the flow is a superposition of both solutions. For a surfer taking the form of a thin circular disk, the minimum dissipation is $16\mu a V^2$, giving a Lighthill efficiency of $1/3$. This places the Marangoni surfers among the hydrodynamically most efficient microswimmers.

1. Introduction

The Marangoni effect describes the motion of liquids due to surface tension gradients, caused, for example, by an uneven distribution of surfactants, and can be harnessed to drive the motion of a microswimmer. The Marangoni propulsion can thereby take place on the swimmer's surface, like in the case of chemically active emulsion droplets (Thutupalli, Seemann & Herminghaus 2011; Herminghaus et al. 2014; Izri et al. 2014; Maass et al. 2016; Schmitt & Stark 2016; Testa et al. 2021; Hokmabad et al. 2022; Michelin 2023). The more common design, however, involves particles embedded in a flat gas–liquid interface that emit surfactants and thereby ‘manipulate’ the surface tension in the interface around them. A classical example are camphor boats, already studied by Tomlinson (1862) and Rayleigh (1890). These surfers create a surface tension gradient by asymmetric emission of a surfactant, but it is also possible for a symmetric swimmer to achieve directed propulsion through spontaneous symmetry breaking (Boniface et al. 2021). Marangoni swimmers have also been fabricated using accessible pen-drawn patterns and powered by fuel in the form of ink, enabling controlled movement and navigation (Song et al. 2023). Finally, some animals like the Microvelia water striders use the Marangoni effect to slide along the surface of water (Bush & Hu 2006).

For certain distributions of surface tension, the propulsion of a thin disk-shaped Marangoni surfer has been exactly solved (Lauga & Davis 2012; Elfring, Leal & Squires 2016; Crowdy 2021). Using the Lorentz reciprocal theorem (Masoud & Stone 2014), the propulsion velocity of a surfer with any shape can in principle be determined provided the solution of the passive particle pulled by an external force is known. Solutions including inertia at finite Reynolds numbers have also been determined (Ender & Kierfeld 2021). Moreover, two-dimensional squirmer models share certain similarities with the Marangoni effect (Matas-Navarro et al. 2014).

The energetic efficiency of a microswimmer has been defined by Lighthill (1952) as the power needed to pull the swimmer through a fluid by an external force, divided by the actual power dissipated when the swimmer actively moves at the same velocity. Although Lighthill's efficiency can theoretically exceed 1 (Leshansky et al. 2007), swimming microorganisms typically achieve values of around $1\,\%$ (Osterman & Vilfan 2011). The theoretical efficiency limit of the three-sphere model swimmer is also of a similar order of magnitude (Nasouri, Vilfan & Golestanian 2019). The efficiency of autophoretic colloidal microswimmers is many orders of magnitude lower than that and is limited by the thickness of the boundary layer (Sabass & Seifert 2010). We have recently derived a minimum dissipation theorem that provides a lower bound on the dissipation by a microswimmer moving through bulk fluid, first for external dissipation alone (Nasouri, Vilfan & Golestanian 2021) and later for the combination of internal and external dissipation (Daddi-Moussa-Ider, Golestanian & Vilfan 2023a). The theorems share the common structure that expresses the dissipation with the reciprocal difference between the drag coefficients of two passive bodies. For example, for the swimmer without internal dissipation, these would be two bodies of the same shape as that of the swimmer, one with the no-slip and one with the perfect-slip boundary (Nasouri et al. 2021). These solutions lead to the question as to whether a similar theorem can be derived for the hydrodynamic efficiency of Marangoni surfers. They differ from the previously solved dissipation problems in two main aspects: first, the propulsive force acts on an infinite plane outside the surfer; and second, the force cannot be optimised freely, but needs to have the form of a surface tension gradient. In this paper, we derive such a theorem – first in a general form and then for a thin surfer with the shape of a circular disk.

2. Minimum dissipation theorem for Marangoni surfers

2.1. Problem formulation

In this section we derive a lower bound for the energy dissipation by a Marangoni surfer of arbitrary shape moving along the surface of a fluid at zero Reynolds number (figure 1a). The surfer comprises a body partially immersed in a gas–liquid interface and suspended by surface tension. The surfer can modify the surface tension in its surroundings, for example by emitting surfactants. For the purpose of this study, we disregard the dynamics of surfactants and treat the surface tension as a given. We assume that the interface remains planar, i.e. it is distorted neither by the induced flows (this corresponds to the low-capillary-number limit) nor by the contact angle at the surfer boundary. The surfer is force-free, driven solely by the surface tension at the contact line and by the flows induced by the Marangoni effect. The submerged part of the surfer forms a no-slip boundary with the fluid. The interface extends across the $x$–$y$ plane, with the $z$ direction oriented perpendicular to this plane.

Figure 1. (a) The optimal Marangoni surfer. The surface colour indicates increased (red) or reduced (green) surface tension leading to the Marangoni effect. (b) A passive object, pulled along the incompressible surface of a fluid with velocity $\boldsymbol {V}_{SI}$. The colours denote the surface tension that builds up in order to maintain the incompressibility condition. (c) A passive object, pulled with velocity $\boldsymbol {V}_{FS}$ in a fluid with a free surface (with a uniform surface tension).

The fluid motion in bulk ($z<0$) is governed by the Stokes equation along with the incompressibility condition:

(2.1a)$$\begin{gather} -\boldsymbol{\nabla}p + \mu \nabla^2 \boldsymbol{v} = \boldsymbol{0}, \end{gather}$$

(2.1b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0, \end{gather}$$

where $\boldsymbol {v}$ and $p$ denote the velocity and pressure fields in the fluid medium and $\mu$ the shear viscosity.

At the submerged part of the swimmer surface $\mathcal {S}$, not necessarily at $z=0$ if the swimmer has a non-flat shape, the fluid is subject to the no-slip boundary $\boldsymbol {v}=\boldsymbol {V}_{A}$, where $\boldsymbol {V}_{A}$ is the translational velocity of the swimmer ($\boldsymbol {V}_{A}\perp \hat {\boldsymbol {e}}_z$). The boundary condition at the gas–liquid interface located at $z=0$ implies $\hat {\boldsymbol {e}}_z \boldsymbol {\cdot } \boldsymbol {v}=0$.

The horizontal force balance on the swimmer states $(\boldsymbol {I}-\hat {\boldsymbol {e}}_z \hat {\boldsymbol {e}}_z) \boldsymbol {\cdot } \int _{\mathcal {S}} {\rm d} S \, \boldsymbol {\sigma } \boldsymbol {\cdot } \hat {\boldsymbol {n}} + \int _\ell {\rm d} s \, \hat {\boldsymbol {n}} \gamma =0$. The first term represents the tractions exerted on the swimmer by the fluid and the second term the effect of the surface tension. Here $\hat {\boldsymbol {n}}$ denotes the surface normal pointing into the fluid and $\ell$ the contact line between the swimmer and the gas–liquid interface (figure 2). In the integral over $\ell$, $\hat {\boldsymbol {n}}$ denotes the in-plane normal to the contact line. The stress tensor is determined as $\boldsymbol {\sigma }=-p \boldsymbol {I}+2\mu \boldsymbol{\mathsf{E}}$, with the strain-rate tensor $\boldsymbol{\mathsf{E}}=(\boldsymbol {\nabla } \boldsymbol {v} + (\boldsymbol {\nabla } \boldsymbol {v})^\top )/2$. The surface tension, reduced by its value in the surface that is not affected by the presence of the surfer, is denoted by $\gamma$. At the gas–liquid interface, the force balance states $(\boldsymbol {I}-\hat {\boldsymbol {e}}_z \hat {\boldsymbol {e}}_z) \boldsymbol {\cdot } \boldsymbol {\sigma }\boldsymbol {\cdot } \hat {\boldsymbol {e}}_z +\boldsymbol {\nabla }_{s} \gamma =0$. Here, $\boldsymbol {\nabla }_{s} = \hat {\boldsymbol {e}}_x \partial _x + \hat {\boldsymbol {e}}_y \partial _y$ represents the gradient within the horizontal plane.

Figure 2. Side view of the surfer. The fluid volume $\mathcal {V}$ is shown in blue, the gas–liquid interface $\mathcal {I}$ as a solid blue line, the submerged surface of the swimmer $\mathcal {S}$ in grey and the contact line $\ell$ in yellow.

The total energy dissipation can be written as either the volume integral of the local dissipation rate or the rate of work exerted by the surface tension, integrated over its surface (Happel & Brenner 1983):

(2.2)\begin{equation} P=\int_{\mathcal{V}} {\rm d} V \, 2\mu \boldsymbol{\mathsf{E}}\boldsymbol{:} \boldsymbol{\mathsf{E}} = \int_{\mathcal{I}} {\rm d} S \, \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\sigma}\boldsymbol{\cdot} \hat{\boldsymbol{e}}_z + \int_{\mathcal{S}} {\rm d} S \, \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\sigma}\boldsymbol{\cdot} (-\hat{\boldsymbol{n}}). \end{equation}

By taking into account the no-slip boundary at $\mathcal {S}$ and the force balance on the swimmer and at the interface, the dissipation can be expressed as

(2.3)\begin{equation} P= \int_{\mathcal{I}} {\rm d} S \, \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}_{s}\gamma + \int_\ell {\rm d} s \, \gamma \hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{V}_{A}. \end{equation}

Finally, we apply the two-dimensional divergence theorem $\int _{\mathcal {I}} {\rm d} S \boldsymbol {\nabla }_{s}\boldsymbol {\cdot }(\gamma \boldsymbol {v})=-\int _\ell {\rm d} s \, \hat {\boldsymbol {n}}\boldsymbol {\cdot } (\gamma \boldsymbol {v})$ to express the dissipation as

(2.4)\begin{equation} P=- \int_{\mathcal{I}} {\rm d} S \, \gamma \boldsymbol{\nabla}_{s}\boldsymbol{\cdot} \boldsymbol{v}, \end{equation}

which represents the rate of change of the total surface energy.

2.2. Derivation of the theorem

In the following we derive a theorem for a lower bound on the dissipation by a Marangoni surfer moving with velocity $\boldsymbol {V}_{A}$. We follow the approach that led to a minimum dissipation theorem for surface-propelled swimmers with external dissipation (Nasouri et al. 2021) and later with combined external and internal dissipation (Daddi-Moussa-Ider et al. 2023a). Using this theorem, the flow field of an optimal nearly spherical swimmer with external dissipation has been obtained using a perturbative analytical method (Daddi-Moussa-Ider et al. 2021a).

The solution consists of two major steps. The first is to find a passive minimum dissipation theorem for flows that satisfy the velocity boundary condition on an object. This passive theorem can be seen as a generalisation of the Helmholtz minimum dissipation theorem that states that among all incompressible flows that satisfy the same fixed-velocity boundary condition, the Stokes flow has the smallest dissipation (Guazzelli & Morris 2009). For example, for surface-driven swimmers with external dissipation, the passive problem consists of a perfect-slip body. In the second step, another passive problem is needed that also satisfies the boundary condition, but is orthogonal to the active problem in the sense that their superposition flow has a dissipation that is the sum of the dissipations of the two flows, each on its own. The minimum dissipation theorem is obtained by applying the inequality from the first condition to this superposition flow.

As a passive minimum dissipation theorem, we use a variant of the Helmholtz minimum dissipation theorem that is derived in Appendix A. The theorem states that among flows that satisfy the no-slip boundary condition on the swimmer moving with a given velocity and zero normal velocity at the fluid–gas interface, the flow without tangential tractions at the interface has the smallest dissipation (figure 1c). Because of the stress-free interface (free surface), we label this solution $\boldsymbol {v}_{FS}$. If the passive body with the shape of the Marangoni surfer is moved with a velocity $\boldsymbol {V}_{FS}$, the drag force acting on it is $\boldsymbol {F}_{FS}= -\boldsymbol {R}_{FS}\boldsymbol {\cdot } \boldsymbol {V}_{FS}$ with the free-surface drag coefficient $\boldsymbol {R}_{FS}$. The total dissipation of any flow, possibly with additional surface forces, then satisfies the inequality

(2.5)\begin{equation} P\ge \boldsymbol{V}_{FS} \boldsymbol{\cdot} \boldsymbol{R}_{FS}\boldsymbol{\cdot} \boldsymbol{V}_{FS}. \end{equation}

The equality is satisfied exactly when there are no additional horizontal forces beyond those pulling the passive body. In the following, we assume that the symmetries of the surfer are such that translational and rotational drag are decoupled. Then a force on the body leads to purely translational motion and a torque to purely rotational motion. The translational drag is then described by a symmetric $2\times 2$ matrix $\boldsymbol {R}_{FS}$.

The second step towards the active minimum dissipation theorem requires finding a problem that is orthogonal to the active problem (Marangoni surfer) in the sense that the dissipation in the superposition of both flows is additive. We will show that this condition is fulfilled for a flow that satisfies the no-slip boundary at the swimmer body and is surrounded by a gas–liquid interface with an incompressible surface, $\boldsymbol {\nabla }_{s} \boldsymbol {\cdot } \boldsymbol {v}_{SI}=0$ at $z=0$ (figure 1b). The incompressibility is achieved by a build-up of a passive surface tension $\gamma _{SI}$, which we also define relative to the unperturbed surface. The incompressible surface occurs in the case of insoluble surfactants in the limit of an infinite Marangoni modulus (Elfring et al. 2016; Manikantan & Squires 2020). The same model also represents a limiting case of the Boussinesq–Scriven model (Scriven 1960; Manikantan & Squires 2017) with an incompressible ‘membrane’, however without additional surface viscosity, i.e. in the limit of vanishing Boussinesq number (Stone & Masoud 2015).

The superposition of the Marangoni surfer moving with velocity $\boldsymbol {V}_{A}$ and the passive object with incompressible surface moving with $\boldsymbol {V}_{SI}$ has a dissipation rate that can be expressed in analogy to (2.3). It corresponds to the sum of the work done by the superposition of external forces (which is just $-\boldsymbol {F}_{SI}$, because $\boldsymbol {F}_{A}=\boldsymbol {0}$), the superposition of both surface tensions at the contact line and the superposition of surface tension gradients on the fluid:

(2.6)\begin{align} P_{A+SI} &= \int_{\mathcal{I}} {\rm d} S \, (\boldsymbol{v}_{A}+\boldsymbol{v}_{SI})\boldsymbol{\cdot} \boldsymbol{\nabla}_{s}(\gamma_{A} +\gamma_{SI}) + \int_\ell {\rm d} s \, (\gamma_{A}+\gamma_{SI}) \hat{\boldsymbol{n}}\boldsymbol{\cdot} (\boldsymbol{V}_{A}+\boldsymbol{V}_{SI}) \nonumber\\ &\quad -\boldsymbol{F}_{SI} \boldsymbol{\cdot} (\boldsymbol{V}_{A} +\boldsymbol{V}_{SI}). \end{align}

Two terms immediately cancel out because of the divergence theorem:

(2.7)\begin{equation} \int_{\mathcal{I}} {\rm d} S \, \boldsymbol{v}_{SI} \boldsymbol{\cdot} \boldsymbol{\nabla}_{s}\gamma_{SI} + \int_\ell {\rm d} s \, \gamma_{SI} \hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{V}_{SI} =- \int_{\mathcal{I}} {\rm d} S \, \gamma_{SI} \boldsymbol{\nabla}_{s}\boldsymbol{\cdot} \boldsymbol{v}_{SI} =0. \end{equation}

We now apply the Lorentz reciprocal theorem (Masoud & Stone 2019) to the two problems:

(2.8)\begin{align} \int_{\mathcal{I}} {\rm d} S \, \boldsymbol{v}_{SI} \boldsymbol{\cdot} \boldsymbol{\sigma}_{A} \boldsymbol{\cdot} \hat{\boldsymbol{e}}_z+ \int_{\mathcal{S}} {\rm d} S\, \boldsymbol{v}_{SI} \boldsymbol{\cdot} \boldsymbol{\sigma}_{A}(-\hat{\boldsymbol{n}})= \int_{\mathcal{I}} {\rm d} S \, \boldsymbol{v}_{A} \boldsymbol{\cdot} \boldsymbol{\sigma}_{SI} \boldsymbol{\cdot} \hat{\boldsymbol{e}}_z+ \int_{\mathcal{S}} {\rm d} S\, \boldsymbol{v}_{A} \boldsymbol{\cdot} \boldsymbol{\sigma}_{SI}(-\hat{\boldsymbol{n}}). \end{align}

By taking into account the force balance at the surface and at both bodies, we obtain

(2.9)\begin{equation} \int_{\mathcal{I}} {\rm d} S \, \boldsymbol{v}_{SI} \boldsymbol{\cdot} \boldsymbol{\nabla}_{s}\gamma_{A}+ \int_\ell {\rm d} s \, \gamma_{A} \hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{V}_{SI} = \int_{\mathcal{I}} {\rm d} S \, \boldsymbol{v}_{A} \boldsymbol{\cdot} \boldsymbol{\nabla}_{s}\gamma_{SI}+ \int_\ell {\rm d} s \, \gamma_{SI} \hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{V}_{A} - \boldsymbol{F}_{SI} \boldsymbol{\cdot} \boldsymbol{V}_{A}. \end{equation}

At the same time, the divergence theorem applied to $\gamma _{A} \boldsymbol {v}_{SI}$ leads to

(2.10)\begin{equation} \int_{\mathcal{I}} {\rm d} S \, (\boldsymbol{v}_{SI} \boldsymbol{\cdot} \boldsymbol{\nabla}_{s}\gamma_{A} + \gamma_{A} \boldsymbol{\nabla}_{s} \boldsymbol{\cdot} \boldsymbol{v}_{SI})+ \int_\ell {\rm d} s \, \gamma_{A} \hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{V}_{SI} = 0. \end{equation}

Because we have chosen $\boldsymbol {\nabla }_{s} \boldsymbol {\cdot } \boldsymbol {v}_{SI}=0$, (2.10) implies that both sides of (2.9) are zero. The total contribution of all mixed terms in (2.6) therefore vanishes and the dissipation of the superposition flow

(2.11)\begin{equation} P_{A+SI}=\int_{\mathcal{I}} {\rm d} S \, \boldsymbol{v}_{A}\boldsymbol{\cdot} \boldsymbol{\nabla}_{s}\gamma_{A} + \int_\ell {\rm d} s \, \gamma_{A} \hat{\boldsymbol{n}}\boldsymbol{\cdot} \boldsymbol{V}_{A}-\boldsymbol{F}_{SI} \boldsymbol{\cdot} \boldsymbol{V}_{SI} = P_{A}+P_{SI} \end{equation}

can be written as the sum of the dissipations in each of the problems. This additivity is the second condition for deriving a minimum dissipation theorem for the active swimmer.

The derivation of the minimum dissipation theorem follows in the same way as for surface-driven swimmers (Nasouri et al. 2021; Daddi-Moussa-Ider et al. 2023a). We apply the inequality (2.5) to the superposition:

(2.12)\begin{equation} P_{A+SI} \ge (\boldsymbol{V}_{A} +\boldsymbol{V}_{SI}) \boldsymbol{\cdot} \boldsymbol{R}_{FS} \boldsymbol{\cdot} (\boldsymbol{V}_{A} +\boldsymbol{V}_{SI}). \end{equation}

This inequality holds for any values of the velocities, but we know that the equality is fulfilled if and only if the superposition represents exactly the flow around the passive body with free surface (FS problem). This condition gives the following set of equations for the superposition: $\boldsymbol {V}_{A}+\boldsymbol {V}_{SI}=\boldsymbol {V}_{FS}$ and $\boldsymbol {F}_{SI}= \boldsymbol {F}_{FS}$. In the absence of translational–rotational coupling, the torque balance is satisfied automatically. At the same time, the drag forces in the FS and SI problem are related to their velocities through through the drag coefficients $\boldsymbol {R}_{FS}$ and $\boldsymbol {R}_{SI}$ as $\boldsymbol {F}_{FS}=- \boldsymbol {R}_{FS}\boldsymbol {\cdot }\boldsymbol {V}_{FS}$ and $\boldsymbol {F}_{SI}=-\boldsymbol {R}_{SI}\boldsymbol {\cdot }\boldsymbol {V}_{SI}$, which closes the equation system. Its solution is

(2.13a)$$\begin{gather} \boldsymbol{V}_{FS}=(\boldsymbol{I} -\boldsymbol{R}_{SI}^{-1} \boldsymbol{\cdot}\boldsymbol{R}_{FS})^{-1}\boldsymbol{\cdot}\boldsymbol{V}_{A}, \end{gather}$$

(2.13b)$$\begin{gather}\boldsymbol{V}_{SI}=(\boldsymbol{R}_{FS}^{-1}\boldsymbol{\cdot} \boldsymbol{R}_{SI}-\boldsymbol{I})^{-1}\boldsymbol{\cdot}\boldsymbol{V}_{A}. \end{gather}$$

By inserting these velocities into the inequality (2.12) we obtain the minimum dissipation theorem:

(2.14)\begin{equation} P_{A} \geq \boldsymbol{V}_{A}\boldsymbol{\cdot}(\boldsymbol{R}_{FS}^{-1} -\boldsymbol{R}_{SI}^{-1})^{-1} \boldsymbol{\cdot}\boldsymbol{V}_{A}, \end{equation}

which allows us to express a lower bound on the hydrodynamic dissipation by a Marangoni surfer with two drag coefficients for the horizontal motion of a body with the same shape in a fluid without surface tension gradients $\boldsymbol {R}_{FS}$ and with surface incompressibility $\boldsymbol {R}_{SI}$. Although we wrote the expressions for translational motion alone, it is straightforward to expand it to rotational motion, too. In this case, the generalised velocities become three-component vectors with two translational and one rotational component and $\boldsymbol {R}_{FS}$ and $\boldsymbol {R}_{SI}$ become the corresponding grand resistance matrices.

From the superposition, we also obtain the flow field of the optimal Marangoni surfer as $\boldsymbol {v}_{A}=\boldsymbol {v}_{FS}-\boldsymbol {v}_{SI}$ and its relative surface tension $\gamma _{A}=-\gamma _{SI}$. Here both passive solutions are evaluated for bodies moving with the velocities as determined by (2.13). Along with the bound on dissipation, the minimum dissipation theorem also provides us with the solution for the surface tension of the active Marangoni surfer that will minimise the hydrodynamic dissipation.

2.3. Minimum dissipation for a circular disk-shaped Marangoni surfer

In the following, we employ the minimum dissipation theorem (2.14) to a thin Marangoni surfer with the shape of a circular disk. For this geometry, both drag coefficients are known from the literature. For an incompressible surface, it can be found as a limiting case in the solutions of Hughes, Pailthorpe & White (1981) and Stone & Ajdari (1998):

(2.15)\begin{equation} R_{SI} = 8\mu a. \end{equation}

The solution with a free planar surface is mathematically equivalent to the edgewise motion of a thin disk in bulk fluid (Happel & Brenner 1983), with half the drag coefficient:

(2.16)\begin{equation} R_{FS} = \tfrac{16}{3} \, \mu a. \end{equation}

With these drag coefficients, the minimum dissipation theorem gives us a lower bound on the hydrodynamic dissipation by a disk-shaped Marangoni surfer as

(2.17)\begin{equation} P_{A} \ge 16\mu a V_{A}^2. \end{equation}

The Lighthill efficiency (Lighthill 1952), defined as the ratio between the dissipation by a passive swimmer steadily towed by an external force and the active swimmer, has the upper limit

(2.18)\begin{equation} \eta_{L}=\frac{R_{FS}V_{A}^2 }{P_{A}} \le \frac 1 3 . \end{equation}

Unlike in surface-slip-driven microswimmers (Leshansky et al. 2007), the Lighthill efficiency is always bounded by 1. A hydrodynamic efficiency of $1/3$ places the Marangoni surfers to the top of realistic swimmer designs. Although an ideal spherical surface-driven swimmer has an efficiency limit of $\eta _{L}\le 1/2$ (Michelin & Lauga 2010; Guo et al. 2021), which is even higher for elongated swimmers (Nasouri et al. 2021), such high efficiency requires a frictionless mechanism of slip generation, which is difficult to realise.

The full solution of the optimal disk-shaped surfer that includes surface tension, fluid velocity and forces can be obtained as a superposition of these fields in the two passive problems. For this purpose, we derive both solutions using a uniform notation in §§ 3 and 4.

2.4. Approximate solution for a surfer of finite depth

If the surfer has a finite depth, it is possible to derive a lower bound on the dissipation in a perturbative way. Specifically, we treat a surfer with the shape of an oblate spheroid with major semi-axis $a$ and minor semi-axis $\varepsilon a$. The spheroid is half-submerged, such that its surface $\mathcal {S}$ is given by the equation $z=-\varepsilon \sqrt {a^2-\rho ^2}$ in cylindrical coordinates.

Perturbative solutions for both drag coefficients, $R_{FS}$ and $R_{SI}$, were obtained by Stone & Masoud (2015). In the case of a force-free surface, the drag coefficient of the half-submerged spheroid is $1/2$ the coefficient for edgewise translation in bulk fluid, which is known exactly (Happel & Brenner 1983). For thin spheroids, its expansion reads

(2.19)\begin{equation} R_{FS} = \frac{16}{3} \mu a \left(1 + \frac{8}{3{\rm \pi}}\varepsilon \right) + O (\varepsilon^2). \end{equation}

The solution in the presence of an incompressible interface was formulated using the Lorentz reciprocal theorem, with the solution for a thin disk serving as the auxiliary problem (Stone & Masoud 2015). An analytical expression can be obtained by integrating (3.6) from Stone & Masoud (2015) in the limit $Bq=0$:

(2.20)\begin{equation} R_{SI} = 8\mu a \left(1 + \frac{4(1-\ln 2)}{\rm \pi}\varepsilon \right) + O (\varepsilon^2). \end{equation}

The minimum dissipation follows from (2.14) as

(2.21)\begin{equation} P_{A} \geq (R_{FS}^{-1}-R_{SI}^{-1})^{-1} = 16\mu a V_{A}^2 \left(1 + \frac{8 \ln 2}{\rm \pi}\varepsilon\right)+ O(\varepsilon^2) \end{equation}

and the upper bound on Lighthill efficiency is

(2.22)\begin{equation} \eta_{L} \leq 1-\frac{R_{FS}}{R_{SI}}= \frac{1}{3} - \frac{8}{9{\rm \pi}} (3\ln 2-1) \varepsilon + O (\varepsilon^2). \end{equation}

The finite depth of the swimmer not only increases the dissipation, but also reduces its Lighthill efficiency. Similar asymptotic solutions are also possible for other swimmer shapes, for example for a partially submerged sphere, with the submerged surface $\mathcal {S}$ taking the shape of a spherical cap (Stone & Masoud 2015).

3. Thin circular disk translating in a fluid with surface incompressibility

In this section, we derive the hydrodynamic field resulting from the translational movement of a passive circular disk positioned at an interface with surface incompressibility. Theoretical studies of the flow past a disk embedded in an incompressible viscous surface, representing a model of a lipid bilayer, date back to Saffman & Delbrück (1975). Saffman (1976) examined the limit ${Bq} \gg 1$ where the Boussinesq number is defined as ${Bq} = \mu _{s} / (a\mu )$ and $\mu _{s}$ represents the surface viscosity. In this limit, the membrane viscosity dominates over the viscosity of the surrounding fluid. Hughes et al. (1981) solved the model for arbitrary viscosities, providing analytical expressions for the viscous flow field and resistance coefficients through the use of dual integral equations. The problem with a finite subphase depth was subsequently solved by Stone & Ajdari (1998), partly through a numerical approach. For an infinite depth and vanishing surface viscosity (${Bq} \to 0$), both solutions give a translational drag coefficient $8\mu a$, which is 50 % larger than the drag coefficient $(16/3)\mu a$ in a fluid with free surface. More recently, Yariv et al. (2023) revisited the translational motion of a disk embedded in a nearly inviscid Langmuir film, i.e. in the limit ${Bq} \ll 1$.

Because the solutions in Hughes et al. (1981) and Stone & Ajdari (1998) are derived for arbitrary $Bq$ and have a more complex form, we rederive the solution without surface viscosity in the following in order to obtain expressions for the velocity and force fields that are more suitable for practical usage. In the following we drop the ‘SI’ subscript as all quantities involved belong to the problem with surface incompressibility.

3.1. Green's functions in Fourier space

We begin by solving the flow equations (2.1) for a given distribution of surface forces. The impermeability condition requires zero normal velocity at the interface, i.e. $v_z = 0$ at $z=0$. The force balance at the interface in the in-plane direction states

(3.1)\begin{equation} \mu \partial_z \boldsymbol{v}_\parallel|_{z=0}= \begin{cases} \boldsymbol{f}_\parallel & \text{for}\ \rho< a \\ \boldsymbol{\nabla}_{s} \gamma & \text{for}\ \rho>a, \end{cases} \end{equation}

wherein $\boldsymbol {f}_\parallel$ is an unknown in-plane surface force density acting on the fluid from the surface of the disk. Finally, we require the in-plane divergence at the interface to be zero:

(3.2)\begin{equation} \boldsymbol{\nabla}_{s} \boldsymbol{\cdot} \boldsymbol{v}_\parallel|_{z=0} = 0. \end{equation}

Together with the incompressibility equation (2.1b), it follows that

(3.3)\begin{equation} \partial_z v_z|_{z=0} = 0. \end{equation}

In order to have a more uniform boundary condition at the surface, we arbitrarily extend the surface tension difference $\gamma$ to the region underneath the disk ($\rho < a$) while requiring continuity at $\rho =a$. At the same time, we introduce a transformed force density:

(3.4)\begin{equation} \boldsymbol{F}_\parallel=\begin{cases} \boldsymbol{f}_\parallel-\boldsymbol{\nabla}_{s} \gamma & \text{for}\ \rho< a\\ \boldsymbol{0} & \text{for}\ \rho>a . \end{cases} \end{equation}

With this redefinition, (3.1) obtains the unified form

(3.5)\begin{equation} \mu \partial_z \boldsymbol{v}_\parallel|_{z=0}= \boldsymbol{F}_\parallel+ \boldsymbol{\nabla}_{s} \gamma. \end{equation}

Equation (3.2) is also trivially satisfied under the disk and therefore in the whole plane.

To determine the solution for the flow velocity and pressure fields, we use a two-dimensional Fourier transform along the $x$ and $y$ directions. We define the forward Fourier transform of a given function $f(\boldsymbol {\rho }, z)$ as

(3.6)\begin{equation} \tilde{f}(\boldsymbol{k},z) = \mathscr{F}\{\,f(\boldsymbol{\rho},z)\} = \int_{\mathbb{R}^2} \mathrm{d}^2 \boldsymbol{\rho} \, f(\boldsymbol{\rho},z) \exp({-{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho}}), \end{equation}

where $\boldsymbol {k}$ represents the wavevector. The $z$ dependence is left unchanged by the transformation.

Similarly, we define the inverse Fourier transform as

(3.7)\begin{equation} f(\boldsymbol{\rho},z) = \mathscr{F}^{-1}\{\,\tilde{f}(\boldsymbol{k},z)\} = \frac{1}{(2{\rm \pi})^2} \int_{\mathbb{R}^2} \mathrm{d}^2 \boldsymbol{k}\, \tilde{f}(\boldsymbol{k}, z) \exp({{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho}}). \end{equation}

In these equations, $\boldsymbol {\rho } = (x, y)$ denotes the position vector along the interface. We also introduce the wavenumber $k = |\boldsymbol {k}|$, which represents the magnitude of the wavevector, and define the unit vector $\hat {\boldsymbol {k}} = \boldsymbol {k}/k$.

3.1.1. Normal velocity

By forming the divergence of (2.1a) and applying the incompressibility condition, it follows that the pressure field is harmonic, $\nabla ^2 p = 0$ (Happel & Brenner 1983). Thus, the velocity field adheres to the biharmonic equation $\nabla ^4 \boldsymbol {v} = 0$, which can be represented in Fourier space as

(3.8)\begin{equation} (\partial_z^2 - k^2)^2 \tilde{\boldsymbol{v}} = \boldsymbol{0}. \end{equation}

This equation represents a homogeneous fourth-order linear differential equation for $\tilde {\boldsymbol {v}}$ and is solved by

(3.9)\begin{equation} \tilde{\boldsymbol{v}} = (\boldsymbol{\alpha}_1 + z \boldsymbol{\alpha}_2) \,{\rm e}^{kz}, \end{equation}

with the unknown wavenumber-dependent vector functions $\boldsymbol {\alpha }_1$ and $\boldsymbol {\alpha }_2$ that will be determined from the boundary conditions. We have retained only the solution with the decaying exponential to meet the regularity condition of finite velocity as $z \to -\infty$.

Given that both the normal velocity at the interface, $\tilde {v}_z$, and its derivative with respect to $z$, $\partial _z \tilde {v}_z$, are zero at the interface, cf. (3.3), we can conclude that $\tilde {v}_z = 0$ throughout. Therefore, the velocity field only has an in-plane component, denoted henceforth as $\tilde {\boldsymbol {v}}_\parallel$.

3.1.2. In-plane velocity

In the following we determine the solution for the in-plane component $\tilde {\boldsymbol {v}}_\parallel$ of the velocity. In Fourier space, the flow equations (2.1) are projected as

(3.10a)$$\begin{gather} - {\rm i}\boldsymbol{k} \tilde{p} + \mu (\partial_z^2 - k^2) \tilde{\boldsymbol{v}}_\parallel= \boldsymbol{0}, \end{gather}$$

(3.10b)$$\begin{gather}\boldsymbol{k} \boldsymbol{\cdot} \tilde{\boldsymbol{v}}_\parallel= 0 \end{gather}$$

and the force balance at the interface (3.5) as

(3.11)\begin{equation} \mu \partial_z \tilde{\boldsymbol{v}}_\parallel|_{z=0} = \tilde{\boldsymbol{F}}_\parallel+ {\rm i}\boldsymbol{k} \tilde{\gamma}. \end{equation}

By multiplying both sides of (3.10a) with $\hat {\boldsymbol {k}}$ and using (3.10b), we see that the pressure vanishes, $\tilde {p} = 0$. Consequently, (3.10a) simplifies to

(3.12)\begin{equation} (\partial_z^2 - k^2) \tilde{\boldsymbol{v}}_\parallel= \boldsymbol{0}, \end{equation}

which is a homogeneous second-order linear differential equation. Its solution that does not diverge for $z\to -\infty$ is given by

(3.13)\begin{equation} \tilde{\boldsymbol{v}}_\parallel= \boldsymbol{\alpha}_3 \,{\rm e}^{kz}, \end{equation}

where $\boldsymbol {\alpha }_3$ is a wavenumber-dependent function that can be determined from the boundary condition (3.11) as

(3.14)\begin{equation} \boldsymbol{\alpha}_3 = \frac{1}{\mu k} (\tilde{\boldsymbol{F}}_\parallel+ {\rm i}\boldsymbol{k} \tilde{\gamma}). \end{equation}

The condition of vanishing in-plane divergence at the interface (3.2) reads $\boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\alpha }_3 = 0$ in Fourier space. It allows us to express the surface tension with the force density as

(3.15)\begin{equation} \tilde{\gamma} = \frac{{\rm i}\boldsymbol{k} \boldsymbol{\cdot} \tilde{\boldsymbol{F}}_\parallel}{k^2}. \end{equation}

Finally, by substituting (3.15) into (3.14) and (3.13), the solution for the velocity is

(3.16)\begin{equation} \tilde{\boldsymbol{v}}_\parallel= \frac{{\rm e}^{kz}}{\mu k} (\boldsymbol{I} - \hat{\boldsymbol{k}} \hat{\boldsymbol{k}})\boldsymbol{\cdot} \tilde{\boldsymbol{F}}_\parallel. \end{equation}

The velocity therefore only has a component that is parallel to the surface and transverse to the wavevector. A more general solution allowing for a finite thickness of the fluid layer and surface viscosity can be found, for instance, in Martínez-Prat et al. (2021).

3.2. Solution for the circular disk

To solve the flow around a circular disk, we use a cylindrical coordinate system and express the unit wavevector $\hat {\boldsymbol {k}}$ as $\hat {k}_x = \cos \phi$ and $\hat {k}_y = \sin \phi$. In addition, we define the unit vector $\hat {\boldsymbol {t}}$ along the direction perpendicular to the wavevector such that $\hat {t}_x = \sin \phi$ and $\hat {t}_y = -\cos \phi$. Accordingly, the Fourier transform of the surface force density can be expressed in the unit vector basis composed of $\hat {\boldsymbol {k}}$ and $\hat {\boldsymbol {t}}$ as

(3.17)\begin{equation} \tilde{\boldsymbol{F}}_\parallel= \tilde{F}_l \hat{\boldsymbol{k}} + \tilde{F}_t \hat{\boldsymbol{t}}, \end{equation}

wherein $\tilde {F}_l$ and $\tilde {F}_t$ stand for the longitudinal and transverse components of the surface force density, respectively (Bickel 2007). Consequently, we can express the in-plane component of the flow velocity in Fourier space by referring to (3.16) as

(3.18)\begin{equation} \tilde{\boldsymbol{v}}_\parallel= \frac{{\rm e}^{kz}}{\mu k} \tilde{F}_t \hat{\boldsymbol{t}}. \end{equation}

Without loss of generality we assume that the disk moves along the positive $x$ direction. We anticipate the angular structure of the Fourier-transformed surface force density to take the following form:

(3.19a)$$\begin{gather} \tilde{F}_l (k,\phi) = \tilde{M}(k) \cos\phi, \end{gather}$$

(3.19b)$$\begin{gather}\tilde{F}_t (k,\phi) = \tilde{N}(k) \sin\phi, \end{gather}$$

where $\tilde {M} (k)$ and $\tilde {N} (k)$ are unknown wavenumber-dependent functions to be determined from the underlying boundary conditions. We made this specific selection because the inverse Fourier transform will result in an angular structure akin to that required for matching in the boundary conditions imposed at the surface of the disk. Inverse Fourier transform of (3.19) yields the radial and tangential components of the surface force density as

(3.20a)$$\begin{gather} F_\rho = \frac{\cos\theta}{4{\rm \pi}} \int_0^\infty k \, \mathrm{d} k [(\tilde{M}(k)+\tilde{N}(k)) \mathrm{J}_0(\rho k) - (\tilde{M}(k)-\tilde{N} (k)) \mathrm{J}_2(\rho k)], \end{gather}$$

(3.20b)$$\begin{gather}F_\theta =-\frac{\sin\theta}{4{\rm \pi}} \int_0^\infty k \, \mathrm{d} k [(\tilde{M}(k)+\tilde{N} (k))\mathrm{J}_0(\rho k) + (\tilde{M}(k)-\tilde{N}(k)) \mathrm{J}_2(\rho k)]. \end{gather}$$

Likewise, the velocities induced by this force distribution transform to real space as

(3.21a)$$\begin{gather} v_\rho = \frac{\cos\theta}{4{\rm \pi}\mu} \int_0^\infty \mathrm{d} k\,\tilde{N} (k) (\mathrm{J}_0(\rho k) + \mathrm{J}_2(\rho k)) \,{\rm e}^{kz}, \end{gather}$$

(3.21b)$$\begin{gather}v_\theta =-\frac{\sin\theta}{4{\rm \pi}\mu} \int_0^\infty \mathrm{d} k \, \tilde{N} (k) (\mathrm{J}_0(\rho k) - \mathrm{J}_2(\rho k)) \,{\rm e}^{kz}. \end{gather}$$

The surface tension, given by (3.15), can be expressed in terms of the longitudinal component of the surface force density as

(3.22)\begin{equation} \tilde{\gamma} (k, \phi) = \frac{{\rm i}}{k} \tilde{F}_l (k, \phi) = \frac{{\rm i}}{k} \tilde{M} (k) \cos\phi, \end{equation}

which can be transformed back to real space through an inverse transformation as

(3.23)\begin{equation} \gamma =-\frac{\cos \theta}{2{\rm \pi}} \int_0^\infty \mathrm{d}k \, \tilde{M} (k) \mathrm{J}_1(\rho k). \end{equation}

3.2.1. Force density

After representing the flow velocity as integrals involving wavenumber-dependent functions, the next step involves the formulation of dual integral equations. The integral equations for the inner domain are obtained from the no-slip boundary conditions at the surface of the disk, $v_\rho = V \cos \theta$ and $v_\theta = -V\sin \theta$ for $z=0$ and $\rho \in [0, a]$:

(3.24a,b)\begin{equation} \int_0^\infty \mathrm{d}k \, \tilde{N}(k) \mathrm{J}_0(\rho k) = 4{\rm \pi}\mu V,\quad \int_0^\infty \mathrm{d}k \, \tilde{N}(k) \mathrm{J}_2(\rho k) = 0. \end{equation}

The integral equations for the outer problem are obtained by imposing a vanishing surface force density (3.20) at the surface outside the disk for $z=0$ and $\rho > a$:

(3.25a,b)\begin{equation} \int_0^\infty k \, \mathrm{d}k (\tilde{M}(k) + \tilde{N}(k)) \mathrm{J}_{0}(\rho k) = 0,\quad \int_0^\infty k \, \mathrm{d}k (\tilde{M}(k) - \tilde{N}(k)) \mathrm{J}_{2}(\rho k) = 0. \end{equation}

Equations (3.24a,b) and (3.25a,b) form a system of dual integral equations on the inner and outer domain boundaries. The solution can be obtained using established techniques as outlined in the works of Sneddon (1960, 1966) and Copson (1947, 1961). The basic idea involves seeking a set of solutions that satisfy the equations for the outer problems (3.25a,b). Typically, the solution for the unknown wavenumber-dependent functions is explored through definite integrals, often resulting in vanishing values in the outer domain when $\rho > a$. For detailed solution approaches, one can refer to Chapter IV of Sneddon (1966) for dual integral equations and Chapter VIII for their applications in electrostatics. When incorporating these solution forms into the equations for the inner problem (3.24a,b), one encounters a system of Fredholm integral equations, the solution of which is not always straightforward to obtain. In our case, we employ the power series method to derive the solution, whereby the unknown function is expanded in terms of unspecified coefficients. Subsequently, we evaluate the resulting integrals analytically. Upon determining the coefficients through comparison with the known right-hand side, the unknown function is then discerned from its series expansion. For further details, the reader is referred to the appendix of Daddi-Moussa-Ider et al. (2020a) where a similar approach was employed. However, in similar scenarios, it is anticipated that the solution may involve combinations of trigonometric and Bessel functions. Therefore, employing a guessed solution can often be advantageous in the given context. For instance, this approach has been applied to the flow generated by different types of force or source singularities near circular interfaces (Daddi-Moussa-Ider et al. 2021b; Daddi-Moussa-Ider, Vilfan & Golestanian 2022; Daddi-Moussa-Ider et al. 2023b). In our case, the dual integral equations are solved by

(3.26a)$$\begin{gather} \tilde{M}(k) = \frac{16\mu V}{k} \mathrm{J}_1(ka), \end{gather}$$

(3.26b)$$\begin{gather}\tilde{N}(k) = \frac{8\mu V}{k} \sin (ka). \end{gather}$$

The solutions we derived can now be inserted into the integrals that stem from the inverse Fourier transforms, (3.20), (3.21) and (3.23). The relative surface tension follows from (3.23):

(3.27)\begin{equation} \gamma =-\frac{4\mu V}{\rm \pi} \cos\theta \times \begin{cases} \dfrac{\rho}{a} & \text{if}\ \rho < a, \\ \dfrac{a}{\rho} & \text{if}\ \rho > a. \end{cases} \end{equation}

The gradient of the surface tension around the disk is therefore proportional to the inverse square of the distance. The expression for the interior, of course, only represents the surface tension arbitrarily extended across the disk with no physical meaning. By inserting it into (3.4), it gives the physical force density exerted on the fluid by the disk:

(3.28)\begin{equation} \boldsymbol{f}_\parallel= \boldsymbol{F}_\parallel-\frac{4}{\rm \pi} \frac{\mu V}{a} \hat{\boldsymbol{e}}_x. \end{equation}

After carrying out the integrals given in (3.20), we can express the resulting force density at the surface of the disk in Cartesian coordinates as

(3.29)\begin{equation} \boldsymbol{f}_\parallel= \frac{4\mu V}{{\rm \pi} a} \{[A(\rho) + B(\rho)\cos (2\theta)] \hat{\boldsymbol{e}}_x+ B(\rho) \sin (2\theta) \hat{\boldsymbol{e}}_y\}, \end{equation}

where we have defined the dimensionless radial functions

(3.30a)$$\begin{gather} A(\rho) = \frac{a}{2 \sqrt{a^2-\rho^2}}, \end{gather}$$

(3.30b)$$\begin{gather}B(\rho) =- \left(\frac{a}{\rho}\right)^2 \biggl(\frac{2a^2-\rho^2}{2a \sqrt{a^2-\rho^2}} - 1\biggr). \end{gather}$$

It can be seen that $A$ is a monotonically increasing function of $\rho$ that starts with the value of $1/2$ at the origin, while $B$ is a monotonically decreasing function of $\rho$, starting from zero at the origin. It is evident that both $A$ and $-B$ exhibit asymptotic divergence, scaling to leading order as $1/\sqrt {a-\rho }$ as $\rho$ approaches $a$. The maximum magnitude of the surface force density is attained at $A - B$ when $\theta \in \{{\rm \pi} /2, 3{\rm \pi} /2\}$, while the minimum magnitude occurs at $A + B$ when $\theta \in \{0, {\rm \pi}\}$.

The drag force on the disk can be calculated as the sum of the shear forces at the bottom surface and the surface tension on the perimeter:

(3.31)\begin{equation} \boldsymbol{F}_{SI}=-\int_{\mathcal{S}} {\rm d} S \, \boldsymbol{f}_\parallel- \int_\ell {\rm d} s \,\gamma \hat{\boldsymbol{n}}. \end{equation}

The contribution of surface tension (3.27) amounts to $-4\mu a V$. The integral of the surface force density (3.29), to which only the term $A(\rho )$ contributes, equally gives $-4\mu a V$. Together, the drag coefficient of a thin circular disk in a fluid with surface incompressibility (2.15) is obtained as

(3.32)\begin{equation} R_{SI} = 8\mu a \end{equation}

and consists of equal contributions by the shear stress and by the surface tension. The drag coefficient agrees with the limiting cases from the expressions by Hughes et al. (1981) and Stone & Ajdari (1998).

3.3. Flow field

To derive a closed analytical expression for the flow field, determined by (3.21), we define the sequence of infinite integrals:

(3.33)\begin{equation} C_n = \int_0^\infty \mathrm{d}k \, \frac{\sin (ka)}{k}\mathrm{J}_n(\rho k) \,{\rm e}^{kz}, \end{equation}

for $z \le 0$ to ensure convergence. It can be shown that $C_0 = \operatorname {arg} \{ Q \}$ and $C_2 = \operatorname {Im} \{ s Q \}$, where we have defined the abbreviation $Q = s + S$ with $s=(z+{\rm i}a)/\rho$ and $S = \sqrt {1+s^2}$. Here, ‘arg’ represents the argument of a complex number, while ‘Im’ signifies its imaginary part. To obtain these results, it is sufficient to represent the sine function using Euler's notation and apply the Laplace transform to the Bessel function (Watson 1922). For a detailed treatment of integrals with similar forms, see the appendix in Daddi-Moussa-Ider et al. (2020b). The radial and azimuthal components of the flow field can then be cast in a compact form as

(3.34a)$$\begin{gather} \frac{v_\rho}{V} = \frac{2}{\rm \pi} (C_0 + C_2) \cos\theta, \end{gather}$$

(3.34b)$$\begin{gather}\frac{v_\theta}{V} =-\frac{2}{\rm \pi} (C_0 - C_2) \sin\theta. \end{gather}$$

At the boundary, where $z=0$, a closed-form expression for $C_n$ can be derived for any arbitrary $n \ge 0$ as

(3.35) \begin{equation} C_n (z=0) = \begin{cases} \dfrac{1}{n} \biggl(\dfrac{\rho}{a + \sqrt{a^2-\rho^2}}\biggr)^n \sin \left(\dfrac{n{\rm \pi}}{2}\right) & \text{if}\ \rho < a, \\ \dfrac{1}{n} \sin \left[n \arcsin \left(\dfrac{a}{\rho}\right)\right] & \text{if}\ \rho > a. \end{cases} \end{equation}

In particular, for $\rho < a$, it can be noticed that $C_{2n}(z=0) = 0$ for $n \ge 1$. For $n=0$, it follows that $C_0(z=0) = {\rm \pi}/2$ if $\rho < a$ and $C_0(z=0) = \arcsin ( a/\rho )$ if $\rho > a$.

The flow field, transformed into the co-moving frame ($(\boldsymbol {v}-\hat {\boldsymbol {e}}_x V)/V$), is shown in figure 3.

Figure 3. The flow field of a circular disk embedded in an incompressible interface (SI problem), shown in the co-moving frame. (a) The top view at the surface and (b) the side view at $y=0$.

4. Thin circular disk translating in a fluid with force-free surface

In the following, we recapitulate the solution for the flow and the drag force on a thin circular disk, moving horizontally at a gas–liquid interface. Because the interface is assumed to be planar and is not affected by the flow, the solution is mathematically equivalent to the edgewise motion of a thin disk in bulk fluid, with half the drag coefficient. The latter has been solved by Ray (1936) using Bessel functions. The flow around a thin disk can also be obtained from the solution for an ellipsoidal particle in the limit of zero thickness (Lamb 1945). A general approach for analysing the arbitrary motion of a circular disk in a Stokes flow has later been formulated using dual integral equation approach (Tanzosh 1994; Tanzosh & Stone 1996). The problem of a disk moving sideways between parallel walls has further been investigated by Davis (1991), who determined the additional drag encountered by the disk by solving a pair of integral equations (Davis 1991). In order to ensure uniform notation needed for the superposition of both solutions, we give a brief summary of the of the derivation by Tanzosh & Stone (1996) using our notation in the following. We will omit the ‘FS’ subscript as all quantities in this section are related to a disk in motion at a free interface.

We express the forces in Fourier space using (3.19). Since this problem has been thoroughly examined in the literature, we refrain from presenting the complete solution steps here. In the case of a free surface, our expectation is that the force density aligns along the direction of motion. Unlike in the case of an incompressible surface, where $\tilde {\boldsymbol {v}}_\parallel$ is found to be dependent solely on $\tilde {F}_t$, see (3.18), the absence of surface tension in the in-plane force balance introduces a different dependence. In this scenario, the in-plane velocity is determined by both $\tilde {F}_t$ and $\tilde {F}_l$. By anticipating that the force is aligned along the $x$ direction, one would posit $\tilde {M}(k) = \tilde {N}(k)$. It is worth noting that a systematic investigation without enforcing this condition leads to the same conclusion. By requiring no-slip boundary conditions at the surface of the disk, the corresponding dual integral integrals for the inner domain for $z=0$ and $\rho \in [0,a]$ are

(4.1a,b)\begin{equation} \int_0^\infty \mathrm{d}k \, \tilde{M}(k) \mathrm{J}_{0}(\rho k) = \frac{8}{3} {\rm \pi}\mu V,\quad \int_0^\infty \mathrm{d}k \, \tilde{M}(k) \mathrm{J}_{2}(\rho k) = 0. \end{equation}

Equations (3.25a,b) on the outer domain remain the same as these are connected to the absence of force beyond the disk. Then, the solution of the resulting dual integral equations is obtained as

(4.2)\begin{equation} \tilde{M}(k) = \tilde{N}(k) = \frac{16 }{3} \frac{\mu V}{k} \sin (ka). \end{equation}

In contrast to the case of incompressible surface, the expression for the surface force density assumes a straightforward form and can be derived from (3.20) as

(4.3)\begin{equation} \boldsymbol{f}_\parallel= \frac{8\mu V}{3{\rm \pi} \sqrt{a^2-\rho^2}} \hat{\boldsymbol{e}}_x. \end{equation}

Integrating this force over the surface of the disk leads to the familiar expression for the drag coefficient (2.16):

(4.4)\begin{equation} R_{FS} = \tfrac{16}{3} \mu a. \end{equation}

The hydrodynamic flow field is determined by carrying out the inverse Fourier transform, resulting in the components of the velocity field as follows:

(4.5a)$$\begin{gather} \frac{v_\rho}{V} = \frac{2}{3{\rm \pi}} \Bigl[ 3 C_0 + C_2 + (D_0 - D_2) \frac{z}{\rho} \Bigr] \cos\theta, \end{gather}$$

(4.5b)$$\begin{gather}\frac{v_\theta}{V} =-\frac{2}{3{\rm \pi}} \Bigl[ 3 C_0 - C_2 + (D_0 + D_2) \frac{z}{\rho} \Bigr] \sin\theta, \end{gather}$$

(4.5c)$$\begin{gather}\frac{v_z}{V} = \frac{4}{3{\rm \pi}} \frac{z}{\rho} D_1 \cos\theta. \end{gather}$$

The sequence $C_n$ was previously defined in terms of an infinite integral by (3.33). Here, $D_n = \rho \partial _z C_n$, also a dimensionless function. It can be shown that $D_0 = \operatorname {Im} \{ 1/S \}$, $D_1 = \operatorname {Im} \{s/ S\}$ and $D_2 = \operatorname {Im} \{Q^2/ S\}$. When $z=0$, a closed-form expression for $D_n$ can likewise be derived for any $n \ge 0$ as

(4.6)\begin{equation} D_n (z=0) = \begin{cases} \dfrac{\rho}{\sqrt{a^2-\rho^2}} \biggl(\dfrac{\rho}{a + \sqrt{a^2-\rho^2}} \biggr)^n \cos \left( \dfrac{n{\rm \pi}}{2} \right) & \text{if}\ \rho < a, \\ \dfrac{\rho}{\sqrt{\rho^2-a^2}} \sin \left[ n \arcsin \left( \dfrac{a}{\rho} \right) \right] & \text{if}\ \rho > a. \end{cases} \end{equation}

In particular, for $z=0$ and $\rho < a$, all odd terms vanish, $D_{2n+1} = 0$.

Figure 4 shows the resulting flow at the surface and in a vertical cross-section. Unlike with an interface with surface incompressibility, where the pressure is zero, the solution for a free surface includes a non-zero pressure field, which is related to the normal velocity as $p = 2\mu v_z/z$.

Figure 4. The flow field of a circular disk embedded in a force-free interface (FS problem), shown in the co-moving frame in top view (a) and side view (b).

5. The optimal propulsion of a disk-shaped Marangoni surfer

With the drag coefficients of the two passive problems, the respective translational velocities in the passive problems follow from (2.13) as

(5.1a,b)\begin{equation} \frac{V_{FS}}{V_{A}} = \frac{R_{SI}}{R_{SI} - R_{FS}}=3,\quad \frac{V_{SI}}{V_{A}} = \frac{R_{FS}}{R_{SI} - R_{FS}}=2. \end{equation}

The full solution for the optimal active surfer can now be obtained as the superposition of the solution for the passive body in an incompressible surface (§ 3) with the negative of the solution for the free surface (§ 4). We thereby replace the body velocity in the former case with $V_{SI}=2V_{A}$ and in the latter case with $V_{FS}=3V_{A}$. The relative surface tension of the optimal active surfer (figure 5a) follows from (3.27) as $\gamma _{A}=-\gamma _{SI}$:

(5.2)\begin{equation} \gamma_{A} = \frac{8\mu V_{A}}{\rm \pi} \frac{a}{\rho} \cos\theta. \end{equation}

The total force on the contact line then equals $8\mu a V_{A}$. Interestingly the force is by a factor of $3/2$ larger than the force that would be needed to pull a passive disk with velocity $V_{A}$, as already stated by Lauga & Davis (2012). Note that the fact that the ratio $3/2$ is the same as the ratio $R_{SI}/R_{FS}$ is specific to the disk geometry and not generally valid.

Figure 5. (a) The surface tension $\gamma _{A}$ of the optimal circular Marangoni surfer (colour scale) relative to the unperturbed surface. The arrows indicate the force density on the fluid, $\boldsymbol {\nabla }_{s}\gamma _{A}$ outside the disk and $\boldsymbol {f}_\parallel$ underneath it. The flow field in the co-moving frame in top view (b) and side view (c).

The force density with which the bottom surface of the surfer acts on the fluid can likewise be obtained as ${\boldsymbol {f}_\parallel }_{FS} - {\boldsymbol {f}_\parallel }_{SI}$. The surface integral of this force over the disk equals the force caused by the surface tension, in accordance with the force-free condition on the active surfer.

Finally, the flow field induced by the optimal Marangoni surfer is obtained as

(5.3a)$$\begin{gather} \frac{v_\rho}{V_{A}} = \frac{2}{\rm \pi} \left( C_-+ D_- \frac{z}{\rho} \right) \cos\theta, \end{gather}$$

(5.3b)$$\begin{gather}\frac{v_\theta}{V_{A}} =-\frac{2}{\rm \pi} \left( C_++ D_+ \frac{z}{\rho} \right) \sin\theta, \end{gather}$$

(5.3c)$$\begin{gather}\frac{v_z}{V_{A}} = \frac{4}{\rm \pi} \frac{z}{\rho} D_1 \cos \theta, \end{gather}$$

where we have introduced the abbreviations $C_\pm = C_0 \pm C_2$ and $D_\pm = D_0 \pm D_2$.

At the interface where $z=0$ and $\rho \ge a$, the velocity can be expressed as

(5.4)\begin{equation} \frac{{\boldsymbol{v}}_\parallel}{V_{A}} =\frac{2}{\rm \pi}(C_- \cos\theta \hat{\boldsymbol{e}}_\rho - C_+ \sin\theta \hat{\boldsymbol{e}}_\theta), \end{equation}

with

(5.5)\begin{equation} C_\pm= \arcsin \left( \frac{a}{\rho} \right) \pm \frac{a \sqrt{\rho^2-a^2}}{\rho^2} . \end{equation}

The resulting flow field is shown in figure 5(b,c). The pressure can likewise be written in a closed form as

(5.6)\begin{equation} p(\rho, z=0)= \frac{8\mu V_{A}}{{\rm \pi}\rho} \frac{a \cos\theta}{\sqrt{\rho^2-a^2}} \end{equation}

for $\rho \ge a$ and is zero underneath the surfer.

6. Discussion

We have derived a minimum dissipation theorem that gives a lower bound on hydrodynamic dissipation by a Marangoni surfer of arbitrary shape in the limit of low Reynolds and capillary numbers. In our calculation we assumed a planar gas–liquid interface, but the expansion to any other interface form, for example, taking into account the meniscus around the contact line, is straightforward as long as the shape is not influenced by the flow. Likewise, the problem can be solved with additional confinement for the liquid phase where an interesting reversal of swimming direction can occur (Vandadi, Kang & Masoud 2017). Another generalisation would be to include the rotational motion of the surfer. In all those cases the only difficulty is to calculate the drag coefficients of the corresponding passive problems.

As a specific example, we calculated the efficiency limit for a surfer taking the form of a thin circular disk. The obtained value $\eta _{L}=1/3$ is significantly higher than that for any other realistic swimmer designs. Interestingly, Lauga & Davis (2012) have already discussed exactly the same solution that we have now proven to be hydrodynamically optimal. Their approach is, however, very different from ours as it is based on the solution of the Laplace equation for the diffusion of surfactants around the surfer in the limit of low Péclet numbers. They have shown that only the first angular mode, corresponding to a dipolar source of surfactants, contributes to propulsion. We have now shown that this mode coincides with the solution of the flow optimisation problem. The two solutions are not entirely comparable, because at low Péclet numbers surfactant diffusion makes the dominant contribution to the total dissipation. On the other hand, the finding that surfactant distributions resulting from diffusion lead to hydrodynamically optimal solutions also has a broader validity, beyond disk-shaped bodies. From the structure of our Green's functions one can see that the optimal distribution of surface tension indeed satisfies a two-dimensional Laplace equation around any flat surfer. We expect, however, that this is not the case for general surfers with a finite vertical extension.

While we are now able to determine the solution with the minimal hydrodynamic dissipation, losses of free energy also take place due to surfactant diffusion. If the surface tension is related to the surface density of insoluble surfactants $c$ as $\gamma =\gamma ^0 - K c$ and the dynamics of the latter is determined by the advection–diffusion equation (Schwartz & Roy 2001)

(6.1)\begin{equation} \partial_t c =- \boldsymbol\nabla_{s}\boldsymbol{\cdot} (\boldsymbol{v} c - D\boldsymbol{\nabla}_{s} c), \end{equation}

then the rate of free energy loss due to diffusion can be written as

(6.2)\begin{equation} P_{D}= K D \int_\mathcal{I} {\rm d} S \frac{(\boldsymbol{\nabla}_{s} c)^2}{c}. \end{equation}

The rate of work done on the fluid (2.4), on the other hand, is $P= K \int _\mathcal {I} {\rm d} S\, c \boldsymbol {\nabla }_{s}\boldsymbol {\cdot } \boldsymbol {v}$. Their ratio is determined by the Péclet number, defined as

(6.3)\begin{equation} Pe=\frac{v a}{D}, \end{equation}

if we assume that the derivatives of $c$ and $\boldsymbol {v}$ are determined by the characteristic length scale $a$. Our theorem therefore gives the total dissipation in the limit of high Péclet numbers. In the general case, it applies to the hydrodynamic dissipation only. At the same time, diffusivity is also needed if, for example, the surfer only emits surfactants at its perimeter, but not elsewhere. This limitation does not apply to all surfers – a counterexample is given by beetles that can secret surfactants from their tail, at a distance from other contact lines (Lang, Seifert & Dettner 2012). Furthermore, our study cannot take into account the energetic cost of synthesising or collecting surfactants. This cost could be significant in the total energy balance of propulsion. Overall, finding optimal solutions for the combined problem including hydrodynamics and diffusion of surfactants remains a challenge for future investigations.