Surface Fluctuating Hydrodynamics Methods for the Drift-Diffusion Dynamics of Particles and Microstructures within Curved Fluid Interfaces

TL;DR

This paper develops fluctuating hydrodynamics methods for curved fluid interfaces to study particle dynamics, revealing unique power-law behaviors and providing theoretical explanations for surface-specific dissipation phenomena.

Contribution

It introduces novel surface fluctuating hydrodynamics approaches, including immersed boundary and stochastic Eulerian-Lagrangian methods, for curved interfaces with thermal fluctuations.

Findings

Power-law velocity autocorrelation scalings depend on geometry and viscosities.

Surface hydrodynamics differ from bulk fluid behavior, showing unique dissipation time-scales.

Methods successfully model passive particles and microswimmers on curved interfaces.

Abstract

We introduce fluctuating hydrodynamics approaches on surfaces for capturing the drift-diffusion dynamics of particles and microstructures immersed within curved fluid interfaces of spherical shape. We take into account the interfacial hydrodynamic coupling, traction coupling with the surrounding bulk fluid, and thermal fluctuations. For fluid-structure interactions, we introduce Immersed Boundary Methods (IBM) and related Stochastic Eulerian-Lagrangian Methods (SELM) for curved surfaces. We use these approaches to investigate the statistics of surface fluctuating hydrodynamics and microstructures. For velocity autocorrelations, we find characteristic power-law scalings $\tau^{-1}$, $\tau^{-2}$, and plateaus can emerge. This depends on the physical regime associated with the geometry, surface viscosity, and bulk viscosity. This differs from the characteristic $\tau^{-3/2}$ scaling for bulk three dimensional fluids. We develop theory explaining these observed power-laws associated with time-scales for dissipation within the fluid interface and coupling to the surrounding fluid. We then use our introduced methods to investigate a few example systems and roles of hydrodynamic coupling and thermal fluctuations including for the kinetics of passive particles and active microswimmers in curved fluid interfaces.

Figures

Tables

Table 1. Table 1: Parameters used for the pair diffusivity study and simulations. We use SI units of atomic-mass-units (amu), nanometers (nm), and picoseconds (ps).

Parameter	Value	Parameter	Value
$R$	15.3 nm	$R_{particle}$	1.5 nm
$k_{B}T$	2.48 amu nm2 / ps	$D$	3e-7 nm2/ps
${\gamma }_{drag}^{-1}$	$D/k_{B}T$	${\gamma }_{slip}$	0.15
${\tau }_D$	$R^2/D$	$\Delta t$	5e-6${\tau }_D$
$L^{+}$	100

Table 2. Table 2: Parameter Values for the Swimmer Simulations. We use SI units of atomic-mass-units (amu), nanometers (nm), and picoseconds (ps).

Variables	Values	Variables	Values
$R$	15.3nm	$R_{particle}$	1.5nm
$k_{B}T$	2.48 amu nm2 / ps2	$D$	9.4577e-1 nm2 / ps
${\tau }_D$	$R^2/D$	$\Delta t$	5e-7${\tau }_D$
$t$	$5e5\Delta t$	${\gamma }_{slip}$	2.6222e-3 ps/amu
resI	100	$\gamma$I	800
$N$	3 particles	$k$	$5e2k_{B}T$
${ϵ}_{wca}$	5$k_{B}T$	${\sigma }_{wca}$	2$R_{particle}$
${\varphi }_0$	0	${\varphi }_2$	$\pi /2$
$A$	$R_{particle}$	$l$	2$R_{particle}$
$\omega$	5e-3

Table 3. Table 3: Estimated swimmer velocity from the simulated swimmer displacements with and without hydrodynamic coupling in Figure  10 . Velocities normalized by v 0 = 4.55 × 10 3 ​ n ​ m / p ​ s subscript 𝑣 0 4.55 superscript 10 3 𝑛 𝑚 𝑝 𝑠 v_{0}=4.55\times 10^{3}nm/ps when L / R = 6.5 𝐿 𝑅 6.5 L/R=6.5 . The non-hydrodynamic swimmer simulations did not exhibit any significant net motions v / v 0 ∼ 5.93 × 10 − 4 ​ n ​ m / p ​ s similar-to 𝑣 subscript 𝑣 0 5.93 superscript 10 4 𝑛 𝑚 𝑝 𝑠 v/v_{0}\sim 5.93\times 10^{-4}nm/ps . We remark that this parameterization in Table  3 is for illustrative purposes of the methods and to obtain more physically realizable systems may require further adjustments.

$L/R$	$v/v_0$ (hydro)	$L/R$	$v/v_0$ (hydro)
0.13	1.167330	6.5	1.000000
0.65	1.138532	52	0.881074

Equations

$$(1)$$

$$\displaystyle\left\{\begin{array}[]{llll}\rho\dot{\mathbf{v}}&=&\overline{\mbox{div}}\left(\boldsymbol{\sigma}\right)+\rho\bar{\mathbf{b}}\\ \dot{\rho}+\rho\left(\overline{\mbox{div}}\left(\mathbf{v}_{\parallel}\right)+\mathbf{v}_{n}H\right)&=&0.\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{llll}\rho\dot{\mathbf{v}}&=&\overline{\mbox{div}}\left(\boldsymbol{\sigma}\right)+\mathbf{b}\\ \overline{\mbox{div}}\left(\mathbf{v}\right)&=&0.\end{array}\right.$$

$$\mathbf{\sigma }={\mu }_m\mathbf{D}-p\mathcal{I}.$$

$$\overline{\text{mbox}div}{\left(\mathbf{\sigma }\right)}^{♭}={\mu }_m\left(-\mathbf{\delta }\mathbf{d}{\mathbf{v}}^{♭}+2K{\mathbf{v}}^{♭}\right)-\mathbf{d}p.$$

$$\rho \frac{d{\mathbf{v}}^{♭}}{dt}$$

$$-\mathbf{\delta }{\mathbf{v}}^{♭}$$

$$m\frac{d\mathbf{V}}{dt}$$

$$\frac{d\mathbf{X}}{dt}$$

$$\mathbf{r}\cdot \nabla \times \mathbf{u}$$

$${\mathbf{u}}^{+}=\sum _{\ell =0}^{\infty }{\mathbf{u}}_{\ell }^{+},{\mathbf{u}}^{-}=\sum _{\ell =1}^{\infty }{\mathbf{u}}_{\ell }^{-}.$$

$${\mathbf{u}}_{\ell }^{+}=\nabla \times \left(\mathbf{r}{\chi }_{-(\ell +1)}\right),{\mathbf{u}}_{\ell }^{-}=\nabla \times \left(\mathbf{r}{\chi }_{\ell }\right),$$

$${\chi }_{\ell }$$

$${\chi }_{-(\ell +1)}$$

$${\mathbf{t}}^{+}={\mathbf{\sigma }}^{+}\cdot {\mathbf{n}}^{+}={\mu }_{+}\frac{\partial {\mathbf{u}}^{+}}{\partial r}+{\mu }_{+}\nabla \left({\mathbf{u}}^{+}\cdot {\mathbf{n}}^{+}\right)$$

$${\mathbf{t}}^{-}={\mathbf{\sigma }}^{-}\cdot {\mathbf{n}}^{-}=-{\mu }_{-}\frac{\partial {\mathbf{u}}^{-}}{\partial r}+{\mu }_{-}\nabla \left({\mathbf{u}}^{-}\cdot {\mathbf{n}}^{-}\right).$$

$${\mathbf{t}}^{+}$$

$${\mathbf{t}}^{-}$$

$$⟨{\mathbf{f}}_{thm}(t){\mathbf{f}}_{thm}(s{)}^T⟩$$

$$⟨{\mathbf{F}}_{thm}(t){\mathbf{F}}_{thm}(s{)}^T⟩$$

$$⟨{\mathbf{F}}_{thm}(t){\mathbf{f}}_{thm}(s{)}^T⟩$$

$${\mathcal{L}}_{ff}$$

$${\mathcal{L}}_f$$

$$⟨{\mathbf{g}}_{thm}(t){\mathbf{g}}_{thm}(s{)}^T⟩$$

$$⟨\mathbf{u}(\mathbf{x})\tilde{\mathbf{u}}(\mathbf{y}{)}^T⟩=(29)$$

$$-\int \int \int \int \mathbf{\alpha }(\mathbf{z},r,\mathbf{x}{)}^{T}2k_{B}T{\mathcal{L}}_{ff}\tilde{\mathbf{\alpha }}(\mathbf{w},q,\mathbf{y})d\mathbf{w}d\mathbf{z}\cdot$$

$$\delta (r-q)drdq.$$

$$⟨\mathbf{U}{\tilde{\mathbf{U}}}^T⟩=\int \int 2k_{B}T\gamma \mathbf{A}(r{)}^T\tilde{\mathbf{A}}(q)\delta (r-q)drdq,$$

$$(31)$$

$$⟨\mathbf{U}\tilde{\mathbf{u}}(\mathbf{y}{)}^T⟩=-\int \int \int 2k_{B}T\gamma \mathbf{A}(r{)}^T\Gamma \tilde{\mathbf{\alpha }}(\mathbf{w},q,\mathbf{y})\cdot$$

$$\delta (r-q)d\mathbf{w}drdq.$$

$$⟨\Gamma \mathbf{v},\mathbf{F}⟩$$

$$⟨\mathbf{v},\Lambda \mathbf{F}⟩$$

$$\Gamma \mathbf{v}$$

$$\Lambda \mathbf{F}$$

$$\mathbf{W}\left[\mathbf{v}\right]$$

$${\mathbf{W}}^{\ast }\left[\mathbf{F}\right]$$

Full text

Surface Fluctuating Hydrodynamics Methods for the Drift-Diffusion

Dynamics of Particles and Microstructures within Curved Fluid Interfaces

David A. Rower*∗*

  

Misha Padidar*∗*

  

and Paul J. Atzberger*∗*

W****e introduce fluctuating hydrodynamics approaches on surfaces for capturing the drift-diffusion dynamics of particles and microstructures immersed within curved fluid interfaces of spherical shape. We take into account the interfacial hydrodynamic coupling, traction coupling with the surrounding bulk fluid, and thermal fluctuations. For fluid-structure interactions, we introduce Immersed Boundary Methods (IBM) and related Stochastic Eulerian-Lagrangian Methods (SELM) for curved surfaces. We use these approaches to investigate the statistics of surface fluctuating hydrodynamics and microstructures. For velocity autocorrelations, we find characteristic power-law scalings $\tau^{-1}$, $\tau^{-2}$, and plateaus can emerge. This depends on the physical regime associated with the geometry, surface viscosity, and bulk viscosity. This differs from the characteristic $\tau^{-3/2}$ scaling for bulk three dimensional fluids. We develop theory explaining these observed power-laws associated with time-scales for dissipation within the fluid interface and coupling to the surrounding fluid. We then use our introduced methods to investigate a few example systems and roles of hydrodynamic coupling and thermal fluctuations including for the kinetics of passive particles and active microswimmers in curved fluid interfaces.

1 Introduction

Soft materials can exhibit rich mechanical responses arising from curved fluid interfaces that mediate interactions between immersed particles and other microstructures [30, 43, 23]. This includes protein interactions within lipid bilayer membranes [2, 43, 73, 4], surfactants and contaminants in bubble interfaces [26, 84], transport in soap films [45, 47], and recent systems with nanoparticles or colloids embedded in fluid interfaces [17, 54, 90, 24, 89, 19, 14, 61]. Related hydrodynamic and curvature mediated phenomena also play an important role in biology and physiology, including transport of surfactants in lung alveoli [42, 57] or in cell mechanics [60, 63, 70]. We develop general approaches to model and simulate the collective drift-diffusion dynamics of particles and other microstructures embedded within curved two-dimensional fluid interfaces. At small length and time scales, our methods also allow for capturing the roles of fluctuations arising from active microstructures and from thermal effects [81, 76, 74, 14, 15].

We develop methods for curved surfaces building on our prior work on stochastic immersed boundary methods and fluctuating hydrodynamics methods [8, 9, 82]. We also draw on our recent work on developing approaches for deterministic incompressible hydrodynamic flows on curved surfaces [36, 37]. Here, we address how to introduce the spontaneous thermal fluctuations and handle the associated drift-diffusive dynamics of microstructures both for hydrodynamics in the inertial regime and the overdamped quasi-steady regime. We develop theory and computational methods that capture for the fluctuations the correlations arising from the hydrodynamic coupling within the curved fluid interface and from the traction stresses with the flows of surrounding bulk fluids.

Many past approaches for investigating hydrodynamic coupling and diffusion have been based on the classic Saffman-Delbrück (SD) hydrodynamics model [79, 78]. These were derived for flat viscous sheets [85, 55, 65, 16, 22, 66]. However, for many problems arising in practice, the geometry plays an important role as a consequence of significant curvature on the SD length-scale or from the surface topology. These effects can significantly change the hydrodynamic responses relative to the flat case [82, 75, 36, 43]. This has motivated recent work going beyond the classic Saffman-Delbrück theory to take into account the role of geometry and additional mechanical effects arising in curved fluid interfaces [41, 5, 75, 80, 24, 82, 36].

We introduce general fluctuating hydrodynamics methods for the drift-diffusion of particles and microstructures immersed within curved fluid interfaces. To demonstrate ideas, we focus particularly on the case of interfaces of spherical shape. To ensure finite diffusivities, as indicated in the SD theory, we formulate hydrodynamic equations for curved fluid interfaces coupled to surrounding bulk fluids. We introduce thermal fluctuations accounting for the fluid-structure interactions and collective dynamics of the particles and microstructures building on our prior work on stochastic immersed boundary methods and related approaches in [8, 9, 82, 36, 37].

For curved fluid interfaces, the geometry and topology pose additional challenges for developing stochastic immersed boundary methods and numerical approaches. This includes formulating and developing methods to solve hydrodynamic equations within curved surfaces and to obtain appropriate fluctuations and fluid-structure coupling operators [82]. We develop stochastic numerical methods and approaches for fluid-structure coupling operators taking the geometry into account. These numerical solvers and coupling operators allow us to obtain mobility tensors for the collective hydrodynamic coupling for modeling both passive and active microstructures within the interface.

The geometry and topology of curved interfaces also pose additional challenges for generating incompressible hydrodynamic fields and thermal fluctuations. We develop techniques to handle these aspects of the fluid mechanics based on formulating generalized vector potentials using the Hodge decomposition for manifolds. For the overdamped quasi-steady regime, we further develop techniques for the mobility tensor by using the embedding space to derive covariance structures that capture consistently hydrodynamic correlations in fluctuating fields. The embedding space expands the dimension of the problem and can create null-spaces for some of the linear operators. To address this issue, we develop algorithms based on stabilizations for computing the stochastic driving forces in the drift-diffusion dynamics of the microstructures. Our introduced approaches provide methods for both inertial and overdamped quasi-steady regimes. Our methods provide approaches for capturing in simulations the interface hydrodynamics, traction coupling with the surrounding bulk fluid, fluid-structure coupling, and thermal fluctuations.

We organize our paper as follows. We formulate the fluctuating hydrodynamic equations and related immersed boundary methods for spherical fluid interfaces in Section 2. Stochastic numerical methods for the drift-diffusion dynamics of microstructures are discussed in Section 3. We use our approaches to investigate the statistical mechanics of surface fluctuating hydrodynamics and to develop theory for explaining observed power-laws in Section 4. We further demonstrate our approaches by investigating the role of hydrodynamic coupling and related diffusive correlations in the kinetics of passive particles and active microswimmers in Section 4. The results show some of the rich phenomena that can arise for microstructures and fluctuating hydrodynamics within curved fluid interfaces.

2 Fluctuating Hydrodynamics for Curved Fluid Interfaces

To formulate the hydrodynamics of curved fluid interfaces, we will first develop the conservation laws with reference to the ambient embedding space. We will then develop descriptions in terms of operators that generalize the techniques of vector calculus used in continuum mechanics [59, 36].

2.1 Conservation Laws for Curved Surfaces

The conservation of mass and momentum of the fluid can be expressed as

[TABLE]

The $\mathbf{v}_{n}$ and $\mathbf{v}_{\parallel}$ denote respectively the components of the fluid velocity normal and tangential to the fluid interface with $\mathbf{v}=\mathbf{v}_{\parallel}+\mathbf{v}_{n}$. The $H$ denotes the local mean curvature of the surface [71]. The $\rho$ is the local mass density, $\boldsymbol{\sigma}$ the internal interfacial stress, and $\overline{\mathbf{b}}$ the body force per unit mass.

V The $\overline{\mbox{div}}\left(\mathbf{t}\right)=\mathbf{t}_{|b}^{b}$ denotes the surface covariant divergence [1, 82, 36]. For a vector field represented as $\mathbf{t}=t^{a}\partial_{x^{a}}$ and $\mathbf{v}=v^{a}\partial_{x^{a}}$, the surface covariant derivative $\nabla_{\mathbf{v}}\mathbf{t}$ gives $\mathbf{w}=\nabla_{\mathbf{v}}\mathbf{t}=t_{|b}^{c}v^{b}\partial_{x^{c}}=w^{c}\partial_{x^{c}}$, with $w^{c}=t_{|b}^{c}v^{b}$. The $t_{|b}^{c}={\partial t^{c}}/{\partial x^{b}}+\Gamma_{ab}^{c}t^{a}$ and $\Gamma_{ab}^{c}$ denotes the Christoffel symbols [1, 71]. The material derivative of the mass on a surface is $\dot{\rho}=\partial\rho/\partial{t}+\rho\left(\overline{\mbox{div}}(\mathbf{v}_{\parallel})+\mathbf{v}_{n}H\right)$, where $H$ is the local mean curvature. The material derivative of the momentum is $\dot{\mathbf{v}}=\partial\mathbf{v}/\partial t+\nabla_{\mathbf{v}}\mathbf{v}$ [1, 36]. In the manifold setting, the material derivative of the vector field $\mathbf{t}$ can be expressed as $\dot{\mathbf{t}}=L_{\mathbf{v}}\mathbf{t}$, where $L_{\mathbf{v}}$ is the Lie derivative of $\mathbf{t}$ under the flow of the velocity field $\mathbf{v}$ [1, 59].

Throughout, we shall consider Newtonian incompressible fluid interfaces of fixed shape. In this case, the hydrodynamic flows are tangential to the surface. As a consequence, $\mathbf{v}_{n}=0$, $\rho=\rho_{0}$, $\dot{\rho}=0$ so that $\overline{\mbox{div}}\left(\mathbf{v}\right)=0$. For notational convenience, we suppress the $\mathbf{v}_{\parallel}$ and $\rho_{0}$ for the velocity and mass denoting them simply as $\mathbf{v}$ and $\rho$. The conservation laws simplify in this setting to

[TABLE]

The $\rho$ is now constant throughout and $\mathbf{b}$ denotes the local force density. The stress tensor $\boldsymbol{\sigma}$ for an incompressible Newtonian fluid for a curved fluid interface can be expressed as

[TABLE]

The $\mu_{m}$ is the surface viscosity of the interfacial fluid, $\mathbf{D}$ is the rate-of-deformation tensor, $p$ is the pressure, and $\mathcal{I}$ is the metric associated identity tensor [59].

Tensors can be expressed in covariant or contravariant form [40, 1]. For vectors in contravariant form we have $\mathbf{t}^{\sharp}=t^{a}\partial_{x^{a}}$ or in covariant form $\mathbf{t}^{\flat}=t_{a}\mathbf{d}x^{a}$ [1]. The $\partial_{x^{a}}$ denotes the $a^{th}$ coordinate basis vector. The $\mathbf{d}x^{a}$ denotes the $a^{th}$ coordinate basis covector (differential $1$-form) with $\mathbf{d}x^{a}[\partial_{x^{b}}]=\delta_{a}^{b}$ [1]. We can convert between vectors and covectors by the mappings $\flat:t^{a}\rightarrow t_{a}=g_{ab}t^{b}$ and $\sharp:t_{a}\rightarrow t^{a}=g^{ab}t_{b}$. The $g_{ab}$ denotes the metric tensor and $g^{ab}$ the inverse metric tensor [1]. The isomorphic maps $\flat$ and $\sharp$ between the tangent space and cotangent space correspond operationally in calculations to lowering and raising indices in the coordinate expressions of the tensors [40, 1].

We find it convenient in our calculations to express tensors in covariant form and use exterior calculus [1]. This allows us to generalize vector calculus and many techniques employed for fluid mechanics to the manifold setting [36, 37, 82]. For an incompressible Newtonian fluid interface, the stress $\boldsymbol{\sigma}$ is given in equation 8. The divergence of the stress tensor on the surface $\overline{\mbox{div}}\left(\boldsymbol{\sigma}\right)$ becomes in covariant form [36, 82, 59]

[TABLE]

The $K$ is the local Gaussian curvature of the surface [71]. The $\mathbf{d}$ is the exterior derivative playing a role similar to the gradient of the vector field on the surface. For a 1-form $\mathbf{v}^{\flat}=v^{i}d\mathbf{x}_{i}$ the exterior derivative is given by $\mathbf{d}\mathbf{v}^{\flat}=\partial v^{i}/\partial x^{j}\,d\mathbf{x}_{j}\wedge d\mathbf{x}_{i}$. The $\boldsymbol{\delta}=-\star\mathbf{d}\star$ is the co-differential playing a role similar to the divergence on the surface [1]. The $\star$ is the Hodge star which for a differential $k$-form $\boldsymbol{\beta}$ gives a complementary $n-k$-form $\star\boldsymbol{\beta}$ so that for any $k$-form $\boldsymbol{\alpha}$ we have $\boldsymbol{\alpha}\wedge\star\boldsymbol{\beta}=\langle\alpha,\beta\rangle\boldsymbol{\omega}$ where $\boldsymbol{\omega}$ is the volume form [1]. This also allows us to generalize vector calculus operators such as the curl and divergence to the surface by $\mbox{curl}_{\mathcal{M}}(\mathbf{v}^{\flat})=-\star\mathbf{d}\mathbf{v}^{\flat}$ and $\mbox{div}_{\mathcal{M}}(\mathbf{v}^{\flat})=\boldsymbol{\delta}\mathbf{v}^{\flat}$. We give coordinate expressions for these operations and additional discussions in Appendix D

2.2 Fluctuating Hydrodynamics with Fluid-Structure Interactions

We develop surface fluctuating hydrodynamics descriptions to account for the drift-diffusion motions of microstructures and their hydrodynamic coupling within curved fluid interfaces [8, 9, 82]. In covariant form, we introduce for curved fluid interfaces fluctuating hydrodynamic equations incorporating fluid-structure interactions. We illustrate our general approach in Figure 1.

The fluid dynamics are modeled by

[TABLE]

The drift-diffusion motions of microstructures are modeled by

[TABLE]

The $\mathbf{X}$ and $\mathbf{V}$ denotes the collective configuration and velocity of the microstructures immersed within the fluid. The $\mathbf{v}^{\flat}$ denotes the fluid velocity and the $-\boldsymbol{\delta}\mathbf{v}^{\flat}=0$ enforces the local incompressibility of the flow. The $\mu_{m}$ is the surface fluid shear viscosity, $p$ is the pressure, $K$ is the local Gaussian curvature, $\boldsymbol{t}^{\flat}$ is the traction stress with the surrounding fluid. The $-\nabla{\phi}$ are the conservative forces acting on the microstructures, and $\mathbf{f}_{thm}$ and $\mathbf{F}_{thm}$ are the stochastic forces accounting for thermal fluctuations of the system.

V The fluid-structure interactions result in two equal-and-opposite forces. The $-\gamma\left(\mathbf{V}-\Gamma\mathbf{v}^{\flat}\right)$ is a drag force of strength $\gamma$ the microstructures experience from local coupling to the fluid. The $\Lambda\left[\gamma\left(\mathbf{V}-\Gamma\mathbf{v}^{\flat}\right)\right]$ is the opposite force density the fluid experiences locally from coupling to the microstructures. The $\Lambda$ is a spreading operator that serves to convert a local force to a local force density acting on the fluid. The $\Gamma$ is an averaging operator that serves to estimate a local reference velocity for a microstructure from the nearby fluid flow. We discuss specific choices for $\Lambda$ and $\Gamma$ in more detail in Section 2.5. V Given the rapid oscillations of the fluid from the thermal fluctuations, we neglect the advection terms which are expected to give lower-order contributions in equation 2.2 [86].

2.3 Traction Stress from Flow of the Surrounding Bulk Fluid

To obtain the traction stresses $\mathbf{t}^{\flat}=\mathcal{T}_{f}\mathbf{v}^{\flat}$ for the spherical geometry, we use Lamb’s solution [39, 49]. The traction arises from the surface flow with velocity $\mathbf{v}$ entraining the surrounding bulk fluid [82]. Our approach makes the assumption that the bulk surrounding flow arises from an incompressible Newtonian fluid that reaches steady-state rapidly when contributing to the surface traction. Let the bulk fluid velocity be denoted by $\mathbf{u}$ with values on the surface $\mathbf{u}=\mathbf{v}+v_{n}\mathbf{n}$, where we shall assume $v_{n}=0$ throughout. The solution can be expressed using the spherical harmonics expansion

[TABLE]

We emphasize here the curl $\nabla\times$ is the usual operator in three dimensional Euclidean space. The $Z_{\ell}$ denotes the combined contributions of all of the solid spherical harmonic terms of degree $\ell$. We also expand the bulk surrounding fluid flows inside the sphere $\mathbf{u}^{-}$ and outside the sphere $\mathbf{u}^{+}$ as

[TABLE]

The $\mathbf{u}_{\ell}^{+}$ and $\mathbf{u}_{\ell}^{-}$ are the solid spherical harmonic expansion terms combined for a given degree $\ell$. The Lamb solutions for the bulk fluid velocity fields are given by [39, 49] expressible as

[TABLE]

where

[TABLE]

The traction stress expressed in contravariant form $\mathbf{t}^{\sharp}=\mathbf{t}^{+}+\mathbf{t}^{-}$ is

[TABLE]

The $\boldsymbol{\sigma}^{\pm}$ is the stress of the bulk surrounding fluid arising at the fluid interface. For the traction stress $\mathbf{t}^{\pm}=\mathbf{t}^{\pm}(\mathbf{x}_{0})$ on the surface, note the $\mathbf{n}^{\pm}=\mathbf{n}^{\pm}(\mathbf{x}_{0})$ is evaluated at the fixed location $\mathbf{x}_{0}$ while $\mathbf{u}^{\pm}=\mathbf{u}^{\pm}(\mathbf{x})$. When taking the gradient the $\mathbf{n}^{\pm}$ does not change, so that $\nabla(\mathbf{u}^{\pm}\cdot\mathbf{n}^{\pm})=\nabla\mathbf{u}^{\pm}\cdot\mathbf{n}^{\pm}$. We take $\mu_{\pm}=\mu_{f}$ throughout. The traction can be expressed in terms of the spherical harmonics expansion as

[TABLE]

The bulk fluid flow has $\mathbf{u}^{+}=\mathbf{u}^{-}=\mathbf{v}$ on the fluid interface and is completely determined by $\mathbf{v}$. In covariant form, we define the traction operators on the surface corresponding to equation 21 and 22 as $\mathbf{t}^{\flat}=\mathcal{T}_{f}\mathbf{v}^{\flat}=\mathcal{T}_{f}^{+}\mathbf{v}^{\flat}+\mathcal{T}_{f}^{-}\mathbf{v}^{\flat}$.

We remark that the approaches we develop can also be extended readily for surrounding fluids in systems that are subject to external shear flows or time-dependent flow responses. This would correspond to $\mathbf{t}^{\flat}$ obtained from equations 19-20 being computed from stresses of such flows from other model hydrodynamic equations or results from numerical solvers [67, 43, 92, 24, 9]. For such surrounding fluids, with flows having time-dependence or flows having non-Newtonian responses, our approaches can also be used capture these dissipitative contributions to the mechanics and develop the associated thermal fluctuations [9, 86].

2.4 Thermal Fluctuations

To determine the associated stochastic driving terms $\mathbf{f}_{thm}(\mathbf{x},t)$ and $\mathbf{F}_{thm}(t)$ that account for thermal fluctuations, we use an approach related to our Stochastic Eulerian Lagrangian Method (SELM) framework [9, 86]. This gives $\mathbf{f}_{thm}(\mathbf{x},t)$ and $\mathbf{F}_{thm}(t)$ as Gaussian processes with $\delta$-correlation in time, mean zero, and covariances

[TABLE]

The notation is to be interpreted as taking an expectation, $\langle Z\rangle=\mathbb{E}[Z]$. The equation 23 involves the dissipative operator $\mathcal{L}_{ff}\mathbf{v}^{\flat}$ associated with the fluid and is defined by

[TABLE]

The terms in equation 24 arise from the dissipative operator $\mathcal{L}_{pp}$ of the microstructure degrees of freedom and is defined by $\mathcal{L}_{pp}\mathbf{V}:=-\gamma\mathcal{I}\mathbf{V}$. The terms in equation 25 gives the cross-correlation that arises from microstructure-fluid coupling giving the operator $\mathcal{L}_{pf}$ and is defined by $\mathcal{L}_{pf}\mathbf{v}^{\flat}:=\gamma\Gamma\mathbf{v}^{\flat}$.

To express and generate efficiently the thermal fluctuations, it is useful to define a term $\mathbf{g}_{thm}$. This term is taken to be independent of $\mathbf{F}_{thm}$ and to have covariance

[TABLE]

We can then express the thermal fluctuations for the hydrodynamics as $\mathbf{f}_{thm}=\mathbf{g}_{thm}-\Lambda[\mathbf{F}_{thm}]$. This gives the correlations in equations 23 – 25.

In our formulation, we use throughout that the fluid-structure coupling operators are adjoints $\Lambda=\Gamma^{T}$, in the sense $\int_{\mathcal{S}}\langle\Lambda\mathbf{U},\mathbf{u}^{\flat}\rangle_{g}dA=\langle\mathbf{U},\Gamma\mathbf{u}^{\flat}\rangle$ for all choices of $\mathbf{U}$ and $\mathbf{u}^{\flat}$ [82, 9]. We discuss choices for the coupling operators and the adjoint conditions in more detail in Section 2.5. For the stochastic driving terms of equations 2.2– 13 and 23– 25, our SELM approach ensures the Gibbs-Boltzmann distribution is invariant under the stochastic dynamics and satisfies detailed-balance [9, 86, 77]. The stochastic equations 2.2– 13 should be interpreted in the sense of Ito Calculus [64, 32].

We remark that the covariance operators in equations 23–25 are to be interpreted in the weak sense [56, 11, 9]. Consider $\mathbf{u}(\mathbf{x})=\int\int\boldsymbol{\alpha}(\mathbf{z},r,\mathbf{x})^{T}\mathbf{f}_{thm}(\mathbf{z},r)d\mathbf{z}dr$ and $\tilde{\mathbf{u}}(\mathbf{y})=\int\int\tilde{\boldsymbol{\alpha}}(\mathbf{w},q,\mathbf{y})^{T}\mathbf{f}_{thm}(\mathbf{w},q)d\mathbf{w}dq$ and $\mathbf{U}=\int\mathbf{A}(r)^{T}\mathbf{F}_{thm}(r)dr$ and $\tilde{\mathbf{U}}=\int\tilde{\mathbf{A}}(q)^{T}\mathbf{F}_{thm}(q)dq$. The $\boldsymbol{\alpha}$, $\tilde{\boldsymbol{\alpha}}$, $\mathbf{A}$, $\tilde{\mathbf{A}}$ are smooth fields and vectors playing the role of test functions. [56]. The associated covariances are

[TABLE]

The differential operator $\mathcal{L}_{ff}$ acts in the parameter $\mathbf{w}$. We also have the covariances

[TABLE]

and

[TABLE]

The operator $\Gamma$ acts in the parameter $\mathbf{w}$. Additional discussions on how to interpret these operator covariances also can be found in [11, 9].

2.5 Fluid-Structure Coupling: Immersed Boundary Methods for Curved

Surfaces

We handle the fluid-structure interactions between the microstructures with the fluid by developing extended immersed boundary methods in the manifold setting, see Figure 1, [69, 8, 9]. Many choices can be made for the operators $\Gamma$ and $\Lambda$ [9]. To ensure that the approximate fluid-structure coupling which converts between the Eulerian and Lagrangian reference frames are non-dissipative, the operators are taken to be adjoints [9, 69, 86].

We require the coupling operators satisfy the following adjoint conditions for any choice of test field $\mathbf{v}$ and vector $\mathbf{F}$,

[TABLE]

where the inner-products are defined as

[TABLE]

The $\mathbf{X}$ denotes the collective vector of particle locations, or more abstract microstructure degrees of freedom [9]. For instance, the $i^{th}$ particle would be at location $\left[\mathbf{X}\right]_{i}$. The $\cdot$ denotes the vector dot-product induced by the ambient physical space. For vector fields on the surface $\mathbf{v}$ and $\mathbf{u}$ we have $\mathbf{v}(\mathbf{x})\cdot\mathbf{u}(\mathbf{x})=\langle\mathbf{v}(\mathbf{x}),\mathbf{v}(\mathbf{x})\rangle_{g}=v^{a}g_{ab}u^{b}$. We use the notation $\Gamma^{T}=\Lambda$ to denote the adjoint condition 32.

In the setting of curved surfaces, the coupling operators $\Gamma$ and $\Lambda$ should be chosen carefully to control the velocity averaging and force spreading. The $\Gamma$ should ensure the velocity averaging takes into account over the surface the different tangential directions at each location. The $\Lambda$ for force densities should have well-controlled components for the surface in the tangential and normal directions. We develop operators of the form

[TABLE]

We use a tensor $\mathbf{W}$ to sample and weight values on the surface to perform velocity averaging. We use the adjoint tensor $\mathbf{W}^{*}$ to produce a corresponding force density field on the surface compatible with our adjoint conditions 32. For the curved surface, we use the geometrically motivated forms

[TABLE]

The sum $i$ runs over the indices of the particle or microstructure and the $\partial_{x^{\alpha}}|_{\mathbf{X}^{[i]}}$ denotes the tangent basis vector in direction $\alpha$ at location $\mathbf{X}^{[i]}$. We refer to the vector field $\mathbf{w}^{[i],\alpha}$ as the probing vector field for direction $\alpha$.

On the sphere, the rotational symmetry can be utilized. For this case with the spherical coordinates $(\theta,\phi)$, we take our probing vector fields to be of the form $\mathbf{w}^{[i],\theta}=\psi(\mathbf{x}-\mathbf{X}^{[i]})\partial_{\theta}$ and $\mathbf{w}^{[i],\phi}=\left(\psi(\mathbf{x}-\mathbf{X}^{[i]})/\cos(\theta)\right)\partial_{\phi}$, where $\psi(r)=C\exp(-r^{2}/2\sigma^{2})$. In practice, we truncate at length $r_{0}=4\sigma$ and use $C$ to normalize so that $\psi(r)$ averages to one on the surface.

2.6 Overdamped Limit

In physical regimes with small Reynolds numbers $Re\ll 1$, we consider the overdamped limit of the fluctuating hydrodynamic equations 2.2– 12  [9, 86]. In this regime, the limiting fluctuating hydrodynamic equations can be expressed as

[TABLE]

where

[TABLE]

This is obtained by taking the overdamped limit while retaining the finite slip term $-\gamma\left(\mathbf{V}-\Gamma\mathbf{v}^{\flat}\right)$ in equation 2.2 and 12. This results in the slip term ${\gamma^{-1}}\mathcal{I}$ in equation 40 [86]. In the strong-coupling limit with $\gamma\rightarrow\infty$ the mobility simplifies to

[TABLE]

The velocity averaging operator $\Gamma$ and force spreading operator $\Lambda$ are as discussed in Section 2.5. The $\mathcal{S}$ denotes the solution operator $\mathbf{v}^{\flat}=\mathcal{S}\mathbf{b}^{\flat}$ for the hydrodynamic equations

[TABLE]

The $\mathcal{T}_{f}$ gives the traction stresses with the surrounding bulk fluid as discussed in Section 2.2. The $\mathbf{v}^{\flat}=\mathcal{S}\Lambda\mathbf{F}$ is the solution in the case with $\mathbf{b}=\Lambda\left[\mathbf{F}\right]$, where in practice we have $\mathbf{F}=-\nabla\phi$.

To simply expressions, we also take for equation 39 the limit of no-slip between the microstructure and the fluid, which corresponds formally to $\gamma/m\rightarrow\infty$. This yields $d\mathbf{X}/dt=\mathbf{V}=\Gamma\mathbf{v}^{\flat}$ throughout. Putting this together we arrive at the first term in equation 39. For the full stochastic system, we perform related detailed dimensional analysis and asymptotics in [9, 86, 82]. The equation 39 should be given the Ito interpretation. The thermal fluctuations involve configuration-dependent correlations resulting in the spontaneous drift term $k_{B}{T}\nabla\cdot\mathbf{M}$, [64, 32, 9, 86]. Additional analysis and reductions for different physical regimes can be found in [9, 86, 82].

2.7 Formulating Fluctuating Hydrodynamics on Surfaces using Vector

Potentials $\Phi$

The hydrodynamic responses both in the inertial fluctuating hydrodynamics of equation 2.2 and in the overdamped regime of equation 42 require computation of the operator $\mathcal{L}_{f}$ of equation 27. Since the hydrodynamic fields are incompressible, we can use the Hodge decomposition of the fluid velocity $\mathbf{v}^{\flat}=\mathbf{d}\tilde{\Phi}+\boldsymbol{\delta}\Psi+h$. The $\tilde{\Phi}$ and $\Psi$ are vector potentials and $h$ is a harmonic function related to the topology of the manifold [59, 1]. The $\Delta_{H}=\mathbf{d}\boldsymbol{\delta}+\boldsymbol{\delta}\mathbf{d}$ is the Hodge Laplacian and the $h$ is harmonic in the sense $\Delta_{H}h=0$. For the case of spherical geometry with tangential velocity fields, we can express incompressible flows as $\mathbf{v}^{\flat}=-\star\mathbf{d}\Phi$, [82]. This uses that incompressibility requires $\boldsymbol{\delta}\mathbf{v}^{\flat}=0$ and from the spherical geometry $h=0$. This simplifies the action of the Hodge Laplacian and the representations. This allows for reformulating the inertial hydrodynamics of equation 2.2 in terms of $\Phi$ as

[TABLE]

The $\mathbf{b}^{\flat}=\mbox{curl}_{\mathcal{M}}(c)=-\star\mathbf{d}c$ denotes the contributions of the forces acting on the fluid. The $\Delta_{LB}=-\boldsymbol{\delta}\mathbf{d}$ denotes the Laplace-Beltrami operator [1]. The $\mathcal{T}_{f}$ denotes the traction stress operator of equation 19 and 20. The $\mathcal{L}_{f}$ denotes the operator in equation 27 for the hydrodynamic response.

We obtained equation 2.7 by substituting for $\Phi$, using $\mathbf{v}^{\flat}=-\star\mathbf{d}\Phi$, and using on both sides the generalized $\mbox{curl}_{\mathcal{M}}(\cdot)=-\star\mathbf{d}(\cdot)$. We also use that the generalized curl commutes with the operator $\mathcal{L}_{f}$ in the sense $(-\star\mathbf{d})\mathcal{L}_{f}(-\star\mathbf{d})=\Delta_{LB}\tilde{\mathcal{L}}_{f}$. The $\mathcal{\tilde{L}}_{f}$ takes on a form similar to equation 27 but now applied to scalar fields. We provide more details for $\mathcal{\tilde{L}}_{f}$ in terms of spherical harmonic coefficients below. These considerations allow for the fluctuating hydrodynamics to be expressed as

[TABLE]

In the overdamped regime, we need to compute the mobility tensor $M(\mathbf{X})$ of the hydrodynamic coupling in equation 40. This requires us to solve the steady-state hydrodynamic equations 42. This can be reformulated as

[TABLE]

We again use that the generalized curl commutes with the operator $\mathcal{L}_{f}$ in the sense $(-\star\mathbf{d})\mathcal{L}_{f}(-\star\mathbf{d})=\Delta_{LB}\tilde{\mathcal{L}}_{f}$. This yields for the steady-state hydrodynamics

[TABLE]

Central to both the inertial fluctuating hydrodynamics and steady-state regimes is computation of the action of the operator $\tilde{\mathcal{L}}_{f}\Phi$ which gives the hydrodynamic responses of the system. We expand $\Phi$ using spherical harmonics to represent the vector potential as $\Phi=\sum_{s}\Phi_{s}=\sum_{s}a_{s}Y_{s}$ where $\Phi_{s}=a_{s}Y_{s}$ and $s=(\ell,m)$ denotes the spherical harmonics index, see Appendix C. For the spherical geometry, we can express solutions using that the spherical harmonics are eigenfunctions of the Laplace-Beltrami operator [82, 36].

The action of the operator can be expressed as

[TABLE]

where

[TABLE]

with $L^{\pm}=\mu_{m}/\mu_{\pm}$. In the case with $\mu_{\pm}=\mu_{f}$ this can be simplified to

[TABLE]

where $L=\mu_{m}/(\mu_{+}+\mu_{-})=\mu_{m}/2\mu_{f}$. We see from this that the hydrodynamic flow on the sphere is characterized by the non-dimensional parameter $\Pi_{1}=L/R$, where $L=\mu_{m}/2\mu_{f}$ is the Saffman-Delbrück length-scale [79] and $R$ the radius of the sphere.

In the overdamped regime, to obtain the steady-state hydrodynamics, we can solve the hydrodynamic equations to obtain directly the spherical harmonics coefficients as

[TABLE]

For the inertial fluctuating hydrodynamics of equation 2.2 and equation 2.7, we can express the dynamics in terms of the spherical harmonics coefficients $a_{s}$ as

[TABLE]

where $c_{s}=\bar{c}_{s}+g_{s}+h_{s}$. The stochastic driving terms $\{g_{s}\}$ are associated with the hydrodynamics. These are complex-valued Gaussian processes $\delta$-correlated in time with mean zero and covariance

[TABLE]

The stochastic driving terms $\{h_{s}\}$ are associated with the fluid-structure coupling discussed in Section 2.5. This involves the contributions from the dissipative terms in equation 2.2. The $h_{s}$ are complex-valued Gaussian processes having correlations with the microstructure stochastic dynamics in equation 12. We can express this as

[TABLE]

The $\langle\cdot,\cdot\rangle_{L^{2}}$ gives the surface $L^{2}$-inner product. The $\lambda_{s}=-\ell(\ell+1)/R^{2}$ denotes the $s^{th}$ eigenvalue of the Laplace-Beltrami operator $\Delta_{LB}=-\boldsymbol{\delta}\mathbf{d}$. The $\delta_{a,b}$ denotes the Kronecker $\delta$-function which is zero for $a\neq b$ and one for $a=b$.

The real-space fluctuating velocity field $\mathbf{v}^{\flat}$ can be recovered using $\mathbf{v}=(-\star\mathbf{d}\Phi)^{\sharp}$. This is derived from the stochastic force $\mathbf{F}_{thm}$ acting on microstructures in equation 12. Together, these stochastic driving terms give in real-space $\mathbf{f}_{thm}^{\flat}=-\star\mathbf{d}\sum_{s}\rho(g_{s}+h_{s})Y_{s}=\mathbf{g}_{thm}^{\flat}-\Lambda[\mathbf{F}_{thm}]$ which satisfy equations 23–25.

In summary, the equations 2.7–  52 give a formulation for surface fluctuating hydrodynamics incorporating fluid-structure interactions in terms of the vector potential $\Phi$. The use of such Hodge decompositions and vector potentials is particularly convenient when representing velocity fields on surfaces. We give additional details on the derivation of these results in Appendix A.

3 Stochastic Numerical Methods

We develop numerical methods for the drift-diffusion dynamics of microstructures both in the inertial and overdamped regimes. To handle the incompressibility constraint for the surface hydrodynamics we make use of the generalized vector potential formulation discussed in Section 2.7.

3.1 Inertial Regime

In the inertial regime, we formulate the time-step integration of the fluctuating hydrodynamics of equation 2.2 in terms of stochastic dynamics of the vector potential $\Phi$ in equation 2.7 and 2.7. We use a time-step update in terms of the spherical harmonics coefficients of $\Phi$ given by

[TABLE]

For the microstructure dynamics, we use an update

[TABLE]

We have $\Phi^{n}=\sum_{s}a_{s}^{n}Y_{s}$ with the fluid velocity obtained by $\mathbf{v}^{n}=(-\star\mathbf{d}\Phi^{n})^{\sharp}$. Similarly, from the fluid-structure coupling $\mathbf{c}^{n}=\Lambda^{n}[\gamma(\mathbf{V}^{n}-\Gamma^{n}\mathbf{v}^{n}]$ we obtain the terms $c_{s}^{n}=\left\langle\Delta_{LB}^{-1}(-\star\mathbf{d})\mathbf{c}^{n},Y_{s}(\mathbf{x})\right\rangle_{L^{2}}$. The hydrodynamic time-step updates in equation 3.1 are related to the Euler-Maruyama method [46] and the microstructure updates in equation LABEL:eqn:verlet_update are related to modified velocity-verlet methods [88, 29, 38]. This gives in the stochastic setting errors with worse-case scalings $O(\sqrt{\Delta{t}})$, [88, 29, 38].

We account for thermal fluctuations of the microstructures over the time-step through the Gaussian term $\mathbf{F}_{thm}^{n}$. This has mean zero and covariance

[TABLE]

The $g_{s}^{n}$ are complex-valued Gaussians with mean zero and covariance

[TABLE]

The $L_{s}$ are given by equation 48 or 2.7 and $\mathcal{C}_{ss^{\prime}}$ are the equilibrium fluctuations of the modes $a_{s}$ given by

[TABLE]

The $\lambda_{s}=-{\ell(\ell+1)}/{R^{2}}$ are the eigenvalues of the Laplace-Beltrami operator $\Delta_{LB}=-\boldsymbol{\delta}\mathbf{d}$ on the surface. For $s=(\ell,m)$ the $\overline{s}=(\ell,-m)$.

The $h_{s}^{n}$ are complex-valued Gaussians which we generate from $\mathbf{F}_{thm}^{n}$ by

[TABLE]

This ensures proper correlations in the system between the microstructure and the fluctuating hydrodynamics. We note that in our derivation of $\mathbf{F}_{thm}$ there were negative cross-correlations with the fluid which we capture in our numerical methods consistent with equation 25. These cross-correlations play the important role of ensuring momentum conservation. They account for the spontaneous fluctuations that exchange momentum back and forth between the fluid and microstructures.

The integrator approach we have introduced works well when the time-scales are comparable between the microstructure evolution and hydrodynamics. When there is a disparity in these time-scales, stiff stochastic numerical time-step integrators can also be developed using the vector potential formulation. For instance, we can develop exponential time-stepping approaches similar to our prior works [8, 93]. As an alternative, an asymptotic analysis of the stochastic dynamics of the fluid-structure system can also be performed as in [86]. This can be used to formulate equations in overdamped regimes and develop stochastic numerical methods [86].

3.2 Overdamped Regime

For the overdamped regime, we develop numerical methods for simulations based on stochastic variants of the velocity-verlet method [88, 29, 38]. We account for the thermal drift term in the stochastic dynamics using an approach related to methods in [21, 27]. We update the collective configuration $\mathbf{X}$ of the particles or microstructures of the system using

[TABLE]

The thermal fluctuations have correlations generated by ${Q}(\mathbf{X}^{n})$, where ${Q}(\mathbf{X}^{n}){Q}(\mathbf{X}^{n})^{T}=2k_{B}{T}M(\mathbf{X}^{n})/\Delta{t}$ $\boldsymbol{\xi}$. The $\boldsymbol{\xi}$ are standard Gaussian random variates with independent components having mean zero and variance one. The thermal drift $k_{B}{T}\nabla\cdot\mathbf{M}$ is approximated by numerically estimating an average by $\langle Z\rangle_{\bar{N}}=1/\bar{N}\sum_{k=1}^{\bar{N}}Z^{[k]}$. For the random variable $Z$, the $Z^{[k]}$ denotes one of the $\bar{N}$ independent samples. This provides a probabilistic estimator for the divergence. The random variables $\mathbf{p}$ and $\mathbf{q}$ satisfy $\langle{p}_{i}{q}_{j}\rangle=\delta_{ij}$. This gives

[TABLE]

This follows an approach related to [21, 27]. The validity follows readily by Taylor expanding $M$ in the variable $\mathbf{X}$ [33].

The stochastic velocity-verlet scheme in equation 59 can be viewed in stages as follows. The $\mathbf{\tilde{X}}^{n+1}$ gives the predictor part of the update of the configuration that is used to evaluate the force $\mathbf{\tilde{F}}^{n+1}$ and mobility $M(\mathbf{\tilde{X}}^{n+1})$. The $\mathbf{X}^{n+1}$ gives the corrector part for the update of the configuration which makes use of the predictor data and additional contributions from the thermal drift term.

3.2.1 Generating Stochastic Forces from the Mobility Tensor

During numerical time-step integration of the drift-diffusion dynamics by equation 59, we must generate the stochastic driving term $\mathbf{h}^{n}=Q(\mathbf{X}^{n})\boldsymbol{\xi}^{n}$. One way to do this is to determine from the mobility the square-root factor $Q$ satisfying $Q(\mathbf{X}^{n})Q(\mathbf{X}^{n})^{T}=2k_{B}{T}M$. Cholesky factorization methods can be used for symmetric positive definite matrices [87]. However, in the current formulation this present challenges given that the mobility tensor for the surface hydrodynamics is singular since it only depends on the tangential component of the force, see equation 2.7. Working only with the surface coordinates also presents challenges since there is no global coordinate chart for spherical topologies [82]. We instead work in the embedding space, but this results in a representation for the mobility tensor that is only positive semi-definite preventing direct use of the Cholesky factorization method.

We remedy this situation by adding a block diagonal tensor to our grand mobility tensor $M$ of equation 41. We do this in a way that preserves the drift-diffusion dynamics in the tangential directions. This is accomplished by using rank-one stabilizations to obtain modified self-mobility blocks of the form $\tilde{M}_{ii}=M_{ii}+\alpha\mathbf{n}_{i}\mathbf{n}_{i}^{T}$. The $\alpha>0$ is any chosen positive weight. The $\mathbf{n}_{i}=\mathbf{n}(\mathbf{X}_{i})$ denotes the normal vector on the sphere at the location of $\mathbf{X}_{i}$. By using this modification, we can obtain a tensor that is positive definite. This allows us to use Cholesky factorizations to obtain $\tilde{Q}$ satisfying $\tilde{Q}(\mathbf{X}^{n})\tilde{Q}(\mathbf{X}^{n})^{T}=2k_{B}{T}\tilde{M}$.

Our introduced stabilizations are constructed orthogonal to the tangent directions of the surface. As a result, this only modifies the mobility responses in the directions normal to the surface and preserves the microstructure drift-diffusion dynamics in the tangential directions. We generate the tangential stochastic forces using $\mathbf{h}=\wp\tilde{\mathbf{h}}^{n}$, where $\tilde{\mathbf{h}}^{n}=Q(\mathbf{X}^{n})\boldsymbol{\xi}^{n}$. The $\boldsymbol{\xi}^{n}$ are Gaussian vectors with independent components having mean zero and variance one. The $\wp$ denotes projection of the stochastic vector $\tilde{\mathbf{h}}^{n}$ to the tangent space of the surface. This allows us to generate the stochastic terms for the thermal fluctuations in the surface drift-diffusion dynamics of the microstructures.

3.3 Computing Hydrodynamic Forces and Responses

To solve for the fluid velocity, we represent the force density term as $c=\mbox{curl}_{\mathcal{M}}(\mathbf{b}^{\flat})=-\star\mathbf{d}\mathbf{b}^{\flat}$ with expansion $c=\sum_{\ell}c_{s}Y_{s}$. We compute this in practice by using the spherical harmonics representation of $\mathbf{b}$ obtained by $\hat{b}^{[a]}_{s}=\langle{b}^{[a]},Y_{s}\rangle$. The $\mathbf{b}={b}^{[a]}\partial_{z^{a}}$ gives the representation of the force density $\mathbf{b}$ in the ambient physical space with coordinates $z^{a}$. The inner-products $\langle\cdot,\cdot\rangle$ are approximated numerically by using Lebedev quadrature [53, 52]. We then compute $c=\mbox{curl}_{\mathcal{M}}(\mathbf{b}^{\flat})=-\star\mathbf{d}\mathbf{b}^{\flat}=\sum_{s}(-\star\mathbf{d})\mathbf{b}_{s}$. This uses the finite spherical harmonics expansion of $\mathbf{b}=\sum_{s}\mathbf{b}_{s}$, where $\mathbf{b}_{s}$ denotes the term of the expansion associated with index $s$. We compute $c_{s}$ using ${c}_{s}=\langle c,Y_{s}\rangle$, where the inner-product is approximated numerically using Lebedev quadrature [53, 52]. From the solution coefficients $a_{s}$, we obtain the fluid velocity on the surface from $\mathbf{v}=\sum_{s}a_{s}(-\star\mathbf{d}Y_{s})^{\sharp}$.

In practice, we numerically approximate $L^{2}$ inner-products by using Lebedev quadratures on the spherical surface [53, 52, 82]. We use $\langle u,v\rangle=\langle u,v\rangle_{Q}=\sum w_{i}u(\mathbf{x}_{i})v(\mathbf{x}_{i})$, where $\langle u,v\rangle_{Q}$ denotes the quadrature computed inner-product with $w_{i}$ the weights and $\mathbf{x}_{i}$ the nodes [53, 52, 82]. This provides a finite spherical harmonics expansion of $\mbox{curl}_{\mathcal{M}}(\mathbf{b})=-\star\mathbf{d}\mathbf{b}$. While other quadrature approaches could be used such as latitude and longitude sampling, the Lebedev quadrature nodes have octahedral symmetry and provide a better distribution of sampling nodes on the spherical surface [36, 82].

We remark that for the hydrodynamics on the sphere most of the representation of the flow responses and evolution is represented analytically using the spherical harmonics expansions. As a result, in the developed numerical methods the primary source of errors is from the truncation of the spherical harmonics expansion and related projections. In these methods, the primary source of computational expense arises from computing the $L^{2}$-inner products using the quadratures and in summing the spherical harmonics expansions in reconstructions.

When using $n$ modes, this part of the numerical algorithms have a computational complexity of $O(n^{2})$. In the inertial regime, the thermal fluctuations requires initially a single off-line computation of the Cholesky factorization of the covariance which costs $O(n^{3})$. Generating the correlated variates then costs $O(n^{2})$ each time-step. In the overdamped case, the primary additional computational expense is from the Cholesky factorization of the mobility tensor. For $m$ microstructure degrees of freedom and $n$ modes, the overdamped case has computational complexity $O(n^{2}+m^{3})$.

The Lebedev quadratures allow for computing the $L^{2}$-inner products to a high level of accuracy. Available Lebedev quadratures allow for $L^{2}$ inner-products for functions expanded up to degree $65$ spherical harmonics to be computed up to round-off errors  [52]. We performed convergence analysis of our spectral numerical methods using Lebedev quadratures for approximating exterior calculus operators and solving hydrodynamic equations in [35, 36].

4 Applications

We discuss a few applications showing how to use our introduced fluctuating hydrodynamics methods for investigating phenomena within curved fluid interfaces. We first discuss the hydrodynamic relaxation of fluctuating fluids and characterize scalings of the velocity autocorrelations. We find different scalings can emerge depending on the physical regimes associated with the interface geometry, surface viscosity, and bulk viscosity. We further investigate the mobilities associated with microstructures embedded within fluid interfaces in the quasi-steady regime. We then demonstrate the computational methods by studying the correlated diffusion of passive particles and the drift-diffusion dynamics of active microswimmers. The results show some of the rich phenomena within curved fluid interfaces that can arise from hydrodynamic coupling, thermal fluctuations, and geometry that our methods can be used to investigate.

4.1 Autocorrelations of the Surface Fluctuating Hydrodynamics

We investigate the autocorrelation of the fluid velocity of interfacial fluctuating hydrodynamics introduced in equation 2.2. We find the velocity autocorrelations can exhibit power-law decay with scalings $\tau^{-1}$ and $\tau^{-2}$ depending on the physical regime. We also find in some regimes a plateau behavior can arise. This differs from bulk three dimensional fluctuating hydrodynamics that would exhibit $\tau^{-3/2}$ power-law decay. This also differs from purely two dimensional fluctuating hydrodynamics that would exhibit only a $\tau^{-1}$ power-law decay [7, 8].

Interfacial fluctuating hydrodynamics involves both dissipation from the propagation of shearing motions within the interfacial fluid surface and from traction coupling with the bulk surrounding fluid. This results in two important time-scales. The first is the time-scale for shear stresses to propagate over the entire spherical surface $\tau_{f}=R^{2}/\mu_{m}$. The second is the time-scale on which rigid-body rotation of the entire spherical interface dissipates energy significantly to the bulk surrounding fluid $\tau_{r}=RL/\mu_{m}$.

We explore these contributions to the velocity autocorrelations by varying the ratio $\tau_{f}/\tau_{r}=R/L$. This involves the sphere radius $R$ and the Saffman-Delbrück length $L$. We show our results in Figure 3. We take as default parameter values $\rho=1$, $\mu_{m}=1$, $K_{B}{T}=1$, $R=1$, $L=1$, and $\ell\leq N=\ell_{*}=50$. We give details on our derivations analyzing these autocorrelations in Appendix B.

The velocity autocorrelations of the surface fluctuating hydrodynamics can be expressed using the following expansion in spherical harmonics

[TABLE]

The $v^{\phi}$ denotes the $\phi$-directional component of the velocity and the $Y_{s}$ denotes the spherical harmonic with mode $s=(\ell,m)$. The metric factor on the sphere is $|g|=R^{2}$. We also use that the velocity is isotropic and that $\langle a_{s}^{2}\rangle=k_{B}TR^{2}/\rho\ell(\ell+1)$ and $\langle a_{s}(0)a_{s}(t)\rangle=\langle a_{s}^{2}\rangle\exp(tL_{s})$, where $L_{s}$ is given in equation 48. We derive scaling laws for different physical regimes using equation 61. We remark that all of our results use a truncated expansion with modes up to degree $N=50$. This introduces a regularization length-scale similar to our immersed boundary coupling approaches for determining responses of particles and microstructures. Details of our derivations can be found in Appendix B.

When $R/L\ll 1$, we find the velocity correlations exhibit an initial decay to a power-law regime with scaling $\tau^{-1}$ on time-scales $\tau\ll\tau_{f}$. This is followed by a plateau regime that persists from $\tau_{f}\ll\tau\ll\tau_{r}$. This eventually gives way to exponential decay for $\tau\gg\tau_{r}$. The initial $\tau^{-1}$ power law decay is associated with the propagation of shear stresses over the interface with negligible dissipation to the bulk. This persists until time-scale $\tau_{f}$, see Figure 3.

We remark that this is similar to the velocity autocorrelations that would be observed in a purely two dimensional viscous Newtonian fluid, which can be computed readily using the methods in [7]. In contrast, for the spherical fluid interface the two dimensional fluid has finite area. This results in finite size effects in the flow responses and correlations.

A plateau arises in the case when $\tau_{r}\gg\tau_{f}$ which creates an intermediate regime. In the intermediate plateau regime, the shear stress has already propagated to the entire surface, but the dissipation into the bulk fluid of the rigid-body rotational motion has not yet become significant. On time-scales $\tau\gg\tau_{r}$, the dissipation into the bulk fluid dominates through the rotational motions and gives exponential decay. As the $R/L$ increases, the plateau regime disappears when $\tau_{r}\sim\tau_{f}$, as seen when moving left to right in Figure 3.

In the $R/L\gg 1$ regime, the time-scales for decay from intra-interface shearing motions becomes reversed and large relative to the time-scale for decay from the coupling to the bulk fluid $\tau_{f}\gg\tau_{r}$. This results in a new regime with power law scaling $\tau^{-2}$ for time-scales $\tau\ll\tau_{r}$. When $\tau\gg\tau_{f}$ this again eventually gives exponential decay from the rotational motions of the entire interface. The $\tau^{-2}$ arises from simultaneous dissipation from the shearing motions within the interface and dissipation from coupling to the bulk fluid. We give more details and derivations of these power-laws and related time-scales in Appendix B.

4.2 Mobility Tensor for Interacting Particles and Microstructures

We compute numerically the mobility $\tilde{M}(\mathbf{X})=\Gamma\tilde{\mathcal{S}}\Lambda$ using equation 41. We obtain a numerical solver $\tilde{\mathcal{S}}$ for the fluid velocity $\tilde{\mathbf{v}}^{\flat}=\tilde{\mathcal{S}}\Lambda[\mathbf{F}]$ with force distribution $\mathbf{f}^{\flat}(\mathbf{x})=\Lambda[\mathbf{F}](\mathbf{x})$ on the spherical surface using equation 2.7 and Lebedev quadratures [53, 52, 82]. We performed convergence analysis of the numerical methods based on Lebedev spectral approaches in [35, 36].

From the symmetry of the sphere the mobility of $N$ particles can be determined from the canonical two particle mobility tensor. The two particle mobility obtained from equation 41 can be expressed as

[TABLE]

with

[TABLE]

The $M_{ij}=\Gamma_{i}\tilde{\mathcal{S}}\Lambda_{j}$ where we use notation $\Gamma_{i}\mathbf{v}^{\flat}=\Gamma[\mathbf{X}_{i}]\mathbf{v}^{\flat}$ and $\Lambda_{j}[\mathbf{F}](\mathbf{x})=\Lambda[\mathbf{X}_{j}][\mathbf{F}](\mathbf{x})$. The $M_{ii}$ denote the self-mobilities and can be computed numerically once and stored.

The $\bar{M}_{12}$ give for a force $\mathbf{F}$ on the first particle the velocity response $\mathbf{V}$ at the second particle. We can express this as $\mathbf{V}=V_{\parallel}\mathbf{e}_{\parallel}+V_{\perp}\mathbf{e}_{\perp}$, where for the two particle displacement we split into the parallel $\parallel$ and perpendicular $\perp$ components. The $V_{\parallel}=M_{\parallel\parallel}\mathbf{F}_{\parallel}+M_{\parallel\perp}\mathbf{F}_{\perp}$ and $V_{\perp}=M_{\perp\parallel}\mathbf{F}_{\parallel}+M_{\perp\perp}\mathbf{F}_{\perp}$. Using the symmetry of the sphere we can tabulate numerically the two particle mobility tensor by using a canonical configuration of the two particles. We rotate the sphere so that the first particle is always situated at the north pole. We then align the geodesic displacement between the two particles in the $xz$-plane with tangent along the positive $x$-axis. We denote the rotation operation by $\mathcal{R}^{T}$ that moves any two particles into this canonical configuration $\tilde{\mathbf{X}}=\mathcal{R}^{T}\mathbf{X}$. We can convert our canonical tabulated two-particle mobilities to the specific mobility of two particles by $M_{12}=\mathcal{R}\tilde{M}_{12}\mathcal{R}^{T}$. We obtain the grand-mobility tensor $M$ for $n$ interacting particles by summing over all pairs of particles the two particle mobility tensors. We show the components of the two-particle mobility $\bar{M}_{12}$ when $L/R=0.13,0.65,6.5,52$ in Figure 4 and 5.

4.3 Equilibrium Fluctuations

We validate our stochastic numerical methods by studying the drift-diffusion dynamics of interacting particles. A strong indication of the validity the methods is provided by consider how particles diffuse over time when subject to a conservative force. This requires the stochastic methods to capture both the drift dynamics accurately while also handling appropriately the thermal fluctuations of the system. From equilibrium statistical mechanics a configuration should have a probability distribution $\rho(\mathbf{X})$ of the Gibbs-Boltzmann form [77, 18]

[TABLE]

The $Z$ denotes the partition function and $k_{B}{T}$ denotes the thermal energy of the system [77].

We consider the drift-diffusion dynamics of two hydrodynamically coupled particles having the non-linear Lennard-Jones interaction

[TABLE]

The $r$ denotes the distance between the two particles. The $\sigma$ denotes the length-scale characterizing the radius of the particles. We find our stochastic numerical methods provide for the drift-diffusion dynamics very good agreement with the distribution predicted from equilibrium statistics mechanics from equation 68. We show these results for a few choices of $\sigma$ in Figure 6.

For parameters, we take throughout the thermal energy $k_{B}T=2.48$ amu$\cdot$nm2/ps2, the strength of the potential $\epsilon=10K_{B}T$, viscosity ratio $L/R=0.65$. We use for the stochastic integrator in equation 59 the time-steps $\Delta{t}=1.3\times 10^{5}$ ps and drift estimator $\delta=10^{-1}$, $\bar{N}=10$.

4.4 Hydrodynamic Correlations in Particle Diffusion

Motivated by proteins in lipid bilayer membranes and recent results for synthetic colloidal systems [20, 31, 17, 61], we consider diffusion limited interactions by particles that interact within a curved fluid interface. We investigate the case with and without hydrodynamic coupling on the collective diffusivity and how this influences the time for particles to come into near contact.

We consider two particles initially at antipodal locations on the sphere (north and south pole) and the amount of time it takes for them to come into contact with each other. We consider the distribution for this meeting time in the case of hydrodynamic coupling with mobility $M$ as in Section 4.1 and in the case without hydrodynamics with a local drag having mobility response $M=-\gamma\mathcal{I}$. We report these results with and without hydrodynamic interactions in Figure 8. We use in our studies the parameter values in Table 1.

The $R$ gives the radius of the spherical fluid interface, the $R_{particle}$ gives the effective size of the particles used to specify the radius of separation for considering contacts, $D$ is the particle diffusivity, and $\gamma_{slip}$ is the slip permitted between the particle and the local fluid. The $L^{\pm}=\mu_{m}/2\mu_{\pm}$ gives the Saffman-Delbrück lengths between the surface fluid viscosity $\mu_{m}$ and the surrounding bulk fluid viscosity $\mu_{-}$ inside and $\mu_{+}$ outside the spherical fluid interface. The $k_{B}{T}$ gives the thermal energy.

We find in both cases that the meeting times remain on average on the same order of magnitude. However, the hydrodynamic coupling introduces significantly more variation in the meeting-time distribution producing a long-tail, see Figure 8. For systems where diffusive kinetics are important, we see the hydrodynamic coupling can significantly influence the distribution of encounters between particles.

4.5 Microscopic Swimmers and Mixing

We investigate hydrodynamic transport and diffusive mixing associated with swimmers at microscopic scales. Behaviors both individually and collectively can differ significantly from macroscopic scales [51, 72]. Swimming in both Newtonian and Non-Newtonian bulk fluids in three dimensional volumes have been investigated in [51, 72, 62]. We investigate here the case of swimmers within two dimensional curved fluid interfaces treated as a Newtonian fluid. Our approaches also could be used to study non-Newtonian fluids either by incorporating explicitly the microstructures in the fluid using approaches of Section 2.5 and [10], or by extending the formulation of the hydrodynamic equations to other constitutive laws in Section 2.1.

As a demonstration of the methods, we consider Golestanian Swimmers [62] that consist of three beads that interact through two oscillating harmonic bonds. The harmonic bonds have time-dependent rest-lengths with energy

[TABLE]

where $i=0,2$. The bond lengths are offset in time to have different phases for $\phi_{0}$ and $\phi_{2}$. To impose excluded volume, both the beads of the microscopic swimmers and the passive particles also interact through the Weeks-Chandler-Andersen (WCA) potential [91]

[TABLE]

The separation distance is denoted by $r=\|\mathbf{X}_{i}-\mathbf{X}_{j}\|$ for two particles $\mathbf{X}_{i}$ and $\mathbf{X}_{j}$. The potential gives an effective particle steric radius of $r_{s}=2^{1/6}\sigma_{wca}$. We give parameters for our model in Table 2. We remark that this parameterization is for illustrative purposes of the methods and to obtain more physically realizable systems may require further adjustments.

The phase differences in the swimming strokes are crucial to break symmetry in time to have the possibility of a net forward motion. For steady-state hydrodynamics, a time-reversible motion would have no net displacement by the Scallop Theorem [72, 50, 44]. We also emphasize that without hydrodynamic coupling between the beads the swimmer would remain stationary. This is a consequence of the equal-and-opposite forces that act on the beads and average out to zero over the periodic strokes. Without the surrounding fluid, such forces can not move the center-of-mass of the swimmer.

The swimmer strokes generate flows that on average pump the fluid. We show the associated stages of the stroke cycle captured by our methods in Figure 9. The hydrodynamics is confined to a surface of spherical topology which results in flows with vortices. We show the average flow that pumps the fluid in the inset of Figure 10.

We study first how the swimmer’s speed depends on the viscosities of the two dimensional interfacial fluid and surrounding bulk fluid. We characterize the viscosities by the ratio $L/R$ of the Saffman-Delbrück length $L$ to the sphere radius $R$ as discussed in Section 2.7. We can consider the swimmer’s angular progression over the sphere from which the swimming speed can be estimated, see Figure 10. As the viscosity ratio $L/R$ increases, the generated flows transition from being relatively localized to enveloping most of the sphere. We find as the viscosity ratio increases from $L/R=0.13$ to $L/R=52$ that the swimming speed drops by approximately $\sim 25\%$, see Table 3.

We next investigate the collective drift-diffusion dynamics of $N_{T}$ passive particles when subjected to mixing by multiple microscopic swimmers $N_{S}$. Both the swimmers and particles are subjected to the hydrodynamic coupling and thermal fluctuations using the approaches we introduced in Section 2. We study as the number of microscopic swimmers increases how the effective diffusivity of the passive tracer particles is influenced.

We characterize the diffusivity by the Mean Squared Displacement (MSD)

[TABLE]

In the spherical geometry the standard Euclidean distance is distorted by the spherical surface. We use as our norm $\|\cdot\|_{g}$ the geodesic distance on the surface between the starting point $X(0)$ and the final point $X(t)$. Unlike bulk three dimensional fluids, the distance between points remains bounded since the surface is a compact manifold. As a consequence, we have that eventually the $MSD(t)\rightarrow m_{0}^{2}$ asymptotes to a limiting value. This corresponds to sampling the $X(t)$ from the stationary distribution over the surface.

To estimate in practice the diffusivity, we consider the $MSD$ over time-scales $\tau$ with $\tau\leq 0.06\tau_{D}$ where $\tau_{D}$ is the time-scale to diffuse the distance $R$. We find the $MSD$ is approximately linear in this regime and we characterize the passive particle motions by the diffusivity $D=\partial{MSD(t)}/{\partial{t}}$, see Figure 12. In practice, we estimate $D$ using the slope of the $MSD$ obtained from a least-squares fit to the data in the range $0\leq t\leq 0.06\tau_{D}$.

We study how the diffusivity $D$ of $N_{T}=80$ passive tracer particles are enhanced by the action of $N_{S}$ swimmers. Both the passive tracers and swimmers undergo the drift-diffusive dynamics of equation 39 with parameters in Table 2. We remark there are some technical challenges in parameterizing consistently for comparisons the non-hydrodynamic and hydrodynamic dynamics. We choose to do this by using the diagonal entry of our mobility tensor $M_{\parallel\parallel}$ as the inverse drag $\gamma_{drag}^{-1}$ for the non-hydrodynamic dynamics. Given symmetries that can arise readily for small numbers of swimmers, we consider cases with $N_{S}\in[5,10]$.

Our initial studies reported here find the swimmers can result in significant enhancement of the tracer diffusivity relative to the non-hydrodynamic case, see Figure 12. In the case with no swimmers $N_{S}=0$, we have diffusivity with hydrodynamic coupling $D_{hy}=1.70295nm^{2}/ps$ and without hydrodynamic coupling $D_{nh}=1.3916nm^{2}/ps$. This has the ratio $D_{hy}/D_{nh}=1.2237$. We see already in the absence of any swimmers an effective enhancement of $\sim 20\%$ in the diffusivity of the particles from the hydrodynamic correlations.

As we introduce more microscopic swimmers $N_{S}$, their collective stroke cycles contribute hydrodynamic flows that mix the tracer particles in addition to the thermal fluctuations. Relative to the non-hydrodynamic case, we find this manifests as an effective diffusivity that is further enhanced in the range of $\sim 40\%-100\%$, see Figure 12. The thermal fluctuations allow for the tracer particles to diffuse between streamlines of the hydrodynamic flows that are transiently generated by the swimmers. The swimmers also can have their own motions driven by the mutual flows and thermal fluctuations that serve both to move their center-of-mass and to rotate their orientation. This combination of effects shows the interplay that can arise between hydrodynamically driven drifts and thermal fluctuations.

The results illustrate some of the hydrodynamic and thermal effects that can be captured using our methods for possible further investigations of passive and active soft materials [58, 81] or mechanics in cell biology [13, 60]. The introduced surface fluctuating hydrodynamics methods can be used to model passive and active spatially extended microstructures to capture both the hydrodynamic coupling and the correlated thermal fluctuations within curved fluid interfaces.

5 Conclusions

We have introduced surface fluctuating hydrodynamics approaches for general investigations of the drift-diffusion dynamics of particles and microstructures immersed within curved fluid interfaces. We introduced computational methods for simulations of fluid-structure interactions and collective dynamics driven by active forces and thermal fluctuations. We studied the velocity autocorrelations for surface fluctuating hydrodynamics in spherical geometries. We found for interfacial fluids there are different scalings that emerge for different physical regimes and depending on the interface geometry, surface viscosity, and bulk viscosities. We also showed how our methods can be used for modeling and simulating the collective drift-diffusion dynamics of both passive and active microstructures. We obtained results investigating the enhanced mixing of particles from active microswimmers. The results show how the introduced surface fluctuating hydrodynamics approaches can be used for investigating some of the rich phenomena that can arise in curved fluid interfaces.

6 Acknowledgments

The authors P.J.A, M.P. and D.A.R. acknowledge support from research grants DOE Grant ASCR PHILMS DE-SC0019246 and NSF Grant DMS - 1616353. We also acknowledge UCSB Center for Scientific Computing NSF MRSEC (DMR-1121053) and UCSB MRL NSF CNS-0960316. P.J.A. would also like to acknowledge a hardware grant from Nvidia.

Appendix A Derivations for Fluctuating Hydrodynamics on Surfaces using Vector

Potentials $\Phi$

We discuss in more detail for fluctuating hydrodynamics on the surface derivation of the stochastic dynamics of the vector potential $\Phi$. We express the dynamics in terms of complex-valued coefficients $a_{s}$ of a spherical harmonics expansion as

[TABLE]

The force expansion terms $c_{s}^{n}$ are obtained from the real-space force $f^{n}$ using equation 2.7. We obtain the stochastic driving fields $g_{s}$ and $h_{s}$ using a fluctuation-dissipation approach for the discretized system in a manner similar to our prior work [8, 9]. We generate Gaussian driving terms with mean zero and covariance

[TABLE]

This requires determination of the covariance $\mathcal{C}_{ss^{\prime}}$ for the equilibrium fluctuations of the modes $a_{s}$.

The Gibbs-Boltzmann distribution $\rho(\mathbf{v}^{\flat})=(1/Z)\exp\left(-E[\mathbf{v}^{\flat}]/k_{B}{T}\right)$ for the equilibrium fluctuations depends on the kinetic energy of the fluid given by

[TABLE]

The $s=(\ell,m)$ where $\ell$ is the degree and $m$ the order. We used here the adjoint property of the co-differential $\boldsymbol{\delta}$ and the exterior derivative $\mathbf{d}$ [1]. The spherical harmonics modes $Y_{s}$ have $L^{2}$-norm given by $\|Y_{s}\|_{2}^{2}$ and are eigenfunctions of the Laplace-Beltrami operator $-\boldsymbol{\delta}\mathbf{d}Y_{s}=\lambda_{s}Y_{s}$ with eigenvalues $\lambda_{s}=-{\ell(\ell+1)}/{R^{2}}$.

The quadratic form of the energy yields that the equilibrium fluctuations of the spherical harmonics coefficients are Gaussian with mean zero and covariance

[TABLE]

The $\delta_{s^{\prime},\overline{s}}$ denotes the Kronecker $\delta$-function where we use notation $\overline{s}=(\ell,-m)$ to denote the conjugate mode index. The spherical harmonics coefficients are complex-valued and the field $\Phi$ must be real-valued. This requires for the coefficients $\overline{a_{s}}=a_{\overline{s}}$.

We generate the stochastic driving terms using $g_{s}=\eta_{s}+i\xi_{s}$ with $i=\sqrt{-1}$ and $\overline{g_{s}}=g_{\overline{s}}$ throughout. From the conditions in equation 74, we have for $m\neq 0$ that

[TABLE]

As a consequence, we must have that $\left\langle\eta_{s}^{2}\right\rangle=\left\langle\xi_{s}^{2}\right\rangle$ and $\left\langle\eta_{s}\xi_{s}\right\rangle=0$. This requires that $\left\langle\eta_{s}^{2}\right\rangle=\left\langle\xi_{s}^{2}\right\rangle=-L_{s}\mathcal{C}_{s,\overline{s}}$. The case with $m=0$ is special since the mode is self-conjugate requiring

[TABLE]

As a consequence, we must have that $\left\langle\xi_{s}^{2}\right\rangle=0$ and $\left\langle\eta_{s}^{2}\right\rangle=-2L_{s}\mathcal{C}_{s,\overline{s}}$.

Algorithmically, for modes $s=(\ell,m)$ these results correspond to generating $g_{s}$ for $m>0$ by computing each of the components $\eta_{s}$ and $\xi_{s}$ as independent Gaussian random variates each having the covariance $-L_{s}\mathcal{C}_{s\overline{s}}$, and for $m<0$ setting $g_{\overline{s}}=\overline{g_{s}}$. The modes with $m=0$ are special since they are self-conjugate and we have $g_{\overline{s}}=\overline{g_{s}}=\eta_{s}$ with covariance $-2L_{s}\mathcal{C}_{s\overline{s}}$.

Appendix B Derivation of Power-Laws for Fluctuating Hydrodynamics on a Sphere

We find the autocorrelation of the fluid velocity on the spherical interface has significantly different behaviors than bulk fluids. For bulk Newtonian fluids occupying a three dimensional volume the velocity autocorrelation function has a well-known characteristic long-tail with scaling $\tau^{-3/2}$. This is supported by continuum theory [7, 12], molecular simulations [3, 25], and experimental evidence [68, 28]. In contrast, we find from equation 2.2 and equation 49 the interfacial fluid velocity exhibits a few different intermediate power-law scalings and exponential decay depending on the considered parameter and temporal regimes, see Figure 3. As we shall show, this arises both from the spherical geometry and from the coupling between the two-dimensional hydrodynamics and bulk surrounding three-dimensional fluid which introduce additional time-scales.

The fluid velocity is obtained from $\Phi$ as $\mathbf{v}=(-\star\mathbf{d}\Phi)^{\sharp}=v^{\theta}\partial_{\theta}+v^{\phi}\partial_{\phi}$ where $v^{\theta}=\sum_{s}\frac{a_{s}}{\sqrt{g}R}\frac{\partial Y_{s}}{\partial\phi}$ and $v^{\phi}=-\sum_{s}\frac{a_{s}}{\sqrt{g}R}\frac{\partial Y_{s}}{\partial\theta}$. The velocity autocorrelation function associated with fluctuating hydrodynamics on the sphere given in equation 2.2 can be expressed as

[TABLE]

For the spherical surface, we use that the metric $|g|=R^{2}$. We also use that $\langle a_{s}^{2}\rangle=k_{B}TR^{2}/\rho\ell(\ell+1)$ and $\langle a_{s}(0)a_{s}(t)\rangle=\langle a_{s}^{2}\rangle\exp(tL_{s})$.

We consider the autocorrelation of the fluid at a point

[TABLE]

Given the $\delta$-spatial correlation of the fluctuating velocity field each of these series diverges in the limit $t\rightarrow 0$ when the full infinite expansion is taken. In practice in techniques such as the Stochastic Immersed Boundary Methods [8] the fluctuating velocity field is spatially averaged to model the dynamics of immersed particles and microstructures. This would result in fluid-structure interactions over the surface only having effectively a responses to fluid fluctuations above some critical length-scale (below some degree $\ell_{*}$) which is related to the object’s geometric size. We can obtain in practice a similar effect by working throughout with truncated series expansions with $\ell\leq\ell_{*}=50$ [7].

We always consider fluctuations at a point $\mathbf{x}_{*}$ on the equator for a given spherical coordinate chart which by symmetry yields

[TABLE]

and

[TABLE]

We can express equation LABEL:equ:expand_autocor as

[TABLE]

where $L_{\ell}$ is given by equation 48.

We approximate these sums asymptotically to estimate significant time-scales governing different behaviors of the autocorrelation functions. We make the ansantz throughout that we can treat the term $\left|Y_{s}\right|^{2}\sim C$. This gives

[TABLE]

We use that $\sum_{|m|\leq\ell}m^{2}=\ell(\ell+1)(2\ell+1)/3$. This very conveniently cancels the term $\ell(\ell+1)$ in the denominator that arose from the eigenvalues of the Laplace-Beltrami operator discussed in Section 2.7. After some rearrangement, we have

[TABLE]

We have absorbed also the additional prefactor constants into $C$. The spherical motion corresponding to rigid-body rotation in the bulk fluid corresponds to the mode $\ell=1$. We see for $L_{\ell}$ only the second term (traction stress term) persists in equation 48. This gives the decay time-scale $\tau_{r}={2RL}/{3\mu_{m}}$ for energy dissipation through the rigid rotation of the entire spherical interface within the bulk fluid. There are two particularly interesting regimes. The first corresponds to when $R/L\gg 1$ indicating localized hydrodynamics on the surface. The second to when $R/L\ll 1$ indicating the hydrodynamics is strongly coupled nearly over the entire surface to give a response that is effectively a rigid-body rotation of the sphere.

We consider first the case when $R/L\gg 1$ and make the approximation

[TABLE]

The $\epsilon$ term includes the higher-order terms. We take $\ell\ll R/L$ so that the $\ell(\ell+1)$ term does not play a significant role. We approximate the sum using integration to obtain

[TABLE]

In this notation, we set $\alpha={\mu_{m}}/{RL}$. We also make the assumption that $t\ll\tau_{r}={2RL}/{3\mu_{m}}$ so we can treat $\exp\left({-\frac{3}{2}\alpha t}\right)\approx 1$. In the case with $R/L\gg 1$ and $t\ll\tau_{r}={2RL}/{3\mu_{m}}$, we have that

[TABLE]

In the case approximating with $R/L\gg 1$ followed by approximating with $\alpha=\mu_{m}/{RL}\gg 1$, we have that

[TABLE]

For the autocorrelation function, these results predict two distinct regimes that will exhibit different power law scalings. The first has power law $t^{-2}$ and the second with power law scaling $t^{-1}$. The $t^{-1}$ power law is consistent with prior studies of pure two dimensional fluid interfaces predicting similar results for the long-time tail and divergence of the diffusion coefficient [12, 3]. Integrating the velocity autocorrelation functions with appropriate truncations given particle size one could obtain effective particle diffusivities within the interface using the Green-Kubo relations [48, 34, 77].

Given the coupling to the bulk surrounding fluid our fluctuating hydrodynamics have additional time-scales mediating these effects. It is interesting that even though some of our parameter regimes exhibit a $t^{-1}$ decay this in fact only persists for a finite amount of time and is eventually mitigated by our coupling to the bulk solvent fluid. Given that finite duration, our exhibited $t^{-1}$ decay would result in finite logarithmic terms in the diffusivity, which is consistent with Saffman-Delbrück theory which considers a similar regime [79]. Our fluctuating hydrodynamics not only capture the classical Saffman-Delbrück results but also extend this to include geometric contributions from the spherical shape and other additional time-scales in regimes where the interfacial hydrodynamics and coupling to the bulk fluid could differ significantly. Further extensions could also be made to include the temporal dynamics of the bulk solvent fluid.

We next consider the regime with $R/L\ll 1$ and make the approximation

[TABLE]

The $\epsilon$ term includes the higher-order terms. We see a key term is $\beta=\mu_{m}/R^{2}$. In this regime we have

[TABLE]

We obtained these results using completion of the square in the exponent.

In the regime with $R/L\ll 1$, we see the velocity autocorrelation has

[TABLE]

This predicts a power law decay with scaling $t^{-1}$.

We see that the relaxation time-scale $\tau_{r}$ for some systems can be quite large relative to the other time-scales. In these regimes, we find something interesting can occur where the velocity autocorrelation function plateaus. This is predicted by our theory to occur in the regime when the time $t$ satisfies $\tau_{a}={R^{2}}/{4\mu_{m}}\ll t\ll\tau_{r}={2RL}/{3\mu_{m}}$. In this regime, the correlations associated with the internal flow of the hydrodynamics within the interfacial decays rapidly to zero. However, the rigid rotational motion of the entire fluid interface can still persist for awhile until the rotation reaches its decay time-scale that dissipates this motion. This leads to the interesting plateaus seen in the velocity autocorrelation function in Figure 3. We also see that for all non-zero parameter choices the autocorrelation function will eventually exhibit an exponential decay when reaching time-scale $t\gg\tau_{r}={2RL}/{3\mu_{m}}$. This is a consequence of the rigid-body rotational mode $\ell=1$ being the longest lived mode and eventually dissipating energy from the interfacial fluid to the bulk surrounding fluid. Our calculations show that surface fluctuating hydrodynamics on quasi two dimensional fluid interfaces can exhibit significantly different phenomena relative to their bulk counter-parts in three dimensional space.

Appendix C Spherical Harmonics

We expand functions $\Phi$ on the surface using the spherical harmonics

[TABLE]

where

[TABLE]

The $m$ denotes the order and $n$ the degree for $n\geq 0$ and $m\in\{-n,\dots,n\}$. The $P^{m}_{n}$ denote the Associated Legendre Polynomials. We denote by $\theta$ the azimuthal angle and by $\phi$ the polar angle of the spherical coordinates [6]. We work with real-valued functions and use that modes are self-conjugate in the sense $Y^{m}_{n}=\overline{Y^{-m}_{n}}$.

We can express the spherical harmonic modes as

[TABLE]

The $X_{n}^{m}$ and $Z_{n}^{m}$ denote the real and imaginary parts. We use this splitting in our numerical methods to construct a purely real set of basis functions on the unit sphere with maximum degree $N$ which consists of $(N+1)^{2}$ basis elements. For the case $N=2$ we have the basis elements

[TABLE]

Similar conventions are used for the basis for the other values of $N$. We take final basis elements $Y_{i}$ that are normalized as $Y_{i}=\tilde{Y}_{i}/\sqrt{\langle\tilde{Y}_{i},\tilde{Y}_{i}\rangle}$.

Derivatives are used within our finite expansions by evaluating analytic formulas whenever possible for the spherical harmonics in order to try to minimize approximation error [6]. Approximation errors are incurred when sampling the values of expressions involving these derivatives at the Lebedev nodes and when performing quadratures [52]. The derivative of the spherical harmonics in the azimuthal coordinate $\theta$ is given by

[TABLE]

We see this has the useful feature that the derivative in $\theta$ of a spherical harmonic of degree $n$ is again a spherical harmonic of degree $n$. As a consequence, we have in our numerics that this derivative can be represented in our finite basis. This allows us to avoid additional $L^{2}$ projections allowing for computation of the derivative in $\theta$ without incurring an approximation error. The derivative of the spherical harmonics in the polar angle $\phi$ is given by

[TABLE]

We see that unlike derivatives in $\theta$ the derivative in $\phi$ can not be represented in general in terms of a finite expansion of spherical harmonics. In our numerics, we use the expression in equation C for $\partial_{\phi}Y^{m}_{n}(\theta,\phi)$ when we need to compute values at the Lebedev quadrature nodes. These analytic results provide a convenient way to compute derivatives of differential forms following the approach discussed in our prior paper [35]. By using these analytic expressions, we have that the subsequent hyperinterpolation of the resulting expressions are where the approximation errors are primarily incurred. Throughout our discussions to simplify the notation we use the convention that $Y^{m}_{n}=0$ when $m\geq n+1$. Further discussion of spherical harmonics can be found [6]. Further discussions about how we use the spherical harmonics in our numerical calculations of exterior calculus operators also can be found in our papers [82, 35].

Appendix D Exterior Calculus: Coordinate Expressions

We use approaches from exterior calculus to generalize operators used in continuum mechanics to the manifold setting. We give here more explicit expressions for these operators in terms of the local surface coordinates. Since there is no global non-singular coordinate system on the sphere, we ensure numerical accuracy by switching between two coordinate charts. In chart $A$ we have coordinates $(\hat{\theta},\hat{\phi})$ with singularities at the north and south poles. In chart $B$ we have coordinates $(\tilde{\theta},\tilde{\phi})$ having singularities at the east and west poles. To avoid issues with singularities when seeking a value at a point $\mathbf{x}$, we evaluate expressions within each chart in the regions with $\pi/4\leq\phi\leq 3\pi/4$ and $\pi/4\leq\tilde{\phi}\leq 3\pi/4$. We give all expressions with generic polar coordinates $(\theta,\phi)$ which we subsequently use in practice in our numerical calculations by choosing the appropriate chart $A$ or chart $B$. More details on our approach can also be found in [82].

The exterior derivative $\mathbf{d}$ for a [math]-form $f$ is given by $\mathbf{d}f=\partial f/\partial x_{j}\,\mathbf{x}_{j}$ and for a 1-form $\boldsymbol{\alpha}=\alpha^{i}d\mathbf{x}_{i}$ is $\mathbf{d}\boldsymbol{\alpha}=\partial\alpha^{i}/\partial x^{j}\,d\mathbf{x}_{j}\wedge d\mathbf{x}_{i}$. The $\boldsymbol{\delta}=-\star\mathbf{d}\star$ is the co-differential playing a role similar to the divergence on the surface [1]. The Hodge $\star$ for a differential $k$-form $\boldsymbol{\beta}$ gives a complementary $n-k$-form $\star\boldsymbol{\beta}$ so that for any $k$-form $\boldsymbol{\alpha}$ we have $\boldsymbol{\alpha}\wedge\star\boldsymbol{\beta}=\langle\alpha,\beta\rangle\boldsymbol{\omega}$ where $\boldsymbol{\omega}$ is the volume form [1]. This allows us to generalize vector calculus operators such as the curl and divergence to the surface by $\mbox{curl}_{\mathcal{M}}(\mathbf{v}^{\flat})=-\star\mathbf{d}\mathbf{v}^{\flat}$ and $\mbox{div}_{\mathcal{M}}(\mathbf{v}^{\flat})=\boldsymbol{\delta}\mathbf{v}^{\flat}$.

We now give coordinate expressions for these operations in terms of the metric tensor and more specialized expressions in the case of two dimensional manifolds (radial $2$-manifolds). For radial $2$-manifolds, the exterior derivatives can be expressed for a [math]-form $f$ and $1$-form $\boldsymbol{\alpha}$ as

[TABLE]

The generalized curl in this setting for [math]-form and $1$-form can be expressed as

[TABLE]

In this notation we have taken the conventions that $f_{j}=\partial_{x^{j}}f$ and $\alpha_{j}$ such that $\boldsymbol{\alpha}=\alpha_{j}\mathbf{d}{x}^{j}$ where $j\in\{\theta,\phi\}$. The isomorphisms $\sharp$ and $\flat$ between vectors and co-vectors can be expressed explicitly as

[TABLE]

We use the notational conventions here that for the embedding map $\boldsymbol{\sigma}$ for spherical coordinates in $\mathbb{R}^{3}$ we have $\boldsymbol{\sigma}_{\theta}=\partial_{\theta}$ and $\boldsymbol{\sigma}_{\phi}=\partial_{\phi}$. Combining the above equations we can express the generalized curl as

[TABLE]

The scalar Laplace-Beltrami operator $\Delta_{LB}=-\boldsymbol{\delta}\mathbf{d}$ that acts on [math]-forms can be expressed in coordinates as

[TABLE]

The $g_{ij}$ denotes the metric tensor, $g^{ij}$ the inverse metric tensor, and $|g|$ the determinant of the metric tensor.

The velocity field of the hydrodynamic flows $\mathbf{v}$ is recovered from the vector potential $\Phi$ as $\mathbf{v}^{\flat}=\mbox{curl}_{\mathcal{M}}(\Phi)=-\star\mathbf{d}\Phi$. The velocity field is obtained from $\mathbf{v}=\mathbf{v}^{\sharp}=\left(-\star\mathbf{d}\Phi\right)^{\sharp}$ using equation LABEL:equ_gen_curl_0form. Similarly, from the force density $\mathbf{b}$ acting on the fluid, we obtain the data $\mbox{curl}_{\mathcal{M}}(\mathbf{b}^{\flat})=-\star\mathbf{d}\mathbf{b}^{\flat}$ for the vector potential formulation of the hydrodynamics using equation LABEL:equ_gen_curl_1form. Additional details and discussions of these operators can be found in our related papers [82, 35] and in [71, 1, 83].