1. Introduction
In this work, we develop a mathematical description of the free-boundary dynamics of a two-dimensional (2-D) incompressible droplet moving atop a bulk Stokes fluid. Following the approach of Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019), this incompressible surface phase can be either active or passive, and is described by the most general linear isotropic fluid model. This model allows a viscous shear stress, an antisymmetric ‘chiral’ stress reflecting the driven rotation of the fluid constituents and an ‘odd’ Hall stress allowed by the consequent loss of time-reversal symmetry at the microscopic level. Given its generality, this model touches upon both classical and emerging areas of fluid dynamics and applied mathematics, including Langmuir films, mixed dimension boundary value problems, Euler and quasigeostrophic vortex systems and active matter systems. We briefly describe connections with these areas before stating our main results.
The interaction of rotating elements in a fluid is a foundational topic in fluid dynamics, going back to the explication of 2-D point vortices of the Euler equations interacting through the Biot–Savart law (Saffman Reference Saffman1995). A 2-D patch of constant vorticity interacts with itself similarly and its dynamics can be reduced to a free-boundary problem (Pullin Reference Pullin1992). The surface quasigeostrophic equations (SQG) of atmospheric physics have their own singular and free-boundary analogues (Held et al. Reference Held, Pierrehumbert, Garner and Swanson1995; Rodrigo & Fefferman Reference Rodrigo and Fefferman2004). Rotational interaction problems also arise, in the guise of so-called active matter, in the zero Reynolds limit of the Stokes equations where solid particles are driven to rotate by an external field, or through internal actuation. Ensembles of such particles will interact through their induced fluid flows, steric interactions and possibly other fields such as magnetic. These systems can show activity-induced phase separation (Yeo, Lushi & Vlahovska Reference Yeo, Lushi and Vlahovska2015), crystallization and hyperuniformity Petroff, Wu & Libchaber Reference Petroff, Wu and Libchaber2015; Oppenheimer, Stein & Shelley Reference Oppenheimer, Stein and Shelley2019; Oppenheimer et al. Reference Oppenheimer, Stein, Yah Ben Zion and Shelley2022), odd surface waves and edge currents (Soni et al. Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019), complex interactions of vortical structures (Bililign et al. Reference Bililign, Balboa Usabiaga, Ganan, Poncet, Soni, Magkiriadou, Shelley, Bartolo and Irvine2021) and forms of active turbulence (Kokot et al. Reference Kokot, Das, Winkler, Gompper, Aranson and Snezhko2017).
When such many-particle systems are modelled as continuous fluidic materials, novel internal stresses can arise. Firstly, the driven rotation of the fluid's constituents gives rise to an anti-symmetric driving stress. Consequent to the microscopic driving of rotation, these out-of-equilibrium fluids do not obey time-reversal symmetry and so can possess an odd or Hall stress which, in its simplest case, is linear in rates of strain and couples longitudinal and transverse flow components. Examples of such systems include quantum Hall fluids (Avron, Seiler & Zograf Reference Avron, Seiler and Zograf1995), vortex fluids (Wiegmann & Abanov Reference Wiegmann and Abanov2014) and electron fluids in graphene (Berdyugin et al. Reference Berdyugin2019). Fluids with an odd viscosity can exhibit rheological properties and exotic flow phenomena markedly different from their Newtonian counterparts such as unidirectional edge currents or topological waves (Soni et al. Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019; Souslov et al. Reference Souslov, Dasbiswas, Fruchart, Vaikuntanathan and Vitelli2019).
Many of the examples above are of rotor assemblies sitting on a 2-D fluid interface either embedded within, or sitting atop, a viscous fluid bulk. Other active matter systems, particularly active nematics formed of microtubule bundles and molecular motors, have been studied in this geometry both experimentally (see, e.g. Sanchez et al. Reference Sanchez, Chen, DeCamp, Heymann and Dogic2012) and theoretically (see, e.g. Gao et al. Reference Gao, Blackwell, Glaser, Betterton and Shelley2015). These systems have generally been modelled as a 2-D incompressible active material covering the entire surface, and the bulk as an incompressible Stokes fluid driven by the surface shear stress. It has been shown that the bulk flows can profoundly modify the active surface dynamics, for example by introducing new length scales of system instability at the onset of active nematic turbulence (Gao et al. Reference Gao, Betterton, Jhang and Shelley2017; Martínez-Prat et al. Reference Martínez-Prat, Ignés-Mullol, Casademunt and Sagués2019). A complementary line of study concerns the turbulent statistics and intermittency of flows within a flat surface that overlay a turbulent and incompressible 3-D bulk (Goldburg et al. Reference Goldburg, Cressman, Vörös, Eckhardt and Schumacher2001; Cressman et al. Reference Cressman, Davoudi, Goldburg and Schumacher2004). Some part of the complexity of surface flows in this case derives from the 2-D compressibility of the 2-D surface velocity in their reflection of 3-D inertial turbulence.
Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019) studied theoretically the dynamics of active chiral surface droplets where the contribution of the underlying bulk fluid was modelled as a simple local drag term, as is appropriate for the dynamics of large droplets near a solid substrate. This yields a homogeneous Brinkman equation, with activity-driven boundary conditions, for the droplets’ in-plane velocity field. The kinematic boundary condition then evolves the droplet domain. Interactions through the bulk fluid, much less with other droplets, are completely screened in this near-wall limit. Here, we allow full coupling between the surface phase and the bulk fluid subphase, allowing the droplet to interact with itself both through its internal stresses, and through induced 3-D fluid motions. In our formulation we use the Neumann-to-Dirichlet map for the 3-D Stokes equations in a finite-depth layer or half-space (Masoud & Shelley Reference Masoud and Shelley2014) to express the surface velocity as a surface convolution of a singular kernel with the surface stress, with that shear stress produced by the surface phase. This relation is quite general and here gives rise to a difficult and novel free-boundary problem.
As first exploratory problems, we restrict our study here to axisymmetric solutions. To study rotational flows we determine the activity-driven flows within circular droplets and in the bulk. To study moving interfaces, we also study domains with holes (annuli) and study the course of hole closure and how it is affected by system parameters. Our analyses and computations thereof show that the bulk surface shear stresses diverge as an inverse square root at the droplet boundary. Nonetheless, despite the divergence, the in-plane velocities remain bounded and continuous. For the disk, this divergence is associated with the rotational drive. Singular flows arise in other rotational systems, such as for an infinitely thin solid disk rotating in a Stokes fluid (Jeffery Reference Jeffery1915) whose edge shear stress diverges similarly, and for the SQG vortex patch problem (Rodrigo & Fefferman Reference Rodrigo and Fefferman2004), which exhibits logarithmic divergences in tangential surface velocity. An identical logarithmic divergence is found, in a continuum limit, for planar assemblies of rotating particles rotating in a Stokes fluid (Yan et al. Reference Yan, Corona, Malhotra, Veerapaneni and Shelley2020).
This work also adds to classic work in applied mathematics on the solution of 3-D mixed boundary value problems in potential theory, which arise when considering elliptic problems on two or more different domains, each of which has a different boundary condition that must be satisfied. Such problems, in axisymmetric settings, are frequently converted into multiple integral equations, and powerful methods have been developed to extract their near-analytical solutions; see Sneddon (Reference Sneddon1966) for an overview. These techniques are not readily applicable as our inhomogeneous forcing does not come from a prescribed stress or velocity field on any domain. Since the boundary includes a 2-D fluid monolayer, there is an additional 2-D elliptic problem with its own higher codimension boundary conditions coupled to the base 3-D one. Other authors have more recently examined related problems in similar geometries but included simplifying assumptions such as the complete absence of a vertical flow component (Stone & McConnell Reference Stone and McConnell1995) or a vanishing monolayer viscosity that reduces the order of the monolayer equation (Alexander et al. Reference Alexander, Bernoff, Mann, Mann and Zou2006). Here, we preserve full generality and formulate a dual integral equation – in addition to the Green's function formulation – to obtain a near-analytical solution for a circular droplet, and demonstrate the formulation as triple integral equations for flow in an annulus.
We begin by giving the mathematical formulation and governing equations in § 2. The solution to the general mathematical problem is stated in terms of a Green's function. We then proceed to specialize the formulation to the axisymmetric case in § 3. In §§ 4 and 5, we demonstrate the solutions to the discal and annular geometries. An appendix reviews the experimental system and parameter values, lists the non-dimensional groups associated with the parameters and considers the solution in the infinite strip geometry.
2. A mathematical model
We consider a surface phase domain 
$\mathcal {D}$ on the upper surface, 
$z=0$, of a layer of passive 3-D Stokes fluid (subphase) of depth 
$H$ and infinite extent in the 
$x$ and 
$y$ directions (figure 1a). Gravitational forces and curvature of the interface are ignored, and we assume that the vertical velocity vanishes at the surface 
$z=0$. At 
$z=-H$, the 3-D fluid subphase is in contact with a wall where it satisfies a no-slip condition. We assume that the 3-D velocity field 
$\boldsymbol {u}$ and pressure field 
$p$ of the subphase satisfy the incompressible Stokes equations
\begin{equation} -\boldsymbol{\nabla}_{3D} p + \mu \boldsymbol{\nabla}_{3D}^2\boldsymbol{u} = 0 \quad \text{and} \quad \boldsymbol{\nabla}_{3D} \boldsymbol{\cdot} \boldsymbol{u}= 0,\end{equation}
where 
$\mu$ is the viscosity of the subphase and 
$\boldsymbol {\nabla }_{3D}$ is the 3-D gradient operator 
$(\partial /\partial x, \partial /\partial y, \partial /\partial z)$.
Figure 1. (a) Schematic of a 2-D monolayer domain 
$\mathcal {D}$ situated on the upper surface of a Stokesian sublayer of depth 
$H$, whose motion generates a shear stress on the monolayer. Outside of 
$\mathcal {D}$, the surface is stress free. The inhomogeneous forcing arises from line tension as well as rotational stresses at the boundary of 
$\mathcal {D}$, represented by coloured triangles. (b) Top-down view of an annular domain with inner radius 
$R_i$ and outer radius 
$R_o$ and the corresponding axisymmetric surface flow field 
$\boldsymbol {U}$. Note that the rotational traction vectors follow the direction of the local tangent vector.
Let 
$\boldsymbol {U}$ be the 2-D fluid velocity field in the 
$z=0$ plane. The surface and bulk velocities are related by continuity: 
$\boldsymbol {u}(x,y,z=0) = (\boldsymbol {U}(x,y),0)$, a notation that captures the condition that the surface remains flat. Following Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019), we take the surface phase in 
$\mathcal {D}$ to be described by a general incompressible and isotropic 2-D fluid with linear viscous and Hall stresses, and driven by an anti-symmetric stress. And so, firstly, we have
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U} = 0, \quad \boldsymbol{x}\in \mathcal{D},\end{equation}
where 
$\boldsymbol {\nabla }=(\partial /\partial x , \partial /\partial y)$ is the 2-D gradient on the surface. Secondly, the stress tensor, 
$\boldsymbol {\sigma }$, of the 2-D active surface phase takes the form from Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019)
\begin{equation} \boldsymbol{\sigma} ={-} P{\boldsymbol{\mathsf{I}}} + \eta (\boldsymbol{\nabla} \boldsymbol{U} + \boldsymbol{\nabla} \boldsymbol{U}^T) + \eta_R(2\varOmega - \omega) {\boldsymbol{\mathsf{J}}} + \eta_O(\boldsymbol{\nabla}^\perp \boldsymbol{U} + \boldsymbol{\nabla} \boldsymbol{U}^\perp),\end{equation}
where 
$P$ is the planar pressure enforcing that 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {U} = 0, \omega = \hat {\boldsymbol {z}} \boldsymbol {\cdot } (\boldsymbol {\nabla }_{3D}\times \boldsymbol {U})$ is the scalar vorticity and 
$\varOmega$ is the rotation frequency of the external magnetic field, taken to be spatially uniform and time independent. The tensors 
${\boldsymbol{\mathsf{I}}}$ and 
${\boldsymbol{\mathsf{J}}}$ are the 2-D identity and anti-symmetric Levi-Civita tensors, respectively, while the operator 
$\perp$ maps a 2-D vector to its rotation by 
${\rm \pi} /2$, i.e. 
$(v_1,v_2)^\perp =(-v_2,v_1)$. There are three different viscous moduli in our model: 
$\eta$ is a standard shear viscosity that arises from overcoming the magnetic attraction between nearby dipoles in relative motion, 
$\eta _R$ is the rotational viscosity that models friction between neighbouring rotating particles and 
$\eta _O$ is the odd viscosity (also known as the Hall viscosity) that gives rise to viscous forces acting transversely to a velocity gradient. As first noted in Avron (Reference Avron1998), the odd viscous stress is a non-dissipative term that is only permissible in 2-D fluids that do not obey local time-reversal symmetry. Its presence is intimately tied to the anti-symmetric driving stress.
The transverse motion of the surface phase droplet generates a shear stress on the bulk fluid below, and so we have for 
$\boldsymbol {u} = (u_1,u_2,u_3)$
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\sigma}=\left.\mu\frac{\partial(u_1,u_2)}{\partial z} \right|_{z=0}=\boldsymbol{f}, \quad \boldsymbol{x}\in \mathcal{D},\end{equation}
which can be interpreted as a boundary condition on the subphase. Outside of the surface phase domain 
$\mathcal {D}$ we have a simple stress-free boundary condition on the subphase, or
\begin{equation} \boldsymbol{f} = \boldsymbol{0}, \quad \boldsymbol{x}\notin \mathcal{D}.\end{equation}
Hence, (2.4) and (2.5) can be combined using a characteristic function 
$\chi$,
\begin{equation} \left.\mu\frac{\partial (u_1,u_2)}{\partial z} \right|_{z=0} = \chi({\mathcal{D}})(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\sigma}).\end{equation}
Here, it has been assumed that the normal stress of the bulk phase at 
$z=0$ is whatever it needs to be to maintain surface flatness. This could be achieved, for example, by having a high surface tension there.
The expression for 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {\sigma }$ has a remarkably simple form. Using the notation of the skew gradient, we note that the scalar vorticity can be written as 
$\omega = \nabla ^\perp \boldsymbol {\cdot } \boldsymbol {U}$. It then follows, using (2.2), that 
$\boldsymbol {\nabla }\boldsymbol {\cdot } (\omega {\boldsymbol{\mathsf{J}}}) = -\nabla ^\perp (\nabla ^\perp \boldsymbol {\cdot } \boldsymbol {U})=-{\rm \Delta} \boldsymbol {U}$. Similar manipulations yield the identities 
$\boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {\nabla }^\perp \boldsymbol {U})=\boldsymbol {0}$ and 
$\boldsymbol {\nabla }\boldsymbol {\cdot } (\boldsymbol {\nabla } \boldsymbol {U}^\perp ) = -\boldsymbol {\nabla } \omega$ for the divergence of the odd viscous tensor, which leads to the important consequence that the effect of the odd viscous stress in the bulk is simply to generate a ‘pressure field’ proportional to vorticity. In fact, it is convenient to define 
$\bar \eta = \eta + \eta _R$ and 
$P \to P + \eta _O\omega$ so that we may write
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\sigma}={-}\boldsymbol{\nabla} P + \bar\eta \boldsymbol{\nabla}^2\boldsymbol{U},\quad \boldsymbol{x}\in \mathcal{D}.\end{equation}
That is, the active phase is described by a 2-D Stokes equation with the viscosity and pressure redefined. The fact that (2.7) does not depend explicitly on 
$\varOmega$ or 
$\eta _O$ implies that the drive and Hall stress can appear only through boundary conditions. Thus, we have
\begin{equation} \chi(\mathcal{D})(-\boldsymbol{\nabla} P + \bar\eta \boldsymbol{\nabla}^2\boldsymbol{U}) =\left.\mu\frac{\partial(u_1,u_2)}{\partial z} \right|_{z=0},\end{equation}coupled to the Stokes equations (2.1a,b) for the bulk phase.
The surface Stokes equation requires additional, transverse boundary conditions. We impose a stress balance condition on the surface phase boundary 
$\partial \mathcal {D}$
\begin{equation} \boldsymbol{\sigma} \boldsymbol{\cdot} \hat{\boldsymbol{n}} = \gamma \kappa \hat{\boldsymbol{n}},\quad \boldsymbol{x} \in \partial \mathcal{D},\end{equation}
where 
$\gamma$ is the line tension, 
$\kappa$ is the local curvature of 
$\partial \mathcal {D}$ and 
$\hat {\boldsymbol {n}}$ is its inward facing normal vector (in the 
$z=0$ plane). In the local Frenet frame of 
$\partial \mathcal {D}$ where 
$\boldsymbol {U} = T\hat {\boldsymbol {t}} + N\hat {\boldsymbol {n}}$, we find that 
$(\boldsymbol {\nabla } \boldsymbol {U} +\boldsymbol {\nabla } \boldsymbol {U}^T)\boldsymbol {\cdot } \hat {\boldsymbol {n}} = 2 \boldsymbol {U}_s^\perp -\omega \hat {\boldsymbol {t}}, {\boldsymbol{\mathsf{J}}}\boldsymbol {\cdot } \hat {\boldsymbol {n}} = \hat {\boldsymbol {t}}$, and 
$(\boldsymbol {\nabla }^\perp \boldsymbol {U} + \boldsymbol {\nabla } \boldsymbol {U}^\perp )\boldsymbol {\cdot } \hat {\boldsymbol {n}} = 2\boldsymbol {U}_s - \omega \hat {\boldsymbol {n}}$ so that (2.9) can be written as
\begin{gather} -P + 2\eta (T_s - \kappa N) + 2\eta_O(N_s + \kappa T) |_{\partial \mathcal{D}} = \gamma \kappa |_{\partial \mathcal{D}}, \end{gather}
\begin{gather}-\bar\eta \omega - 2\eta (N_s +\kappa T) + 2\eta_O (T_s-\kappa N) |_{\partial \mathcal{D}}={-} 2\eta_R \varOmega |_{\partial\mathcal{D}}, \end{gather}
where the subscript 
$s$ denotes differentiation in the direction of 
$\hat {\boldsymbol {t}}$. This form of the boundary condition makes it clear that the role of the odd viscosity is to ‘complement’ the shear viscosity by producing a tangential boundary stress that depends on the local normal velocity and vice versa, thus coupling the flows in the two directions. We remark that the 
$\varOmega$ term in (2.11) provides the inhomogeneous forcing through which a non-trivial solution arises.
The dynamics of the active surface phase is a free-boundary problem for 
$\partial \mathcal {D}$. If its velocity is 
$\boldsymbol {V}$ then we evolve 
$\partial \mathcal {D}$ through the kinematic boundary condition 
$\lim _{x\to \partial \mathcal {D}} (\boldsymbol {V} - \boldsymbol {U}) \boldsymbol {\cdot } \hat {\boldsymbol {n}} = 0$. Here, we assume the limit is taken from the interior of the monolayer, so that no assumptions are made about continuity of velocity along 
$\partial \mathcal {D}$. Equations ((2.1a,b)–(2.3)), (2.6) and (2.9), together with the conditions 
$\hat {\boldsymbol {z}}\boldsymbol {\cdot } \boldsymbol {u}|_{z=0}=0$ and continuity of subphase and surface phase horizontal velocities at 
$z=0$, yield a complex, but complete, formulation for the determination of 
$\boldsymbol {U}$.
Some length and time scales: the many parameters of the model (Appendix A) give rise to a number of relevant length and time scales, a detailed analysis of which is provided in Appendix B. Here, we mention of few. One important length scale, independent of geometry and activity, is the Saffman–Delbrück length 
$\ell _{SD} = \eta /\mu$ (Saffman & Delbrück Reference Saffman and Delbrück1975). On length scales smaller than 
$\ell _{SD}$, momentum travels primarily in the plane of the surface phase, while for length scales larger than 
$\ell _{SD}$, momentum travels through the subphase as well. For the system at hand, 
$\ell _{SD}$ can be between 0.1 and 
$100\ {\rm \mu}{\rm m}$, depending on the viscosity of the subphase. For problems in which the characteristic size is variable, such as the edge tension-driven closure of a cavity punctured in the monolayer, the dynamics may take place in both regimes (Jia & Shelley Reference Jia and Shelley2022). A related length scale for surface phase droplets very close to the bottom wall is the penetration depth 
${\bar \delta }=(H{\bar \eta }/\mu )^{1/2}$. This is the length scale of edge currents driven by the rotational drive (Soni et al. Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019), and typically measures around 
$5\ {\rm \mu}{\rm m}$ (roughly 3 particle lengths). One obvious time scale is that of the rotational drive 
$\tau _1=\varOmega ^{-1}$, while another is the relaxational time scale 
$\tau _2=\eta R/\gamma$ driven by line tension; here, 
$R$ is the characteristic size of the monolayer. Both arise in the solutions constructed herein, even though the time scale for rotation is approximately one hundredth of that for relaxation (that is, for the experiments of Soni et al. Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019).
2.1. A Green's function formulation for the free-boundary problem
The infinite horizontal extent of the domain makes the problem amenable to classical Fourier transform methods. Masoud & Shelley (Reference Masoud and Shelley2014) showed that in this geometry, 2-D Fourier transforming (2.4) in 
$x$ and 
$y$ and applying boundary conditions yields the relation
\begin{equation} \hat{\boldsymbol{f}} = \mu k[A(kH) {\boldsymbol{\mathsf{I}}} + B(kH)\hat{\boldsymbol{k}} \hat{\boldsymbol{k}} ]\hat{\boldsymbol{U}},\end{equation}with
\begin{equation} A(\alpha) = \coth \alpha, \end{equation}and
\begin{equation} B(\alpha) = \frac{\alpha^2 \coth \alpha -2\alpha + \sinh \alpha \cosh \alpha}{\sinh^2 \alpha - \alpha^2}, \end{equation}
where 
$\boldsymbol {k}$ is the 2-D wavenumber with magnitude 
$k$ and unit vector 
$\hat {\boldsymbol {k}}$. If the surface velocity further satisfies 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {U} = 0$ everywhere on 
$z=0$, as has been assumed in other works (e.g. Stone & McConnell Reference Stone and McConnell1995; Lubensky & Goldstein Reference Lubensky and Goldstein1996; Alexander et al. Reference Alexander, Bernoff, Mann, Mann and Wintersmith2007), (2.12) then simplifies dramatically to
\begin{equation} \hat{\boldsymbol{f}} = \mu k A(kH) \hat{\boldsymbol{U}} = \mu k \coth (kH)\hat{\boldsymbol{U}}. \end{equation}
Such an assumption confers the advantage of making 
$w=0$ and 
$p=0$ in the bulk fluid, which in turn makes some problems analytically tractable. However, we emphasize that this simplification does not follow from assuming that 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {U} = 0$ in 
$\mathcal {D}$ alone, even in the axisymmetric case, so we will not make that assumption here.
Equation (2.12) can be rewritten as
\begin{equation} \hat{\boldsymbol{U}} =\frac{1}{\mu k} \left[ \frac{1}{A(kH)} {\boldsymbol{\mathsf{I}}} - \frac{B(kH)}{A(kH)[A(kH) + B(kH)]} \hat{\boldsymbol{k}} \hat{{\boldsymbol{k}}} \right] \hat{\boldsymbol{f}}, \end{equation}
which can be interpreted as a statement that 
$\boldsymbol {U}$ is a convolution of 
$\boldsymbol {f}$ against some tensorial Green's function 
${\boldsymbol{\mathsf{G}}}_{{\boldsymbol{\mathsf{H}}}}$, expressible as an inverse 2-D Fourier transform
\begin{equation} {\boldsymbol{\mathsf{G}}}_{{\boldsymbol{\mathsf{H}}}} = \frac{1}{(2{\rm \pi})^2} \int_0^{2{\rm \pi}}\, \mathrm{d}\phi \int_0^\infty\, \mathrm{d}k\, k \frac{1}{\mu k} \left[ \frac{1}{A(kH)} {\boldsymbol{\mathsf{I}}} - \frac{B(kH)}{A(kH)[A(kH) + B(kH)]}\hat{\boldsymbol{k}} \hat{\boldsymbol{k}} \right] \, {\rm e}^{\mathrm{i} kr\cos(\theta-\phi)}, \end{equation}
where 
$r$ and 
$\theta$ are the polar coordinates of real space and 
$k$ and 
$\phi$ are the polar coordinates of Fourier space. We define 
$\beta = \theta - \phi$, change variables, and apply periodicity in 
$\beta$ to obtain
\begin{equation} {\boldsymbol{\mathsf{G}}}_{{\boldsymbol{\mathsf{H}}}} = \frac{1}{(2{\rm \pi})^2} \int_0^{2{\rm \pi}}\, \mathrm{d}\beta \int_0^\infty\, \mathrm{d}k \frac{1}{\mu}\left[\frac{1}{A(kH)}{\boldsymbol{\mathsf{I}}} - \frac{B(kH)}{A(kH)[A(kH) + B(kH)]} \hat{\boldsymbol{Z}}\hat{\boldsymbol{Z}}\right] \, {\rm e}^{\mathrm{i} kr\cos \beta}, \end{equation}
where 
$\hat {\boldsymbol {Z}} = \cos \beta \hat {\boldsymbol {x}} -\sin \beta \hat {\boldsymbol {x}}^\perp$ comes from expressing 
$\hat {\boldsymbol {k}}$ in the 
$\{\hat {\boldsymbol {x}},\hat {\boldsymbol {x}}^\perp \}$ basis. Performing a Jacobi–Anger expansion (formula 8.511.1 of Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007) to evaluate the integrals over 
$\beta$, we arrive at
\begin{equation} {\boldsymbol{\mathsf{G}}}_{{\boldsymbol{\mathsf{H}}}} = \frac{1}{2{\rm \pi} \mu} \int_0^\infty \mathrm{d}k \left\{\frac{J_0(kr)}{A(kH)}{\boldsymbol{\mathsf{I}}}- \frac{B(kH)}{A(kH)[A(kH)+B(kH)]} \left[{J_1'(kr)}\hat{\boldsymbol{x}}\hat{\boldsymbol{x}} + \frac{J_1(kr)}{kr} \hat{\boldsymbol{x}}^\perp \hat{\boldsymbol{x}}^\perp\right]\right\}, \end{equation}
where 
$J_0$ and 
$J_1$ are Bessel functions of the first kind. By the convolution theorem and (2.6),
\begin{equation} \boldsymbol{U}={\boldsymbol{\mathsf{G}}}_{{\boldsymbol{\mathsf{H}}}}*[\chi(\mathcal{D})\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\sigma}] ={\boldsymbol{\mathsf{G}}}_{{\boldsymbol{\mathsf{H}}}}*[\chi(\mathcal{D})(-\boldsymbol{\nabla} P + \bar\eta \boldsymbol{\nabla}^2 \boldsymbol{U})]. \end{equation}
This equation holds for all points in the 
$z=0$ plane and, if evaluated inside of 
$\mathcal {D}$, generates an integral relation for the unknown 
$\boldsymbol {U}$. No standard method exists for its solution. The first expression involving 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\sigma }$ serves only to emphasize its generality – other surface stress tensors could be considered. For example, if 
$\boldsymbol {\sigma }$ is Newtonian, (2.20) could model a classical incompressible Langmuir film (see Alexander et al. Reference Alexander, Bernoff, Mann, Mann and Zou2006). Equation (2.20) is likewise applicable to other materials such as active nematic films (e.g. Gao et al. Reference Gao, Blackwell, Glaser, Betterton and Shelley2015).
Thus, (2.20), together with the incompressibility condition (2.2), the stress boundary condition (2.9) and the kinematic boundary condition, 
$\lim _{x\to \partial \mathcal {D}} (\boldsymbol {V} - \boldsymbol {U}) \boldsymbol {\cdot } \hat {\boldsymbol {n}}$, where the limit is taken from the interior of 𝒟, completely specify the initial value problem for 
$\partial \mathcal {D}$. Note (2.20) is easily extended to a domain composed of multiple monolayers.
There are two limits where 
${\boldsymbol{\mathsf{G}}}_{{\boldsymbol{\mathsf{H}}}}$ takes a particularly simple closed form.
(i) In the limit of infinite 
$H, A(kH)$ and 
$B(kH)$ both approach unity with zero slope, and 
${\boldsymbol{\mathsf{G}}}_{{\boldsymbol{\mathsf{H}}}}$ approaches
\begin{equation} {\boldsymbol{\mathsf{G}}}_{\infty} = \frac{1}{2{\rm \pi} \mu r} \left({\boldsymbol{\mathsf{I}}} - \frac{1}{2} {\hat{\boldsymbol{x}}^\perp\hat{\boldsymbol{x}}^\perp}\right). \end{equation}
Using the fact that 
${\boldsymbol{\mathsf{I}}} =\hat {\boldsymbol {x}} \hat {\boldsymbol {x}}+\hat {\boldsymbol {x}}^\perp \hat {\boldsymbol {x}}^\perp$, we can rewrite this expression as
\begin{equation} {\boldsymbol{\mathsf{G}}}_{\infty} = \frac{{\boldsymbol{\mathsf{I}}} + \hat{\boldsymbol{x}}\hat{\boldsymbol{x}}}{4{\rm \pi} \mu r}. \end{equation}That is, in the case of infinite subphase depth, the Green's function is proportional to the classical Stokeslet.
(ii) In the limit of small depth, 
$H\to 0, A(kH) \to (kH)^{-1}$ and 
$B(kH) \to (kH/3)^{-1}$. In this case, we interpret the formal limit
\begin{equation} {\boldsymbol{\mathsf{G}}}_{\boldsymbol 0}=\frac{H}{2{\rm \pi} \mu} \int_0^\infty\, \mathrm{d}k \left\{kJ_0(kr){\boldsymbol{\mathsf{I}}} - \frac{3}{4}k \left[J_1'(kr) \hat{\boldsymbol{x}}\hat{\boldsymbol{x}} + \frac{J_1(kr)}{kr}\hat{\boldsymbol{x}}^\perp \hat{\boldsymbol{x}}^\perp\right]\right\}, \end{equation}using generalized functions. Starting with the orthogonality relation (formula 6.512.8 of Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007)
\begin{equation} \int_0^\infty \, \mathrm{d}k\,kJ_0(kr)J_0(kr') = \frac{\delta(r-r')}{r}, \end{equation}
and letting 
$r'\to 0$, we find that
\begin{equation} \int_0^\infty \, \mathrm{d}k\, kJ_0(kr) = \frac{\delta(r)}{r}. \end{equation}This can be combined with the convergent integral
\begin{equation} \int_0^\infty \, \mathrm{d}k\, k\frac{J_1(kr)}{kr}=\frac{1}{r^2}, \end{equation}to show
\begin{equation} \int_0^\infty \, \mathrm{d}k\, kJ_1'(kr) = \int_0^\infty \, \mathrm{d}k \left[J_0(kr)-\frac{J_1(kr)}{kr}\right] = \frac{\delta(r)}{r} - \frac{1}{r^2}. \end{equation}Thus,
\begin{equation} {\boldsymbol{\mathsf{G}}}_{\boldsymbol 0}= \frac{H}{2{\rm \pi} \mu} \left\{\frac{\delta(r)}{r} {\boldsymbol{\mathsf{I}}} - \frac{3}{4} \left[\left(\frac{\delta(r)}{r} - \frac{1}{r^2}\right)\hat{\boldsymbol{x}}\hat{\boldsymbol{x}} + \frac{1}{r^2} \hat{\boldsymbol{x}}^\perp\hat{\boldsymbol{x}}^\perp\right]\right\}.\end{equation}
We can rewrite the quantity in square brackets by observing that its Fourier transform upon convolution with 
$\boldsymbol {f}$ is
\begin{equation} 2{\rm \pi}(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}) \boldsymbol{\cdot} \hat{\boldsymbol{f}} ={-}2{\rm \pi}\mathrm{i}\frac{\hat{\boldsymbol{k}}}{k} (\mathrm{i} {\boldsymbol{k}}\boldsymbol{\cdot} \hat{\boldsymbol{f}}).\end{equation}
Since 
$-1/k^2$ is the Fourier transform of the fundamental solution to the Laplace equation, 
$-\mathrm {i}\hat {\boldsymbol {k}}/k$ is the Fourier transform of its gradient, namely 
$\hat {\boldsymbol {x}}/(2{\rm \pi} r)$ in two dimensions. Applying the convolution theorem to the quantity on the right-hand side of (2.29) thus establishes the relation
\begin{equation} \left[\left(\frac{\delta(r)}{r} - \frac{1}{r^2}\right)\hat{\boldsymbol{x}}\hat{\boldsymbol{x}} + \frac{1}{r^2} \hat{\boldsymbol{x}}^\perp\hat{\boldsymbol{x}}^\perp \right]*\boldsymbol{f}= \frac{\hat{\boldsymbol{x}}}{r} * (\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{f}}), \end{equation}so that
\begin{equation} \boldsymbol{U} = {\boldsymbol{\mathsf{G}}}_{\boldsymbol 0} * \boldsymbol{f}= \frac{H}{2{\rm \pi} \mu} \left[\frac{\delta(r)}{r}{\boldsymbol{\mathsf{I}}} * \boldsymbol{f} - \frac{3}{4} \frac{\hat{\boldsymbol{x}}}{r}* (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{f})\right]=\frac{H}{\mu}\left[ \boldsymbol{f} - \frac{3}{4} \frac{\hat{\boldsymbol{x}}}{2{\rm \pi} r} * (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{f})\right]. \end{equation}
Taking a divergence and using the fact that 
$\boldsymbol {\nabla } \boldsymbol {\cdot } (\hat {\boldsymbol {x}}/r) = 2{\rm \pi} \delta _2(\boldsymbol {x})$, where 
$\delta _2(\boldsymbol {x})$ is the 2-D Dirac delta distribution, we arrive at the simple relation
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{U} = \frac{H}{4\mu}(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{f}). \end{equation}
Note that this result is also directly obtainable by letting 
$H\to 0$ in (2.12), which gives
\begin{equation} \hat{\boldsymbol{f}} = \frac{\mu}{H} ({\boldsymbol{\mathsf{I}}} + 3\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}) \hat{\boldsymbol{U}}, \end{equation}
multiplying by 
$\mathrm {i}\hat {\boldsymbol {k}}$, and taking an inverse Fourier transform.
Equation (2.32), in addition to the assumptions that 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {U}= 0$ in 
$\mathcal {D}$ and 
$\boldsymbol {f} = 0$ in 
$\mathcal {D}^C$, is sufficient to guarantee that 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {U} = \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {f}= 0$ on the entire surface so that ultimately
\begin{equation} {\boldsymbol{\mathsf{G}}}_{\boldsymbol 0}=\frac{\delta_2(\boldsymbol{x}){\boldsymbol{\mathsf{I}}}}{ \varGamma }, \end{equation}
where 
$\varGamma = \mu /H$. Equation (2.20) then becomes
\begin{equation} \chi(\mathcal{D})(-\boldsymbol{\nabla} P + \bar\eta \boldsymbol{\nabla}^2\boldsymbol{U}) =\varGamma\boldsymbol{U}. \end{equation}
That is, at leading order, the forcing is not only linear but local in 
$\boldsymbol {U}$, and the equation of motion reduces to a Brinkman equation. Equation (2.35) can also be derived by taking a thin film lubrication approximation in the subphase (Barentin et al. Reference Barentin, Ybert, di Meglio and Joanny1999; Elfring, Leal & Squires Reference Elfring, Leal and Squires2016). The Brinkman equation arises frequently in the modelling of porous materials and especially Newtonian mono- and bilayers (Evans & Sackmann Reference Evans and Sackmann1988). It is readily seen that in this limit, the pressure inside the monolayer is harmonic, and the fluid velocity outside of the monolayer is zero. Expanding (2.15) further in small 
$H$ yields
\begin{equation} \hat{\boldsymbol{f}} = \mu k \left(\frac{1}{kH} + \frac{kH}{3} + \cdots\right)\hat{\boldsymbol{U}} = \varGamma \hat{\boldsymbol{U}} + \frac{\mu H}{3} k^2 \hat{\boldsymbol{U}} + \cdots . \end{equation}
Thus, in real space the 
$O(H)$ correction enters as a term proportional to 
$-{\rm \Delta} \boldsymbol {U}$, which has the effect of slightly increasing the monolayer viscosity (Lubensky & Goldstein Reference Lubensky and Goldstein1996).
In summary, the small 
$H$ limit of the low friction case is the high friction case with substrate drag coefficient 
$\varGamma$.
3. The axisymmetric case
We proceed to specialize the integral equation derived above to axisymmetric domains of either disks or annuli. It is advantageous to consider the radial and azimuthal equations separately. To this end, let 
$\boldsymbol {U} = U(r) \hat {\boldsymbol {x}} + V(r)\hat {\boldsymbol {x}}^\perp$ and 
$\boldsymbol {f} = f(r)\hat {\boldsymbol {x}} + g(r) \hat {\boldsymbol {x}}^\perp$. Substituting into (2.8), we obtain the momentum equations
\begin{gather} \chi(\mathcal{D}) \left(-\frac{\mathrm{d} P}{\mathrm{d} r} + \bar\eta \mathcal{L}[U] \right)= f, \end{gather}
\begin{gather}\chi(\mathcal{D}) \bar\eta{\mathcal{L}}[V] = g, \end{gather}where we have defined the differential operator
\begin{equation} \mathcal{L} = \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial }{\partial r} - \frac{1}{r^2}. \end{equation}For axisymmetric quantities, 2-D Fourier transforms reduce to Hankel transforms
\begin{gather} \hat{\boldsymbol{U}} ={-}2{\rm \pi} \mathrm{i}\hat{\boldsymbol{k}} \int_0^\infty \, \mathrm{d}r\, r U(r) J_1(kr) - 2{\rm \pi} \mathrm{i} \hat{\boldsymbol{k}}^\perp\int_0^\infty \, \mathrm{d}r\, r V(r) J_1(kr), \end{gather}
\begin{gather}\hat{\boldsymbol{f}} ={-}2{\rm \pi} \mathrm{i}\hat{\boldsymbol{k}} \int_0^\infty \, \mathrm{d}r\, r f(r) J_1(kr) - 2{\rm \pi} \mathrm{i} \hat{\boldsymbol{k}}^\perp\int_0^\infty \, \mathrm{d}r\, r g(r) J_1(kr). \end{gather}Substituting these expressions into (2.16) yields the relations
\begin{gather} \int_0^\infty \, \mathrm{d}r\, r U(r) J_1(kr) = \frac{1}{\mu k [A(kH) + B(kH)]} \int_0^\infty \, \mathrm{d}r\, rf(r) J_1(kr), \end{gather}
\begin{gather}\int_0^\infty \, \mathrm{d}r\, rV(r) J_1(kr) = \frac{1}{\mu k A(kH)} \int_0^\infty \, \mathrm{d}r\, rg(r) J_1(kr). \end{gather}
These equations show that the azimuthal and radial dynamics are decoupled in the bulk, but as we will see, this is not always the case at the boundary 
$\partial \mathcal {D}$, where the two interact through the odd viscosity. By taking another Hankel transform to solve for the velocity components, we obtain the axisymmetric analogues of (2.20)
\begin{gather} U(r) = \frac{1}{\mu}\int_0^\infty \, \mathrm{d}k \frac{J_1(kr)}{A(kH)+B(kH)} \int \, \mathrm{d}r'\, r'f(r') J_1(kr')= \frac{1}{\mu} \int_{R_i}^{R_o} \, \mathrm{d}r'\, r' M(r,r') f(r'), \end{gather}
\begin{gather}V(r) = \frac{1}{\mu}\int_0^\infty \, \mathrm{d}k \frac{J_1(kr)}{A(kH)} \int \, \mathrm{d}r'\, r'g(r') J_1(kr')=\frac{1}{\mu} \int_{R_i}^{R_o} \, \mathrm{d}r'\, r' L(r,r') g(r'), \end{gather}
where 
$R_i< R_o, R_i$ can be zero, and 
$R_o$ can be positive infinity. We have also defined the kernels
\begin{gather} M(r,r') = \int_0^\infty \, \mathrm{d}k \frac{J_1(kr)J_1(kr')}{A(kH) + B(kH)}, \end{gather}
\begin{gather}L(r,r') = \int_0^\infty \, \mathrm{d}k \frac{J_1(kr)J_1(kr')}{A(kH)}. \end{gather} In the limit 
$H\to \infty$, we obtain closed form expressions for 
$M$ and 
$L$, which we denote with a bar, from Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007) formulas 6.512.1, 8.126.3 and 8.126.4
\begin{align} \bar{L}(r,r')= 2\bar{M}(r,r') &= \int_0^\infty \, \mathrm{d}k\, J_1(kr)J_1(kr') \end{align}
\begin{align} &= \left\{\begin{array}{@{}cc} \dfrac{2}{{\rm \pi} r}({-}E[r^2/r'^2] + K[r^2/r'^2]) & \text{if }r< r'\\ \dfrac{2}{{\rm \pi} r'}({-}E[r'^2/r^2] + K[r'^2/r^2]) & \text{if }r'< r \end{array}\right. \end{align}
\begin{align} &={-}\frac{1}{{\rm \pi} rr'(r+r')}[(r+r')^2E[\xi] -(r^2+r'^2)K[\xi]], \end{align}
where 
$\xi = 4rr'/(r+r')^2$, and 
$K$ and 
$E$ are complete elliptic integrals of the first and second kind, respectively. The reader is cautioned that our notational convention for elliptic integrals differs from that of Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007) by a square root in the argument. Since 
$K$ has a logarithmic singularity when its argument approaches unity, these expressions illustrate that 
$L(r,r')$ and 
$M(r,r')$ are both logarithmically singular when 
$r\sim r'$, and it becomes useful to isolate the most singular behaviour. We write
\begin{equation} \bar{L}(r,r') ={-}\frac{\log|r-r'|}{{\rm \pi} r'} + \tilde{L}(r,r'), \end{equation}
where 
$\tilde {L}$ has a removable singularity at 
$r=r'$, with an analogous expansion for 
$\bar {M}$.
In the opposite limit of 
$H\to 0, A(kH) \sim (kH)^{-1}$ and 
$L$ is proportional to a Dirac delta
\begin{equation} \int_0^\infty \, \mathrm{d}k\, kH\, J_1(kr)J_1(kr') = H \frac{\delta(r-r')}{r'}, \end{equation}
using orthogonality properties of 
$J_1$ (formula 6.512.8 from Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007). The convolution integral thus reduces to a Brinkman ordinary differential equation whose solution is discussed in §§ 4.3 and 5.3. For the radial direction, we similarly find 
$M\to 4H\delta (r-r')/r'$.
As noted by Yan & Sloan (Reference Yan and Sloan1988), integral equations with logarithmically singular kernels generally have unique solutions that diverge like an inverse square root of the distance to the boundary. For our domains, we can therefore expect the surface shear stresses of the bulk fluid 
$f(r)$ and 
$g(r)$ to take the form
\begin{equation} f(r) = \frac{\tilde{f}(r)}{\sqrt{r-R_i}\sqrt{R_o-r}}\chi(R_i< r< R_o), \end{equation}and
\begin{equation} g(r) = \frac{\tilde{g}(r)}{\sqrt{r-R_i}\sqrt{R_o-r}}\chi(R_i< r< R_o), \end{equation}
for smooth functions 
$\tilde {f}$ and 
$\tilde {g}$ that do not vanish at 
$r=R_i$ and 
$r=R_o$. As such, 
$V''(r)$, as well as 
$P'(r)$ in the case of an annulus, similarly diverge like an inverse square root of the distance to the boundaries when the limit is taken from the domain interior. (The exterior velocity fields will also exhibit singularities in their derivatives as the boundary is approached; see § 5.2 for one analytical example.) Nonetheless, in spite of these divergences, both 
$V(r)$ and 
$V'(r)$ remain bounded. This phenomenon is consistent with other axisymmetric systems with vanishingly thin domains immersed in a continuous 3-D medium, such as the rotating solid disk submerged in a Stokes fluid analysed by Jeffery (Reference Jeffery1915) or the penny-shaped crack in a 3-D elastic medium analysed by Sneddon (Reference Sneddon1946). For the rotating disk, Sherwood (Reference Sherwood2013) has found that relaxing the no-slip condition by allowing for a finite slip length can regularize this singularity, but such a condition is not present in our model. A consequence of the divergence at the boundary is that even the problem of the linearly perturbed disc requires a great deal of subtlety, as linearization requires further differentiation of the boundary terms. Linear stability analysis of a related system was previously treated by Stone & McConnell (Reference Stone and McConnell1995), although they included a simplifying global incompressibility assumption as well as a non-zero surface viscosity everywhere, which circumvents the issue at hand. A linear stability analysis for the 
$H\to 0$ case was previously completed in Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019) since the integral kernel is delta singular rather than logarithmically singular and velocity gradients do not diverge in this case.
In the following sections, we describe several ways to derive and numerically solve the singular integral equations (3.8) and (3.9) for two important axisymmetric geometries: the disc and the annulus.
4. A disc-shaped domain
We begin with the simplest non-trivial axisymmetric geometry and take the domain 
$\mathcal {D}$ to be the disc of radius 
$R$ centred at the origin. The incompressibility of the monolayer disc centred at the origin automatically implies that there is no radial component to the axisymmetric flow; consequently, the entire surface flow field is incompressible and (2.15) holds. The pressure gradient also vanishes. Hence, we can formulate the problem of finding the flow field on the surface as a scalar mixed boundary value problem using (2.6). Let 
$v(r,z)$ be the azimuthal component of 
$\boldsymbol {u}(r,z)$ and let 
$V(r) = v(r,z=0)$. The azimuthal momentum equation can be written as
\begin{equation} \left. \begin{array}{ll} \displaystyle\left. \mu \dfrac{\partial v}{\partial z}\right|_{z=0}=\bar\eta \left(\dfrac{\mathrm{d}^2 V}{\mathrm{d} r^2} + \dfrac{1}{r} \dfrac{\mathrm{d} V}{\mathrm{d} r} - \dfrac{V}{r^2} \right) & \mbox{for}\; r < R\\ \displaystyle\left. \mu \dfrac{\partial v}{\partial z}\right|_{z=0}= 0 & \mbox{for}\; r>R \end{array}\right\},\end{equation}or, more compactly,
\begin{equation} g(r) = \bar\eta \mathcal{L}[V(r)]\chi(r< R), \end{equation}
where the 2-D vector Laplacian operator 
$\mathcal {L}$ was defined in (3.3). If the 2-D surface pressure outside of the monolayer is taken to be zero, the stress boundary condition (2.9) in this geometry reduces to the Robin boundary condition
\begin{equation} V'(R^-) - \frac{\eta-\eta_R}{\bar\eta R} V(R) = \frac{2\eta_R\varOmega}{\bar\eta}, \end{equation}in the azimuthal direction; it is unnecessary to analyse the radial direction as it ultimately simply sets the pressure difference. Note that in this situation, the absence of a radial velocity means the odd stress only produces a transverse stress which is balanced by the pressure. Also note this boundary condition is what introduces an inhomogeneous forcing into the problem. We demonstrate two techniques to solve for the droplet flow field: a direct inversion of the singular integral equation in convolution form (3.9) and an equivalent formulation as a dual integral equation for the Hankel transformed velocity field that admits a semi-analytic solution.
4.1. Solution via Green's function
We begin by noting that the solution to (4.1) for 
$r< R$ can be written as
\begin{equation} V(r) = V_p(r) + V_h(r), \end{equation}where the particular solution is of the form
\begin{equation} V_p(r) = \int \, \mathrm{d}r'\, G(r,r') g(r'), \end{equation}
with the Green's function 
$G$ of the operator 
$\bar \eta \mathcal {L}$ satisfying
\begin{equation} \bar\eta \left( \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial }{ \partial r} - \frac{1}{r^2}\right) {G(r,r')} = \delta(r-r'), \end{equation}and the homogeneous solution is of the form
\begin{equation} V_h(r) = A_hr + \frac{B_h}{r}, \end{equation}
with constants 
$A_h$ and 
$B_h$ that are chosen to satisfy the boundary condition (4.3) as well as the condition 
$V(0) = 0$ that is a consequence of axisymmetry.
The Green's function 
$G(r,r')$ is easily calculated via any of a multitude of methods (Hankel transforms for instance). We find
\begin{equation} G(r,r') ={-}\frac{1}{\bar\eta} \int_0^\infty \, \mathrm{d}k \frac{J_1(kr)J_1(kr')}{k} ={-}\frac{1}{\bar\eta} \left[\frac{r}{2} \chi({r< r'}) + \frac{r'^2}{2r} \chi({r'< r})\right], \end{equation}so that
\begin{equation} V_p(r) = \int_0^R \, \mathrm{d}r'\, G(r,r') g(r') ={-} \frac{1}{2\bar\eta r} \int_0^r \, \mathrm{d}r'\, (r')^2 g(r') - \frac{r}{2\bar\eta} \int_r^R \, \mathrm{d}r'\, g(r'). \end{equation}
It remains to solve for 
$A_h$ and 
$B_h$. The condition 
$V(0)=0$ implies that 
$B_h=0$. The azimuthal velocity 
$V= V_p + A_hr$ is continuous and piecewise smooth, so upon substitution into (4.3), we obtain
\begin{equation} A_h = \varOmega + \frac{\eta-\eta_R}{2\eta_R R} V_p(R) - \frac{\bar\eta}{2\eta_R} V_p'(R^-), \end{equation}
where 
$R^-$ denotes the limit as 
$r$ approaches 
$R$ from the left. Using (4.9), this can be written as
\begin{equation} A_h = \varOmega - \frac{\eta}{2\eta_R \bar\eta R^2} \int_0^R \, \mathrm{d}r'\, (r')^2 g(r'). \end{equation}
Finally, substituting everything into (3.9) yields the following integral equation for 
$g$ in the interval 
$0< r< R$:
\begin{equation} \int_0^R \, \mathrm{d}r'\,G(r,r') g(r') + r \left[\varOmega - \frac{\eta}{2\eta_R \bar\eta R^2} \int_0^R \, \mathrm{d}r'\,(r')^2 g(r')\right] = \frac{1}{\mu} \int_0^R \, \mathrm{d}r'\, r' L(r,r') g(r'). \end{equation}
The discretization of the singular kernel 
$L(r,r')$ is handled as follows. First, we expand in 
$H$ so that
\begin{equation} L(r,r') = \bar{L}(r,r') + L^*(r,r'), \end{equation}where
\begin{equation} L^*(r,r') = \int_0^\infty \, \mathrm{d}k\, J_1(kr)J_1(kr') (\tanh kH - 1). \end{equation}
Because the integrand of 
$L^*$ decays to zero exponentially, this integral is well approximated by taking a finite cutoff, say by replacing the upper limit of 
$\infty$ by 
$10/H$, and numerically integrating. Next, we focus on the diagonal terms by writing the integral as
\begin{equation} \frac{1}{\mu} \int_0^R \, \mathrm{d}r'\, r' \bar{L}(r,r') g(r') = \frac{1}{\mu}\int_0^R \, \mathrm{d}r'\, r'\bar{L}(r,r') [g(r')-g(r)] + \frac{1}{\mu}g(r)\int_0^R \, \mathrm{d}r'\, r' \bar{L}(r,r'). \end{equation}
The first term on the right-hand side vanishes as 
$r$ approaches 
$r'$ if 
$g$ is smooth enough. Formula 6.561.13 from Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007) allows us to evaluate the second integral on the right-hand side, which we denote 
$L_d$, in terms of hypergeometric functions 
${}_pF_q$
\begin{align} L_d(r) &= \int_0^R \, \mathrm{d}r'\, r'\bar{L}(r,r') \nonumber\\ &= \int_0^\infty \, \mathrm{d}k\, J_1(kr) \int_0^R\mathrm{d}r'\, r'J_1(kr') \nonumber\\ &=\frac{r}{32}\left({-}8 - \frac{3r^2}{R^2} {}_4F_3\left[\left\{1,1,\frac{3}{2},\frac{5}{2}\right\}, \{2,2,3\},\frac{r^2}{R^2}\right]+ 16 \log \frac{4R}{r}\right). \end{align}
Once 
$g$ is known, (4.9) and (4.11) can be used to calculate 
$V$ inside the monolayer. When 
$r>R$, the integral equation is non-singular away from 
$r=R$ and the same 
$g$ can be substituted into (3.9) to calculate 
$V$ directly.
Figure 2(a) shows 
$V(r)$ calculated for various subphase thicknesses 
$H$. Within 
$\mathcal {D}, V(r)$ is upwardly convex and at larger values of 
$H$ (well described by the infinite 
$H$ case) shows nearly solid-body rotation near the droplet centre with a delocalized, faster edge current. This agrees qualitatively with the observations of Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019). We find that 
$V$ is continuous and smooth everywhere except at the interface 
$r=R$. As the limit is taken from either side of the interface, 
$V'$ is found to be finite valued but 
$V''$ diverges like 
$|r-R|^{-1/2}$, as previously discussed. The maximal 
$V$ is always found at 
$r=R$; outside of the disc, 
$V(r)$ decays exponentially with increasing 
$r$ if 
$H$ is finite. (For infinite 
$H$, the rate of decay is an inverse quadratic.) Decreasing 
$H$ has the effect of generally reducing the motion both in the monolayer bulk and in the exterior. In the limit of 
$H\ll R$, the motion becomes largely localized to a boundary layer at 
$r=R$ whose thickness scales with the penetration depth 
$\bar {\delta } = \sqrt {\bar \eta H/\mu }$, and 
$V$ is well-approximated by (4.44). This agrees with the analytical prediction of an edge current in the high friction case, as discussed in § 4.4 and in Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019).
Figure 2. (a) The surface azimuthal velocity 
$V$ across an active disc-shaped monolayer of radius 
$R$ as subphase depth 
$H$ is varied. As 
$H\to 0$, a boundary layer of width 
$\bar {\delta } = \sqrt {\bar \eta H/\mu }$ becomes visible and 
$V$ is well approximated by (4.44). Parameters: 
$\eta _R/\eta =1.875, \mu R/\eta = 30$. Here, 
$R$ is large compared with the Saffman–Delbrück length 
$\ell _{SD}=\eta /\mu$, and 
$V$ scales like 
$\ell _{SD}\varOmega$ for fixed 
$H$. Since the radial velocity 
$U=0$ for the disc, the flow fields are independent of line tension and odd viscosity. (b) For the same parameters, with 
$H/R=1$, the azimuthal velocity at different sublayer depths.
The above results can be framed in terms of the Saffman–Delbrück length, 
$\ell _{SD} = \eta /\mu$. Non-dimensionalizing the monolayer momentum equation in the simplest case where 
$\eta = \eta _R$ and 
$H\to \infty$, we obtain
\begin{equation} \bar{\beta} \mathcal{L}[V] = \left.\frac{\partial v}{\partial z}\right|_{z=0}, \end{equation}
where all lengths have been scaled by 
$R$ and 
$\bar {\beta } = 2\ell _{SD}/R$ is a dimensionless parameter formed from a ratio of the only remaining length scales. If 
$\bar {\beta } \gg 1$, then within the region where 
$r< R\ll \ell _{{SD}}$, momentum is dissipated primarily through the monolayer, and 
$\mathcal {L}[V]$ is necessarily small. Taken with the boundary condition (4.3), this implies 
$V \approx \varOmega r$ inside of the monolayer. Thus, we recover the rotating rigid disc on a half-space of Stokes fluid with vanishing surface viscosity analysed by Goodrich (Reference Goodrich1969) (see also the solution of Jeffery (Reference Jeffery1915) to the problem of a rotating solid disc in an infinite Stokes medium). In this case, the bulk shear stress at the surface is
\begin{equation} \left.\mu \frac{\partial v}{\partial z}\right|_{z=0} = \frac{4\mu\varOmega}{\rm \pi}\frac{r}{\sqrt{R^2-r^2}}\chi(r< R), \end{equation}and the surface velocity field is given by
\begin{equation} V(r) = \left\{\begin{array}{ll} \varOmega r & 0< r< R\\ \dfrac{2\varOmega}{{\rm \pi} r}\left[r^2\sin^{{-}1}\left(\dfrac{R}{r}\right) - R\sqrt{r^2 -R^2}\right] & r>R \end{array}\right. , \end{equation}so that the aforementioned regularity properties are all explicitly observable.
In the opposite limit of 
$\bar {\beta } \ll 1$, momentum is dissipated primarily through the bulk subphase, and the dimensionless velocity profile is nonlinear. Now, the dimensionless 
$(\partial v/\partial z)|_{z=0}$ scales like 
$\bar \beta$ so that the dimensional velocity scales like 
$\ell _{SD}\varOmega$ in this regime. This is the 
$\bar \beta$ regime depicted in figure 2.
Once the surface velocity has been found, the subphase velocity and pressure can be calculated by specializing the formulas that appear in Masoud & Shelley (Reference Masoud and Shelley2014). Define the Hankel transformed velocity component
\begin{equation} \mathcal{H}[U(r')](k) = \int_0^\infty \, \mathrm{d}r'\, r'J_1(kr')U(r'), \end{equation}
and similarly for 
$V(r')$. For axisymmetric geometries
\begin{equation} v(r,z) = \int_0^\infty \, \mathrm{d}k\, k J_1(kr)\mathcal{H}[V(r')] \frac{\sinh k(z+H)}{\sinh kH}, \end{equation}
\begin{align} &u(r,z) = \int_0^\infty \, \mathrm{d}k\, k J_1(kr)\mathcal{H}[U(r')]\nonumber\\ &\quad \times \left\{\dfrac{\sinh k(z+H)}{\sinh kH} + \dfrac{kz\sinh kH \coth k(z+H) + kH [(kH \coth kH-1)\sinh kz - kz\cosh kz]}{\sinh^2 kH - k^2H^2}\right\}, \end{align}
\begin{gather} w(r,z) = \int_0^\infty \, \mathrm{d}k \, kJ_0(kr) \mathcal{H}[U(r')] \frac{k^2H (z+H) \sinh kz - kz \sinh kH \sinh k(z+H)}{\sinh^2 kH -k^2 H^2}, \end{gather}
\begin{gather}p(r,z) = 2\mu \int_0^\infty \, \mathrm{d}k\, k^2 J_0(kr) \mathcal{H}[U(r')] \frac{\sinh kH \sinh k(z+H) - kH \sinh kz}{\sinh^2 kH -k^2 H^2}. \end{gather}
For the disc, 
$u, w$ and 
$p$ are clearly seen to vanish because 
$U=0$. The azimuthal velocity in the subphase 
$v$ is plotted for various depths 
$z$ in figure 2(b). It satisfies the no-slip boundary condition at 
$z=-H$ and the continuity condition 
$v=V$ at 
$z=0$. For 
$-H< z<0$, the exponential decay of the integrand as 
$k\to \infty$ in (4.21) ensures that the subsurface velocity field is smooth.
4.2. A formulation as dual integral equations
As a check, and as an extension of more classical approaches, we now reframe the problem of finding 
$V$ as solving a pair of integral equations, one each on the adjoining intervals 
$(0,R)$ and 
$(R,\infty )$. Define 
$a(k)$ to be the Hankel transform of 
$g(r)$
\begin{equation} g(r) = \int_0^\infty \, \mathrm{d}k\, k a(k)J_1(kr), \end{equation}which, by virtue of (3.7), implies
\begin{equation} V(r) =\frac{1}{\mu} \int_0^\infty \, \mathrm{d}k \tanh kH\,a(k) J_1(kr). \end{equation}Substitute into (4.1) to obtain
\begin{gather} \int_0^\infty \, \mathrm{d}k\, k \left(\frac{\bar\eta}{\mu}k \tanh kH + 1\right) a(k)J_1(kr) = 0,\quad 0< r< R, \end{gather}
\begin{gather}\int_0^\infty \, \mathrm{d}k\, k a(k)J_1(kr)=0,\quad r>R. \end{gather} Equations (4.27) and (4.28) constitute a set of dual integral equations with Bessel-type kernel in the unknown 
$a(k)$ that is homogeneous for 
$r>R$. Several powerful methods, such as those developed by Busbridge (Reference Busbridge1938), Cooke (Reference Cooke1956) or Sneddon (Reference Sneddon1975), have been introduced over the years in order to solve problems of this form. Perhaps the most computationally straightforward of these is the method of Tranter (Reference Tranter1954), who found an explicit countable basis for 
$a(k)$ and cast the problem as an infinite linear system, thus reducing the problem of finding the azimuthal velocity field to the problem of finding a small number of coefficients. Tranter's method has previously been used in works such as Stone (Reference Stone1995), Henle & Levine (Reference Henle and Levine2009) and Martin & Smith (Reference Martin and Smith2011) to solve Cartesian or axisymmetric mixed boundary problems in fluid mechanics where the flow or a stress is prescribed on an inner region, but the procedure here requires a modification since no such information is provided in our system – in the problem at hand, the inhomogeneous forcing comes from the boundary condition (4.3). In the following section, we show how to modify Tranter's procedure in such a way that the boundary data naturally enter into the problem.
4.2.1. Solution via Tranter's method
We begin by deriving the momentum equation in Fourier space. In the case of an axisymmetric monolayer, 2-D Fourier transforms reduce to Hankel transforms
\begin{equation} \hat{\boldsymbol{f}} = \int_\mathcal{D} \, \mathrm{d}A\, \bar\eta {\rm \Delta} \boldsymbol{U} \, {\rm e}^{-\mathrm{i} \boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{x}} ={-}2{\rm \pi} \mathrm{i} \hat{\boldsymbol{k}}^\perp\int_0^R \, \mathrm{d}r\, r \bar\eta \mathcal{L}[V](r)J_1(kr), \end{equation}and
\begin{equation} \hat{\boldsymbol{U}} = \int_\mathcal{D} \, \mathrm{d}A\, \boldsymbol{U} \, {\rm e}^{-\mathrm{i} \boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{x}} ={-}2{\rm \pi} \mathrm{i}\hat{\boldsymbol{k}}^\perp \int_0^\infty \, \mathrm{d}r\, r V(r)J_1(kr). \end{equation}
Here, 
$\mathcal {L}$ is the differential operator defined in (3.3). It is advantageous to integrate by parts
\begin{align} \hat{\boldsymbol{f}} ={-}2{\rm \pi} \mathrm{i} \bar\eta \hat{\boldsymbol{k}}^\perp \left[V'(R^-) R J_1(kR) + V(R) [J_1(kR)- kR J_0(kR)] -k^2\int_0^R \, \mathrm{d}r\, r V(r)J_1(kr)\right], \end{align}and substitute in the stress boundary condition (4.3), to arrive at the integral version of the inhomogeneous momentum equation in Fourier space
\begin{equation} -2\eta_R \varOmega R J_1(kR) = V(R) [2\eta J_1(kR) - \bar\eta kR J_0(kR)] - \bar\eta k^2 \int_0^R \, \mathrm{d}r\, r V(r) J_1(kr) - a(k).\end{equation}
Here, we have defined 
$a(k)$ to be the Hankel transform of 
$g(r)$ as in (4.25).
Tranter (Reference Tranter1954) observed that, without loss of generality, we may take our 
$a(k)$ to be of the form
\begin{equation} a(k) = k^{-\beta} \sum_{n=0}^\infty a_n J_{2n+1+\beta}(kR), \end{equation}
where 
$\beta > 0$ is arbitrary and the coefficients 
$a_n$ are unknown; such a form for 
$a(k)$ automatically satisfies the condition 
$\boldsymbol {f} = 0$ for 
$r>R$. This expression is substituted into the integral equation and projected back onto the basis to yield an infinite system of linear equations for the 
$\{a_n\}$, which are then truncated and solved to obtain 
$V$. By (3.7) and the Hankel inversion theorem
\begin{equation} V(r) =\frac{1}{\mu} \int_0^\infty \, \mathrm{d}k\, a(k) \tanh kH\, J_1(kr) = \frac{1}{\mu}\sum_{n=0}^\infty a_n \int_0^\infty \, \mathrm{d}k\, k^{-\beta} \tanh kH\, J_{2n+1+\beta}(kR) J_1(kr). \end{equation}Upon substitution into (4.32),
\begin{align} &-2\eta_R \varOmega R \mu J_1(kR) \nonumber\\ &\quad = [2\eta J_1(kR) - \bar\eta kR J_0(kR)]\sum_{n=0}^\infty a_n \int_0^\infty \, \mathrm{d}k'\, k'^{-\beta} \tanh k'H\, J_{2n+1+\beta}(k'R) J_1(k'R) \nonumber\\ &\qquad - \bar\eta k^2 \sum_{n=0}^\infty a_n\int_0^R \, \mathrm{d}r\, r J_1(kr) \int_0^\infty \, \mathrm{d}k'\, k'^{-\beta} \tanh kH\, J_{2n+1+\beta}(k'R) J_1(k'r) \nonumber\\ &\qquad - \mu k^{-\beta} \sum_{n=0}^\infty a_n J_{2n+1+\beta}(kR). \end{align}
We now multiply the equation by 
$k^{-1-\beta } J_{2m+1+\beta }(kR)$, where 
$m$ is a non-negative integer, and integrate 
$k$ from 
$0$ to 
$\infty$. This yields the infinite system of equations
\begin{equation} -2\eta_R \varOmega R \mu g_n = \sum_{m=0}^\infty a_m [(2\eta g_n - \bar\eta\xi_n) \varLambda_m - \bar\eta M_{mn} - \mu \varDelta_{mn}],\end{equation}where, using formulas from Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007),
\begin{gather} g_n = \int_0^\infty \, \mathrm{d}k\, J_1(kR) J_{2n+1+\beta}(kR) k^{{-}1-\beta} = \frac{R^\beta}{2^{1+\beta} (1+\beta)!}\delta_{n0}, \end{gather}
\begin{gather}\xi_n = \int_0^\infty \, \mathrm{d}k\, k J_0(kR) J_{2n+1+\beta}(kR) k^{-\beta} = 0, \end{gather}
\begin{gather}\varLambda_n = \int_0^\infty \, \mathrm{d}k'\, (k')^{-\beta} \tanh k'H\, J_{2n+1+\beta}(k'R) J_1(k'R) , \end{gather}
\begin{gather}M_{mn} = \int_0^\infty \, \mathrm{d}k'\, (k')^{{-}2\beta}\tanh k'H\, J_{2n+1+\beta}(k'R) J_{2m+1+\beta}(k'R), \end{gather}
\begin{align} \varDelta_{mn} &= \int_0^\infty \, \mathrm{d}k\, k^{{-}1-2\beta} J_{2n+1+\beta}(kR)J_{2m+1+\beta}(kR) \nonumber\\ &= \frac{ \beta R^{2\beta} (2\beta-1)! (m+n)!}{4^{\beta}(\beta+m-n)!(\beta-m+n)!(1+2\beta +m + n)!} \end{align}
and 
$\delta_{n0} = 1$ if 
$n = 0$ and 0 otherwise. Here, non-integer factorials assume their usual definition via the gamma function (formula 8.310.1 of Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007). In cases where a negative integer factorial is being taken in the denominator, the factorial is interpreted as infinity so that the integral vanishes.
At this stage, we choose 
$\beta =1/2$ so that the matrix 
$M_{mn}$ is nearly diagonal for large 
$k$. If 
$H\to \infty$, then 
$\beta =1/2$ ensures 
$M_{mn}$ is exactly diagonal. This choice of 
$\beta$ aids the convergence of the numerical routine by capturing the anticipated inverse square root divergence of 
$g(r)$ at the boundary and is hence optimal, although we emphasize that the routine converges to the same solution regardless of 
$\beta$.
The infinite system of (4.36) is truncated and inverted to solve for the coefficients 
$a_m$. We find that keeping the first twenty terms is generally sufficient, so the computation is fast, and the flow field at the surface can be reconstructed from the solution via (4.34). Analytical formulas exist for 
$\varLambda _n$ and 
$M_{mn}$ in the case where 
$H$ is either vanishingly small or infinitely large. In the remaining cases, these integrals must be computed numerically, for instance with an asymptotic expansion that exploits the rapid decay of 
$\tanh kH$ to unity or a specialized computational package for oscillatory integrals like the IIPBF package developed by Ratnanather et al. (Reference Ratnanather, Kim, Zhang, Davis and Lucas2014). The results obtained using Tranter's method are in excellent agreement with those obtained by using the Green's function approach, as demonstrated in figure 3.
Figure 3. Comparison of the solution of the integral equation (3.9) computed using Tranter's method (20 terms, dashed line) with the solution using the Green's function formulation (800 Chebyshev nodes on the interval 
$0< r< R$, solid line) for a disc-shaped monolayer on a subphase of infinite depth. Parameters: 
$\eta _R/\eta = 1.875, \mu R/\eta = 30$.
4.3. Asymptotic solution when 
$H\to 0$
In the limit of vanishing subphase thickness, the substrate drag dominates. This can be seen by letting 
$H\to 0$ in (2.15), which reduces to the simple condition
\begin{equation} \hat{\boldsymbol{f}} = \varGamma \hat{\boldsymbol{U}}, \end{equation}
where 
$\varGamma = \mu / H$ is the substrate drag coefficient. Inverting the Fourier transform yields a Brinkman equation inside the disc
\begin{equation} - \boldsymbol{\nabla} P + \bar\eta\boldsymbol{\nabla}^2 \boldsymbol{U} = \varGamma \boldsymbol{U},\end{equation}
whose general solution in the axisymmetric case is of the form 
$V(r) = C K_1(r/\bar {\delta }) + D I_1(r/\bar {\delta })$, where 
$\bar {\delta }^2 = \bar \eta /\varGamma, I_1$ and 
$K_1$ are modified Bessel functions and 
$C$ and 
$D$ are constants. In order to avoid a blowup at the origin, 
$C=0$; the constant 
$D$ is then found by applying the stress boundary condition (4.3) at 
$r=R$, ultimately yielding the solution
\begin{equation} V(r) = 2\varOmega \bar{\delta} \frac{\eta_R I_1(r/\bar{\delta})}{\eta I_2(R/\bar{\delta}) + \eta_R I_0(R/\bar{\delta})}\chi( r < R). \end{equation}
Note this function is discontinuous at 
$r=R$. Figure 2 shows (4.44) is the asymptotic limit of the solution to (4.1) as 
$H\to 0$. In this ‘high friction case’, the flow is largely confined to a boundary layer of width 
$\bar {\delta }$.
5. An annular domain
We now take 
$\mathcal {D}$ to be the annulus with radii 
$0 < R_i(t) < R_o(t)$. Unlike the disc case, there will be a radial component to the flow field in addition to a azimuthal one, and the odd viscosity will play a non-trivial role in this geometry (figure 1b).
The divergence of the surface velocity field is assuredly non-zero since the interface moves inward; however, for 
$R_i< r< R_o$, the divergence is still zero by assumption. This condition restricts the radial flow to something of the form
\begin{equation} U(r) = \frac{F}{r}, \quad R_i< r< R_o, \end{equation}
for some constant 
$F<0$, to be determined. The operator 
$\mathcal {L}$ annihilates 
$U$ on 
$R_i< r< R_o$, and the radial momentum equation simply reduces to
\begin{equation} -\frac{\mathrm{d}P}{\mathrm{d}r} = f,\end{equation}where
\begin{equation} f =\left.\mu \frac{\partial u}{ \partial z} \right|_{z=0}, \end{equation}
obeys (3.6) when Hankel transformed. While (5.2) seems to imply the radial bulk dynamics is completely independent of the three monolayer viscosities in this geometry, we recall the definition of 
$P$ contains 
$\eta _O$. If we take the surface pressure outside of the annulus to be zero, the radial and azimuthal stress boundary conditions for the annulus are found to be
\begin{gather} \left.-P(r)-\frac{2\eta F}{r^2} + 2\eta_O \frac{V(r)}{r} \mp \frac{\gamma}{r}\right|_{r=R_i,R_o} = 0, \end{gather}
\begin{gather}\left.-\bar\eta V'(r) +2\eta_R \varOmega + (\eta-\eta_R) \frac{V(r)}{r} + 2\eta_O \frac{F}{r^2} \right|_{r=R_i,R_o} = 0, \end{gather}
where in the radial boundary conditions, the negative sign corresponds to 
$r=R_i$ and the positive sign corresponds to 
$r=R_o$. Note that when 
$\eta _O\neq 0$, the radial and azimuthal flows are coupled through the boundary conditions.
5.1. Solution via Green's function
The solution to the annulus problem using the Green's function formulation mirrors that of the disc problem. We begin with the azimuthal velocity which is again decomposed as
\begin{equation} V(r) = V_p(r) + V_h(r), \end{equation}
where the definition of 
$V_p$ is slightly modified to account for the new limits of integration
\begin{equation} V_p(r) = \int_{R_i}^{R_o} \, \mathrm{d}r'\, G(r,r') g(r') ={-} \frac{1}{2\bar\eta r} \int_{R_i}^r \, \mathrm{d}r' \, r'^2 g(r') - \frac{r}{2\bar\eta} \int_r^{R_o} \, \mathrm{d}r'\, g(r'), \end{equation}
and the definition of 
$V_h$ remains intact. For the annulus, the constants 
$A_h$ and 
$B_h$ satisfy
\begin{gather} V_p'(R_i^+) + A_h -\frac{B_h}{R_i^2} - \frac{\eta-\eta_R}{\bar\eta R_i} \left[V_p(R_i) + A_hR_i +\frac{B_h}{R_i}\right]= \frac{2\eta_R\varOmega}{\bar\eta}, \end{gather}
\begin{gather}V_p'(R_o^-) + A_h -\frac{B_h}{R_o^2} - \frac{\eta-\eta_R}{\bar\eta R_o} \left[V_p(R_o) + A_hR_o +\frac{B_h}{R_o}\right]= \frac{2\eta_R\varOmega}{\bar\eta}, \end{gather}
where 
$R_i^+$ and 
$R_o^-$ are right and left limits, respectively. Solving for 
$A_h$ and 
$B_h$
\begin{gather} A_h = {\varOmega} +\frac{ (\eta-\eta_R)[R_iV_p(R_i) - R_oV_p(R_o)] + \bar\eta [R_o^2 V_p'(R_o^-) - R_i^2 V_p'(R_i^+)]}{2\eta_R (R_i^2-R_o^2)}, \end{gather}
\begin{gather}B_h ={-}\frac{2\eta_OF}{\eta} \!+\! \left[\frac{\eta\!-\!\eta_R}{2\eta} \left(\frac{V_p(R_i)}{R_i} \!-\! \frac{V_p(R_o)}{R_o}\right) \!+\! \frac{\bar\eta}{2\eta}[V_p'(R_o^-)-V_p'(R_i^+)]\right] \left(\frac{1}{R_o^2}-\frac{1}{R_i^2}\right)^{{-}1}. \end{gather}From (5.7), we obtain the identities
\begin{gather} R_iV_p(R_i) - R_oV_p(R_o) = \frac{1}{2\bar\eta} \int_{R_i}^{R_o} \, \mathrm{d}r'\, (r'^2 -R_i^2 )g(r'), \end{gather}
\begin{gather}R_o^2V'_p(R_o^-) - R_i^2V_p'(R_i^+) = \frac{1}{2\bar\eta} \int_{R_i}^{R_o} \, \mathrm{d}r'\, (r'^2 +R_i^2) g(r'), \end{gather}
\begin{gather}\frac{V_p(R_i)}{R_i} - \frac{V_p(R_o)}{R_o} = \frac{1}{2\bar\eta} \int_{R_i}^{R_o} \, \mathrm{d}r'\, \left({-}1 + \frac{r'^2}{R_o^2}\right) g(r'), \end{gather}
\begin{gather}V_p'(R_o^-) - V_p'(R_i^+) = \frac{1}{2\bar\eta} \int_{R_i}^{R_o} \, \mathrm{d}r'\, \left(1 + \frac{r'^2}{R_o^2} \right)g(r'). \end{gather}Substituting everything into (3.9), we finally have
\begin{align}
&\int_{R_i}^{R_o} \, \mathrm{d}r'\, G(r,r')
g(r') + r\left[{\varOmega} +\frac{\eta}{2\eta_R\bar\eta(R_i^2-R_o^2)} \int_{R_i}^{R_o}
\, \mathrm{d}r'\, r'^2 g(r') + \frac{R_i^2}{2\bar\eta
(R_i^2-R_o^2)} \int_{R_i}^{R_o} \, \mathrm{d}r'\, g(r')
\right]\nonumber\\
&\quad +\frac{1}{r} \left[-\frac{2\eta_O F}{\eta} + \frac{\rho}{2\bar\eta R_o^2}\int_{R_i}^{R_o} \,
\mathrm{d}r'\, r'^2 g(r') + \frac{\eta_R \rho}{2\bar\eta
\eta} \int_{R_i}^{R_o} \, \mathrm{d}r'\, g(r')\right] =
\frac{1}{\mu}\int_{R_i}^{R_o} dr'~r' L(r,r') g(r'),
\end{align}
where 
$\rho = [(1/R_o)^2 -(1/R_i)^2]^{-1}$. As before, this integral equation can be inverted for 
$g$ as a function of 
$F$, which in turn gives 
$V$ in terms of 
$F$. Note that the definition of 
$L_d$ changes due to the new limits of integration
\begin{align} L_d(r) &= \int_{R_i}^{R_o} \, \mathrm{d}r'~r'\bar{L}(r,r') \nonumber\\ &= \int_0^\infty \, \mathrm{d}k\, J_1(kr) \int_{R_i}^{R_o}\mathrm{d}r'\,r'J_1(kr') \nonumber\\ &={-}\frac{R_i^3}{6r^2}{}_3F_2\left[\left\{\frac{1}{2},\frac{3}{2}, \frac{3}{2}\right\},\left\{2,\frac{5}{2}\right\}, \frac{R_i^2}{r^2}\right] \nonumber\\ &\quad + \frac{r}{32}\left(16\log \frac{4R_o}{r}-8-\frac{3r^2}{R_o^2} {}_4F_3\left[\left\{1,1, \frac{3}{2},\frac{5}{2}\right\},\{2,2,3\}, \frac{r^2}{R_o^2}\right]\right). \end{align}
Also note that if 
$R_i$ and 
$F$ are allowed to approach zero, the integral equation reduces to that of the disc case, (4.12), with 
$R=R_o$.
With the azimuthal velocity essentially solved, we turn to the radial velocity. Since 
$\mathcal {L}[U]$ vanishes in (3.8), we are left with
\begin{equation} \frac{F}{r} = \frac{1}{\mu} \int_{R_i}^{R_o}\mathrm{d}r' \, r' M(r,r') f(r'). \end{equation}
By expanding 
$M(r,r')$ in the same way as 
$L(r,r')$ above, this integral equation can be numerically inverted to solve for 
$f/F$. Since 
$V$ is known (up to 
$F$), the constant 
$F$ can be determined from the fact that
\begin{equation} P(R_o^-) - P(R_i^+) ={-}\int_{R_i}^{R_o} \, \mathrm{d}r'\, f(r'), \end{equation}
along with (5.4). This completes the solution of the instantaneous surface flow field. Representative solutions for different values of 
$H$ are shown in figure 4. The scaling arguments established in the disc case carry over: in the 
$\bar \beta \gg 1$ limit, the appropriate length scale for the velocity is 
$R$, while it is 
$\ell _{SD}$ in the small 
$\bar \beta$ limit. However, the radial and azimuthal components have different time scales. For instance, in figure 4, 
$U$ and 
$V$ exhibit a disparity in scale, with 
$U$ being approximately a hundred times smaller in magnitude than 
$V$. This can be traced back to (2.10) and (2.11), which show that in the 
$\eta _O=0$ case, the azimuthal motion originates from the rotational drive while the radial motion originates from the line tension. The corresponding time scales are 
$\tau _1 = \varOmega ^{-1}$ and 
$\tau _2=\eta R_i/\gamma$; for the curves plotted in figure 4, the ratio of these time scales was taken to be 
$\tau _1/\tau _2 = \gamma /(\eta R_i \varOmega ) = 0.01$. Thus, for 
$\bar \beta \gg 1, V\sim R_i/\tau _1=R_i\varOmega$ and 
$U\sim R_i/\tau _2=\gamma /\eta$, while for 
$\bar \beta \ll 1, V\sim \ell _{SD}/\tau _1= \eta \varOmega /\mu$ and 
$U\sim \ell _{SD}/\tau _2 = \gamma /(\mu R_i)$. Figure 5 illustrates the subphase velocity field for the 
$H=R_i$ case calculated from (4.21)-(4.23). Dynamics follows from advection of the domain boundary according to the kinematic boundary condition, which we describe in more detail in § 5.4.
Figure 4. The instantaneous surface azimuthal velocity 
$V$ (a) and radial velocity 
$U$ (b) as well as the monolayer pressure 
$P$ (c) due to an annular monolayer of radius 
$R_o/R_i = 2$ as subphase depth 
$H$ is varied. As 
$H\to 0$, boundary layers of width 
$\bar {\delta } = \sqrt {\bar \eta H/\mu }$ becomes visible in 
$V$. Parameters: 
$\eta _R/\eta =1.875 , \mu R/\eta = 30, \eta _O/\eta = 0, \tau _1/\tau _2=\gamma /(\eta R_{i}\varOmega ) = 0.01$. Since 
$\ell _{SD}\ll R, V$ scales like 
$\ell _{SD}/\tau _1 = \eta \varOmega /\mu$ and 
$U$ scales like 
$\ell _{SD}/\tau _2=\gamma /(\mu R_i)$.
Figure 5. The azimuthal component 
$v$ (a), radial component 
$u$ (b) and vertical component 
$w$ (c) of the velocity field 
$\boldsymbol {u}$ in the subphase at various depths for an annular monolayer of 
$R_o/R_i=2$ and 
$H=1$. Parameters: 
$\eta _R/\eta =1.875 , \mu R/\eta = 30, \eta _O/\eta = 0, \gamma /(\eta R_{i}\varOmega ) = 0.01$.
5.2. Formulation as triple integral equations
Similar to the disc case, we can convert the problem into two sets of triple integral equations. Keeping the definition of 
$a(k)$ from (4.25), the azimuthal set is
\begin{gather} \int_0^\infty \, \mathrm{d}k \, k a(k) J_1(kr)=0, \quad 0< r< R_i, \end{gather}
\begin{gather}\int_0^\infty \, \mathrm{d}k\, k \left(\frac{\bar\eta}{\mu} k \tanh kH + 1\right) a(k)J_1(kr)=0 , \quad R_i< r< R_o, \end{gather}
\begin{gather}\int_0^\infty \, \mathrm{d}k \, k a(k) J_1(kr)=0,\quad r>R_o, \end{gather}and the radial set is
\begin{gather} \int_0^\infty \, \mathrm{d}k\, k b(k) J_1(kr)=0, \quad 0< r< R_i, \end{gather}
\begin{gather}\int_0^\infty \, \mathrm{d}k\, b(k)[A(kH) + B(kH)]^{{-}1}J_1(kr)=\frac{\mu F}{r} ,\quad R_i< r< R_o, \end{gather}
\begin{gather}\int_0^\infty \, \mathrm{d}k \, k b(k) J_1(kr)=0 ,\quad r>R_o, \end{gather}where we have defined
\begin{equation} b(k) = \int_0^\infty \, \mathrm{d}r\, r f(r) J_1(kr), \end{equation}and used (3.6).
Difficulties similar to the dual integral equations of the previous section plague the azimuthal triple integral equations. In particular, ((5.20)–(5.22)) again erroneously appear to be homogeneous because the problem as stated above is not closed. Unfortunately, no convenient basis analogous to Tranter's for the disc appears to resolve the problem; attempting to use the Tranter basis as before will lead to a result that is discontinuous at 
$r=R_i$, as the basis is ‘unaware’ of the divergence there. One possible workaround is to Fourier transform the equations to include boundary conditions and discretize the equations directly on the interval 
$(0,\infty )$ to yield a large linear system of equations. However, this method converges very slowly and performing it repeatedly to time step the kinematic boundary condition is not practical. Thus, we content ourselves with solving the equations using the aforementioned singular integral formulation.
On the other hand, several methods exist for solving the radial triple integral equations. One approach due to Cooke (Reference Cooke1965) involving Erdélyi–Kober operators (generalized fractional derivatives) can be used to solve the radial equations, up to the unknown constant 
$F$. In the simpler case of 
$H\to \infty$, a different method also devised by Cooke (Reference Cooke1963) involving rewriting the equation as a composition of Abel transforms (see Noble Reference Noble1958) can be used as well. These methods reduce the set of triple integral equations to a single Fredholm integral equation, which must be numerically solved, ultimately making them more work than the solution we have presented.
It is notable that the radial set of triple integral equations has a simple exact solution in the limit 
$H,R_o\to \infty$, as found by Alexander et al. (Reference Alexander, Bernoff, Mann, Mann and Zou2006). In this limit
\begin{equation} b(k)= \mu F\frac{\sin kR_i}{k R_i},\end{equation}which, via the Hankel inversion theorem, leads to (cf. Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007) formula (6.693.1) and differentiate under the integral)
\begin{equation} -\frac{\mathrm{d}P}{\mathrm{d} r} =\frac{\mu F}{r\sqrt{r^2-R_i^2}}\chi({r>R_i}), \end{equation}so that the pressure (with constant of integration zero) is
\begin{equation} P ={-}\mu F \frac{\cos^{{-}1}(R_i/r)}{R_i}, \end{equation}
for 
$r>R_i$, and
\begin{equation} U = \frac{F}{r}\left(1-\sqrt{1-\frac{r^2}{R_i^2}}\right), \end{equation}
for 
$0< r< R_i$. For 
$r>R_i, U=F/r$ as previously.
5.3. Asymptotic solution when $H\to 0$
As with the disc-shaped domain, the small 
$H$ case reduces to a Brinkman equation; we will impose 
$\boldsymbol {U} = \boldsymbol {0}$ on 
$\mathcal {D}^C$ so that 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {U} = 0$ everywhere and the drag is isotropic. The general solution in the axisymmetric case is
\begin{gather} U(r) = \frac{F}{r}, \end{gather}
\begin{gather}V(r) = CK_1\left(\frac{r}{\bar{\delta}}\right) + DI_1\left(\frac{r}{\bar{\delta}}\right), \end{gather}
\begin{gather}P(r) ={-}\varGamma F \log r + G. \end{gather}
The four unknowns 
$C,D,F$ and 
$G$ are found by substituting these expressions into the four boundary conditions (5.4) and (5.5). Exact but cumbersome expressions for these constants can be found; in the interest of brevity, we omit them. However, in the simple case of zero odd viscosity, the constant 
$F$ is found to be
\begin{equation} F ={-}\frac{\gamma [(1/R_i) + (1/R_o)]}{\varGamma \log (R_o/R_i) +2\eta [(1/R_i^{2}) -(1/R_o^{2})]},\end{equation}
which, through the kinematic boundary condition, determines the size of the cavity. Note that 
$R_o$ is related to 
$R_i$ through the monolayer incompressibility constraint: 
$A = {\rm \pi}(R_o^2 -R_i^2)$. We remark that as in the disc case, boundary layers of width 
$\bar {\delta } = \sqrt {\bar \eta /\varGamma }$ are visible in the azimuthal velocity field in this limit (figure 4). Figure 6 shows the parameter 
$F$ as a function of time as 
$\eta _O$ is varied. The impact of this parameter on the closing of the cavity is discussed in the next section.
Figure 6. The parameter 
$F$ as a function of time for various odd viscosities in the high friction case. The initial slope depends on 
$\gamma$ and 
$\eta _O$, where line tension and odd viscosity are the dominant contributors to the radial dynamics. As the hole closes and the hole radius becomes comparable to 
$\bar {\delta }, F$ reaches a local minimum and then rapidly approaches zero in a small shear viscosity dominated regime so that 
$F(t)$ appears nearly vertical. Initially, 
$R_{o,0}/R_{i,0} = 5$. Curves from right to left: 
$\eta _O/\eta = 0,0.5,1.0,1.5,2.0$. Parameters: 
$\eta _R/\eta = 1.857\times 10^{-2}, \varGamma R_{i,0}^2/\eta = 4.074\times 10^{2}, \gamma /(\eta R_{i,0} \varOmega ) = 5.214\times 10^{-4}$.
5.4. Hole closure dynamics and the effect of odd viscosity
For the disc-shaped domain, the absence of a radial velocity means that the boundary never moves. The odd viscous stresses are in the radial direction but are offset by the pressure and hence have no effect on the domain shape or flow field.
For the annular domain, the kinematic boundary condition states that the radii of the circular boundaries will change according to the local radial surface velocity
\begin{equation} \frac{\mathrm{d}R_i}{\mathrm{d}t} = U(R_i) =\frac{F}{R_i},\end{equation}
and similarly for the outer radius. Equivalently, the outer radius can be found by applying the constraint that the monolayer area 
${\rm \pi} (R_o^2 -R_i^2)$ is constant. Equation (5.35) can be numerically integrated to find 
$R_i(t)$, up until 
$R_i=0$, at which point the circular disc case is recovered. At each time step, we must solve for 
$F$ by solving for the flow field using the procedure outlined in the previous section.
We first obtain some simple limits for the high friction case. Following the discussion in Appendix B, the Saffman–Delbrück length for the high friction case is 
$\bar \delta$. For 
$\bar \delta \gg R_i$, the expected radial velocity scale is 
$R_i/\tau _2 = \gamma /\eta$, while for 
$\bar {\delta } \ll R_i$, it is 
$\bar {\delta }/\tau _2 = \gamma \bar {\delta }/(\eta R_i)$. In the initial phase of the experiment, the cavity radius is large compared with 
$\bar \delta$, so 
$F\sim \gamma \bar {\delta }/\eta$. However, as the hole is just about to close, 
$F\sim \gamma R_i/\eta$. These two regimes are illustrated in figure 6, where 
$F$ is initially nearly constant in the 
$\eta _O=0$ case and transitions to a linear regime with large slope when 
$R_i$ becomes comparable to 
$\bar \delta$.
An asymptotic analysis of the small 
$R_i$ limit reveals the closure time is finite. The argument proceeds as follows: for simplicity, assume that 
$H\to \infty$ and 
$V(R_i)\to 0$. If 
$R_i \ll R_o$, we can use the exact solution (5.27) to calculate the pressure difference 
$P(R_o) - P(R_i) = -\mu {\rm \pi}F/(2R_i)$, which is negligible compared with the shear viscous stress which scales like 
$F/R_i^{2}$. The boundary condition (5.4) shows that this stress must be balanced by the line tension, from which we find 
$F \sim -\gamma R_i/(2\eta )$. The kinematic boundary condition (5.35) then implies 
$\mathrm {d}R_i/\mathrm {d}t$ is constant in this limit so that there is no blowup. Prior to this regime, the cavity area decreases at a nearly constant rate, and is well approximated by 
$A_C = {\rm \pi}R_{i,0}^2 + 2{\rm \pi} t F_0$, where 
$R_{i,0}$ is the initial cavity radius and 
$F_0$ is the value of 
$F$ at time 
$t=0$. Figure 7 illustrates these different regimes and compares the cavity radius and area as functions of time for different odd viscosities. Note that in the high friction case, the 
$t\to t^*$ shear viscosity-dominated regime occurs at a cavity radius comparable to 
$\bar \delta$, which is itself comparable to the particle size, and is thus not expected to be experimentally detectable.
Figure 7. (a) Cavity radius as a function of time for various odd viscosity values in the high friction case. Initially, 
$R_{o,0}/R_{i,0} = 5$. Curves from right to left: 
$\eta _O/\eta = 0,0.5,1.0,1.5,2.0$. Parameters: 
$\eta _R/\eta = 1.857\times 10^{-2}, \varGamma /(\eta R_{i,0}^2) = 4.074\times 10^{2}, \gamma /(\eta R_{i,0} \varOmega ) = 5.214\times 10^{-4}$. (b) Corresponding area of cavity 
$A_C$ as a function of time using the same parameters and colour scheme. Increasing odd viscosity changes the apparent concavity of the cavity area vs time curve. Inset: larger version of the 
$\eta _O=0$ curve near 
$t=t^*$, where the hole size is comparable to the penetration depth 
$\bar \delta$. In this limit, the closing is dominated by shear viscosity and the area decreases quadratically.
The most noticeable effect of odd viscosity is that it decreases the time it takes for the hole to close. If 
$\eta _O = 0$, the closure time is independent of 
$\varOmega$ and 
$\eta _R$. On the other hand, a non-zero odd viscosity couples the azimuthal and radial velocities via the stress boundary conditions. At both the inner and outer boundaries, the line tension forces point radially inward toward the origin. The odd stress is oriented inward at the outer boundary but outward at the inner boundary; naively, one may think this causes the hole to close slower. However, this argument does not account for the pressure. As a consequence of domain incompressibility, the rate of change of the inner radius must be larger than that of the outer radius
\begin{equation} \frac{\mathrm{d} R_i}{\mathrm{d} t } = \frac{R_o}{R_i} \frac{\mathrm{d}R_o}{\mathrm{d}t}. \end{equation}
This restriction implies that if the outer boundary is moving in faster with non-zero 
$\eta _O$, so too must the inner boundary. The odd viscosity also changes the concavity of the area vs time curve. If 
$\eta _O\ll \eta$, the curve is concave (except for the shear viscosity-dominated region when the hole is very small, where it is always convex) but if 
$\eta _O$ is sufficiently large, it experiences regions of convexity as well. This behaviour is not observed in the low friction case (figure 8).
Figure 8. Cavity radius (a) and area (b) as functions of time for the low friction 
$H\to \infty$ case. Increasing the odd viscosity decreases the time needed to close the hole. The quadratic behaviour of the cavity area in the shear viscosity-dominated regime is visible at larger radii compared with the high friction case. Parameters: 
$R_{o,0}/R_{i,0} = 5, \eta _R/\eta = 1.857\times 10^{-2}, \gamma /(\eta R_{i,0}\varOmega ) = 5.214\times 10^{-4}, \mu R_{i,0}/\eta = 8.79$. In the figure on the right, the transition from linear to quadratic behaviour takes place when 
$R_i\sim \ell _{SD}=0.11R_{i,0}$.
For the low friction case, we find the same general trends (figure 9). For zero odd viscosity, infinite substrate depth and asymptotically large outer radius 
$R_o/R_i\to \infty$, a complete analytical description is possible using (5.27). The cavity area as a function of time in this case is given by 
$A_C={\rm \pi} R_i(t)^2$, where
\begin{equation} R_i(t) = \frac{2\eta}{{\rm \pi} \mu} \left[ \sqrt{1 +\frac{{\rm \pi} \gamma \mu (t^*-t)}{2\eta^2}}-1\right], \end{equation}
and 
$t^* = (2\eta R_{i,0}+{\rm \pi} \mu R_{i,0}^2/4)/\gamma$ is the closing time of the cavity written in terms of the initial radius 
$R_{i,0}$ at time 
$t=0$. When 
$R_i$ is large compared with 
$\ell _{SD}$, the area changes linearly in time with a constant 
$2{\rm \pi} F\sim -4\gamma /\mu$. In the limit where 
$R_i\ll \ell _{SD}$, we find that 
$F\sim -\gamma R_i/(2\eta )$ just as in the high friction case (in fact, this limit is independent of 
$H$). Figure 8 shows the analogous radius and area curves as a function of time for the low friction case. The small radius viscosity-dominated regime is more clearly visible. The reader is referred to Jia & Shelley (Reference Jia and Shelley2022) for details about this analytically tractable case.
Figure 9. Increasing substrate height drastically decreases the time it takes for the hole to close. (a) Parameters: 
$R_{o,0}/R_{i,0} = 5, \eta _R/\eta = 1.857\times 10^{-2}, \eta _O/\eta = 0, \mu R_{i,0}/\eta = 8.79, \gamma /(\eta R_{i,0}\varOmega ) = 5.214\times 10^{-4}$. (b) Same parameters except 
$\eta _O/\eta = 0.5$.
5.5. Odd viscosity drives hole closure
We next turn to the case where line tension is small compared with 
$\mu R_i^2\varOmega$. As seen by letting 
$\gamma \to 0$ in (5.34), an ordinary Brinkman fluid (i.e. one with 
$\eta _O=0$) with no line tension has 
$F=0$ so that no radial dynamics occurs. However, for an odd fluid with a non-trivial tangential velocity, it is possible for odd viscosity to serve as a driver of normal motion, even when there is no line tension. In this case, we might expect the constant 
$F$ to be proportional to 
$\eta _O\varOmega /\varGamma$, on the basis of dimensional analysis. While the exact expression for 
$F$ in this scenario is still too cumbersome to reproduce here, in the limit of 
$\bar {\delta } \ll R_i < R_o$, we find
\begin{equation} F\sim{-}\frac{2\eta_O\varOmega}{\varGamma} \frac{ \eta_R\bar{\delta} [(1/R_i)+(1/R_o)]}{ \bar{\eta}\log(R_o/R_i)}. \end{equation}
There are prominent similarities between this expression and (5.34), highlighting the interpretation of odd viscosity as an additional line tension-like stress that in this geometry drives radial motion. In the opposite limit of 
$R_o \gg \bar {\delta }\gg R_i$, the scaling is quadratic in 
$R_i$
\begin{equation} F\sim{-}\frac{\eta_O\eta_R \varOmega R_i^2}{\eta^2+\eta_O^2},\end{equation}
which, using (5.35), implies that the radius decays exponentially with time scale 
$[\eta _O\eta _R\varOmega /(\eta ^2+\eta _O^2)]^{-1}$. Consequently, in the zero line tension case, the cavity does not close up in finite time, as figure 10 shows. However, this limit is not expected to apply to experiments as it requires 
$R_i$ to be much smaller than 
$\bar \delta$, which is itself comparable to the particle size.
Figure 10. (a) The parameter 
$F$ as a function of time for various odd viscosities in the high friction case, with no line tension. Parameters same as in figure 6, with the exception of 
$\gamma$, which has been set to zero. When 
$\eta _O=0$, inhomogeneous forces are absent in the system, and 
$F=0$, corresponding to no motion. Increasing 
$\eta _O$ increases the radial stress, showing that odd viscosity alone can drive the radial dynamics. Dashed black line represents the asymptotic values of 
$F$ for the 
$\eta _O/\eta = 0.5$ line. Inset: larger version of the 
$\eta _O/\eta =0.5$ curve and its asymptotic approximation on a semilogarithmic scale, showing long term exponential decay of 
$F$ as a function of 
$t$. The plateau corresponds to the smallest positive floating-point number representable in MATLAB. (b) The corresponding cavity radius as a function of time, also exhibiting an exponential decay as 
$R_i$ becomes comparable to 
$\bar \delta$. Without line tension, the hole does not close in finite time. (c) The corresponding cavity area as a function of time.
While a full analytical description is not available for the low friction case with odd viscosity in the 
$R_o/R_i\to \infty$ limit, (5.27) is still valid and (5.4) thus provides a relation between 
$F$ and 
$V(R_i)$
\begin{equation} F={-}\frac{R_i(\gamma -2\eta_OV(R_i))}{2\eta + {\rm \pi}\mu R_i}. \end{equation}
As for the high friction case, if 
$\gamma \to 0$, the odd viscous stress becomes the driver of radial motion (recall that generally 
$V(R_i)<0$ for the annulus). This expression shows that, similar to the 
$\eta _O=0$ case, the viscous dissipation is mostly due to 
$\mu$ when 
$R_i \gg \ell _{SD}$ and due to 
$\eta$ when 
$R_i\ll \ell _{SD}$. However, if 
$V(R_i)>0$, which is possible in certain parameter regimes such as the large 
$\eta$ limit, this expression actually implies that the hole expands. Figure 11 shows the cavity area as a function of time for the low friction case with near infinite subphase depth. As with the low friction case, the area decays exponentially, with 
$F$ scaling as in (5.39).
5.6. Conclusion and future work
We have developed a formulation for the dynamics of an active, chiral surface phase coupled to a passive fluid underneath. We showed how to formulate the problem as calculating the surface velocity, given the surface stress, using a Green's function; this formulation is highly general and could be used to model other types of active (or passive) surface phases, or to study the dynamics in more complicated, multi-connected domains with little to no modification. Using analytical and numerical methods, we proceeded to calculate the velocity fields for a disc-shaped and an annular monolayer. For the case of a disc-shaped monolayer, a modification of Tranter's method allowed for a semi-analytical description and efficient numerical solution. For the case of an annulus, we thoroughly explored the effects of odd viscosity on the closing of a 2-D circular cavity. Our main results include a decrease in the cavity closure time in the presence of odd viscosity and a change in concavity of the cavity area vs time curve as 
$\eta _O$ is increased in the high friction case. These results may provide another way to experimentally estimate the odd viscosity coefficient. We have also seen that in the absence of edge tension, odd viscosity alone can drive cavity closure.
Ongoing work is focused in several different directions. Firstly, a boundary integral formulation for the high friction case to handle non-axisymmetric shapes is under development. A full numerical formulation of the low friction case, much less its linear stability theory, is particularly challenging. Great care is needed to accommodate the divergent surface stresses, which are a fundamental part of the basic model. As a basic problem in applied mathematics, it would be interesting to find an analogue to the countable Tranter basis for the annular case and in that way develop a near analytical solution for its dynamics.
Acknowledgements
The authors are grateful to F. Balboa-Usabiaga, E. Bililign and Y. Ganan for helpful discussions.
Funding
W.T.M.I. acknowledges support from the National Science Foundation under awards DMR-2011854 (University of Chicago MRSEC) and DMR-1905974. M.J.S. acknowledges support by the National Science Foundation under awards DMR-1420073 (NYU MRSEC) and DMR-2004469.
Declaration of interests
The authors report no conflict of interest.
Appendix A. The physical system and experimental parameter values
The active chiral fluid we consider is a monolayer composed of thousands to millions of hematite particles, each roughly 1.6 μm in size and equipped with a magnetic dipole moment. The colloids are suspended in water and sedimented onto either a glass slide or an air/water interface; we refer to the former as the ‘high friction case’ and the latter as the ‘low friction case’. Note that the particles are denser than water so that in the low friction case, the monolayer is found at the bottom of the water ‘subphase’, which is a top-down reflection of what is depicted in the schematic in figure 1. For the monolayer sizes considered here (
$R\lesssim 500\ {\rm \mu}{\rm m}$), the interface is well approximated by an infinite plane, so for convenience, we may take the reflected configuration as our model without affecting any of our results. The depth of the water subphase, 
$H$, is typically comparable to 
$R$ in low friction experiments.
Under the application of an external rotating magnetic field, the particles spin; for frequencies in the range of roughly 
$\varOmega = 2$ to 12 Hz, the particles’ rotational inertia is negligible so that the dipole moments are effectively always aligned with the external magnetic field. Since the average magnetic interaction is attractive, the system experiences effective surface and line tensions that form a cohesive 2-D incompressible fluid. Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019) showed experimental examples of fluidic behaviour and put forth a descriptive zero Reynolds number hydrodynamic theory accounting for three kinds of bulk viscous interparticle stresses: a shear viscous stress arising from attractions between neighbouring dipoles, a rotational stress arising from rotor–rotor friction and an odd stress possibly arising from the collisions of rotating particles.
Rheological tests by Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019) suggest that the shear viscosity 
$\eta$ of the colloidal fluid is around fifty times greater than its rotational viscosity 
$\eta _R$, while fitting the dispersion relation of low friction edge waves suggests that the odd viscosity 
$\eta _O$ is comparable in magnitude to the shear viscosity: 
$\eta = 4.9\pm 0.2 \times 10^{-8}$ Pa m s, 
$\eta _R = 9.1\pm 0.1 \times 10^{-10}$ Pa m s and 
$\eta _O = 1.5\pm 0.1 \times 10^{-8}$ Pa m s. In the high friction case, these stresses are balanced against an external substrate friction 
$\varGamma = 2.49\pm 0.03 \times 10^3$ Pa s m
$^{-1}$ that is generally found to be isotropic and proportional to the monolayer velocity. In the more complicated low friction case, the external forcing comes from the shear stress due to the motion of the fluid subphase with viscosity 
$\mu$, which is intimately coupled to that of the monolayer. At the boundary, the internal stresses are balanced by an edge tension 
$\gamma = 2.3\pm 0.2 \times 10^{-13}$ N. The theoretical analysis in Soni et al. (Reference Soni, Bililign, Magkiriadou, Sacanna, Bartolo, Shelley and Irvine2019) is restricted to the simpler high friction case of a monolayer situated on glass substrate; here, we will focus on the more general low friction case of a fluid subphase.
Appendix B. Non-dimensional groups
It is instructive to consider the dimensionless versions of the monolayer momentum equation and associated boundary conditions. In the case of a monolayer with length scale 
$R$, we take velocity scale to be 
$\bar {U} = \varOmega R$ and the pressure scale to be 
$\mu \bar {U} = \mu \varOmega R$. Note that there are other possible choices for the velocity scale associated with the line tension, such as 
$\bar {V} = \gamma /(\mu R)$. In keeping with the experiments, where the time scale of rotation is much smaller than that of relaxation due to line tension, we elect to scale velocities by 
$\bar {U}$ instead of 
$\bar {V}$. Temporarily identifying dimensionless quantities with their dimensional counterparts, the momentum equation for the monolayer in arbitrary coordinates that are consistent with the Frenet frame at the boundary becomes
\begin{equation} -\boldsymbol{\nabla} P + \frac{\eta + \eta_R}{\mu R}{\rm \Delta} \boldsymbol{U} = \left. \frac{\partial \boldsymbol{u}}{\partial z} \right|_{z=0}, \end{equation}which reveals the ratios of two types of Saffman–Delbrück length to the monolayer size as two dimensionless parameters
\begin{gather} \beta_S = \frac{\eta/\mu }{ R}, \end{gather}
\begin{gather}\beta_R = \frac{\eta_R/\mu}{ R}. \end{gather}Rescaling the stress boundary conditions (2.10) and (2.11) in the same manner yields
\begin{gather} -P + 2\beta_S ( T_s -\kappa N) + 2\beta_O(N_s+\kappa T) |_{\partial\mathcal{D}}= \alpha \kappa|_{\partial \mathcal{D}}, \end{gather}
\begin{gather}-(\beta_S + \beta_R)\omega - 2 \beta_S( N_s + \kappa T) + 2\beta_O(T_s - \kappa N) |_{\partial \mathcal{D}}={-} 2\beta_R |_{\partial\mathcal{D}}, \end{gather}where
\begin{equation} \alpha= \frac{\gamma}{\mu R^2 \varOmega} = \frac{\bar{V}}{\bar{U}}, \end{equation}is the ratio of the two aforementioned velocity scales and
\begin{equation} \beta_O = \frac{\eta_O/\mu}{R}. \end{equation}
Thus, the five dimensionless parameters for the low friction problem are 
$\alpha, \beta _S, \beta _R, \beta _O$ and 
$\zeta = H/R$.
The analysis proceeds nearly identically for the high friction Brinkman equation, with one modification: the velocity scale 
$\bar {V}$ is written in terms of the substrate friction, becoming 
$\gamma /(\varGamma R^2)$. We find the corresponding definitions of 
$\beta _S, \beta _R$ and 
$\beta _O$
\begin{equation} \beta_S = \frac{\eta }{\varGamma R^2},\quad\beta_R = \frac{\eta_R }{\varGamma R^2},\quad\beta_O = \frac{\eta_O }{\varGamma R^2}. \end{equation}
Note that, in terms of the penetration depth 
$\bar \delta$,
\begin{equation} \beta_S + \beta_R = \frac{\bar\delta^2}{R^2}. \end{equation}Appendix C. An infinite strip
Here, we consider the flow field when 
$\mathcal {D}$ is an infinite strip of half-width 
$R$ oriented axially along the 
$y$-axis. The flow is assumed to be steady and unidirectional. For this section, we will use Cartesian coordinates so that the flow field may be expressed as 
$\boldsymbol {u}(x,y,z) = v(x,z)\hat {\boldsymbol {y}}$, with 
$\boldsymbol {U}(x,y) = \boldsymbol {u}(x,y,0) = V(x)\hat {\boldsymbol {y}}$. This type of flow field satisfies 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {U} = 0$ on the entire surface, so that (2.15) applies. Defining 
$g(x) = \mu \partial v/\partial z|_{z=0}$, the 
$x$-component of the momentum equation inside the monolayer is
\begin{equation} -\frac{\mathrm{d}P}{\mathrm{d}\kern0.7pt x}+ \bar\eta \frac{\mathrm{d}^2V}{\mathrm{d}\kern0.7pt x^2} = g, \end{equation}
when 
$|x|< R$. On the other hand, when 
$|x|>R$, we have 
$g=0$. At the boundary, 
$\kappa =0$ and 
$\hat {\boldsymbol {n}} = \pm \hat {\boldsymbol {x}}$, so the boundary conditions (2.9) for this geometry are
\begin{equation} -P|_{x={\pm} R} = 0 \quad \text{and} \quad \left.\bar\eta \frac{\mathrm{d}V}{\mathrm{d}\kern0.7pt x} - 2\eta_R \varOmega\right|_{x ={\pm} R} =0,\end{equation}
where the outer pressure has been taken to be zero. Analogous to the disc case, the odd viscosity does not enter explicitly in this strip geometry. Note that because the flow is unidirectional, the pressure 
$P$ inside the domain is harmonic. Since (C2a,b) shows 
$P$ vanishes along its boundary, 
$P$ must be zero everywhere. The momentum equation is then simply
\begin{equation} \bar\eta \frac{\mathrm{d}^2V}{\mathrm{d}\kern0.7pt x^2}\chi(|x|< R) = g. \end{equation}
Since the flow field has an odd symmetry, we define 
$a(k)$ to be the Fourier sine transform of 
$g(x)$
\begin{equation} a(k) = \int_{0}^\infty \, \mathrm{d}\kern0.7pt x\, g(x) \sin kx \Longleftrightarrow g(x) = \frac{2}{\rm \pi}\int_0^\infty\mathrm{d}k\, a(k)\sin kx. \end{equation}
We take the Fourier sine transform of the momentum equation in 
$x$ to obtain
\begin{equation} \bar\eta\int_{0}^R \, \mathrm{d}\kern0.7pt x \frac{\mathrm{d}^2V}{\mathrm{d}\kern0.7pt x^2} \sin kx = a(k). \end{equation}Integration by parts yields
\begin{equation} 2\eta_R\varOmega \sin kR - \bar\eta V(R) k \cos kR - \bar\eta k^2 \int_0^R\mathrm{d}\kern0.7pt x\, V(x)\sin kx = a(k), \end{equation}
where the fact that 
$V(x)$ is odd and (C2a,b) have been used to simplify boundary terms. Since 
$\sin z = \sqrt {{\rm \pi} z/2} J_{1/2}(z)$, this equation is amenable to Tranter's method, which we now demonstrate; the prescription is nearly identical to that of the disc geometry given in § 4.2.1. Naturally, adapting the Green's function formulation as in § 4.1 yields an identical answer.
Adapting Tranter (Reference Tranter1954), we let
\begin{equation} a(k) = k^{1/2-\beta} \sum_{n=0}^\infty a_n J_{2n+1/2+\beta}(kR), \end{equation}
where 
$\beta >0$ is arbitrary and the coefficients 
$\{a_n\}$ are unknown. Note that (2.15) combined with (C4) and (C7) implies
\begin{equation} V(x) = \frac{2}{{\rm \pi}\mu} \sum_{n=0}^\infty a_n \int_0^{\infty} \, \mathrm{d}k\, k^{{-}1/2-\beta} \tanh kH\, J_{2n+1/2+\beta}(kR) \sin kx, \end{equation}for the geometry at hand. Substituting this into (C6) gives
\begin{align} &2\eta_R\varOmega \sin kR-k \cos kR\frac{2\bar\eta}{{\rm \pi} \mu} \sum_{n=0}^\infty a_n \int_0^\infty \, \mathrm{d}k'~(k')^{{-}1/2-\beta} \tanh k'H\, J_{2n+1/2+\beta}(k'R) \sin k'R \nonumber\\ &\qquad - \frac{2\bar\eta}{{\rm \pi} \mu} \sum_{n=0}^\infty a_n k^2\int_0^R \, \mathrm{d}\kern0.7pt x~ \sin kx \int_0^\infty \, \mathrm{d}k'~(k')^{{-}1/2-\beta} \tanh k'H\, J_{2n+1/2+\beta}(k'R)\sin k'x \nonumber\\ &\quad = k^{1/2-\beta} \sum_{n=0}^\infty a_n J_{2n+1/2+\beta}(kR). \end{align}
Multiplying both sides by 
$k^{-3/2-\beta }J_{2m+1+\beta }(kR)$, where 
$m$ is a non-negative integer, integrating from 
$0$ to 
$\infty$ in 
$k$ and interchanging integrals yields the system
\begin{equation} 2\eta_R\varOmega\mu g^{(s)}_m = \sum_{n=0}^\infty a_n\left[\frac{2\bar\eta}{{\rm \pi} }\xi_m^{(s)} \varLambda_n^{(s)} + \frac{2\bar\eta}{\rm \pi} M_{nm}^{(s)} + \mu \varDelta_{nm}^{(s)}\right], \end{equation}where
\begin{gather} g^{(s)}_m = \int_0^\infty \, \mathrm{d}k\, k^{{-}3/2-\beta} J_{2m+1/2+\beta}(kR) \sin kR = \frac{{\rm \pi} (R/2)^{\beta+1/2}}{2(\beta+1/2)!} \delta_{n0}, \end{gather}
\begin{gather}\xi^{(s)}_m = \int_0^\infty \, \mathrm{d}k\, k^{{-}3/2-\beta} J_{2m+1/2+\beta}(kR) k\cos kR = 0, \end{gather}
\begin{gather}\varLambda^{(s)}_n = \int_0^\infty \, \mathrm{d}k'~(k')^{{-}1/2-\beta} \tanh k'H\, J_{2n+1/2+\beta}(k'R)\sin k'R, \end{gather}
\begin{gather}M^{(s)}_{mn} = \int_0^\infty \, \mathrm{d}k'~(k')^{{-}2\beta} \tanh k'H\, J_{2n+1/2+\beta}(k'R)J_{2m+1/2+\beta}(k'R) \end{gather}
\begin{align} {\rm \Delta}^{(s)}_{mn} &= \int_0^\infty \, \mathrm{d}k\, k^{{-}1-2\beta} J_{2n+1/2+\beta}(kR)J_{2m+1/2+\beta}(kR) \nonumber\\ &=\frac{\beta(2\beta-1)!(m+n-1/2)!}{4^\beta (\beta+m-n)!(\beta+n-m)!(1/2+2\beta+m+n)!}, \end{align}
and factorials assume their usual definition via the gamma function (formula 8.310.1 of Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007). Note the strong resemblance to (4.36) to (4.41) for the disc. Following the discussion in § 4.2.1, we choose 
$\beta = 1/2$ and numerically evaluate 
$M^{(s)}_{mn}$ and 
$\varLambda ^{(s)}_{n}$ for 
$m$ and 
$n< 20$. The truncated linear system is quickly solved for the coefficients 
$\{a_n\}$. Figure 12 depicts the resulting surface flow field 
$V(x)$ for different values of 
$H/R$, which is found by evaluating equation (C8).
Figure 12. The surface azimuthal velocity 
$V$ due to an infinite strip of half-width 
$R$ as subphase depth 
$H$ is varied. As 
$H\to 0$, a boundary layer of width 
$\bar {\delta } = \sqrt {\bar \eta H/\mu }$ becomes visible and 
$V$ is well approximated by (C17). Parameters: 
$\eta _R/\eta =1.875, \mu R/\eta = 30$.
Finally, we turn to the high friction (
$H\to 0$) case, where 
$V$ satisfies the Brinkman equation
\begin{equation} \bar\eta \frac{\mathrm{d}^2 V}{\mathrm{d}\kern0.7pt x^2}\chi(|x|< R) = \varGamma V, \end{equation}
with 
$\varGamma = \mu /H$. Using the boundary condition (C2a,b), the solution is
\begin{equation} V(x) = \frac{ 2\varOmega \bar{\delta}\eta_R}{\bar\eta}\frac{\sinh (x/\bar{\delta})}{\cosh (R/\bar{\delta})}\chi(|x|< R), \end{equation}
where 
$\bar {\delta } = \sqrt {\bar \eta /\varGamma }$ is the penetration depth of the edge current. The convergence to this solution as 
$H$ is decreased can be seen in figure 12.
