# Slender body theories for rotating filaments

Ondrej Maxian and Aleksandar Donev

###### Abstract

Slender fibers are ubiquitous in biology, physics, and engineering, with
prominent examples including bacterial flagella and cytoskeletal fibers. In
this setting, slender body theories (SBTs), which give the resistance on the
fiber asymptotically in its slenderness $\epsilon$, are useful tools for both
analysis and computations. However, a difficulty arises when accounting for
twist and cross-sectional rotation: because the angular velocity of a filament
can vary depending on the order of magnitude of the applied torque, asymptotic
theories must give accurate results for rotational dynamics over a range of
angular velocities. In this paper, we first survey the challenges in applying
existing SBTs, which are based on either singularity or full boundary integral
representations, to rotating filaments, showing in particular that they fail
to consistently treat rotation-translation coupling in curved filaments. We
then provide an alternative approach which uses the Rotne-Prager-Yamakawa
tensor as a regularized singularity to represent the fiber centerline. Unlike
existing SBTs, this approach gives a grand mobility with symmetric rotation-
translation and translation-rotation coupling, making it applicable to a broad
range of angular velocities. To restore fidelity to the three-dimensional
filament geometry, we use our regularized singularity model to inform a simple
empirical equation which relates the mean force and torque along the filament
centerline to the translational and rotational velocity of the cross section.
The single unknown coefficient in the model is estimated numerically from
three-dimensional boundary integral calculations on a rotating, curved
filament.

##  1 Introduction

Whether an organism has one or millions of cells, slender filaments play a
vital role in its motion and maintenance. For instance, these types of
filaments make up the flagellar appendages that bacteria use to propel
themselves in run and tumble motions [4, 23]. In animal cells, slender
filaments are a key component in the cell cytoskeleton, which provides a
scaffolding for the cell [1], and adapts to promote cellular division,
migration, and resistance to deformation [11].

Given the vital role of slender filaments in biology, it is not surprising
that their interaction with the fluid medium in which they are immersed has
been a subject of study in physics and applied mathematics. Because most of
the applications of interest involve small lengthscales, the relevant fluid
equations are the Stokes equations, which are linear and therefore possible to
solve analytically through the use of the Stokes Green’s function. However,
because the body geometry is slender, it is expensive and often intractable to
solve the Stokes equations on the surface of the slender filament. This
intractability led Batchelor to develop the first so-called “slender body
theory” (SBT) [3], in which the body is approximated via a one-dimensional
line of Stokeslet singularities. Subsequent work by Keller and Rubinow [19]
and Johnson [15] placed Batchelor’s work on more rigorous asymptotic footing;
in particular, Johnson [15, 16] showed that adding a doublet to the line
integral produces a velocity that is guaranteed to be constant to
$\mathcal{O}(\epsilon)$ on the fiber cross section, where $\epsilon$ is the
filament slenderness.

Recent work has focused on placing these theories in the context of three-
dimensional well-posed partial differential equations (PDEs) [33, 31, 32, 22].
In particular, it has been shown that SBT can be derived from a three-
dimensional boundary integral equation [22, 13], and that its solution is an
$\mathcal{O}(\epsilon)$ approximation to a Stokes PDE with mixed Dirichlet-
Neumann boundary data and a boundary condition that the fiber maintain the
integrity of its cross section [33, 31, 32]. A second class of recent work has
developed SBTs for different geometries, including ribbons [21], cables [20],
and bodies with non-circular cross section [5].

Despite the breadth of literature on SBT and the novel theories for complex
geometries, there remains an open problem which is simple to formulate, but
difficult to solve, having to do with the _rotation_ of a filament. Consider a
filament with one clamped end and one free end, where a motor spins the
clamped end with frequency $\omega$. At some critical frequency $\omega_{c}$,
the applied torque overcomes the bending moment, the steady twirling state
becomes unstable, and large-scale filament deformations are introduced. This
instability, which has been studied extensively in [49, 25, 45, 24, 30] has
applications in the motion of bacterial flagella [7, 50] and supercoiling
during DNA transcription [26]. It has been shown that the critical frquency
$\omega_{c}$ scales like $1/\epsilon^{2}$, and thus that the torque required
to induce it is approximately independent of the fiber aspect ratio [38, Eq.
(62)]. Because of this, fibers being driven at the critical frequency have
angular velocities which are sufficiently large as to also affect the
translational dynamics. But as of yet, we are not aware of any slender body
theory that can consistently account for the full hydrodynamics of a slender,
translating, _and_ rotating filament, and in particular for the coupling
between the two. In fact, we recently showed that the failure to account for
rotation-translation coupling causes the SBT-based analysis of [49] to
overestimate $\omega_{c}$ by as much as 20% [30, Sec. 5.3]. To understand the
reason for this, in Section 2 we present a survey of an old and new approach
to SBT, and an explanation of the challenges involved when we try to extend
these approaches to account for rotation.

In Section 3 of this paper, we offer a line integral of regularized
singularities as a potential alternative to traditional SBT. Mathematically,
the approach of the original slender body pioneers [3, 19, 15] was to
distribute singularities along the centerline, then evaluate the flow from
these singularities on a fiber cross section. The flow is nonsingular since it
is _evaluated_ a fiber radius away from the location of the singularities. In
regularized singularities, we instead distribute the force and torque over a
“cloud” or “blob” region of characteristic size ${\hat{a}}\sim a$, which is
defined by some regularized delta function $\delta_{\hat{a}}$. In this case,
the velocity field is nonsingular because the force and torque are
_distributed_ over a region of finite size approximately equal to the fiber
radius. We can then convolve the velocity field (or half the vorticity field)
with the regularized delta function again to obtain the velocity (or angular
velocity) of the fiber centerline. In this way, we derive fundamental
regularized singularities for all components of the mobility operator, then
perform asymptotics on line integrals of the regularized singularities.

Representing filaments using regularized singularities is not a novel idea, as
it was first employed for numerical purposes in the immersed boundary (IB)
method of the early 1970s to model the flow patters inside the heart [36].
Later variations, including the force coupling method [27], method of
regularized Stokeslets [8, 9], and Rotne-Prager-Yamakawa (RPY) tensor [40,
46], focused on choosing the regularized delta function such that the Stokes
equations are solvable analytically. The accuracy of using regularized
singularities to approximate slender _translating_ filaments has been studied
numerically [6] and analytically [10, 29, 34, 51]. The conclusion of these
studies is that regularization methods can effectively match slender body
theory with a judicious choice of regularization radius ${\hat{a}}$. There has
yet to be an analysis, however, of how regularized singularities perform in
the context of _rotation_. For instance, do they reproduce the relationship
between torque and rotational velocity for a straight cylinder? And what about
the feedback of rotational velocity on translational motion for curved
filaments?

These are the main questions we try to answer in this paper. To do so, we
first look at how far we can get with traditional SBT in Section 2. We show
that the singularity methods of Keller and Rubinow and Johnson yield well-
defined mobilities for translation from force (trans-trans) and rotation from
torque (rot-rot) terms, but that both leave an unknown $\mathcal{O}(1)$
constant in the rotation-translation (rot-trans) coupling term which is
relevant when the fiber is driven by a large torque. Then, we perform
asymptotic expansions on line integrals of regularized RPY singularities to
derive a consistent SBT for filaments rotating and translating with speeds of
arbitrary order. We show that the regularized singularity radius ${\hat{a}}$
can be chosen to match either the trans-trans _or_ rot-rot mobility, but not
both. We also show that, unlike SBTs that are more faithful to the filament
geometry, the regularized RPY singularity approach gives symmetric rot-trans
coupling terms, and can easily be evaluated at the fiber endpoints without
having to modify the radius function.

The main drawback of our regularized singularity approach is that we lose
fidelity to the true tubular geometry of the filament. However, since we
expect the regularized singularity approach to capture the essential physics
of the interaction between the filament and fluid, we can utilize it as a
starting point to model the three-dimensional dynamics. To this end, in
Section 4 we use our results from Section 3 to postulate an asymptotic
relationship between the force and torque on a filament’s cross section and
its translational and rotational velocity. This equation has a single unknown
constant, which we estimate through three-dimensional boundary integral
simulations, thereby developing an empirical SBT for tubular rotating
filaments.

##  2 Three-dimensional SBTs for twist

In this section, we discuss how previous SBTs for translation might be
extended to also account for rotation. This section is not meant to be a
comprehensive list of previous SBTs (see [42] for a such a review), but rather
a focus on two different techniques and the challenges involved in applying
them to a (rapidly) _rotating_ filament. Our first technique is that used by
Keller and Rubinow [19] and Johnson [15, 16], which is to represent the flow
field by an integral of singularities positioned on the filament centerline.
The approach of Keller and Rubinow is technically distinct from that of
Johnson since they only use the singularity representation in the outer
expansion, but the end result for is the same, and so we discuss both in
Section 2.1. We then discuss the more recent approach of Koens and Lauga [22],
which evaluates the three-dimensional boundary integral representation of the
surface velocity asymptotically in $\epsilon$. The asymptotics of this
approach are more involved, since it maps traction on a surface to velocity on
a surface, rather than singularities on a line to velocity on the surface.

To establish some notation, we let $\boldsymbol{X}(s)$ be the fiber
centerline, where $s\in[0,L]$ is an arclength parameterization of the curve,
and thus $\boldsymbol{\tau}(s)=\partial_{s}\boldsymbol{X}(s)$ is the unit-
length tangent vector. The fiber surface $\widehat{\boldsymbol{X}}$ is then
defined as

$\widehat{\boldsymbol{X}}(s,\theta)=\boldsymbol{X}(s)+a\rho(s)\left(\boldsymbol{e}_{n}(s)\cos{\theta}+\boldsymbol{e}_{b}(s)\sin{\theta}\right):=\boldsymbol{X}(s)+a\rho(s)\widehat{\boldsymbol{r}}(s,\theta),$
(1)

where $a\rho(s)$, $0\leq\rho(s)\leq 1$, is the radius of the (circular) cross
section at $s$, $\boldsymbol{e}_{n}(s)=\boldsymbol{X}_{ss}(s)/\kappa(s)$ is
the fiber normal, $\kappa(s)$ is the centerline curvature at $s$, and
$\boldsymbol{e}_{b}=\boldsymbol{\tau}\times\boldsymbol{e}_{n}$ is the
binormal. The maximum fiber radius is $a$, which gives an aspect ratio
$\epsilon:=a/L$. Unless otherwise specified, we will assume that we are
working in a region of the fiber where $\rho(s)\equiv 1$; thus $a\rho=a$.

We will assume that the fiber centerline is translating with velocity
$\boldsymbol{U}(s)$ while its cross section rigidly rotates with angular
velocity $\Psi^{\parallel}(s)$ around the tangent vector
$\boldsymbol{\tau}(s)$, so that the velocity on the surface of the fiber is
given by

$\widehat{\boldsymbol{U}}(s,\theta)=\boldsymbol{U}(s)+a\Psi^{\parallel}(s)\left(\boldsymbol{\tau}(s)\times\widehat{\boldsymbol{r}}(s,\theta)\right)$
(2)

The goal of our analysis is to find an equation relating this motion to the
total force density $\boldsymbol{f}(s)$ and parallel torque density
$n^{\parallel}(s)$ on the fiber cross section. We assume here that the no-slip
condition applies to the fluid velocity $\boldsymbol{u}(\boldsymbol{x})$, so
that

$\boldsymbol{u}\left(\widehat{\boldsymbol{X}}(s,\theta)\right)=\widehat{\boldsymbol{U}}(s,\theta).$
(3)

###  2.1 Singularity representations

We first try to apply an older SBT approach to the spinning filament problem:
the singularity approach of Johnson [15], and later Götz [14]. The approach of
Johnson is to look for a series of singularities along the fiber centerline to
represent the flow field globally (everywhere outside the fiber). The
singularities are chosen so that the boundary condition (3) is satisfied with
an error $\mathcal{O}(\epsilon)$.

Using superposition, we can break the boundary velocity (2) into a velocity
with $\boldsymbol{U}=\boldsymbol{0}$ and another with $\Psi^{\parallel}=0$.
For $\Psi^{\parallel}=0$, Johnson shows that the boundary condition (3) can be
satisfied to $\mathcal{O}(\epsilon)$ by setting the global flow outside of the
fiber to be an integral of Stokeslets and doublets,

$\displaystyle\boldsymbol{u}^{(u)}(\boldsymbol{x})=\frac{1}{8\pi\mu}\int_{0}^{L}\left(\mathbb{S}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)+\frac{\left(a\rho(s)\right)^{2}}{2}\mathbb{D}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)\right)\boldsymbol{f}(s)\,ds$
(4)
$\displaystyle\text{where}\qquad\mathbb{S}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)=\left(\frac{\boldsymbol{I}}{r}+\frac{\boldsymbol{r}\boldsymbol{r}^{T}}{r^{3}}\right),\qquad\qquad\mathbb{D}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)=\left(\frac{\boldsymbol{I}}{r^{3}}-3\frac{\boldsymbol{r}\boldsymbol{r}^{T}}{r^{5}}\right),$

and $\boldsymbol{r}=\boldsymbol{x}-\boldsymbol{X}(s)$ with
$r=\left\|\boldsymbol{r}\right\|$. The doublet strength is necessary to cancel
the $\mathcal{O}(1)$ angular-dependent flow induced by the Stokeslet on a
filament cross section. Since the doublet has strength
$\mathcal{O}(\epsilon^{2})$, it makes a negligible contribution to the outer
expansion, but corrects the Stokeslet flow in the inner expansion. The
resulting asymptotic evaluation of $\boldsymbol{u}^{(u)}$ on the filament
cross section $\widehat{\boldsymbol{X}}(s,\theta)$ is independent of $\theta$
to $\mathcal{O}(\epsilon)$, and is given by

$\displaystyle
8\pi\mu\boldsymbol{u}^{(u)}\left(\widehat{\boldsymbol{X}}(s,\theta)\right)=\boldsymbol{U}_{\text{f}}^{\text{(LD)}}(s)+\boldsymbol{U}_{\text{f}}^{\text{(FP)}}+\mathcal{O}(\epsilon),$
(5)
$\displaystyle\boldsymbol{U}_{\text{f}}^{\text{(LD)}}(s)=\ln{\left(\frac{4s(L-s)}{a^{2}}\right)}\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)+\left(\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)$
(6)
$\displaystyle\boldsymbol{U}_{\text{f}}^{\text{(FP)}}(s)=\int_{0}^{L}\left(\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(t)\right)\boldsymbol{f}(t)-\left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|t-s|}\right)\boldsymbol{f}(s)\right)$
(7)

where the doublet contributes the
$\left(\boldsymbol{I}-3\boldsymbol{\tau}\boldsymbol{\tau}\right)$ term [15,
Eqs. (10–12)], and the term $\boldsymbol{U}_{\text{f}}^{\text{(FP)}}$ is a
finite part integral which is only well-defined as a difference of two terms
(the outer and inner expansions). The expression (5) has an error on the cross
section of order $\epsilon$, which means there will be a nonzero
$\mathcal{O}(1)$ rotational velocity of the cross section. To zero out the
rotational velocity, Johnson shows that it is necessary to introduce a rotlet,
a source, two stresslets, and two quadrupoles, with strengths given as a
function of the $\mathcal{O}(\epsilon)$ velocity generated by the surface flow
(5). While these singularities zero out the $\mathcal{O}(\epsilon)$ term, they
_do not_ impart any additional translational velocity on the cross section,
meaning that (5) gives the _translational_ velocity on the cross section to
order $\epsilon^{2}\log{\epsilon}$.

Let us now consider how the fluid flow outside the fiber relates to its
angular velocity. Johnson first sets

$n^{\parallel}=4\pi\mu a^{2}\Psi^{\parallel},$ (8)

which is the result for an infinite, straight cylinder [38, Eq. (62)]. The
singularity representation Johnson uses to match the angular velocity
everywhere is then given by a line integral of rotlets,

$\displaystyle
8\pi\mu\boldsymbol{u}^{(\Psi)}\left(\boldsymbol{x}\right)=\int_{0}^{L}\mathbb{R}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)n^{\parallel}(s)\,ds,\qquad\mathbb{R}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)=\frac{\boldsymbol{\tau}(s)\times\boldsymbol{r}}{r^{3}}$
(9)

Asymptotically, this velocity is equal to

$\displaystyle
8\pi\mu\boldsymbol{u}^{(\Psi)}\left(\widehat{\boldsymbol{X}}\left(s,\theta\right)\right)\approx\frac{2n^{\parallel}}{a}\left(\boldsymbol{\tau}\times\hat{\boldsymbol{r}}\right)+\frac{n^{\parallel}}{2}\left(\ln{\left(\frac{4s(L-s)}{a^{2}}\right)}-2\right)\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)+\boldsymbol{U}_{\text{n}}^{\text{(FP)}}$
(10)
$\displaystyle+n^{\parallel}\kappa\cos{\theta}\left(\boldsymbol{\tau}\times\hat{\boldsymbol{r}}\right)+\mathcal{O}\left(\epsilon\log{\epsilon}\right),$
$\displaystyle\text{where}\qquad\boldsymbol{U}_{\text{n}}^{\text{(FP)}}(s)=\int_{0}^{L}\left(\mathbb{R}\left(\boldsymbol{X}(s),\boldsymbol{X}(t)\right)n^{\parallel}(t)-\frac{1}{2}\left(\frac{\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)}{|t-s|}\right)n^{\parallel}(s)\right)$
(11)

on the fiber surface. Substituting (8) into (10), we see that the first term
is exactly the angular velocity of (2), while the second and third terms on
the first line of (10) give an additional constant velocity on the cross
section (the rot-trans coupling term). However, we still have an angular-
dependent term on the second line of (10): since
$\boldsymbol{\tau}\times\hat{\boldsymbol{r}}=\boldsymbol{e}_{b}\cos{\theta}-\boldsymbol{e}_{n}\sin{\theta}$,
the term in the second line has the angle dependence

$\displaystyle
n^{\parallel}\kappa\cos{\theta}\left(\boldsymbol{\tau}\times\hat{\boldsymbol{r}}\right)$
$\displaystyle=\frac{n^{\parallel}\kappa}{2}\left(\boldsymbol{e}_{b}+\boldsymbol{e}_{b}\cos{2\theta}-\boldsymbol{e}_{n}\sin{2\theta}\right).$
(12)

In Johnson’s analysis, the angular velocity $\Psi^{\parallel}$ and
translational velocity $\boldsymbol{U}$ are of the same order in $\epsilon$,
and so the contribution of all subleading terms in (10) to translational
velocity on the cross section is still $\epsilon^{2}$ relative to the
contribution of forcing. However, if force and torque are of the same order in
$\epsilon$ (as when a fiber is driven by a torque $n^{\parallel}$ independent
of $\epsilon$), the second line will make an $\mathcal{O}(1)$ contribution to
the cross-sectional velocity, violating the boundary condition (3). Thus,
there is a $\mathcal{O}(1)$ angle dependence coming from the rotlet term which
is unresolved. This dependence is the reason why the final integral equation
of Keller and Rubinow [19, Eq. (28)] for twisting filaments contains terms
which depend on the angle $\theta$ at which the inner expansion (of a
rotating, translating cylinder) is matched with the outer expansion (of
Stokeslets and rotlets), meaning there is no general solution for the
Stokeslet and rotlet strength which is independent of the matching angle.

In addition, there is also a lack of symmetry in the final result of Johnson,
since torque makes an $\mathcal{O}(\log{\epsilon})$ contribution to
$\boldsymbol{U}$ in (10), but force does not contribute to $\Psi^{\parallel}$
in (8). In our prototypical example of an unstable twirling fiber, this is not
an issue, since force makes a negligible contribution to rotational velocity
anyway.111Specifically, if the motor is driven by a constant torque
$n^{\parallel}$, the angular velocity $\Psi^{\parallel}\sim
n^{\parallel}/\epsilon^{2}$. Symmetry tells us that the contribution of force
to $\Psi^{\parallel}$ is $\mathcal{O}(\log{\epsilon})$, which means that the
rotational velocity from force is $\Psi_{\text{f}}^{\parallel}\sim
f\log{\epsilon}$, which is much smaller than the rotational velocity from
torque $\Psi_{\text{n}}^{\parallel}\sim n/\epsilon^{2}$, if $f$ and $n$ scale
similarly with $\epsilon$. In the case when $\Psi^{\parallel}$ does not scale
with $\epsilon$, however, the torque $n^{\parallel}$ is order $\epsilon^{2}$
and therefore makes a negligible contribution to $\boldsymbol{U}$, but the
force might make a nontrivial contribution to $\Psi^{\parallel}$. Thus there
are two main issues when applying singularity representations to twisting
filaments: first, the representation of the rotation in terms of rotlets is
insufficient to satisfy the boundary condition (3) to $\mathcal{O}(\epsilon)$,
and, second, there is a lack of symmetry in the rot-trans coupling terms.

###  2.2 Boundary integral representations

Recently, effort has been made to place the asymptotic theories of Keller and
Rubinow and Johnson on a more rigorous three-dimensional footing. Foremost
among these are Mori et al. [32, 33, 31], who show that slender body theory is
an $\mathcal{O}(\epsilon)$ approximation of a three-dimensional PDE with non-
standard boundary conditions, and Koens and Lauga [22], who show that the
translational SBT of Johnson and Keller and Rubinow can be obtained from
asymptotic expansion of the single layer222Koens and Lauga begin with the full
boundary integral representation, but the single layer representation is
sufficient as long as there is no change in the fiber volume. However, it is
important to point out that $\widehat{\boldsymbol{f}}$ in (13) is not in
general the force density on the fiber surface. Notable exceptions are a
filament undergoing rigid body motion or a “filament” whose interior is a
fluid of viscosity equal to that of the fiber exterior. boundary integral
equation [39, Sec. 4.1]

$8\pi\mu\,\widehat{\boldsymbol{U}}\left(s,\theta\right)=\int_{0}^{2\pi}\int_{0}^{L}\mathbb{S}\left(\widehat{\boldsymbol{X}}(s,\theta),\widehat{\boldsymbol{X}}(s^{\prime},\theta^{\prime})\right)\widehat{\boldsymbol{f}}\left(s^{\prime},\theta^{\prime}\right)\,ds^{\prime}d\theta^{\prime}.$
(13)

Here we assume the fiber is a thin shell filled with fluid, so that
$\widehat{\boldsymbol{f}}(s,\theta)$, which is the jump in surface traction
times the surface area element [39, p. 105], is also equal to the force
density on the fiber surface. Koens and Lauga show that the constant cross-
sectional $\boldsymbol{U}(s)$ in (2) is related to the total force on a cross
section,

$\boldsymbol{f}(s)=\int_{0}^{2\pi}\widehat{\boldsymbol{f}}\left(s,\theta\right)\,d\theta$
(14)

by the classical SBT (5), and that, for a _straight_ cylinder, the torque-
angular velocity relationship (8) is recovered [22, Eq. (4.11)] for the total
parallel torque on a cross section,

$n^{\parallel}(s)=\boldsymbol{\tau}(s)\cdot\int_{0}^{2\pi}a\rho(s)\left(\widehat{\boldsymbol{r}}(s,\theta)\times\widehat{\boldsymbol{f}}(s,\theta)\right)\,d\theta.$
(15)

Here we briefly discuss how the approach of Koens and Lauga might be extended
to account for rotation in curved filaments. The main issue is that, if the
boundary integral equation (13) is only expanded to leading order (i.e., with
an error of $\mathcal{O}(\epsilon)$), it is impossible to make observations
about the rotational velocity, since the velocity on the cross section induced
by an angular velocity $\Psi^{\parallel}$ is $\sim\epsilon\Psi^{\parallel}$.
Thus more terms are necessary in the asymptotics. Unfortunately, because the
expansions must be conducted around both $\theta$ and $s$ in (13), even
expanding a single layer formulation to the next order is complex.

To get a flavor for the challenges involved, we consider an asymptotic
expansion of the _isotropic part_ of the single-layer potential (c.f. [22, Eq.
(5.1)]

$\boldsymbol{U}_{i1}(s,\theta):=\int_{0}^{2\pi}\int_{0}^{L}\frac{\boldsymbol{f}(s^{\prime},\theta^{\prime})}{R\left(\widehat{\boldsymbol{X}}(s,\theta),\widehat{\boldsymbol{X}}(s^{\prime},\theta^{\prime})\right)}\,d\theta^{\prime}\,ds^{\prime},$
(16)

in an inner region where $s-s^{\prime}=\mathcal{O}(a)$. The details of the
expansion up to integration over $\theta$ are given in Appendix A. In order to
integrate over $\theta^{\prime}$, we need to assume a Fourier-series
representation of $\boldsymbol{f}$. For simplicity, we will look only at the
first two terms,

$\widehat{\boldsymbol{f}}(s,\theta)=\frac{1}{2\pi}\boldsymbol{f}(s)+\frac{1}{2\pi}\left(\boldsymbol{f}_{c,1}(s)\cos{\theta}+\boldsymbol{f}_{s,1}(s)\sin{\theta}\right),$
(17)

which we substitute into (A.5). After integrating over $\theta^{\prime}$, the
final inner expansion for the surface velocity is333We use the variable
$\theta$ here for notational convenience; more precisely, $\theta$ should be
measured relative to the curve torsion angle $\theta_{i}$(s), which is defined
in Appendix A. That is, $\theta$ in (17) and (18) stands for
$\theta-\theta_{i}(s)$.

$\displaystyle\boldsymbol{U}_{i1}(s,\theta)=\ln{\left(\frac{4s(L-s)}{a^{2}}\right)}\boldsymbol{f}(s)+(L-2s)\partial_{s}\boldsymbol{f}(s)+\boldsymbol{f}_{c,1}(s)\cos{\theta}$
(18)
$\displaystyle-\frac{a\kappa}{2}\left(\boldsymbol{f}(s)\cos{\theta}+\left(\frac{1}{4}\cos{2\theta}+\frac{1}{2}\left(\ln{\left(\frac{4s(L-s)}{a^{2}}\right)}-2\right)\right)\boldsymbol{f}_{c,1}(s)-\frac{1}{2}\boldsymbol{f}_{s,1}(s)\cos{\theta}\sin{\theta}\right)+\mathcal{O}(\epsilon^{2}).$

We can make three important observations from (18): first, the term
$\boldsymbol{f}(s)\kappa\cos{\theta}$ demonstrates rot-trans coupling
(rotational motion from constant forcing) which results from including
additional terms in the SBT expansion, or more specifically from accounting
for the curvature of the fiber in the inner expansion. Likewise, the _torque_
(15) generated from the traction scales as $a\boldsymbol{f}_{c/s,1}$, and thus
(18) tells us that torque makes an $\mathcal{O}(\log{\epsilon})$ contribution
to $\boldsymbol{U}(s)$, as we expect from the result (10) of singularity
methods.

While the result of Koens and Lauga is vital to understanding translational
SBT, there are clearly limitations to extending their approach to rotation.
Even for the simplest possible term in the boundary integral formulation
(isotropic part of the single layer), there are many terms that have to be
tracked if we expand to $\mathcal{O}(\epsilon)$, i.e., to an error of
$\mathcal{O}\left(\epsilon^{2}\right)$. More importantly, it is not even clear
if we can solve (18), since the Fourier modes no longer decouple as they do
for translation.

Summing up this section, we have seen that a slender body theory which
properly accounts for twist dynamics _and_ is faithful to the three-
dimensional fiber geometry remains elusive. In order to make some progress, we
make an approximation of the fiber geometry and define the fiber as a series
of infinitely many regularized singularities. Then, we use our regularized
singularity SBT to inform an empirical SBT for the full three-dimensional
geometry.

##  3 Slender body theory from regularized singularities

Our survey in Section 2 reveals the challenging nature of developing a slender
body theory for twist. Because of the separation of scales between the angular
and translational velocity, more terms are required in the asymptotic
analysis, which makes these SBTs tedious to derive and cumbersome to work with
in practice.

An alternative approach, which we consider in this section, is to relax the
fidelity to a true cylindrical geometry and instead represent the cylinder
using _regularized_ singularities distributed along the centerline. The idea
of the regularized singularity approach is to replace the delta function in
the force and torque with a regularized delta function $\delta_{\hat{a}}$,
which is smoothed over some “regularization radius” ${\hat{a}}$. The pointwise
fluid velocity due to the regularized singularities is then obtained by
solving the Stokes equations [30, Sec. 2.1]

$\nabla\pi(\boldsymbol{x})=\mu\nabla^{2}\boldsymbol{u}(\boldsymbol{x})+\int_{0}^{L}\left(\left(\boldsymbol{f}(s)+\frac{\nabla}{2}\times
n^{\parallel}(s)\boldsymbol{\tau}(s)\right)\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{X}(s)\right)\right)\,ds,$
(19)

and $\nabla\cdot\boldsymbol{u}=0$. The translational and parallel rotational
velocity of the fiber centerline are then given by convolving the pointwise
velocity $\boldsymbol{u}(\boldsymbol{x})$ and half vorticity once more with
the regularized delta function,

$\displaystyle\boldsymbol{U}\left(s\right)=\int\boldsymbol{u}(\boldsymbol{x})\,\delta_{{\hat{a}}}\left(\boldsymbol{X}(s)-\boldsymbol{x}\right)\,d\boldsymbol{x}$
(20)
$\displaystyle\Psi^{\parallel}\left(s\right)=\boldsymbol{\tau}(s)\cdot\int\delta_{{\hat{a}}}\left(\boldsymbol{X}(s)-\boldsymbol{x}\right)\frac{\nabla}{2}\times\boldsymbol{u}\left(\boldsymbol{x}\right)\,d\boldsymbol{x}.$
(21)

This yields a mobility that is guaranteed to be symmetric positive definite
[17, Footnote 3]. Here we only consider the parallel component of rotational
velocity, since, for unshearable filaments, the other components can be
derived from the translational velocity $\boldsymbol{U}(s)$ [30, Sec. 2].

Using superposition, the translational and rotational velocity of the fiber
centerline can be defined as integrals of the fundamental _regularized_
singularities

$\displaystyle\boldsymbol{U}(s)$
$\displaystyle=\int_{0}^{L}\left(\mathbb{M}_{\text{tt}}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}(s^{\prime})+\mathbb{M}_{\text{tr}}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)n^{\parallel}(s^{\prime})\boldsymbol{\tau}(s^{\prime})\right)\,ds^{\prime}$
(22) $\displaystyle\Psi^{\parallel}(s)$
$\displaystyle=\boldsymbol{\tau}(s)\cdot\int_{0}^{L}\left(\mathbb{M}_{\text{rt}}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}(s^{\prime})+\mathbb{M}_{\text{rr}}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)n^{\parallel}(s^{\prime})\boldsymbol{\tau}(s^{\prime})\right)\,ds^{\prime}$

where the fundamental regularized singularities
$\mathbb{M}_{\text{tt}},\mathbb{M}_{\text{tr}},\mathbb{M}_{\text{rt}},$ and
$\mathbb{M}_{\text{rr}}$ are obtained by convolving the Stokeslet with
$\delta_{\hat{a}}$, with a half curl taken each time rotation comes in,

$\displaystyle\mathbb{M}_{\text{tt}}\left(\boldsymbol{x},\boldsymbol{y}\right)$
$\displaystyle=\int\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{z}\right)\int\mathbb{S}\left(\boldsymbol{z},\boldsymbol{w}\right)\delta_{{\hat{a}}}\left(\boldsymbol{w}-\boldsymbol{y}\right)\,d\boldsymbol{w}\,d\boldsymbol{z},$
(23)
$\displaystyle\mathbb{M}_{\text{tr}}\left(\boldsymbol{x},\boldsymbol{y}\right)$
$\displaystyle=\int\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{z}\right)\int\mathbb{S}\left(\boldsymbol{z},\boldsymbol{w}\right)\frac{\nabla_{\boldsymbol{w}}}{2}\times\delta_{{\hat{a}}}\left(\boldsymbol{w}-\boldsymbol{y}\right)\,d\boldsymbol{w}\,d\boldsymbol{z},$
$\displaystyle\mathbb{M}_{\text{rt}}\left(\boldsymbol{x},\boldsymbol{y}\right)$
$\displaystyle=\int\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{z}\right)\frac{\nabla_{\boldsymbol{z}}}{2}\times\int\mathbb{S}\left(\boldsymbol{z},\boldsymbol{w}\right)\delta_{{\hat{a}}}\left(\boldsymbol{w}-\boldsymbol{y}\right)\,d\boldsymbol{w}\,d\boldsymbol{z}=\mathbb{M}_{\text{tr}}^{T}\left(\boldsymbol{y},\boldsymbol{x}\right),$
$\displaystyle\mathbb{M}_{\text{rr}}\left(\boldsymbol{x},\boldsymbol{y}\right)$
$\displaystyle=\int\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{z}\right)\frac{\nabla_{\boldsymbol{z}}}{2}\times\int\mathbb{S}\left(\boldsymbol{z},\boldsymbol{w}\right)\frac{\nabla_{\boldsymbol{w}}}{2}\times\delta_{{\hat{a}}}\left(\boldsymbol{w}-\boldsymbol{y}\right)\,d\boldsymbol{w}\,d\boldsymbol{z}.$

Using this approach, the four components of the mobility operator relating
force and torque to angular and rotational velocity can be written
_separately_ , giving four different integrals for the four different pieces
of the velocity in (22). Thus, we can simply do asymptotics on each of these
four integrals to get an approximation to each term to some order in
$\epsilon$. This is in contrast to the SBTs of Section 2, where the
translational velocity on the cross section was some combination of
contributions from force and torque, and we were left with trying to determine
the order of each term from the desired order of $\boldsymbol{U}$ and
$\Psi^{\parallel}$ on the cross section.

Once we define the velocity in (22), it remains to specify the regularized
singularity kernels and perform the integration. Our choice of kernel is the
Rotne-Prager-Yamakawa (RPY) tensor [40, 46], for which
$\delta_{\hat{a}}(\boldsymbol{r})=\delta(r-{\hat{a}})/(4\pi{\hat{a}}^{2})$,
i.e., the regularized delta function is a surface delta function on a sphere
of radius ${\hat{a}}$. This choice is more faithful to the original
cylindrical geometry than delta functions with infinite support [8, 9, 28] or
those whose support is tied to a numerical fluid grid [37]. In particular,
since each regularized singularity is a sphere, it is not hard to imagine this
resulting in a cylinder of constant radius ${\hat{a}}$ along a length $L$,
with hemispherical caps at each end. Still, it is an approximation, since the
regularized delta function is the same in the axial and radial directions, and
also because the RPY tensor does not include stresslet terms required to keep
the spheres rigid in flow [12], even though stresslets are of the same
multipole order as torques, which are included.

In a companion paper [30], we develop efficient quadrature schemes to evaluate
the integrals (22) to numerical, rather than asymptotic, precision. The RPY
kernels change their form when the distance between the spheres is less than
$2{\hat{a}}$, i.e., when the spheres overlap, thus regularizing the
singularity that occurs in the interactions between point forces. Still, the
kernels are nearly singular, and in the limit ${\hat{a}}/L:=\hat{\epsilon}\ll
1$ it is more efficient to evaluate them asymptotically, since some of the
singularities can safely be assumed to only contribute to the “inner”
expansion.

In this section, we therefore perform asymptotic expansions on the integrals
(22) to obtain an SBT for a twisting filament. Our strategy is standard
matched asymptotics and similar to the approaches of Johnson [15] and Gotz
[14] for SBT. Because we use regularized singularities, however, all of the
integrals take place on the fiber centerline, and so the asymptotics are in
one dimension, rather than on the fiber cross section (two dimensions). It
therefore becomes much simpler to obtain expressions for the velocity at the
endpoints (see Appendix B), where we only have to redefine the domain of
integration in the inner expansion.

For each integral in (22), we first compute an outer expansion by considering
the region where $|s-s^{\prime}|$ is $\mathcal{O}(1)$. Then in the inner
expansion, we consider the part of the integrals where $|s-s^{\prime}|$ is
$\mathcal{O}({\hat{a}})$. In this case, we follow [14] and introduce the
rescaled variable

$\xi=\frac{s^{\prime}-s}{{\hat{a}}},$ (24)

so that $\xi$ is $\mathcal{O}(1)$ in the inner expansion, and the domain of
$\xi$ is $[-s/{\hat{a}},(L-s)/{\hat{a}}]$. With this definition, it is
straightforward to integrate over $\xi$ and obtain an inner expansion, which
is done in two pieces because the RPY tensor changes at $s=2{\hat{a}}$ (the
case of overlapping spheres). Adding the inner and outer solutions together
and subtracting the common part then gives the final matched asymptotic
solution. Unless otherwise specified, all expansions will be carried out to
$\mathcal{O}(\hat{\epsilon})$, i.e., in this section $\approx$ denotes
equality up to $\mathcal{O}(\hat{\epsilon})$.

###  3.1 Translation from force

We begin by expanding the translational velocity from force, which we define
as

$\boldsymbol{U}_{\text{f}}=\int_{0}^{L}\mathbb{M}_{\text{tt}}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}(s^{\prime})\,ds^{\prime}.$
(25)

The following derivation is a repeat of that given in [29, Appendix A], but it
is included here for completeness, and because the set of steps for the other
3 blocks of the mobility operator is identical. When we substitute the
definition of $\mathbb{M}_{\text{tt}}$ for RPY singularities from [46, Eq.
(3.12)] into (25), we have the integral

$\displaystyle 8\pi\mu\boldsymbol{U}_{\text{f}}\left(s\right)=$
$\displaystyle\int_{R>2{\hat{a}}}\left(\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)+\frac{2{\hat{a}}^{2}}{3}\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\right)\boldsymbol{f}(s^{\prime})\,ds^{\prime}$
(26) $\displaystyle+$ $\displaystyle\int_{R\leq
2{\hat{a}}}\left(\left(\frac{4}{3{\hat{a}}}-\frac{3R\left(s^{\prime}\right)}{8{\hat{a}}^{2}}\right)\boldsymbol{I}+\frac{1}{8{\hat{a}}^{2}R\left(s^{\prime}\right)}\left(\boldsymbol{R}\boldsymbol{R}\right)\left(s^{\prime}\right)\right)\boldsymbol{f}\left(s^{\prime}\right)\,ds^{\prime},$

where
$\boldsymbol{R}\left(s^{\prime}\right)=\boldsymbol{X}(s)-\boldsymbol{X}\left(s^{\prime}\right)$,
and $R=\left\|\boldsymbol{R}\right\|$. The separation of the integrals
captures the change in the RPY tensor when $R\leq 2{\hat{a}}$. Because $s$ is
an arclength parameterization, we will replace the bound $R>2{\hat{a}}$ with
$|s-s^{\prime}|>2{\hat{a}}$ going forward. This incurs a relative error on the
order $\hat{\epsilon}$, which is the same order to which we carry the
asymptotic expansions.

In the outer expansion, we consider the part of the integral (26) where
$|s-s^{\prime}|$ is $\mathcal{O}(1)$. In this case, the doublet term in (26)
is insignificant and we obtain the outer velocity by integrating the Stokeslet
over the fiber centerline,

$8\pi\mu\boldsymbol{U}_{\text{f}}^{(\text{outer})}(s)=\int_{|s-s^{\prime}|>2{\hat{a}}}\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}\left(s^{\prime}\right)\,ds^{\prime}.$
(27)

The part of the kernel (26) for $R\leq 2{\hat{a}}$ makes no contribution to
the outer expansion since $|s-s^{\prime}|$ is $\mathcal{O}({\hat{a}})$ there.

The inner expansion of (26) is derived in a straightforward way by
substituting the Taylor expansions (B.1) into the RPY velocity (26) and
integrating over $\xi$. The details are given in Appendix B, with the final
result that

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{f}}^{(\text{inner})}(s)=\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)\ln{\left(\frac{(L-s)s}{4{\hat{a}}^{2}}\right)}$
(28)
$\displaystyle+\left(\left(4-\frac{{\hat{a}}^{2}}{3s^{2}}-\frac{{\hat{a}}^{2}}{3(L-s)^{2}}\right)\boldsymbol{I}+\left(\frac{{\hat{a}}^{2}}{s^{2}}+\frac{{\hat{a}}^{2}}{(L-s)^{2}}\right)\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s),$

in the fiber interior. The complete formulas including endpoint modifications
are given in (B.7); the $\mathcal{O}\left({\hat{a}}^{2}\right)$ terms are
included in (28) because they are necessary to give $\mathcal{C}^{2}$
continuity with the endpoint formulas.

The common part is the outer velocity written in terms of the inner variables,

$8\pi\mu\boldsymbol{U}_{\text{f}}^{(\text{common})}(s)=\int_{|s-s^{\prime}|>2{\hat{a}}}\left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s-s^{\prime}|}\right)\boldsymbol{f}(s)\,ds^{\prime}.$
(29)

The total velocity is then the sum of the inner and outer expansions, with the
common part subtracted,

$\displaystyle\boldsymbol{U}_{\text{f}}(s)=\boldsymbol{U}_{\text{f}}^{(\text{inner})}(s)+\boldsymbol{U}_{\text{f}}^{(\text{outer})}(s)-\boldsymbol{U}_{\text{f}}^{(\text{common})}(s).$
(30)

This can be written as

$\displaystyle\boldsymbol{U}_{\text{f}}(s)=\boldsymbol{U}_{\text{f}}^{(\text{inner})}(s)+\frac{1}{8\pi\mu}\int_{|s-s^{\prime}|>2{\hat{a}}}\left(\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}\left(s^{\prime}\right)-\left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s-s^{\prime}|}\right)\boldsymbol{f}(s)\right)\,ds^{\prime},$
(31)

where $\boldsymbol{U}_{\text{f}}^{(\text{inner})}$ is defined in (B.7). The
velocity (31) has a form similar to that of SBT, but with different bounds on
the integral. Unlike SBT, the integral (31) is not a finite part integral, but
a nearly singular integral that makes sense for each term separately.
Numerical evidence suggests that, unlike traditional SBT, our asymptotic
formula (31) keeps the eigenvalues of the mobility above zero [30, Sec 4.4.1].

To compare with SBT, however, it is advantageous to observe that the integrand
in (31) is $\mathcal{O}(\hat{\epsilon})$ when $R\leq 2{\hat{a}}$, and so we
can add the excluded part back in without changing the asymptotic accuracy.
This gives a velocity of the exact same form as SBT,

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{f}}(s)=\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)\ln{\left(\frac{(L-s)s}{4{\hat{a}}^{2}}\right)}+4\boldsymbol{f}(s)$
(32)
$\displaystyle+\int_{0}^{L}\left(\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}\left(s^{\prime}\right)-\left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s-s^{\prime}|}\right)\boldsymbol{f}(s)\right)\,ds^{\prime},$

in the fiber interior, where we have dropped $\mathcal{O}({\hat{a}}^{2})$
terms from (28). The comparison of this formula to translational SBT will be
given in Section 3.5.

###  3.2 Translation from torque

To complete our formulation for translational velocity, we next consider an
asymptotic expansion of the velocity $\boldsymbol{U}_{\text{n}}(s)$ due to a
torque density $n^{\parallel}(s)$ on the fiber centerline,

$\boldsymbol{U}_{\text{n}}(s)=\int_{0}^{L}\mathbb{M}_{\text{tr}}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)n^{\parallel}(s^{\prime})\boldsymbol{\tau}(s^{\prime})\,ds^{\prime}.$
(33)

When we substitute the definition of $\mathbb{M}_{\text{tr}}$ for RPY
singularities from [46, Eq. (3.16)] into (33), we have the integral

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{n}}(s)=\int_{R>2{\hat{a}}}\frac{\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})}{R(s^{\prime})^{3}}n^{\parallel}(s^{\prime})\,ds^{\prime}$
(34) $\displaystyle+\frac{1}{2{\hat{a}}^{2}}\int_{R\leq
2{\hat{a}}}\left(\frac{1}{{\hat{a}}}-\frac{3R(s^{\prime})}{8{\hat{a}}^{2}}\right)\left(\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)n^{\parallel}(s^{\prime})ds^{\prime},$

where
$\boldsymbol{R}(s^{\prime})=\boldsymbol{X}(s)-\boldsymbol{X}(s^{\prime})$. The
first term in (34) is the Rotlet, while the second term is the RPY tensor when
$R\leq 2{\hat{a}}$. As before, we will use $R\approx|s-s^{\prime}|$ to modify
the bounds on the integrals.

In the outer expansion, we consider the part of the integral (34) where
$|s^{\prime}-s|$ is $\mathcal{O}(1)$. In this case the correction term is
insignificant and we get the outer velocity

$8\pi\mu\boldsymbol{U}_{\text{n}}^{\text{(outer)}}(s)=\int_{|s-s^{\prime}|>2{\hat{a}}}\frac{\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})}{R(s^{\prime})^{3}}n^{\parallel}(s^{\prime})\,ds^{\prime}.$
(35)

In the inner expansion, we use the inner asymptotics (B.1) to write the inner
expansion as

$\displaystyle 8\pi\mu\boldsymbol{U}_{\text{n}}^{\text{(inner)}}(s)$
$\displaystyle=\frac{1}{2}\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)\left(\ln{\left(\frac{s(L-s)}{4{\hat{a}}^{2}}\right)}+\frac{7}{6}\right)$
(36)

in the fiber interior, with the corresponding endpoint formula given in
(B.11). The details of this inner expansion can be found in Appendix B.

The common part for trans-rot coupling is the inner expansion written in terms
of the outer variables,

$8\pi\mu\boldsymbol{U}_{\text{n}}^{\text{(common)}}(s)=\int_{|s-s^{\prime}|>2{\hat{a}}}\left(\frac{\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)}{2|s-s^{\prime}|}\right)n^{\parallel}(s)\,ds^{\prime}.$
(37)

To $\mathcal{O}(\hat{\epsilon})$, the velocity due to torque
$\boldsymbol{n}(s)=n(s)\boldsymbol{\tau}(s)$ on the fiber centerline is given
by

$\displaystyle\boldsymbol{U}_{\text{n}}(s)=\boldsymbol{U}_{\text{n}}^{\text{(inner)}}(s)+\boldsymbol{U}_{\text{n}}^{\text{(outer)}}(s)-\boldsymbol{U}_{\text{n}}^{\text{(common)}}(s)$
(38)
$\displaystyle\boldsymbol{U}_{\text{n}}(s)=\boldsymbol{U}_{\text{n}}^{\text{(inner)}}(s)+\frac{1}{8\pi\mu}\int_{R>2{\hat{a}}}\left(\frac{\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})}{R(s^{\prime})^{3}}n^{\parallel}(s^{\prime})-\frac{\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)}{2|s-s^{\prime}|}n^{\parallel}(s)\right)\,ds^{\prime},$
(39)

We can make this velocity (39) more SBT-like by changing the bounds on the
integral. Because the outer and common part match to
$\mathcal{O}(\hat{\epsilon})$ when $R\leq 2{\hat{a}}$, we can change the limit
of integration in the integral to obtain a finite part integral, which gives
the final result

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{n}}(s)=\frac{1}{2}\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)\left(\ln{\left(\frac{s(L-s)}{4{\hat{a}}^{2}}\right)}+\frac{7}{6}\right)$
(40)
$\displaystyle+\int_{0}^{L}\left(\left(\frac{\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})}{R(s^{\prime})^{3}}\right)n^{\parallel}(s^{\prime})-\left(\frac{\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)}{2|s-s^{\prime}|}\right)n^{\parallel}(s)\right)\,ds^{\prime},$

for the velocity in the fiber interior. This is of exactly the same form as
the first line of the singularity representation (10), but without the
additional angle-dependent terms in the second line of (10).

###  3.3 Rotation from force

We next calculate the parallel rotational velocity
$\Psi_{\text{f}}^{\parallel}(s)$ due to the force density $\boldsymbol{f}(s)$
on the fiber centerline,

$\Psi_{\text{f}}^{\parallel}(s)=\boldsymbol{\tau}(s)\cdot\int_{0}^{L}\mathbb{M}_{\text{rt}}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}(s^{\prime})\,ds^{\prime}.$
(41)

When we substitute the definition of $\mathbb{M}_{\text{tr}}$ from [46, Eq.
(3.16)] into (41), we have

$\displaystyle
8\pi\mu\Psi_{\text{f}}^{\parallel}(s)=\int_{R>2{\hat{a}}}\frac{\left(\boldsymbol{f}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)}{R(s^{\prime})^{3}}\cdot\boldsymbol{\tau}(s)\,ds^{\prime}$
(42)
$\displaystyle+\frac{1}{2{\hat{a}}^{2}}\int_{R<2{\hat{a}}}\left(\frac{1}{{\hat{a}}}-\frac{3R(s^{\prime})}{8{\hat{a}}^{2}}\right)\left(\boldsymbol{f}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)\cdot\boldsymbol{\tau}(s)\,ds^{\prime}.$

After using the triple cross product identity to rewrite
$\left(\boldsymbol{f}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)\cdot\boldsymbol{\tau}(s)=\left(\boldsymbol{R}(s^{\prime})\times\boldsymbol{\tau}(s)\right)\cdot\boldsymbol{f}(s^{\prime})$,
the symmetry of the kernels $\mathbb{M}_{\text{rt}}$ and
$\mathbb{M}_{\text{tr}}$ makes the asymptotics in the inner, outer, and common
expansions much the same as the previous section. The final result is that

$\displaystyle
8\pi\mu\Psi_{\text{f}}^{\parallel}(s)=\frac{1}{2}\left(\ln{\left(\dfrac{s(L-s)}{4{\hat{a}}^{2}}\right)}+\dfrac{7}{6}\right)\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)\cdot\boldsymbol{f}(s)+\Psi_{\text{f}}^{\text{(FP)}}(s)$
(43)
$\displaystyle\Psi_{\text{f}}^{\text{(FP)}}(s)=\int_{0}^{L}\left(\frac{\boldsymbol{R}(s^{\prime})\times\boldsymbol{\tau}(s)}{R(s^{\prime})^{3}}\cdot\boldsymbol{f}(s^{\prime})-\frac{\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)}{2|s-s^{\prime}|}\cdot\boldsymbol{f}(s)\right)\,ds^{\prime}.$
(44)

in the fiber interior. To obtain a formula accurate to
$\mathcal{O}(\hat{\epsilon})$ up to and including the fiber endpoints, we
replace the first line in (43) with $\Psi_{\text{f}}^{\parallel,\text{inner}}$
defined in (B.12).

Using the endpoint formulas in Appendix B, it is straightforward to show the
symmetry (self-adjointness in $L^{2}$) of the mobility operator,

$\int_{0}^{L}\Psi_{\text{f}}^{\parallel}(s)n^{\parallel}(s)\,ds=\int_{0}^{L}\boldsymbol{U}_{\text{n}}(s)\cdot\boldsymbol{f}(s)\,ds.$
(45)

The only nontrivial part of this calculation is the rotlet term in the finite
part integrals, in which the order of integration in $s$ and $s^{\prime}$ must
be swapped.

###  3.4 Rotation from torque

The final asymptotic expansion we need is the rotational velocity
$\boldsymbol{\Psi}(s)$ due to the torque $n^{\parallel}(s)$. Because the
doublet singularity has already been expanded in the translational case in
(B.4), it will be convenient to work with a full vector rotational velocity
and torque

$\boldsymbol{\Psi}_{\text{n}}(s)=\int_{0}^{L}\mathbb{M}_{\text{rr}}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{n}(s^{\prime})\,ds^{\prime},$
(46)

and specialize to the case of parallel velocity and torque later. We
substitute the definition of $\mathbb{M}_{\text{rr}}$ from [46, Eq. (3.14)]
into (46) to obtain the vector rotational velocity,

$\displaystyle
8\pi\mu\boldsymbol{\Psi}_{\text{n}}(s)=-\frac{1}{2}\int_{R>2{\hat{a}}}\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{n}(s^{\prime})\,ds^{\prime}$
(47) $\displaystyle+\frac{1}{{\hat{a}}^{3}}\int_{R\leq
2{\hat{a}}}\left(\left(1-\frac{27R(s^{\prime})}{32{\hat{a}}}+\frac{5R(s^{\prime})^{3}}{64{\hat{a}}^{3}}\right)\boldsymbol{I}+\left(\frac{9}{32{\hat{a}}R(s^{\prime})}-\frac{3R(s^{\prime})}{64{\hat{a}}^{3}}\right)\left(\boldsymbol{R}\boldsymbol{R}\right)(s^{\prime})\right)\boldsymbol{n}(s^{\prime})\,ds^{\prime}.$

The first term in (47) is the doublet, while the second term is the RPY tensor
for $R\leq 2{\hat{a}}$.

In the outer expansion, we consider the part of the integral (47) where
$|s^{\prime}-s|$ is $\mathcal{O}(1)$,

$8\pi\mu\boldsymbol{\Psi}_{\text{n}}^{\text{(outer)}}(s)=-\frac{1}{2}\int_{|s-s^{\prime}|>2{\hat{a}}}\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{n}(s^{\prime})\,ds^{\prime}.$
(48)

Importantly, the outer expansion (and the common expansion) is
$\mathcal{O}(1)$. The inner expansion will be
$\mathcal{O}\left(\hat{\epsilon}^{-2}\right)$, so we will wind up dropping the
entire outer expansion.

In Appendix B, we use the inner asymptotics (B.1) to derive the inner angular
velocity

$\displaystyle
8\pi\mu\boldsymbol{\Psi}_{\text{n}}^{\text{(inner)}}(s)=\left(p_{I}(s)\boldsymbol{I}+p_{\tau}(s)\boldsymbol{\boldsymbol{\tau}}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{n}(s),\qquad\text{where}$
(49) $\displaystyle
p_{I}(s)=\dfrac{1}{{\hat{a}}^{2}}\left(\dfrac{9}{8}+\dfrac{{\hat{a}}^{2}}{4}\left(\dfrac{1}{s^{2}}+\dfrac{1}{(L-s)^{2}}\right)\right)$
$\displaystyle
p_{\tau}(s)=\dfrac{1}{{\hat{a}}^{2}}\left(\dfrac{9}{8}-\dfrac{3{\hat{a}}^{2}}{4}\left(\dfrac{1}{s^{2}}+\dfrac{1}{(L-s)^{2}}\right)\right),$

in the fiber interior. The functions $p_{I}$ and $p_{\tau}$ are different at
the fiber endpoints, and are given in (B.15). As in the translation case, the
$\mathcal{O}\left({\hat{a}}^{2}\right)$ terms in (49) are negligible in the
fiber interior, but are necessary to give $\mathcal{C}^{2}$ continuity with
the endpoint formulas in (B.15).

Because the outer expansion is $\mathcal{O}(\hat{\epsilon}^{2})$ smaller than
the inner expansion, the total asymptotic velocity is just given by the inner
expansion,

$\boldsymbol{\Psi}_{\text{n}}=\boldsymbol{\Psi}_{\text{n}}^{\text{(inner)}}.$
(50)

Specializing now to the case of parallel angular velocity from parallel
torque, $\boldsymbol{n}=n^{\parallel}\boldsymbol{\tau}$, and
$\Psi_{\text{n}}^{\parallel}=\boldsymbol{\Psi}_{\text{n}}\cdot\boldsymbol{\tau}$,
and dropping the negligible terms in the fiber interior from (49), we obtain
the simple formula

$8\pi\mu\Psi_{\text{n}}^{\parallel}(s)=\frac{9n^{\parallel}(s)}{4{\hat{a}}^{2}}$
(51)

in the fiber interior.

###  3.5 Comparison to Keller-Rubinow-Johnson

Let us compare our RPY-based slender body theory to the results we have from
Section 2. We first compare the theories when we have a well-defined SBT,
i.e., for translation-translation and rotation-rotation. Then, we use the form
of the translation-rotation RPY coupling and our results in Section 2 to
postulate an SBT equation for trans-rot coupling with a single unknown
parameter. The parameter will be determined in the next section.

####  3.5.1 Trans-trans and rot-rot

All three-dimensional slender body theories agree that the translational
velocity in the case when $\Psi^{\parallel}=0$ is given by

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{f}}=\ln{\left(\frac{4s(L-s)}{a^{2}}\right)}\left(\boldsymbol{I}+\boldsymbol{\tau}\boldsymbol{\tau}\right)\boldsymbol{f}+\left(\boldsymbol{I}-3\boldsymbol{\tau}\boldsymbol{\tau}\right)\boldsymbol{f}+\boldsymbol{U}_{\text{f}}^{\text{(FP)}}+\mathcal{O}(\epsilon),$
(52)

while our RPY-based theory (32) gives

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{f}}=\ln{\left(\frac{s(L-s)}{4{\hat{a}}^{2}}\right)}\left(\boldsymbol{I}+\boldsymbol{\tau}\boldsymbol{\tau}\right)\boldsymbol{f}+4\boldsymbol{f}+\boldsymbol{U}_{\text{f}}^{\text{(FP)}}+\mathcal{O}(\epsilon).$
(53)

For the case of rotation from parallel torque, all theories of Section 2 agree
that

$8\pi\mu\Psi_{\text{n}}^{\parallel}=\frac{2n^{\parallel}}{a^{2}},$ (54)

while our RPY-based theory gives (51).

While we can set ${\hat{a}}$ to match the two translational theories or the
two rotational theories, there is no choice of ${\hat{a}}$ that will match
both. There are thus two possible candidates for ${\hat{a}}$. If we take

${\hat{a}}={\hat{a}}_{\text{tt}}=a\frac{e^{3/2}}{4}\approx 1.12a,$ (55)

then the two translational theories match [29, Appendix A], but for rotation
we obtain

$\displaystyle
8\pi\mu{\Psi}_{\text{n}}^{\parallel}=\frac{9n^{\parallel}(s)}{4{\hat{a}}_{\text{tt}}^{2}}=\frac{36n^{\parallel}(s)}{e^{3}a^{2}}\approx
1.79\frac{n^{\parallel}}{a^{2}},$ (56)

which is a 10% error from the known formula for a cylinder (50).

On the other hand, if we take

${\hat{a}}={\hat{a}}_{\text{rr}}=\sqrt{9/8}a\approx 1.06a,$ (57)

the rotational mobilities match, and the resulting RPY translational mobility
becomes

$\displaystyle 8\pi\mu\boldsymbol{U}$
$\displaystyle=\ln{\left(\frac{s(L-s)}{4{\hat{a}}_{\text{rr}}^{2}}\right)}\left(\boldsymbol{I}+\boldsymbol{\tau}\boldsymbol{\tau}\right)\boldsymbol{f}+4\boldsymbol{f}+\boldsymbol{U}_{\text{f}}^{\text{(FP)}}$
$\displaystyle\approx\ln{\left(\frac{4s(L-s)}{a^{2}}\right)}\left(\boldsymbol{I}+\boldsymbol{\tau}\boldsymbol{\tau}\right)\boldsymbol{f}+\left(1.1\boldsymbol{I}-2.9\boldsymbol{\tau}\boldsymbol{\tau}\right)\boldsymbol{f}+\boldsymbol{U}_{\text{f}}^{\text{(FP)}},$

which gives a small $\mathcal{O}(1)$ error relative to translational SBT (52).

####  3.5.2 Trans-rot coupling

The more compelling question is whether we can use our findings to gain
insights into the unknown $\mathcal{O}(1)$ coefficient in translation-rotation
SBT. Our RPY-based result is

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{n}}=\left(\ln{\left(\frac{s(L-s)}{4{\hat{a}}^{2}}\right)}+\frac{7}{6}\right)\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)\frac{n^{\parallel}}{2}+\boldsymbol{U}_{\text{f}}^{\text{(FP)}},$
(58)

with the corresponding symmetric result for $\Psi_{\text{f}}^{\parallel}$ in
(43). Since Section 3.5.1 shows that the RPY models capture the essential
physics of the interaction of the fiber with the fluid, but with different
coefficients, we can postulate that the rotation-translation coupling in a
true tubular geometry is also of the form (58), but with different
coefficients. Furthermore, since there is an $\mathcal{O}(1)$ error in the
boundary condition (3) in the result (10) of Johnson, we can expect that any
modifications to (10) will be $\mathcal{O}(1)$ (not $\log{\epsilon}$).
Combining these two assumptions, we postulate a rotation-translation coupling
term of the form

$8\pi\mu\boldsymbol{U}_{\text{n}}=\left(\ln{\left(\frac{4s(L-s)}{a^{2}}\right)}-k\right)\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)\frac{n^{\parallel}}{2}+\boldsymbol{U}_{\text{n}}^{\text{(FP)}}+\mathcal{O}(\epsilon)$
(59)

for some unknown constant $k$, with the corresponding symmetric form for
$\Psi_{\text{f}}^{\parallel}$. In the next section, we fit the constant $k$ to
three dimensional boundary integral results. This gives us an empirical three-
dimensional SBT for rot-trans coupling to which we can compare the RPY theory.

##  4 Comparing to three-dimensional boundary integral

In this section, we compare our slender body theories of Sections 2 and 3 to a
three-dimensional boundary integral (BI) calculation. Because there is one
uncertain coefficient, the $\mathcal{O}(1)$ coefficient $k$ in the rot-trans
coupling term (59), the goal of these calculations is to find a $k$ that can
predict the full BI results to $\mathcal{O}(\epsilon)$ accuracy. Because of
this, we can fix the fiber geometry, vary the aspect ratio $\epsilon$, and
check that the error in SBT decreases with $\epsilon$. It will be convenient
to fix $L=2$ and work on $s\in[-1,1]$ with an ellipsoidally-decaying radius
function

$\rho(s)=\sqrt{1-s^{2}}.$ (60)

The radius at each point $s$ is $a\rho(s)$, where we will vary $a$ to study
the limit $\epsilon=a/L\rightarrow 0$.

As discussed in Section 2.2, under the assumption that the fiber is an
infinitely thin membrane containing fluid of viscosity equal to that of the
exterior fluid, the force density multiplied by the surface area element on
the tubular boundary can be obtained by solving the integral equation

$\displaystyle
8\pi\mu\widehat{\boldsymbol{U}}(s,\theta)=\int_{-1}^{1}\int_{0}^{2\pi}\mathbb{S}\left(\widehat{\boldsymbol{X}}(s,\theta),\widehat{\boldsymbol{X}}\left(s^{\prime},\theta^{\prime}\right)\right)\widehat{\boldsymbol{f}}\left(s^{\prime},\theta^{\prime}\right)\,d\theta^{\prime}ds^{\prime},$
(61)

for $\widehat{\boldsymbol{f}}(s,\theta)$, where $\widehat{\boldsymbol{X}}$ and
$\widehat{\boldsymbol{U}}$ are defined in (1) and (2).

In order to compare the BI output to SBT, we use the following sequence of
steps:

  1. 1.

Given the surface velocity
$\widehat{\boldsymbol{U}}(s,\theta)=\boldsymbol{U}(s)+a\Psi^{\parallel}(s)\left(\boldsymbol{\tau}(s)\times\widehat{\boldsymbol{r}}(s,\theta)\right)$,
solve (61) using the numerical method of Section 4.1 to obtain the surface
force density $\widehat{\boldsymbol{f}}(s,\theta)$.

  2. 2.

Compute the centerline force $\boldsymbol{f}(s)$ and parallel torque
$n^{\parallel}(s)$ from $\widehat{\boldsymbol{f}}(s,\theta)$ via (14) and
(15).

  3. 3.

Compute the SBT velocity induced from the BI force as

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{SB}}=\ln{\left(\epsilon^{-2}\right)}\left(\boldsymbol{I}+\boldsymbol{\tau}\boldsymbol{\tau}\right)\boldsymbol{f}+\left(\boldsymbol{I}-3\boldsymbol{\tau}\boldsymbol{\tau}\right)\boldsymbol{f}+\boldsymbol{U}_{\text{f}}^{\text{(FP)}}$
(62)
$\displaystyle+\left(\ln{\left(\epsilon^{-2}\right)}-k\right)\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)\frac{n^{\parallel}}{2}+\boldsymbol{U}_{\text{n}}^{\text{(FP)}},$
$\displaystyle
8\pi\mu\Psi^{\parallel}_{\text{SB}}=\frac{2n^{\parallel}}{a^{2}\rho^{2}}+\frac{1}{2}\left(\ln{\left(\epsilon^{-2}\right)}-k\right)\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)\cdot\boldsymbol{f}+\Psi_{\text{f}}^{\text{(FP)}},$
(63)

where the finite part integrals are defined in (7), (11), and (44), and we
have used the local drag coefficients for ellipsoidally-decaying radius (60).

  4. 4.

Calculate the $L^{2}$ norm of the error
$\boldsymbol{U}_{\text{SB}}-\boldsymbol{U}$ and
$\Psi^{\parallel}_{\text{SB}}-\Psi^{\parallel}$, and plot it as a function of
$\epsilon$.

Note that we apply the forward SBT operator, which gives velocity from
force/torque, in (62) and (63). This avoids inverting the slender body
mobility, which is an ill-posed problem [14, 44]. To evaluate the finite part
integrals, we use the singular quadrature scheme developed in [43, 29, 30].

We have verified that the error in $\boldsymbol{U}_{\text{SB}}$ and
$\Psi^{\parallel}_{\text{SB}}$ decreases with $\epsilon$ in cases when rot-
trans coupling is negligible (i.e., cases when
$\Psi^{\parallel}_{\text{SB}}=0$ or the fiber is straight). Having done this,
our focus will be on the case of curved filaments with parallel rotational
velocity $\Psi^{\parallel}\sim 1/\epsilon^{2}$, so that the magnitude of the
parallel torque $n^{\parallel}$ is unchanged to leading order as we change
$\epsilon$. In this case, which is the one that exposed the flaws in the SBTs
of Section 2, the rot-trans coupling term in (62) is important, while the one
in (63) is negligible. Our goal is to find a $k$ in (62) such that the
formulas are accurate to $\mathcal{O}(\epsilon)$; that is, the relative errors
in both $\left\|\boldsymbol{U}_{\text{SB}}-\boldsymbol{U}\right\|$ and
$\Psi^{\parallel}_{\text{SB}}-\Psi^{\parallel}$ are $\mathcal{O}(\epsilon)$.

###  4.1 Numerical method

To solve (61), we need a way to integrate the singularity at
$\left(s^{\prime},\theta^{\prime}\right)=\left(s,\theta\right)$. We follow the
idea of Batchelor [2] and Koens [20] in adding and subtracting the flow from a
translating, rigid prolate ellipsoid, which has constant
$\widehat{\boldsymbol{f}}$ along its surface. The geometry of the prolate
ellipsoid is chosen to exactly cancel the singularity in (61). Specifically,
we introduce an ellipsoid with geometry

$\widehat{\boldsymbol{X}}_{e}\left(s^{\prime},\theta^{\prime},s,\theta\right)=\boldsymbol{X}(s)+\ell(s,\theta)\left(s^{\prime}-s_{c}(s)\right)\boldsymbol{\tau}(s)+a_{e}(s)\widehat{\boldsymbol{r}}(s,\theta^{\prime})\sqrt{1-\left(s^{\prime}\right)^{2}}.$
(64)

This ellipsoid has half axis length $\ell$ and radius $a_{e}$. The point
$s_{c}(s)$ is the coordinate on the ellipsoid at which we want to cancel the
singularity at $\widehat{\boldsymbol{X}}\left(s,\theta\right)$. Following
Koens [20] again, we cancel the singularity by imposing the conditions

$\displaystyle\widehat{\boldsymbol{X}}_{e}\left(s_{c},\theta^{\prime},s,\theta\right)=\widehat{\boldsymbol{X}}\left(s,\theta\right)$
(65)
$\displaystyle\partial_{s^{\prime}}\widehat{\boldsymbol{X}}_{e}\left(s_{c},\theta,s,\theta\right)=\partial_{s}\widehat{\boldsymbol{X}}(s,\theta)\qquad\partial_{\theta^{\prime}}\widehat{\boldsymbol{X}}_{e}\left(s_{c},\theta,s,\theta\right)=\partial_{\theta}\widehat{\boldsymbol{X}}(s,\theta).$

Substituting the definition of $\widehat{\boldsymbol{X}}$ from (1) and
$\widehat{\boldsymbol{X}}_{e}$ from (64), the matching conditions (65) give
three equations for the unknown parameters $\ell,a_{e}$, and $s_{c}$ in (64),

$\displaystyle a\rho=a_{e}\sqrt{1-s_{c}^{2}}$ (66) $\displaystyle
1+a\rho\left(\boldsymbol{\tau}\cdot\partial_{s}\widehat{\boldsymbol{r}}\right)=\ell$
(67) $\displaystyle a\partial_{s}\rho=-\frac{a_{e}s_{c}}{\sqrt{1-s_{c}^{2}}}.$
(68)

The first of these equations matches the positions and automatically matches
the $\theta$ derivatives, while the second two are obtained by projecting the
$s$ derivative equation onto the $\boldsymbol{\tau}$ and
$\widehat{\boldsymbol{r}}$ directions, respectively. To solve for
$s_{c},a_{e}$, and $\ell$, we first observe that (67) gives $\ell$ in terms of
the fiber geometry, and so we are left with a $2\times 2$ nonlinear system for
$a_{e}$ and $s_{c}$. This system has two solutions, only one of which gives
$s_{c}\in[-1,1]$,

$\displaystyle
s_{c}=\frac{\rho-\sqrt{\rho^{2}+4\left(\partial_{s}\rho\right)^{2}}}{2\rho},$
(69)

with $a_{e}$ determined from $s_{c}$ using (66). When the fiber is a straight
ellipsoid, $\rho(s)=\sqrt{1-s^{2}}$ and
$\partial_{s}\widehat{\boldsymbol{r}}=0$, so $a_{e}=a$ and $s_{c}=s$, and
therefore the effective ellipsoid reduces to the ellipsoid itself.

Because the traction multiplied by the area element on a translating prolate
ellipsoid is constant, the singularity subtraction scheme for (61) can now be
written as [20, Eq. 5]

$\displaystyle\widehat{\boldsymbol{U}}(s,\theta)=\boldsymbol{M}_{e}(s,\theta)\boldsymbol{f}(s,\theta)+\frac{1}{8\pi\mu}\int_{-1}^{1}\int_{0}^{2\pi}\Bigg{(}\mathbb{S}\left(\widehat{\boldsymbol{X}}(s,\theta),\widehat{\boldsymbol{X}}\left(s^{\prime},\theta^{\prime}\right)\right)\boldsymbol{f}\left(s^{\prime},\theta^{\prime}\right)$
(70)
$\displaystyle-\mathbb{S}\left(\widehat{\boldsymbol{X}}_{e}(s_{c},\theta,s,\theta),\widehat{\boldsymbol{X}}_{e}\left(s^{\prime},\theta^{\prime},s,\theta\right)\right)\boldsymbol{f}\left(s,\theta\right)\Bigg{)}\,d\theta^{\prime}ds^{\prime},$

where the matrix $\boldsymbol{M}_{e}(s,\theta)$ can be determined from
analytical formulas for the resistance on a prolate ellipsoid. In particular,
defining the eccentricity

$e=\sqrt{1-\frac{a_{e}^{2}}{\ell^{2}}},\qquad
q=\log{\left(\frac{1+e}{1-e}\right)},$ (71)

the resistance on a translating prolate ellipsoid is [18, Eqs. (44–45)] (see
also [35])

$\displaystyle\xi_{\parallel}=\frac{16\pi\ell
e^{3}\mu}{(1+e^{2})q-2e}\qquad\xi_{\perp}=\frac{32\pi\ell
e^{3}\mu}{(3e^{2}-1)q+2e}$ (72)

These give the relationship $F_{\parallel}=U_{\parallel}\xi_{\parallel}$ and
$F_{\perp}=U_{\perp}\xi_{\perp}$. To obtain the traction times area element on
the surface, we divide the resistance by $4\pi$ and invert the result, thus
obtaining an expression for the matrix $\boldsymbol{M}_{e}$,

$\displaystyle\boldsymbol{M}_{e}=\frac{4\pi}{\xi_{\parallel}}\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)+\frac{4\pi}{\xi_{\perp}}\left(\boldsymbol{I}-\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right).$
(73)

Because the singularity in (70) is canceled, we can now discretize the
integral by standard quadrature, with the singular point given a weight of
zero (i.e., skipped). We will use $N_{\theta}$ uniformly spaced points in
$\theta$ and $N$ Clenshaw-Curtis quadrature points in $s$. By assigning the
singular point a weight of zero, we reduce the accuracy of this quadrature
scheme, which is spectrally accurate for smooth integrands, to first order, as
demonstrated in Appendix C. Our Matlab implementation of this quadrature
scheme is available at
https://github.com/stochasticHydroTools/SlenderBody/tree/master/Matlab/NumericalSBT/SingleLayer.

###  4.2 Rotating curved fiber

To extract the unknown coefficient $k$ in the trans-rot coupling term in (62),
we consider the half-turn helix

$\boldsymbol{X}(s)=\frac{1}{\sqrt{2}}\left(\frac{1}{\pi}\cos{\left(\pi
s\right)},\frac{1}{\pi}\sin{\left(\pi s\right)},s\right),$ (74)

which spins with angular velocity $\Psi^{\parallel}\equiv 1/a^{2}$ and without
moving, $\boldsymbol{U}=0$ in (2). Because the fiber is not straight, the
torque required to generate $\Psi^{\parallel}$ will induce a translational
velocity, and a nonzero average force $\boldsymbol{f}$ must cancel this
translational velocity to satisfy the boundary condition (3). This means that
the translation-rotation coupling term in (62) will play a key role in the
accuracy of our SBT. Because of the large torque, the rotational velocity in
(63) is dominated by the rotation-rotation mobility, and so to study
translation-rotation coupling we will focus only on the translational velocity
(62).

In Appendix C, we verify first-order accuracy in obtaining the cross sectional
force $\boldsymbol{f}(s)$ in (14) from the traction jump
$\widehat{\boldsymbol{f}}(s,\theta)$ in (70) as we refine $N_{\theta}$ and $N$
_together_. We also show how the errors
$\left\|\boldsymbol{U}_{\text{SB}}-\boldsymbol{U}\right\|_{L^{2}}$, which are
obtained from the traction jump via (14) and then (62), change as we refine
$N_{\theta}$ and $N$ for three different values of $k$ ranging from $k=2$
(which corresponds to throwing out the second line of the rotlet line integral
(10)), to $k=3$. Since the values of
$\left\|\boldsymbol{U}_{\text{SB}}-\boldsymbol{U}\right\|_{L^{2}}$ do change
as we continue to refine our discretization, we report error bars in Fig. 1
that give the difference in $\boldsymbol{U}_{\text{SB}}$ between our last two
levels of refinement. Based on Fig. 4, the numerical error in the BI
calculations likely causes an underestimate of the SBT error for $k=2$ and
$k=2.5$, while for $k=2.85$ and $k=3$ we likely overestimate the SBT error.

Figure 1: Error in the rot-trans SBT relative to the boundary integral
calculation for various $\epsilon$ and $k$. We use the helical fiber (74) and
prescribe a rotational velocity $\Psi^{\parallel}\equiv 1/a^{2}$ with
translational velocity $\boldsymbol{U}=\boldsymbol{0}$. We then solve the
boundary integral equation (61) for the surface traction, use (14) and (15) to
obtain the centerline force density and parallel torque density, and then (62)
to obtain an SBT velocity $\boldsymbol{U}_{\text{SB}}$ from these densities.
We plot the $L^{2}$ norm of $\boldsymbol{U}_{\text{SB}}$ (which should be
zero) for several different values of $k$. For $k\neq 2.85$, the error is
clearly $\mathcal{O}(1)$ with respect to $\epsilon$, while for $k=2.85$ the
error decreases with $\epsilon$ at a rate of roughly $\epsilon^{1}$. The data
points here come from the most refined discretization in Fig. 2, while the
error bars are the difference between the most refined and second-most refined
calculation.

Looking at the SBT errors in Fig. 1, $k=2.85$ is the clear winner among the
possible choices of $k$. For $k=2$ (which is our first guess from the first
line of the rotlet integral expansion (10)), and $k=2.5$, we observe an error
which saturates to a constant as we decrease $\epsilon$, meaning that the
error in the SBT (62) when $k=2$ or $2.5$ is asymptotically $\mathcal{O}(1)$.
Increasing to $k=2.85$, we observe saturated errors that are roughly
proportional to $\epsilon$, which is exactly what we want to achieve in SBT.
When we increase again to $k=3$, we see errors which initially decrease with
$\epsilon$, but appear to plateau to an $\mathcal{O}(1)$ value. It appears
then, that there is an optimal $k$ around $k=2.85$ which makes the SBT formula
(62) accurate to $\mathcal{O}(\epsilon)$.

While this is by no means a proof that the correct rot-trans coupling
coefficient444Indeed, there is no _a priori_ reason to believe (62) or a
similarly simple equation holds rigorously from (2) and (13)–(15). in (62) is
$k=2.85$, it does provide us with a good educated guess that can be used to
cheaply approximate boundary integral simulations with small error, and get an
idea of how close our RPY formulas are to the true three-dimensional dynamics.
In particular, using an RPY radius of ${\hat{a}}\approx 1.86a$ in (58) gives
the formula (59) with $k=2.85$. This value of the RPY radius is comparable,
but fairly larger than, what we obtained by matching the trans-trans
(${\hat{a}}_{\text{tt}}\approx 1.12a$) and rot-rot
(${\hat{a}}_{\text{rr}}\approx 1.06a$) mobilities in Section 3.5.

##  5 Conclusion

We examined the deficiencies of existing slender body theories (SBTs) for
rotating filaments, and then proposed a semi-empirical SBT for tubular
filaments derived from line integrals of Rotne-Prager-Yamakawa (RPY)
singularities. We showed that the singularity-based SBTs of Keller+Rubinow and
Johnson [19, 15] both suffer from the same issues related to rotation-
translation coupling: when the fiber is driven by an applied torque at its
base, their theories imply that torque induces an $\mathcal{O}(1)$ angular-
dependent translational velocity on the fiber cross section, which violates
the SBT boundary condition that the fiber only translates and rotates. In
addition, these previous SBTs exhibit a lack of symmetry in the $2\times 2$
mobility operator because of assumptions about the relative order of magnitude
of the fiber translational and rotational velocities. As we showed previously
[30, Sec. 5.3], neglecting rotation-translation coupling can lead to
overestimation of the torque required to induce large-scale filament
deformations. It is therefore important to derive an SBT which properly treats
rotation-translation coupling.

As a workaround to the true problem of obtaining an SBT for a three-
dimensional cylindrical (or prolate-spherioidal) filament, we used a different
definition of a filament as a line of infinitely many regularized
singularities, similar to what has been done previously by Cortez and Nicholas
[10]. Each component of the $2\times 2$ mobility operator (translation-
translation, rotation-translation, etc.) is then defined as a line integral of
the corresponding fundamental regularized singularity, which are defined in
(23). Because these regularized singularities, which define the mobility for a
pair of “blobs,” give a symmetric $2\times 2$ mobility, our mobility for
fibers is automatically symmetric, which alleviates one of the two major
problems with previous SBTs. The other issue, the $\mathcal{O}(1)$ angular
dependence of the translational velocity on a cross section, is sidestepped
because we use an isotropic regularized singularity function (a delta function
on a sphere). This comes at a price: by using a series of spheres to represent
the fiber centerline, we lose fidelity to the true three-dimensional geometry.
While this is far from ideal, it gives a mobility which is physical and simple
to compute, since all that results is a local drag term plus a finite part
integral for the trans-trans, rot-trans, and trans-rot terms, for which
special quadrature schemes have already been developed [43].

Using matched asymptotics on the RPY kernels, we found a $2\times 2$ mobility
relationship, which we then compared to SBT for a three-dimensional tubular
fiber. For translation-translation and rotation-rotation, we found that the
two mobilities take the same form, but with an $\mathcal{O}(1)$ difference in
the constants. The discrepancy in the constants can be removed via a judicious
choice of regularized singularity radius, which is different for each
component of the mobility; ${\hat{a}}\approx 1.06a$ gives a perfect match in
the rot-rot mobilities, while ${\hat{a}}\approx 1.12a$ matches the trans-trans
mobilities. We used the correspondence between RPY asymptotics and SBT to
posit the relationship (59) for the off-diagonal rotation-translation
components in SBT. The relationship has a single unknown coefficient $k\approx
2.85$, which we obtained in Section 4 from boundary integral calculations. The
value of $k$ we found corresponds to a regularized singularity radius
${\hat{a}}\approx 1.86a$.

Using the RPY kernels simplified a number of potential issues that have
plagued the development of SBT. First, we were able to easily derive mobility
formulas near the fiber endpoints by modifying the domain of integration in
the inner expansion (see Appendix B). The resulting endpoint mobilities given
in Appendix B are _different_ from classical SBT, where the fiber radius
function _must_ decay ellipsoidally at the endpoints, with the mobility
transitioning to a constant near the fiber boundaries. In other words, the RPY
approach gives a physical regularization of the local drag coefficients near
the fiber endpoints, which is in contrast to other _ad hoc_ approaches [29,
Sec. 2.1]. Furthermore, the RPY mobility immediately generalizes to multiple
filaments, since the flow induced on filament $j$ can also be computed by an
integral of the RPY kernel over filament $i$. In SBT, while the pointwise
velocity induced by singularities might be constant on a cross section of
fiber $i$, the same might not be true on a nearby fiber $j$, which leads to an
ambiguity in how the induced flow on the centerline of fiber $j$ is assigned.

Our boundary integral calculations in Section 4 can be viewed as both a check
on the accuracy of the rotation-translation coupling terms in our RPY theory,
and also as a semi-empirical SBT for translation-rotation coupling. Indeed,
these calculations are an attempt to find an SBT by reversing the process
Koens and Lauga [22] used to show that translational SBT can be derived from
the BI equation. While we showed that this approach is unworkable analytically
for the rotation-translation coupling problem, numerical calculation allowed
us to extract the relationship

$8\pi\mu\boldsymbol{U}_{\text{n}}\approx\left(\ln{\left(\epsilon^{-2}\right)}-2.85\right)\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)\frac{n^{\parallel}}{2}+\boldsymbol{U}_{\text{n}}^{\text{(FP)}}+\mathcal{O}(\epsilon)$
(75)

for the translational velocity due to parallel torque on an ellipsoidal
filament, with the finite part integral
$\boldsymbol{U}_{\text{n}}^{\text{(FP)}}$ defined in (11). Of course, this
result is an estimate for a particular filament, and there is no way to show
it is universal, although tests on a one-turn and two-turn helix gave similar
results. Still, our guess for the rotation-translation coupling tensor is
informed by the asymptotics of the RPY tensor, and so, if it is incorrect or
non-universal, the RPY mobility (58) must miss some of the essential physics
in the modeling of slender filaments. Given the number of studies that have
used the RPY tensor to model filament hydrodynamics (e.g., [47, 41, 17, 48]),
this would be concerning.

Our work here is by no means the final say on SBT for twisting filament, but
is rather intended as a starting point to provoke future research in three
main directions. First, there is the possibility of developing a singularity
representation for the flow around a twisting filament that properly matches
the velocity on the fiber cross section. For an ellipsoidal filament, the
velocity on a cross section due to a rotlet (10) can be expanded using (12) as

$\displaystyle
8\pi\mu\boldsymbol{u}^{(\Psi)}\left(\widehat{\boldsymbol{X}}\left(s,\theta\right)\right)=\frac{2n^{\parallel}}{a}\left(\boldsymbol{\tau}\times\hat{\boldsymbol{r}}\right)+\frac{n^{\parallel}}{2}\left(\ln{\left(\epsilon^{-2}\right)}-3\right)\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)+\boldsymbol{U}_{\text{n}}^{\text{(FP)}}$
(76)
$\displaystyle+\frac{n^{\parallel}\kappa}{2}\left(\boldsymbol{e}_{b}\cos{2\theta}-\boldsymbol{e}_{n}\sin{2\theta}\right)+\mathcal{O}\left(\epsilon\right).$

Thus, the task is to find another singularity (or set of singularities) that
can eliminate the $\mathcal{O}(1)$ dependence on $2\theta$ in the second line,
thereby making the cross-sectional velocity a pure rotation (the
$\boldsymbol{\tau}\times\widehat{\boldsymbol{r}}$ term) plus a pure
translation (the next two terms in the first line). Based on our empirical SBT
(75), we expect such a singularity, if it exists, to eliminate the angular
dependence while introducing a small correction to the $\mathcal{O}(1)$ term
multiplying $\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}$.

Given that it is difficult to find such a singularity representation, if it
exists at all, our empirical approach in Section 4 offers some promise, if
only we could reach even smaller $\epsilon$. In the case when
$\Psi^{\parallel}\sim 1/\epsilon^{2}$, a more accurate numerical approach is
necessary to do so, since in that case the surface tractions scale like
$1/\epsilon$, and it becomes difficult to confirm $\mathcal{O}(\epsilon)$
scaling of asymptotic expressions. And finally, there is an opportunity to
analyze what Stokes-related partial differential equation such a slender body
theory is the solution of (this is an extension of the work of Mori, Ohm, and
others for translation [32, 31, 33]).

## Acknowledgments

We thank Ajeet Kumar and Mohit Garg for useful discussions that led to this
work. This work was supported by the National Science Foundation through
Research Training Group in Modeling and Simulation under award RTG/DMS-1646339
and through the Division of Mathematical Sciences award DMS-2052515.

##  Appendix A Isotropic single layer expansion

To expand the isotropic single layer (18), we begin with an inner expansion of
the displacement
$\boldsymbol{R}(s,s^{\prime},\theta,\theta^{\prime})=\widehat{\boldsymbol{X}}(s,\theta)-\widehat{\boldsymbol{X}}(s^{\prime},\theta^{\prime})$,

$\displaystyle\boldsymbol{R}(s,s^{\prime},\theta,\theta^{\prime})=\widehat{\boldsymbol{X}}(s,\theta)-\widehat{\boldsymbol{X}}(s^{\prime},\theta^{\prime})=\boldsymbol{X}(s)-\boldsymbol{X}(s^{\prime})+\hat{\boldsymbol{e}}_{\rho}(s,\theta)-\hat{\boldsymbol{e}}_{\rho}(s^{\prime},\theta^{\prime}),$
(A.1)

where $\hat{\boldsymbol{e}}_{\rho}(s,\theta)$ is the unit vector going from
$\boldsymbol{X}(s)$ to $\widehat{\boldsymbol{X}}(s,\theta)$. Since the
geometry is more complicated in this true three-dimensional case, it is
helpful to parameterize $\hat{\boldsymbol{e}}_{\rho}$ by an angle
$\theta_{i}(s)$ satisfying
$\partial_{s}\theta_{i}=\boldsymbol{e}_{b}\cdot\partial_{s}\boldsymbol{e}_{n}$
(the curve torsion) [22, Eq. (2.2)]. This gives [22, Sec. 2]

$\displaystyle\hat{\boldsymbol{e}}_{\rho}(s,\theta)=\boldsymbol{e}_{n}\cos{(\theta-\theta_{i}(s))}+\boldsymbol{e}_{b}\sin{(\theta-\theta_{i}(s))}$
(A.2)
$\displaystyle\partial_{s}\hat{\boldsymbol{e}}_{\rho}(s,\theta)=-\kappa(s)\cos{(\theta-\theta_{i})}\boldsymbol{\tau}(s)$

Substituting into (A.1), we can now write the first two terms in the Taylor
expansion of $\boldsymbol{R}$ and $1/R$ as

$\displaystyle\boldsymbol{R}\left(\widehat{\boldsymbol{X}}(s,\theta),\widehat{\boldsymbol{X}}(s^{\prime},\theta^{\prime})\right)\approx
a\left(-\xi\boldsymbol{\tau}+\Delta\hat{\boldsymbol{e}}_{\rho}(s,\theta,\theta^{\prime})\right)+a^{2}\left(\frac{\xi^{2}}{2}\partial_{s}\boldsymbol{\tau}-\xi\kappa\cos{(\theta^{\prime}-\theta_{i})}\boldsymbol{\tau}\right)$
$\displaystyle\frac{1}{R\left(\widehat{\boldsymbol{X}}(s,\theta),\widehat{\boldsymbol{X}}(s^{\prime},\theta^{\prime})\right)}=\frac{1}{a\sqrt{\xi^{2}+\Delta\hat{\boldsymbol{e}}_{\rho}\cdot\Delta\hat{\boldsymbol{e}}_{\rho}}}+\frac{-\xi^{2}\kappa\cos{(\theta^{\prime}-\theta_{i})}+\xi^{2}/2\boldsymbol{\tau}^{\prime}\cdot\Delta\hat{\boldsymbol{e}}_{\rho}}{2\sqrt{\xi^{2}+\Delta\hat{\boldsymbol{e}}_{\rho}\cdot\Delta\hat{\boldsymbol{e}}_{\rho}}^{3}}$
(A.3)

where
$\Delta\hat{\boldsymbol{e}}_{\rho}(s,\theta,\theta^{\prime})=\hat{\boldsymbol{e}}_{\rho}(s,\theta)-\hat{\boldsymbol{e}}_{\rho}(s,\theta^{\prime})$
and $\xi=(s^{\prime}-s)/a$ is order 1. Note the appearance of the centerline
curvature (scalar $\kappa$ and vector $\partial_{s}\boldsymbol{\tau}$) at this
order, which is vital to yield proper rot-trans coupling. Combining (A.3) with
the Taylor series representation
$\boldsymbol{f}(s^{\prime},\theta^{\prime})=\boldsymbol{f}(s,\theta^{\prime})+a\xi\partial_{s}\boldsymbol{f}(s,\theta^{\prime})$,
we have the inner expansion of the integrand in (16) as

$\displaystyle\frac{\boldsymbol{f}(s^{\prime},\theta^{\prime})}{R}=\frac{\boldsymbol{f}(s,\theta^{\prime})}{a\sqrt{\xi^{2}+\Delta\hat{\boldsymbol{e}}_{\rho}\Delta\hat{\boldsymbol{e}}_{\rho}}}+\frac{\xi\partial_{s}\boldsymbol{f}(s,\theta^{\prime})}{\sqrt{\xi^{2}+\Delta\hat{\boldsymbol{e}}_{\rho}\Delta\hat{\boldsymbol{e}}_{\rho}}}$
(A.4)
$\displaystyle+\frac{-\xi^{2}\kappa\cos{(\theta^{\prime}-\theta_{i})}+\xi^{2}/2(\partial_{s}\boldsymbol{\tau}\cdot\Delta\hat{\boldsymbol{e}}_{\rho})}{2\sqrt{\xi^{2}+\Delta\hat{\boldsymbol{e}}_{\rho}\Delta\hat{\boldsymbol{e}}_{\rho}}^{3}}\boldsymbol{f}(s,\theta^{\prime})+\mathcal{O}(\epsilon).$

The second line already makes it clear how rot-trans coupling is going to
arise. If $\boldsymbol{f}(s,\theta^{\prime})=\boldsymbol{f}(s)$ (no angular
dependence), then there will still be a term in the velocity proportional to
$\cos{\theta}\boldsymbol{f}(s)$, which is rotational velocity coming from
translational force, or rot-trans coupling.

To complete the expansion, we proceed with the integration of (A.4) from
$\xi=-s/a$ to $\xi=(L-s)/a$. The integration can be done exactly, and then
expanded in $\epsilon$ to yield simpler expressions [22, Table 1]. After using
(A.2) to determine
$\partial_{s}\boldsymbol{\tau}\cdot\Delta\hat{\boldsymbol{e}}_{\rho}=0$, then
integrating what remains, we have

$\displaystyle\boldsymbol{U}_{i1}(s,\theta)=\int_{0}^{2\pi}\Bigg{(}\ln{\left(\frac{4s(L-s)}{a^{2}(\Delta\hat{\boldsymbol{e}}_{\rho}\cdot\Delta\hat{\boldsymbol{e}}_{\rho})}\right)}\boldsymbol{f}(s,\theta^{\prime})+\left(L-2s\right)\partial_{s}\boldsymbol{f}(s,\theta)$
(A.5)
$\displaystyle-\frac{\epsilon\kappa\cos{\left(\theta^{\prime}-\theta_{i}(s)\right)}\boldsymbol{f}(s,\theta^{\prime})}{2}\left(\ln{\left(\frac{4s(L-s)}{a^{2}(\Delta\hat{\boldsymbol{e}}_{\rho}\cdot\Delta\hat{\boldsymbol{e}}_{\rho})}\right)}-2\right)\Bigg{)}\,d\theta^{\prime}+\mathcal{O}(\epsilon^{2}).$

We substitute the one-term Fourier series representation (17) into (A.5) and
integrate over $\theta^{\prime}$ to obtain (18). To integrate, we use
$\Delta\hat{\boldsymbol{e}}_{\rho}\cdot\Delta\hat{\boldsymbol{e}}_{\rho}=2\left(1-\cos{\left(\theta-\theta^{\prime}\right)}\right)$.

##  Appendix B Inner expansions for RPY asymptotics

This appendix is devoted to the inner expansions of the RPY velocities (22).
Using the rescaled variable $\xi$ defined in (24), we write the series
expansions

$\displaystyle\boldsymbol{R}(s^{\prime})=\boldsymbol{X}(s)-\boldsymbol{X}(s^{\prime})=-\xi{\hat{a}}\boldsymbol{\tau}(s)-\frac{\xi^{2}{\hat{a}}^{2}}{2}\partial_{s}\boldsymbol{\tau}(s)+\mathcal{O}\left({\hat{a}}^{3}\right)$
(B.1)
$\displaystyle\boldsymbol{R}\boldsymbol{R}=\xi^{2}{\hat{a}}^{2}\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)+\mathcal{O}\left({\hat{a}}^{3}\right),$
$\displaystyle
R^{2}=\boldsymbol{R}\cdot\boldsymbol{R}=\xi^{2}{\hat{a}}^{2}+\mathcal{O}\left({\hat{a}}^{3}\right),\quad
R=|\xi|{\hat{a}}+\mathcal{O}\left({\hat{a}}^{2}\right),$ $\displaystyle
R^{-1}=\frac{1}{|\xi|{\hat{a}}}+\mathcal{O}(1),\quad
R^{-3}=\frac{1}{|\xi|^{3}{\hat{a}}^{3}}+\mathcal{O}\left({\hat{a}}^{-2}\right),\quad
R^{-5}=\frac{1}{|\xi|^{5}{\hat{a}}^{5}}+\mathcal{O}\left({\hat{a}}^{-4}\right).$
$\displaystyle\boldsymbol{f}\left(s^{\prime}\right)=\boldsymbol{f}(s)+\xi{\hat{a}}\partial_{s}\boldsymbol{f}(s)+\mathcal{O}({\hat{a}}^{2})\qquad
n^{\parallel}\left(s^{\prime}\right)=n^{\parallel}(s)+\xi{\hat{a}}\partial_{s}n^{\parallel}(s)+\mathcal{O}({\hat{a}}^{2}),$

which we substitute into the corresponding RPY velocity formulas (26), (34),
(42), and (47). What results is two separate integrals over $\xi$: one for
$|\xi|>2$ (when the RPY spheres do not overlap), and another for $|\xi|\leq 2$
(when they do). The bounds on the integrals are

$\displaystyle\int_{|\xi|>2}f(\xi)\,d\xi=\begin{cases}\int_{-s/{\hat{a}}}^{-2}f(\xi)\,d\xi+\int_{2}^{(L-s)/{\hat{a}}}f(\xi)\,d\xi&2{\hat{a}}<s<L-2{\hat{a}}\\\
\int_{2}^{(L-s)/{\hat{a}}}f(\xi)\,d\xi&s\leq 2{\hat{a}}\\\
\int_{-s/{\hat{a}}}^{-2}f(\xi)\,d\xi&s\geq L-2{\hat{a}}\\\ \end{cases}$ (B.2)
$\displaystyle\int_{|\xi|\leq
2}f(\xi)\,d\xi=\begin{cases}\int_{-2}^{-2}f(\xi)\,d\xi&2{\hat{a}}<s<L-2{\hat{a}}\\\
\int_{-s/{\hat{a}}}^{2}f(\xi)\,d\xi&s\leq 2{\hat{a}}\\\
\int_{-2}^{(L-s)/{\hat{a}}}f(\xi)\,d\xi&s\geq L-2{\hat{a}}\\\ \end{cases}$

With these integration bounds, the inner velocity can be computed for any $s$
on the filament, as we do next.

###  B.1 Translation from force

Beginning with the trans-trans mobility (26), we start with the part of the
integral in the region $R>2{\hat{a}}$. For this we first need to integrate the
inner expansion of the Stokeslet

$\displaystyle\int_{R>2{\hat{a}}}\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)$
$\displaystyle\boldsymbol{f}\left(s^{\prime}\right)\,ds^{\prime}\approx\int_{|\xi|>2}\left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|\xi|{\hat{a}}}\right)\boldsymbol{f}(s){\hat{a}}\,d\xi$
(B.3)
$\displaystyle=\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)\int_{|\xi|>2}\frac{1}{|\xi|}d\xi.$

Unlike the Stokeslet, the doublet term is only included in the inner
expansion. Its expansion is given by

$\displaystyle\frac{2{\hat{a}}^{2}}{3}$
$\displaystyle\int_{R>2{\hat{a}}}\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\,ds^{\prime}\approx\frac{2{\hat{a}}^{2}}{3}\int_{|\xi|>2}\left(\frac{\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|\xi|^{3}{\hat{a}}^{3}}\right)\boldsymbol{f}(s){\hat{a}}\,d\xi$
(B.4)
$\displaystyle=\frac{2\left(\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)}{3}\int_{|\xi|>2}\frac{1}{|\xi|^{3}}\,d\xi.$

It still remains to include in the inner expansion the term for $R\leq
2{\hat{a}}$. For this we have the inner expansion of the first term

$\displaystyle\int_{R\leq
2{\hat{a}}}\left(\frac{4}{3{\hat{a}}}-\frac{3R\left(s^{\prime}\right)}{8{\hat{a}}^{2}}\right)\boldsymbol{f}\left(s^{\prime}\right)\,ds^{\prime}$
$\displaystyle\approx\boldsymbol{f}(s)\int_{|\xi|\leq
2}\left(\frac{4}{3}-\frac{3|\xi|}{8}\right)\,d\xi.$ (B.5)

Likewise the second term for $R\leq 2{\hat{a}}$ has the inner expansion

$\displaystyle\int_{R\leq
2{\hat{a}}}\frac{1}{8{\hat{a}}^{2}R\left(s^{\prime}\right)}\left(\boldsymbol{R}\boldsymbol{R}\right)\left(s^{\prime}\right)\boldsymbol{f}\left(s^{\prime}\right)\,ds^{\prime}$
$\displaystyle\approx\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\boldsymbol{f}\left(s\right)\int_{|\xi|\leq
2}\frac{|\xi|}{8}\,d\xi.$ (B.6)

Adding the terms (B.3), (B.4), (B.5), and (B.6), and computing integrals using
the bounds (B.2), we obtain the complete inner expansion

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{f}}^{(\text{inner})}(s)=\left(a_{L}(s)\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)+a_{I}(s)\boldsymbol{I}+a_{\tau}(s)\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)$
(B.7) $\displaystyle
a_{L}(s)=\begin{cases}\ln{\left(\dfrac{(L-s)s}{4{\hat{a}}^{2}}\right)}&2{\hat{a}}<s<L-2{\hat{a}}\\\\[4.0pt]
\ln{\left(\dfrac{(L-s)}{2{\hat{a}}}\right)}&s\leq 2{\hat{a}}\\\\[4.0pt]
\ln{\left(\dfrac{s}{2{\hat{a}}}\right)}&s\geq L-2{\hat{a}}\end{cases}$
$\displaystyle
a_{I}(s)=\begin{cases}4-\dfrac{{\hat{a}}^{2}}{3s^{2}}-\dfrac{{\hat{a}}^{2}}{3(L-s)^{2}}&2{\hat{a}}<s<L-2{\hat{a}}\\\\[4.0pt]
2+\dfrac{4s}{3{\hat{a}}}-\dfrac{3s^{2}}{16{\hat{a}}^{2}}-\dfrac{{\hat{a}}^{2}}{3(L-s)^{2}}&s\leq
2{\hat{a}}\\\\[4.0pt]
2+\dfrac{4(L-s)}{3{\hat{a}}}-\dfrac{3(L-s)^{2}}{16{\hat{a}}^{2}}-\dfrac{{\hat{a}}^{2}}{3s^{2}}&s\geq
L-2{\hat{a}}\end{cases}a_{\tau}(s)=\begin{cases}\dfrac{{\hat{a}}^{2}}{s^{2}}+\dfrac{{\hat{a}}^{2}}{(L-s)^{2}}&2{\hat{a}}<s<L-2{\hat{a}}\\\\[4.0pt]
\dfrac{s^{2}}{16{\hat{a}}^{2}}+\dfrac{{\hat{a}}^{2}}{(L-s)^{2}}&s\leq
2{\hat{a}}\\\\[4.0pt]
\dfrac{(L-s)^{2}}{16{\hat{a}}^{2}}+\dfrac{{\hat{a}}^{2}}{s^{2}}&s\geq
L-2{\hat{a}}\end{cases}$

The velocity $\boldsymbol{U}_{\text{f}}^{(\text{inner})}$ is twice
continuously differentiable, since it is the integral of a function which is
once continuously differentiable.

###  B.2 Rot-trans and trans-rot coupling

We now move on to the inner expansion in the case of translation-rotation
coupling (34). Beginning with the rotlet term for $R>2{\hat{a}}$, we write the
rotlet in terms of the inner variables as

$\displaystyle\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})$
$\displaystyle\approx\frac{\xi^{2}{\hat{a}}^{2}}{2}\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)$
(B.8)
$\displaystyle\int_{R>2{\hat{a}}}\frac{\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})}{R(s^{\prime})^{3}}n^{\parallel}(s^{\prime})\,ds^{\prime}$
$\displaystyle\approx\int_{|\xi|>2}\frac{\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)}{2|\xi|{\hat{a}}}n^{\parallel}(s){\hat{a}}\,d\xi.$
(B.9)

We next need to account for the term for $R\leq 2{\hat{a}}$, or the second
term in (34). Using (B.8), we have that

$\displaystyle\frac{1}{2{\hat{a}}^{2}}\int_{R<2{\hat{a}}}\left(\frac{1}{{\hat{a}}}-\frac{3R}{8{\hat{a}}^{2}}\right)\left(\boldsymbol{n}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)\,ds^{\prime}$
$\displaystyle\approx\frac{\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)}{4{\hat{a}}}\int_{|\xi|<2}\left(1-\frac{3}{8}|\xi|\right)\xi^{2}{\hat{a}}\,d\xi.$
(B.10)

Adding (B.9) and (B.10) and using (B.2) to compute the inner expansion
integrals, we obtain

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{n}}^{\text{(inner)}}(s)=\frac{n^{\parallel}(s)}{2}\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)c_{c}(s)$
(B.11) $\displaystyle
c_{c}(s)=\begin{cases}\ln{\left(\dfrac{(L-s)s}{4{\hat{a}}^{2}}\right)}+\dfrac{7}{6}&2{\hat{a}}<s<L-2{\hat{a}}\\\\[4.0pt]
\ln{\left(\dfrac{(L-s)}{2{\hat{a}}}\right)}+\dfrac{7}{12}+\dfrac{s^{3}}{6{\hat{a}}^{3}}-\dfrac{3s^{4}}{64{\hat{a}}^{4}}&s\leq
2{\hat{a}}\\\\[4.0pt]
\ln{\left(\dfrac{s}{2{\hat{a}}}\right)}+\dfrac{7}{12}+\dfrac{(L-s)^{3}}{6{\hat{a}}^{3}}-\dfrac{3(L-s)^{4}}{64{\hat{a}}^{4}}&s\geq
L-2{\hat{a}}\end{cases}$

which is a twice continuously differentiable function.

The inner expansion of angular velocity from force is given symmetrically as

$\displaystyle
8\pi\mu\Psi_{\text{f}}^{\parallel,\text{(inner)}}(s)=\frac{1}{2}c_{c}(s)\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)\cdot\boldsymbol{f}(s).$
(B.12)

###  B.3 Rotation from torque

The inner expansion for (47) contains two parts: the doublet, and terms for
$R\leq 2{\hat{a}}$. We have already seen the expansion of the doublet in
(B.4). Multiplying the result in (B.4) by the appropriate prefactors, we
obtain

$\displaystyle-\frac{1}{2}\int_{R>2{\hat{a}}}\left(\frac{\boldsymbol{n}(s^{\prime})}{R(s^{\prime})^{3}}-3\frac{\boldsymbol{R}(s^{\prime})\left(\boldsymbol{R}(s^{\prime})\cdot\boldsymbol{n}(s^{\prime})\right)}{R(s^{\prime})^{5}}\right)\,ds^{\prime}\approx-\frac{\left(\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)}{2{\hat{a}}^{2}}\int_{|\xi|>2}\frac{1}{|\xi|^{3}}\,d\xi,$
(B.13)

which which is accurate to $\mathcal{O}\left(\hat{\epsilon}^{-1}\right)$,
although in the fiber interior that term actually cancels, leaving an
$\mathcal{O}\left(\log{\hat{\epsilon}^{-1}}\right)$ dependence.

In a similar way, the inner expansions of the $R\leq 2{\hat{a}}$ terms in the
fiber interior are

$\displaystyle\frac{1}{{\hat{a}}^{3}}\int_{R<2{\hat{a}}}\left(1-\frac{27R(s^{\prime})}{32{\hat{a}}}+\frac{5R(s^{\prime})^{3}}{64{\hat{a}}^{3}}\right)\boldsymbol{n}(s^{\prime})\,ds^{\prime}\approx\frac{\boldsymbol{n}(s)}{{\hat{a}}^{2}}\int_{|\xi|\leq
2}\left(1-\frac{27}{32}|\xi|+\frac{5}{64}|\xi|^{3}\right)\,d\xi$ (B.14)
$\displaystyle\frac{1}{{\hat{a}}^{3}}\int_{R<2{\hat{a}}}\left(\left(\frac{9}{32aR(s^{\prime})}-\frac{3R(s^{\prime})}{64{\hat{a}}^{3}}\right)\left(\boldsymbol{R}\boldsymbol{R}\right)(s^{\prime})\right)\boldsymbol{n}(s^{\prime})\,ds^{\prime}\approx\frac{\boldsymbol{\tau}(s)\left(\boldsymbol{\tau}(s)\cdot\boldsymbol{n}(s)\right)}{{\hat{a}}^{2}}\int_{|\xi|\leq
2}\left(\frac{18-3|\xi|^{2}}{64}\right)|\xi|\,d\xi$

with an error of $\mathcal{O}(\ln{\hat{\epsilon}})$ in the fiber interior and
$\mathcal{O}(\hat{\epsilon}^{-1})$ at the endpoints.

Using the integral bounds (B.2) to integrate and sum (B.13) and (B.14), the
complete inner expansion for $\boldsymbol{\Psi}_{\text{n}}$ is given by

$\displaystyle
8\pi\mu\boldsymbol{\Psi}_{\text{n}}^{\text{(inner)}}(s)=\left(p_{I}(s)\boldsymbol{I}+p_{\tau}(s)\boldsymbol{\boldsymbol{\tau}}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{n}(s),\qquad\text{where}$
(B.15) $\displaystyle
p_{I}(s)=\dfrac{1}{{\hat{a}}^{2}}\begin{cases}\dfrac{9}{8}+\dfrac{{\hat{a}}^{2}}{4}\left(\dfrac{1}{s^{2}}+\dfrac{1}{(L-s)^{2}}\right)&2{\hat{a}}<s<L-2{\hat{a}}\\\\[6.0pt]
\dfrac{9}{16}+\dfrac{s}{{\hat{a}}}-\dfrac{27s^{2}}{64{\hat{a}}^{2}}+\dfrac{5s^{4}}{256{\hat{a}}^{4}}+\dfrac{{\hat{a}}^{2}}{4(L-s)^{2}}&s\leq
2{\hat{a}}\\\\[6.0pt]
\dfrac{9}{16}+\dfrac{L-s}{{\hat{a}}}-\dfrac{27(L-s)^{2}}{64{\hat{a}}^{2}}+\dfrac{5(L-s)^{4}}{256{\hat{a}}^{4}}+\dfrac{{\hat{a}}^{2}}{4s^{2}}&s\geq
L-2{\hat{a}}\end{cases}$ $\displaystyle
p_{\tau}(s)=\dfrac{1}{{\hat{a}}^{2}}\begin{cases}\dfrac{9}{8}-\dfrac{3{\hat{a}}^{2}}{4}\left(\dfrac{1}{s^{2}}+\dfrac{1}{(L-s)^{2}}\right)&2{\hat{a}}<s<L-2{\hat{a}}\\\\[6.0pt]
\dfrac{9}{16}+\dfrac{9s^{2}}{64{\hat{a}}^{2}}-\dfrac{3s^{4}}{256{\hat{a}}^{4}}-\dfrac{3{\hat{a}}^{2}}{4(L-s)^{2}}&s\leq
2{\hat{a}}\\\\[6.0pt]
\dfrac{9}{16}+\dfrac{9(L-s)^{2}}{64{\hat{a}}^{2}}-\dfrac{3(L-s)^{4}}{256{\hat{a}}^{4}}-\dfrac{3{\hat{a}}^{2}}{4s^{2}}&s\geq
L-2{\hat{a}}\end{cases}$

which is twice continuously differentiable.

##  Appendix C Convergence of quadrature scheme

In this appendix, we discuss the convergence of the quadrature scheme for (70)
for the rotating curved fiber in Section 4.2. We report convergence as a
function of $N_{t}=NN_{\theta}$ the total number of points on the filament
boundary. In order to obtain smooth convergence, both $N$ and $N_{\theta}$
must be increased together. To maintain a similar spacing in the axial and
circumferential directions, we use $N_{\theta}=16,20,\dots 32$, and $N\approx
0.25N_{\theta}/\epsilon$.

We first perform a self-convergence study in Fig. 2, where we look at the
$L^{2}$ errors in $\boldsymbol{f}$ and $n^{\parallel}$, computed from (14) and
(15), under refinement. For the force, the error at a given $N_{t}$ is defined
as the $L^{2}$ difference of $\boldsymbol{f}$ when compared to the next
largest $N_{t}$, and similarly for $n^{\parallel}$. As $N_{t}$ increases, the
force is asymptotically first-order accurate for every $\epsilon$. The torque
converges at a faster rate because, in this example, torque is dominated by
local rotation-rotation dynamics.

Figure 2:  Convergence plot for the curved rotating filament studied in
Section 4.2. We show the self-convergence of our numerical method for (left)
the force $\boldsymbol{f}$ and (right) the parallel torque $n^{\parallel}$ on
each cross section. The “error” is the $L^{2}$ difference of each quantity
relative to the next level of refinement, normalized by the $L^{2}$ norm of
the most refined solution.  Figure 3:  SBT asymptotic errors for rotational
velocity in the fully coupled problem. We show how the relative error
$a^{2}\left(\Psi^{\parallel}_{\text{SB}}-\Psi^{\parallel}\right)$ changes as
we change $\epsilon$ and the total number of discretization points $N_{t}$.
The SBT rotational velocity $\Psi^{\parallel}_{\text{SB}}$ is computed from
(63) and is dominated by rotation-rotation dynamics. It therefore does not
depend on $k$, and so we show only a single $k=2.85$.

In Figs. 3 and 4, we show how the convergence of the force and torque affects
the convergence of the SBT errors, $\boldsymbol{U}_{\text{SB}}-\boldsymbol{U}$
and $\Psi^{\parallel}_{\text{SB}}-\Psi^{\parallel}$ for various $\epsilon$ and
$k$. As we increase $N_{t}$, our hope is that the errors will saturate as we
isolate the SBT asymptotic error from the boundary integral numerical error.
For rotational velocity, which is dominated by torque and therefore
independent of the $k$ used in (63), this is indeed the case, as we see
constant errors under refinement in Fig. 3. Furthermore, the relative error in
$\Psi^{\parallel}_{\text{SB}}$ scales roughly as $\epsilon^{2}\log{\epsilon}$,
which is expected since the rotational velocity is dominated by the rot-rot
term.

For translational velocity, we do not observe saturation of
$\left\|\boldsymbol{U}_{\text{SB}}-\boldsymbol{U}\right\|$ in Fig. 4, but we
can still obtain a good estimate of the SBT error by looking at our finest
discretizations. In particular, the errors for $k=2$ and $k=2.5$ appear to be
approaching a constant of about 0.2 and 0.1, respectively, irrespective of the
particular value of $\epsilon$. This indicates an $\mathcal{O}(1)$ error in
the SBT asymptotic equation (62) for those values of $k$. However, when we
increase $k$ further to 2.85 and later 3, the errors are much smaller,
indicating that the SBT theory gets more accurate. To account for our BI
solutions not being fully converged, we report error bars in Fig. 1 that give
the difference between our most refined and second-most refined errors in
$\boldsymbol{U}_{\text{SB}}$ (the difference between the rightmost and second-
rightmost point on each colored curve in Fig. 4).

Figure 4:  SBT asymptotic errors for translational velocity in the fully
coupled problem. We show how the error in the slender body velocity
$\boldsymbol{U}_{\text{SB}}$ changes as we refine the number of discretization
points and the fiber aspect ratio. Here the velocity
$\boldsymbol{U}_{\text{SB}}$ is computed using (62) with the value of $k$
indicated in the legend. The bottom plot shows the error
$a^{2}\left(\Psi^{\parallel}_{\text{SB}}-\Psi^{\parallel}\right)$, which is
computed from (63) and is dominated by rotation-rotation dynamics. It
therefore does not depend on $k$, and so we show only a single $k=2.85$ at the
bottom.

## References

  * [1] Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, and Peter Walter.  Molecular biology of the cell.  Garland Science, 2002.
  * [2] GK Batchelor.  Slender-body theory for particles of arbitrary cross-section in stokes flow.  Journal of Fluid Mechanics, 44(3):419–440, 1970.
  * [3] GK Batchelor.  The stress system in a suspension of force-free particles.  Journal of fluid mechanics, 41(3):545–570, 1970.
  * [4] Howard C Berg.  E. coli in Motion.  Springer, 2004.
  * [5] Neeraj S Borker and Donald L Koch.  Slender body theory for particles with non-circular cross-sections with application to particle dynamics in shear flows.  Journal of Fluid Mechanics, 877:1098–1133, 2019.
  * [6] Thomas T Bringley and Charles S Peskin.  Validation of a simple method for representing spheres and slender bodies in an immersed boundary method for stokes flow on an unbounded domain.  Journal of Computational Physics, 227(11):5397–5425, 2008.
  * [7] Xiaobing Chen and Howard C Berg.  Torque-speed relationship of the flagellar rotary motor of escherichia coli.  Biophysical journal, 78(2):1036–1041, 2000.
  * [8] Ricardo Cortez.  The method of regularized stokeslets.  SIAM Journal on Scientific Computing, 23(4):1204–1225, 2001.
  * [9] Ricardo Cortez, Lisa Fauci, and Alexei Medovikov.  The method of regularized stokeslets in three dimensions: analysis, validation, and application to helical swimming.  Physics of Fluids, 17(3):031504, 2005.
  * [10] Ricardo Cortez and Michael Nicholas.  Slender body theory for stokes flows with regularized forces.  Communications in Applied Mathematics and Computational Science, 7(1):33–62, 2012.
  * [11] Frédéric Eghiaian, Annafrancesca Rigato, and Simon Scheuring.  Structural, mechanical, and dynamical variability of the actin cortex in living cells.  Biophysical journal, 108(6):1330–1340, 2015.
  * [12] Andrew M Fiore and James W Swan.  Rapid sampling of stochastic displacements in brownian dynamics simulations with stresslet constraints.  The Journal of chemical physics, 148(4):044114, 2018.
  * [13] Mohit Garg and Ajeet Kumar.  A slender body theory for the motion of special cosserat filaments in stokes flow.
  * [14] Thomas Götz.  Interactions of fibers and flow: asymptotics, theory and numerics.  dissertation. de, 2001.
  * [15] Robert E Johnson.  An improved slender-body theory for stokes flow.  Journal of Fluid Mechanics, 99(2):411–431, 1980.
  * [16] Robert Edward Johnson.  Slender-body theory for Stokes flow and flagellar hydrodynamics.  PhD thesis, California Institute of Technology, 1977.
  * [17] Bakytzhan Kallemov, Amneet Bhalla, Boyce Griffith, and Aleksandar Donev.  An immersed boundary method for rigid bodies.  Communications in Applied Mathematics and Computational Science, 11(1):79–141, 2016.
  * [18] Eric E Keaveny and Michael J Shelley.  Applying a second-kind boundary integral equation for surface tractions in stokes flow.  Journal of Computational Physics, 230(5):2141–2159, 2011.
  * [19] Joseph B Keller and Sol I Rubinow.  Slender-body theory for slow viscous flow.  Journal of Fluid Mechanics, 75(4):705–714, 1976.
  * [20] Lyndon Koens.  Tubular-body theory for viscous flows.  Physical Review Fluids, 7(3):034101, 2022.
  * [21] Lyndon Koens and Eric Lauga.  Slender-ribbon theory.  Physics of Fluids, 28(1):013101, 2016.
  * [22] Lyndon Koens and Eric Lauga.  The boundary integral formulation of stokes flows includes slender-body theory.  Journal of Fluid Mechanics, 850, 2018.
  * [23] Eric Lauga and Thomas R Powers.  The hydrodynamics of swimming microorganisms.  Reports on Progress in Physics, 72(9):096601, 2009.
  * [24] Wanho Lee, Yongsam Kim, Sarah D Olson, and Sookkyung Lim.  Nonlinear dynamics of a rotating elastic rod in a viscous fluid.  Physical Review E, 90(3):033012, 2014.
  * [25] Sookkyung Lim and Charles S Peskin.  Simulations of the whirling instability by the immersed boundary method.  SIAM Journal on Scientific Computing, 25(6):2066–2083, 2004.
  * [26] Leroy F Liu and James C Wang.  Supercoiling of the dna template during transcription.  Proceedings of the National Academy of Sciences, 84(20):7024–7027, 1987.
  * [27] Sune Lomholt and Martin R Maxey.  Force-coupling method for particulate two-phase flow: Stokes flow.  J. Comput. Phys., 184(2):381–405, 2003.
  * [28] MR Maxey and BK Patel.  Localized force representations for particles sedimenting in stokes flow.  International journal of multiphase flow, 27(9):1603–1626, 2001\.
  * [29] Ondrej Maxian, Alex Mogilner, and Aleksandar Donev.  Integral-based spectral method for inextensible slender fibers in stokes flow.  Physical Review Fluids, 6(1):014102, 2021.
  * [30] Ondrej Maxian, Brennan Sprinkle, Charles S Peskin, and Aleksandar Donev.  The hydrodynamics of a twisting, bending, inextensible fiber in stokes flow.  arXiv preprint arXiv:2201.04187, 2022.
  * [31] Yoichiro Mori and Laurel Ohm.  Accuracy of slender body theory in approximating force exerted by thin fiber on viscous fluid.  Studies in Applied Mathematics, 2021.
  * [32] Yoichiro Mori, Laurel Ohm, and Daniel Spirn.  Theoretical justification and error analysis for slender body theory.  Communications on Pure and Applied Mathematics, 73(6):1245–1314, 2020.
  * [33] Yoichiro Mori, Laurel Ohm, and Daniel Spirn.  Theoretical justification and error analysis for slender body theory with free ends.  Archive for Rational Mechanics and Analysis, 235(3):1905–1978, 2020\.
  * [34] Laurel Ohm.  Remarks on regularized stokeslets in slender body theory.  Fluids, 6(8):283, 2021.
  * [35] Juan C Padrino, James E Sprittles, and Duncan A Lockerby.  Comment on “applying a second-kind boundary integral equation for surface tractions in stokes flow”.  Journal of Computational Physics, 401:109007, 2020.
  * [36] Charles S Peskin.  Flow patterns around heart valves: a numerical method.  J. Comput. Phys., 10(2):252–271, 1972.
  * [37] Charles S Peskin.  The immersed boundary method.  Acta Numer., 11:479–517, 2002.
  * [38] Thomas R Powers.  Dynamics of filaments and membranes in a viscous fluid.  Reviews of Modern Physics, 82(2):1607, 2010.
  * [39] Constantine Pozrikidis et al.  Boundary integral and singularity methods for linearized viscous flow.  Cambridge university press, 1992.
  * [40] Jens Rotne and Stephen Prager.  Variational treatment of hydrodynamic interaction in polymers.  J. Chem. Phys., 50(11):4831–4837, 1969.
  * [41] Simon F Schoeller, Adam K Townsend, Timothy A Westwood, and Eric E Keaveny.  Methods for suspensions of passive and active filaments.  Journal of Computational Physics, 424:109846, 2021.
  * [42] DK Srivastava.  Slender body theory for stokes flow past axisymmetric bodies: a review article.  International Journal of Applied Mathematics and Mechanics, 8(15):14–39, 2012.
  * [43] Anna-Karin Tornberg.  Accurate evaluation of integrals in slender-body formulations for fibers in viscous flow.  arXiv preprint arXiv:2012.12585, 2020.
  * [44] Anna-Karin Tornberg and Michael J Shelley.  Simulating the dynamics and interactions of flexible fibers in stokes flows.  Journal of Computational Physics, 196(1):8–40, 2004.
  * [45] Hirofumi Wada and Roland R Netz.  Non-equilibrium hydrodynamics of a rotating filament.  EPL (Europhysics Letters), 75(4):645, 2006.
  * [46] Eligiusz Wajnryb, Krzysztof A Mizerski, Pawel J Zuk, and Piotr Szymczak.  Generalization of the rotne–prager–yamakawa mobility and shear disturbance tensors.  Journal of Fluid Mechanics, 731, 2013.
  * [47] Jianghui Wang, Emilio J Tozzi, Michael D Graham, and Daniel J Klingenberg.  Flipping, scooping, and spinning: Drift of rigid curved nonchiral fibers in simple shear flow.  Physics of Fluids, 24(12):123304, 2012.
  * [48] Timothy A Westwood and Eric E Keaveny.  Coordinated motion of active filaments on spherical surfaces.  Physical Review Fluids, 6(12):L121101, 2021.
  * [49] Charles W Wolgemuth, Thomas R Powers, and Raymond E Goldstein.  Twirling and whirling: Viscous dynamics of rotating elastic filaments.  Physical Review Letters, 84(7):1623, 2000.
  * [50] Junhua Yuan, Karen A Fahrner, Linda Turner, and Howard C Berg.  Asymmetry in the clockwise and counterclockwise rotation of the bacterial flagellar motor.  Proceedings of the National Academy of Sciences, 107(29):12846–12849, 2010.
  * [51] Boan Zhao and Lyndon Koens.  Regularized stokeslets lines suitable for slender bodies in viscous flow.  Fluids, 6(9):335, 2021.