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 Reference Happel and Brenner1983):

(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. Reference Nasouri, Vilfan and Golestanian2021) and later with combined external and internal dissipation (Daddi-Moussa-Ider et al. Reference Daddi-Moussa-Ider, Golestanian and Vilfan2023a). 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. Reference Daddi-Moussa-Ider, Nasouri, Vilfan and Golestanian2021a).

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 Reference Guazzelli and Morris2009). 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. Reference Elfring, Leal and Squires2016; Manikantan & Squires Reference Manikantan and Squires2020). The same model also represents a limiting case of the Boussinesq–Scriven model (Scriven Reference Scriven1960; Manikantan & Squires Reference Manikantan and Squires2017) with an incompressible ‘membrane’, however without additional surface viscosity, i.e. in the limit of vanishing Boussinesq number (Stone & Masoud Reference Stone and Masoud2015).

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 Reference Masoud and Stone2019) 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. Reference Nasouri, Vilfan and Golestanian2021; Daddi-Moussa-Ider et al. Reference Daddi-Moussa-Ider, Golestanian and Vilfan2023a). 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 (Reference Hughes, Pailthorpe and White1981) and Stone & Ajdari (Reference Stone and Ajdari1998):

(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 Reference Happel and Brenner1983), 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 Reference Lighthill1952), 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. Reference Leshansky, Kenneth, Gat and Avron2007), 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 Reference Michelin and Lauga2010; Guo et al. Reference Guo, Zhu, Liu, Bonnet and Veerapaneni2021), which is even higher for elongated swimmers (Nasouri et al. Reference Nasouri, Vilfan and Golestanian2021), 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 (Reference Stone and Masoud2015). 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 Reference Happel and Brenner1983). 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 Reference Stone and Masoud2015). An analytical expression can be obtained by integrating (3.6) from Stone & Masoud (Reference Stone and Masoud2015) 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 Reference Stone and Masoud2015).

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 Reference Happel and Brenner1983). 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. (Reference Martínez-Prat, Alert, Meng, Ignés-Mullol, Joanny, Casademunt, Golestanian and Sagués2021).

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 Reference Bickel2007). 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 (Reference Sneddon1960, Reference Sneddon1966) and Copson (Reference Copson1947, Reference Copson1961). 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 (Reference Sneddon1966) 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. (Reference Daddi-Moussa-Ider, Lisicki, Löwen and Menzel2020a) 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. Reference Daddi-Moussa-Ider, Sprenger, Richter, Löwen and Menzel2021b; Daddi-Moussa-Ider, Vilfan & Golestanian Reference Daddi-Moussa-Ider, Vilfan and Golestanian2022; Daddi-Moussa-Ider et al. Reference Daddi-Moussa-Ider, Hosaka, Vilfan and Golestanian2023b). 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. (Reference Hughes, Pailthorpe and White1981) and Stone & Ajdari (Reference Stone and Ajdari1998).

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 Reference Watson1922). For a detailed treatment of integrals with similar forms, see the appendix in Daddi-Moussa-Ider et al. (Reference Daddi-Moussa-Ider, Sprenger, Amarouchene, Salez, Schönecker, Richter, Löwen and Menzel2020b). 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$.