# Hydrodynamic Aggregation of Membrane Inclusions due to Non-Newtonian Surface
Rheology

Vishnu Vig  Harishankar Manikantan  Department of Chemical Engineering, UC
Davis, Davis, CA<EMAIL_ADDRESS>
###### Abstract

Biological membranes are self-assembled complex fluid interfaces that host
proteins, molecular motors and other macromolecules essential for cellular
function. These membranes have a distinct in-plane fluid response with a
surface viscosity that has been well characterized. The resulting quasi-2D
fluid dynamical problem describes the motion of embedded proteins or
particles. However, the viscous response of biological membranes is often non-
Newtonian: in particular, the surface shear viscosity of phospholipids that
comprise the membrane depends strongly on the surface pressure. We use the
Lorentz reciprocal theorem to extract the effective long-ranged hydrodynamic
interaction among membrane inclusions that arises due to such non-trivial
rheology. We show that the corrective force that emerges ties back to the
interplay between membrane flow and non-constant viscosity, which suggests a
mechanism for biologically favorable protein aggregation within membranes. We
quantify and describe the mechanism for such a large-scale concentration
instability using a mean-field model. Finally, we employ numerical simulations
to demonstrate the formation of hexatic crystals due to the effective
hydrodynamic interactions within the membrane.

††preprint: AIP/123-QED

##  I Introduction

Lipids and other surface-active macromolecules play critical roles in
biological interfaces. Membranes of eukaryotic cells are self-assembled
phospholipid bilayers [1, 2]. These complex systems display rich dynamics due
to their inherent quasi-2D viscous nature that emerges from the coupling to
the surrounding bulk phases [3, 4]. Structural and rheological properties of
fluid membranes have been well characterized in the past decades, conceptually
driven by connecting in-plane protein diffusion to membrane fluidity. Surface
viscosities are now known for single and multi-component lipid monolayers [5,
6, 7], phospholipid bilayers [8, 9, 10], protein-laden interfaces [11], and
polymeric multilayers [12]. Simultaneously, continuum fluid dynamical models
have been widely successful in describing the isolated and cooperative motion
of disk-like [13, 14, 15], rod-like [16, 17, 18], active [19, 20, 21], and
polymeric [22] inclusions in simple 2D Newtonian membranes.

However, the rheology of real membranes is rarely Newtonian. Insoluble
surfactant species often phase separate to form 2D dispersions or ‘rafts’ of
condensed liquid crystalline phases suspended in a liquid expanded disordered
phase [23, 24, 25]. Unlike 3D fluids, surfactants are much easier to compress
and even a facile in-plane deformation triggers significant changes in the
rheological signature of the lipid interface. The transport of surfactant
molecules between coexisting phases as well as the morphology and interlocking
of stiffer domains contribute to a strongly nonlinear surface rheological
response. Insoluble surfactant layers exhibit jamming and yield stress
behavior [5], surface viscoelasticity [26], surface shear thinning [27], and
surface-concentration-dependent surface rheology [28, 7]. While a lot of the
intuition of 3D fluid mechanics can be ported to the analysis of 2D viscous
manifolds, the nature of momentum transport in interfacial systems and the
novel rheological relations present unique challenges.

In this work, we aim to theoretically investigate and numerically demonstrate
the role of a class of non-Newtonian membrane behavior on driven or anchored
trans-membrane proteins or protein-associated domains that play critical roles
in cell motility and structure [29, 30]. These membrane inclusions are
anchored to and driven relative to the membrane by the intracellular network
of actin and microtubules, or by external motor proteins such as kinesin or
dynein [1]. We will specifically target the nontrivial hydrodynamic
interactions expected to arise due to surface-pressure–dependent surface
viscosity of lipids. Despite this strong nonlinearity, recent mathematical
efforts have demonstrated the qualitatively new phenomena that emerge in thin
interfacial gaps [31], in pore-spanning monolayers or membranes [32], and in
the resulting stability and dynamics of drops containing such surfactants [33,
34]. These past works point to kinematic symmetry breaking within the 2D
interfacial layer, leading to trajectories of probe particles, lipid rafts, or
naturally occuring membrane inclusions that are not expected in Newtonian
systems. Yet, the exact nature of the hydrodynamic interactions between pairs
of driven inclusions representing cytoskeletal anchors, membrane motors, and
adhesion junction proteins remain a mystery. We derive analytic results for
such hydrodynamic interactions in this work, shedding light on the effective
inter-particle attraction or repulsion between interfacially driven particles.
Collectively, such interactions can lead to large-scale aggregation, which
might be biologically favorable in immune response, locomotion, and tubulation
[30]. Past work on membrane inclusions have demonstrated long-ranged
interactions and assembly purely due to protein activity [21, 20], membrane
curvature [35, 36], and inclusion size mismatch [30]. Building on our analysis
of pair dynamics, we demonstrate here that non-Newtonian surface rheological
response can also lead to aggregation and packing of ‘2D suspensions’ of
driven or anchored inclusions on the membrane.

This paper is organized as follow: we first evaluate the motion of a single
particle driven in a rheologically complex membrane in §II. In this section,
we formulate the fluid dynamical problem of a disk-like particle embedded in
the membrane and driven by an external force while subject to the 2D flow in
the plane of the membrane. We will develop a perturbative solution for weakly
non-Newtonian behavior, derive the corrected velocity of the particle, and
physically rationalize the result. In §III, we extend the study to collective
behavior via membrane hydrodynamic interactions due to ‘2D suspensions’ of
particles driven in the plane of a membrane. We propose a mean-field
description to quantify the instability that leads to large-scale aggregation,
develop a Langevin description that allows us to simulate multi-particle
dynamics, and explore the hexatic order that emerges from crystalline
aggregation due to attractive hydrodynamic interactions. We conclude with a
discussion of the relevance of these results to real biological systems, and
connections to other problems where viscosity depends on pressure.

##  II Single particle dynamics in non-Newtonian membranes

###  II.1 Problem formulation

Figure 1: Illustration of a disk-like inclusion representing a transmembrane
protein embedded in a phospholipid membrane. The inclusion may be driven by
intracellular forces exerted by the cytoskeleton, or by extracellular forces
due to motor protein activity during transport or at adhesion junctions. The
fluid dynamical problem is equivalent to an anchored immobile inclusion within
a flowing membrane. Forces driving such inclusions can generate in-plane
hydrodynamic disturbance flows and interactions.

The geometry of our fluid dynamical system is shown in Fig. 1. A bulk fluid
phase of viscosity $\eta$ representing intracellular space underlies a
membrane of potentially non-constant surface shear viscosity $\eta_{s}$. Fluid
flow in the bulk phase with velocity $\mathbf{v}$ is governed by the Navier-
Stokes equation, with a no-slip condition to match the membrane velocity
$\mathbf{u}$. Additionally, the stress jump at the interface describes
conservation of 2D momentum [37, 38]:

$\rho_{s}\frac{D\mathbf{u}_{s}}{Dt}=\bm{\nabla}_{s}\bm{\cdot}\bm{\sigma}_{s}-\eta\partialderivative{\mathbf{v}}{z}\bigg{|}_{\text{z}=0},$
(1)

where we have assumed a flat membrane along the $x$–$y$ plane for simplicity,
and $D/Dt$ is the material derivative.
$\bm{\nabla}_{s}=\mathbf{I}_{s}\bm{\cdot}\bm{\nabla}$ is the surface gradient
operator with $\mathbf{I}_{s}=\mathbf{I}-\hat{\mathbf{n}}\hat{\mathbf{n}}$ the
surface identity tensor on a plane with unit normal $\hat{\mathbf{n}}$.
Equation (1) is readily generalized to bulk fluids on either side of the
membrane: the results that follow are qualitatively unchanged and are only
modified by a prefactor in that case.

Equation (1) is a 2D Cauchy equation for a species of density $\rho_{s}$ and
2D stress tensor $\bm{\sigma}_{s}$ that is forced via viscous tractions from
the adjacent bulk phase(s). Alternatively, we may interpret Eq. (1) as a
boundary condition for the Navier-Stokes equation governing 3D fluid field
$\mathbf{v}$. The 2D stress tensor $\bm{\sigma}_{s}$ may be decomposed into
isotropic and deviatoric parts:

$\bm{\sigma}_{s}=-\Pi\,\mathbf{I}_{s}+\bm{\tau}_{s},$ (2)

where $\Pi$ is the 2D pressure and the deviatoric stress tensor
$\bm{\tau}_{s}$ is prescribed by a constitutive relation. For instance,

$\bm{\tau}_{s}=\eta_{s}\left[\bm{\nabla}_{s}\mathbf{u}\bm{\cdot}\mathbf{I}_{s}+\mathbf{I}_{s}\bm{\cdot}(\bm{\nabla}_{s}\mathbf{u})^{\text{T}}\right]$
(3)

in the simple case of a 2D Newtonian membrane. Underlying such a constitutive
relation is the approximation that the interface is 2D incompressible. In
other words, we assume that the surface flow rearranges to that of a 2D
incompressible material so long as Marangoni flows (driven by gradients in the
surface pressure $\Pi$) are generated faster than the rate of in-plane
compression. This is almost always the case for insoluble phospholipids that
make up most biological membranes [38, 39, 40]. Then, neglecting inertia in
these overdamped systems, Eq. (1) simplifies to

$\eta\partialderivative{\mathbf{v}}{z}\bigg{|}_{\text{z}=0}=-\bm{\nabla}_{s}\Pi+\eta_{s}\nabla^{2}\mathbf{u}_{s},\quad\bm{\nabla}_{s}\bm{\cdot}\mathbf{u}_{s}=0.$
(4)

Equation (4) governs momentum and mass conservation of an insoluble
incompressible Newtonian membrane with constant surface shear viscosity. A
wide range of theoretical studies spanning decades [3, 4, 13, 14] describing
in-plane motion of particles in Newtonian membranes start with Eq. (4).

However, we intend to explore the role of non-constant $\eta_{s}$ on membrane
inclusion dynamics. Specifically, we wish to account for surface-pressure-
dependent surface viscosity of lipids. In what follows, we assume that surface
pressure enforces incompressibility and $\Pi$-dependent viscosities are
incorporated by treating the the interfacial momentum equation (4) as a
generalized Newtonian model. This is analogous to dropping the equation of
state for density in the incompressible Navier-Stokes equation: the
thermodynamic quantity $\Pi$ plays the role of the ‘mechanical’ pressure
within this approximation. The validity of such an approximation is diskussed
in detail in past works [31] and this model has been widely successful in
studying these complex systems [32, 33, 34].

Such a description of the non-constant surface shear viscosity also relates
directly to typical experiments where the surface viscosity is measured as a
function of surface pressure. $\eta_{s}$ changes order of magnitude over a
range of $\Pi$ commonly accessible in experiments, unlike 3D fluids whose
viscosity only changes under extreme pressures, if at all. For example, the
phospholipid dipalmitoylphosphatidylcholine (DPPC) is a major constituent of
cell membranes and forms stable monolayers at an air-water interface that
undergo a liquid expanded (LE) to liquid condensed (LC) phase transition at
$\Pi\sim 8\,\text{mN/ m}$ at room temperature [23]. Above this critical
surface-pressure, the monolayer viscosity increases exponentially with surface
pressure [28]. Similar ‘$\Pi$-thickening’ behavior has been observed in
certain fatty acids like nonadaconic (C19) acid, heneicosanoic (C21) acid and
behenic (C22) acid [41]. Conversely, phase transitions associated with chain-
tilt lead to molecules shearing past each other: as a result $\eta_{s}$
decreases exponentially with increasing $\Pi$ from $10$ to $20\,\text{mN/m}$
in the ‘$\Pi$-thinning’ surfactant eicosanol [42].

The simplest constitutive relation for $\eta_{s}(\Pi)$ follows free-area
models of surface viscosity [41, 43] to give an exponential relation:

$\eta_{s}(\Pi)=\eta^{0}_{s}e^{\Pi-\Pi_{\infty}/\Pi_{c}},$ (5)

where $\Pi_{c}$ is the characteristic surface pressure change required to
produce a noticeable change in $\eta_{s}$, and $\eta^{0}_{s}$ is a reference
or unperturbed viscosity at a reference pressure $\Pi_{\infty}$. Setting
$\Pi_{c}\rightarrow\infty$ retrieves the Newtonian limit of constant surface
viscosity. Equation (5) accommodates both $\Pi$-thickening ($\Pi_{c}>0$) and
$\Pi$-thinning ($\Pi_{c}<0$) surfactants. Plugging Eq. (5) into the Cauchy
momentum equation and simplifying in the overdamped limit then gives

$-{\bm{\nabla}}_{s}\Pi+{\bm{\nabla}}_{s}\bm{\cdot}\left[\eta_{s}(\Pi)\left({\bm{\nabla}}_{s}\mathbf{u}_{s}+{\bm{\nabla}}_{s}\mathbf{u}_{s}^{T}\right)\right]=\eta\partialderivative{\mathbf{v}}{z}\bigg{|}_{\text{z}=0},$
(6)

which along with the incompressibility condition
$\bm{\nabla}_{s}\bm{\cdot}\mathbf{u}_{s}=0$ governs our system.

###  II.2 Non-dimensionalization and perturbation expansion

Scaling velocities and lengths over characteristic values $U$ and $a$ gives a
characteristic surface pressure $\Pi_{0}=\eta_{s}^{0}U/a$, and the
dimensionless surface momentum equation becomes:

$-\tilde{\bm{\nabla}}_{s}\tilde{\Pi}+\tilde{\bm{\nabla}}_{s}\bm{\cdot}\left[\tilde{\eta}_{s}(\tilde{\Pi})\left(\tilde{\bm{\nabla}}_{s}\mathbf{\tilde{u}}_{s}+\tilde{\bm{\nabla}}_{s}\mathbf{\tilde{u}}_{s}^{T}\right)\right]=\frac{1}{Bq}\partialderivative{\mathbf{\tilde{v}}}{\tilde{z}}\bigg{|}_{\text{z}=0}.$
(7)

In problems where a characteristic external force $F$ (instead of a velocity
$U$) is specified, we can equivalently define $U=F/\eta_{s}^{0}$ or
$\Pi_{0}=F/a$. The dimensionless surface viscosity is
$\tilde{\eta}_{s}=\eta_{s}/\eta_{s}^{0}=e^{\Pi/\Pi_{c}}$ and the Boussinesq
number

$Bq=\frac{\eta_{s}}{\eta a}$ (8)

compares surface viscous stresses to subphase viscous stresses. Equation (7)
is essentially a 2D Stokes equation with non-constant viscosity forced by the
traction from the bulk phase. Indeed, the solution of the corresponding
Newtonian problem at large $Bq$ limits to 2D Stokes flow. Such a solution is
classically known to be singular at large distances (due to the ‘Stokes
paradox’ [44, 38]). In the case of membranes or other viscous interfaces,
traction from the bulk ultimately catches up to interfacial viscous stresses
over the Saffman-Delbrück length [3, 38] and regularizes the singularity. To
keep the analysis tractable, we will work under the large-$Bq$ limit, and
assume that the Saffman-Delbrück length is sufficiently large so that we can
safely ignore the forcing term on the RHS of Eq. (7). We will appeal later to
the physically realistic situation of a finite-size membrane to regularize the
2D Stokes description as described by Saffman & Delbrück [3].

In what follows, we drop the hats on dimensionless variables. The
dimensionless governing equations are

$\bm{\nabla}_{s}\bm{\cdot}{\mathbf{\bm{\sigma}}}=-\bm{\nabla}_{s}\Pi+\bm{\nabla}_{s}\bm{\cdot}\left[\eta_{s}(\Pi)\mathbf{S}\right]=0,\quad\mathbf{\bm{\nabla}_{s}\cdot{u}}=0,$
(9)

where
$\mathbf{S}=\bm{\nabla}_{s}\mathbf{u}_{s}+\bm{\nabla}_{s}\mathbf{u}_{s}^{T}$
is twice the typical notation of the rate-of-strain tensor. Equation (9) is
heavily nonlinear due to the dependence on $\eta(\Pi)$. To make analytical
progress, we take a perturbative approach for large $\Pi_{c}$ or small
departures of $\eta_{s}(\Pi)$ from $\eta^{0}_{s}$:

$\eta_{s}(\Pi)=1+\frac{\Pi}{\;\Pi_{c}}+\mathcal{O}(\Pi_{c}^{-2})=1+\beta\frac{\Pi}{\;\Pi_{0}}+\mathcal{O}(\beta^{2}),$
(10)

where

$\beta=\frac{\Pi_{0}}{\Pi_{c}}=\frac{\eta_{s}^{0}U}{\Pi_{c}a}=\frac{F}{\Pi_{c}a},$
(11)

is a dimensionless parameter that is small for weak pressure dependence of
viscosity. Note that $\beta$ is positive when $\Pi$-thickening and negative
when $\Pi$-thinning.

The momentum equation then expands to

$\bm{\nabla}_{s}\bm{\cdot}{\mathbf{\bm{\sigma}}}=-\bm{\nabla}_{s}\Pi+\bm{\nabla}_{s}\bm{\cdot}\mathbf{S}+\beta\bm{\nabla}_{s}\bm{\cdot}\left[\Pi\mathbf{S}\right]+\mathcal{O}(\beta^{2}).$
(12)

In the weakly non-linear regime when $\beta\ll 1$, we asymptotically treat the
interfacial velocity and surface pressure fields as regular expansions in
$\beta$:

$\displaystyle\mathbf{u}$
$\displaystyle=\mathbf{u}^{(0)}+\beta\mathbf{u}^{(1)}+\mathcal{O}(\beta^{2}),$
(13a) $\displaystyle\Pi$
$\displaystyle=\Pi^{(0)}+\beta\Pi^{(1)}+\mathcal{O}(\beta^{2}),$ (13b)

so that the corresponding dimensionless stress tensor can be written as:

$\displaystyle\bm{\sigma}\hskip 10.0pt$
$\displaystyle=\bm{\sigma}^{(0)}+\beta\bm{\sigma}^{(1)}+\mathcal{O}(\beta^{2}),$
(14a) $\displaystyle\bm{\sigma}^{(0)}$
$\displaystyle=-\Pi^{(0)}\mathbf{I+S}^{(0)},$ (14b)
$\displaystyle\bm{\sigma}^{(1)}$
$\displaystyle=-\Pi^{(1)}\mathbf{I+S}^{(1)}+\Pi^{(0)}\mathbf{S}^{(0)}.$ (14c)

We note that the non-Newtonian term in the $\mathcal{O}(\beta)$ stress tensor
$\bm{\sigma}^{(1)}$ depends on the leading-order solution. The momentum and
mass conservation equations become:

$\displaystyle\bm{\nabla}_{s}\bm{\cdot}{\bm{\sigma}^{(0)}}$
$\displaystyle=\bm{\nabla}_{s}\bm{\cdot}\mathbf{u}^{(0)}=0$ (15a)
$\displaystyle\bm{\nabla}_{s}\bm{\cdot}{\bm{\sigma}^{(1)}}$
$\displaystyle=\bm{\nabla}_{s}\bm{\cdot}\mathbf{u}^{(1)}=0$ (15b)

We will first solve the leading-order Newtonian problem for a disk placed in
arbitrary background surface flow. Then, we use the resulting stress field as
the heterogeneity in the governing equation in $\mathcal{O}(\beta)$ to obtain
the correction at $\mathcal{O}(\beta)$.

###  II.3 Isolated disk driven in a background flow

In building towards a description of the coupled dynamics of multiple disks,
we first examine the Newtonian response of a single disk driven by an external
force while placed in a background flow. We hope to then port the intuition as
well as mathematical results to the case of a dilute 2D suspension of
interacting disks by summing over the pair-wise hydrodynamic interactions.
Such a Newtonian solution describes the leading order flow corresponding to
${\bm{\sigma}^{(0)}}$, $\mathbf{u}^{(0)}$, and $\Pi^{(0)}$.

We analyze the translation of an isolated particle driven by a force
$\mathbf{F}_{p}$ in the non-Newtonian interface with an imposed linear
velocity field given by:

$\mathbf{V(x)=V}_{0}+\mathbf{x}\bm{\cdot}\bm{\Gamma},$ (16)

where $\bm{\Gamma}=\bm{\nabla}_{s}\mathbf{V}$ is the imposed background
velocity gradient. Let $q(x)$ and
$\mathbf{A}=\bm{\nabla}_{s}\mathbf{V}+(\bm{\nabla}_{s}\mathbf{V})^{T}$ be the
surface pressure and surface rate-of-strain fields associated with the imposed
background surface flow. Meanwhile, we set $\Pi$ and $\mathbf{S}$ to be the
pressure and rate of strain due to forced translation of the disk due to an
applied external force $\mathbf{F}_{p}$.

The stress fields at the leading and perturbed order, following Eq. (14), then
take the forms

$\displaystyle\bm{\sigma}^{(0)}=-(\Pi^{(0)}+q)\mathbf{I}_{s}+(\mathbf{S}^{(0)}+\mathbf{A}),$
(17)
$\displaystyle\bm{\sigma}^{(1)}=-\Pi^{(1)}\mathbf{I}_{s}+\mathbf{S}^{(1)}+(\Pi^{(0)}+q)(\mathbf{S}^{(0)}+\mathbf{A}).$
(18)

It is advantageous in what follows to have disturbance fields decay to zero
far from the particle. So, we define disturbance flow variables
$\mathbf{\hat{u}}$, $\hat{\Pi}$, and $\hat{\bm{\sigma}}$ by subtracting
imposed background flow variables from the total flow variables:

$\displaystyle\mathbf{\hat{u}(x)=u(x)-v(x)},$ (19a)
$\displaystyle\hat{\Pi}(\mathbf{x})=\Pi_{t}-q(\mathbf{x}),$ (19b)
$\displaystyle\hat{\bm{\sigma}}(\mathbf{x})=\bm{\sigma}(\mathbf{x})-\bm{\tau}(\mathbf{x}),$
(19c)

where $\bm{\tau}=-q\mathbf{I}_{s}+\mathbf{A}+\beta q\mathbf{A}$. These
disturbance variables will also satisfy the momentum conservation and
continuity equations and can be asymptotically expanded as a regular expansion
in $\beta$:

$\displaystyle\mathbf{\hat{u}}$
$\displaystyle=\mathbf{\hat{u}}^{(0)}+\beta\mathbf{\hat{u}}^{(1)}+\ldots$
(20a) $\displaystyle\hat{\Pi}$
$\displaystyle=\hat{\Pi}^{(0)}+\beta\hat{\Pi}^{(1)}+\ldots$ (20b)
$\displaystyle\hat{\bm{\sigma}}$
$\displaystyle=\hat{\bm{\sigma}}^{(0)}+\beta\hat{\bm{\sigma}}^{(1)}+\ldots$
(20c)

The $\mathcal{O}(1)$ disturbance problem corresponds to the linear Newtonian
solution and can be written as the sum of disturbance velocity due to a
translating disk and disturbance velocity due to the presence of a disk in an
ambient flow. Using standard singularity methods of 2D Stokes flow [45, 46],
this velocity field can be shown to be

$\begin{split}\mathbf{\hat{u}}^{(0)}(\mathbf{x})=\left[2\left(-\ln{(r)}\mathbf{I}+\frac{\mathbf{xx}}{r^{2}}\right)+\left(\frac{\mathbf{I}}{r^{2}}-\frac{2\mathbf{xx}}{r^{4}}\right)\right]\bm{\cdot}\frac{\mathbf{F}_{p}}{4\pi}\\\
-\left[\frac{\mathbf{xxx}}{r^{4}}+\frac{1}{2}\left(\frac{\mathbf{Ix}}{r^{4}}-\frac{2\mathbf{xxx}}{r^{6}}\right)\right]\bm{:}\mathbf{A}.\end{split}$
(21)

The terms proportional to $\mathbf{F}_{p}$ represent flow due to forced
translation, whereas those proportional to the imposed velocity gradient via
$\mathbf{A}$ are disturbances to background flow due to the presence of the
disk. Note that the leading term due to the external force, corresponding to
the 2D stokeslet, has a logarithmic singularity that is tied to the Stokes
paradox mentioned previously. In reality, flow in the plane of the membrane is
governed by 3D hydrodynamics ($\mathbf{\hat{u}}^{(0)}\propto 1/r$) beyond the
Saffman-Delbrück length, decaying to zero as the viscous traction from the
subphase becomes dominant in the momentum balance. We will safely assume that
the membrane is finite and within the momentum crossover or Saffman-Delbrück
length, and this justifies the use of Eq. (21) in capturing membrane
hydrodynamics when $Bq\gg 1$.

The leading-order pressure disturbances also enter the $\mathcal{O}(\beta)$
governing equation. Like with the velocity, the Newtonian response can be
written as a sum of contributions due to external force and imposed background
velocity:

$\hat{\Pi}^{(0)}(\mathbf{x})=\left[\frac{4\mathbf{x}}{r^{2}}\right]\bm{\cdot}\frac{\mathbf{F}_{p}}{4\pi}-\left[\frac{2\mathbf{xx}}{r^{4}}\right]\bm{:}\mathbf{A}.$
(22)

Force and torque balance ($\mathbf{F}_{\text{net}}=\mathbf{F}_{p}$ and
$\mathbf{T}_{\text{net}}=\mathbf{0}$) with these Newtonian fields gives the
translational velocity $\mathbf{U}_{p}=\mathbf{V}_{0}+\mathbf{F}_{p}/4\pi$ and
rotational velocity $\bm{\Omega}_{p}=\bm{\nabla}_{s}\times\mathbf{V}$ at
leading order.

###  II.4 Reciprocal theorem: non-Newtonian problem

With the leading-order solution now known, we can write the non-Newtonian
stress tensor at $\mathcal{O}(\beta)$ as

$\hat{\bm{\sigma}}^{(1)}=-\Pi^{(1)}\mathbf{I}+\mathbf{S}^{(1)}+(\Pi^{(0)}+q)(\mathbf{S}^{(0)}+\mathbf{A})-(q\mathbf{A}),$
(23)

The effect of pressure-dependent surface viscosity could yield a non-zero
$\mathbf{\hat{F}}\textsuperscript{(1)}$ at $\mathcal{O}(\beta)$, given by:

$\mathbf{\hat{F}\textsuperscript{(1)}}=\int\hat{\mathbf{n}}\bm{\cdot}\hat{\bm{\sigma}}^{(1)}dl$
(24)

where, $\mathbf{\hat{n}}$ is the normal vector pointing from the disk into the
2D fluid and $l$ denotes the boundary of the domain. As we have set up
velocity field to decay far from the disk, this boundary is simply the
perimeter of the disk. Without loss of generality, we choose to constrain the
particle to its leading-order motion and determine the force or torque
generated at $\mathcal{O}(\beta)$, i.e.,
$\mathbf{U}_{p}^{(1)}=\bm{\Omega}_{p}^{(1)}=\mathbf{0}$. The governing
equations at this order, Eq. (15b), thus satisfy the boundary conditions

$\displaystyle\mathbf{\hat{u}}^{(1)}$
$\displaystyle=\mathbf{0}\quad\text{for}\quad\absolutevalue{\mathbf{x}}\leq
1,$ (25) $\displaystyle\mathbf{\hat{u}}^{(1)}$
$\displaystyle\rightarrow\mathbf{0}\quad\text{as}\quad\absolutevalue{\mathbf{x}}\rightarrow\infty.$
(26)

Evaluating the non-linear effects on the disk involves solving the
$\mathcal{O}(\beta)$ problem for $\mathbf{\hat{u}}^{(1)}$. This is an
inhomogeneous Stokes equation that is extremely difficult to solve
analytically. Instead, we make use of the Lorentz reciprocal theorem [47] to
get around solving for the full flow and pressure fields using standard
methods. This framework is well-established, and has previously been employed
to evaluate perturbative solutions to several classes of weakly non-Newtonian
problems involving viscoelasticity [48, 32, 49].

To set up the reciprocal theorem, we first need an auxiliary solution
$\\{\Pi\textsubscript{aux},\;\mathbf{u\textsubscript{aux}}\;\bm{\sigma}\textsubscript{aux}\\}$
that satisfies the homogeneous Stokes problem. We choose the solution to the
translation of a disk in a 2D Newtonian interface with a velocity
$\mathbf{{U}\textsubscript{aux}}$. This solution is standard and is readily
given by the singularity solutions corresponding to steady translation (i.e.,
the parts containing $\mathbf{F}_{p}$) in Eqs. (21) and (22), with
$\mathbf{{U}\textsubscript{aux}}$ replacing $\mathbf{F}_{p}/4\pi$ in both
cases. This Newtonian problem satisfies
$\bm{\nabla}_{s}\bm{\cdot}{\bm{\sigma}\textsubscript{aux}}=\mathbf{0}$.

The reciprocal relation can be derived starting with the following vector
relations:

$\displaystyle\bm{\nabla}_{s}\bm{\cdot}\left({\bm{\sigma}\textsubscript{aux}}\bm{\cdot}\mathbf{\hat{u}}^{(1)}\right)={\bm{\sigma}\textsubscript{aux}}\bm{:}\bm{\nabla}\mathbf{\hat{u}}^{(1)},$
(27a)
$\displaystyle\bm{\nabla}_{s}\bm{\cdot}\left(\hat{\bm{\sigma}}^{(1)}\bm{\cdot}\mathbf{u}\textsubscript{aux}\right)=\hat{\bm{\sigma}}^{(1)}\bm{:}\bm{\nabla}\mathbf{u}\textsubscript{aux}.$
(27b)

Subtracting Eq. (27b) from Eq. (27a) and integrating over the entire 2D fluid
domain gives:

$\begin{split}\int\bm{\nabla}_{s}\bm{\cdot}\left({\bm{\sigma}\textsubscript{aux}}\bm{\cdot}\mathbf{\hat{u}}^{(1)}\right)-\bm{\nabla}_{s}\bm{\cdot}\left(\hat{\bm{\sigma}}^{(1)}\bm{\cdot}\mathbf{u}\textsubscript{aux}\right)dS\\\
=\int{\bm{\sigma}\textsubscript{aux}}\bm{:}\bm{\nabla}\mathbf{\hat{u}}^{(1)}-\hat{\bm{\sigma}}^{(1)}\bm{:}\bm{\nabla}\mathbf{u}\textsubscript{aux}~{}dS.\end{split}$
(28)

The LHS of Eq. (28) can be simplified using the divergence theorem and force
balance to give

$\begin{split}-\int\hat{\mathbf{n}}\bm{\cdot}\bm{\sigma}\textsubscript{aux}\bm{\cdot}\mathbf{\hat{u}}^{(1)}~{}dl+\int\hat{\mathbf{n}}\bm{\cdot}\hat{\bm{\sigma}}^{(1)}\bm{\cdot}\mathbf{u}\textsubscript{aux}~{}dl\\\
=-\mathbf{F}\textsubscript{aux}\bm{\cdot}\mathbf{U\textsuperscript{(1)}}+\mathbf{F\textsuperscript{(1)}\bm{\cdot}U\textsubscript{aux}}.\end{split}$
(29)

Since we constrain the particle to its leading-order motion, the particle
velocity at $\mathcal{O}(\beta)$ vanishes:
$\mathbf{U\textsuperscript{(1)}}=\mathbf{0}$. Then, using Eq. (29) for the LHS
in Eq. (28) and simplifying gives

$\begin{split}\mathbf{F\textsuperscript{(1)}\bm{\cdot}U\textsubscript{aux}}=-\int\left[(\Pi^{(0)}+q)(\mathbf{S}^{(0)}+\mathbf{A})-q\mathbf{A}\right]\bm{:}\bm{\nabla}_{s}\mathbf{u\textsubscript{aux}}~{}dS.\end{split}$
(30)

Importantly, Eq. (30) uses only the zeroth order (Newtonian) solution to solve
a first order (non-Newtonian) problem.

As we are interested in far-field effects in dilute suspensions, we seek to
incorporate leading order effects of the the disturbance velocity due to other
disks as the imposed background flow. Over the lengthscale of a particle, this
disturbance is primarily a unidirectional shear flow (so that
$q(\mathbf{x})=0$). With this assumption, the integration in Eq. (30) is still
cumbersome but can be be performed analytically. Additionally, integrals of
the leading order flow field involve the logarithm that diverges at large
distances as noted before. We therefore restrict our fluid domain to a
(dimensionless) lateral membrane size $R$ that is smaller that or comparable
to the Saffman-Delbrück length: this corresponds to the finite membrane size
resolution of the Stokes paradox Saffman and Delbruck [3]. With the velocity
and pressure fields of the auxiliary Newtonian problem fully known, we finally
get:

$\mathbf{F\textsuperscript{(1)}}=\zeta\mathbf{F}_{p}\bm{\cdot}\mathbf{A},$
(31a)
$\zeta=4\left[\ln\left(R\right)-\frac{1}{4}-\frac{3}{R^{2}}+\frac{11}{4R^{4}}-\frac{1}{R^{6}}\right],$
(31b)

where all quantities are dimensionless. The logarithmic dependence on membrane
size that is characteristic of a 2D calculation of a driven object (unless
otherwise regularized by bulk fluid momentum at large distances) shows up here
and dominates the nonlinear response.

The force that arises due to surface-pressure-dependent surface viscosity
becomes evident if we revert to dimensional variables. Putting everything
together, the leading-order dimensional force on a driven particle in an
ambient linear flow field in this 2D fluid is

$\mathbf{F}_{\text{tot}}=\mathbf{F}_{p}+\frac{\zeta\eta_{s}^{0}}{\Pi_{c}}\mathbf{F}_{p}\bm{\cdot}\mathbf{A}+\mathcal{O}(\beta^{2}).$
(32)

This is the main result of this paper. As expected, the non-Newtonian
correction vanishes as $\Pi_{c}\rightarrow\infty$. More useful in many-
particle simulations or continuum descriptions (§ III) is the net
translational velocity resulting from such a force. Using the leading-order
translational mobility
$\mathbf{M}=\mathbf{I}_{s}/4\pi\eta_{s}^{0}=M_{0}\mathbf{I}_{s}$, and
acknowledging the uniform component $\mathbf{V}_{0}$ of the background
velocity, we find

$\mathbf{U}_{p}=\mathbf{V}_{0}+M_{0}\left[\mathbf{I}_{s}+\frac{\zeta\eta_{s}^{0}}{\Pi_{c}}\mathbf{A}\right]\bm{\cdot}\mathbf{F}_{p}+\mathcal{O}(\beta^{2}).$
(33)

Notably, the correction proportional to $\mathbf{A}\bm{\cdot}\mathbf{F}_{p}$
is not necessarily along the direction of the driving force $\mathbf{F}_{p}$,
thus setting the stage for lateral motion and concentration fluctuations in
large-scale ‘suspensions’ of trans-membrane inclusions.

###  II.5 Physical interpretation and direction of motion

Figure 2: Physical mechanism leading to transverse migration in a driven
particle placed in shear flow within the membrane as shown in (a). Background
shear as shown in (b) sets up an extensional flow as shown (c). The external
force ${\bf F}_{p}$ on the particle also sets up pressure fields that lead to
increase (in red) in viscosity ahead of the particle and decrease (in green)
behind the particle. The resulting in-plane viscous stresses are proportional
to the local viscosity and are asymmetric as indicated with the solid arrows.
This results in a net force ${\bf F}^{(1)}$ to the right in this particular
configuration.

We will first aim to mechanistically understand the direction of motion due to
the non-Newtonian effect. A simple case is a background shear flow that is
along the direction of the driving force (Fig. 2). This geometry is
particularly relevant as it depicts the dominant effect of disturbance fields
from neighboring particles, and is a common setup in classic studies of 3D
sedimenting suspension mechanics [48, 50]. In our case, we will take
$\hat{\mathbf{e}}_{x}$ to be the direction along which a disk embedded within
a membrane is driven with an external horizontal force, so that
$\mathbf{F}_{p}=F_{p}\hat{\mathbf{e}}_{x}$. Simultaneously, the 2D fluid is
subject to an in-plane shear flow
$\mathbf{v(x)}=\dot{\gamma}y\,\hat{\mathbf{e}}_{x}$ centered at the particle,
so that the background strain tensor is

$\mathbf{A}=\bm{\nabla}_{s}\mathbf{v}+\bm{\nabla}_{s}\mathbf{v}^{T}=\dot{\gamma}(\hat{\mathbf{e}}_{x}\hat{\mathbf{e}}_{y}+\hat{\mathbf{e}}_{y}\hat{\mathbf{e}}_{x}).$
(34)

Then, the resulting net velocity following Eq. (33) is

$\mathbf{U}_{p}=M_{0}F_{p}\hat{\mathbf{e}}_{x}+\frac{\zeta\eta_{s}^{0}\dot{\gamma}M_{0}F_{p}}{\Pi_{c}}\,\hat{\mathbf{e}}_{y}.$
(35)

For $\Pi$-thickening suspensions, $\Pi_{c}>0$, and so the disk drifts in the
$y$ direction in addition to being driven as expected in the $x$ direction by
the applied force.

The direction of the non-Newtonian force can be gleaned by examining the
nature of the pressure perturbation and disturbance flow around the disk. A
shear flow can be locally decomposed into an extension and a rotation. The
extension associated with the supplied shear is shown in Fig. 2(c). In a
Newtonian medium, such an extensional flow does not generate a net force on
the particle as surface viscous stresses around the perimeter of the particle
cancel out. In a pressure-thickening fluid, however, the dipolar surface
pressure field generates an associated surface viscosity field $\eta_{s}(\Pi)$
that makes the 2D medium more viscous ahead of the disk in the direction of
$\mathbf{F}_{p}$ and less viscous behind it. The resulting surface viscous
stresses are asymmetric around the particle, leading to a net force due to the
extensional disturbance flow. The local surface viscous stresses are
qualitatively illustrated in Fig. 2(c), which suggests that the resulting net
force is to the right in this configuration.

This intuitive understanding of the interplay between extension due to
background flow and local viscosity variations helps build up to more
sophisticated models for coordinated behavior of collections of proteins or
particles embedded within a membrane. In the following section, we take the
mathematical result for the corrective force $\mathbf{F\textsuperscript{(1)}}$
and the mechanistic picture drawn here to quantify such collective dynamics
using mean-field kinetic models and particle simulations.

##  III Collective dynamics: Concentration Instability and Crystallization

###  III.1 Mean-field model for dilute suspensions

While the model for a single particle in a weakly non-Newtonian membrane is
analytically tractable, the collective behavior of large-scale ‘2D
suspensions’ is much more complicated. Nevertheless, the direction of the
corrective force derived above suggests a mechanism for large-scale
aggregation of particles driven in a membrane. Simply put, the disturbance
flow associated with the motion of each inclusion drives a local shear flow
around every other inclusion. Cross-streamline motion within such a flow must
directly lead to effective hydrodynamic drift towards regions of higher
concentration of particles (Fig. 3). We intend to describe the concentration
instability that results from such interactions.

We will first turn to a continuum description to glean physical features and
stabilizing mechanisms of such an instability. We will follow the classic
approach of Koch and Shaqfeh [50] in studying 3D suspensions of sedimenting
rods, which has since been widely adapted to examine viscoelastic fluids [48],
active suspensions [51], and flexible fibers [52].

We begin with a continuous variable $c(\mathbf{x},t)$ that describes the
concentration of particles embedded within this 2D medium at position
$\mathbf{x}$ at time $t$. We will use $n$ to describe the mean number density
(per unit area) on a interface of area $A$ such that

$\frac{1}{A}\int c(\mathbf{x},t)\,dA=n.$ (36)

Then, the evolution of local concentration is described by a conservation
equation

$\frac{\partial c}{\partial
t}+\bm{\nabla}_{s}\bm{\cdot}(\dot{\mathbf{x}}\,c)=D\nabla_{s}^{2}c,$ (37)

where $\dot{\mathbf{x}}$ is the flux velocity and $D$ is a diffusivity that
may represent Brownian fluctuations or hydrodynamic diffusion. We assume that
the diffusivity is constant and isotropic, although such an assumption can be
relaxed without qualitative changes to the conclusions drawn below. The flux
velocity can be written as a combination of the background velocity field and
the response of the particle to the applied force following Eq. (33):

$\dot{\mathbf{x}}=\mathbf{u}(\mathbf{x})+M_{0}\left[\mathbf{I}_{s}+\frac{\zeta\eta_{s}^{0}}{\Pi_{c}}\left(\bm{\nabla}_{s}\mathbf{u}+\bm{\nabla}_{s}\mathbf{u}^{T}\right)\right]\bm{\cdot}\mathbf{F}_{p}.$
(38)

Written this way, the conservation equation (37) captures the Newtonian
response of each particle to the external force (via the term
$M_{0}\mathbf{F}$, the non-Newtonian correction due to weak pressure-
dependence (via the term proportional to $\zeta$), and the disturbance field
due to the presence of neighboring particles at the location of each particle
(via the background velocity field $\mathbf{u}(\mathbf{x})$). This disturbed
fluid velocity field is still unknown, and obeys

$-\bm{\nabla}_{s}\Pi+\eta_{s}^{0}\nabla_{s}^{2}\mathbf{u}+\frac{\eta_{s}^{0}}{\Pi_{c}}\bm{\nabla}_{s}\bm{\cdot}(\Pi\mathbf{S})+\mathbf{F}_{p}[c(\mathbf{x})-n]=\mathbf{0},$
(39)

which is the 2D momentum equation perturbed by the external force at the
location of the particles.

We will follow past works [48] to first argue that the nonlinear term in Eq.
(39) is sub-dominant in determining the flux velocity if the suspension is
dilute and concentration fluctuations are small. We will stay in the regime of
small perturbations in local concentration away from the uniform value $n$, so
that $c(\mathbf{x},t)=n+\epsilon c^{\prime}(\mathbf{x},t)$, with
$|\epsilon|\ll 1$ and $c^{\prime}(\mathbf{x},t)=\mathcal{O}(n)$. We will also
expand $\mathbf{u}$ as a perturbation in $\beta$, so
$\mathbf{u}=\mathbf{u}^{(0)}+\beta\mathbf{u}^{(1)}+\mathcal{O}(\beta^{2})$
like in previous sections. The leading-order disturbance velocity
$\mathbf{u}^{(0)}$ then has characteristic magnitude $U_{0}\sim\epsilon
F_{p}nL^{2}/\eta_{s}^{0}$ where $L$ is the suspension length scale (typically,
the extent of the membrane or interface). The perturbative correction
$\mathbf{u}^{(1)}$, by design, has characteristic scale $U_{1}\sim\beta
U_{0}$. We can compare this disturbance velocity field correction to the
direct contribution due to the non-Newtonian effect obtained in Eq. (33) which
scales like $U_{c}\sim\beta F_{p}/\eta_{s}^{0}$. So, the nonlinear term in Eq.
(39) is relevant only if $U_{1}$ is comparable to $U_{c}$. However,
$U_{1}/U_{c}\sim U_{0}\eta_{S}^{0}/F_{p}\sim\epsilon nL^{2}$. In a dilute 2D
suspension, the area fraction of particles is small ($nL^{2}\ll 1$) and the
associated perturbation to linear stability is also small ($\epsilon\ll 1$).
The non-Newtonian correction in equation (39) can therefore be safely
neglected, and the disturbance field due to concentration fluctuations follow

$\bm{\nabla}_{s}\Pi=\eta_{s}^{0}\nabla_{s}^{2}\mathbf{u}+\mathbf{F}_{p}[c(\mathbf{x})-n],$
(40)

accurate to $\mathcal{O}(\epsilon\beta)$.

Figure 3: Mechanism for suspension instability. A local increase in particle
concentration drives a shear flow that interacts with the non-uniform
viscosity field around each particle in a manner that draws particles towards
regions of high concentration. This magnifies the concentration fluctuation
and the instability grows.

We wish to examine the linear stability of such a system to perturbations in
concentrations away from the uniform state of $c(\mathbf{x})=n$. The
corresponding base state fluid velocity and surface pressure fields are
$\mathbf{u}=\mathbf{0}$ are $\Pi(\mathbf{x})=\Pi_{0}$, respectively.
Perturbing the concentration as $c=n+\epsilon c^{\prime}(\mathbf{x},t)$ where
$c^{\prime}=\mathcal{O}(n)$ generates associated fluid fields
$\mathbf{u}=\epsilon\mathbf{u}^{\prime}(\mathbf{x},t)$ and
$\Pi=\Pi_{0}+\epsilon\Pi^{\prime}(\mathbf{x},t)$. Using these perturbed fields
in the conservation equation (37) and linearizing gives the following form at
$\mathcal{O}(\epsilon)$:

$\frac{\partial c^{\prime}}{\partial
t}+\frac{nM_{0}\zeta\eta_{s}^{0}}{\Pi_{c}}\nabla_{s}^{2}\mathbf{u}^{\prime}\bm{\cdot}\mathbf{F}_{p}+M_{0}\mathbf{F}_{p}\bm{\cdot}\bm{\nabla}_{s}c^{\prime}=D\nabla_{s}^{2}c^{\prime}.$
(41)

For stability analysis, we will consider normal modes of the form
$c^{\prime}=\tilde{c}(\mathbf{k})\exp[i\mathbf{k}\bm{\cdot}\mathbf{u}+\sigma
t]$ and
$\mathbf{u}^{\prime}=\tilde{\mathbf{u}}(\mathbf{k})\exp[i\mathbf{k}\bm{\cdot}\mathbf{u}+\sigma
t]$ with a 2D wavevector $\mathbf{k}$. The forced Stokes equation can be
solved by applying Fourier transforms and using standard methods of projecting
perpendicular to the pressure term [48, 53, 38] to eliminate
$\tilde{\mathbf{u}}(\mathbf{k})$ in favor of $\tilde{c}(\mathbf{k})$:

$\tilde{\mathbf{u}}(\mathbf{k})=\frac{1}{\eta_{s}^{0}k^{2}}\left(\mathbf{I}_{s}-\hat{\mathbf{k}}\hat{\mathbf{k}}\right)\bm{\cdot}\mathbf{F}_{p}\tilde{c}(\mathbf{k}),$
(42)

where $k=|\mathbf{k}|$, and $\hat{\mathbf{k}}=\mathbf{k}/k^{2}$ is a unit
wavevector. Using these Fourier coefficients in the linearized conservation
equation (41) and simplifying gives

$\begin{split}\sigma\,\tilde{c}(\mathbf{k})-\frac{nM_{0}\zeta}{\Pi_{c}}\mathbf{F}_{p}\bm{\cdot}\left(\mathbf{I}_{s}-\hat{\mathbf{k}}\hat{\mathbf{k}}\right)\bm{\cdot}\mathbf{F}_{p}\,\tilde{c}(\mathbf{k})\\\
+ikM_{0}\mathbf{F}_{p}\bm{\cdot}\mathbf{k}\,\tilde{c}(\mathbf{k})+k^{2}D\,\tilde{c}(\mathbf{k})=0.\end{split}$
(43)

Eliminating $\tilde{c}(\mathbf{k})$ and defining
$\theta=\cos^{-1}(\hat{\mathbf{e}}_{x}\bm{\cdot}\hat{\mathbf{k}})$ as the
angle between the applied external force (taken to be along the $x$ direction,
without loss of generality) and the wave vector, we obtain the real part of
$\sigma$ representing the growth rate of concentration fluctuations:

$\sigma_{R}=\mathrm{Re}[\sigma]=\frac{nM_{0}\zeta
F_{p}^{2}}{\Pi_{c}}\sin^{2}\theta-k^{2}D.$ (44)

Positive $\sigma_{R}$ indicates an exponentially growing perturbation based on
the form of the normal modes. As $\zeta$, $\Pi_{c}$ and $M_{0}$ are all
positive, the system is linearly unstable to concentration fluctuations so
long as $\theta$ is non-zero.

Figure 4: (a) Growth rate of concentration instability, nondimensionalized as
$\sigma_{R}^{*}=\sigma\Pi_{c}/nM_{0}\zeta F_{p}^{2})$, as a function of
disturbance wavenumber $k^{*}=k/\sqrt{n}$. Curves represent different
orientations of the wavevector relative to
$\mathbf{F}_{p}\sim\hat{\mathbf{e}}_{x}$, with
$\cos\theta=\hat{\mathbf{e}}_{x}\bm{\cdot}\hat{\mathbf{k}}$. (b) Diffusion
acts to stabilize large wavenumber fluctuations, and a higher diffusivity
suppresses a wider range of modes: here, $D^{*}=D\Pi_{c}/M_{0}\zeta
F_{p}^{2}$.

Figure 4 shows the dispersion relation from Eq. (44) made dimensionless by
scaling the growth rate over $nM_{0}\zeta F_{p}^{2}/\Pi_{c}$, the wave number
over $\sqrt{n}$ and the diffusion constant over $M_{0}\zeta
F_{p}^{2}/\Pi_{c}$. Such a driven suspension is evidently unstable to long-
wavelength perturbations. Indeed, the growth rate is always largest at $k=0$,
corresponding to concentration fluctuations that span the system size.
Perturbations in the direction perpendicular to the driving force (so that
$\theta=\pi/2$) are the most destabilizing, analogous to classic theories of
sedimentation stability in 3D suspensions [50, 51, 52]. The only stabilizing
mechanism is diffusion, that suppresses large wave numbers (as shown in Fig.
4b). As expected, the instability does not occur in the Newtonian limit as
$\Pi_{c}\rightarrow\infty$.

The associated mechanism is illustrated in Fig. 3. A perturbation that
increases local particle concentration drives a shear flow that interacts with
the non-uniform viscosity field around each particle to generate a non-
Newtonian correction to the drift velocity as shown in previous sections (see
Fig. 2). In a suspension, this drift draws particles towards regions of high
concentration. This magnifies the concentration fluctuation and the
instability grows.

###  III.2 Langevin description

The mean-field model above describes short-term linear growth of concentration
fluctuations and suggests a mechanism for aggregation. In what follow, we will
demonstrate the long-term effect of the nonlinear force on collections of
membrane inclusions using a Langevin description. We will consider $N$ disks
embedded in a membrane with an infinite sub-phase, where each disk translates
with a self velocity due to an applied force $\mathbf{F}_{p}$. This could be
due to anchoring proteins on a membrane, or equivalently due to externally
applied electrical, thermal, or optical forces on particles embedded within a
complex interface. We assume that the suspension is sufficiently dilute such
that pairwise interactions can be added to obtain net force on each disk. The
trajectory of each disk is then determined by the Newtonian membrane
hydrodynamic disturbance field due to the other $N-1$ disks, the correction
due to non-Newtonian nature of the membrane as derived above, and Brownian
fluctuations.

The corresponding Langevin equation for disk $i$ reads:

$\frac{\partial\mathbf{x}_{i}}{\partial
t}=M_{0}F_{p}\hat{\mathbf{e}}_{x}+\mathbf{u}_{i}^{\text{Br}}+\sum_{i\neq
j}^{N}\left[\mathbf{u}^{d}_{j}(\mathbf{x}_{i})+\mathbf{u}^{\text{St}}_{ij}+M_{0}\beta\mathbf{F}_{ij}\right].$
(45)

The first term on the RHS describes the translational response of disk $i$ to
the applied force. We will take $\mathbf{F}_{p}=F_{p}\hat{\mathbf{e}}_{x}$ to
be constant and along the $x$ direction. The second term describes Brownian
fluctuations, which in general depend on the collective hydrodynamic motility
of the system [14]. Owing to the diluteness of our system, we only account for
the local mobilities for each disk. This simplification decouples the
fluctuation-dissipation relation of each disk from the rest so that the
translational Brownian velocity satisfies:

$\displaystyle\langle\mathbf{u}^{\text{Br}}(t)\rangle=0,\quad\langle\mathbf{u}^{\text{Br}}(t)\mathbf{u}^{\text{Br}}(t^{\prime})\rangle=2k_{B}T\mathbf{M}\delta(t-t^{\prime}),$
(46)

where $\mathbf{M}=M_{0}\mathbf{I}_{s}=\mathbf{I_{s}}/4\pi\eta_{s}^{0}$ is the
Newtonian translational mobility of a disk. Using $U=F_{p}/\eta_{s}^{0}$ and
$a/U$ as the characteristic scales for velocity and time, we obtain the
dimensionless Brownian velocity from Eq. (46) as

$\mathbf{u}_{i}^{\text{Br}}=\sqrt{\frac{\tilde{T}}{2\pi\Delta\tilde{t}}}~{}\mathbf{w},\qquad\tilde{T}=\frac{k_{B}T}{F_{p}a},$
(47)

where $\Delta\tilde{t}$ is the dimensionless time step and $\tilde{T}$ is a
dimensionless temperature. The white noise vector $\mathbf{w}$ is populated by
random numbers sampled from a normal distribution of mean 0 and variance 1
such that the fluctuation-dissipation relation Eq. (46) is satisfied.

The next term in the Langevin equation describes leading-order disturbance
field at the location of disk $i$ due to an external force on neighboring disk
$j$:

$\mathbf{u}^{d}_{j}(\mathbf{x}_{i})=\mathbf{G}_{ij}(\mathbf{x}_{i}-\mathbf{x}_{j})\bm{\cdot}\mathbf{F}_{p}.$
(48)

Here, $\mathbf{G}_{ij}$ is the Green’s function or 2D stokeslet for the
Newtonian problem. This tensor is readily obtained from the term corresponding
to $\mathbf{F}_{p}$ in Eq. (21).

Although the suspension is dilute and particles are initially far from each
other, fluctuations due to pair-wise addition of attractive hydrodynamic
forces can bring disks within contact distance of each other. To prevent
overlap, we add a soft repulsion between disks at short range. These repulsive
excluded volume interactions are accounted via a steric velocity

$\mathbf{u}^{\text{St}}_{ij}=U_{s}\frac{e^{-\alpha(r-a)}}{1+e^{-\alpha(r-a)}}\,\hat{\mathbf{r}}_{ij},$
(49)

where $r=|\mathbf{x}_{i}-\mathbf{x}_{j}|$, and which acts along the unit
vector $\hat{\mathbf{r}}_{ij}=(\mathbf{x}_{i}-\mathbf{x}_{j})/r$ that connects
disk $i$ to disk $j$. This soft repulsion decays exponentially over a length
scale $\alpha^{-1}=\mathcal{O}(a)$ and $U_{s}$ is the contact velocity.

Figure 5: (a) Local extension field around particles due to the flow generated
by a reference particle at the center. (b) Effective non-Newtonian interaction
due to the interplay between local extension and non-uniform viscosity reveals
a propensity to chain particles along the direction of the applied force. (c)
Snapshots at equal intervals of time of a pair of particles translating in
response to a downward force, simulated using the Langevin description. Dashed
unfilled particles are the Newtonian reference case, which maintain relative
separation and orientation. Solid filled circles account for pressure-
dependent viscosity, which draws them together to chain up along the direction
of the force.

The final term in the Langevin description is the non-Newtonian correction
from Eq. (33) obtained via the reciprocal calculation in section II. This term
also acts pair-wise as the local background velocity gradient felt by disk $i$
is that generated by the disturbance velocity associated with disk $j$. We
time integrate the Langevin equation numerically using standard Brownian
dynamics methods established for similar systems [20, 21].

###  III.3 Effective hydrodynamic interactions and chaining

Before discussing long-time features and hydrodynamic crystallization, we
first take a closer look at the effective hydrodynamic interactions in large-
scale suspensions to understand the mechanism of particle aggregation. The
interplay between local fluid extension and non-uniform viscosity field around
each particle dictates the effective hydrodynamic interactions that drive
relative particle motion. Fig. 5(a) borrows from the classic work of Koch and
Shaqfeh [50] to show local extension around a particle that is placed near a
point force (corresponding to a reference particle) at the center. Each
surrounding particle represents a possible location of a membrane anchor
protein in the disturbance field induced by the reference particle. The
orientation of the extensional field changes based on the location of the
particle, but always acts to draw fluid towards the reference point force in
the upper half and to draw fluid away from it in the lower half.

Recalling that inclusions are also driven downward by external forces, the
local viscosity field around each inclusion is simultaneously weakly
perturbed. The particles experience a relatively higher viscosity ahead of
them. The effective hydrodynamic force, or the resulting drift velocity, can
then be determined by evaluating the asymmetric surface viscous stresses
around each particle. Again, such extensional flows would not generate a net
force on disks in a Newtonian system (akin to spheres in 3D fluids [50]): the
non-uniform viscosity around each inclusion breaks this symmetry.

We can readily evaluate the form of such an interaction. The flow field set up
by the reference particle (the dashed lines in Fig. 5a) is given by the term
proportional to the force in Eq. (21). For illustration, we consider only the
leading-order disturbance field corresponding to a 2D stokeslet or a point
force. The correction due to the interplay between the external forces and the
strain rate corresponding to such a flow at a location $\mathbf{x}$ is given,
following Eq. (31), by

$\mathbf{F\textsuperscript{(1)}}=\zeta\mathbf{F}_{p}\bm{\cdot}\mathbf{A}=\zeta\mathbf{F}_{p}\bm{\cdot}\left[\frac{\mathbf{I}_{s}\mathbf{x}}{r^{2}}-\frac{2\mathbf{x}\mathbf{x}\mathbf{x}}{r^{4}}\right]\bm{\cdot}\mathbf{F}_{p},$
(50)

where $r=|\mathbf{x}|$. Taking
$\mathbf{F\textsuperscript{(1)}}=F_{p}\hat{\mathbf{e}}_{x}$ and denoting the
radial unit vector originating at the reference particle as
$\hat{\mathbf{e}}_{r}$, Eq. (50) becomes

$\mathbf{F\textsuperscript{(1)}}=\zeta
F_{p}^{2}\left(\frac{\cos\theta}{r}\hat{\mathbf{e}}_{x}-\frac{2\cos^{2}\theta}{r}\hat{\mathbf{e}}_{r}\right),$
(51)

where now $\cos\theta=\hat{\mathbf{e}}_{x}\bm{\cdot}\hat{\mathbf{e}}_{r}$.

This force leads to a velocity $M_{0}\mathbf{F\textsuperscript{(1)}}$ at the
location $\mathbf{x}$, which can be interpreted as an effective hydrodynamic
interaction due to pressure-dependent viscosity. $M_{0}$ is the isotropic
mobility of a disk, and the direction of the velocity follows that of
$\mathbf{F\textsuperscript{(1)}}$. This correction to the Newtonian-flow field
around each particle, following Eq. (51), is shown in Fig. 5(b). As expected,
this correction vanishes in the horizontal line where $\theta=\pi/2$, where
the extensional field due to a stokeslet is zero. The correction is non-zero
at all other orientations and maintains a front-back symmetry relative to
$\hat{\mathbf{e}}_{x}$. Notably, the direction of
$\mathbf{F\textsuperscript{(1)}}$ is such that it acts to draw neighboring
particles into chains along the direction of external force.

Figure 6: Time series of a 2D suspension of disks driven by an external
(downward) force leading to hydrodynamic aggregation. Chains form initially
which eventually give way to large scale hexatically ordered crystals. Here,
$\beta=0.1$ and $\tilde{T}=0.001$.

Such a mechanism of chaining is clearly seen in Fig. 5(c) in two-particle
simulations following the Langevin description from the preceding section.
While a pair of particles in a 2D Newtonian fluid would maintain relative
separation and orientation as a they translate in response to an external
force, the non-Newtonian correction derived here causes them to line up along
the direction of the external force. We therefore expect particles to form
chains. By analogy with molecular systems that display a long-ranged
attraction and short-range repulsion, we expect clustering and aggregation
when a large number of particles are present.

Indeed, our simulations show that chains form initially but eventually give
way to larger-scale stable clusters (Fig. 6) with increasing $\beta$.
Aggregates are stable and do not disintegrate for long simulation times (as
observed till $\tilde{t}=120$ in units of time made dimensionless over
characteristic time $a\eta_{s}^{0}/F_{p}$, Fig. 6). We see crystalline order
emerge as expected from analogous problems with ‘hydrosteric’ interactions
[21], i.e, long-ranged attraction and short-rainge repulsive interactions.
Previous works have examined aggregation of driven transmembrane proteins due
membrane curvature [35, 36], and inclusion size mismatch [30]. Our simulations
show that stable hexatic clusters form as a result of the hydrodynamic
attraction, physically representing a collection of membrane anchor points
crystallizing for a potential biological advantage.

###  III.4 Aggregation and hexatic Order

We will use particle pair proximity and resulting packing order to quantify
the degree of aggregation in this system. A straightforward metric is the
$n$-th packing order parameter:

$\Psi^{j}_{n}=\frac{1}{n_{j}}\sum_{k}e^{in\theta_{kj}}.$ (52)

$\Psi^{j}_{n}$ measures the orientation and packing order around particle $j$,
where $n_{j}$ is the number of nearest neighbors and $\theta_{kj}$ is the
angle between the line joining disks $k,j$ and the x-axis. The average local
$\langle\absolutevalue{\Psi_{n}}\rangle$ order parameter quantifies the
packing order, with 0 corresponding to an unordered aggregate and 1
representing a perfect $n$-th order aggregate. The local order parameter tells
us whether each particle in an aggregate forms a $n$-th order lattice with its
nearest neighbor. Note that the _global_ hexatic order parameter
$|\langle\Psi_{n}\rangle|$ can indicate system-spanning crystalline
structure[21]. Our simulations are non-periodic and have a finite number of
particles, and we confirm that no qualitative difference exists between the
local and global order parameters for systems at the scales shown.

Given the 2D nature of our systems and the short-ranged repulsion, we always
see hexatic ($n=6$) order formation for non-zero $\beta$. This is comparable
to other works on 2D nonlinear hydrodynamics [21, 54, 55] that have shown that
rigid inclusions rotating on a membrane form stable hexagonal crystals. The
effective attractive potential in our case is orientation-dependent (Fig. 5b)
and arises due to surface-pressure dependent rheology. As expected, disks
self-assemble faster with increasing $\beta$ (Fig. 7) corresponding to a
stronger role of pressure-dependent surface viscosity.

Figure 7: A stronger pressure-dependence of viscosity (smaller $\Pi_{c}$ or
larger $\beta$) leads to faster aggregation and crystallization. Snapshots
taken at $\tilde{t}=25$.

We can also use the evolution of $\langle\absolutevalue{\Psi_{6}}\rangle$ to
quantitatively describe the rate of assembly, and qualitatively describe the
transition towards hexatic order. Hexagonal order of the disk crystals
plateaus with time to attain a quasi-steady value (Fig. 8a). The rate at which
this plateau is reached, as well the magnitude of the plateau are clear
indicators of departure from Newtonian behavior. The transition from disorder
to crystalline order due to $\Pi$-dependent $\eta_{s}$ is evident in Fig.
8(b).

Finally, note that the strong thermal fluctuations can melt crystals and
restabilize the suspension towards a uniform concentration. We have set the
effective temperature $\tilde{T}=k_{B}T/F_{p}a=0.001$ in the results shown
here, which corresponds to $F_{p}=\mathcal{O}(\text{pN})$ for particles or
membrane inclusions in the size range of $a=\mathcal{O}(100\text{ nm})$. This
is in the range of forces exerted on membranes by biological activity such as
intracellular pinning and cytoskeleton-associated growth and locomotive forces
[29, 56]. Larger forces or bigger particles would only accentuate the effect,
and we have tested that crystallization transition is qualitatively unchanged
upon increasing thermal noise upto $\tilde{T}=0.01$. Increasing thermal
fluctuations within this biologically relevant range slightly decreases the
plateau value of $\langle\absolutevalue{\Psi_{6}}\rangle$: the rate at which
the plateau value is reached for all $\beta$, however, remains unchanged.

Figure 8: (a) Evolution of the hexatic order parameter with time. Darker hue
indicates a larger $\beta$, and the curves correspond to $\beta=0.01$ to
$\beta=0.1$ in steps of $0.01$. (b) Plateau value of
$\langle|\Psi_{6}|\rangle$ at $\tilde{t}=30$.

##  IV Conclusion

Collective motility, aggregation, and crystallization of membrane proteins has
distinct biological advantage in immune response, locomotion, and tubulation
[30, 29, 35]. Previous works have examined aggregation that arises due to
protein activity [21, 20], membrane curvature [35, 36], and inclusion size
mismatch [30]. We have shown that a commonly measured class of non-Newtonian
membrane viscosity can have the same effect, and that these non-trivial fluid
dynamics can play critical roles in membrane organization principles [24, 25].

We have made several approximations to keep this complicated problem
tractable. Neglecting the coupling to the bulk phase will have quantitative
changes to the corrective force, especially for large membranes where 3D fluid
stresses catch up over the Saffman-Delbrück crossover length. We expect this
change to be minimal for finite membrane sizes or crowded system, as the bulk
flow effect is screened at much shorter length scales, and so the qualitative
picture remains unchanged. A more careful analysis of the coupling of a non-
Newtonian membrane to a 3D fluid phase is left as future work. We have also
chosen to develop a minimal point-particle-based Langevin approach for
simulations. Our objective was to demonstrate the long-time dynamics with a
minimal simulation model. We expect that high-fidelity computations will
confirm these dynamics, and also go beyond the weakly nonlinear regime. We
also recognize that not all inclusions within a membrane are driven: a real
membrane hosts a collection of anchored, driven, active, and passive
inclusions. Nor are inclusions disk-like. While shape does not significantly
alter long-ranged membrane hydrodynamics, the near-field steric effects and
packing order will shift to accommodate inclusion shape. Nonetheless, we
expect that the physical intuition gained in this study will port over to
future analyses of combinations of mobile and immobile inclusions

We end by noting that the mathematical framework and physical picture that
emerges in this work is generalizable to a broader class of problems beyond
lipids with pressure-dependent viscosity. The formulation using the reciprocal
theorem is fairly general and can be extended to other common non-Newtonian
constitutive relations. Force corrections due to shear-thinning [49] and
viscoelastic effects [48] follow their 3D analogs. The mean-field model can
also be readily adapted to any form of the corrective force, opening a rich
avenue of problems that can be viewed through the lens of this framework.
Further, this model can be applied to a broader class of 3D fluids with
pressure-dependent viscosity, such as those that appear in piezoviscous
problems [57], high-pressure polymer melts [58], geophysical flows [59], and
crude oil mixtures [60].

Conflict of Interest: The authors have no conflicts to disclose.

Acknowledgments: We thank Ronald J. Phillips and Gregory H. Miller for
insightful feedback on this work.

## References

  * Phillips  _et al._ [2012] R. Phillips, J. Kondev, J. Theriot, and H. Garcia,  _Physical biology of the cell_ (Garland Science, 2012).
  * Van Meer, Voelker, and Feigenson [2008] G. Van Meer, D. R. Voelker, and G. W. Feigenson, “Membrane lipids: where they are and how they behave,” Nature reviews Molecular cell biology 9, 112–124 (2008).
  * Saffman and Delbruck [1975] P. G. Saffman and M. Delbruck, “Brownian motion in biological membranes.” Proceedings of the National Academy of Sciences 72, 3111–3113 (1975).
  * Hughes, Pailthorpe, and White [1981] B. D. Hughes, B. A. Pailthorpe, and L. R. White, “The translational and rotational drag on a cylinder moving in a membrane,” Journal of Fluid Mechanics 110, 349 (1981).
  * Choi  _et al._ [2011] S. Q. Choi, S. Steltenkamp, J. A. Zasadzinski, and T. M. Squires, “Active microrheology and simultaneous visualization of sheared phospholipid monolayers.” Nature communications 2, 312 (2011).
  * Espinosa  _et al._ [2011] G. Espinosa, I. López-Montero, F. Monroya, and D. Langevin, “Shear rheology of lipid monolayers and insights on membrane fluidity,” Proceedings of the National Academy of Sciences of the United States of America 108, 6008–6013 (2011).
  * Samaniuk and Vermant [2014] J. R. Samaniuk and J. Vermant, “Micro and macrorheology at fluid–fluid interfaces,” Soft Matter 10, 7023–7033 (2014).
  * Gambin  _et al._ [2006] Y. Gambin, R. Lopez-Esparza, M. Reffay, E. Sierecki, N. S. Gov, M. Genest, R. S. Hodges, and W. Urbach, “Lateral mobility of proteins in liquid membranes revisited,” Proceedings of the National Academy of Sciences 103, 2098–2102 (2006).
  * Cicuta, Keller, and Veatch [2007] P. Cicuta, S. L. Keller, and S. L. Veatch, “Diffusion of liquid domains in lipid bilayer membranes,” Journal of Physical Chemistry B 111, 3328–3331 (2007).
  * Schütte  _et al._ [2017] O. M. Schütte, I. Mey, J. Enderlein, F. Savić, B. Geil, A. Janshoff, and C. Steinem, “Size and mobility of lipid domains tuned by geometrical constraints,” Proceedings of the National Academy of Sciences 114, E6064–E6071 (2017).
  * Prasad, Koehler, and Weeks [2006] V. Prasad, S. Koehler, and E. R. Weeks, “Two-particle microrheology of quasi-2d viscous systems,” Physical Review Letters 97, 176001 (2006).
  * Pepicelli  _et al._ [2019] M. Pepicelli, N. Jaensson, C. Tregouët, B. Schroyen, A. Alicke, T. Tervoort, C. Monteux, and J. Vermant, “Surface viscoelasticity in model polymer multilayers: From planar interfaces to rising bubbles,” Journal of Rheology 63, 815–828 (2019).
  * Stone and Ajdari [1998] H. A. Stone and A. Ajdari, “Hydrodynamics of particles embedded in a flat surfactant layer overlying a subphase of finite depth,” Journal of Fluid Mechanics 369, 151–173 (1998).
  * Oppenheimer and Diamant [2009] N. Oppenheimer and H. Diamant, “Correlated Diffusion of Membrane Proteins and Their Effect on Membrane Viscosity,” Biophysical Journal 96, 3041–3049 (2009).
  * Camley and Brown [2013] B. A. Camley and F. L. Brown, “Diffusion of complex objects embedded in free and supported lipid bilayer membranes: role of shape anisotropy and leaflet structure,” Soft Matter 9, 4767–4779 (2013).
  * Fischer [2004a] T. M. Fischer, “The drag on needles moving in a Langmuir monolayer,” Journal of Fluid Mechanics 498, 123–137 (2004a).
  * Levine, Liverpool, and MacKintosh [2004] A. J. Levine, T. Liverpool, and F. C. MacKintosh, “Dynamics of rigid and flexible extended bodies in viscous films and membranes,” Physical review letters 93, 038102 (2004).
  * Shi, Moradi, and Nazockdast [2022] W. Shi, M. Moradi, and E. Nazockdast, “Hydrodynamics of a single filament moving in a spherical membrane,” Physical Review Fluids 7, 084004 (2022).
  * Mikhailov and Kapral [2015] A. S. Mikhailov and R. Kapral, “Hydrodynamic collective effects of active protein machines in solution and lipid bilayers,” Proceedings of the National Academy of Sciences 112, E3639–E3644 (2015).
  * Manikantan [2020] H. Manikantan, “Tunable Collective Dynamics of Active Inclusions in Viscous Membranes,” Physical Review Letters 125, 268101 (2020).
  * Oppenheimer, Stein, and Shelley [2019] N. Oppenheimer, D. B. Stein, and M. J. Shelley, “Rotating membrane inclusions crystallize through hydrodynamic and steric interactions,” Physical Review Letters 123, 148101 (2019).
  * Ramachandran  _et al._ [2011] S. Ramachandran, S. Komura, K. Seki, and G. Gompper, “Dynamics of a polymer chain confined in a membrane,” The European Physical Journal E 34, 46 (2011).
  * Kaganer, Möhwald, and Dutta [1999] V. M. Kaganer, H. Möhwald, and P. Dutta, “Structure and phase transitions in langmuir monolayers,” Rev. Mod. Phys. 71, 779–819 (1999).
  * Lingwood and Simons [2010] D. Lingwood and K. Simons, “Lipid Rafts As a Membrane-Organizing Principle,” Science 327, 46–50 (2010).
  * Sezgin  _et al._ [2017] E. Sezgin, I. Levental, S. Mayor, and C. Eggeling, “The mystery of membrane organization: Composition, regulation and roles of lipid rafts,” Nature Reviews Molecular Cell Biology 18, 361–374 (2017).
  * Arriaga  _et al._ [2010] L. R. Arriaga, I. López-Montero, J. Ignés-Mullol, and F. Monroy, “Domain-growth kinetic origin of nonhorizontal phase coexistence plateaux in langmuir monolayers: compression rigidity of a raft-like lipid distribution,” The Journal of Physical Chemistry B 114, 4509–4520 (2010).
  * Kurnaz and Schwartz [1997] M. L. Kurnaz and D. K. Schwartz, “Channel flow in a Langmuir monolayer: Unusual velocity profiles in a liquid-crystalline mesophase,” Physical Review E 56, 3378–3384 (1997).
  * Kim  _et al._ [2011] K. Kim, S. Q. Choi, J. A. Zasadzinski, and T. M. Squires, “Interfacial microrheology of DPPC monolayers at the air–water interface,” Soft Matter 7, 7782 (2011).
  * Fogelson and Mogilner [2014] B. Fogelson and A. Mogilner, “Computational estimates of membrane flow and tension gradient in motile cells,” PloS One 9, e84524 (2014).
  * Bruinsma and Pincus [1996] R. Bruinsma and P. Pincus, “Protein aggregation in membranes,” Current Opinion in Solid State and Materials Science 1, 401–406 (1996).
  * Manikantan and Squires [2017a] H. Manikantan and T. M. Squires, “Pressure-dependent surface viscosity and its surprising consequences in interfacial lubrication flows,” Physical Review Fluids 2, 023301 (2017a).
  * Manikantan and Squires [2017b] H. Manikantan and T. M. Squires, “Irreversible particle motion in surfactant-laden interfaces due to pressure-dependent surface viscosity,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 473, 20170346 (2017b).
  * Singh and Narsimhan [2021] N. Singh and V. Narsimhan, “Impact of surface viscosity on the stability of a droplet translating through a stagnant fluid,” Journal of Fluid Mechanics 927, A44 (2021).
  * Herrada  _et al._ [2022] M. A. Herrada, A. Ponce-Torres, M. Rubio, J. Eggers, and J. M. Montanero, “Stability and tip streaming of a surfactant-loaded drop in an extensional flow. influence of surface viscosity,” Journal of Fluid Mechanics 934, A26 (2022).
  * Reynwar  _et al._ [2007] B. J. Reynwar, G. Illya, V. A. Harmandaris, M. M. Müller, K. Kremer, and M. Deserno, “Aggregation and vesiculation of membrane proteins by curvature-mediated interactions,” Nature 447, 461–464 (2007).
  * Dasgupta, Auth, and Gompper [2017] S. Dasgupta, T. Auth, and G. Gompper, “Nano- and microparticles at fluid and biological interfaces,” Journal of Physics: Condensed Matter 29, 373003 (2017).
  * Slattery and Sagis [2007] J. C. Slattery and L. Sagis,  _Interfacial Transport Phenomena._ (Springer, 2007).
  * Manikantan and Squires [2020] H. Manikantan and T. M. Squires, “Surfactant dynamics: hidden variables controlling fluid flows,” Journal of Fluid Mechanics 892, P1 (2020).
  * Fischer [2004b] T. M. Fischer, “Comment on “shear viscosity of langmuir monolayers in the low-density limit”,” Physical Review Letters 92, 139603 (2004b).
  * Stone and Masoud [2015] H. A. Stone and H. Masoud, “Mobility of membrane-trapped particles,” Journal of Fluid Mechanics 781, 494–505 (2015).
  * Alonso and Zasadzinski [2006] C. Alonso and J. A. Zasadzinski, “A brief review of the relationships between monolayer viscosity, phase behavior, surface pressure, and temperature using a simple monolayer viscometer,” The Journal of Physical Chemistry B 110, 22185–22191 (2006).
  * Zell  _et al._ [2014] Z. A. Zell, A. Nowbahar, V. Mansard, L. G. Leal, S. S. Deshmukh, J. M. Mecca, C. J. Tucker, and T. M. Squires, “Surface shear inviscidity of soluble surfactants,” Proceedings of the National Academy of Sciences 111, 3677–3682 (2014).
  * Kim  _et al._ [2013] K. Kim, S. Q. Choi, Z. A. Zell, T. M. Squires, and J. A. Zasadzinski, “Effect of cholesterol nanodomains on monolayer morphology and dynamics,” Proceedings of the National Academy of Sciences 110, E3054–E3060 (2013).
  * Leal [2007] L. G. Leal,  _Cambridge Series in Chemical Engineering_ (Cambridge University Press, Cambridge, 2007).
  * Pozrikidis [1992] C. Pozrikidis,  _Boundary integral and singularity methods for linearized viscous flow_ (Cambridge university press, 1992).
  * Kim and Karrila [1991] S. Kim and S. J. Karrila,  _Microhydrodynamics: principles and selected applications_ (Butterworth-Heinemann, 1991).
  * Masoud and Stone [2019] H. Masoud and H. A. Stone, “The reciprocal theorem in fluid dynamics and transport phenomena,” Journal of Fluid Mechanics 879, P1 (2019).
  * Vishnampet and Saintillan [2012] R. Vishnampet and D. Saintillan, “Concentration instability of sedimenting spheres in a second-order fluid,” Physics of Fluids 24, 073302 (2012).
  * Chen  _et al._ [2021] Y. Chen, E. Demir, W. Gao, Y.-N. Young, and O. S. Pak, “Wall-induced translation of a rotating particle in a shear-thinning fluid,” Journal of Fluid Mechanics 927, R2 (2021).
  * Koch and Shaqfeh [1989] D. L. Koch and E. S. Shaqfeh, “The instability of a dispersion of sedimenting spheroids,” Journal of Fluid Mechanics 209, 521–542 (1989).
  * Saintillan and Shelley [2008] D. Saintillan and M. J. Shelley, “Instabilities and pattern formation in active particle suspensions: kinetic theory and continuum simulations,” Physical Review Letters 100, 178103 (2008).
  * Manikantan  _et al._ [2014] H. Manikantan, L. Li, S. E. Spagnolie, and D. Saintillan, “The instability of a sedimenting suspension of weakly flexible fibres,” Journal of Fluid Mechanics 756, 935–964 (2014).
  * Hasimoto [1959] H. Hasimoto, “On the periodic fundamental solutions of the stokes equations and their application to viscous flow past a cubic array of spheres,” Journal of Fluid Mechanics 5, 317–328 (1959).
  * Grzybowski, Stone, and Whitesides [2000] B. A. Grzybowski, H. A. Stone, and G. M. Whitesides, “Dynamic self-assembly of magnetized, millimetre-sized objects rotating at a liquid–air interface,” Nature 405, 1033–1036 (2000).
  * Goto and Tanaka [2015] Y. Goto and H. Tanaka, “Purely hydrodynamic ordering of rotating disks at a finite reynolds number,” Nature communications 6, 5994 (2015).
  * Cohen and Shi [2020] A. E. Cohen and Z. Shi, “Do cell membranes flow like honey or jiggle like jello?” BioEssays 42, 1900142 (2020).
  * Hron, Malek, and Rajagopal [2001] J. Hron, J. Malek, and K. R. Rajagopal, “Simple flows of fluids with pressure-dependent viscosities,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 457, 1603–1622 (2001).
  * Denn [1981] M. M. Denn, “Pressure drop-flow rate equation for adiabatic capillary flow with a pressure- and temperature-dependent viscosity,” Polymer Engineering and Science 21, 65–68 (1981).
  * Schoof [2007] C. Schoof, “Pressure-dependent viscosity and interfacial instability in coupled ice–sediment flow,” Journal of Fluid Mechanics 570, 227 (2007).
  * Srinivasan, Bonito, and Rajagopal [2013] S. Srinivasan, A. Bonito, and K. R. Rajagopal, “Flow of a fluid through a porous solid due to high pressure gradients,” Journal of Porous Media 16, 193–203 (2013).