Statistical Mobility of Multicellular Colonies of Flagellated Swimming Cells

Yonatan L Ashenafi, Peter R Kramer

1 Department of Mathematical and Statistical Sciences, University of Alberta,
Edmonton, AB, Canada
2 Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy,
NY, USA

*<EMAIL_ADDRESS>
## Abstract

We study the stochastic hydrodynamics of colonies of flagellated swimming
cells, typified by multicellular choanoflagellates, which can form both
rosette and chainlike shapes. The objective is to link cell-scale dynamics to
colony-scale dynamics for various colonial morphologies. Via autoregressive
stochastic models for the cycle-averaged flagellar force dynamics and
statistical models for demographic cell-to-cell variability in flagellar
properties and placement, we derive effective transport properties of the
colonies, including cell-to-cell variability. We provide the most quantitative
detail on disclike geometries to model rosettes, but also present formulas for
the dynamics of general planar colony morphologies, which includes planar
chain-like configurations.

##  1 Introduction

Various forms of eukaryotic cells ranging from protozoa to spermatozoa employ
flagellar beating to navigate in their fluid environment, for example in
search of food, chemicals, and other cells while also avoiding their
predators. Experimental work has investigated vital components of flagella-
driven swimming motion such as those concerning fluid rheology (Rafaï et al.,
2010; Matsui et al., 2020), strength of the propulsive forcing due to flagella
(Bayly et al., 2011), the statistics of the orientation of the flagella from
its base (Kirkegaard et al., 2016a), and the ability to sense and react to
various environmental cues (Kirkegaard et al., 2016b; DeMarco et al., 2020;
Foster and Smyth, 1980). Theoretical models and analyses have made predictions
on the flows and mobility mechanisms of the swimming cells, including
properties in suspension (Roper et al., 2013; Blake, 1974; Dean, 1996;
Sparacino et al., 2020).

Some recent experimental studies on swimming microorganisms have turned to
multicellular colonies driven by multiple flagella, for example in protozoa
(Kirkegaard et al., 2016a; de Maleprade et al., 2020) and green algae
(Goldstein, 2015). The present work is largely motivated by the theoretical
analysis in Kirkegaard et al. (2016a) for the observed dynamics of protozoan
choanoflagellate colonies. These are colonial eukaryotic organisms which are
close relatives of animals, and are being used to understand the pathways to
emergence of mutlicelluarity in animals (Brunet and King, 2017; King, 2004;
Koehl, 2020; Kumler et al., 2020). Choanoflagellate colonies mimic their
distant mobile predecessors and choanoflagellate cells are also shown to be
quite similar to sponges found in the ocean (Leadbeater, 2015).

In Kirkegaard et al. (2016a), the theoretical connection between the
individual cellular behavior and the effective “aggregate random walker” model
of translational and rotational dynamics for the colony is made through
phenomenology and scaling arguments. We endeavor here to derive more
systematically the colony dynamics beginning from a formulation of specific
statistically independent dynamical models for each flagellum together with a
rigid-body model for their geometric configuration. A key aspect of the
swimming performance of colonies is that in both rosette and chain
configurations, the flagella of the constituent cells are oriented in various
directions and therefore partially or even largely cancel each other other’s
propulsive force on the colony (Koehl, 2020). Irregularities in the flagellar
positioning introduce qualitative changes to the swimming of a colony by
mitigating the cancellation of flagellar forces and inducing a net torque
causing the colony to rotate (Koehl, 2020; Laundon et al., 2019). Roper et al.
(2013) studied feeding fluxes in regular models of colony shapes by
consideration of the flows generated by the geometric arrangement of flagellar
forces, with brief consideration of irregular flagellar placement. One of our
central objectives is to formulate stochastic models in the flagellar dynamics
and statistical models for the flagellar geometry, and to examine in a precise
manner how these two sources of randomnness together influence the mobility of
the colony.

In Section 2, we develop the model of a disclike colony that is constructed in
the image of choanoflagellate rosettes. The cycle-averaged flagellar forces
and geometry are modeled stochastically to account for both observed dynamical
variations in time and demographic variations between flagella. We discuss as
well how we estimate the model parameters from the experimental literature.
Next in Section 3 we summarize the results of our calculations for some key
mobility characterizations of the colony: mean-square speed as well as the
effective translational and rotational drift and diffusity. We present
formulas both for an individual colony with its peculiar irregular
configuration of flagella as well as for population statistics over an
ensemble of colonies with a prescribed demographic distribution of variability
in flagellar properties and placement on the colony. The derivation of these
results, which amount to probabilistic calculations and some asymptotic
approximations, are explained in Section 4. In Section 5, we show how to
generalize our characterization of the mobility of a colony of microswimmers
of arbitrarily planar geometry in terms of the properties of the individual
microswimmers and their geometric arrangement, provided the dynamical noise in
the cycle-averaged flagellar dynamics can be taken as relatively small. We
discuss how our conclusions for the simplified model may be expected to apply
in more realistic models in Section 6.

##  2 Disclike Colony Model

In order to obtain analytical relationships between the effective swimming
dynamics of the colony in terms of those of the individual cells, we begin
with a rather minimal modeling framework. A simple disclike geometric model
for the colony is described in Subsection 2.1, followed in Subsection 2.2 by
the stochastic dynamical model for the propulsion induced by each cell’s
flagellum. The resulting dynamical model for the colony is presented in
Subsection 2.3. Suitable values for the various biophysical model parameters
are discussed in Subsection 2.4, which will inform some simplifications in the
calculations to follow.

###  2.1 Colony Geometry

Our primary geometric model for a colony in a rosette shape is a two-
dimensional disc, composed of a fixed number $N$ of swimming cells. We in
particular neglect dynamical changes in cell number due to cell division or
other morphology changes, as these take place on a time scale of hours
(Fairclough et al., 2010), much longer than the basic dynamical time scales of
seconds to minutes of interest. The two-dimensionality of our model is chosen
simply to avoid the complications of three-dimensional rotational dynamics.
More precisely, we envision the colony as a thin disc in a three-dimensional
fluid, whose translational and rotational dynamics are constrained to remain
in the plane defined by the disc. The disclike approximation may actually be
reasonable for small _Salpingoeca rosetta_ colonies of 4-7 cells; larger
colonies exhibit a more three-dimensional structure (Larson et al., 2020). But
even with a truly flat morphology, we impose a simplification by constraining
the colony to only move within the plane defined by its geometry.

In our models, the cells are assumed to be of equal size and symmetrically
situated about the disc, with each cell being exposed to the ambient fluid
over an arc length $l$ of the disc, as well as along the flat transverse faces
of the thin disc. In particular, this implies the disc representing the colony
has a circumference $Nl$ and thus a radius $a=\frac{Nl}{2\pi}$. This
proportionality between radius and cell number is consistent with a two-
dimensional version of the argument in Solari et al. (2006). In reality, cells
in _S. rosettas_ colonies do exhibit some size and shape disparities (Naumann
and Burkhardt, 2019), possibly influenced by colony size (Laundon et al.,
2019; Larson et al., 2020), but the available data at this point seems
insufficient for a meaningful mathematical model of this cellular variability.
We in particular treat the physical colony in an entirely smooth and symmetric
fashion, neglecting the natural protrusion of each cell with respect to the
nominal disc shape of colony, not to mention the microvillar collar (Nguyen et
al., 2019). Asymmetries in the flagellar properties and placement will however
be represented in our model, as these appear to be more fundamental to the
determination of the colony dynamics. We view these simplifying geometric
contrivances, illustrated in Figure 1 as still capable of preserving the
substance of the interactions between the cellular components in determining
the colonial swimming properties.

The flagellum for each cell $i=1,\ldots,N$ is assumed to be within the plane
of the colonial disc, with its base attached to the associated arc at a
displacement $S_{i}\sim U(-\delta_{S},\delta_{S})$ independently and uniformly
randomly distributed within an arclength $\delta_{S}\leq l/2$ of the center of
the arc. Each flagellum $i=1,\ldots,N$ is moreover allowed demographic
stochasticity in terms of its mean propulsion force $F_{i}^{(0)}$ and its
“relaxed” angle $\Theta_{i}^{(0)}$ which it makes with the normal at its basal
point of attachment to the cell. These are each taken to be independent random
variables with $\Theta_{i}^{(0)}\sim N(0,\delta_{\Theta}^{2})$, Gaussian
distributed with mean $0$ and standard deviation $\delta_{\Theta}$. The
flagellar force $F_{i}^{(0)}$ only appears linearly or quadratically in the
transport expressions, so we will only characterize it in terms of its first
four moments or cumulants.

Figure 1: Rendition of a colony of ten cells (left) and a model schematic of
eight cells (right). The cell has the shape of a disc of radius $a$, with each
cell represented by an equal arc of length $l$. The attachment point of the
flagellum associated to the $i$th cell is displaced by an arclength $S_{i}$
from the center of its arc. The direction of the flagellar force is
dynamically directed along a direction $\Theta_{i}(t)$ with respect to the
normal. The colony’s dynamical orientation is denoted by an angle
$\Theta^{(c)}(t)$ with respect to a fixed reference direction which we take to
be along the rightward horizontal. As usual, all quantities associated to
rotational measures are taken to be positive in the counterclockwise
direction.

###  2.2 Cell Propulsive Forcing Model

Experimental and theoretical studies have characterized the approximately
periodic dynamics of flagellar beating, including the coupling between
flagellar dynamics in various cells and multicellular organisms (Brumley et
al., 2012; Guirao and Joanny, 2007). For at least choanoflagellates, however,
the flagella appear to be sufficiently distant and irregularly positioned to
make correlations between their dynamics unobservable (Kirkegaard et al.,
2016a). In our model, we will therefore neglect correlations between the
flagellar beating dynamics. The detailed representation of the periodic
beating of each flagellum $i=1,\ldots,N$ does not then appear to be necessary,
and we represent the flagellar dynamics instead by a locally cycle-averaged
force $F_{i}(t)$ and cycle-averaged orientation $\Theta_{i}(t)$ of the
flagellum with respect to the normal at its basal attachment to the cell. We
treat these variables as continuous in time; they could be thought of as
windowed time averages over a period width centered at the current time $t$.
Since flagellar beating is known to have effectively random disturbances due
to the noise in the dynein activity driving them as well as environmental
effects, we model the dynamics of the cycle-averaged flagellar force magnitude
and orientation for cell $i$ via autoregressive stochastic processes with
relaxed states $F_{i}^{(0)}$ and $\Theta_{i}^{(0)}$, dynamical variances
$\sigma_{F}^{2}$ and $\sigma_{\Theta}^{2}$, and damping parameters (or inverse
correlation times) $\gamma_{F}$ and $\gamma_{\Theta}$ respectively. The SDEs
for these processes are given as:

$\displaystyle dF_{i}(t)=$
$\displaystyle-\gamma_{F}(F_{i}(t)-F_{i}^{(0)})dt+\sqrt{2\sigma_{F}^{2}\gamma_{F}}\,\mathrm{d}W_{i}^{F}(t)$
(1) $\displaystyle d\Theta_{i}(t)=$
$\displaystyle-\gamma_{\Theta}(\Theta_{i}(t)-\Theta_{i}^{(0)})dt+\sqrt{2\sigma_{\Theta}^{2}\gamma_{\Theta}}\,\mathrm{d}W_{i}^{\Theta}(t)$
(2)

Here $\\{W_{i}^{F}\\}_{i=1}^{N}$ and $\\{W_{i}^{\Theta}\\}_{i=1}^{N}$ stand
for independent standard Wiener processes. Note that we do not associate the
noise of the flagellar dynamics to temperature, as the noise from active
driving processes is apparently dominant. Our model does not account for a
possible jumpy change in the swimming direction of a cell, as seen for fast
swimming _S. rosetta_ cells in Sparacino et al. (2020), but it is slow
swimmers that form _S. rosetta_ colonies (Dayel et al., 2011). In our
analytical and numerical computations, we take the initial conditions for the
flagella to be independently distributed according to their stationary
distributions, so $F_{i}(0)$ and $\Theta_{i}(0)$, are normally distributed
with mean $F_{i}^{(0)}$ and standard deviation $\sigma_{F}$ (respectively,
$\Theta_{i}^{(0)}$ and $\sigma_{\Theta}$).

###  2.3 Colony Dynamics Model

We assume the colony is sufficiently isolated that we can neglect external
interaction or hydrodynamic forces. The Reynolds number for colonial
microswimmers is approximately $10^{-4}-10^{-3}$ so we may use low Reynolds
number hydrodynamics based on the linear Stokes equation to represent the
hydrodynamic drag on the colony motion. Using the drag formulas from Kim and
Karrila (1991, Ch. 3) for a spheroid in the oblate limit, we take for,
respectively, the translational and rotational drag coefficient:

$\displaystyle\gamma_{\mathrm{t}}=\frac{32}{3}a\eta,\hskip
17.07182pt\gamma_{\mathrm{r}}=\frac{32}{3}a^{3}\eta$ (3)

where $\eta$ is the dynamic viscosity of the fluid.

Balancing the active and drag forces and torque, and incorporating thermal
noise, we have the following overdamped Langevin equations for the center of
mass ($\mathbf{X}^{(\mathrm{c})}$) and angular orientation ($\Theta^{(c)}$):

$\displaystyle
d\mathbf{X}^{(\mathrm{c})}(t)=\frac{\mathbf{F}^{(c)}(t)}{\gamma_{\mathrm{t}}}dt+\sqrt{2D_{\mathrm{t}}}d\mathbf{W}^{x,c}(t)$
(4) $\displaystyle
d\Theta^{(c)}(t)=\frac{T^{(c)}(t)}{\gamma_{\mathrm{r}}}dt+\sqrt{2D_{\mathrm{r}}}dW^{\Theta,c}(t).$
(5)

Here $D_{\mathrm{t}}=k_{B}T/\gamma_{\mathrm{t}}$ and
$D_{\mathrm{r}}=k_{B}T/\gamma_{\mathrm{r}}$ are, respectively, the thermally
induced translational and rotational diffusivity, $k_{B}$ is the Boltzmann
constant and T is temperature, and $\mathbf{W}^{X,c}$ and
$\mathbf{W}^{\Theta,c}$ are independent Wiener processes representing thermal
noise. The aggregate active force on the colony is obtained by adding up the
vectorial contributions from each flagellum:

$\displaystyle\mathbf{F}^{(c)}=\sum_{i=1}^{N}-F_{i}(t)\begin{pmatrix}\cos(2\pi[\frac{i-\frac{1}{2}+\frac{S_{i}}{l}}{N}]+\Theta_{i}(t)+\Theta^{(c)}(t))\\\
\sin(2\pi[\frac{i-\frac{1}{2}+\frac{S_{i}}{l}}{N}]+\Theta_{i}(t)+\Theta^{(c)}(t))\end{pmatrix}.$
(6)

The aggregate active torque is similarly obtained by adding up the torque
contributions from each flagellum:

$\displaystyle T^{(c)}=-\sum_{i=1}^{N}F_{i}(t)a\sin\Theta_{i}(t)$ (7)

To obtain the representation (6) for the total force on the colony, we have
defined the colony orientation angle $\Theta^{(\mathrm{c})}(t)$ to be the
angle with respect to direction of the horizontal made by the radius from the
center of the colony to the beginning of the arc defining cell $1$, with all
angular orientations taken in the usual counterclockwise sense (see Figure 1).
Then the center of the arc associated to cell $i$ will, at time $t$, be
situated at angular position
$\Theta^{(\mathrm{c})}(t)+\frac{2\pi}{N}\left(i-\frac{1}{2}\right)$, and the
associated flagellum base is at angular position
$\Theta^{(\mathrm{c})}(t)+\frac{2\pi}{N}\left(i-\frac{1}{2}+\frac{S_{i}}{l}\right)$.
The force applied by flagellum $i$ will be directed in the opposite direction
(into the cell body), after correction for the current cycle-averaged
displacement $\Theta_{i}(t)$ of the flagellar orientation from the normal
direction to the cell.

Because the environment in our model is spatially homogenous and isotropic,
without loss, we will take the initial configuration of the colony to satisfy
$\mathbf{X}^{(\mathrm{c})}(t=0)=\bm{0}$ and $\Theta^{(\mathrm{c})}(t=0)=0$.

###  2.4 Parameter Estimation

Estimates for the various biophysical parameters in our model are presented in
Table 1. We next explain how they were obtained from the experimental
literature.

Table 1: Estimates of model parameters

Parameter Estimated value Number of cells (N) 2-10 Arclength of cell exposure
($l$) $2\pi\,$\mathrm{\SIUnitSymbolMicro m}$$ Flagellar force ($F^{(0)}_{i}$)
$2\pm 1\,$\mathrm{pN}$$ Relaxed flagellar orientation with respect to normal
($\Theta^{(0)}_{i}$) $\pm 0.02$ Damping parameter of force orientation
($\gamma_{\Theta}$) $10\,$\sec$^{-1}$ Damping parameter of force magnitude
($\gamma_{F}$) $10\,$\sec$^{-1}$ Dynamical variance of force orientation
($\sigma_{\Theta}^{2}$) $0.002$ Dynamical variance of force magnitude
($\sigma_{F}^{2}$) $0.1\,$\mathrm{pN}$^{2}$ Temperature(T) 300 K Dynamic
viscosity ($\eta$) $0.01\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}\sec$

####  2.4.1 Estimation of Geometric Parameters

Images in Dayel et al. (2011); Kirkegaard et al. (2016a)) indicate the
individual cells in an _S. Rosetta_ colony have a diameter on the order of
$5\text{\,}\mathrm{\SIUnitSymbolMicro m}$, and Fairclough et al. (2013) report
intercellular bridges between cells have lengths of about
$0.15\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The experimental observations in
Kirkegaard et al. (2016a), which are the target of our modeling efforts,
report colonial projected areas of
$20-250$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}$$. If we assume the colony
cross sections to be approximately circular, this would correspond to
circumferences of $10-60$\mathrm{\SIUnitSymbolMicro m}$$. This would
correspond to about $N\sim 2-10$ cells under our two-dimensional disc model if
we partition the circumference of the colony into segments for each cell and
their intercellular bridges to their neighbors, and approximate the arc length
occupied by a cell by its diameter. Larson et al. (2020) find in fact _S.
rosetta_ to proceed through an approximately two-dimensional phase of growth
of 4-7 cells to a three-dimensional phase of growth from 8-12 cells. Kumler et
al. (2020) report $2-13$ cells in _Salpingoeca helianthica_ colonies. We note
that Larson et al. (2020) report shrinking cell sizes as the colony increases,
but the reported data does not seem to permit a clean scaling representation.
One could of course allow the exposed cell arclength $l$ to decrease as
desired with the number of cells $N$; our analysis is only affected in terms
of the dependence of the asymptotic errors on the colony size. We take
$l=2\pi\,$\mathrm{\SIUnitSymbolMicro m}$$ for modeling convenience since it is
roughly consistent with observations.

####  2.4.2 Estimation of Flagellar Force Parameters

The typical flagellar force magnitude $F_{i}^{(0)}$ is taken from Roper et al.
(2013). To estimate the dynamical variance of the flagellar force, we assume
the force is proportional to the frequency (Guirao and Joanny, 2007) and
consider the periodogram estimates of the frequency spectrum (Figure 2) using
the time series from Figures 2a and 2b in Kirkegaard et al. (2016a). Gaussian
fits to each resulting frequency distribution yielded estimates for the ratio
of the frequency spread (standard deviation) to the mean, which we took as an
estimate for the ratio of the dynamical standard deviation of the force to its
mean. The peak frequency of the two flagella differ by about a factor of $2$,
suggesting a rather broad demographic distribution of the flagellar force,
consistent with the spread of flagellar frequencies reported in Kirkegaard et
al. (2016a, Fig. 2e). By contrast, the polydispersivity of beating frequencies
in _Volvox carteri_ , whose partial metachronal synchronization is out of the
scope of our study, is observed to be only about 10%, or 3% when controlled
for spatial location (Brumley et al., 2012).

Figure 2: The spectral distributions of the time series for the orientation of
two flagellaa in Kirkegaard et al. (2016a, Fig. 2a. Fig. 2b)

For estimating the cycle-averaged statistics of the flagellar orientation, we
applied a moving average filter over a window equal to one period as estimated
by the frequency peak to the sample time series of the orientation of two
flagella in Figures 2a and 2b of Kirkegaard et al. (2016a), obtaining the
plots in Figure 3. The average values were both found to be $0.02$, so we used
that value as a typical magnitude (with either sign) of the relaxed or
dynamical mean $\Theta_{i}^{(0)}$ of the flagellar orientation with respect to
the normal. The dynamical variance $\sigma_{\Theta}^{2}$ of the flagellar
orientation was then estimated as the variance of these cycle-averaged time
courses of the flagellar orientation.

Figure 3: The cycle-averaged smoothing of the time series for the flagellar
orientation $\Theta$ reported in Figures 2a (blue) and 2b (red) in (Kirkegaard
et al. (2016a)

The damping parameter $\gamma_{\Theta}$ for the cycle-averaged force
orientation is easily inferred from the decay times of the orientation
correlation presented in Kirkegaard et al. (2016a, Fig. 2f). This decay time
was seen to have a broad variation across flagella, but we neglect the
demographic stochasticity of this parameter in our model. We make the
plausible choice $\gamma_{F}=\gamma_{\Theta}$ since we are not aware of time
series data on the flagellar force.

##  3 Statistics of Mobility

We present here a summary of the effective statistics of motion of the
disclike colony model described in Section 2. In Subsection 3.1 we report the
effective rotational drift and diffusivity of the colony, while in Subsection
3.2 we report the mean-square speed and translational diffusivity of the
colony. Our model has two sorts of randomness: 1) dynamical stochasticity
arising from the noisy driving terms on the flagellar force and orientation in
Subsection 2.2, and 2) demographic stochasticity, meaning random variations
between cells of the static features $F_{i}^{(0)}$, $\Theta_{i}^{(0)}$, and
$S_{i}$ of their associated flagella as defined in Subsection 2.1. For each
mobility descriptor, we first report the formulas for a particular colony as a
function of its peculiar values of the static flagellar features, but averaged
over the dynamical noise. We use $\langle\cdot\rangle$ to denote such averages
over the dynamical stochasticity, with a fixed realization of the static
flagellar properties. After the mobility quantifier formula for an individual
colony is presented, we then express its average and variance with respect to
the demographic distribution of flagellar properties described in Subsection
2.1. In this summary section, we reserve the notation $\mathbb{E}$ and Var for
this demographic average and variance. The demographic moments of the relaxed
flagellar force will be denoted
$\mu^{(0)}_{F,k}\equiv\mathbb{E}\left[(F^{(0)})^{k}\right]$, and their
corresponding cumulants by $\kappa^{(0)}_{F,k}$.

The derivation of these results will follow in Section 4. We provide some
visualization of the formulas and comparison with Monte Carlo simulations; the
method for conducting the simulations is discussed in Appendix A. We assume
the flagella are initialized with their statistically stationary distribution
for their force and orientation. If they are not, the formulas to follow only
apply at times $t\gg\gamma_{F}^{-1},\gamma_{\Theta}^{-1}$ over which the
flagella relax to their stationary distribution. The diffusivity formulas are
inherently long-time so unaffected by the initial conditions.

###  3.1 Rotational Mobility Statistics

####  3.1.1 Rotational Drift

:

$\displaystyle\Omega^{\ast}_{\mathrm{r}}\equiv\frac{\mathrm{d}\langle\Theta^{(\mathrm{c})}(t)\rangle}{\mathrm{d}t}=-\frac{a}{\gamma_{\mathrm{r}}}e^{-\frac{1}{2}\sigma_{\Theta}^{2}}\sum_{i=1}^{N}F_{i}^{(0)}\sin(\Theta_{i}^{(0)})$
(8)

The rotational drift is driven by the natural relaxed torque diminished by the
factor $\mathrm{e}^{-\frac{1}{2}\sigma_{\Theta}^{2}}$ arising from the
dynamical variation of the angular orientation of the flagella of the cells,
while being resisted by rotational drag due to the fluid environment.

Because we assume $\Theta_{i}^{(0)}$ is symmetrically distributed about $0$,
the demographic average of the rotational rate of the colony is zero.
Individual colonies will break the symmetry and have some rotational drift. We
can characterize the typical magnitude of the rotational drift by its
demographic standard deviation:

$\sigma_{\Omega^{\ast}_{\mathrm{r}}}=\sqrt{\text{Var}[\Omega^{\ast}_{\mathrm{r}}]}=\frac{a}{\gamma_{\mathrm{r}}}e^{-\frac{1}{2}\sigma_{\Theta}^{2}}\sqrt{N\mu^{(0)}_{F,2}\left(\frac{1-\mathrm{e}^{-2\delta_{\Theta}^{2}}}{2}\right)}$
(9)

We see that the magnitude of the rotational drift of a colony increases
monotonically with an increase in the demographic variability
$\delta_{\Theta}$ of the flagellar basal attachment angles, with linear
proportionality at small demographic variations. A comparison of this
demographic spread of rotational drift with Monte Carlo simulations at various
dynamical variances $\sigma_{\Theta}^{2}$ is presented in Figure 4. As
illustrated again by comparison with Monte Carlo simulations in Figure 5, the
magnitude of the rotational drift will tend to decrease with colony size as
$N^{-3/2}$ since the rotational drag coefficient $\gamma_{\mathrm{r}}$ scales
cubically with the number of cells $N$ while the radius $a$ scales linearly.

Figure 4:  Comparison of theoretical (blue) and Monte Carlo simulation (red)
results of rotational drift over a long time ($t=10000$\mathrm{s}$$) for a
colony of $N=10$ cells. The analytical result has zero mean with demographic
standard deviation from (9) indicated with error bars. Simulation results are
shown separately for 1 simulation of each of 10 colony samples at each value
of $\sigma_{\Theta}^{2}$, with the flagellar displacements modeled as
$S_{i}\sim U(\frac{-l}{2},\frac{l}{2})$. The rest of the parameters are
specified in Appendix A and Table 1. Figure 5:  Comparison of theoretical
(blue) and Monte Carlo simulation (red) results of rotational drift over a
long time ($t=10000$\mathrm{s}$$) for colonies of various sizes. The
analytical result has zero mean with demographic standard deviation from Eq.
(9) indicated with error bars. Simulation results are shown separately for 1
simulation of each of 10 colony samples at each colony size, with the
flagellar displacements modeled as $S_{i}\sim U(\frac{-l}{2},\frac{l}{2})$.
The rest of the parameters are specified in Appendix A and Table 1.

####  3.1.2 Rotational Diffusion

$\displaystyle D_{\mathrm{r}}^{\ast}$
$\displaystyle\equiv\lim\limits_{t\rightarrow\infty}\frac{(\Theta^{(c)}-\Omega^{\ast}_{\mathrm{r}}t)^{2}}{2t}=D_{\mathrm{r}}$
(10)
$\displaystyle\qquad\qquad+\frac{1}{2\gamma_{\mathrm{r}}^{2}}\Big{(}\sum_{i=1}^{N}\big{[}{(aF_{i}^{(0)})^{2}}e^{-\sigma_{\Theta}^{2}}\gamma_{\Theta}^{-1}(-\mathrm{Ein}(-\sigma_{\Theta}^{2})+\cos(2\Theta_{i}^{(0)})\mathrm{Ein}(\sigma_{\Theta}^{2}))$
$\displaystyle\qquad\qquad+a^{2}\sigma_{F}^{2}\gamma_{F}^{-1}\big{(}{}_{1}F_{1}[1,\frac{\gamma_{F}}{\gamma_{\theta}}+1,-\sigma_{\Theta}^{2}]-\cos(2\Theta_{i}^{(0)})e^{-2\sigma_{\Theta}^{2}}{}_{1}F_{1}[1,\frac{\gamma_{F}}{\gamma_{\theta}}+1,\sigma_{\Theta}^{2}]\big{)}\big{]}\Big{)}$

$\mathrm{Ein}(z)\equiv\int_{0}^{z}(1-\mathrm{e}^{-z^{\prime}})/z^{\prime}\,\mathrm{d}z^{\prime}$
is the complementary exponential integral and ${}_{1}F_{1}$ is Kummer’s
confluent hypergeometric function (Abramowitz and Stegun, 1964).

Taking the mean and variance of this rotational diffusion formula over the
demographic distribution of flagella, we obtain:

$\displaystyle\mathbb{E}\left[D_{\mathrm{r}}^{\ast}\right]$
$\displaystyle=D_{\mathrm{r}}+\frac{Na^{2}}{2\gamma_{\mathrm{r}}^{2}}\left(\frac{\mathrm{e}^{-\sigma_{\Theta}^{2}}}{\gamma_{\Theta}}\mu^{(0)}_{F,2}\Big{(}-\mathrm{Ein}(-\sigma_{\Theta}^{2})+\mathrm{e}^{-2\delta_{\Theta}^{2}}\mathrm{Ein}(\sigma_{\Theta}^{2})\Big{)}\right.$
(11)
$\displaystyle\left.\qquad\qquad+\frac{\sigma_{F}^{2}}{\gamma_{F}}({}_{1}F_{1}[1,\frac{\gamma_{F}}{\gamma_{\theta}}+1,-\sigma_{\Theta}^{2}]-\mathrm{e}^{-2\delta_{\Theta}^{2}-2\sigma_{\Theta}^{2}}{}_{1}F_{1}[1,\frac{\gamma_{F}}{\gamma_{\theta}}+1,\sigma_{\Theta}^{2}])\Big{)}\right),$
$\displaystyle\text{Var}\left[D_{\mathrm{r}}^{\ast}\right]$
$\displaystyle=\frac{Na^{4}}{4\gamma_{\mathrm{r}}^{4}}$
$\displaystyle\times\left(\frac{\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{\gamma_{\theta}^{2}}\mu^{(0)}_{F,4}\Big{(}\mathrm{Ein}^{2}(-\sigma_{\Theta}^{2})-2\operatorname{Ein}(-\sigma_{\Theta}^{2})\operatorname{Ein}(\sigma_{\Theta}^{2})\mathrm{e}^{-2\delta_{\Theta}^{2}}+\frac{1+\mathrm{e}^{-8\delta_{\Theta}^{2}}}{2}\operatorname{Ein}^{2}(\sigma_{\Theta}^{2})\Big{)}\right.$
$\displaystyle\qquad\qquad-\frac{\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{\gamma_{\theta}^{2}}(\mu^{(0)}_{F,2})^{2}\Big{(}-\mathrm{Ein}(-\sigma_{\Theta}^{2})+\mathrm{e}^{-2\delta_{\Theta}^{2}}\operatorname{Ein}(\sigma_{\Theta}^{2})\Big{)}^{2}$
$\displaystyle\qquad\qquad\left.+\frac{\sigma_{F}^{4}\mathrm{e}^{-4\sigma_{\Theta}^{2}}}{2\gamma_{F}^{2}}(1-\mathrm{e}^{-4\delta_{\Theta}^{2}})^{2}({}_{1}F_{1}^{2}[1,\frac{\gamma_{F}}{\gamma_{\theta}}+1,\sigma_{\Theta}^{2}])\right),$

These expressions can be expressed more simply for the relevant regime of
small flagellar angular displacements (dropping the fourth order terms in the
asymptotic expansion):

$\displaystyle D_{\mathrm{r}}^{\ast}$
$\displaystyle=D_{\mathrm{r}}+\frac{a^{2}}{\gamma_{\mathrm{r}}^{2}}\left(\sum_{i=1}^{N}\left[(F_{i}^{(0)})^{2}\frac{\sigma_{\Theta}^{2}}{\gamma_{\Theta}}+\sigma_{F}^{2}\left(\frac{\left(\Theta_{i}^{(0)}\right)^{2}}{\gamma_{F}}+\frac{\sigma_{\Theta}^{2}}{\gamma_{F}+\gamma_{\Theta}}\right)\right]\right)$
(12)

with demographic mean and standard deviation:

$\displaystyle\mathbb{E}[D_{\mathrm{r}}^{\ast}]$
$\displaystyle=D_{\mathrm{r}}+\frac{Na^{2}}{\gamma_{\mathrm{r}}^{2}}\Big{(}\mathbb{E}[(F^{(0)})^{2}]\frac{\sigma_{\Theta}^{2}}{\gamma_{\Theta}}+\sigma_{F}^{2}\left(\frac{\delta_{\Theta}^{2}}{\gamma_{F}}+\frac{\sigma_{\Theta}^{2}}{\gamma_{F}+\gamma_{\Theta}}\right)\Big{)},$
(13) $\displaystyle\text{Var}[D_{\mathrm{r}}^{\ast}]$
$\displaystyle=\frac{Na^{4}}{\gamma_{\mathrm{r}}^{4}}\Big{(}\text{Var}[(F^{(0)})^{2}]\frac{\sigma_{\Theta}^{4}}{\gamma_{\Theta}^{2}}+\frac{2\sigma_{F}^{4}\delta_{\Theta}^{4}}{\gamma_{F}^{2}}\big{)}\Big{)}$

These formulas are illustrated as a function of dynamical variance
$\sigma_{\Theta}^{2}$ in comparison with Monte Carlo simulations in Figure 6.
Using the biophysical parameter values from Table 1, we verify that the
contribution from rotational diffusion due to the active flagellar driving is
significantly stronger than the thermal one. This can be seen in the
comparison of the formulas for various colony sizes against Monte Carlo
simulations in Figure 7. Since the rotational drag is proportional to the cube
of the colony size the long time rotational diffusion enhancement will
decrease with colony size as $\frac{1}{N^{3}}$, and its demographic standard
deviation as $\frac{1}{N^{7/2}}$. As can be expected on physical grounds, the
active contribution to the rotational diffusion relies on dynamical torque
fluctuations from the flagella induced both by orientational fluctuations and
force fluctuations along a quenched angle deflected from the normal by the
cell’s orientation within the colony.

On comparing Figures 4 and 6, we see that the rotational drift is typically
small compared to the rotational diffusivity over most of the potential range
of dynamical angular variances plotted. But when the flagellar orientations
are small deviations from the normal, as it appears to be for
choanoflagellates (Table 1), the typical rotational drift and diffusivity for
a colony are comparable (Figures 5 and 7). Indeed, taking the relaxed
orientations as symmetrically distributed about the normal with small
demographic standard deviation $\delta_{\Theta}\ll 1$ and neglecting the
thermal contribution $D_{\mathrm{r}}$ to the rotational diffusivity, the mean
rotation drift rate (8) should scale as
$\Omega^{\ast}_{\mathrm{r}}\sim\sqrt{N}\bar{F}a\delta_{\Theta}/\gamma_{\mathrm{r}}$,
while the effective rotational diffusivity should scale as
$D_{\mathrm{r}}^{\ast}\sim
Na^{2}/\gamma_{\mathrm{r}}^{2}[(\bar{F}^{2}\sigma_{\Theta}^{2}/\gamma_{\Theta}+\sigma_{F}^{2}(\delta_{\Theta}^{2}/\gamma_{F}+\sigma_{\Theta}^{2}/(\gamma_{F}+\gamma_{\Theta})]$,
where $\bar{F}$ denotes a typical magnitude of the force $F_{i}^{(0)}$ of a
flagellum. The ratio of the effective rotational drift and effective
rotational diffusivity is thus

$\frac{\Omega^{\ast}_{\mathrm{r}}}{D_{\mathrm{r}}^{\ast}}\sim
N^{-1/2}\left(\frac{\gamma_{\mathrm{r}}}{\bar{F}a}\right)\frac{\delta_{\Theta}}{\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1}+(\sigma_{F}/\bar{F})^{2}(\delta_{\Theta}^{2}\gamma_{F}^{-1}+\sigma_{\Theta}^{2}(\gamma_{F}+\gamma_{\Theta})^{-1})}$
(14)

Taking $\sigma_{\Theta}^{2}$ and $\delta_{\Theta}^{2}$ as comparable, and
$\sigma_{F}/\bar{F}$ as small, consistent with Table 1, we can simplify this
to:

$\frac{\Omega^{\ast}_{\mathrm{r}}}{D_{\mathrm{r}}^{\ast}}\sim
N^{-1/2}\frac{\gamma_{\mathrm{r}}\gamma_{\Theta}}{\bar{F}a}\frac{\delta_{\Theta}}{\sigma_{\Theta}^{2}}.$
(15)

We notice the ratio of the physical parameters can be viewed as the ratio of
the relaxation rate $\gamma_{\Theta}$ of the flagellar orientation and the
nominal angular frequency scale $\bar{F}a/\gamma_{\mathrm{r}}$ a flagellum
would impart to the colony, if its force were applied normally to the colony
surface (which it is not). From the numerical parameters in Table 1:

$\frac{\gamma_{\mathrm{r}}\gamma_{\Theta}}{\bar{F}a}\sim
0.05N^{2},\qquad\frac{\delta_{\Theta}}{\sigma_{\Theta}^{2}}=10$ (16)

so $\Omega^{\ast}_{\mathrm{r}}/D_{\mathrm{r}}^{\ast}\sim 0.5N^{3/2}$. Thus,
for the biophysical parameters in Table 1, the effective rotational drift and
diffusivity are comparable for small colonies but the effective rotational
diffusivity is reduced at a faster $N^{-3}$ rate relative to the $N^{-3/2}$
dependence of the rotational drift.

The relation between the effective rotational drift and diffusivity can also
be expressed somewhat more naturally, under just the assumption that the
relaxed flagellar orientations $\Theta_{j}^{(0)}$ are symmetrically
distributed about the normal,

$\frac{D_{\mathrm{r}}^{\ast}}{\Omega^{\ast}_{\mathrm{r}}{}^{2}}\sim\frac{\sigma_{\Theta}^{2}}{\delta_{\Theta}^{2}\gamma_{\Theta}}+\frac{\sigma_{F}^{2}}{\gamma_{F}\bar{F}^{2}}\left(1+\frac{\sigma_{\Theta}^{2}}{\delta_{\Theta}^{2}}\right),$
(17)

These relative sizes of rotational drift and diffusivity have substantial
impact on the translational diffusivity as will be seen in Subsection 3.2.2.

Figure 6: Comparison of theoretical (blue) and Monte Carlo simulation (red)
results of rotational diffusivity over a long time ($t=10000$\mathrm{s}$$) for
a colony of $N=10$ cells. The analytical result is plotted as error bars
centered at the mean with width equal to one standard deviation using the
formulas from (11). The parameters and simulations are the same as in Figure
4. Figure 7:  Comparison of theoretical (blue) and Monte Carlo simulation
(red) results of rotational diffusivity over a long time
($t=10000$\mathrm{s}$$) for colonies of various sizes. The analytical result
is plotted as error bars centered at the mean with width equal to one standard
deviation using the formulas from (13), which are visually indistinguishable
from the more precise expressions in Eq. (11). The parameters and simulations
are the same as in Figure 5.

###  3.2 Translational Mobility Statistics

Because the environment is isotropic, any initial mean velocity of the colony
will dissipate on the time scale of rotational diffusion, and therefore the
colony does not exhibit sustained translational drift.

####  3.2.1 Mean-Square Speed

We therefore look at the instantaneous colony velocity induced by active
flagellar forces,
$\mathbf{V}^{(\mathrm{c})}\equiv\frac{\mathbf{F}^{(c)}}{\gamma_{\mathrm{t}}}$.
The thermal Brownian motion component to the colony velocity is not well-
defined in our model (4) which coarse-grains inertial dynamics, but on
physical grounds would have mean-square value $\frac{2k_{B}T}{m}$ where $m$ is
the mass of the colony. We will neglect this thermal component in the results
we are about to present, and verify afterward that its physical contribution
would be subdominant for the biophysical parameter values in Table 1.

An individual colony’s mean-square speed (from active flagellar-driven
contributions) is:

$\displaystyle\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle$
$\displaystyle=\gamma_{\mathrm{t}}^{-2}\left[N\sigma_{F}^{2}+e^{-\sigma_{\Theta}^{2}}\left|\sum_{j=1}^{N}F_{j}^{(0)}\mathrm{e}^{\mathrm{i}\Theta_{j}^{(0)}+\frac{2\pi\mathrm{i}}{N}\left(j+\frac{S_{j}}{l}\right)}\right|^{2}+(1-\mathrm{e}^{-\sigma_{\Theta}^{2}})\sum_{j=1}^{N}\left(F_{j}^{(0)}\right)^{2}\right]$
(18)

The squared modulus of the complex sum is simply the squared norm of the total
force on the colony when each flagellum $j$ is exerting its mean force
$F_{j}^{(0)}$ and its relaxed orientation $\Theta_{j}^{(0)}$. The low variance
$\sigma_{\Theta}^{2}$ of the angular orientation of the flagellar forcing
implies a substantial impact of the asymmetrical contributions of forcing from
the cells on the instantaneous velocity of the colony.

We remark that the mean-square displacement of the colony adjusts from an
initial ballistic phase to a second ballistic phase over time scales between
those of the flagellar dynamics and the colony dynamics:

$\displaystyle\mathbb{E}|\mathbf{X}^{(\mathrm{c})}(t+\tau)-\mathbf{X}^{(\mathrm{c})}(t)|^{2}\sim
4D_{\mathrm{t}}\tau+V^{\ast}{}^{2}\tau^{2}\text{ for
}\gamma_{F}^{-1},\gamma_{\Theta}^{-1}\ll\tau\ll\Omega^{\ast}_{\mathrm{r}}{}^{-1},D_{\mathrm{r}}^{\ast}{}^{-1},$
$\displaystyle
V^{\ast}{}^{2}\sim\gamma_{\mathrm{t}}^{-2}e^{-\sigma_{\Theta}^{2}}\left|\sum_{j=1}^{N}F_{j}^{(0)}\mathrm{e}^{\mathrm{i}\Theta_{j}^{(0)}+\frac{2\pi\mathrm{i}}{N}\left(j+\frac{S_{j}}{l}\right)}\right|^{2}=\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle-\frac{N\sigma_{F}^{2}}{\gamma_{\mathrm{t}}^{2}}-\frac{(1-\mathrm{e}^{-\sigma_{\Theta}^{2}})}{\gamma_{\mathrm{t}}^{2}}\sum_{j=1}^{N}\left(F_{j}^{(0)}\right)^{2}$
(19)

where the expression for the mean-square velocity inferred over these
intermediate time scales is approximated through second order in the dynamical
flagellar fluctuations. We note this mean-square velocity, which might be more
relevant for comparison to experiments, differs from the instantaneous mean-
square velocity (18) by terms second order in the dynamical flagellar
fluctuations. It is simply the square of the colony speed $V^{\ast}$ induced
by each flagellum $j$ applying its mean force $F_{j}^{(0)}$ about its mean
orientation $\Theta_{j}^{(0)}$ relative to the normal from the attachment
point, multiplied by $\mathrm{e}^{-\sigma_{\Theta}^{2}}$ to account for the
averaging over dynamical flagellar fluctuations. We continue our discussion in
terms of the instantaneous mean-square velocity (18); rather similar
conclusions with slight adjustments apply for the mean-square velocity
observed over longer intervals in Eq. (19).

The demographic average and standard deviation, derived in Appendix C, are:

$\displaystyle\mathbb{E}[\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle]=N\gamma_{\mathrm{t}}^{-2}\left(\sigma_{F}^{2}+\mu^{(0)}_{F,2}-\mathrm{e}^{-\delta_{\Theta}^{2}-\sigma_{\Theta}^{2}}\operatorname{sinc}(\delta_{S}/a)(\mu^{(0)}_{F,1})^{2}\right)$
(20)
$\displaystyle\text{Var}(\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle)=N^{2}\mathrm{e}^{-2\sigma_{\Theta}^{2}}\gamma_{\mathrm{t}}^{-4}\left(\mu^{(0)}_{F,2}-\mathrm{e}^{-\delta_{\Theta}^{2}}\operatorname{sinc}^{2}(\delta_{S}/a)(\mu^{(0)}_{F,1})^{2})\right)^{2}$
(21)
$\displaystyle\qquad\qquad+4N\gamma_{\mathrm{t}}^{-4}(1-\mathrm{e}^{-\sigma_{\Theta}^{2}}\mathrm{e}^{-\delta_{\Theta}^{2}}\operatorname{sinc}^{2}(\delta_{S}/a))\kappa^{(0)}_{F,3}\mu^{(0)}_{F,1}$
$\displaystyle\qquad\qquad+N\gamma_{\mathrm{t}}^{-4}\kappa^{(0)}_{F,4}+N\gamma_{\mathrm{t}}^{-4}(2-\mathrm{e}^{-2\sigma_{\Theta}^{2}}-\mathrm{e}^{-2\sigma_{\Theta}^{2}-4\delta_{\Theta}^{2}}\operatorname{sinc}^{2}(2\delta_{S}/a))(\mu^{(0)}_{F,2})^{2}$
$\displaystyle\qquad\qquad-2N\gamma_{\mathrm{t}}^{-4}(1-4\mathrm{e}^{-\sigma_{\Theta}^{2}-\delta_{\Theta}^{2}}\operatorname{sinc}^{2}(\delta_{S}/a)+3\mathrm{e}^{-2\sigma_{\Theta}^{2}-2\delta_{\Theta}^{2}}\operatorname{sinc}^{4}(\delta_{S}/a))(\mu^{(0)}_{F,1})^{4}$
$\displaystyle\qquad\qquad+4N\gamma_{\mathrm{t}}^{-4}\operatorname{sinc}^{2}(\delta_{S}/a)(-2\mathrm{e}^{-\sigma_{\Theta}^{2}-\delta_{\Theta}^{2}}+\mathrm{e}^{-2\sigma_{\Theta}^{2}-3\delta_{\Theta}^{2}}\operatorname{sinc}(2\delta_{S}/a)+\mathrm{e}^{-2\sigma_{\Theta}^{2}-\delta_{\Theta}^{2}})(\mu^{(0)}_{F,1})^{2}\mu^{(0)}_{F,2}.$

The representation of the variance fully in terms of cumulants, by the
substitution $\mu^{(0)}_{F,2}=\kappa^{(0)}_{F,2}+(\mu^{(0)}_{F,1})^{2}$ and
$\mu^{(0)}_{F,1}=\kappa^{(0)}_{F,1}$ leads to a more cumbersome expression. We
present the demographic statistics for the mean-square speed rather than the
root-mean-square speed because they have a less complicated analytic
expression. We can simplify the expressions somewhat by neglecting the
demographic and dynamical fluctuation in flagellar orientation relative to
that in the flagellar force ($\delta_{\Theta},\sigma_{\Theta}\downarrow 0$) to
obtain:

$\displaystyle\mathbb{E}[\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle]$
$\displaystyle=N\gamma_{\mathrm{t}}^{-2}\left(\sigma_{F}^{2}+\text{Var}[F^{(0)}]+\left(1-\operatorname{sinc}\left(\frac{\delta_{S}}{a}\right)\right)(\mathbb{E}[F^{(0)}])^{2}\right)$
(22) $\displaystyle\text{Var}(\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle)$
$\displaystyle=N^{2}\gamma_{\mathrm{t}}^{-4}\left[\text{Var}[F^{(0)}]+\left(1-\operatorname{sinc}^{2}\left(\frac{\delta_{S}}{a}\right)\right)(\mathbb{E}[F^{(0)}])^{2}\right]^{2}$
$\displaystyle\qquad\qquad+4N\gamma_{\mathrm{t}}^{-4}(1-\operatorname{sinc}^{2}(\delta_{S}/a))\kappa^{(0)}_{F,3}\mu^{(0)}_{F,1}$
$\displaystyle\qquad\qquad+N\gamma_{\mathrm{t}}^{-4}\kappa^{(0)}_{F,4}+N\gamma_{\mathrm{t}}^{-4}(1-\operatorname{sinc}^{2}(2\delta_{S}/a))(\mu^{(0)}_{F,2})^{2}$
$\displaystyle\qquad\qquad-2N\gamma_{\mathrm{t}}^{-4}(1-4\operatorname{sinc}^{2}(\delta_{S}/a)+3\operatorname{sinc}^{4}(\delta_{S}/a))(\mu^{(0)}_{F,1})^{4}$
$\displaystyle\qquad\qquad+4N\gamma_{\mathrm{t}}^{-4}\operatorname{sinc}^{2}(\delta_{S}/a)(-1+\operatorname{sinc}(2\delta_{S}/a))(\mu^{(0)}_{F,1})^{2}\mu^{(0)}_{F,2}$

Further neglecting the demographic fluctuation in flagellar placement, so that
only variations in flagellar force magnitude are taken into account, we take
$\delta_{S}/a,\delta_{\Theta},\sigma_{\Theta}\downarrow 0$ to obtain:

$\displaystyle\mathbb{E}[\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle]$
$\displaystyle=N\gamma_{\mathrm{t}}^{-2}\left(\sigma_{F}^{2}+\text{Var}[F^{(0)}]\right)$
(23) $\displaystyle\text{Var}(\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle)$
$\displaystyle=\gamma_{\mathrm{t}}^{-4}\left[N^{2}(\text{Var}[F^{(0)}])^{2}+N\kappa^{(0)}_{F,4}\right].$

Figure 8: Comparison of theoretical (blue) and simulated (red) results of
colony root mean square speed plotted against population size for a colony
with $S_{i}\sim U(\frac{-l}{2},\frac{l}{2})$ for each cell. The colony
population size ranges from 2 to 10 cells, with 10 independently sampled
colonies for each cell size. The theoretical values are plotted as
$\sqrt{\mathbb{E}[\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle]}$ with error
bars extending from
$\sqrt{\mathbb{E}[\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle]\pm\sqrt{\text{Var}(\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle)}}$,
using the formulas in Eq. (22). For most cell sizes, in fact the theoretical
standard deviation of the mean-square speed exceeds its mean, in which case
the error bars are terminated below at $0$ speed. The rest of the parameters
are specified in Appendix A and Table 1.

The approximation used in the theoretical formulas (22) seem to be well
satisfied for the biophysical parameter values in Table 1, and indeed are
shown in Figure 8 to agree well with Monte Carlo simulations. The even simpler
formula (23) for the statistics of the root-mean-square velocity, which
neglects variability in flagellar placement, is not accurate for small colony
sizes and underpredicts within about 20% of the more precise formula (22) for
$N\geq 6$ for our parameter choices. To check on the neglect of the thermal
contribution to the colony root-mean-square speed, we take a single cell to be
a sphere of diameter $5\text{\,}\mathrm{\SIUnitSymbolMicro m}$, the diameter
of a cell (Dayel et al., 2011; Kirkegaard et al., 2016a), and have mass
density on the order of $1$\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$.
Then we estimate the mass of a single cell as $M\approx 60$\mathrm{pg}$$, and
thus the root-mean-square thermal speed scale of a single cell as
$\sqrt{\frac{2k_{B}T}{M}}\sim$10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{/}\mathrm{s}$$.
The thermal speed of the colony will decrease with size, and we confirm from
Figure 8 that the colony speed contribution from active flagellar forces is
about an order of magnitude greater than that of thermal effects.

On comparison with the experimental data for the choanoflagellate _S. rosetta_
from Kirkegaard et al. (2016a, Fig. 4), the speeds from our theoretical model
are an order unity multiple too large, but this may be attributed to an
underrepresentation of the theoretical drag by our treatment of the swimming
cells as wedges rather than osculating discs or spheres of the observed
diameter. More importantly, our theoretical model predicts a general
$N^{-1/2}$ decrease of the root-mean-square-speed with colony size (as
$\gamma_{\mathrm{t}}\propto N$), while Kirkegaard et al. (2016a) observe a
slow increase of speed with colony size. This was attributed, in a three-
dimensional model, to the number of flagella scaling as the square of the
colony radius while the drag scales only linearly. This effect would be lost
in our two-dimensional disclike model, where the number of flagella scales
only linearly with colony radius. The observed trend of a speed increase with
colony size is still a bit mysterious because it would require all the
flagella in the three-dimensional model to coordinate their orientation to a
common direction, which would seem a bit awkward for the cells situated away
from the rear of the direction of motion. Our model treats each flagellum as
autonomous and indifferent to the direction of colony propulsion so would have
no such coordination effect. One other feature observed in Kirkegaard et al.
(2016a, Fig. 4) is the wide variability in the speed of colonies of a certain
size, with the standard deviation being about half the average swimming speed.
Figure 8 shows similar behavior in our model, with greater variability at
smaller colony sizes. For large colony sizes, our theoretical results (20) and
(21) give

$\frac{\sqrt{\text{Var}(\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle)}}{\mathbb{E}[\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle]}\sim\frac{\mathrm{e}^{-\sigma_{\Theta}^{2}}\mu^{(0)}_{F,2}}{\sigma_{F}^{2}+\mu^{(0)}_{F,2}}\text{
as }N\rightarrow\infty,$

which should be around $1$ when the demographic variability in the flagellar
force $\text{Var}[F^{(0)}]$ is substantially greater than the dynamical
variance in flagellar force $\sigma_{F}^{2}$ and the dynamical orientation
variance $\sigma_{\Theta}^{2}$ is small, as is the case for the biophysically
relevant parameter values inferred in Table 1.

Kumler et al. (2020) compare the swimming speeds of unicellular and
multicellular forms of _S. helianthica_ , another choanoflagellate with a
similar morphology to _S. rosetta_ , but presumably different biophysical
parameters. One of their observations is the swimming speed of multicellular
colonies (of up to $N=13$ cells) is comparable to that of unicellular
colonies, and that the standard deviation of their speeds are no more than
about a third of the mean. Our model does not even qualitatively explain these
features. Among various reasons for the disagreement, one could be that our
model is for an isolated swimming colony while the experiments generally
involve suspensions.

####  3.2.2 Translational Diffusivity

The translational diffusivity of the colony becomes unwieldy to calculate for
general noise parameters, so we will specialize to the case where the
dynamical fluctuations $\sigma_{\Theta}$ of flagellar orientation are small,
and the magnitude $\sigma_{F}$ of dynamical force variations is small relative
to the magnitude of the relaxed flagellar forces $F_{j}^{(0)}$. Our results
are presented perturbatively through second order in these parameters being
assumed small. Note we are allowing substantial demographic variations in
flagellar force magnitudes $F_{j}^{(0)}$ and relaxed orientations
$\Theta_{j}^{(0)}$, and are stipulating only that the _dynamical_ variations
in flagellar force magnitude and orientation are relatively small. This is
important for the generalization of our calculation to other geometries in
Section 5.

We also assume that the effective rotational drift and diffusivity of the
colony, including the active contributions from the flagella, occur slowly
relative to the flagellar fluctuation time scale:

$\displaystyle\Omega^{\ast}_{\mathrm{r}},D_{\mathrm{r}}^{\ast}\ll\gamma_{F},\gamma_{\Theta}.$
(24)

The colony force autocorrelation function will then decay on a longer time
scale than the time scale of adjustment of the individual flagellar forces.
This in particular requires the thermal rotational diffusion rate
$D_{\mathrm{r}}$ to be small compared to $\gamma_{F}$ and $\gamma_{\Theta}$,
which amounts to the well-satisfied relation $0.4N^{-3}$\sec$^{-1}\ll
10$\sec$^{-1}$ for the biophysical parameters from Table 1. Figure 7 show that
the inequality (24), when including active contributions, is still satisfied
by about a factor of 10, with a wider ratio for larger colonies, for the same
specified biophysical parameter values.

The above assumptions are fundamental to our approximate calculation, and a
further mild technical assumption allows us to avoid tedious consideration of
uninteresting cases:

  * •

The time scales of relaxation of flagellar force and orientation are
comparable:

$\gamma_{F}/\gamma_{\Theta}\sim\operatorname{ord}(1).$ (25)

Indeed, we find this condition satisfied by the parameters we inferred from
experimental observations in Table 1.

Under the assumptions described above, we obtain the following approximation
for the translational diffusivity:

$\displaystyle\lim\limits_{t\rightarrow\infty}\frac{\langle\mathbf{X}^{(\mathrm{c})}(t)\otimes\mathbf{X}^{(\mathrm{c})}(t)\rangle}{2t}$
$\displaystyle=D_{\mathrm{t}}^{\ast}\mathit{\mathsf{I}},$ $\displaystyle
D_{\mathrm{t}}^{\ast}$ $\displaystyle\equiv
D_{\mathrm{t}}+\frac{V^{\ast}{}^{2}}{2}\frac{D_{\mathrm{r}}^{\ast}}{D_{\mathrm{r}}^{\ast}{}^{2}+\Omega^{\ast}_{\mathrm{r}}{}^{2}}+\tilde{D_{\mathrm{t}}}$
(26)

where $V^{\ast}{}^{2}$ is the mean-square velocity coarse-grained over the
flagellar time scale (19), and

$\displaystyle\tilde{D_{\mathrm{t}}}\equiv\frac{a\Omega^{\ast}_{\mathrm{r}}}{2\gamma_{\mathrm{t}}^{2}\gamma_{\mathrm{r}}(\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2})}\left[\frac{2\sigma_{F}^{2}}{\gamma_{F}}\sum_{j,j^{\prime}=1}^{N}F_{j}^{(0)}\sin\Theta_{j^{\prime}}^{(0)}\cos\left(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)}+\frac{2\pi}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)\right)\right.$
$\displaystyle\left.\qquad+\frac{\sigma_{\Theta}^{2}}{\gamma_{\Theta}}\sum_{j,j^{\prime}=1}^{N}F_{j}^{(0)}F_{j^{\prime}}^{(0)}(F_{j^{\prime}}^{(0)}\cos(\Theta_{j^{\prime}}^{(0)})-F_{j}^{(0)}\cos(\Theta_{j}^{(0)}))\sin\left(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)}+\frac{2\pi}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)\right)\right]$

Here $\Omega^{\ast}_{\mathrm{r}}$ is the effective rotational drift
coefficient and $D_{\mathrm{r}}^{\ast}$ is the effective rotational diffusion
coefficient given in (8) and (10) respectively. As shown in Subsection 4.2,
the error in our calculation is essentially fourth order in the small
parameters $\sigma_{\Theta}$ and $\sigma_{F}$.

The first two terms of the expression in Eq. (26) are exactly the
translational diffusion for a colony with thermal diffusivity
$D_{\mathrm{t}}$, moving at a constant speed $V^{\ast}$ along an orientation
with constant rotational drift $\Omega^{\ast}_{\mathrm{r}}$ and constant
rotational diffusion $D_{\mathrm{r}}^{\ast}$. As we shall shortly explain, the
complicated terms in $\tilde{D_{\mathrm{t}}}$ are included for completeness
and consistency with the results presented for general geometries in
Subsection 5.2, but should typically be considerably smaller than the simpler
terms presented in Eq. (26) for disclike geometries.

Due to the nonpolynomial dependence of the translational diffusion on the
effective rotational statistics, a rigorous computation of the demographic
mean is complicated. We therefore content ourselves with an estimate for the
demographic mean by averaging separately the numerator and denominator,
neglecting demographic correlations between the colony speed and rotational
characteristics, and dropping the small contribution from
$\tilde{D_{\mathrm{t}}}$:

$\mathbb{E}[D_{\mathrm{t}}^{\ast}]\approx\mathbb{E}[V^{\ast}{}^{2}]\frac{\mathbb{E}[D_{\mathrm{r}}^{\ast}]}{2(\text{Var}[D_{\mathrm{r}}^{\ast}]+(\mathbb{E}[D_{\mathrm{r}}^{\ast}])^{2}+\text{Var}[\Omega^{\ast}_{\mathrm{r}}])}+D_{\mathrm{t}}.$
(27)

All component demographic averages have been reported above, other than
$\mathbb{E}[V^{\ast}{}^{2}]$ which is clearly representable as a linear
combination of $\mathbb{E}[\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle]$ and
$\mathbb{E}[(F_{j}^{(0)})^{2}]$. We do not report a demographic standard
deviation because it is cumbersome, and in any case is simply a compounding of
the demographic fluctuations from the mean-square speed (21) and the the rate
of velocity decorrelation induced by the effective rotational drift (9) and
rotational diffusion (13) of the colony.

Figure 9:  Comparison of theoretical (blue) and simulated (red) results of the
effective translational diffusivity plotted against colony size. For each
colony size, we generate 10 colonies with flagellar displacements modeled as
$S_{i}\sim U(\frac{-l}{2},\frac{l}{2})$ for each cell. For each of these
colonies, the theoretical diffusivity computed from (26) using the peculiar
flagellar parameters of that colony is plotted as a blue circle, while the
demographic mean formula (27) is plotted as a blue dashed line. The
contribution from terms in $\tilde{D_{\mathrm{t}}}$ are neglected in the
theoretical calculation. The red circles, offset horizontally for clarity,
represent the translational diffusivity estimated from simulations. The rest
of the parameters are specified using Section A and Table 1.

As noted in Subsection 3.2.1, the mean-square colony speed scales inversely
with colony size $N$.

The second factor in the expressions (26) and (27) is essentially a rotational
decorrelation time scale, with this rotational time being determined by
rotational diffusivity when it dominates rotational drift, and by a shorter
time scale when rotational drift is stronger. The reason for the reduction in
the translational diffusivity by rotational drift is that the steady
rotational drift causes (in absence of rotational diffusion) a perfect
cancellation of sustained translational motion. From Figures 5 and 7, we see
the typical rotational drift is somewhat larger than the rotational
diffusivity using the biophysical parameters from Table 1, noting moreover
from the discussion in Subsection 3.1.1 and Subsection 3.1.2 that rotational
drift decays more slowly ($\sim N^{-3/2}$) with colony size than rotational
diffusion ($\sim N^{-3}$). Thus, when the rotational drift dominates in Eq.
(27), the rotational decorrelation time scale (second factor) is not
substantially varying with $N$ and the demographically averaged translational
diffusivity should decrease inversely with colony size $N$.

We show in Figure 9 a plot of the theoretical formulas for translational
diffusivity compared against simulations. The actual magnitude is presumably
an order of magnitude too large because of our disc model’s underestimate of
the translational drag. Focusing on the qualitative features, we see the
theoretical predictions based on the properties of specific colonies (blue
circles) agree roughly with the simulated results (red circles), while the
approximation (27) for the demographic average of the translational
diffusivity (dashed line) does not seem so accurate. In particular, the
approximation (27) for the demographic average predicts an order of magnitude
decrease in the theoretical translational diffusivity as the colony size
increases from $N=2$ to $10$, while neither the simulations nor the
distribution of theoretical predictions for specific colonies show evidence of
a decrease with colony size. We also observe the demographic variability over
two orders of magnitude for the theoretical translational diffusivity for a
given colony size, while the simulations show a somewhat more pronounced
variability, with a few colonies tending to show considerably higher
translational diffusivity than theoretically predicted. The strong
variabilities in the theoretical translational diffusivity, and therefore the
inadequacy of the averaging approximation in (27), can be traced to the
variability in the rotational drift contribution
$\Omega^{\ast}_{\mathrm{r}}{}^{2}$, which becomes typically more dominant in
the denominator as the colony size increases. Recall from Subsection 3.1.1
that the rotational drift has mean zero due to symmetry but a substantial
demographic variation. For colonies with a low rotational drift, the
rotational decorrelation time is anomalously increased and the theoretical
translational diffusivity well above the demographic mean.

Referring to Eq. (14), we see that the ratio of the more complex expression
proportional to $\Omega^{\ast}_{\mathrm{r}}$ to the expression proportional to
$D_{\mathrm{r}}^{\ast}$ in Eq. (26) is no larger than

$\displaystyle\frac{2\tilde{D_{\mathrm{t}}}(D_{\mathrm{r}}^{\ast}{}^{2}+\Omega^{\ast}_{\mathrm{r}}{}^{2})}{D_{\mathrm{r}}^{\ast}V^{\ast}{}^{2}}$
$\displaystyle\sim
O\left(\frac{\Omega^{\ast}_{\mathrm{r}}}{D_{\mathrm{r}}^{\ast}}\frac{a\bar{F}}{\gamma_{\mathrm{r}}}\left[\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1}+(\sigma_{F}^{2}/\bar{F}^{2})\gamma_{F}^{-1}\delta_{\Theta}\right]\right)\sim
O\left(N^{-1/2}\delta_{\Theta}\frac{\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1}+(\sigma_{F}^{2}/\bar{F}^{2})\gamma_{F}^{-1}\delta_{\Theta}}{\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1}+(\sigma_{F}^{2}/\bar{F}^{2})\gamma_{F}^{-1}\delta_{\Theta}^{2}}\right)$
(28) $\displaystyle\sim\begin{cases}O(N^{-1/2}\delta_{\Theta}^{2})&\text{ if
}(\sigma_{F}^{2}/\bar{F}^{2})\gamma_{F}^{-1}\delta_{\Theta}\lesssim\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1},\\\
O\left(N^{-1/2}\delta_{\Theta}^{2}\frac{(\sigma_{F}^{2}/\bar{F}^{2})\gamma_{F}^{-1}}{\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1}}\right)&\text{
if
}(\sigma_{F}^{2}/\bar{F}^{2})\gamma_{F}^{-1}\delta_{\Theta}^{2}\ll\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1}\ll(\sigma_{F}^{2}/\bar{F}^{2})\gamma_{F}^{-1}\delta_{\Theta}\\\
O(N^{-1/2})&\text{ if
}(\sigma_{F}^{2}/\bar{F}^{2})\gamma_{F}^{-1}\delta_{\Theta}^{2}\gg\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1}.\end{cases}$

This is small when the colony size is large, and particularly small if
$\delta_{\Theta}$ is small and the relative strength of force magnitude
fluctuations is not much larger that those of the force orientation
fluctuations.

For typical situations within the paradigm of small stochastic variations in
the flagellar forces and orientation variations from the normal to the
surface, $\tilde{D_{\mathrm{t}}}$ should not contribute substantially relative
to the simple terms in (26). Indeed, for the biophysical parameters listed in
Table 1, the ratio in Eq. (28) is approximately $0.04N^{-1/2}$. This
correction term is not sufficiently potent to rectify the discrepancies seen
in Figure 9.

Note the terms $\tilde{D_{\mathrm{t}}}$ arise from correlations between
individual flagellar dynamics and the rotation of the colony as a whole. It
appears to be related to the wobbling of swimming choanoflagellate cells
observed in the computational simulations of Nguyen et al. (2019), though we’d
expect the wobbling to be reduced for colonies due to the increased rotational
drag. The argument for the smallness of $\tilde{D_{\mathrm{t}}}$ relies
largely on the disclike colony configuration, with nearly normal flagella, and
does not apply to general colony configurations. Indeed, the rotational drift
should become even stronger relative to the rotational diffusivity in less
symmetric colonies. Thus the terms which generalize $\tilde{D_{\mathrm{t}}}$
to more general nearly planar colony configurations are retained in Eq. (64).

##  4 Methods

In this section we will present the detailed derivations for the transport
statistics reported in Section 3 for the disc colony described in Section 2.3.
Rotational statistics are computed in Subsection 4.1, and translational
statistics in Subsection 4.2. All statistical averaging operations here,
denoted by angle brackets $\langle\cdot\rangle$, are over dynamical noise
only, for fixed colony parameters. Thus, in contrast to Section 3, we will use
Cov and Var to denote covariance and variance of random variables with respect
to the dynamical noise. The subsequent demographic averages are all
straightforward except for the mean-square speed; the demographic statistic
for that quantity is derived in detail in Appendix C.

###  4.1 Rotational Statistics

Under the model described in Subsection 2.2, the statistics for the flagellar
force magnitude $F_{i}(t)$ and orientation $\Theta_{i}(t)$ are all
independent, statistically stationary Gaussian random processes with means
$\langle F_{i}\rangle=F_{i}^{(0)}$ and
$\langle\Theta_{i}\rangle=\Theta_{i}^{(0)}$ and correlation functions

$\text{Cov}(\Theta_{i}(t),\Theta_{i}(t^{\prime}))=(\sigma_{\Theta}^{2})e^{-\gamma_{\Theta}|t-t^{\prime}|},\qquad\qquad\text{Cov}(F_{i}(t),F_{i}(t^{\prime}))=(\sigma_{F}^{2})e^{-\gamma_{F}|t-t^{\prime}|}$
(29)

In particular, the correlation time of $\Theta_{i}$ is
$\frac{1}{\gamma_{\Theta}}$ and for $F_{i}$ is $\frac{1}{\gamma_{F}}$.

We proceed next to compute the first and second order statistics of the total
torque (7) on the colony. Beginning with the mean, we use

the independence of $F_{i}$ and $\Theta_{i}$ and Eq. (78) to write:

$\displaystyle\langle T^{(c)}(t)\rangle=-\sum_{i=1}^{N}\langle F_{i}(t)\rangle
a\langle\sin(\Theta_{i}(t))\rangle=-\sum_{i=1}^{N}F_{i}^{(0)}a\sin(\Theta_{i}^{(0)})e^{-\frac{1}{2}\sigma_{\Theta}^{2}}$
(30)

This mean torque induces the mean rotational velocity for a given colony (see
Eq. (5)):

$\Omega^{\ast}_{\mathrm{r}}\equiv\lim_{t\rightarrow\infty}\frac{\langle\Theta^{(\mathrm{c})}(t)-\Theta^{(\mathrm{c})}(0)\rangle}{t}=-\frac{a}{\gamma_{\mathrm{r}}}e^{-\frac{1}{2}\sigma_{\Theta}^{2}}\sum_{i=1}^{N}F_{i}^{(0)}\sin(\Theta_{i}^{(0)})$

To compute the variance of the colony orientation, we require the
autocorrelation function of the colony torque. We begin by using independence
of all component stochastic processes to write:

$\displaystyle\text{Cov}[T^{(c)}(t),T^{(c)}(t^{\prime})]$
$\displaystyle=\text{Cov}[-\sum_{i=1}^{N}F_{i}(t)a\sin(\Theta_{i}(t)),-\sum_{j=1}^{N}F_{j}(t^{\prime})a\sin(\Theta_{j}(t^{\prime}))]$
(31)
$\displaystyle=\sum_{j=1}^{N}\text{Cov}[F_{j}(t)a\sin(\Theta_{j}(t)),F_{j}(t^{\prime})a\sin(\Theta_{j}(t^{\prime}))]$
$\displaystyle=\sum_{j=1}^{N}\left[\langle(F_{j}(t)a\sin\Theta_{j}(t))(F_{j}(t^{\prime})a\sin\Theta_{j}(t^{\prime}))\rangle-\langle
F_{j}(t)a\sin\Theta_{j}(t)\rangle\langle
F_{j}(t^{\prime})a\sin\Theta_{j}(t^{\prime})\rangle\right]$
$\displaystyle=a^{2}\sum_{j=1}^{N}\langle
F_{j}(t)F_{j}(t^{\prime})\rangle\frac{1}{2}\langle\cos(\Theta_{j}(t)-\Theta_{j}(t^{\prime}))-\cos(\Theta_{j}(t)+\Theta_{j}(t^{\prime}))\rangle$
$\displaystyle\qquad\qquad-a^{2}\sum_{j=1}^{N}(F_{j}^{(0)})^{2}\sin^{2}\Theta_{j}^{(0)}\mathrm{e}^{-\sigma_{\Theta}^{2}}$

Next, we compute each summand using the formula (78) for the average of the
sinusoid of a Gaussian random variable:

$\displaystyle\text{Cov}[F_{j}(t)a\sin(\Theta_{j}(t)),F_{j}(t^{\prime})a\sin(\Theta_{j}(t^{\prime}))]$
$\displaystyle\qquad\qquad=\frac{1}{2}a^{2}\left(\sigma_{F}^{2}\mathrm{e}^{-\gamma_{F}|t-t^{\prime}|}+(F_{j}^{(0)})^{2}\right)\left(\mathrm{e}^{-\frac{1}{2}\text{Var}[\Theta_{j}(t)-\Theta_{j}(t^{\prime})]}-\cos(2\Theta_{j}^{(0)})\mathrm{e}^{-\frac{1}{2}\text{Var}[\Theta_{j}(t)+\Theta_{j}(t^{\prime})]}\right)$
$\displaystyle\qquad\qquad\qquad\qquad-a^{2}(F_{j}^{(0)})^{2}\sin^{2}\Theta_{j}^{(0)}\mathrm{e}^{-\sigma_{\Theta}^{2}}$
$\displaystyle\qquad\qquad=\frac{a^{2}}{2}\left((\sigma_{F}^{2})e^{-\gamma_{F}|t-t^{\prime}|}+(F_{j}^{(0)})^{2}\right)(e^{[-\sigma_{\Theta}^{2}+\sigma_{\Theta}^{2}e^{-\gamma_{\Theta}|t-t^{\prime}|}]}-\cos(2\Theta_{j}^{(0)})e^{[-\sigma_{\Theta}^{2}-\sigma_{\Theta}^{2}e^{-\gamma_{\Theta}|t-t^{\prime}|}]})$
$\displaystyle\qquad\qquad\qquad\qquad-a^{2}(F_{j}^{(0)})^{2}\frac{1}{2}(1-\cos(2\Theta_{i}^{(0)})))e^{-\sigma_{\Theta}^{2}}$

We thereby obtain

$\displaystyle\text{Cov}[T^{(c)}(t),T^{(c)}(t^{\prime})]$ (32)
$\displaystyle\qquad\sum_{j=1}^{N}\frac{a^{2}}{2}((\sigma_{F}^{2})e^{-\gamma_{F}|t-t^{\prime}|}+(F_{j}^{(0)})^{2})(e^{[-\sigma_{\Theta}^{2}+\sigma_{\Theta}^{2}e^{-\gamma_{\Theta}|t-t^{\prime}|}]}-\cos(2\Theta_{j}^{(0)})e^{[-\sigma_{\Theta}^{2}-\sigma_{\Theta}^{2}e^{-\gamma_{\Theta}|t-t^{\prime}|}]})$
$\displaystyle\qquad\qquad-\frac{a^{2}}{2}(F_{j}^{(0)})^{2}(1-\cos(2\Theta_{j}^{(0)})))e^{-\sigma_{\Theta}^{2}}$

We can then express the evolution of the variance of the colony orientation by
noting the independence of the active torque and thermally driven rotational
diffusion:

$\displaystyle\text{Var}(\Theta^{(c)}(t))$
$\displaystyle=\text{Var}(\frac{1}{\gamma_{\mathrm{r}}^{2}}\int_{0}^{t}T^{(c)}(t^{\prime})dt^{\prime})+\text{Var}(\sqrt{2D_{\mathrm{r}}}W^{\Theta,c}(t))$
(33)
$\displaystyle=\left(\frac{1}{\gamma_{\mathrm{r}}}\right)^{2}\int_{0}^{t}\int_{0}^{t}\text{Cov}(T^{(c)}(t^{\prime\prime}),T^{(c)}(t^{\prime}))dt^{\prime\prime}dt^{\prime}+2D_{\mathrm{r}}t$

We are primarily interested in the long-time statistical behavior of the
orientation. To this end, we observe that
$\text{Cov}(T^{(c)}(t^{\prime\prime}),T^{(c)}(t^{\prime}))$ is a function only
of $|t^{\prime\prime}-t^{\prime}|$, and that for integrable functions $g$, a
change of variables to $u=t^{\prime\prime}-t^{\prime}$ and
$u^{\prime}=t^{\prime\prime}+t^{\prime}$ gives:

$\displaystyle\int\limits_{0}^{t}\int\limits_{0}^{t}g(|t^{\prime\prime}-t^{\prime}|)dt^{\prime}dt^{\prime\prime}$
$\displaystyle=\frac{1}{2}\int\limits_{-t}^{t}\int\limits_{|u|}^{2t-|u|}g(|u|)\mathrm{d}u^{\prime}\,\mathrm{d}u=2\int\limits_{0}^{t}(t-u)g(u)\,\mathrm{d}u\sim
2t\int_{0}^{\infty}g(u)\,\mathrm{d}u\text{ as }t\rightarrow\infty.$ (34)

To evaluate the time integral, we use the following formulas, which involve a
change of integration variable to $u=\mathrm{e}^{-kt}$

$\displaystyle\int_{0}^{\infty}(\mathrm{e}^{c\mathrm{e}^{-kt}}-1)\,\mathrm{d}t$
$\displaystyle=\int_{0}^{1}\frac{\mathrm{e}^{cu}-1}{ku}\,\mathrm{d}u=-k^{-1}\operatorname{Ein}(-c)\text{
for }c,k\in\mathbb{R},k>0,$
$\displaystyle\int_{0}^{\infty}\mathrm{e}^{-bt+c\mathrm{e}^{-kt}}\,\mathrm{d}t$
$\displaystyle=k^{-1}\int_{0}^{1}u^{b/k-1}\mathrm{e}^{cu}\,\mathrm{d}u=\frac{1}{b}{}_{1}F_{1}(b/k,b/k+1,c)$
$\displaystyle=\frac{1}{b}\mathrm{e}^{c}{}_{1}F_{1}(1,b/k+1,-c)\text{ for
}b,c,k\in\mathbb{R},b,k>0$

where we have used integral representations and identities for the entire
exponential integral $\operatorname{Ein}$ and the confluent hypergeometric
function of the first kind ${}_{1}F_{1}$. Applying the asymptotic integral
formula (34) with these integral identities to the integral of the active
torque correlation function (32) yields the result (10).

###  4.2 Translational Statistics Derivations

To compute the translational mobility of the colony, we next examine the force
on the colony (6). We will use complex notation for expressing vector
calculations, with $\mathrm{i}$ in these calculations representing the
imaginary unit.

$\displaystyle\mathbf{F}^{(c)}(t)$
$\displaystyle=\sum_{j=1}^{N}-F_{j}(t)\begin{pmatrix}\cos\left(2\pi\left[\frac{j-\frac{1}{2}+\frac{S_{j}}{l}}{N}\right]+\Theta_{j}(t)+\Theta^{(c)}(t)\right)\\\
\sin\left(2\pi\left[\frac{j-\frac{1}{2}+\frac{S_{j}}{l}}{N}\right]+\Theta_{j}(t)+\Theta^{(c)}(t)\right)\end{pmatrix}$
$\displaystyle\sim-\sum_{j=1}^{N}F_{j}(t)e^{\mathrm{i}\left[2\pi\left(\frac{j-\frac{1}{2}+\frac{S_{j}}{l}}{N}\right)+\Theta_{j}(t)+\Theta^{(c)}(t)\right]}$
(35)

While individual colonies will have asymmetrically arranged flagella, the
isotropic enviroment will ensure that the colony orientation
$\Theta^{(\mathrm{c})}(t)$ will eventually become uniformly distributed,
implying an eventual zero mean to the total force vector on the colony. In
more detail, as the calculations in Ashenafi (2021) show, the mean force on
the colony relaxes on a time scale
$\max(\gamma_{F}^{-1},\gamma_{\Theta}^{-1})$ from the initial flagellar forces
and orientations to a flagellar-averaged form, which then rotates with drift
$\Omega^{\ast}_{\mathrm{r}}$ (8) and dissipates at a rate
$D_{\mathrm{r}}^{\ast}$ due to rotational diffusion. The flagellar-averaged
root-mean-square force will be recovered in the calculation of the colony
force covariance to follow. The translational drift of the colony is
proportional to the colony force, and thus will eventually tend to zero due to
the randomization of its direction.

The homogeneity and isotropy of the environment implies the long-time dynamics
will also be statistically isotropic, so we have

$\text{Cov}(\mathbf{F}^{(\mathrm{c})}(t),\mathbf{F}^{(\mathrm{c})}(t^{\prime}))\sim\frac{1}{2}\langle\mathbf{F}^{(\mathrm{c})}(t)\cdot\mathbf{F}^{(\mathrm{c})}(t^{\prime})\rangle\mathit{\mathsf{I}}\text{
for }t,t^{\prime}\rightarrow\infty.$ (36)

where $\mathit{\mathsf{I}}$ denotes the identity matrix. Using the complex
representation of the vectorial colony force (4.2), we have:

$\langle|\mathbf{F}^{(\mathrm{c})}(t)|^{2}\rangle=\sum_{j=1}^{N}\sum_{j^{\prime}=1}^{N}\mathrm{e}^{\frac{2\pi\mathrm{i}}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)}\langle
F_{j}(t)F_{j^{\prime}}(t)\mathrm{e}^{\mathrm{i}(\Theta_{j}(t)-\Theta_{j^{\prime}}(t))}\rangle.$

Now the flagellar forces and orientations are assumed independently
distributed with respect to their Gaussian stationary distributions. Thus, we
can evaluate the average over the orientation fluctuations using the formula
(74) for the average of the exponential of a Gaussian random variable.

$\displaystyle\langle|\mathbf{F}^{(\mathrm{c})}(t)|^{2}\rangle$
$\displaystyle=\sum_{j=1}^{N}\sum_{j^{\prime}=1}^{N}\mathrm{e}^{\frac{2\pi\mathrm{i}}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)}\langle
F_{j}(t)F_{j^{\prime}}(t)\rangle\langle\mathrm{e}^{\mathrm{i}(\Theta_{j}(t)-\Theta_{j^{\prime}}(t))}\rangle$
$\displaystyle\qquad\qquad=\sum_{j=1}^{N}(\sigma_{F}^{2}+(F_{j}^{(0)})^{2})+\sum_{j=1}^{N}\sum_{\begin{subarray}{c}j^{\prime}=1\\\
j^{\prime}\neq
j\end{subarray}}^{N}\mathrm{e}^{\frac{2\pi\mathrm{i}}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)}F_{j}^{(0)}F_{j^{\prime}}^{(0)}\mathrm{e}^{\mathrm{i}(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)})}\mathrm{e}^{-\sigma_{\Theta}^{2}}$
$\displaystyle=N\sigma_{F}^{2}+\left[e^{-\sigma_{\Theta}^{2}}\left|\sum_{j=1}^{N}F_{j}^{(0)}\mathrm{e}^{\mathrm{i}\Theta_{j}^{(0)}+\frac{2\pi\mathrm{i}}{N}\left(j+\frac{S_{j}}{l}\right)}\right|^{2}+(1-\mathrm{e}^{-\sigma_{\Theta}^{2}})\sum_{j=1}^{N}\left(F_{j}^{(0)}\right)^{2}\right]$

Using the proportionality of the colony velocity (neglecting the thermal
contribution) to the colony force via the translational drag, we obtain the
instantaneous mean-square velocity (4.2) for a colony.

To obtain the effective translational diffusivity of the colony requires
integrating the autocorrelation function of the colony velocity, which in turn
requires us to consider the colony force correlation function at unequal times
$t,t^{\prime}$:

$\langle\mathbf{F}^{(\mathrm{c})}(t)\cdot\mathbf{F}^{(\mathrm{c})}(t^{\prime})\rangle=\operatorname{\mathrm{Re}}\sum_{j=1}^{N}\sum_{j^{\prime}=1}^{N}\mathrm{e}^{\frac{2\pi\mathrm{i}}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)}\langle
F_{j}(t)F_{j^{\prime}}(t^{\prime})\mathrm{e}^{\mathrm{i}(\Theta_{j}(t)+\Theta^{(\mathrm{c})}(t)-\Theta_{j^{\prime}}(t^{\prime})-\Theta^{(\mathrm{c})}(t^{\prime}))}\rangle$
(37)

Here the calculation in general is complicated relative to our previous ones
due to the correlations between the colony rotation
$\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t)$ and the flagellar
forces and orientations. Moreover, the calculation of the average is
complicated by the fact that $\Theta^{(\mathrm{c})}(t)$ is non-Gaussian. We
proceed by assuming the dynamical force fluctuation standard deviation
$\sigma_{F}$ is small relative to the magnitude of the flagellar forces
$F_{j}^{(0)}$, which we denote by $\bar{F}$, and the dynamical orientation
standard deviation $\sigma_{\Theta}$ is also small. Our attempts to split off
higher terms by Taylor expansion led, when substituted into the integral (50),
to pessimistic error bounds that compete with significant terms. We thus take
a different approximation route by treating
$\Theta^{(\mathrm{c})}(t)-\Theta^{(\mathrm{c})}(t^{\prime}),\Theta_{j}(t),\Theta_{j^{\prime}}(t^{\prime}),F_{j}(t),F_{(}j^{\prime})(t^{\prime})$
as jointly Gaussian random variables. We use the notation
$\langle\cdot\rangle_{\mathrm{G}}$ to denote an average computed under this
approximation. The assumption here is on the change in colony orientation
$\Theta^{(\mathrm{c})}(t)-\Theta^{(\mathrm{c})}(t^{\prime})$, which is non-
Gaussian because the colony torque (7) is not Gaussian. But they are Gaussian
through first order in the flagellar force and orientation fluctuations, which
we are assuming to be small. Moreover, by the functional central limit
theorem, $\Theta^{(\mathrm{c})}(t)-\Theta^{(\mathrm{c})}(t^{\prime})$ can be
treated as approximately Gaussian when
$|t-t^{\prime}|\gg\max(\gamma_{F}^{-1},\gamma_{\Theta}^{-1})$, the correlation
time of the colony torque (see Eq. (32)). Actually, this does not tell us that
$\Theta^{(\mathrm{c})}(t)-\Theta^{(\mathrm{c})}(t^{\prime})$ is _jointly_
Gaussian with the other random variables, but as we will see, they are weakly
correlated so this concern should not be of much import.

Proceeding with this plausible jointly Gaussian approximation, we have from
the formulas in Appendix B:

$\displaystyle\langle
F_{j}(t)F_{j^{\prime}}(t^{\prime})\mathrm{e}^{\mathrm{i}(\Theta_{j}(t)+\Theta^{(\mathrm{c})}(t)-\Theta_{j^{\prime}}(t^{\prime})-\Theta^{(\mathrm{c})}(t^{\prime}))}\rangle_{G}=\mathrm{e}^{\mathrm{i}\left(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)}+\langle\Theta^{(\mathrm{c})}(t)-\Theta^{(\mathrm{c})}(t^{\prime})\rangle\right)}$
(38)
$\displaystyle\qquad\qquad=\left(F_{j}^{(0)}-\mathrm{i}\text{Cov}(F_{j}(t),\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t))\right)\left(F_{j^{\prime}}^{(0)}-\mathrm{i}\text{Cov}(F_{j^{\prime}}(t^{\prime}),\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t)\right)$
$\displaystyle\qquad\qquad\qquad\qquad\cdot\mathrm{e}^{-\frac{1}{2}(2\sigma_{\Theta}^{2}+\text{Var}(\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t)))-2\text{Cov}(\Theta_{j}(t)-\Theta_{j^{\prime}}(t^{\prime}),\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t))}\text{
for }j\neq j^{\prime}$

We next express the colony rotation in terms of its mean (first term) and
fluctuating behavior:

$\displaystyle\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t)=\Omega^{\ast}_{\mathrm{r}}(t^{\prime}-t)-\int_{t}^{t^{\prime}}\frac{a}{\gamma_{\mathrm{r}}}\sum_{k=1}^{N}[F_{k}(t^{\prime\prime})\sin(\Theta_{k}(t^{\prime\prime}))-F_{k}^{(0)}\sin(\Theta_{k}^{(0)})\mathrm{e}^{-\sigma_{\Theta}^{2}/2}]\,\mathrm{d}t^{\prime\prime}$
$\displaystyle\qquad\qquad\qquad\qquad+\sqrt{2D_{\mathrm{r}}}(W^{\Theta,c}(t^{\prime})-W^{\Theta,c}(t))$

and then decompose this expression into a Gaussian component and a smaller
non-Gaussian component:

$\displaystyle\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t)=\Delta\Theta^{(\mathrm{c})}_{\mathrm{G}}(t^{\prime},t)+\Delta\Theta^{(\mathrm{c})}_{\mathrm{NG}}(t^{\prime},t),$
(39)
$\displaystyle\Delta\Theta^{(\mathrm{c})}_{\mathrm{G}}(t^{\prime},t)=\Omega^{\ast}_{\mathrm{r}}(t^{\prime}-t)+\sqrt{2D_{\mathrm{r}}}(W^{\Theta,c}(t^{\prime})-W^{\Theta,c}(t))$
(40)
$\displaystyle\qquad\qquad\qquad\qquad-\frac{a}{\gamma_{\mathrm{r}}}\int_{t}^{t^{\prime}}\sum_{k=1}^{N}\left[F_{k}^{(0)}\cos(\Theta_{k}^{(0)})(\Theta_{k}(t^{\prime\prime})-\Theta_{k}^{(0)})+(F_{k}(t^{\prime\prime})-F_{k}^{(0)})\sin(\Theta_{k}^{(0)})\mathrm{e}^{-\sigma_{\Theta}^{2}/2}\right]\,\mathrm{d}t^{\prime\prime}$
$\displaystyle\Delta\Theta^{(\mathrm{c})}_{\mathrm{NG}}(t^{\prime},t)=-\sum_{k=1}^{N}\int_{t}^{t^{\prime}}\frac{a}{\gamma_{\mathrm{r}}}F_{k}(t^{\prime\prime})\left(\sin\Theta_{k}(t^{\prime\prime})-\sin(\Theta_{k}^{(0)})\mathrm{e}^{-\sigma_{\Theta}^{2}/2}-\cos(\Theta_{k}^{(0)})(\Theta_{k}(t^{\prime\prime})-\Theta_{k}^{(0)})\right)\,\mathrm{d}t^{\prime\prime}$
$\displaystyle\qquad\qquad\qquad\qquad-\sum_{k=1}^{N}\int_{t}^{t^{\prime}}\frac{a}{\gamma_{\mathrm{r}}}(F_{k}(t^{\prime\prime})-F_{k}^{(0)})\cos(\Theta_{k}^{(0)})(\Theta_{k}(t^{\prime\prime})-\Theta_{k}^{(0)})\,\mathrm{d}t^{\prime\prime}$

The covariances in Eq. (38) can now be seen to be clearly dominated by just
the covariance with the Gaussian component so we compute:

$\displaystyle\text{Cov}(F_{j}(t),\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t))=\text{Cov}(F_{j}(t),\Delta\Theta^{(\mathrm{c})}_{\mathrm{G}}(t^{\prime},t))+\text{Cov}(F_{j}(t),\Delta\Theta^{(\mathrm{c})}_{\mathrm{NG}}(t^{\prime},t))$
$\displaystyle\qquad\qquad=-\frac{a}{\gamma_{\mathrm{r}}}\int_{t}^{t^{\prime}}\left[\sigma_{F}^{2}\sin\Theta_{j}^{(0)}+O(\sigma_{F}^{2}\sigma_{\Theta}^{4})\right]\mathrm{e}^{-\sigma_{\Theta}^{2}/2}\mathrm{e}^{-\gamma_{F}|t^{\prime\prime}-t|}\,\mathrm{d}t^{\prime\prime}$
$\displaystyle\qquad\qquad=\frac{a\sigma_{F}^{2}\sin\Theta_{j}^{(0)}}{\gamma_{\mathrm{r}}}\frac{\operatorname{sgn}(t^{\prime}-t)(\mathrm{e}^{-\gamma_{F}|t-t^{\prime}|}-1)}{\gamma_{F}}\left(1+O(\sigma_{\Theta}^{2})\right),$
(41)
$\displaystyle\text{Cov}(\Theta_{j}(t)-\Theta_{j^{\prime}}(t^{\prime}),\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t))$
$\displaystyle\qquad\qquad=\text{Cov}(\Theta_{j}(t)-\Theta_{j^{\prime}}(t^{\prime}),\Delta\Theta^{(\mathrm{c})}_{\mathrm{G}}(t^{\prime},t))+\text{Cov}(\Theta_{j}(t)-\Theta_{j^{\prime}}(t^{\prime}),\Delta\Theta^{(\mathrm{c})}_{\mathrm{NG}}(t^{\prime},t))$
$\displaystyle\qquad\qquad=-\int_{t}^{t^{\prime}}a\sigma_{\Theta}^{2}\frac{F_{j}^{(0)}\cos(\Theta_{j}^{(0)})\mathrm{e}^{-\gamma_{\Theta}|t-t^{\prime\prime}|}-F_{j^{\prime}}^{(0)}\cos(\Theta_{j^{\prime}}^{(0)})\mathrm{e}^{-\gamma_{\Theta}|t^{\prime}-t^{\prime\prime}|}}{\gamma_{\mathrm{r}}}\left(1+O(\sigma_{\Theta}^{2})\right)$
$\displaystyle\qquad\qquad=-\frac{a\sigma_{\Theta}^{2}\left(F_{j}^{(0)}\cos(\Theta_{j}^{(0)})-F_{j^{\prime}}^{(0)}\cos(\Theta_{j^{\prime}}^{(0)})\right)}{\gamma_{\mathrm{r}}}\frac{\operatorname{sgn}(t^{\prime}-t)(1-\mathrm{e}^{-\gamma_{\Theta}|t-t^{\prime}|})}{\gamma_{\Theta}}\left(1+O(\sigma_{\Theta}^{2})\right)$
(42)

and similarly

$\text{Cov}(F_{j^{\prime}}(t^{\prime}),\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t))=\frac{a\sigma_{F}^{2}\sin\Theta_{j^{\prime}}^{(0)}}{\gamma_{\mathrm{r}}}\frac{\operatorname{sgn}(t^{\prime}-t)(\mathrm{e}^{-\gamma_{F}|t-t^{\prime}|}-1)}{\gamma_{F}}\left(1+O(\sigma_{\Theta}^{2})\right)$
(43)

We estimate the variance of the colony rotation in terms of its long term-
asymptotics using a generalization of Eq. (33):

$\displaystyle\text{Var}(\Theta^{(\mathrm{c})}(t^{\prime})-\Theta^{(\mathrm{c})}(t))$
$\displaystyle=\frac{2}{\gamma_{\mathrm{r}}^{2}}\int_{0}^{|t^{\prime}-t|}(|t^{\prime}-t|-u)\text{Cov}(T^{(c)}(t),T^{(c)}(t+u))\,\mathrm{d}u+2D_{\mathrm{r}}|t^{\prime}-t|$
$\displaystyle=2D_{\mathrm{r}}^{\ast}|t^{\prime}-t|+V_{\Delta\Theta,\mathrm{G}}(|t^{\prime}-t|)$

with the error term estimated via Taylor expansion of the torque correlation
function (32) with respect to the small parameters $\sigma_{\Theta}^{2}$ and
$\sigma_{F}^{2}/\bar{F}^{2}$

$\displaystyle V_{\Delta\Theta,\mathrm{G}}(\tau)$
$\displaystyle\equiv-\frac{2\tau}{\gamma_{\mathrm{r}}^{2}}\int_{\tau}^{\infty}\text{Cov}(T^{(c)}(t),T^{(c)}(t+u))\,\mathrm{d}u-\frac{2}{\gamma_{\mathrm{r}}^{2}}\int_{0}^{\tau}u\text{Cov}(T^{(c)}(t),T^{(c)}(t+u))\,\mathrm{d}u$
$\displaystyle\sim
O\left(\frac{Na^{2}}{\gamma_{\mathrm{r}}^{2}\gamma_{-}}(\sigma_{F}^{2}+\bar{F}^{2}\sigma_{\Theta}^{2})\left[\tau\mathrm{e}^{-\gamma_{-}\tau}+\gamma_{-}^{-1}\right]\right)\sim
O\left(\frac{D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}}}{\gamma_{-}}\left[\tau\mathrm{e}^{-\gamma_{-}\tau}+\gamma_{-}^{-1}\right]\right)$
$\displaystyle V_{\Delta\Theta,\mathrm{G}}^{\prime}(\tau)$ $\displaystyle\sim
O\left(\frac{Na^{2}}{\gamma_{\mathrm{r}}^{2}\gamma_{-}}(\sigma_{F}^{2}+\bar{F}^{2}\sigma_{\Theta}^{2})\left[(1+\gamma_{-}\tau)\mathrm{e}^{-\gamma_{-}\tau}\right]\right)\sim
O(\left(\frac{D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}}}{\gamma_{-}}\left[(1+\gamma_{-}\tau)\mathrm{e}^{-\gamma_{-}\tau}\right]\right)$
(44)

Substituting Eqs. (41), (42), and (43) into Eq. (38), we obtain

$\displaystyle\langle
F_{j}(t)F_{j^{\prime}}(t^{\prime})\mathrm{e}^{\mathrm{i}(\Theta_{j}(t)+\Theta^{(\mathrm{c})}(t)-\Theta_{j^{\prime}}(t^{\prime})-\Theta^{(\mathrm{c})}(t^{\prime}))}\rangle_{G}$
(45)
$\displaystyle=\left(F_{j}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin(\Theta_{j}^{(0)})\operatorname{sgn}(t^{\prime}-t)(1-\mathrm{e}^{-\gamma_{F}|t^{\prime}-t|})(1+O(\sigma_{\Theta}^{2}))\right)$
$\displaystyle\qquad\qquad\cdot\left(F_{j^{\prime}}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin(\Theta_{j^{\prime}}^{(0)})\operatorname{sgn}(t^{\prime}-t)(1-\mathrm{e}^{-\gamma_{F}|t^{\prime}-t|})(1+O(\sigma_{\Theta}^{2}))\right)$
$\displaystyle\qquad\qquad\cdot
e^{\mathrm{i}[\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)}]-\sigma_{\Theta}^{2}-\mathrm{i}\Omega^{\ast}_{\mathrm{r}}(t^{\prime}-t)-(D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}})|t^{\prime}-t|-\frac{1}{2}V_{\Delta\Theta,\mathrm{G}}(|t^{\prime}-t|)}$
$\displaystyle\qquad\qquad\cdot\mathrm{e}^{-\frac{a\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{\Theta}}(F_{j}^{(0)}\cos(\Theta_{j}^{(0)})-F_{j^{\prime}}^{(0)}\cos(\Theta_{j^{\prime}}^{(0)}))\operatorname{sgn}(t^{\prime}-t)(1-\mathrm{e}^{-\gamma_{\Theta}(t^{\prime}-t)})(1+O(\sigma_{\Theta}^{2}))}\text{
for }j\neq j^{\prime}.$

A similar calculation gives:

$\displaystyle\langle
F_{j}(t)F_{j^{\prime}}(t^{\prime})\mathrm{e}^{\mathrm{i}(\Theta_{j}(t)+\Theta^{(\mathrm{c})}(t)-\Theta_{j^{\prime}}(t^{\prime})-\Theta^{(\mathrm{c})}(t^{\prime}))}\rangle_{G}e^{-\sigma_{\Theta}^{2}(1-\mathrm{e}^{-\gamma_{\Theta}(t^{\prime}-t)})-\mathrm{i}\Omega^{\ast}_{\mathrm{r}}(t^{\prime}-t)-(D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}})|t^{\prime}-t|-\frac{1}{2}V_{\Delta\Theta,\mathrm{G}}|t^{\prime}-t|}$
(46)
$\displaystyle\qquad\cdot\left(\left(F_{j}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin(\Theta_{j}^{(0)})\operatorname{sgn}(t^{\prime}-t)(1-\mathrm{e}^{-\gamma_{F}(t^{\prime}-t)})(1+O(\sigma_{\Theta}^{2}))\right)^{2}+\sigma_{F}^{2}\mathrm{e}^{-\gamma_{F}(t^{\prime}-t)}\right)\text{
for }j=j^{\prime}$

The expressions in (46) and (45) can be simplified when the time lag
$|t^{\prime}-t|$ is large compared to the time scale of flagellar force and
orientation relaxation:

$\displaystyle\langle
F_{j}(t)F_{j^{\prime}}(t^{\prime})\mathrm{e}^{\mathrm{i}(\Theta_{j}(t)+\Theta^{(\mathrm{c})}(t)-\Theta_{j^{\prime}}(t^{\prime})-\Theta^{(\mathrm{c})}(t^{\prime}))}\rangle_{G}=e^{-\sigma_{\Theta}^{2}-\mathrm{i}\Omega^{\ast}_{\mathrm{r}}(t^{\prime}-t)-D_{\mathrm{r}}^{\ast}|t^{\prime}-t|-\frac{1}{2}V_{\Delta\Theta,\mathrm{G}}(|t^{\prime}-t|)}$
(47)
$\displaystyle\qquad\cdot\left(F_{j}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin\Theta_{j}^{(0)}\operatorname{sgn}(t^{\prime}-t)(1+O(\sigma_{\Theta}^{2}))\right)\left(F_{j^{\prime}}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin\Theta_{j^{\prime}}^{(0)}\operatorname{sgn}(t^{\prime}-t)(1+O(\sigma_{\Theta}^{2}))\right)$
$\displaystyle\qquad\cdot\mathrm{e}^{\mathrm{i}(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)})}\mathrm{e}^{-\frac{a\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{\Theta}}\operatorname{sgn}(t^{\prime}-t)(F_{j}^{(0)}\cos\Theta_{j}^{(0)}-F_{j^{\prime}}^{(0)}\cos\Theta_{j^{\prime}}^{(0)})(1+O(\sigma_{\Theta}^{2}))}$
$\displaystyle\qquad\qquad\qquad\qquad\text{ for
}|t^{\prime}-t|\gg\gamma_{-}^{-1}.$

We note for what follows that the error in using this expression also for
short time lags can be estimated as:

$O\left(\left(1+\frac{a\bar{F}}{\gamma_{\mathrm{r}}\gamma_{-}}\right)\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]\right).$
(48)

Continuing with the Gaussian approximation and decorating averages with a
subscript G to denote this approximation, we examine the mean-square
displacement over observation intervals large compared to the flagellar
dynamics but short compared to the colony dynamics:

$\mathbb{E}|\mathbf{X}^{(\mathrm{c})}(t+\tau)-\mathbf{X}^{(\mathrm{c})}(t)|^{2}=\Big{(}\frac{1}{\gamma_{\mathrm{t}}}\Big{)}^{2}\int_{t}^{t+\tau}\int_{t}^{t+\tau}\mathbf{F}^{(\mathrm{c})}(t^{\prime})\cdot\mathbf{F}^{(\mathrm{c})}(t^{\prime\prime}))dt^{\prime}\,dt^{\prime\prime}+4D_{\mathrm{t}}\tau.$

Provided $\tau$ is short compared to the time scale on which the colony as a
whole rotates, the force correlation function in the integrand simply adjusts
from its zero lag value to a different constant value inferred from
substituting Eq. (47) into Eq. (37). Approximating the integrand by the second
constant and noting the error in this approximation from Eq. (48), we obtain:

$\displaystyle\mathbb{E}|\mathbf{X}^{(\mathrm{c})}(t+\tau)-\mathbf{X}^{(\mathrm{c})}(t)|^{2}\sim
4D_{\mathrm{t}}\tau$ (49)
$\displaystyle\qquad+\frac{e^{-\sigma_{\Theta}^{2}}\tau^{2}}{2\gamma_{\mathrm{t}}^{2}}\operatorname{\mathrm{Re}}\sum_{j=1}^{N}\sum_{j^{\prime}=1}^{N}\left[\left(F_{j}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin\Theta_{j}^{(0)}(1+O(\sigma_{\Theta}^{2}))\right)\left(F_{j^{\prime}}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin\Theta_{j^{\prime}}^{(0)}(1+O(\sigma_{\Theta}^{2}))\right)\right.$
$\displaystyle\qquad\qquad\cdot\mathrm{e}^{\mathrm{i}(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)})+\frac{2\pi\mathrm{i}}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)}\mathrm{e}^{-\frac{a\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{\Theta}}(F_{j}^{(0)}\cos\Theta_{j}^{(0)}-F_{j^{\prime}}^{(0)}\cos\Theta_{j^{\prime}}^{(0)})(1+O(\sigma_{\Theta}^{2}))}$
$\displaystyle\qquad\qquad+O\left(\left(1+\frac{a\bar{F}}{\gamma_{\mathrm{r}}\gamma_{-}}\right)\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]\gamma_{\mathrm{t}}^{-2}\gamma_{-}^{-1}\tau\right)\text{
for
}\gamma_{-}^{-1}\ll\tau\ll\Omega^{\ast}_{\mathrm{r}}{}^{-1},D_{\mathrm{r}}^{\ast}{}^{-1}$

The term proportional to $\tau^{2}$ term corresponds to a second ballistic
phase, whose coefficient gives the mean-square speed that would be estimated
from observations taken at a resolution longer than that of the flagellar
dynamics. Expanding this expression with respect to $\sigma_{F}$ and
$\sigma_{\Theta}$, the complicated terms at second order all cancel out in the
double sum, and we arrive at the value reported in Eq. (19), with a relative
error

$\sim
O\left(\left(1+\frac{a\bar{F}}{\gamma_{\mathrm{r}}\gamma_{-}}\right)^{2}\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]^{2}\right)$

Indeed, the expression in Eq. (19) would follow by a simple calculation
starting from Eq. (4.2) and freezing the colony rotation. Note the ballistic
term should dominate the error term in Eq. (49) for $\tau\gg\gamma_{-}^{-1}$.

Finally, we can compute the long-time covariance of the displacement of the
colony from integrating the velocity autocorrelation function, which is
directly proportional to the force autocorrelation function, and adding the
translational diffusion component:

$\displaystyle\text{Cov}_{\mathrm{G}}(\mathbf{X}^{(\mathrm{c})}(t),\mathbf{X}^{(\mathrm{c})}(t))$
$\displaystyle=\Big{(}\frac{1}{\gamma_{\mathrm{t}}}\Big{)}^{2}\int_{0}^{t}\int_{0}^{t}\text{Cov}_{\mathrm{G}}(\mathbf{F}^{(\mathrm{c})}(t^{\prime}),\mathbf{F}^{(\mathrm{c})}(t^{\prime\prime}))dt^{\prime}dt^{\prime\prime}+2D_{\mathrm{t}}\mathit{\mathsf{I}}t$
(50)

We use now the expressions (36), (37), and (47) for the covariance of the
colony forcing at large time lag,

We note that the expression (47) exhibits complex conjugation symmetry under
the joint interchange $t\leftrightarrow t^{\prime}$, $j\leftrightarrow
j^{\prime}$, which in turn makes the force autocorrelation function (37)
depend only on the absolute value of the time lag $|t-t^{\prime}|$, as it
must. Thus, we use again Eq. (34) to get:

$\displaystyle\lim_{t\rightarrow\infty}\frac{\text{Cov}_{\mathrm{G}}(\mathbf{X}^{(\mathrm{c})}(t),\mathbf{X}^{(\mathrm{c})}(t))}{2t}$
(51)
$\displaystyle=D_{\mathrm{t}}\mathit{\mathsf{I}}+\frac{1}{2\gamma_{\mathrm{t}}^{2}}\mathit{\mathsf{I}}\int_{0}^{\infty}\operatorname{\mathrm{Re}}\sum_{j=1}^{N}\sum_{j^{\prime}=1}^{N}\mathrm{e}^{\frac{2\pi\mathrm{i}}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)}\langle
F_{j}(0)F_{j^{\prime}}(u)\mathrm{e}^{\mathrm{i}(\Theta_{j}(0)+\Theta^{(\mathrm{c})}(0)-\Theta_{j^{\prime}}(u)-\Theta^{(\mathrm{c})}(u)}\rangle_{G}\,\mathrm{d}u$
$\displaystyle=D_{\mathrm{t}}\mathit{\mathsf{I}}+\operatorname{\mathrm{Re}}\mathit{\mathsf{I}}\frac{e^{-\sigma_{\Theta}^{2}}}{2\gamma_{\mathrm{t}}^{2}}\sum_{j=1}^{N}\sum_{j^{\prime}=1}^{N}\left[\left(F_{j}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin\Theta_{j}^{(0)}(1+O(\sigma_{\Theta}^{2}))\right)\left(F_{j^{\prime}}^{(0)}+\mathrm{i}\frac{a\sigma_{F}^{2}}{\gamma_{\mathrm{r}}\gamma_{F}}\sin\Theta_{j^{\prime}}^{(0)}(1+O(\sigma_{\Theta}^{2}))\right)\right.$
$\displaystyle\qquad\qquad\cdot\mathrm{e}^{\mathrm{i}(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)})+\frac{2\pi\mathrm{i}}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)}\mathrm{e}^{-\frac{a\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{\Theta}}(F_{j}^{(0)}\cos\Theta_{j}^{(0)}-F_{j^{\prime}}^{(0)}\cos\Theta_{j^{\prime}}^{(0)})(1+O(\sigma_{\Theta}^{2}))}$
$\displaystyle\qquad\qquad\cdot\int_{0}^{\infty}e^{-\mathrm{i}\Omega^{\ast}_{\mathrm{r}}u-D_{\mathrm{r}}^{\ast}u-\frac{1}{2}V_{\Delta\Theta,\mathrm{G}}(u)}\,\,\mathrm{d}u$
$\displaystyle\left.\qquad\qquad\qquad\qquad+O\left(\int_{0}^{\gamma_{-}^{-1}}\left(1+\frac{a\bar{F}}{\gamma_{\mathrm{r}}\gamma_{-}}\right)\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]\,\mathrm{d}u\right)\right]$

Integrating by parts in the main integral, noting from Eq. (44) that
$V_{\Delta\Theta,\mathrm{G}}(u)$ is subdominant to $D_{\mathrm{r}}^{\ast}u$
for large $u$ and $V_{\Delta\Theta,\mathrm{G}}(0)=0$:

$\int_{0}^{\infty}e^{-\mathrm{i}\Omega^{\ast}_{\mathrm{r}}u-D_{\mathrm{r}}^{\ast}u-\frac{1}{2}V_{\Delta\Theta,\mathrm{G}}(u)}\,\mathrm{d}u=\frac{1}{\mathrm{i}\Omega^{\ast}_{\mathrm{r}}+D_{\mathrm{r}}^{\ast}}\left[1-\frac{1}{2}\int_{0}^{\infty}V_{\Delta\Theta,\mathrm{G}}^{\prime}(u)e^{-\mathrm{i}\Omega^{\ast}_{\mathrm{r}}u-D_{\mathrm{r}}^{\ast}u-V_{\Delta\Theta,\mathrm{G}}(u)}\,\mathrm{d}u\right].$

Applying the estimate on $V_{\Delta\Theta,\mathrm{G}}^{\prime}(u)$ from Eq.
(44), we have:

$\int_{0}^{\infty}e^{-\mathrm{i}\Omega^{\ast}_{\mathrm{r}}u-D_{\mathrm{r}}^{\ast}u-V_{\Delta\Theta,\mathrm{G}}(u)}\,\mathrm{d}u=\frac{1}{\mathrm{i}\Omega^{\ast}_{\mathrm{r}}+D_{\mathrm{r}}^{\ast}}\left[1+O\left(\frac{D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}}}{\gamma_{-}}\right)\right].$

For the correction term in Eq. (51), we estimate:

$\displaystyle\int_{0}^{\gamma_{-}^{-1}}\left(1+\frac{a\bar{F}}{\gamma_{\mathrm{r}}\gamma_{-}}\right)\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]\,\mathrm{d}u$
$\displaystyle\sim
O\left(\left(1+\frac{a\bar{F}}{\gamma_{\mathrm{r}}\gamma_{-}}\right)\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]\gamma_{-}^{-1}\right)$

Next computing the real part of the resulting complex expression in Eq. (51),
we obtain:

$\displaystyle\lim_{t\rightarrow\infty}\frac{\text{Cov}_{\mathrm{G}}(\mathbf{X}^{(\mathrm{c})}(t),\mathbf{X}^{(\mathrm{c})}(t))}{2t}=D_{\mathrm{t}}\mathit{\mathsf{I}}+\frac{e^{-\sigma_{\Theta}^{2}}}{2\gamma_{\mathrm{t}}^{2}(\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2})}\mathit{\mathsf{I}}$
(52)
$\displaystyle\qquad\qquad\cdot\left[\sum_{j,j^{\prime}=1}^{N}F_{j}^{(0)}\left(F_{j^{\prime}}^{(0)}(D_{\mathrm{r}}^{\ast}+O(\omega)))+\frac{2a\sigma_{F}^{2}\sin\Theta_{j^{\prime}}^{(0)}(\Omega^{\ast}_{\mathrm{r}}+O(\omega))}{\gamma_{\mathrm{r}}\gamma_{F}}\right)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\cdot\cosh\left(\frac{a\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{\Theta}}(F_{j}^{(0)}\cos(\Theta_{j}^{(0)})-F_{j^{\prime}}^{(0)}\cos(\Theta_{j^{\prime}}^{(0)}))\right)$
$\displaystyle\qquad\qquad\qquad\qquad\cdot\cos\left(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)}+\frac{2\pi}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)\right)$
$\displaystyle\qquad\qquad+\sum_{j,j^{\prime}=1}^{N}F_{j}^{(0)}\left(F_{j^{\prime}}^{(0)}(\Omega^{\ast}_{\mathrm{r}}+O(\omega))-\frac{2a\sigma_{F}^{2}\sin\Theta_{j^{\prime}}^{(0)}(D_{\mathrm{r}}^{\ast}+O(\omega))}{\gamma_{\mathrm{r}}\gamma_{F}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\cdot\sinh\left(\frac{a\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{\Theta}}(F_{j^{\prime}}^{(0)}\cos(\Theta_{j^{\prime}}^{(0)})-F_{j}^{(0)}\cos(\Theta_{j}^{(0)}))\right)$
$\displaystyle\qquad\qquad\qquad\qquad\left.\cdot\sin\left(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)}+\frac{2\pi}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)\right)\right]$
$\displaystyle\qquad\qquad\qquad\cdot\left[1+O\left(\frac{D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}}}{\gamma_{-}}\right)+O\left(\frac{a\bar{F}((\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2})\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{-}}\right)+O\left(\frac{a^{2}\sigma_{F}^{4}}{\gamma_{\mathrm{r}}^{2}\gamma_{F}^{2}\bar{F}^{2}}\right)\right]$
(53)

where

$\omega\equiv\left(1+\frac{a\bar{F}}{\gamma_{\mathrm{r}}\gamma_{-}}\right)\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]\frac{\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2}}{\gamma_{-}},$

a quantity with dimensions of rate, denotes the error arising from the
integration over times $\lesssim\gamma_{-}^{-1}$.

In $\omega$ we have separated out a nondimensional factor of
$\frac{\bar{F}a}{\gamma_{\mathrm{r}}\gamma_{-}}$ which from Eq. (16), is $\sim
20/N^{2}$ and thus order unity for the biophysical parameters in Table 1 and
reasonable colony sizes. But to allow for situations in which this factor may
be large and thus somewhat counteract the asymptotically small factors of
$\sigma_{\Theta}$ and $\sigma_{F}/\bar{F}$, we recast errors involving this
factor in terms of the effective rotational diffusion via the expresssion (10)
Taylor expanded with respect to these small factors:

$\displaystyle\frac{a\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{\Theta}}F_{j}^{(0)}$
$\displaystyle\sim
O(N^{-1/2}(D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}})^{1/2}\gamma_{+}^{1/2}\gamma_{-}^{-1}\sigma_{\Theta}),$
(54)
$\displaystyle\frac{a\sigma_{F}^{2}}{F_{j}^{(0)}\gamma_{\mathrm{r}}\gamma_{F}}$
$\displaystyle\sim
O(N^{-1/2}(D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}})^{1/2}\gamma_{+}^{1/2}\gamma_{-}^{-1}(\sigma_{F}/\bar{F})),$
$\displaystyle\frac{a\bar{F}}{\gamma_{\mathrm{r}}\gamma_{-}}\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]$
$\displaystyle\sim
O\left(N^{-1/2}(D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}})^{1/2}\gamma_{+}^{1/2}\gamma_{-}^{-1}(\sigma_{\Theta}+(\sigma_{F}/\bar{F}))\right).$

Applying now the assumption (25) that $\gamma_{+}/\gamma_{-}$ is order unity,
we see that the squares of the first two expressions in Eq. (54) are small
compared to the $O((D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}})\gamma_{-}^{-1})$
integration error, and the functions of them in Eq. (52) should be self-
consistently represented by their first order Taylor approximations. The last
expression implies that

$\omega\sim
O\left(\left(\left[(\sigma_{F}/\bar{F})^{2}+\sigma_{\Theta}^{2}\right]+\left(\frac{D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}}}{N\gamma_{-}}\right)^{1/2}\left[(\sigma_{F}/\bar{F})+\sigma_{\Theta}\right]\right)\frac{\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2}}{\gamma_{-}}\right)$

With these simplifications and reductions we arrive at:

$\displaystyle\lim_{t\rightarrow\infty}\frac{\text{Cov}_{\mathrm{G}}(\mathbf{X}^{(\mathrm{c})}(t),\mathbf{X}^{(\mathrm{c})}(t))}{2t}=D_{\mathrm{t}}\mathit{\mathsf{I}}+\frac{e^{-\sigma_{\Theta}^{2}}}{2\gamma_{\mathrm{t}}^{2}(\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2})}\mathit{\mathsf{I}}$
(55)
$\displaystyle\qquad\qquad\cdot\left\\{\left[\sum_{j,j^{\prime}=1}^{N}F_{j}^{(0)}\left(F_{j^{\prime}}^{(0)}(D_{\mathrm{r}}^{\ast}+O(\omega)))+\frac{2a\sigma_{F}^{2}\sin\Theta_{j^{\prime}}^{(0)}(\Omega^{\ast}_{\mathrm{r}}+O(\omega))}{\gamma_{\mathrm{r}}\gamma_{F}}\right)\right.\right.$
$\displaystyle\qquad\qquad\qquad\qquad\cdot\cos\left(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)}+\frac{2\pi}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)\right)$
$\displaystyle\qquad\qquad+\frac{a(\Omega^{\ast}_{\mathrm{r}}+O(\omega))\sigma_{\Theta}^{2}}{\gamma_{\mathrm{r}}\gamma_{\Theta}}\sum_{j,j^{\prime}=1}^{N}F_{j}^{(0)}F_{j^{\prime}}^{(0)}(F_{j^{\prime}}^{(0)}\cos(\Theta_{j^{\prime}}^{(0)})-F_{j}^{(0)}\cos(\Theta_{j}^{(0)}))$
$\displaystyle\qquad\qquad\qquad\qquad\left.\cdot\sin\left(\Theta_{j}^{(0)}-\Theta_{j^{\prime}}^{(0)}+\frac{2\pi}{N}\left(j-j^{\prime}+\frac{S_{j}-S_{j^{\prime}}}{l}\right)\right)\right]$
$\displaystyle\qquad\qquad\qquad\cdot\left[1+O\left(\frac{D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}}}{\gamma_{-}}\right)+O\left(\left(\frac{D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}}}{N\gamma_{-}}\right)^{1/2}\left[(\sigma_{F}/\bar{F})+\sigma_{\Theta}\right]\sigma_{\Theta}^{2}\right)\right]$
(56)

The error terms are all formally second order in the small parameters
$\sigma_{F}/\bar{F}$ and $\sigma_{\Theta}$, but we’ve also noted when these
small parameters are therefore better expressed in terms of the effective
rotational diffusivity, which Eq. (12) shows is itself a second order quantity
in these small parameters. We now explain why we retain some second order
terms in explicit form. First of all, we note that $\omega$ has, beyond the
second order factors of smallness in dynamical flagellar variation, an
additional factor of
$(\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2})/\gamma_{-}$
which is small compared to the larger of $\Omega^{\ast}_{\mathrm{r}}$ and
$D_{\mathrm{r}}^{\ast}$ under the self-consistency assumption (24) that the
effective colony dynamics are slow relative to the flagellar dynamics.
Secondly, the explicitly presented second order terms are proportional to
$\Omega^{\ast}_{\mathrm{r}}$, which we have seen in Subsection 3.1 to
typically be larger than $D_{\mathrm{r}}^{\ast}$, particularly for larger
colonies. So even though their ratio to $\Omega^{\ast}_{\mathrm{r}}$ is
$O(N^{-1}(D_{\mathrm{r}}^{\ast}-D_{\mathrm{r}})^{1/2}((\sigma_{F}/\bar{F})+\sigma_{\Theta}))$
and thus negligible to the level of our approximation, their ratio to
$D_{\mathrm{r}}^{\ast}$ won’t be as small.

In the summary discussion in Subsection 3.2.2, we report the translational
diffusivity based on the expression (55), dropping the error terms. As
discussed at the beginning of this calcluation, we did not attempt to quantify
the error in computing averages as if the various quantities are jointly
Gaussian. Informal considerations would suggest these errors would arise at
fourth order, which would be negligible with respect to the results of our
more precise estimates within the Gaussian approximation framework.

##  5 Generalized Colony Geometries

The disclike model from Section 2 may be an appropriate simplified model for a
rosette colony, but choanoflagettes are known also to form chain colonies
(Dayel et al., 2011). We show in this section how much of our analysis of the
mobility of disclike colonies can be extended to more general thin two-
dimensional colony shapes. We do continue to consider the colonies as rigid,
which is reasonable even for choanoflagellate chain colonies due to their
intercellular bridges (Roper et al., 2013). The flagellar force model remains
of the same Langevin type as in Subsection 2.2; we only generalize in
Subsection 5.1 the representation of the relaxed locations and orientations of
the flagella and the hydrodynamic resistance of the colony. Within this
generalized model, we present mobility statistics for the colony in Subsection
5.2 under the general assumption of small dynamical fluctuations of the
flagellar force magnitudes and angles. We do not present demographic mobility
statistics, as these would be too complex to develop in a general geometric
setting. In Subsection 5.3, we sketch how these results can be derived using
the methods of Section 4.

###  5.1 Model

We first appeal to low Reynolds number hydrodynamics to define the resistance
supermatrix (Happel and Brenner, 1965, Sec. 5.4) for the joint translation and
rotational dynamics of the rigid two-dimensional colony within its own plane:

$\mathit{\mathsf{\Gamma}}=\begin{pmatrix}\gamma_{\mathrm{t}1}&0&\gamma_{\mathrm{tr}1}\cr
0&\gamma_{\mathrm{t}2}&\gamma_{\mathrm{tr}2}\cr\gamma_{\mathrm{tr}1}&\gamma_{\mathrm{tr}2}&\gamma_{\mathrm{r}}\end{pmatrix}$
(57)

We have chosen a body frame coordinate system with the Cartesian axes aligned
with the principal axes of the translational component of the friction matrix
in the upper left $2\times 2$ block. If the translational friction
coefficients $\gamma_{\mathrm{t}1}\leq\gamma_{\mathrm{t}2}$ along the two axes
are unequal, we take the convention of choosing the first Cartesian axis to
lie along the direction with a lower friction coefficient. If the colony had a
thin ellipsoidal shape, this would be equivalent to choosing the body frame
coordinates based on the major and minor axes of the elliptical cross section.
The origin of the body frame is chosen to be the center of mass. Then
$\gamma_{\mathrm{r}}$ is the rotational drag coefficient, and
$\gamma_{\mathrm{tr}1}$, $\gamma_{\mathrm{tr}2}$ are the frictional coupling
coefficients between translational motion along the principal body axes and
rotation.

For each of the $j=1,\ldots,N$ flagella, we specify the two-dimensional
position $\mathbf{B}_{j}$ of the base of each flagellum, a unit vector
$\hat{\mathbf{e}}_{j}$ describing the mean direction of force applied by that
flagellum, and a force magnitude $F_{j}^{(0)}$ of that flagellum. A schematic
for a chain colony shape is given in Figure 10. For the disc colony model from
Subsection 2.3, the body frame orientation is arbitrary,

$\mathbf{B}_{j}=a\begin{pmatrix}\cos(2\pi[\frac{j-\frac{1}{2}+\frac{S_{j}}{l}}{N}])\\\
\sin(2\pi[\frac{j-\frac{1}{2}+\frac{S_{j}}{l}}{N}]\end{pmatrix},\hat{\mathbf{e}}_{j}=\begin{pmatrix}-\cos\left(2\pi[\frac{j-\frac{1}{2}+\frac{S_{j}}{l}}{N}]+\Theta_{j}^{(0)}\right)\\\
-\sin\left(2\pi[\frac{j-\frac{1}{2}+\frac{S_{j}}{l}}{N}]+\Theta_{j}^{(0)}\right)\end{pmatrix},$
(58)

Figure 10: Model Schematic: A colony of five cells arranged in a chain with
with colony orientation $\Theta^{(c)}$ describing the rotation, about the
center of mass at $\mathbf{C}$, of the major axis of the hydrodynamic friction
matrix from the horizontal $\mathbf{B}_{j}$ is the location of the juncture of
flagellum $j$ with the colony surface, and the relaxed orientation of its
force is $\hat{\mathbf{e}}_{j}$. $\Theta_{j}$ is the angle between the
momentary flagellar force orientation relative to its relaxed orientation.

The force summation rules of the generalized colony in the body frame are:

$\displaystyle\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}}$
$\displaystyle=\sum_{j=1}^{N}F_{j}(t)\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j},$
(59) $\displaystyle T^{(\mathrm{c})}$
$\displaystyle=\sum_{j=1}^{N}F_{j}(t)(\mathbf{B}_{j}-\mathbf{C})\cdot\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}^{\perp}$

where $\Theta_{j}(t)$ is now the force orientation for flagellum $j$ measured
with respect to its natural orientation $\hat{\mathbf{e}}_{j}$. For the
disclike case considered in the previous sections, we introduced an angle
$\Theta_{j}^{(0)}$ for the angle between the flagellum’s natural orientation
$\hat{\mathbf{e}}_{j}$ and the surface normal. For an asymmetric colony, there
is no natural reference angular position so there is no point in introducing
such a demographic parameter; it is embedded in $\hat{\mathbf{e}}_{j}$ (as can
be seen for the disclike case in Eq. (58)). Thus, we carry over the flagellar
dynamic models from Subsection 2.2 by simply setting $\Theta_{j}^{(0)}\equiv
0$ in Eq. (2) for the generalized geometry model under current discussion.

###  5.2 Colony Transport Statistics

In this subsection, without recurrent comment, all results are presented only
to second order in the force fluctuation relative to its mean value
$\sigma_{F}/\bar{F}$ and angular fluctuation magnitude $\sigma_{\Theta}$, both
assumed small. We express the results with $\sim$ to indicate they are
asymptotic approximations in this sense. We find it useful to represent some
results using the rotation matrix

$\mathit{\mathsf{R}}_{\theta}\equiv\begin{pmatrix}\cos\theta&-\sin\theta\cr\sin\theta&\cos\theta\end{pmatrix}.$

and the perpendicular direction
$\hat{\mathbf{e}}_{j}^{\perp}=\mathit{\mathsf{R}}_{-\pi/2}\hat{\mathbf{e}}_{j}$
to the orientation of flagellum $j$, directed to give the appropriate sign for
the torque in Eq. (59). The colony transport statistics are naturally
expressed in terms of sums of quantities associated with the geometric
orientation of each flagellum $j$ in relation to the colony body:

  * •

the vectors (with units of time per mass) describing the ratio of the colony
velocity induced per force directed along the flagellum’s relaxed orientation
$\hat{\mathbf{e}}_{j}$:

$\vec{\nu}_{j}\equiv\left[\mathit{\mathsf{\Gamma}}^{-1}\begin{pmatrix}\hat{\mathbf{e}}_{j}\cr(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}^{\perp}\end{pmatrix}\right]_{\mathbf{x}}$

and per unit force directed transversely (along
$\hat{\mathbf{e}}_{j}^{\perp}$) to the flagellum’s relaxed orientation:

$\vec{\nu}_{\perp,j}\equiv\left[\mathit{\mathsf{\Gamma}}^{-1}\begin{pmatrix}\hat{\mathbf{e}}_{j}^{\perp}\cr-(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}\end{pmatrix}\right]_{\mathbf{x}}$

with $[\cdot]_{\mathbf{x}}$ denoting the projection of the three-dimensional
supermatrix onto the translational components (first two indices).

  * •

the scalar (with units of time per mass per length) describing the ratio of
the the colony angular velocity induced per unit force directed along the
flagellum’s relaxed orientation $\hat{\mathbf{e}}_{j}$:

$\rho_{j}\equiv\left[\mathit{\mathsf{\Gamma}}^{-1}\begin{pmatrix}\hat{\mathbf{e}}_{j}\cr(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}^{\perp}\end{pmatrix}\right]_{\theta}$

and per unit force directed transversely (along
$\hat{\mathbf{e}}_{j}^{\perp}$) to the flagellum’s relaxed orientation:

$\rho_{\perp,j}\equiv\left[\mathit{\mathsf{\Gamma}}^{-1}\begin{pmatrix}\hat{\mathbf{e}}_{j}^{\perp}\cr-(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}\end{pmatrix}\right]_{\theta}$

The results naturally are not dependent on simultaneous sign reversal of
$\vec{\nu}_{\perp,j}$ and $\rho_{\perp,j}$, corresponding to choosing the
opposite “transverse” direction. For the disc configuration developed in
Section 2, we’d have
$\vec{\nu}_{j}=\gamma_{\mathrm{t}}^{-1}\hat{\mathbf{e}}_{j}$,
$\vec{\nu}_{\perp,j}=\gamma_{\mathrm{t}}^{-1}\hat{\mathbf{e}}_{j}^{\perp}$,
$\rho_{j}=-a\gamma_{\mathrm{r}}^{-1}\sin\Theta_{j}^{(0)}$, and
$\rho_{\perp,j}=a\gamma_{\mathrm{r}}^{-1}\cos\Theta_{j}^{(0)}$.

The fundamental transport statistics for a colony under our flagellar forcing
model can then be expressed in terms of the geometric configuration:

  * •

Rotational Drift:

$\Omega^{\ast}_{\mathrm{r}}=\sum_{j=1}^{N}F_{j}^{(0)}\mathrm{e}^{-\frac{1}{2}\sigma_{\Theta}^{2}}\rho_{j}$
(60)

where $[\cdot]_{\theta}$ the angular (third index) component of the three-
dimensional supervector in the brackets.

  * •

Effective Rotational Diffusion:

$\displaystyle D_{\mathrm{r}}^{\ast}$ $\displaystyle\sim
D_{\mathrm{r}}+\sigma_{F}^{2}\gamma_{F}^{-1}\sum_{j=1}^{N}\rho_{j}^{2}+\sigma_{\Theta}^{2}\gamma_{\Theta}^{-1}\sum_{j=1}^{N}(F_{j}^{(0)})^{2}\rho_{\perp,j}^{2}$
(61)

The first term is the thermal contribution

$D_{\mathrm{r}}=\frac{k_{B}T}{\gamma_{\mathrm{r}}-\gamma_{\mathrm{tr}1}^{2}/\gamma_{\mathrm{t}1}-\gamma_{\mathrm{tr}2}^{2}/\gamma_{\mathrm{t}2}}$

while the second term arises from the active flagellar force fluctuations and
the third term from the orientational fluctuations in the active flagellar
forces.

  * •

Mean-Square Instantaneous Speed:

$\displaystyle\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle$
$\displaystyle\sim\mathrm{e}^{-\sigma_{\Theta}^{2}}\left|\sum_{j=1}^{N}F_{j}^{(0)}\vec{\nu}_{j}\right|^{2}+\sigma_{F}^{2}\sum_{j=1}^{N}\left|\vec{\nu}_{j}\right|^{2}+\sigma_{\Theta}^{2}\sum_{j=1}^{N}(F_{j}^{(0)})^{2}\left|\vec{\nu}_{\perp,j}\right|^{2}$
(62)

The first term is the squared speed when all fluctuations are neglected,
multiplied by $\mathrm{e}^{-\sigma_{\Theta}^{2}}$ to account for the
mitigation in the mean flagellar force applied due to angular fluctuations.
The second and third term are contributions measuring the magnitudes of
velocities induced by rapid dynamical flagellar fluctuations, and we show next
that they disappear when viewed over a time scale long compared to these
fluctuations.

  * •

Mean-Square Speed Over Coarse Time Intervals The mean-square displacement
exhibits a ballistic regime over time scales long compared to the flagellar
time scales but short compared to the rotational dynamics of colony:

$\displaystyle\mathbb{E}|\mathbf{X}^{(\mathrm{c})}(t+\tau)-\mathbf{X}^{(\mathrm{c})}(t)|^{2}$
$\displaystyle\sim 4D_{\mathrm{t}}\tau+V^{\ast}{}^{2}\tau^{2}\text{ for
}\gamma_{F}^{-1},\gamma_{\Theta}^{-1}\ll\tau\ll\Omega^{\ast}_{\mathrm{r}}{}^{-1},D_{\mathrm{r}}^{\ast}{}^{-1},$
$\displaystyle V^{\ast}{}^{2}$
$\displaystyle\sim\mathrm{e}^{-\sigma_{\Theta}^{2}}\left|\sum_{j=1}^{N}F_{j}^{(0)}\vec{\nu}_{j}\right|^{2}$
(63)

  * •

Effective Translational Diffusion: The colony will have no sustained drift due
to statistical isotropy, so its mean position will remain bounded. The long-
time diffusivity matrix will similarly be a simple multiple of the identity:

$\lim_{t\rightarrow\infty}\frac{\text{Cov}(\mathbf{X}^{(\mathrm{c})}(t),\mathbf{X}^{(\mathrm{c})}(t))}{2t}=D_{\mathrm{t}}^{\ast}\mathit{\mathsf{I}}$

with

$\displaystyle D_{\mathrm{t}}^{\ast}\sim
D_{\mathrm{t}}+\frac{D_{\mathrm{r}}^{\ast}V^{\ast}{}^{2}}{2(D_{\mathrm{r}}^{\ast}{}^{2}+\Omega^{\ast}_{\mathrm{r}}{}^{2})}-\frac{\Omega^{\ast}_{\mathrm{r}}k_{B}T}{(\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2})}\frac{\sigma_{F}^{2}}{\gamma_{F}}\sum_{j=1}^{N}F_{j}^{(0)}\vec{\nu}_{j}\cdot\left[\mathit{\mathsf{\Gamma}}^{-1}\right]_{\mathbf{x},\theta}$
(64)
$\displaystyle\qquad-\frac{\Omega^{\ast}_{\mathrm{r}}}{(\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2})}\frac{\sigma_{F}^{2}}{\gamma_{F}}\sum_{j=1}^{N}F_{j}^{(0)}\vec{\nu}_{j}\cdot\sum_{j^{\prime}=1}^{N}\rho_{j^{\prime}}\vec{\nu}_{j^{\prime}}-\frac{\Omega^{\ast}_{\mathrm{r}}}{(\Omega^{\ast}_{\mathrm{r}}{}^{2}+D_{\mathrm{r}}^{\ast}{}^{2})}\frac{\sigma_{\Theta}^{2}}{\gamma_{\Theta}}\sum_{j=1}^{N}F_{j}^{(0)}\vec{\nu}_{j}\cdot\sum_{j^{\prime}=1}^{N}(F_{j^{\prime}}^{(0)})^{2}\rho_{\perp,j^{\prime}}\vec{\nu}_{\perp,j^{\prime}}$

where $D_{\mathrm{t}}$ is the thermal diffusivity, $V^{\ast}{}^{2}$ is the
mean-square velocity observed on scales long compared to flagellar time scales
(63), and $[\mathit{\mathsf{\Gamma}}^{-1}]_{\mathbf{x},\theta}$ is the vector
describing the translational-rotational coupling of the inverse resistance
matrix in the body frame of the colony. This translational diffusivity result
requires the additional assumption (24) that the effective rotational
diffusivity is small compared with the flagellar fluctuation rates together
with the assumption (25) that the relaxation rates of flagellar force and
orientation are comparable.

The first two terms in Eq. (64) correspond to an isotropic particle moving at
constant speed $V^{\ast}{}^{2}$ in a direction governed by rotation at rate
$\Omega^{\ast}_{\mathrm{r}}$ together with rotational diffusion
$D_{\mathrm{r}}^{\ast}$. For the disc-like case, our argument for the
subdominance of the terms beyond the first two in Eq. (64) relied on
relationships derived in Subsection 3.1.2 between the effective rotational
drift and diffusivity for a nearly symmetric disc arrangement. These cannot be
expected to hold for general geometries, and those formally second order terms
might be significant when the rotational velocity $\Omega^{\ast}_{\mathrm{r}}$
is substantially greater than the effective rotational diffusivity
$D_{\mathrm{r}}^{\ast}$ of the colony.

All these results can be verified to agree with the results obtained for the
disc-like model when specialized to that geometry. Computing demographic
averages of these quantities analytically is considerably more challenging
than for the disc-like model from Section 2.1. For example, demographic
variability across colonies will affect, in a correlated way, the resistance
matrix $\mathit{\mathsf{\Gamma}}$ in (57) and other geometric features of the
colony which are multiplied together in all mobility expressions. But these
transport quantities are easily computed (without the need for running long
Monte Carlo simulations) for any given realization of a colony’s physical
makeup, and thus we can easily generate a statistical distribution of the
mobility statistics for a specified demographic colony distribution model by
direct sampling, just as we did for translational diffusivity in the disclike
colony model (blue circles in Figure 9).

###  5.3 Justification of Generalized Geometry Formulas

The derivation of the formulas for the effective translation and rotation of
colonies with a general configuration follows the same mathematical procedure
as for the disc shape from Section 4, just with more notationally complex
versions of the formulas. We will only discuss some key components of the
argument and add details only where they are novel to the generalized colony
geometry setting. All calculations are only accurate to second order in
$\sigma_{F}/\bar{F}$ and $\bar{\Theta}$; we do not report formal error
estimates here.

We first compute the first and second order statistics of the total body force
and torque on the colony:

  * •

the mean total force on the colony, in the body frame:

$\langle\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}}\rangle=\sum_{j=1}^{N}F_{j}^{(0)}\hat{\mathbf{e}}_{j}\mathrm{e}^{-\frac{1}{2}\sigma_{\Theta}^{2}}$
(65)

  * •

the mean torque on the colony:

$\langle
T^{(\mathrm{c})}\rangle=\sum_{j=1}^{N}F_{j}^{(0)}\mathrm{e}^{-\frac{1}{2}\sigma_{\Theta}^{2}}(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}^{\perp}$
(66)

  * •

the covariance matrix of the total force on the colony, in the body frame

$\text{Cov}(\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}},\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}})=\sigma_{F}^{2}\sum_{j=1}^{N}\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}+\sigma_{\Theta}^{2}\sum_{j=1}^{N}(F_{j}^{(0)})^{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}^{\perp}$

  * •

the variance of the torque on the colony:

$\text{Var}(T^{(\mathrm{c})})\sim\sum_{j=1}^{N}(F_{j}^{(0)})^{2}\sigma_{\Theta}^{2}\left|(\mathbf{B}_{j}-\mathbf{C})\bm{\cdot}\hat{\mathbf{e}}_{j}\right|^{2}+\sigma_{F}^{2}\left|(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}^{\perp}\right|^{2}$

  * •

the covariance of the torque and vectorial colony force in the body frame

$\text{Cov}(T^{(\mathrm{c})},\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}})=\sum_{j=1}^{N}\sigma_{F}^{2}\left[(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}^{\perp}\right]\hat{\mathbf{e}}_{j}-(F_{j}^{(0)})^{2}\sigma_{\Theta}^{2}\left[(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}\right]\hat{\mathbf{e}}_{j}^{\perp})$
(67)

These calculations are quite similar to those described in Subsection 4.1,
averaging over the stationary distribution of the flagellar variables
$\\{\Theta_{j}(t)\\}_{j=1}^{N}$, $\\{F_{j}(t)\\}_{j=1}^{N}$. The main
complication is the need to introduce a more general friction matrix (57),
which moreover couples rotational and translational motion in general (Happel
and Brenner, 1965, Sec. 5.4). As a consequence, we need to account for
correlations between rotations and translations in computing the effective
rotational diffusivity. The covariance (67) was obtained by the following
calculation together with an expansion with respect to small
$\sigma_{\Theta}$:

$\displaystyle\text{Cov}(T^{(\mathrm{c})},\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}})$
$\displaystyle=\sum_{j=1}^{N}\sigma_{F}^{2}(\mathbf{B}_{j}-\mathbf{C})\cdot\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}^{\perp}\otimes\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}\rangle$
$\displaystyle\qquad\qquad+(F_{j}^{(0)})^{2}(\mathbf{B}_{j}-\mathbf{C})\cdot\text{Cov}(\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}^{\perp},\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}),$
$\displaystyle\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}^{\perp}\otimes\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}\rangle$
$\displaystyle=\langle\cos^{2}\Theta_{j}(t)\rangle\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}-\langle\sin^{2}\Theta_{j}(t)\rangle(\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}^{\perp})$
$\displaystyle\qquad\qquad+\langle\sin\Theta_{j}(t)\cos\Theta_{j}(t)\rangle(\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}-\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}^{\perp})$
$\displaystyle=\frac{1+\langle\cos
2\Theta_{j}(t)\rangle}{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}-\frac{1-\langle\cos
2\Theta_{j}(t)\rangle}{2}\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}^{\perp}+0(\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}-\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}^{\perp})$
$\displaystyle=\frac{1+\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}-\frac{1-\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}^{\perp}$
$\displaystyle\text{Cov}(\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}^{\perp},\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j})$
$\displaystyle=\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}^{\perp}\otimes\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}\rangle-(\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\rangle\hat{\mathbf{e}}_{j}^{\perp})\otimes(\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\rangle\hat{\mathbf{e}}_{j})$
$\displaystyle=\frac{1+\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}-\frac{1-\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}(\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}^{\perp})-(\mathrm{e}^{-\sigma_{\Theta}^{2}/2}\hat{\mathbf{e}}_{j}^{\perp})\otimes(\mathrm{e}^{-\sigma_{\Theta}^{2}/2}\hat{\mathbf{e}}_{j})$
$\displaystyle=\frac{(1-\mathrm{e}^{-\sigma_{\Theta}^{2}})^{2}}{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}-\frac{1-\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}(\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}^{\perp})$

Similarly, to compute the covariance matrix of the body force (which
contributes to the enhanced rotational diffusion through the rotational-
translational coupling)

$\displaystyle\text{Cov}(\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}},\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}})$
$\displaystyle=\sum_{j=1}^{N}\sigma_{F}^{2}\langle(\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j})\otimes(\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j})\rangle+\sum_{j=1}^{N}(F_{j}^{(0)})^{2}\text{Cov}(\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j},\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}),$
$\displaystyle\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}\otimes\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}\rangle$
$\displaystyle=\langle\cos^{2}\Theta_{j}(t)\rangle\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}+\langle\sin^{2}\Theta_{j}(t)\rangle(\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}^{\perp})$
$\displaystyle\qquad\qquad-\langle\sin\Theta_{j}(t)\cos\Theta_{j}(t)\rangle(\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}^{\perp}+\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j})$
$\displaystyle=\frac{1+\langle\cos
2\Theta_{j}(t)\rangle}{2}\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}+\frac{1-\langle\cos
2\Theta_{j}(t)\rangle}{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}^{\perp}+0(\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}^{\perp}+\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j})$
$\displaystyle=\frac{1+\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}+\frac{1-\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}^{\perp},$
$\displaystyle\text{Cov}(\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j},\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j})$
$\displaystyle=\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}\otimes\mathit{\mathsf{R}}_{\Theta_{j}(t)}\hat{\mathbf{e}}_{j}\rangle-(\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\rangle\hat{\mathbf{e}}_{j})\otimes(\langle\mathit{\mathsf{R}}_{\Theta_{j}(t)}\rangle\hat{\mathbf{e}}_{j})$
$\displaystyle=\frac{1+\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}+\frac{1-\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}^{\perp}-(\mathrm{e}^{-\sigma_{\Theta}^{2}/2}\hat{\mathbf{e}}_{j})\otimes(\mathrm{e}^{-\sigma_{\Theta}^{2}/2}\hat{\mathbf{e}}_{j})$
$\displaystyle=\frac{(1-\mathrm{e}^{-\sigma_{\Theta}^{2}})^{2}}{2}\hat{\mathbf{e}}_{j}\otimes\hat{\mathbf{e}}_{j}+\frac{1-\mathrm{e}^{-2\sigma_{\Theta}^{2}}}{2}\hat{\mathbf{e}}_{j}^{\perp}\otimes\hat{\mathbf{e}}_{j}^{\perp}$

The instantaneous mean-square speed (62) follows from the second-order
statistics via the relation:

$\mathbf{V}^{(\mathrm{c})}(t)=\mathit{\mathsf{R}}_{\Theta^{(\mathrm{c})}(t)}\left[\mathit{\mathsf{\Gamma}}^{-1}\begin{pmatrix}\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}}(t)\cr
T^{(\mathrm{c})}(t)\end{pmatrix}\right]_{\mathbf{x}}$ (68)

which implies:

$\langle|\mathbf{V}^{(\mathrm{c})}|^{2}\rangle=\left|\left[\mathit{\mathsf{\Gamma}}^{-1}\begin{pmatrix}\langle\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}}\rangle\cr\langle
T^{(\mathrm{c})}\rangle\end{pmatrix}\right]_{\mathbf{x}}\right|^{2}+\operatorname{Tr}\left[\mathit{\mathsf{\Gamma}}^{-1}\begin{pmatrix}\text{Cov}(\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}},\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}})&\text{Cov}(\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}},T^{(\mathrm{c})})\cr\text{Cov}(T^{(\mathrm{c})},\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}})&\text{Var}(T^{(\mathrm{c})})\end{pmatrix}\mathit{\mathsf{\Gamma}}^{-1}\right]_{\mathbf{x}},$

The matrix appearing here takes the form of a sum of dyads from each flagellum
which allows the somewhat simplified representation in Eq. (62).

To compute the long-run drift and diffusion, we turn to the equations of
microhydrodynamics which balance the active forces and torques, drag force and
torque, and the Brownian forces and torques, enhanced by the active flagellar
noise. The flagellar force fluctuations (which contribute to the terms
proportional to $\sigma_{F}^{2}$) decay exponentially with time scale
$\gamma_{F}^{-1}$, and the force orientation fluctuations decay exponentially
with time scale $\gamma_{\Theta}^{-1}$. By standard results on white-noise
approximations of Ornstein-Uhlenbeck noise (Gardiner, 2009, Sec. 8.1), the
square of the coefficient of the white noise should be chosen as twice the
product of the variance and correlation time of the Ornstein-Uhlenbeck process
to have an equivalent effect on slower variables. Indeed the derivation of
this result relies on the formula (34). Applying this coarse-graining in
conjunction with the single-time statistics of the colony force and torque
listed above, the microhydrodynamics equations for the colony in its body
frame reads (Brady and Bossis, 1988)

$\displaystyle\begin{pmatrix}\langle\mathbf{F}^{(\mathrm{c})}_{\mathrm{b}}\rangle\cr\langle
T^{(\mathrm{c})}\rangle\end{pmatrix}+\sqrt{2\mathit{\mathsf{\Sigma}}}\frac{\mathrm{d}\mathbf{W}^{(\mathrm{A})}}{\mathrm{d}t}+\sqrt{2k_{B}T\mathit{\mathsf{\Gamma}}}\frac{\mathrm{d}\mathbf{W}^{(\mathrm{T})}(t)}{\mathrm{d}t}$
$\displaystyle=\mathit{\mathsf{\Gamma}}\begin{pmatrix}\mathbf{V}^{(\mathrm{c)}}_{\mathrm{b}}\cr\Omega\end{pmatrix},$
(69)

We have represented generic three-dimensional white noise representing thermal
fluctuations by the formal expression
$\frac{\mathrm{d}\mathbf{W}^{(\mathrm{T})}(t)}{\mathrm{d}t}$ to facilitate the
re-expression of the microhydrodynamic equations in the language of stochastic
differential equations for the colony position and orientation in the lab
frame below. We have another independent copy of white noise
$\frac{\mathrm{d}\mathbf{W}^{(\mathrm{T})}(t)}{\mathrm{d}t}$ representing the
coarse-graining of the active flagellar fluctuations on the colony’s body
force and torque; the active noise matrix is given by

$\displaystyle\mathit{\mathsf{\Sigma}}$
$\displaystyle=\sigma_{F}^{2}\gamma_{F}^{-1}\sum_{j=1}^{N}\begin{pmatrix}\hat{\mathbf{e}}_{j}\cr(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}^{\perp}\end{pmatrix}\begin{pmatrix}\hat{\mathbf{e}}_{j}&(\mathbf{B}_{j}-\mathbf{C})\cdot\hat{\mathbf{e}}_{j}^{\perp}\end{pmatrix}$