Vortex Flows and Streamline Topology in Curved Biological Membranes
VVortex Flows and Streamline Topology in Curved Biological Membranes
R. Samanta
1, 2 and N. Oppenheimer
1, 2 Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel The Center for Physics and Chemistry of Living Systems, Tel Aviv University (Dated: December 23, 2020)When considering flows in biological membranes, they are usually treated as flat, though moreoften than not, they are curved surfaces, even extremely curved, as in the case of the endoplasmicreticulum. Here, we study the topological effects of curvature on flows in membranes. Focusing ona system of many point vortical defects, we are able to cast the viscous dynamics of the defects interms of a geometric Hamiltonian. In contrast to the planar situation, the flows generate additionaldefects of positive index. For the simpler situation of two vortices, we analytically predict the loca-tion of these stagnation points. At the low curvature limit, the dynamics resemble that of vorticesin an ideal fluid, but considerable deviations occur at high curvatures. The geometric formulationallows us to construct the spatio-temporal evolution of streamline topology of the flows resultingfrom hydrodynamic interactions between the vortices. The streamlines reveal novel dynamical bi-furcations leading to spontaneous defect-pair creation and fusion. Further, we find that membranecurvature mediates defect binding and imparts a global rotation to the many-vortex system, withthe individual vortices still interacting locally.
Contents
I. Introduction II. Setup : Viscous hydrodynamics in curved membranes coupled to external solvents.
III. Streamline Topology for Membrane Vortices
IV. Many vortices V. Membranes of general curvature and other possible extensions VI. Conclusion VII. Acknowledgements VIII. Appendices s l : Poles of the stream function in Legendre basis on the sphere. 33 References a r X i v : . [ phy s i c s . f l u - dyn ] D ec I. INTRODUCTION
We study two dimensional (2D) flows in curved biological membranes arising from the dynamics of rotatingembedded particles. In particular, we explore the spatio-temporal evolution of topological features of such 2D flows.The study is relevant in the context of biological systems featuring vortical defects. In this work we are primarilymotivated by proteins embedded in lipid membranes, in particular the rotating ATP synthase proteins [1], abundantin the endoplasmic reticulum. A detailed knowledge of such flows will shed light on possible mechanisms of mixingin biological viscous fluids [2], as well as serve as a guiding principle to engineer molecular rotors [3], artificiallycontrolled microswimmers and nano-carriers in targetted drug delivery [4], or in wound detection and healing [5].Moreover, the rotating inclusions that we consider in much of our analysis can be realized in experiments e.g. byparamagnetic microscopic particles in a rotating magnetic field [6], birefringent particles rotated by laser tweezers [7]and biological swimmers such as bacteria, Volvox algae and diatoms [8, 9].The dynamics of physical systems in the presence of topological defects and curvature is currently an activearea of research. Topological defects can play a major role in key macroscopic properties of the system — be itin driving phase transitions, creating fluid flow patterns, or the emergence of turbulence. Mostly, defect dynamicsare investigated in planar 2D systems, a few prime example are vortices in superfluids [10, 11], Abrikasov vortexlattices in superconductors [12, 13], and vortex driven Berezinski Kosterlitz Thouless (BKT) transition [14, 15]. Thedynamics of vortices in ideal fluids, including integrability, chaos, and stability analysis is also a topic of intenseresearch [16–21]. The natural world often features motion of defects on curved surfaces. Point vortex dynamics ina spherical geometry may be used as an approximation for air flow in the earth’s atmosphere and oceans [22]. Inrecent years, experimental advances in condensed matter (Bose-Einstein condensates in particular) have also openedup the possibility to explore superfluid vortices in a curved substrate or in optical traps [23, 24]. These systems thusallow a rich interplay between curvature and dynamics of topological defects.In the biological world, curvature and defects feature in the vast majority of living systems e.g. cell membranes,tumor growth, and morphogenesis [25–27]. For planar systems, topological defects play a vital role in BKT-like phasetransitions in active nematics [28]. Such defects lead to turbulent flow patterns even in such highly viscous fluids[29–35]. Recently, motivated by biological examples, defect dynamics in nematic films have been explored on curvedsurfaces as well [36, 37]. Polar active fluids confined to curved surfaces show flocking and topologically protectedsound modes [38].In this work we perform a detailed study of vortical defects in curved biological membranes. The investigation ofin-plane and out-of-plane dynamics of membranes separating two viscous fluids has been a subject of much interest[39–42]. A detailed analysis of different modes in such membranes, including force response, mobility calculations,and many-body interactions was performed [43–49]. In particular, a study of 2D flows and dynamics of rotors wascarried out in detail in [50, 51]. The quasi-2D nature of membranes leads to a new length-scale (here termed theSaffman length) given by the ratio of membrane and solvent viscosities. This length scale acts as a natural cutoff forthe logarithmic divergence of two dimensional flows. Beyond the Saffman length, the external solvent contributes tothe in-plane dynamics, regulating the divergence.A detailed study of biophysical transport applicable to curved membranes was carried out in a few recent works[52–58]. In particular [52, 53] generalized the pioneering works of Saffman and Delbr¨uck [39, 40] for curved surfacesof static geometry. For a spherical membrane, particle mobility was computed. In the limit of high curvature (smallradius), one finds a reduced mobility, with the sphere radius playing the role of the Saffman length. On the otherhand, at low curvature, the Saffman length still continues to regulate the logarithmic divergence. Further, the studyreveals the existence of a zero mode due to curvature, that imparts a global motion to the system. The initial partof our analysis here is a direct follow-up of these works [52, 53].The work presented here outlines the following aspects of 2D flows in curved membranes:1. At the single particle level, we extend the works of Henle and Levine [52, 53] to account for rotational flows,such as those generated by point vortices and torque dipoles, as applicable to ATP synthase proteins.2. We compute the relevant Green’s function in real-space in closed form using Appell Hypergeometric Functions.These functions, arising from the geometry and topology of the membrane prove extremely useful in analysingmany-particle dynamics. For two particles, we analytically predict the location of singularities in the flow fieldsusing this approach.3. For low curvatures, we find a surprising structural similarity between the equations of viscous membrane hydro-dynamics sourced by rotating inclusions and the equations of point vortices in ideal fluids on curved surfaces.This analogy suggests that ideal point vortex models may be a useful tool to gain a basic understanding ofdefect mediated biological turbulent flows observed in the viscous low Reynolds regime, see [35].4. At high curvatures, we find that the dynamics in membranes deviate from the ideal fluid case. The soft modedue to curvature imparts a global rotation to the many-body system, while the rotating inclusions individuallyinteract in a manner similar to local point vortex like interactions.5. We provide explicit formulas for the dynamical equations and flows (Eqs. 18, 19) and the rotation rates (Eq. 25and Eq. 36) in the full parameter space of the biological model.6. We cast the viscous dynamics of the rotating inclusions in the spherical membrane in terms of a geometricHamiltonian, with associated conservation laws. The analytic structure of the Hamiltonian allows us to constructthe spatio-temporal evolution of streamlines resulting from hydrodynamic interactions between the point rotors.In contrast to the planar situation, we find the existence of new vortical defects of positive index (centers) inthe flow-fields where the velocity field vanishes. The number of such new stagnation points is strictly governedby the Euler Characteristic of the surface, consistent with Poincare Index Theorem.7. For many vortices with varying circulations, the temporal evolution of the streamline topologies reveal noveldynamical bifurcations leading to defect-pair fission and fusion. We are able to demonstrate all these effectswith a relatively small number of vortices. Further, we observe that the membrane curvature imparts a globalrotation to the system and can drive the binding of defects with opposite index, similar to activity driven defectbinding and unbinding phenomena observed in 2D nematic fluids [28].The paper is organized as follows : In Section II we present a short review of the basic equations for viscoushydrodynamics on curved membranes. In Section II A, II B and II C, we provide typical examples of the 2D fluid flowsin spherical membranes due to 3 types of sources: a point force, a point torque and a torque dipole. The detailedcalculations are presented in the Appendix VIII A, VIII B, VIII C. In particular, Section II B explores the connectionsto equations arising in vortex dynamics in ideal fluids on curved surfaces. The rotating solutions allow us to construct aHamiltonian description for a system of rotating inclusions embedded in the membrane. This Hamiltonian descriptionis presented in Section III B along with basic equations III A to explore the streamline topology of the hydrodynamicflow fields. These equations are used to construct the spatio-temporal evolution of streamline topology of the in-planeflow fields. We explore the streamline flows for different vortex circulations along with an analytic understanding ofthe associated stagnation points. Next in Section IV we present some interesting scenarios of spontaneous creationof defect pairs and defect fusion that arise in such systems in the chaotic regime of many interacting point rotors.Finally in section V and VI we conclude with possible generalizations.The Appendix contains many details of the calculations and formulas used in the main text. Appendix VIII A, VIII BVIII C describe the full structure of the real space Green’s functions used in our study, while Appendix VIII E discussesthe pole structure of the Green’s function in Legendre basis. Appendix VIII D supplements an analytic investigationof stagnation points and streamline topology carried out in the main text for the situation of two vortices.
II. SETUP : VISCOUS HYDRODYNAMICS IN CURVED MEMBRANES COUPLED TO EXTERNALSOLVENTS.
Let us start by describing the hydrodynamic equations for curved membranes. We use the pioneering works ofSaffman and Delbr¨uck [39, 40] as adapted to a spherical membrane [52, 53]. We approximate the membrane as a two-dimensional viscous fluid surrounded above and below by three-dimensional viscous fluids. We also assume strictlytangential flows within the membrane. In such situations, the appropriate generalization of the Stokes equationsdescribing 2D flows is D α v α = 0 , σ extα = − η d (cid:0) K ( (cid:126)x ) v α + D µ D µ v α (cid:1) + D α p + (cid:16) σ dαz | z → − − σ dαz | z → + (cid:17) , (1)where x represents a general coordinate on the surface, v α is the in-plane 2D fluid velocity ( α runs over surfaceco-ordinates) and η D denotes the viscosity of the 2D membrane fluid. D is the two dimensional covariant derivativewhich generalizes the partial derivative of flat space, K ( (cid:126)x ) is the local Gaussian curvature, p is the local membranepressure, σ D denotes the bulk fluid stress tensor while z denotes a generalized co-ordinate in the normal directionto the surface. The first of the two equations in Eq. 1 ensures incompressibility of the membrane fluid while thesecond equation is a stress balance condition on the membrane surface. The point source embedded in the membraneprovides σ ext . This is balanced by the stress provided by the 2D membrane fluid and the external solvents aboveand below the membrane. In the limit of vanishing curvature, K = 0, one recovers the usual Stokes equations in thepresence of a source σ ext . In curved surfaces, the covariant derivatives fail to commute and this is the origin of thecurvature term K ( (cid:126)x ). For more details, see Appendix VIII A.These equations need to be supplemented by the appropriate Stokes equations for the 3D outer fluid, η ± ∇ v ± = ∇ ± p ± , ∇ · v ± = 0 , (2)where v + ( v − ) is the fluid velocity above (below) the membrane, with similar notation for pressure p ± and viscosities η ± . One can define two length scales given by the ratio of membrane and solvent viscosities λ ± = η D η ± . (3)The curvature introduces a new scale in the problem. In spherical membranes for example, this will be the radius R .The coupling between the 2D membrane flows and the 3D external solvents is mediated via the no-slip boundarycondition and by the stress balance on the membrane surface (the last term in Eq. 1). For membranes of arbitraryshape, these equations can be solved numerically, e.g. see [57]. However, in the simpler situation where the curvatureis constant, one can analytically extract flows, which we describe next.The incompressibility requirement D α v α = 0 allows us to express the flow field in terms of a stream function asfollows: v α ( x ) = (cid:15) αγ D γ φ ( (cid:126)x ) , (4)where (cid:15) αβ is the antisymmetric Levi-Civita symbol. One thus needs to solve for φ ( (cid:126)x ), given a point source σ ext ,taking into account the membrane curvature and boundary conditions. As shown in detail in Refs. [52, 53] (seealso our appendix VIII A, VIII B and VIII C), such a response calculation is conceptually simple. One needs toinvert the curved surface Laplace operator in the presence of the curvature and traction terms . For a non-trivialspatially varying curvature K ( x ), a Fourier decomposition can be done numerically but requires knowledge of thespectrum. For surfaces of constant curvature, the spectrum is often known. One can use the known eigenfunctionsto perform the inversion in Fourier Space. For example, for a sphere, one decomposes the above equations in thebasis of spherical harmonics, taking into account the stick and stress boundary conditions. Henle and Levine [52, 53]express the final solution of the stream function in terms of such eigenmodes. Using analytically continued AppellHypergeometric Functions, we are able to perform the inverse Fourier transform and find closed-form expressions forthe stream function in real space. We present a detailed description in Appendix VIII A, VIII B and VIII C for eachof the sources: point force, point torque and a torque dipole respectively.Before proceeding, let us briefly mention some general topological constraints that the flow fields on the sphericalmembrane must satisfy. First, the hairy ball theorem implies that flow fields on the spherical membrane must featurestagnation points where the velocity field vanishes. Second, each of the singular points of the flow field can be assignedan index which keeps track of the winding of the flow field around the core of the singularity. The sum of these indices must equate to 2, the Euler Characteristic of the sphere (Poincare Index Theorem). In all the examples involvingsingle inclusions that we are about to study, we will observe these topological features in the flow fields.We now describe the flow fields resulting from a point force, point torque, and torque-dipole, one by one: A. Velocity field due to a point force on the spherical membrane.
The main body of this work concerns vortices in a membrane, but for completeness and consistency, in this sectionwe reproduce the results of Henle and Levine for a point force acting on a spherical membrane. The velocity field at an Without the traction contribution, the curvature term gives rise to a zero mode. However, in the presence of traction, the zero mode isremoved and the operator can be inverted without issues. We will present a more general criterion for the situation of many embedded particles, see Eq. 21 in Section III B
FIG. 1: Streamline plot of the velocity field in the low (top) and high (bottom) curvature regimes, in response to a forcelocalized on the spherical membrane. On the left, the flow field is shown in a θ, φ chart while on the right, the flow field iswrapped on a spherical membrane. Location of the force is marked in red. Note the creation of two vortical defects around theforce. In the top row (low curvature regime) we show the flow field for a point force localized at θ = 1 .
5. In the bottom row(high curvature regime) for a point force localized at θ = 0, the vortices migrate to the equatorial regions. arbitrary point ( θ, φ ) on the sphere due to a point force localized on the membrane surface at ( θ , φ ) is summarizedby an Oseen tensor on S given by v = G ( θ, θ , φ, φ ) F , where the different components of the Green’s function canbe expressed in terms of double derivatives of a function SG θθ = csc θ csc θ πη D ∂ φ ∂ φ S, G θφ = − csc θ πη D ∂ φ ∂ θ S, G φθ = − csc θ πη D ∂ θ ∂ φ S, G φφ = πη D ∂ θ ∂ θ S, (5)and the function S is defined in the basis of Legendre Polynomials, S := (cid:80) ∞ l =1 (2 l +1) s l l ( l +1) P l (cos γ ), where s l = l ( l + 1) − Rλ − ( l −
1) + Rλ + ( l + 2), and cos γ is the cosine of the geodesic angle between the source at ( θ , φ ) and response at( θ, φ ) cos γ = sin θ sin θ cos( φ − φ ) + cos θ cos θ . (6)The function S thus varies with the geodesic angle and the physical parameters, namely the sphere radius, and themembrane and solvent viscosities. The full structure of the function S in real space is presented in Appendix VIII Ain terms of analytically continued Appell Hypergeometric Functions. With the knowledge of the real space Green’sfunction in the full parameter space at hand, we now plot the resulting flows due to a Stokeslet (point force) localizedon the north pole of the sphere. In these plots we have chosen η + = η − = η d in Eq. 3. Thus, one can compare theradius of the sphere R with respect to the unique Saffman Length λ = η d η d . There are two distinct regimes R > λ (low curvature) or
R < λ (high curvature).In the low curvature regime, the velocity field exhibits a dipole-like structure around the point of application ofthe force. The dipole has a topological index +2 which agrees with the Euler Characteristic of the sphere. As thecurvature is increased, the dipole structure breaks into two +1 vortices which migrate away to diametrically oppositepoints. These features were predicted first in [52, 53] and generalized to lipid bilayers with slip velocity in [56].We observe that our real-space Green’s function (Appendix VIII A) also reproduces these effects. This provides aconsistency check of our summation procedure explained in Appendix VIII A.
B. Velocity field due to a point torque on the spherical membrane
The velocity field at a point ( θ, φ ) produced by a point torque of circulation τ localized at ( θ , φ ) on the spherecan be expressed as v = τη D [ ∇ S ⊥ ] ψ , (7)where [ ∇ S ⊥ ] = (cid:16) ˆ θ R sin θ ∂ φ − ˆ φ R ∂ θ (cid:17) and ψ represents the dimensionless stream function. In terms of Legendre modes, ψ is given by ψ [ θ, φ, θ , φ ] = (cid:88) l (2 l + 1)4 π s l P l (cos γ ) . (8)where cos γ = sin θ sin θ cos( φ − φ ) + cos θ cos θ is the cosine of the geodesic angle between the source and responselocations. The real space representation of ψ in the full parameter space is presented in Appendix VIII B. There aretwo different representations of the stream function, one valid at low curvatures (Eq. 104) and the other valid at highcurvature (Eq. 105).We now focus on the associated topology of the flow-field due to the rotating inclusion localized at the north pole,see Fig. 2. We find that a new vortical defect of positive index (center) develops at the south pole to make the totalindex +2 as required by the topology of the sphere.Let us briefly comment here on an interesting connection between the point vortex flows we study in this viscous FIG. 2: Flow fields due to a vortex localized at north pole. Top row shows the low-curvature case (Eq. 104) and bottom rowshows the high-curvature one (Eq. 105). On the left, the magnitude of the azimuthal velocity v φ is shown in the θ direction.In the high curvature case the flow develops a local maximum due to the effect of the global rotation, shown as a gray line (seefurther discussion in Sec. II D). Right panels show the flow field wrapped on a spherical membrane. set-up and the ideal point vortex problem on curved surfaces. For simplicity, let us consider the equation of stressbalance for viscous hydrodynamics of the membrane fluid (with no external solvent) in the presence of a point rotorof unit strength. η D (cid:0) K ( x ) + ∆ (cid:1) v α = (cid:15) αβ D β δ ( θ − θ ) δ ( φ − φ ) R sin θ (cid:124) (cid:123)(cid:122) (cid:125) ω , (9)where we have defined ω as suggested in the above equation.In terms of the scalar stream function defined as v α = (cid:15) αβ D β φ , the above equation reads η D (2 K ( x ) (cid:15) αγ D γ φ + (cid:15) αγ D γ ∆ φ ) = (cid:15) αγ D γ ω (10)In the limit of low curvature where K ( x ) can be ignored, we are left with η D ∆ φ = ω. (11)This equation is identical to the that of a point vortex in an ideal fluid on a curved surface, where ∆ is the surfaceLaplace Beltrami Operator. Due to this equivalence, we expect that at low curvatures (and hence on the plane inparticular) the response to rotating inclusions in a viscous fluid is similar to point vortex flows in an ideal fluid. Thesituation in the biological model we consider here departs from the ideal vortex problem once the curvature termbecomes important. Indeed, we find that at low curvatures where the radius of the sphere is much larger comparedto the Saffman length, the flows resemble those of ideal point vortex problem on the sphere. However, as the radiusdecreases, the curvature contributes to a zero mode that imparts a global rotation to the many body system, whilethe rotating inclusions individually still continue to interact in a manner similar to local point-vortex like interactions.At high curvatures, due to the effect of the global rotation, the flow is no longer monotonically decreasing, see Fig. 2. C. Velocity field due to a torque dipole on the spherical membrane
Rotating inclusions that arise in biological examples feature no external torque e.g. rotor proteins such as ATPsynthase. To account for this additional structure, we also construct a model of counter-rotating torque-dipole [50, 51].On a spherical membrane the solution is constructed in Appendix VIII C. The dynamics are very similar, mainly givingrise to a faster spatial decay. The solutions are identical in terms of their topology, and thus, in the rest of this paperwe will focus on the solutions due to point rotors, but the Appendix outlines the results for torque-dipoles as well.
D. Emergence of Global Rotation at high curvature
Having understood the flow fields due to simple source terms on the spherical membrane, we now illustrate theemergence of the global rotation in the high curvature regime for all three situations. This global rotation was firstreported in [52, 53] for a Stokeslet. In all the situations we studied, the dimensionless stream function on the sphericalmembrane has the generic structure in the basis of Legendre polynomials: ψ [ θ, φ, θ , φ ] = (cid:88) l f l π s l g l P l (cos γ ) , (12)where f l and g l are polynomials in Legendre modes denoted by l and s l = l ( l + 1) − Rλ − ( l −
1) + Rλ + ( l + 2). Thegeodesic angle between the source and response locations is denoted by γ . One can understand the emergence ofthe global rotation from the common denominator s l arising in all three situations. For simplicity, let us assume theSaffman lengths associated with external and internal solvents to be the same and denote it by λ . The classification Please see Appendix VIII B where we show that the membrane pressure vanishes in this situation. For the spherical membrane of constant curvature, the complete equation reads w = η D ( R φ + ∆ φ ) of high and low curvature regimes is then simply determined by the ratio λ/R .In the high curvature limit of λ/R (cid:29) s l ∼ l ( l + 1) − . (13)The zero mode l = 1 dominates the Legendre sum in the stream function of Eq. 12 and generates the global rotation.To see this, let us consider the situation where we place a rotating inclusion at the north pole. The stream functionin Eq. 12 leads to a velocity field with flows only in the azimuthal direction given by v φ = τ πη d R (cid:88) l (2 l + 1) P l (cos θ ) s l , (14)where P l (cos θ ) denotes the associated Legendre function of first order. In the limit of high curvature, the zero mode l = 1 dominates the sum and we get v φ = τ sin θ πR η + , (15)which corresponds to a global rotation of the flow withΩ = τ πR η + . (16)Let us note that the global rotation rate is purely regulated by the external solvent.One can also use the asymmetry of external and internal solvents at a fixed radius R to generate this global rotation.This happens when λ + (cid:29) λ − . s l = l ( l + 1) − Rλ − ( l −
1) + Rλ + ( l + 2) ∼ l ( l + 1) − Rλ − ( l − . (17)Let us note that l = 1 continues to be the zero mode in this situation leading to the global rotation. This is unlikethe opposite limit when λ + (cid:28) λ − .It is interesting to note how our real space representation of the stream function Eq. 7 in terms of Appell Hypergeo-metric functions capture the global rotation. As shown in AppendixVIII B, VIII C, depending on the roots of s l = 0( the roots are analyzed in Appendix VIII E) the Legendre sum defined in Eq. 8 leads to two different representationsof the stream function in real space. One of these representations Eq. 104 is valid in the regime of low curvature,while the other Eq. 105 is valid in the high curvature regime. The appropriate stream function in the high curvatureregime indeed shows a dominance of the global rotation term, see Fig. 13 in AppendixVIII B. III. STREAMLINE TOPOLOGY FOR MEMBRANE VORTICES
In this section we set up the equations needed to explore the streamline topology of flows on the biological membranedue to rotating inclusions.
A. Dynamic Equations for an ensemble of membrane vortices
Let us consider N rotating inclusions embedded in a spherical membrane with viscosity η d , surrounded by externalsolvents with viscosities η ± . The evolution equations for purely hydrodynamic interactions between the vortices aregiven by ˙ θ i = 1 η d R N (cid:88) j (cid:54) = i τ j sin θ i ψ (cid:48) [ γ ij ] ∂ φ i [ γ ij ] , ˙ φ i = − η d R N (cid:88) j (cid:54) = i τ j sin θ i ψ (cid:48) [ γ ij ] ∂ θ i [ γ ij ] . (18)where ψ (cid:48) is the derivative of the stream function defined in Eq. 7 for point rotors, with the explicit structure in termsof Appell Hypergeometric Functions given in Appendix VIII B by Eqs. 104 and 105 for the low and high curvaturesrespectively (see Eq. 127, Eq. 128 in Appendix VIII C for the corresponding expressions for torque dipoles). Let usnote that the stream function has two different representations that are dictated by curvature. While performing thedynamical simulations, one needs to insert the appropriate representation of ψ into Eq. 18. Finally, γ ij is defined inEq. 6.Let us add some comments on the absence of a self-drive term in Eq. 18. This can be argued from symmetryconsiderations. Due to spherical symmetry, there is no preferred direction and a single vortex does not move. Thisargument works for the planar situation as well. However, in generic surfaces with no (or restricted) symmetry, oneneeds to treat the self interaction term with a proper regularization procedure. This will in general lead to a self driveterm. For the spherical membrane however, such a regularization leads to a constant (due to symmetry) and has noeffect on the dynamics.We will be interested in the flow fields resulting from the above dynamics as well. For this purpose, one constructsthe hydrodynamic velocity field at any given point p via superposition (taking into account the contributions from allpoint rotors). v θ p = 1 η d R N (cid:88) j τ j sin θ p ψ (cid:48) [ γ pj ] ∂ φ p [ γ pj ] ,v φ p = − η d R N (cid:88) j τ j ψ (cid:48) [ γ pj ] ∂ θ p [ γ pj ] . (19)where γ pj denotes the geodesic angle between ( θ p , φ p ) and ( θ j , φ j ) and now the sum runs over all vortices. B. Hamiltonian Formulation
The dynamical equations (Eq. 18) can be cast in terms of a geometric Hamiltonian. In terms of canonicalcoordinates, Q i = (cid:112) | τ i | φ i , P i = (cid:112) | τ i | cos θ i , Eq. 18 can be re-written as ˙ Q i = ∂ P i H, ˙ P i = − ∂ Q i H,H = 1 R η D (cid:88) i The theorem forbids the existence of a nowhere vanishing vector field on the sphere,i.e. there is at least one point where the flow field vanishes on the sphere. Let us mention that the so defined canonical momentum P in our Hamilton description is thus a function of the coordinate θ . ThePoisson bracket relation [ Q i , P j ] = δ ij τ i thus implies that position coordinates θ, φ do not commute in this geometric formulation. Thismay seem surprising at first sight, however this is not uncommon. For example, consider motion of electrons in a magnetic field (say inthe ˆ z direction), where the canonical momentum is given by p x = m ˙ x − eBmc y . The canonical Poisson brackets — [ x, p x ] = 1 in the large B limit reads [ x, y ] ∼ − mceB (cid:54) = 0. Indeed this has been a major motivation to explore if vortices in superfluids can form quantum Hallstates [59–61]. Poincare Index Theorem. The topological defects in the flow field on a closed, compact surface can beassigned an index corresponding to the winding of the field around the singular core. The theorem implies thatthe sum of the indices over all singularities is equal to the Euler Characteristic of the surface. For a sphericalmembrane (or membranes deformable to a sphere), the Euler characteristic is 2. As we will see later, thetopology of the spherical membrane leads to creation of N c new centers (vortex defects of index +1 where thevelocity field vanishes, i.e. a stagnation point of the flow-field), a phenomenon not observed on the plane [18].Such centers have unit positive index. In addition, one has N S saddles of unit negative index (an anti-vortex).Further, each of the N rotating inclusions that we consider contribute a positive index of +1. The index countingthus demands N + N c − N s = 2 . (21)The circulation τ of the vortices is independent of the index. In particular, both positive and negativecirculations have the same index +1.3. Integrability and the Liouville-Arnold theorem. A Hamiltonian system with 2 N dimensional phase-spaceis integrable if there exist N independent integrals of motion which are all mutually Poisson commuting, i.e.they are in mutual involution. The Hamiltonian we constructed in Eq. 20 has the same set of symmetries asideal point vortices on a sphere which have 3 mutually commuting conserved quantities. Thus, the system losesintegrability for N ≥ N = 4 situation is integrable if the total circulation of vortices is zero [19].We are now ready to discuss the dynamics of the vortices and spatio-temporal evolution of vortical defects in theflows within the biological membrane. The basic methodology we adopt here is very simple: we first simulate thedynamics of the vortices as given in Eq. 18 and feed the dynamic locations of the vortices into Eq. 19 to get thehydrodynamic vector field ( v θ , v φ ). Using Mathematica [62], we next plot the streamlines associated with this flowfield in the ( θ, φ ) chart. In order to understand the dynamics and flows better, we restrict our discussions in thissection to two vortices and discuss the many vortex situation in the next section. We separate the discussion of twovortices into regimes of low and high curvature. In each regime, we first consider the simple situation of two vorticeswith same (opposite) circulations τ , discuss the dynamics and the associated streamline topology of the flow field. Inappropriate places, we comment on important distinctions from the planar dynamics and flows. From our discussionsnear Eq. 11, we expect the dynamics and flows to be similar to ideal vortices in the low curvature regime. The flowdeparts from ideal vortices once curvature becomes important. C. Low Curvature Regime A single vortex does not move due to spherical symmetry. The dynamics becomes interesting once we have twovortices or more. To understand the dynamics of two vortices better, we first note that the Hamiltonian Eq. 20preserves the chord distance C = | (cid:126)X − (cid:126)X | between the vortices and they orbit around each other. The rotationrate can be easily estimated by converting the dynamical equations for two vortices (Eq. 18) into Cartesian form: ddt (cid:126)X = τ Rη D ψ (cid:48) [ γ [ C ]]sin[ γ [ C ]] (cid:126)X × (cid:126)X R ,ddt (cid:126)X = τ Rη D ψ (cid:48) [ γ [ C ]]sin[ γ [ C ]] (cid:126)X × (cid:126)X R . (22)where γ [ C ] is related to the chord distance C = | (cid:126)X − (cid:126)X | via γ [ C ] = 2 arcsin (cid:104) C R (cid:105) .The center of vorticity vector (cid:126)M in this situation is given by (cid:126)M = τ (cid:126)X + τ (cid:126)X τ , (23) Drawing the streamlines on the ( θ, φ ) chart is convenient to keep track of the evolution of vortical defects over the entire sphericaldomain. However one can equally well wrap the flow field on the spherical membrane. τ = τ + τ is the total circulation. It is easy to see that Eq. 22 can be written using the above vector (cid:126)M as ddt (cid:126)X = 1 Rη D ψ (cid:48) [ γ [ C ]]sin[ γ [ C ]] (cid:126)X × τ (cid:126)MR ,ddt (cid:126)X = 1 Rη D ψ (cid:48) [ γ [ C ]]sin[ γ [ C ]] (cid:126)X × τ (cid:126)MR . (24)From Eq. 24, we can read the rotation rate at low curvature ω vorticesLC asΩ vorticesLC = τ | (cid:126)M | R η D ψ (cid:48) [ γ [ C ]]sin[ γ [ C ]] . (25)Let us note that the same formula also holds for the torque-dipole case with appropriate ψ given in Appendix VIII C.We now elaborate on the dynamics for two vortices in Fig. 3. Same Circulation : In this situation, the vortices orbit each other with an angular frequency ω vorticesLC given byEq. 25. In Fig. 3 we show an example of an orbit (the trajectory of the vortices) and the associated flow fields attwo instants of time. One observes that in addition to the original centers created by the vortices themselves, theflow field exhibits a new center and a saddle where the velocity vanishes. Overall the index adds up to 2 , consistentwith Poincare Index Theorem, as expected from the topology of the spherical membrane. In Fig. 5 , we show howthe location of stagnation points change as we change the distance between the vortices (of same circulation). Inparticular, let us note that when the distance is π then there is a continuum of stagnation points formed along themid-line between them (left most figure in lower panel of Fig. 5). Let us note that the saddle is always formed betweenthe vortices, as can be seen also from continuity of the flow. Opposite Circulation : In this case the vortices move together, such that the perpendicular bisector of the linejoining them follows a geodesic. Interestingly, the flow fields exhibit no new center or saddle. The Poincare IndexTheorem is still satisfied, since the index contribution from the two vortices of opposite circulation is 2. Let us mentionthat Kimura [17] predicted that a vortex dipole (vortices of equal and opposite strength placed close to each other)traces a geodesic in all closed Riemann Surfaces deformable to the sphere. For our spherical membrane, we indeedfind this property holds true, see Fig. 4. Calculation of Stagnation Points. We now proceed to a calculation of the locations of the stagnation points.As shown in Appendix VIII D, one can project the dynamical equations Eq. 18 via stereographic projection on theplane. Using complex coordinates to denote the locations of the vortices, the equation of a tracer particle in thepresence of the vortices can be cast in complex notation ddt ¯ z p = iη D R (1 + | z p | ) ∂ z p H p , (26)where z p denotes the location of the tracer particle, H p is defined using the streamfunction ψ (with the structurepresented in appendix VIII B , see appendix VIII C for the torque dipole case) H p = N (cid:88) j τ j ψ [ γ pj ] , (27)where the geodesic distance in complex notation is given by γ pj = arccos (cid:32) (1 − | z p | )(1 − | z j | ) + 4 Re [ z p ¯ z j ](1 + | z p | )(1 + | z j | ) (cid:33) . (28)It follows that solving for the stagnation points amounts to finding solutions to ddt ¯ z p = 0 . (29)As shown in Appendix VIII D, this amounts to solving an equation of the general form N (cid:88) j τ j F [ z p , z j ] G [ z p , z j ] = 0 , (30)2 FIG. 3: Streamline plot for two vortices at the low curvature regime, both released at the equator with one at φ = 0 . φ = 1 . where the factors F and G arise from the derivative of the stream function i.e. ∂ z p ψ = ∂ ψ ∂ cos γ ∂ z p cos γ := F × G . Letus note that although F is dependent on the choice of parameters, the factor G is purely geometric. To proceed further,we need to compute F from the appropriate stream function ψ . The full structure of the stream function in AppendixVIII B makes the analysis somewhat complicated, however one can choose a set of parameters for the model to simplifythe stream function. As explained in Appendix VIII D, for a particular choice η D = 3 / , η − = 1 , η + = 2 , R = 1 forwhich λ/R = 1 / 2, the stream function is given a relatively simple expression (Eq. 141). Specializing to the case oftwo vortices on the spherical membrane, let their positions in the complex plane be denoted by z and z . Since thestagnation points are always constrained to lie on the great circle joining the two locations [19], we can essentially3 FIG. 5: Location of staganation points for various relative distances between two vortices of the same strength. Plots are shownfor ∆ φ = 2 , , π in the low curvature regime in the θ − φ plane. A saddle stagnation point (anti-vortex) is shown between thetwo positive vortex defects. A continuum of stagnation points occur at φ = π shown in the rightmost figure. map the dynamics to the unit circle on the complex plane . Using polar representation z = re i ˜ θ , we choose withoutloss of generality z = 1 , z = e iφ , τ = 1 , τ = τ. (31)Plugging Eq. 31 into Eq .30 , we convert it to an effective two parameter problem where the stagnation point z p = e iθ p has to be solved as a function of the relative circulation τ and the location of the second vortex parametrized by φ ,from the equation (see Appendix VIII D for a complete derivation) I = f [cos θ p ] 14 (1 − e − iθ p ) + τ f [cos( θ p − φ )] 12 (cid:16) e − iφ − cos( θ p − φ ) e − iθ p (cid:17) = 0 , (32)where f = 10 − √ − x + x (cid:16) − √ − x + 15( − √ − x ) x + 30 x ( x − 1) arcCoth(1 + √ − x ) (cid:17) π ( x − . (33)We can now systematically search for the location of the new stagnation points θ p as a function of the relative vortexcirculation τ and relative distance between the two original vortices φ . As an example we plot the results for φ = 2 . τ = 1. We plot the real and imaginary parts of I in Eq. 32 as a function of θ . The common zeros of the plots arein good agreement with Fig. 6. FIG. 6: Zeros of the Real and Imaginary parts of the L.H.S. of Eq. 32 as a function of θ p . The new stagnation points arelocated at around the common zeroes i.e. θ p = 1 . , . 5, marked with a red cross in the streamline plot on the right. Via stereographic projection, the azimuthal angle on the spherical membrane maps to the polar angle on the plane. D. High Curvature Regime As we saw in Sec. II, the curvature term in Eq. 1 imparts a global rotation to the system. For a single vortex, wepredicted this rate in Eq. 16. In the generic situation of more than one interacting vortices with varied circulations,one can proceed as follows:In the high curvature regime, we saw in section II D that the l = 1 term dominates the Legendre sum in Eq. 12. Thisterm, as we saw in earlier sections, leads to a global rotation. The global rotation rate for a system of many rotatinginclusions be easily extracted by noting that for R (cid:28) λ , as far as global effects are concerned, one can ignore the localhydrodynamic interactions and keep only the global term in the stream function appearing in the dynamical equationEq. 18 i.e. ψ (cid:48) [ C ij ]sin[ γ [ C ij ]] ∼ η D πRη + ∀ ( i, j ) . (34)Using this approximation in the high curvature regime gives ddt (cid:126)X i = (cid:88) j (cid:54) = i τ j Rη D ψ (cid:48) [ C ij ]sin[ γ [ C ij ]] (cid:126)X i × (cid:126)X j R ∼ (cid:88) j τ j Rη D η D πRη + (cid:126)X i × (cid:126)X j R = (cid:126)X i × (cid:80) j τ j (cid:126)X j πR η + . (35)Thus the global rotation rate is ω NvorticesHC ∼ τ | (cid:126)M | πR η + , (36)where (cid:126)M = (cid:80) Nj τ j (cid:126)X j τ denotes the conserved center of vorticity vector, τ is the total circulation τ = (cid:80) j τ j . Let usnote that the above rate agrees with the one found for a single vortex Eq. 16. The same formula holds for the caseof torque-dipoles, with an additional suppression factor d/R coming from the stream function Eq. 124 .The essential features for two vortices of same (opposite) circulation are detailed below. Same Circulation : The dynamics of two vortices is similar to the low curvature case, with the rate given by Eq. 25,but with the appropriate stream function for the high curvature regime Eq. 105. In the limit of high curvature, thisrate can be approximated by our estimate Eq. 36. In addition, the flow field develops two new centers due to theglobal rotation. The two vortices orbit around one of the global centers. Thus, compared to the low curvature regime,there is now an extra center and saddle appearing in the high curvature regime. The location of these new globalcenters is universal and will be calculated soon. As curvature is increased, the saddles (anti-vortex) move towards theoriginal vortices. This is reminiscent of a binding event of activity driven defects in nematic fluids [28]. For membranevortices, the binding is mediated by curvature versus the nematic case where it is driven by elasticity. Figure 8)shows the new stagnation points. Zoomed images show that curvature drives the newly formed saddles (of index -1)closer to the original vortices (index +1). Opposite Circulation : Here as well the global rotation creates two new centers. Unlike the case of same circulation,the two vortices now orbit around different global centers. Note that the new saddles are formed between the originalvortices and the new global centers, as required by continuity of vector fields. The location of these new global defectsis again universal and independent of the details of the model (see below). Here curvature also drives the newlyformed saddles towards the original vortices, each of which now orbits a different global center. Location of global centers : In the regime of high curvature, the stagnation points are still given by Eq. 29, andEq. 30 as in the low-curvature regime, only the streamfunction used is the one appropriate for high curvature, seeEq. 105 in Appendix VIII B. However, unlike the low curvature regime, there is no simple choice of parameters whichsimplifies the functional form of the stream function. One can still determine the location of the global defects, sincein the limit where the stream function is dominated purely by the global term, one can approximate F in Eq. 30 asfollows : F = d ψ d (cos γ ) ∼ − η D πRη + , (37)while G is a purely geometric factor, same as the low curvature regime (see Appendix VIII D ). G [ z p , z j ] = (1 − | z j | )( − z p ) + 4 (cid:16) (1 + | z p | ) ¯ z j − Re [ z p ¯ z j ]¯ z p (cid:17) (1 + | z j | )(1 + | z p | ) . (38)5 FIG. 7: Streamline plot for two vortices, both at the equator in the high curvature regime, top row with same strength(positioned at φ = 0 . φ = 3 . 3) and bottom row with opposite strength (positioned at φ = 0 . φ = 1 . t = 0 for two vortices, one at θ = π/ , φ = 0 . θ = π/ , φ = 3 . R/λ ∼ − for vortices of same circulation. On the left, we see the appearance of two global centers as expected (maintext). The original vortices are orbiting around the global center formed on the lower half of the leftmost figure. The newlyformed saddles are tightly bound to the original vortices, as shown in the zoomed images on the right, showing each of thevortices separately near φ = 0 . φ = 3 . Using these F and G , the equation for stagnation points z p , given by Eq. 30 simplifies considerably and is purelydetermined by the geometric function G . G ( z p , 1) + τ G ( z p , e iφ ) ∼ ⇒ 14 (1 − ¯ z p ) + τ (cid:16) e − iφ − Re [ z p e − iφ ]¯ z p (cid:17) ∼ z p = e iθ p in the above we get14 (1 − e − iθ ) + τ (cid:16) e − iφ − cos( θ − φ ) e − iθ (cid:17) = 0 , (40)which has the following solutions θ p = ± arccos (cid:20) ± τ cos φ √ τ +2 τ cos φ (cid:21) . (41)For example, for τ = − φ = 1, this yields θ p = 2 . , . , − . , − . θ p = 2 . θ p = − . . FIG. 9: For vortex locations φ = 0 and φ = 1 with opposite circulation, the global centers at the high curvature limit appearat θ p = 2 and θ p = 5 . IV. MANY VORTICES In this section we briefly discuss the situation of many rotating inclusions with varied circulations, with the dynamicsand flow fields described by Eq. 18 and Eq. 19 respectively. For many vortices the dynamics is ergodic and it isin general difficult to perform an analytic investigation. However, one can still simulate the dynamics, numericallysolving Eq. 18 and using the solution to keep track of spatio-temporal evolution of vortical defects via Eq. 19. Onecan build some general intuition in specific situations, as explained in Fig. 10. For example, as shown in top row ofFig. 10, in the low curvature regime, 12 closely spaced vortices of the same circulation tend to rotate together as asingle effective center. Together with the creation of an isolated center, the flow fields furnish a coarse grained versionof Poincare Index Theorem. This also follows from the fact that the symmetries of the Hamiltonian we constructed(Eq. 20) preserve the second moment M = (cid:88) i (cid:54) = j τ i τ j C ij , (42)where C ij denotes the chord distance between the vortices. In this situation since all circulations are same, thisimplies that the vortices will remain geometrically confined within a region of the membrane.With alternating circulations (second row), this is no longer the case, with the 12 centers breaking into smallergroups and spreading across the whole membrane, while still conserving M . In the third and fourth rows of Fig. 10 we Integrability is lost beyond N = 3 vortices. For N = 4 vortices, the system is still integrable if the sum of the circulations vanish. FIG. 10: Dynamics and flows for 12 vortices. The first two rows are at low curvature (large radius): with same (row 1) andalternating (row 2) circulation strengths. Crosses (points) mark clockwise (counter-clockwise) circulation. The last two rowsare at high curvature (small radius): with same (row 3) and alternating (row 4) circulation strengths. In all the rows, thevortices have the same initial locations centred in a small square. The left column shows the trajectories traced with time, themid and right column show the streamlines at initial and final time respectively. Note how in the case of alternating circulationsat the low curvature regime the dynamics span the entire spherical domain, whereas for same circulation vortices the dynamicsstay bounded as resulting from the conservation of the second moment. Eq. 42. Pair creation and fusion Typically, collapse of vortices happens under very special initial conditions [19].As we have seen, unlike the situation on the plane, the topology of the spherical membrane generically leadsto the creation of new vortical defects in the flow fields. This creates the possibility for a spontaneous creationof vortical defect pairs as well as fusion events on the spherical membrane. We observe that this is indeed thecase. The dynamics of vortices drives the system from one configuration of defects to another, with a differentnumber of defects (still satisfying the Poincare Index theorem before and after the bifurcation). We are able todemonstrate these effects with a small number of vortices. Fig. 11 shows a temporal evolution, exhibiting spon-taneous creation and subsequent disappearance of a pair of vortical defects of opposite index ( vortex anti-vortex pair).Within the biological context, one may also incorporate the finite size of the rotating inclusions by introducinga soft repulsion between vortices, in addition to hydrodynamic interactions that we have considered so far. Thismakes the dynamics and flow fields very interesting, see Fig. 12 where we demonstrate a fusion event betweenan original membrane rotor and a newly created defect arising from the spherical topology. The fusion happensvia a bridging saddle of negative index. We expect the number of such events to rapidly proliferate in the sit-uation of large number of inclusions. The many rotor system will be explored in more detail in an upcoming work [63]. FIG. 11: Left to right : an example of dynamical bifurcation for vortex strengths (1,1,-1) where a defect pair is createdspontaneously in the third figure and then disappears as the dipole returns to its original position. Flow fields are shown in( θ, φ ) chart.FIG. 12: Dynamical fusion with five rotors. Black arrows mark the two centers which in the last frame combine into one. Flowfields are shown in ( θ, φ ) chart. V. MEMBRANES OF GENERAL CURVATURE AND OTHER POSSIBLE EXTENSIONS It is straightforward (although computationally intensive) to generalize the calculations presented here to arbitrarilycurved static geometries. Let us denote the coordinates in 3D ambient space by x α and arbitrary coordinates on themembrane by y a . A natural choice for the section of the ambient frame bundle for arbitrarily curved geometries(embedded in R ) is the set ˆ n, (cid:126)e , (cid:126)e where ˆ n = ∇ φ |∇ φ | , e αa = ∂x α ∂y a (43)where the surface is described by the equation φ = 0, a runs over 1 , α runs over all three coordinates on the ambient space. It is a simple matter to repeat the calculations with the normalvector ˆ n now playing the role of the radial vector in the boundary conditions. For a given embedding of the surface,9one first decomposes Lamb’s solution [67] for the 3D fluid (which we constructed in the spherical co-ordinate basis,see Appendix VIII A) in the above new basis constructed out of the hypersurface embedding. Then we apply thestick and stress boundary conditions in this basis . The in-plane membrane velocity field can be extracted by solvingthe spectral eigenfunctions ψ ( y, s ) and eigenvalues λ s from the corresponding equations of membrane hydrodynamics(see Appendix VIII A for a derivation) written in the basis (cid:126)e , (cid:126)e (defined in Eq. 43): η D (cid:0) K ( y ) D a ψ ( y, s ) + D a ∆ ψ ( y, s ) (cid:1) = λ s D a ψ ( y, s ) , (44)for an arbitrary local curvature K ( y ) and where D is associated with the pullback metric on the hypersurface h ab = g αβ e αa e βb . (45)For simple geometries, the spectrum is known and one can derive analytic expressions for the stream functions. Itwill be interesting to see, for example, how the anisotropies of an ellipsoid or of the negative curvature of hyperbolicspaces affect the streamline topologies. We also plan to investigate the dynamics and streamline topology of flowfields in fluctuating biological membranes with bending rigidity and surface tension, see for example [55] for anexperimental perspective. One can also use the flows we constructed in this paper to study the optimal time for naviga-tion of microswimmers in spherical membranes, see [4]. We plan to report on these investigations in subsequent works. VI. CONCLUSION To summarize, in this work we explored in detail the topological aspects of 2D flows resulting from the dynamicsof inclusions embedded in curved biological membranes. To get a first estimate of such flows, we considered a wellknown model [39, 40] adapted to spherical geometry, in the presence of external solvents. In the examples of flows dueto point sources (point force, torque and torque-dipole), we presented new closed form expressions for the respectiveGreen’s function in real space, using analytically continued Appell Hypergeometric Functions. We investigated thetopological features of the flow fields in some detail. Such solutions can be used to model hydrodynamic interactionsof proteins embedded in positively curved surfaces. The topology of the membrane creates many additional vorticaldefects. For the simpler situation of two vortices, we were able to analytically predict the location of such defectswith reasonable precision. The point rotor solutions allowed us to construct a many-body geometric Hamiltonian thatgenerates the dynamics of the vortices on the spherical membrane. We studied the spatio-temporal evolution of defectmediated 2D flows in the spherical membrane. In particular, we found that at low curvature, the flows generated bythe rotating inclusions are similar to flows generated by vortices in an ideal fluid . High curvature imparts a globalrotation to the many-body system, with the individual vortices interacting locally. Already in this simple model atlow Reynolds, we saw surprisingly rich dynamics and flows mediated by the curvature and topology of the sphericalmembrane. The spatio-temporal evolution of streamlines revealed spontaneous creation and fusion of vortical defects,not present in the planar versions of the model. Some of our key formulas in this biological model are the dynamicalequations Eq. 18, 19 and the rotation rates Eq. 25 and Eq. 36, which motivates experiments along the lines of Refs. [6],[7] and [9].With the advent of many advanced imaging techniques and fluorescent rotor probes [3, 5], it is now possible toperform accurate local transport measurements in curved biological membranes. Moreover, accurate particle trackingand velocity field imaging are also being conducted in many artificial setups [54, 55, 64]. Such investigations help usunderstand and characterize the biophysical properties of membranes and their impact on various cellular processes.Although such simple models overlook many complex details of a living cell, they are essential to get a first estimateof biophysical transport processes that routinely take place in living cells [65, 66]. A good understanding of flowshappening in a membrane also has great potential to aid in efficient drug delivery [4]. We hope our results in thiswork will motivate more studies and experiments in these exciting directions. VII. ACKNOWLEDGEMENTS We would like to thank Haim Diamant for many insightful discussions. This research was supported by the ISRAELSCIENCE FOUNDATION (grant No. 1752/20).0 VIII. APPENDICESA. More details on the Point Force. In this appendix section, we briefly review the computation of in-plane velocity-force response function, closelyfollowing [52, 53]. We approximate the membrane as a two-dimensional viscous incompressible fluid surroundedabove and below by three-dimensional viscous fluids, with same notations for parameters as in the main text. Greekindices are used to denote in-plane 2D objects while latin indices will be used for objects living in 3D ambient space.The incompressibility and Stokes equations for the curved membrane thus have the general form : D α v α = 0 , D β Π αβ = 0 (46)where Π αβ = pg αβ − η αβµγ D µ v γ (47)and η αβµγ = η D (cid:0) g αµ g βγ + g αγ g µβ (cid:1) + ( ξ − η d ) g αβ g µγ (48)where v denotes the in-plane 2D fluid velocity, D is the 2D covariant derivative which generalizes the partial derivativeof flat space, p is the local pressure and g µν is the 2D metric. η D is the shear viscosity of the 2D membrane fluidwhile ξ is the bulk viscosity.We now simplify the momentum equation of Eq. 46 using the incompressibility property D α v α = 0 and themetrinilic properties of the surface covariant derivative ie. D α g µν = 0, so that the metric g µν can freely pass in andout of the covariant derivative. The surface metric g µν is also used to raise/ lower appropriate indices. D β Π αβ = 0 ⇒ D β (cid:16) pg αβ − (cid:0) η D ( g αµ g βγ + g αγ g µβ ) + ( ξ − η D ) g αβ g µγ (cid:1) D µ v γ (cid:17) = 0 ⇒ D α p − η D (cid:0) D γ D α v γ + D µ D µ v α (cid:1) + ( ξ − η D ) D α (cid:24)(cid:24)(cid:24)(cid:24)(cid:58) ( D γ v γ ) = 0 ⇒ D α p − η D (cid:32) D γ D α v γ (cid:124) (cid:123)(cid:122) (cid:125) + D µ D µ v α (cid:33) = 0 (49)Let us consider the term shown in braces in the above equation. In flat space , the derivatives D γ and D α arejust partial derivatives and they may be interchanged (flat space derivatives commute), then using incompressibility ∂ γ v γ = 0 the term vanishes. However they no longer commute in curved surfaces and the commutator is proportionalto the local Gaussian curvature K ( x ). Thus, D γ D α v γ = D γ D α v γ = [ D γ ,D α ] v γ + D α (cid:8)(cid:8)(cid:8)(cid:42) D γ v γ = K ( x ) v α (50)Plugging this back to Eq. 49 we get D α p − η D (cid:0) K ( x ) v α + D µ D µ v α (cid:1) = 0 (51)Abbreviating the curved laplacian D µ D µ by ∆ , we have the final form of the Stokes equation written in terms of theGaussian curvature. η D (cid:0) K ( x ) + ∆ (cid:1) v α = D α p (52)In the limit of zero curvature, one thus recovers the usual 2D Stokes equations.We now turn to the analysis of the spectrum of the operator η D (cid:0) K ( x ) + ∆ (cid:1) . η D (cid:0) K ( x ) + ∆ (cid:1) v α ( (cid:126)x, s ) = λ s v α ( (cid:126)x, s ) (53)1where we use s to label the eigenvalues λ s and eigenfunctions v α ( (cid:126)x, s ). Any arbitrary velocity field on the curvedsurface can be decomposed in terms of the eigenfunctions v α ( (cid:126)x, s ) as follows: v α ( x ) = (cid:88) s A s v α ( (cid:126)x, s ) (54)It helps to write the velocity eigenfunctions v α ( (cid:126)x, s ) in terms of a stream function eigenmodes φ ( (cid:126)x, s ), satisfying theincompressibility requirement D α v α = 0, as follows: v α ( (cid:126)x, s ) = (cid:15) αγ D γ φ ( (cid:126)x, s ) (55)where (cid:15) is proportional to the totally antisymmetric permutation symbol e . (cid:15) αγ = √ g e αγ (56)and √ g denotes the determinant of the surface metric g µν .Plugging Eq. 55 in Eq. 54 we get the total stream function φ ( (cid:126)x ) ie. v α ( x ) = (cid:88) s A s (cid:15) αγ D γ φ ( (cid:126)x, s ) := (cid:15) αγ D γ φ ( (cid:126)x ) (57)where the total stream function is given by a decomposition in the eigenmodes φ ( (cid:126)x ) = (cid:88) s A s φ ( (cid:126)x, s ) (58)We now rewrite the eigenvalue equation Eq. 53 in terms of the stream function φ . η D (cid:0) K ( x ) D γ φ ( (cid:126)x, s ) + D γ ∆ φ ( (cid:126)x, s ) (cid:1) = λ s D γ φ ( (cid:126)x, s ) (59)For a general local Gaussian curvature, one can solve this equation numerically , however for surfaces of constantcurvature, there is an additional simplification :∆ φ ( (cid:126)x, s ) = λ s − Kη D η D φ ( (cid:126)x, s ) (60)We now turn to a discussion of the external solvents and boundary conditions. External fluids : We have ignored the role of external fluids in the fairly generic treatment above. We nowincorporate the effects of the external solvent. The external fluids are described by the usual Stokes equations in theflat ambient space. Let us denote the velocity, pressure, viscosity for r > R by + and same quantities for r < R by − . (cid:126) ∇ · (cid:126)v ± = 0 , η ± ∇ v ± = (cid:126) ∇ p ± (61)We have two boundary conditions :a) Stick boundary condition : velocities must coincide on the membrane surface. v ± | r = R = v (62)b) Stress Balance condition on the membrane σ extα = D α p − η D (cid:0) K ( x ) + ∆ (cid:1) v α + T α (63)where σ extα is the external point force or point torque applied to the membrane and T α is the traction due to couplingof the membrane fluid with the external solvents. T α = σ − αr − σ + αr | r = R , σ ± ij = η ± (cid:16) D i v ± j + D j v ± i (cid:17) − g ij p ± (64)2Rewriting the stress balance equation Eq. 63 in terms of stream function using Eq. 55 and eliminating the membranepressure by taking the antisymmetric derivative one arrives at, (cid:15) αβ D β σ extα = − (cid:88) s A s λ s ∆ φ s + (cid:15) αβ D β T α (65)where the traction T is given by Eq. 64. So far our discussions are applicable to any curved static membrane geometry.We now specialize to spherical membrane. Spherical membrane The 2D sphere metric and Levi civita are listed below along with their inverses : g αβ = R R sin θ (66) g αβ = R R sin θ (67) (cid:15) αβ = R sin θ − R sin θ (68) (cid:15) αβ = csc θR − csc θR (69)The spectrum of the Laplace Beltrami operator on a sphere is known, Eq. 60 immediately yields the eigenvalues λ s in terms of the constant curvature K and η D . For the sphere of radius R , the eigenvalue label s is given by thespherical harmonic mode index ( l, m ) . Further, K = 1 /R and the known spectrum for spherical laplacian is givenby ∆ φ lm = − l ( l + 1) R φ lm . (70)Comparing Eq. 60 and Eq. 70 gives λ l = 2 − l ( l + 1) R η D (71)and the eigenfunctions are given by φ lm = Y lm ( θ, φ ) (72)The mode decomposition for the velocity field on the sphere is given by Eq. 54 with the role of s played by ( l, m ). v α = (cid:88) lm A lm (cid:15) αγ D γ Y lm (73)where we used Eq. 72. The stress balance condition Eq. 63 will help us determine the unknown coefficients A lm inEq. 73 for the membrane velocity on the sphere in terms of the applied force. For this, we need to compute thetraction T α appearing in the stress balance condition using the knowledge of the known Lamb’s solution [67] for theexternal solvent Eq. 74. We carry out the steps below :3 External solvent in spherical co-ordinate- Lamb’s solution : Let us denote the velocities for r > R by v + and for r < R by v − . The solution of Stokes equations in 3D flat space [] is given by (cid:126)v − = ∞ (cid:88) l =1 v − l , v − l = (cid:126) ∇ × ( (cid:126)rq − l ) + (cid:126) ∇ w − l + 1 η − ( l + 1)(2 l + 3) (cid:18) 12 ( l + 3) r (cid:126) ∇ p − l − l(cid:126)rp − l (cid:19) (74)where q − l , w − l , p − l are harmonic functions of ( r, θ, φ ) ie. ∇ q − l = 0 , ∇ w − l = 0 , ∇ p − l = 0. q − l = m = l (cid:88) m = − l q − l,m r l Y lm ( θ, φ ) w − l = m = l (cid:88) m = − l w − l,m r l Y lm ( θ, φ ) p − l = m = l (cid:88) m = − l p − l,m r l Y lm ( θ, φ ) (75)Similarly for r > R the solution is obtained by the replacement l → − l − 1. However, the stick boundary conditionsdemand w − l = 0 , p − l = 0 , w + l = 0 , p + l = 0 (76)and q − lm = A lm R l +1 , q + lm = R l A lm (77)Meanwhile, the traction can be computes using the definition Eq. 64. T α = (cid:88) lm (cid:18) η − R ( l − 1) + η + R ( l + 2) (cid:19) A lm (cid:15) αβ D β Y lm ( θ, φ ) (78)Finally using the stress balance condition Eq. 65 and decomposing the point force localized at ( θ , φ ) via σ extα = F α R ∞ (cid:88) l =0 l (cid:88) m = − l Y lm ( θ, φ ) Y ∗ lm ( θ , φ ) (cid:124) (cid:123)(cid:122) (cid:125) θ δ ( θ − θ ) δ ( φ − φ ) , (79)we can solve for the unknown membrane velocity coefficients A lm in terms of force components F α A lm = csc θ η D s l l ( l + 1) (cid:0) F θ ∂ φ Y ∗ lm ( θ , φ ) − F φ ∂ θ Y ∗ lm ( θ , φ ) (cid:1) (80)where s l = l ( l + 1) − Rλ − ( l − 1) + Rλ + ( l + 2) and λ ± = η D η ± . Let us note that the traction contribution kills thezero mode l = 1. Plugging the A lm into the expression for the velocity field , we finally arrive at the stream functioncorresponding to the velocity field on the membrane surface for the point force: ψ = (cid:88) lm A lm Y lm ( θ, φ ) = (cid:88) lm csc θ η D s l l ( l + 1) (cid:0) F θ ∂ φ Y ∗ lm [ θ , φ ] − F φ ∂ θ Y ∗ lm [ θ , φ ] (cid:1) Y lm ( θ, φ ) (81)where s l = l ( l + 1) − Rλ − ( l − 1) + Rλ + ( l + 2) and λ ± = η D η ± .Performing the sum over m, this yields ψ = csc θ πη D (cid:0) F θ ∂ φ S − F φ ∂ θ S (cid:1) (82)4where S = (cid:88) l l + 1 s l l ( l + 1) P l [cos γ ] (83)where cos γ = cos ∠ (cid:0) ( θ , φ ) , ( θ, φ ) (cid:1) = sin θ sin θ cos( φ − φ ) + cos θ cos θ .The physical velocity field on the membrane surface is summarized by an Oseen tensor on S given by v θ = G θθ F θ + G θφ F φ v φ = G φθ F θ + G φφ F φ These are physical components of velocity and force as opposed to components in a covariantbasis. G θθ = csc θ csc θ πη D ∂ φ ∂ φ SG θφ = − csc θ πη D ∂ φ ∂ θ SG φθ = − csc θ πη D ∂ θ ∂ φ SG φφ = πη D ∂ θ ∂ θ S (84)We now perform the sum over Legendre Polynomials to obtain a closed form expression for S, Eq. 83. We proceedfirst by noting that the roots of s l = 0 are given by l p = − ( η d + Rη − + Rη + ) + (cid:112) η d + 6 Rη d ( η − − η + ) + R ( η − + η + ) η d l m = − ( η d + Rη − + Rη + ) − (cid:112) η d + 6 Rη d ( η − − η + ) + R ( η − + η + ) η d (85)We analyze the structure of the roots in detail in Appendix VIII E. Here we just import those results. Depending onthe parameters of the model, l p lies in the range −∞ < l p ≤ l m is always negative. BreakingS (Eq. 83) into partial fractions and summing the individual parts , one gets two different real space representationsof S depending on the sign of l p , as we show in Eq. 86 and Eq. 89 below.Further, in the situation when the internal and external Saffman lengths are same ie. λ − = λ + := λ , the analysisin Appendix VIII E shows that for large radius R (cid:29) λ , the root l p is negative. In the opposite situation of highcurvature (small radius) ie. R (cid:28) λ , the root l p is positive. Thus one has different representations of S depending onthe sign of the root l p . We use the appropriate one for our simulations. We list them below : Case 1 : −∞ < l p < S l p < = 1 l m l p (cid:16) log[2] − log( − cos γ + (cid:112) − γ + 1) (cid:17) + 1(1 + l m )(1 + l p ) log[ cos γ − √ − γ − γ − l m l m (1 + l m )( l m − l p ) A [ l m ] + 1 + 2 l p l p (1 + l p )( l p − l m ) A [ l p ] (86)where the function A [ l m ] is defined by a combination of Appell Hypergeometric function. A [ l m ] = ( − l m ) l m A [2 − l m , , , − l m ,e iγ ,e − iγ ] − ( − l m ) ( ( − l m ) A [ − l m , − , − , − l m ,e iγ ,e − iγ ]+2 l m A [1 − l m , , , − l m ,e iγ ,e − iγ ] Cos γ ) ( − l m )( − l m ) l m (87)and a similar relation for A [ l p ].The function A becomes simpler for integer values of the negative root. We list some of them below. A | l m =0 = log 2 − log (cid:104) − x + √ − x + 1 (cid:105) | x =cos γ A | l m = − = log (cid:34) x − √ − x − x − (cid:35) | x =cos γ A | l m = − = √ − x + 2 x coth − ( √ − x + 1) | x =cos γ (88)5 Case 2 : 0 < l p < S l p > = 1 l m l p (cid:16) log[2] − log( − cos γ + (cid:112) − γ + 1) (cid:17) + 1(1 + l m )(1 + l p ) log[ cos γ − √ − γ − γ − l m l m (1 + l m )( l m − l p ) A [ l m ] + 1 + 2 l p l p (1 + l p )( l p − l m ) B [ l p ] (89)where B [ l p ] = − l p + 1 − A [ − l p , , , − l p , e iγ , e − iγ ] l p (90)There are some special points in the parameter space where the above representations need to be supplemented bythe following: Case 3 : one of the roots is zeroThis enforces the other root , let us call it ˜ l = η − +4 η + η − − η + and η d = R (2 η + − η − ). Note that for η d > η + > η − , which implies ˜ l has to be negative. S = (cid:88) l l + 1( l − ˜ l ) l ( l + 1) P l [cos γ ] = − − ˜ l ˜ l (cid:16) log[2] − log( − cos γ + (cid:112) − γ + 1) (cid:17) − l S + 11 + ˜ l A [ − 1] + 1 + 2˜ l ˜ l (1 + ˜ l ) A [˜ l ] (91)where S = ∞ (cid:88) l =1 P l [cos γ ] l (92)is convergent and can be evaluated numerically. Case 4 : one of the roots is -1This enforces the other root , let us call it ˜ l = η − + η + )2 η − − η + and η d = R ( η + − η − ). Note that for η d > η + > η − , which implies ˜ l has to be negative. S = (cid:88) l l + 1( l − l p ) l ( l + 1) P l [cos γ ] = − l (cid:16) log[2] − log( − cos γ + (cid:112) − γ + 1) (cid:17) + 1 − ˜ l − S + ˜ l (˜ l + 1) A [ − 1] + 1 + 2˜ l ˜ l (1 + ˜ l ) A [˜ l ] (93)where ˜ S = ∞ (cid:88) l =1 P l [cos γ ]( l + 1) (94)is convergent and can be evaluated numerically. B. More details on the Point Torque Keeping the same notations as the force-velocity response calculation, the equation for stress balance Eq. 65 in thesituation of a rotor embedded in the spherical membrane takes the following form in the basis of spherical harmonics: − (cid:15) αγ D γ [ τ (cid:15) αβ D β ] | θ ,φ R (cid:88) l,m Y lm ( θ, φ ) Y ∗ lm ( θ , φ ) (cid:124) (cid:123)(cid:122) (cid:125) θ δ ( θ − θ ) δ ( φ − φ ) = (cid:88) lm η D l ( l + 1) R s l (cid:124)(cid:123)(cid:122)(cid:125) membrane stress +traction A lm Y lm ( θ, φ ) ,s l = l ( l + 1) − Rl − ( l − 1) + Rl + ( l + 2) . (95)6where the in-plane membrane velocity field is decomposed as Eq. 73, ie. v α = (cid:80) lm A lm (cid:15) αγ D γ Y lm and tau denotesthe rotor circulation. Solving for the unknown coefficients A lm from the above equation yields A lm = τ Y ∗ lm ( θ , φ ) η D s l (96)Plugging this into the mode expansion for the velocity field Eq. 73, we get v α = (cid:88) lm τ Y ∗ lm ( θ , φ ) η D s l (cid:15) αγ D γ | θ,φ Y lm ( θ, φ ) (97)Performing the sum over m, (cid:80) m = lm = − l Y lm ( θ , φ ) Y ∗ lm ( θ , φ ) = l +14 π P l (cos γ ) we finally have v α = (cid:88) l τ (2 l + 1)4 πη D s l (cid:15) αγ D γ | θ,φ P l (cos γ ) (98)where γ is the geodesic angle between ( θ, φ ) and ( θ , φ ).Introducing the operator [ ∇ S ⊥ ] = − (cid:16) ˆ θ R sin θ ∂ φ − ˆ φ R ∂ θ (cid:17) we find that the physical velocity field can be expressed as v = τη D [ ∇ S ⊥ ] ψ (99)where the dimensionless stream function ψ is given by ψ [ θ, φ, θ , φ ] = (cid:88) l (2 l + 1)4 π s l P l (cos γ ) (100)where s l = l ( l + 1) − Rl − ( l − 1) + Rl + ( l + 2) and λ ± = η D η ± and cos γ = sin θ sin θ cos( φ − φ ) + cos θ cos θ . Vanishing of membrane pressure in rotors : Taking the symmetric derivative D α of the stress balancecondition σ extα = D α p − (cid:80) lm A lm λ l (cid:15) αγ D γ φ lm + T α we get∆ p = 0 (101)for the case of a rotor ie. σ extα = τ (cid:15) αβ D β δ ( θ − θ , φ − φ ) and traction T α given by Eq. 78. Since ∆ is a LaplaceOperator on a compact manifold S (see [68] for details), this means p can only be a harmonic function with eigenvaluezero which forces it to be a constant and drops out of the hydrodynamic equation because it appears as a gradi-ent. Let us note that this does not happen for the external point force where we will get a non zero membrane pressure. Flat space limit : The flat space limit of Eq. 100 can be understood by introducing a momentum variable q = lR and converting the sum into an integral in the limit of large radius. ψ [ θ, φ, θ , φ ] = (cid:80) l max l =1 R (cid:124)(cid:123)(cid:122)(cid:125) dq = d q πq (2 l + 1) R π s l (cid:124) (cid:123)(cid:122) (cid:125) π ( q + λ − P l (cos γ ) (cid:124) (cid:123)(cid:122) (cid:125) e iq.r → (cid:82) d qq ( q + λ − ) e iq.r (102)where we used (2 l +1) R π s l = (2 qR +1) R π (cid:18) qR ( qR +1) − Rλ − ( qR − Rλ + ( qR +2) (cid:19) ∼ qR π (cid:18) q R + Rλ − ( qR )+ Rλ + ( qR ) (cid:19) ∼ (cid:18) q + λ − + λ + (cid:19) (103)7 FIG. 13: Comparison of the velocity component v φ vs θ (using our analytic expression Eq. 105) with the contribution to velocitypurely from global rotation term (Eq. 15 of main text) shown in gray, for the point torque positioned at the north pole. Topleft to bottom right, the curvature is progressively increased. The non-monotonicity (which appears first on the second plotfrom the left) is due to the dominance of the global rotation term as curvature is increased. Performing the sum over Legendre Polynomials In order to find closed form expression for ψ given by Eq. 100,we proceed exactly as the situation of the point force. Depending on the nature of the roots of the equation s l = 0,we again have different representations of ψ . We list them below : Case 1 : −∞ < l p < ψ l p < = 14 π (cid:32) l m + 1 l m − l p A [ l m ] − l p + 1 l m − l p A [ l p ] (cid:33) (104)with A by Eq. 87. Case 2 : 0 < l p < ψ l p > = 14 π (cid:32) l m + 1 l m − l p A [ l m ] − l p + 1 l m − l p B [ l p ] (cid:33) (105)with A and B given by Eq. 87 and 90. Case 3 : l p = 0In this situation, ψ l p =0 = (cid:80) l l +1 l ( l − l m ) P l [cos γ ]= − l m log[ − cos γ + √ − γ +1 ] + l m l m ( A [ l m ] + l m ) (106)where A is defined in Eq. 87. C. More details on the Counter Rotating Torque We now consider the situation of counter rotating torque dipole, one situated at r = R with torque τ and the otherpositioned at r = R + d with torque − τ . This is a useful first approximation of rotor proteins in a membrane [51].Following Lamb’s solution [67] in spherical co-ordinates, we first consider the solutions of the Stokes equations in 3D, as given in Eq. 74 with Eq. 76 . We list the profiles in the appropriate domains ( v − , v int , v + denote the innermost, In our notation the physical components of the 3D velocity field are given by (cid:126)v = v r ˆ r + v θ ˆ θ + v φ ˆ φ and similarly for the membrane (cid:126)v = v θ ˆ θ + v φ ˆ φ r < R we have [ v ] − r = 0 , [ v ] − θ = (cid:80) lm θ q − lm r l ∂ φ Y lm ( θ, φ )[ v ] − φ = − (cid:80) lm q − lm r l ∂ θ Y lm ( θ, φ ) (107)and for r > R + d are given by [ v ] + r = 0 , [ v ] + θ = (cid:80) lm θ q + lm r − l − ∂ φ Y lm ( θ, φ )[ v ] + φ = − (cid:80) lm q + lm r − l − ∂ θ Y lm ( θ, φ ) (108)For R < r < R + d we have both the rising and falling solutions in the intermediate velocity field v int .[ v int ] − r = 0 , [ v int ] − θ = (cid:80) lm (cid:16) θ ˜ q − lm r l ∂ φ Y lm ( θ, φ ) + θ ˜ q + lm r − l − ∂ φ Y lm ( θ, φ ) (cid:17) [ v int ] − φ = − (cid:80) lm (cid:16) ˜ q − lm r l ∂ θ Y lm ( θ, φ ) + (cid:80) lm ˜ q + lm r − l − ∂ θ Y lm ( θ, φ ) (cid:17) (109)The physical components of the membrane velocity field are as follows : v θ == R sin θ (cid:80) lm A lm ∂ φ Y lm ( θ, φ ) v φ = − R (cid:80) lm A lm ∂ θ Y lm ( θ, φ ) . (110)Thus altogether we have 5 undetermined coefficients A lm , q − lm , q + lm , ˜ q + lm , ˜ q − lm which are to be determined via 5 equations: 3 equations from velocity matching and 2 stress balance equations at r = R and r = R + d . Velocity matching at r = R gives 2 equations since the membrane velocity field has to match with v int and v − simultaneously at r = R .This gives A lm R = q − lm R l = ˜ q − lm R l + ˜ q + lm R − l − (111)Further v int has to match with v + at r = R + d . This gives us˜ q − lm ( R + d ) l + ˜ q + lm R + d ) l +1 = q + lm ( R + d ) l +1 (112) Traction at r = R + d : ˜ T α = η + (cid:88) lm ˜ q − lm ( R + d ) l ( l − 1) + (˜ q + lm − q + lm )( R + d ) − l − ( − l − (cid:15) αβ D β Y lm (113)Stress balance at r = R + d : The equation for stress balance in this surface(plugged ˜ T from Eq. 113) gives τ ( R + d ) Y ∗ lm ( θ , φ ) = η + (cid:16) ˜ q − lm ( R + d ) l ( l − 1) + (˜ q + lm − q + lm )( R + d ) − l − ( − l − (cid:17) (114) Traction at r = R : T α = (cid:104) η − q − lm ( l − R l − η + ˜ q − lm ( l − R l − η + ˜ q + lm ( − l − R − l − (cid:105) (cid:15) αβ D β Y lm (115)9The stress balance at r = R now becomes − τR Y ∗ lm ( θ , φ ) = − A lm (2 − l ( l +1)) R η D + (cid:16) η − q − lm ( l − R l − η + ˜ q − lm ( l − R l − η + ˜ q + lm ( − l − R − l − (cid:17) (116)Thus overall we have the following set of 5 equations : A lm R = q − lm R l = ˜ q − lm R l + ˜ q + lm R − l − ˜ q − lm ( R + d ) l + ˜ q + lm R + d ) l +1 = q + lm ( R + d ) l +1 τ ( R + d ) Y ∗ lm ( θ , φ ) = η + (cid:16) ˜ q − lm ( R + d ) l ( l − 1) + (˜ q + lm − q + lm )( R + d ) − l − ( − l − (cid:17) − τR Y ∗ lm ( θ , φ ) = − A lm (2 − l ( l +1)) R η D + (cid:16) η − q − lm ( l − R l − η + ˜ q − lm ( l − R l − η + ˜ q + lm ( − l − R − l − (cid:17) (117)The above system of five equations can be solved for the five unknowns A lm , q − lm , q + lm , ˜ q + lm , ˜ q − lm .The solution is given by A lm = τ ( R + d ) − − l ( R l − d ( d + R ) l − dR ( d + R ) l − R ( d + R ) l ) ( − l + l ) η D + R [ η − ( l − η + ( l +2) ] Y ∗ lm ( θ , φ )= τη D ( R + d ) − − l ( R l − d ( d + R ) l − dR ( d + R ) l − R ( d + R ) l ) s l Y ∗ lm ( θ , φ ) (118) q − lm = τη D R − l − ( d + R ) − − l ( R l − d ( d + R ) l − dR ( d + R ) l − R ( d + R ) l ) Y ∗ lm ( θ ,φ ) s l q + lm = τη D ( d + R ) − − l ( (1 − l ) R l ((2+ l ) η d + R ( η − − η + )) − (1+2 l ) R l ( d + R ) l η p +( d + R ) l (( − l )((2+ l ) η d + Rη − )+(2+ l ) Rη + ) ) (1+2 l ) η + s l Y ∗ lm ( θ , φ )˜ q + lm = τη d R l ((1 − l ) R l ( d + R ) − − l ((2+ l ) η d + R ( η − − η + )) − (1+2 l ) η + (1+2 l ) η + s l Y ∗ lm ( θ , φ )˜ q − lm = ( d + R ) − − l τY ∗ lm ( θ ,φ )(2 l +1) η + (119) v α = (cid:80) lm A lm (cid:15) αγ D γ Y lm = (cid:80) lm τη D ( R + d ) − − l ( R l − d ( d + R ) l − dR ( d + R ) l − R ( d + R ) l ) s l Y ∗ lm ( θ , φ ) (cid:15) αγ D γ Y lm ( θ, φ ) (120)Performing the sum over m, we get v α = (cid:88) l τη D C l s l l + 14 π (cid:15) αγ D γ P l (cos γ ) (121)where C l := ( R + d ) − − l (cid:16) R l − d ( d + R ) l − dR ( d + R ) l − R ( d + R ) l (cid:17) = d ( − − l ) R + O ( d ) (122)Introducing the operator [ ∇ S ⊥ ] = − (cid:16) ˆ θ R sin θ ∂ φ − ˆ φ R ∂ θ (cid:17) we find that the physical velocity field can be expressed as v = τη D [ ∇ S ⊥ ] ψ (123)where the dimensionless stream function ψ is given by ψ [ θ, φ, θ , φ ] = (cid:88) l (2 l + 1) C l π s l P l (cos γ ) (124)0where s l = l ( l + 1) − Rl − ( l − 1) + Rl + ( l + 2) and l ± = η D η ± and cos γ = sin θ sin θ cos( φ − φ ) + cos θ cos θ and C l isdefined in Eq. 122.Let us also note from Eq. 122 we note that the velocity field vanishes when the distance d betweenthe counter rotating torques go to zero. Flat space limit : The flat space limit of Eq. 124 can be understood by again introducing a momentumvariable q = lR and converting the sum into an integral in the limit of large radius. ψ [ θ, φ, θ , φ ] = (cid:80) l max l =1 R (cid:124)(cid:123)(cid:122)(cid:125) dq = d q πq (2 l + 1) C l R π s l (cid:124) (cid:123)(cid:122) (cid:125) qd π ( q + λ − + O ( d ) P l (cos γ ) (cid:124) (cid:123)(cid:122) (cid:125) e iq.r → d × (cid:82) d q ( q + λ − ) e iq.r (125)where we used (2 l +1) C l R π s l = (2 qR +1) dR ( − − qR ) R π (cid:18) qR ( qR +1) − Rλ − ( qR − Rλ + ( qR +2) (cid:19) ∼ dq R π (cid:18) q R + Rλ − ( qR )+ Rλ + ( qR ) (cid:19) ∼ q (cid:18) q + λ − + λ + (cid:19) (126)Thus the velocity has dimensions τdη D × λ .We need to perform the sum ( to lowest order in d ). Once again, the roots of the equation s l = 0 are given by Eq. 85. Case 1 : −∞ < l p < ψ l p < = − d πR (cid:20) √ − γ + l m +2 l m l m − l p A [ l m ] + − − l p − l p l m − l p A [ l p ] (cid:21) (127)where A is defined by Eq. 87. Case 2 : If 0 < l p < ψ l p > = − d πR (cid:20) √ − γ + l m +2 l m l m − l p A [ l m ] + − − l p − l p l m − l p B [ l p ] (cid:21) (128)where A and B are defined by Eq. 87 and Eq. 90. Case 3 : l p = 0In this situation, ψ l p =0 = (cid:80) l (2 l +1)( − l − l ( l − l m ) P l [cos γ ]= − √ − γ + l m log[ − cos γ + √ − γ +1 ] + − − l m − l m l m ( A [ l m ] + l m ) (129)where A is defined in Eq. 87. D. Analytical Investigations of Streamline Topologies In this section, we provide details of the derivation of Eq. 32 of main text. We initially consider N rotors and laterspecialize to two rotors. One can project the dynamical equations Eq. 18 of main text via stereographic projectionon the plane. If we denote the planar polar co-ordinates by ( r, ˜ θ ), then the stereographic map relates ( r, ˜ θ ) tocoordinates ( θ, φ ) on the sphere via the relations ˜ θ = φr = tan θ (130)1Using this mapping, the hydrodynamic evolution equations take the form ddt r i = η D R (cid:80) Nj (cid:54) = i τ j (1+ r i ) ∂ ˜ θ i ψ [ γ ij ] ddt ˜ θ i = η D R (cid:80) Nj (cid:54) = i − τ j (1+ r i ) ∂ r i ψ [ γ ij ] (131)where γ ij = arccos (cid:32) (1 − r i )(1 − r j ) + 4 r i r j cos(˜ θ i − ˜ θ j )(1 + r i )(1 + r j ) (cid:33) (132)Similarly, the equation of a tracer particle (denoted by suffix p) moving in the presence of N rotors can be written interms of H p as follows : ddt r p = η D R (1+ r p ) ∂ ˜ θ p H pddt ˜ θ p = η D R − (1+ r p ) ∂ r p H p (133) H p = N (cid:88) j τ j ψ [ γ pj ] (134)and γ pj = arccos (cid:32) (1 − r p )(1 − r j ) + 4 r p r j cos(˜ θ p − ˜ θ j )(1 + r p )(1 + r j ) (cid:33) (135)Introducing complex coordinates on the plane z = re i ˜ θ we can write Eq. 133 in complex notation ddt z ∗ p = iη D R (1 + | z p | ) ∂ z p H p (136)where H p is the same as defined in Eq. 20 with the geodesic distance in complex notation given by γ pj = arccos (cid:32) (1 − | z p | )(1 − | z j | ) + 4 Re [ z p ¯ z j ](1 + | z p | )(1 + | z j | ) (cid:33) (138)In this section, we perform an analytical treatment of the location of stagnation points on the spherical membrane inthe regimes of low and high curvature separately. 1. Regime of low curvature In general, it follows from Eq. 136 that solving for stagnation points amounts to solving for solutions to ddt ¯ z p = iη D R (1 + | z p | ) ∂ z p N (cid:88) j τ j ψ [ γ pj ] = 0 (139) Let us note that upon substituting the stream function for ideal vortices given by ψ [ γ pj ] = log(1 − cos γ pj ) into Eq. 20, one gets thestandard Hamiltonian for ideal vortices on the sphere given by H idealp = N (cid:88) j τ j log (cid:32) | z p − z j | (1 + | z p | )(1 + | z j | ) (cid:33) . (137) γ pj = arccos (cid:32) (1 − | z p | )(1 − | z j | ) + 4 Re [ z p ¯ z j ](1 + | z p | )(1 + | z j | ) (cid:33) (140)and the stream function ψ given by Eq. 100. Because the stream function is complicated in structure after performingthe Legendre sum, we propose here to choose a set of parameters that enables us to simplify the stream function andsubsequent analysis of stagnation points.To be concrete, let us choose η D = 3 / , η − = 1 , η + = 2 , R = 1 for which λ/R = 1 / 2. This yields the tworoots of s l = 0 to be l m = − , l p = 0Using Case 3 of summed up versions of Eq. 100, we get ψ [ γ ] = 112 π (cid:20) (cid:2) (6 cos γ − 2) arcCoth( (cid:112) − γ + 1) + 3 cos γ ( (cid:112) − γ − 1) + (cid:112) − γ (cid:19) − log( − cos γ + (cid:112) − γ + 1) (cid:3) (141)Plugging Eq. 141 into Eq. 139 iη D R (1 + | z p | ) N (cid:88) j τ j F [ z p , z j ] G [ z p , z j ] = 0 (142)where the factors F and G arise from the derivative of the stream function ie. ∂ z p ψ = ∂ ψ ∂ cos γ ∂ z p cos γ := F × G . Letus note that although F is dependent on the choice of parameters, the factor G is essentially purely geometric. Forour choice of parameters, F [ z p , z j ] = − √ − γ pj +cos γ pj (cid:16) − √ − γ pj +15( − √ − γ pj ) cos γ pj − 30 cos γ pj sin γ pj arcCoth(1+ √ − γ pj ) (cid:17) − π sin γ pj G [ z p , z j ] = (1 −| z j | )( − z p )+4 (cid:16) (1+ | z p | ) ¯ zj − Re [ z p ¯ z j ]¯ z p (cid:17) (1+ | z j | )(1+ | z p | ) (143)where in the expression of F we havecos γ pj = (cid:32) (1 − | z p | )(1 − | z j | ) + 4 Re [ z p ¯ z j ](1 + | z p | )(1 + | z j | ) (cid:33) (144)We now specialize to the case of two rotors on the spherical membrane. Let their positions in the complex planebe denoted by z and z . Since the stagnation points are always constrained to lie on the great circle joining thetwo locations, we can essentially map the dynamics to the unit circle on the complex plane. We choose coordinatessuch that the location of first rotor is at ˜ θ = 0 and the second rotor at ˜ θ = φ . We further choose, without loss ofgenerality, the strength of the first rotor to be 1 and relative strength between the rotors be denoted by τ . Thus z = 1 , z = e iφ , τ = 1 , τ = τ (145)Plugging in Eq. 145 into Eq. 142 and Eq. 143,we convert it to a effective two parameter problem where the stagnationpoint z p has to be solved as a function of the relative vortex strength τ and the location of the second vortexparametrized by φ , from the equation F [ z p , G [ z p , 1] + τ F [ z p , e iφ ] G [ z p , e iφ ] = 0 ⇒ f [ Re ( z p )] (1 − ¯ z p ) + τ f [ Re ( z p e − iφ )] (cid:0) e − iφ − Re ( z p e − iφ )¯ z p (cid:1) = 0 (146)where f = 10 − √ − x + x (cid:16) − √ − x + 15( − √ − x ) x + 30 x ( x − 1) arcCoth(1 + √ − x ) (cid:17) π ( x − 1) (147)Substituting ansatz z p = e iθ p into Eq. 148, we get f [cos θ p ] 14 (1 − e − iθ p ) + τ f [cos( θ p − φ )] 12 (cid:16) e − iφ − cos( θ p − φ ) e − iθ p (cid:17) = 0 (148)3 2. Regime of high curvature In the regime of high curvature, ignoring all local corrections, F = d ψ d (cos γ ) = − η D πRη + (149)This follows directly from the fact that the global rotation term is ddγ ψ = η D πRη + sin γ (150)Using this F, the equation for stagnation points simplifies considerably. G ( z p , 1) + τ G ( z p , e iφ ) ∼ ⇒ 14 (1 − ¯ z p ) + τ (cid:16) e − iφ − Re [ z p e − iφ ]¯ z p (cid:17) = 0 (151)Substituting z p = e iθ p in the above we get14 (1 − e − iθ ) + τ (cid:16) e − iφ − cos( θ − φ ) e − iθ (cid:17) = 0 (152)This is the equation discussed in main text Eq. 40. E. Roots of s l : Poles of the stream function in Legendre basis on the sphere. In all the examples we studied in the main text, the dimensionless stream function on the spherical membrane hasthe following generic structure in the basis of Legendre polynomials: ψ [ θ, φ, θ , φ ] = (cid:88) l f l π s l g l P l (cos γ ) (153)where f l and g l are some polynomials in Legendre modes denoted by l and s l = l ( l + 1) − Rλ − ( l − 1) + Rλ + ( l + 2).The geodesic angle between the source and response locations is denoted by γ . In order to find the real space Greensfunction, one is thus left with the task of performing the sum Eq. 153. As mentioned in Appendix VIII A, VIII B,VIII C, the real space representation of the stream function crucially depends on the root structure of the equation s l = 0.In this appendix we discuss the nature of the roots of the equation s l = 0. s l = l ( l + 1) − R η − η D ( l − 1) + Rη + η D ( l + 2) = 0 ⇒ l + l (cid:18) Rη − η d + Rη + η d (cid:19) + (cid:32) − − Rη − η d + 2 Rη + η d (cid:33) = 0 (154) l p = − ( η d + Rη − + Rη + ) + (cid:112) η d + 6 Rη d ( η − − η + ) + R ( η − + η + ) η d l m = − ( η d + Rη − + Rη + ) − (cid:112) η d + 6 Rη d ( η − − η + ) + R ( η − + η + ) η d (155)Let us now discuss the nature of the roots in the space of parameters : Nature of the root l m : always negative. Nature of the root l p : The range of this root is −∞ < l p ≤ l p .In order to understand this better, let us first consider the simpler situation η + = η − := η d . Defining the uniqueSaffman length as λ := η d η d , we first note from Eq. 154 that the product of the two roots is l m l p = R λ − . (156)Since l m is always negative, it is clear that for large radius R > λ , the root l p is negative. In the opposite situationof high curvature (small radius) ie. R < λ , the root l p is positive.To explore the more generic situation where η + (cid:54) = η − , let us consider expansions of l p in terms of radius R.For small radius, l p = 1 − Rη + η d + O ( R ) (157)Thus shows that l p < l p < 0. This demands l p < , l m l p = (cid:32) − − Rη − η d + 2 Rη + η d (cid:33) > η + − η − > , η d < R η + − η − ) (159)For large radius , l p = − η + − η − η + + η − + O (1 /R ) (160)One notes that now the sign of l p is more subtle, positive if η − > η + and negative for η − < η + . [1] H. Ueno, T. Suzuki, K. Kinosita, and M. Yoshida, ATP-driven stepwise rotation of FoF1-ATP synthase , Proc. Natl. Acad.Sci. , 1333 (2005).[2] C. Duprat, H. A. Shore. Fluid-Structure Interactions in Low-Reynolds-Number Flows Royal Society of Chemistry, 2015.[3] Kuimova et al, Molecular Rotor Measures Viscosity of Live Cells via Fluorescence Lifetime Imaging , JACS , ,6672(2008).[4] Piro, Tang and Golestanian, Optimal navigation strategies for microswimmers on curved manifolds , arXiv 2010.07580.[5] Wang et al, Probing cell membrane damage using a molecular rotor probe with membrane-to-nucleus translocation , MaterialsHorizons, The Royal Society of Chemistry, 2020.[6] Grzybowski, Stone and Whitesides, Dynamic self-assembly of magnetized, millimetre-sized objects rotating at a liquid-airinterface , Nature ,1033 (2000).[7] Terray, Oakey and Marr , Microfluidic control using colloidal devices , Science , 1841 (2002).[8] E. Lauga, W. R. DiLuzio, G. M. Whitesides, and H. A. Stone. Swimming in circles: motion of bacteria near solid boundaries Biophys. J. , (2006).[9] Drescher, Leptos, Tuval, Ishikawa, Pedley and Goldstein, Dancing Volvox : Hydrodynamic Bound States of SwimmingAlgae , Phys. Rev. Lett. , 168101 (2009).[10] N.D. Mermin, The topological theory of defects in ordered media , Reviews of Modern Physics , 591 (1979).[11] David. R. Nelson, Defects and geometry in condensed matter physics , Cambridge University Press, 2002.[12] A. A. Abrikasov, Type II superconductors and the vortex lattice , Reviews of Modern Physics, (2004).[13] M. Tinkham, Introduction to Superconductivity , Dover publications, New York.[14] V. Berezinskii, Destruction of long-range order in one-dimensional and two dimensional systems having a continuoussymmetry group: classical and quantum systems , Soviet Journal of Experimental and Theoretical Physics , 493 (1971)and , 610 (1972).[15] J.M. Kosterlitz, D.J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems , Journal of PhysicsC : Solid State Physics, , 1181, (1972). [16] Bogomolov, Dynamics of vorticity at a sphere , 1977.[17] Kimura, Vortex motion on a sphere , Journal of Physical Society of Japan, 1987.[18] H. Aref, M. Brons, On stagnation points and streamline topology in vortex flows , J. Fluid Mech. , 1-27 (1998).[19] Kidambi and Newton, Streamline Topologies for integrable vortex motion on a sphere , Physica D: Nonlinear Phenomena , 95-125 (2000)[20] Boatto and Koiller, Book chapter: Geometry, Mechanics and Dynamics , Fields Institute Communications, , 185,Springer, New York (2015).[21] Dritschel and Boatto, The motion of point vortices on closed surfaces , Proc. R. Soc. A : 20140890 (2015).[22] I. M. Held, R. T. Pierrehumbert, S. T. Garner, and K. L. Swanson. Surface quasi-geostrophic dynamics J. Fluid Mech. Anomalous coupling between topological defects and curvature , PRL, , 215301 (2004).[24] Turner, Vitelli and Nelson, Vortices on curved surfaces , Rev. Mod. Phys. , 1301 (2010).[25] Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren, Topological defects in the nematicorder of actin fibers as organization centers of Hydra morphogenesis bioRxiv (2020).[26] Shankar, Souslov, Bowick, Marchetti, Vitelli, Topological Active Matter , arXiv 2010.00364.[27] Sone and Ashida, Anomalous Topological Active Matter , PRL , 205502(2019).[28] Shankar, Ramaswamy, Marchetti and Bowick, Defect unbinding in active nematics , PRL , 108002 (2018).[29] Giomi, Bowick, Mishra, Sknepnek and Marchetti, Defect dynamics in active nematics , Philosophical Transactions of theRoyal Society A, , 20130365 (2014).[30] Giomi, Geometry and Topology of Turbulence in active nematics , PRX, , 031003 (2015).[31] Giomi, Kos, Ravnik and Sengupta, Cross-talk between topological defects in different fields revealed by nematic microfluidics ,PNAS, 2017.[32] Opathalage et al, Self-Organized dynamics and the transition to turbulance of confined active nematics , PNAS (2019).[33] Duclos et al, Topological structure and dynamics of three dimensional active nematics , Science, , 1120, 2020.[34] T. Sanchez, D. TN. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic. Spontaneous motion in hierarchically assembledactive matter Nature , (2012).[35] Tan, Liu, Miller, Tekant, Dunkel and Fakhri, Topological turbulence in the membrane of a living cell , Nature Physics, ,657 (2020)[36] Shin, Bowick and Xing, Topological Defects in Spherical Nematics , PRL 101.037802 (2008).[37] Henkes, Marchetti and Sknepnek, Dynamical patterns in active nematics on a sphere , Phys. Rev. E , 042605(2018).[38] Shankar, Bowick and Marchetti, Topological Sound and Flocking on Curved Surfaces , Phys. Rev. X , 031039(2017).[39] P. G. Saffman, Brownian motion in thin sheets of viscous fluid , J. Fluid Mech. 73:593–602, (1975).[40] P. G. Saffman and M. Delbr¨uck, Brownian motion in biological membranes , Proc. Natl. Acad. Sci. USA. 72:3111–3113(1975).[41] B.D. Hughes, B. A. Pailthorpe, and L. R. White. The translational and rotational drag on a cylinder moving in a membrane ,J. Fluid Mech. :349–372 (1981).[42] Evans E., Sackmann E. Translational and rotational drag coefficients for a disk moving in a liquid membrane associatedwith a rigid substrate , J. Fluid Mech. :553–561 (1988).[43] Levine A.J., MacKintosh F.C. Dynamics of viscoelastic membranes , Phys. Rev. E , Stat. Nonlin. Soft Matter Phys. :061606 (2002)[44] Levine A.J., Liverpool T.B., MacKintosh F.C. Dynamics of rigid and flexible extended bodies in viscous films and mem-branes , Phys. Rev. Lett. :038102 (2004).[45] N. Oppenheimer and H. Diamant, Correlated diffusion of membrane proteins and their effect on membrane viscosity ,Biophysical Journal , 3041–3049 (2009).[46] N. Oppenheimer and H. Diamant Correlated dynamics of inclusions in a supported membrane , Physical Review E 82,041912 (2010).[47] N. Oppenheimer and H. Diamant, Dynamics of membranes with immobile inclusions , Physical Review Letters 107, 258102(2011).[48] H. A. Stone, and A. Ajdari. Hydrodynamics of particles embedded in a flat surfactant layer overlying a subphase of finitedepth , Journal of Fluid Mechanics (1998).[49] K. Seki, S. Mogre, S. Komura, Diffusion coefficients in leaflets of bilayer membranes , Phys. Rev. E. , 022713 (2014).[50] Lenz, Joanny, Julicher and Prost, Membranes with rotating motors , PRL91, 108104 (2003).[51] N.Oppenheimer, D.B. Stein, M.J.Shelley, Rotating Membrane Inclusions Crystallize Through Hydrodynamic and StericInteractions , Physical Review Letters ,148101, (2019).[52] Mark L Henle, R McGorty, AB Schofield, AD Dinsmore, AJ Levine, The effect of curvature and topology on membranehydrodynamics , EPL (Europhysics Letters), , 48001 (2008 ).[53] Mark L Henle, Alex J Levine, Hydrodynamics in curved membranes: The effect of geometry on particulate mobility , PhysicalReview E, , 011905 (2010).[54] Smith et al, Membrane Viscosity Determined from Shear-Driven Flow in Giant Vesicles , PRL , 038103(2013).[55] F. Quemeneura, J.K. Sigurdsson, M. Rennerf, P.J. Atzberger, P. Bassereaua, and D. Lacoste, Shape matters in proteinmobility within membranes , Proceedings of the National Academy of Sciences (PNAS), , No. 14, pg. 5083–5087, (2014).[56] J.K. Sigurdsson and P.J. Atzberger, Hydrodynamic Coupling of Particle Inclusions Embedded in Curved Lipid BilayerMembranes , , 6685-6707, Soft Matter, The Royal Society of Chemistry, (2016).[57] Gross and Atzberger, Hydrodynamic Flows on Curved Surfaces : Spectral Numerical Methods for radial manifold shapes , Journal of Computational Physics, , 663 (2018).[58] Rower, Padidar and Atzberger, Surface Fluctuating Hydrodynamics methods for the drift-diffusion dynamics of particlesand microstructures within curved fluid interfaces , arXiv 1906.01146.[59] S. Viefers, Quantum Hall physics in rotating Bose-Einstein condensates , Journal of Physics:Condensed Matter, (2008)[60] M.C. Beeler et al, The Spin Hall effect in a quantum gas , Nature , 201(2013)[61] T. Byrnes, J. Dowling, Quantum Hall Effect with small number of vortices in Bose Einstein condensates , Phys. Rev. A , 023629 (2015).[62] Wolfram Research, Inc., Mathematica, Champaign, IL.[63] R. Samanta and N. Oppenheimer, in preparation.[64] Li, Santos-Otte, Au, Rentsch, Block and Ewers, Directed manipulation of membrane proteins by fluorescent magneticnanoparticles , Nature Communications, , 4259 (2020).[65] Enkavi, Javanainen, Kulig, Rog and Vattulinen, Multiscale Simulations of Biological Membranes: The Challenge To Un-derstand Biological Phenomena in a Living Substance , Chemical Reviews, , 5607 (2019).[66] Sarkis and Vie, Biomimetic Models to Investigate Membrane Biophysics Affecting Lipid–Protein Interaction , Frontiers inBioengineering and Biotechnology, , 270 (2020).[67] H. Lamb, Hydrodynamics (University Press, 1924)[68] Heiko Gimperlein,