fiber begins a steady crankshafting rigid rotation around the $y$ axis, such
that $Y(s)$ and $\sqrt{X(s)^{2}+Z(s)^{2}}$ are constant for all $s$. In Fig.
8, we show how the amplitude $R=\sqrt{X(L)^{2}+Z(L)^{2}}$ (Fig. 8(b)) and
frequency $\chi$ (inset of Fig. 8(a)) of the crankshafting motion depend on
the driving motor frequency. For $\omega/\omega_{c}\lesssim 1.6$, the
amplitude contracts rapidly to about 5% of the fiber length, while the
crankshafting frequency is roughly constant at $1-2$% of the motor frequency.
On $\omega/\omega_{c}\gtrsim 1.6$, however, the fiber amplitude function is
roughly constant; indeed, Fig. 8(c) shows that the steady state overwhirling
shapes approach a fiber that is essentially twirling upside down in this case.
Instead of bringing the endpoint closer to the $y$ axis, the extra work done
by the motor is used to increase the frequency of the crankshafting motion at
a higher rate, so that it transitions to about 10% of the motor frequency. The
mean $\chi/\omega$ of about 0.05 is similar to what Wada and Netz found for a
thermally fluctuating fiber [70, Fig. 1(c)]. Our frequency results are also
similar to those presented previously in [39, Fig. 6(a)], where nonlinear
behavior is observed for $\omega/\omega_{c}\gtrsim 1.6$.

Interestingly, while we find that the overwhirling state can be stable for
frequencies as large as $3\omega_{c}$, we also find that the ability to reach
it depends on the initial fiber configuration. For instance, when
$\omega/\omega_{c}$ is as small as 2.1, the fiber intersects itself on its way
to the overwhirling state. An animation of this intersection, which is more
severe for larger $\omega/\omega_{c}$, is shown when $\omega/\omega_{c}=2.8$
in one of the supplementary movies. This suggests a mechanism for the
formation of plectonemes that were observed at high $\omega$ in [13], and that
the full nonlinear dynamics are likely more complex than previously thought.

##  6 Conclusions

In this paper, we studied a single, inextensible filament that resists bending
and twisting while immersed in a viscous fluid. Our focus was on adding twist
elasticity to the model that we used in previous work [48] for inextensible
filaments with bending resistance. We did this in two parts: first, we
compared two possible models for twist elasticity in the presence of a fluid,
which we called the Kirchhoff and Euler models. While the two models have the
same elastic constituitive laws, in the Kirchhoff model the angular velocity
of the material frame is constrained to be equal to the angular fluid
velocity, while in the Euler model the perpendicular components of the angular
fluid velocity are discarded. The contribution of the first part of this paper
was to reformulate the Kirchhoff model in a form similar to the Euler model,
and use this to show that these two models produce dynamics (parallel and
perpendicular angular velocities) that differ by $\mathcal{O}(\epsilon^{2})$
in the fiber interior, where $\epsilon$ is the fiber aspect ratio. We chose to
utilize the Euler model in the remainder of the paper, since the extra degree
of rotational freedom (parallel torque and parallel rotational velocity) can
be easily incorporated into the formulation we introduced in [48]. To compute
the additional force and parallel torque due to twist, we used the Bishop
(twist-free) frame and evolved the angle of twist relative to that frame [8].

We defined the hydrodynamic mobility operators for translation-translation,
translation-rotation, rotation-translation, and rotation-rotation coupling in
terms of integrals of the corresponding Rotne-Prager-Yamakawa (RPY) kernels,
then evaluated these integrals both asymptotically and numerically. Our
preference for trans-trans, trans-rot, and rot-trans was to use the slender
body _quadrature_ schemes described in Appendix G to evaluate the integrals;
we showed that this removes the spurious negative eigenvalues that result when
the trans-trans integral is evaluated asymptotically (as in slender body
_theory_). The only remaining uncontrolled approximation in our mobility
evaluation is to use local drag, and not the full RPY integral, for the rot-
rot mobility. We did this because the $\mathcal{O}(\epsilon^{2})$ accuracy of
this approximation was sufficient for our purposes. If necessary, the special
quadrature scheme we developed for the doublet (in the context of trans-trans
mobility) could also be applied to the rot-rot kernel to account for the
nonlocal rot-rot hydrodynamics to higher accuracy. Using the RPY kernels to
account for the fiber thickness undoubtedly gives some error relative to the
true three dimensional geometry. Nevertheless, the equations we wrote are
intended as a mathematical _model_ of a thin elastic rod in fluid, in much the
same spirit as the Kirchhoff rod model itself, which is not derived by solving
any three-dimensional elasticity problem but which aims to capture the essence
of a thin rod’s elastic behavior in a one-dimensional self-consistent
continuum model. In the same way, our formulation of the fluid-structure
interaction equations for a thin rod in fluid does not exactly correspond to
the solution of any 3D Stokes problem with an internal moving tubular boundary
[36], but instead we work with the variables of the Kirchhoff rod model and
couple them to the fluid in a physically reasonable way.

We applied our numerical method to two common single-filament problems in
which twist elasticity has been traditionally neglected and included,
respectively. We considered a relaxing bent filament and showed that twist
elasticity has a negligible $\mathcal{O}(\epsilon^{2})$ effect, which
justifies the assumption in many studies [51, 48] to neglect it. By contrast,
the spinning filament problem requires finite twist elasticity to drive the
filament into an unstable trajectory for spinning frequencies
$\omega>\omega_{c}$. We showed that our numerical method reproduces exactly
the semi-numerical prediction of [75] for $\omega_{c}$ when using
ellipsoidally-tapered filaments with local drag, and that the shape of the
filament and the addition of nonlocal trans-trans hydrodynamics modify the
critical frequency by at most 5% [70]. We showed that rot-trans coupling is
more significant, lowering the frequency by another 15% to about 80% of the
value predicted for ellipsoidally-tapered local drag, which explains about
half of the deviation observed experimentally [13].

Above the critical frequency $\omega_{c}$, the fiber transitions to a steady
crankshafting rotation in the shape of a shepherd’s crook (overwhirling) [41].
While we found this state to be stable up to the highest frequencies we
investigated, getting to that state from an initially straight filament
requires complex dynamics that can lead to self-intersections, suggesting that
in practice plectonemes may be observed instead [13]. Our numerical method
cannot handle self intersections, both because of the quadrature scheme and
because we do not include lubrication, frictional, or steric forces [10].

Despite the length and complexity of this paper, we only considered a single
filament! Yet it is quite interesting how many modeling and numerical choices
we had to make for this seemingly simple problem. For instance, we chose _not_
to regularize the shape of the filament, keeping the radius constant all the
way up to the endpoints. Because we use integrals of the RPY tensor, the
mobility for constant regularized radius ${\hat{a}}$ _is_ well-defined at the
fiber endpoints, although its lack of smoothness causes slow convergence in
our spectral method. This is in contrast to SBT, in which the local drag
coefficient is not well-defined for cylindrical fibers, and some kind of
ellipsoidally-tapered radius function has to be assumed [48, Sec. 2.1]. Future
work could develop a hybrid approach that uses the RPY kernel for _unequally-
sized spheres_ [77] to simulate fibers with radius functions that decay
smoothly to zero at the endpoints. This should both smooth out the mobility
and keep the eigenvalues of the trans-trans mobility away from zero. The
integrals involved are more complex than those for equally sized spheres
because the integrand includes the radius function ${\hat{a}}(s)$. Still,
because the singularities (Stokeslet, doublet, rotlet, etc.) in the
polydisperse integral kernels are unchanged, it should be possible to account
for non-constant ${\hat{a}}(s)$ using the slender body quadrature schemes
developed here.

Another choice we made for the single filament was to stick with a first kind
integral equation for the force-velocity relationship, i.e., the trans-trans
mobility. While this yielded sufficient accuracy for our purposes, the
eigenvalues of the first kind mobility cluster around zero, leading to ill-
defined solutions near the fiber endpoints. Of course, one argument around
this (which we made in [48]) is that both the ill-conditioning and negative
eigenvalues of the SBT mobility occur when $N\gtrsim 1/\hat{\epsilon}$, which
is too large for slender body approaches to be efficient. Still, a better-
conditioned formulation (for trans-trans mobility) would be to use a second-
kind integral equation of the form given in [5, Eq. (4)]. To enable rapid
dynamic simulation of the model in [5], it is necessary to develop slender
body quadratures to evaluate integrals of regularized Stokeslets and doublets
along the centerline, as we did for the first kind RPY mobility in Appendix G.

An important future extension of this work is to consider the dynamics of
multiple filaments in a suspension with twist, as we did without twist in [48,
49]. The mobility of Section 2 directly generalizes to multiple filaments,
since the same RPY kernel can be used between distinct filaments and a
filament and itself. The challenge lies in efficiently numerically computing
the solution. Focusing on the static case first, we need to develop an
efficient numerical method to sum the rot-trans and trans-rot hydrodynamic
interactions over many bodies, since our previous work [48, Sec. 4] already
treated the trans-trans kernel with a combination of the spectral Ewald method
[23] (for well-separated fibers)666We chose an Ewald approach because we are
interested in triply-periodic systems; in free space a fast multipole method
[27] is likely the best choice, and systems with mixed periodicity are an
avenue of active research [42, 76]. and corrections for nearly-touching fibers
[1]. Extending our many-body-summation approach to twist is straightforward,
but arduous. In particular, spectral Ewald methods have been developed for the
rot-trans and trans-rot singularities of Stokes flow [24, 2], and the nearly-
singular quadratures of [1] readily generalize to the other singularities. For
rot-rot coupling, the fast decay of the doublet singularity makes it a good
approximation to only include the flows generated by nearly contacting fibers.

Then there is the question of computing the dynamics of many twisting
filaments. In our previous SBT-based work [48], we treated terms involving
nonlocal hydrodynamic interactions explicitly, and only treated the local drag
terms implicitly. This allowed us to solve essentially diagonal linear
systems, and to evaluate long-ranged hydrodynamics only once or at most a few
times per time step. In the slender body quadrature approach followed here,
there is not, at first sight, an explicit splitting into a “local drag”
contribution and a nonlocal contribution, as there is in traditional SBT. But
in the trans-trans mobility, we used singularity subtraction to factor out the
leading order singularity in the RPY integral. This, combined with the
integrals for $R<2{\hat{a}}$, give a “local drag” contribution that can be
treated implicitly. The rot-rot mobility is even more dominated by local
behavior, and so we can treat the angular velocity coming from torque on any
other fibers explicitly and only treat the self contribution implicitly. A
challenge for future exploration is how to efficiently handle pairs of
filaments (or pieces of a single filament) in near contact in a semi-implicit
temporal integrator, without expensive GMRES iterations [48, Secs. 4.5.2, 5.2]
at each time step.

The spectral numerical method developed in this work is efficient, but lacks
the simplicity, flexibility, and robustness of second-order “blob-based”
methods that use $\mathcal{O}(\hat{\epsilon}^{-1})$ discrete regularized
singularities along the centerline. Spectral methods excel especially when
blob-based discretizations are expensive (because of small aspect ratios) or
unnecessary (because of smooth fiber shapes). They are more challenging,
however, when the fiber shapes are nonsmooth, as it was only possible to model
overwhirling with a great effort in anti-aliasing. A related even greater
challenge for the future is to account for thermal fluctuations of elastic
filaments in spectral approaches.

## Acknowledgments

We thank Bjorn Dunkel and Vishal Prakash Patil for helpful discussions
involving the Bishop frame. 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 RPY mobility

When the regularized delta function $\delta_{\hat{a}}$ is a surface Delta
function on a sphere of radius ${\hat{a}}$ (i.e.,
$\delta_{\hat{a}}(\boldsymbol{r})=1/(4\pi{\hat{a}}^{2})\delta_{\hat{a}}(r-{\hat{a}})$),
the kernels
$\mathbb{M}_{\text{tt}},\mathbb{M}_{\text{tr}},\mathbb{M}_{\text{rt}}$, and
$\mathbb{M}_{\text{rr}}$ in (6) and (9) can be written analytically, and are
collectively known as the Rotne-Prager-Yamakawa (RPY) tensor. For the
translational velocity (6) at position $\boldsymbol{x}$, they take the form
[72]

$\displaystyle\left(8\pi\mu\right)\,\mathbb{M}_{\text{tt}}\left(\boldsymbol{x},\boldsymbol{y}\right)\boldsymbol{F}=\begin{cases}\left(\dfrac{\boldsymbol{I}+\widehat{\boldsymbol{R}}\widehat{\boldsymbol{R}}}{R}+\dfrac{2{\hat{a}}^{2}}{3}\dfrac{\boldsymbol{I}-3\widehat{\boldsymbol{R}}\widehat{\boldsymbol{R}}}{R^{3}}\right)\boldsymbol{F}&R>2{\hat{a}}\\\\[8.0pt]
\left(\left(\dfrac{4}{3{\hat{a}}}-\dfrac{3R}{8{\hat{a}}^{2}}\right)\boldsymbol{I}+\dfrac{R}{8{\hat{a}}^{2}}\widehat{\boldsymbol{R}}\widehat{\boldsymbol{R}}\right)\boldsymbol{F}&R\leq
2{\hat{a}}\end{cases}$ (A.1)
$\left(8\pi\mu\right)\,\mathbb{M}_{\text{tr}}\left(\boldsymbol{x},\boldsymbol{y}\right)\boldsymbol{N}=\begin{cases}\dfrac{\boldsymbol{N}\times\widehat{\boldsymbol{R}}}{R^{2}}&R>2{\hat{a}}\\\\[8.0pt]
\dfrac{1}{2{\hat{a}}^{2}}\left(\dfrac{R}{{\hat{a}}}-\dfrac{3R^{2}}{8{\hat{a}}^{2}}\right)\left(\boldsymbol{N}\times\widehat{\boldsymbol{R}}\right)&R\leq
2{\hat{a}}\end{cases}$ (A.2)

where $\boldsymbol{F}$ and $\boldsymbol{N}$ are the strength of a regularized
force and torque, respectively, applied at position $\boldsymbol{y}$, the
displacement vector $\boldsymbol{R}=\boldsymbol{x}-\boldsymbol{y}$,
$R=\left\|\boldsymbol{R}\right\|$,
$\widehat{\boldsymbol{R}}=\boldsymbol{R}/\left\|\boldsymbol{R}\right\|$, and
$\widehat{\boldsymbol{R}}\widehat{\boldsymbol{R}}$ denotes the outer product
of $\widehat{\boldsymbol{R}}$ with itself. Note that the RPY tensor changes
form for the case of overlapping spheres ($R\leq 2{\hat{a}}$). At
$R=2{\hat{a}}$, the mobilities are once continuously differentiable.

Likewise, for the angular velocity, the double convolutions in (9) can be
computed analytically, which yields the formulas [72]

$\displaystyle\left(8\pi\mu\right)\,\mathbb{M}_{\text{rt}}\left(\boldsymbol{x},\boldsymbol{y}\right)\boldsymbol{F}=\begin{cases}\dfrac{\boldsymbol{F}\times\widehat{\boldsymbol{R}}}{R^{2}}&R>2{\hat{a}}\\\\[8.0pt]
\dfrac{1}{2{\hat{a}}^{2}}\left(\dfrac{R}{{\hat{a}}}-\dfrac{3R^{2}}{8{\hat{a}}^{2}}\right)\left(\boldsymbol{F}\times\widehat{\boldsymbol{R}}\right)&R\leq
2{\hat{a}}\end{cases}$ (A.3)
$\displaystyle\left(8\pi\mu\right)\,\mathbb{M}_{\text{rr}}\left(\boldsymbol{x},\boldsymbol{y}\right)\boldsymbol{N}=\begin{cases}-\dfrac{1}{2}\left(\dfrac{\boldsymbol{I}-3\widehat{\boldsymbol{R}}\widehat{\boldsymbol{R}}}{R^{3}}\right)\boldsymbol{N}&R>2{\hat{a}}\\\\[2.0pt]
\dfrac{1}{{\hat{a}}^{3}}\left(\left(1-\dfrac{27R}{32{\hat{a}}}+\dfrac{5R^{3}}{64{\hat{a}}^{3}}\right)\boldsymbol{I}+\left(\dfrac{9R}{32{\hat{a}}}-\dfrac{3R^{3}}{64{\hat{a}}^{3}}\right)\widehat{\boldsymbol{R}}\widehat{\boldsymbol{R}}\right)\boldsymbol{N}&R\leq
2{\hat{a}}.\end{cases}$ (A.4)

##  Appendix B Derivation of twist evolution equation

In this appendix, we derive the evolution equation (24) for the twist $\psi$.
We present this derivation because it uses straightforward calculus, and not
the variational approach in the review article [58].

Because the material frame
$\left(\boldsymbol{D}^{(1)}(s),\boldsymbol{D}^{(2)}(s),\boldsymbol{D}^{(3)}(s)\right)$
is orthonormal, it evolves by the parallel transport equation

$\partial_{s}\boldsymbol{D}^{(i)}=\boldsymbol{\omega}\times\boldsymbol{D}^{(i)},$
(B.1)

which also gives

$\displaystyle\psi$
$\displaystyle=\partial_{s}\boldsymbol{D}^{(1)}\cdot\boldsymbol{D}^{(2)}$
$\displaystyle=\left(\boldsymbol{\omega}\times\boldsymbol{D}^{(1)}\right)\cdot\boldsymbol{D}^{(2)}=\left(\boldsymbol{D}^{(1)}\times\boldsymbol{D}^{(2)}\right)\cdot\boldsymbol{\omega}=\boldsymbol{\tau}\cdot\boldsymbol{\omega}$
(B.2)

So $\omega^{\parallel}=\psi$. Likewise, setting $i=3$ in (B.1) and using the
constraint $\boldsymbol{\tau}=\boldsymbol{D}^{(3)}$ from (13) gives
$\partial_{s}\boldsymbol{\tau}=\boldsymbol{\omega}\times\boldsymbol{\tau}$.
Crossing both sides with $\boldsymbol{\tau}$, we obtain the total
representation of $\boldsymbol{\omega}$ as [8]

$\boldsymbol{\omega}=\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}+\boldsymbol{\tau}\psi$
(B.3)

We now equate the $t$ derivative of (B.1) with the $s$ derivative of (22),

$\displaystyle\partial_{t}\left(\boldsymbol{\omega}\times\boldsymbol{D}^{(i)}\right)=\partial_{s}\left(\boldsymbol{\Omega}\times\boldsymbol{D}^{(i)}\right)$
$\displaystyle\partial_{t}\boldsymbol{\omega}\times\boldsymbol{D}^{(i)}+\boldsymbol{\omega}\times\partial_{t}\boldsymbol{D}^{(i)}=\partial_{s}\boldsymbol{\Omega}\times\boldsymbol{D}^{(i)}+\boldsymbol{\Omega}\times\partial_{s}\boldsymbol{D}^{(i)}$
$\displaystyle\partial_{t}\boldsymbol{\omega}\times\boldsymbol{D}^{(i)}=\partial_{s}\boldsymbol{\Omega}\times\boldsymbol{D}^{(i)}+\boldsymbol{\Omega}\times\left(\boldsymbol{\omega}\times\boldsymbol{D}^{(i)}\right)-\boldsymbol{\omega}\times\left(\boldsymbol{\Omega}\times\boldsymbol{D}^{(i)}\right).$
(B.4)

The last equality used the evolution equations (B.1) and (22) for
$\boldsymbol{D}^{(i)}$. Using the triple cross product identity
$\boldsymbol{a}\times\left(\boldsymbol{b}\times\boldsymbol{c}\right)+\boldsymbol{c}\times\left(\boldsymbol{a}\times\boldsymbol{b}\right)+\boldsymbol{b}\times\left(\boldsymbol{c}\times\boldsymbol{a}\right)=\boldsymbol{0}$,
we can write (B.4) as

$\displaystyle\partial_{t}\boldsymbol{\omega}\times\boldsymbol{D}^{(i)}=\partial_{s}\boldsymbol{\Omega}\times\boldsymbol{D}^{(i)}+\left(\boldsymbol{\Omega}\times\boldsymbol{\omega}\right)\times\boldsymbol{D}^{(i)},$
(B.5)

which must hold for all $i$. This gives [58, Eq. (54)]

$\partial_{t}\boldsymbol{\omega}=\partial_{s}\boldsymbol{\Omega}+\left(\boldsymbol{\Omega}\times\boldsymbol{\omega}\right).$
(B.6)

To derive the evolution equation for the twist angle, we take the
$\boldsymbol{\tau}$ component of (B.6) and use (22)

$\displaystyle\partial_{t}\left(\boldsymbol{\omega}\cdot\boldsymbol{\tau}\right)-\boldsymbol{\omega}\cdot\partial_{t}\boldsymbol{\tau}=\partial_{s}\left(\boldsymbol{\Omega}\cdot\boldsymbol{\tau}\right)-\boldsymbol{\Omega}\cdot\partial_{s}\boldsymbol{\tau}+\left(\boldsymbol{\Omega}\times\boldsymbol{\omega}\right)\cdot\boldsymbol{\tau}$
(B.7)
$\displaystyle\partial_{t}\left(\boldsymbol{\omega}\cdot\boldsymbol{\tau}\right)=\boldsymbol{\omega}\cdot\left(\boldsymbol{\Omega}\times\boldsymbol{\tau}\right)+\partial_{s}\left(\boldsymbol{\Omega}\cdot\boldsymbol{\tau}\right)-\boldsymbol{\Omega}\cdot\partial_{s}\boldsymbol{\tau}+\left(\boldsymbol{\Omega}\times\boldsymbol{\omega}\right)\cdot\boldsymbol{\tau}.$

The triple products cancel, and using (B.2), we get the twist evolution
equation (24),

$\displaystyle\partial_{t}\psi=\partial_{s}\Omega^{\parallel}-\left(\boldsymbol{\Omega}\cdot\partial_{s}\boldsymbol{\tau}\right).$

##  Appendix C Euler model for intrinsic curvature and twist

In Sections 2.3 and 2.4, we saw that the Euler model can be derived from the
Kirchhoff rod model by first setting
$\boldsymbol{F}=\boldsymbol{F}^{(\kappa)}+\boldsymbol{F}^{(\gamma)}+\boldsymbol{\Lambda}$
in (33), where $\boldsymbol{F}^{(\kappa)}$ and $\boldsymbol{F}^{(\gamma)}$ are
defined to eliminate the perpendicular torque in $\partial_{s}\boldsymbol{N}$,
then setting $\boldsymbol{\Lambda}=T\boldsymbol{\tau}$, which eliminates the
perpendicular torque Lagrange multipliers that are the distinguishing feature
of the Kirchhoff rod model. These two steps eliminate _all_ perpendicular
torques from the formulation, which leads to the Euler model in Section 2.4.

In this appendix, we repeat this procedure on a fiber with intrinsic curvature
and twist (e.g., a flagellum) to derive the corresponding Euler model. In this
case, we have the moment [54]

$\displaystyle\boldsymbol{N}$
$\displaystyle=\kappa\left(\left(\partial_{s}\boldsymbol{D}^{(2)}\cdot\boldsymbol{D}^{(3)}-\kappa_{1}\right)\boldsymbol{D}^{(1)}+\left(\partial_{s}\boldsymbol{D}^{(3)}\cdot\boldsymbol{D}^{(1)}-\kappa_{2}\right)\boldsymbol{D}^{(2)}\right)+\gamma\left(\partial_{s}\boldsymbol{D}^{(1)}\cdot\boldsymbol{D}^{(2)}-\phi\right)\boldsymbol{D}^{(3)}$
$\displaystyle=\kappa\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}-\kappa_{1}\boldsymbol{D}^{(1)}-\kappa_{2}\boldsymbol{D}^{(2)}\right)+\gamma\left(\psi-\phi\right)\boldsymbol{\tau},$
(C.1)

where $\kappa_{1}$ and $\kappa_{2}$ are the preferred curvatures, and $\phi$
is the preferred twist angle.

We illustrate the key steps here assuming constant $\kappa_{1}$, $\kappa_{2}$,
and $\phi$. Taking the $s$ derivative of $\boldsymbol{N}$ and separating the
result into tangential and perpendicular parts, we have

$\displaystyle\partial_{s}\boldsymbol{N}$
$\displaystyle=\kappa\left(\boldsymbol{\tau}\times\partial_{s}^{2}{\boldsymbol{\tau}}-\kappa_{1}\partial_{s}\boldsymbol{D}^{(1)}-\kappa\partial_{s}\boldsymbol{D}^{(2)}\right)+\gamma\left(\left(\partial_{s}\psi\right)\boldsymbol{\tau}+(\psi-\phi)\partial_{s}\boldsymbol{\tau}\right)$
(C.2)
$\displaystyle=\boldsymbol{\tau}\times\left(\kappa\left(\partial_{s}^{2}{\boldsymbol{\tau}}+\kappa_{1}\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{D}^{(1)}\right)+\kappa_{2}\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{D}^{(2)}\right)\right)-\gamma\left(\psi-\phi\right)\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)\right)$
(C.3)
$\displaystyle+\left(\kappa\left(-\kappa_{1}\left(\partial_{s}\boldsymbol{D}^{(1)}\cdot\boldsymbol{\tau}\right)-\kappa_{2}\left(\partial_{s}\boldsymbol{D}^{(2)}\cdot\boldsymbol{\tau}\right)\right)+\gamma\partial_{s}\psi\right)\boldsymbol{\tau}$

To derive an Euler model, we now amend (32) and (33) to

$\displaystyle\boldsymbol{F}^{(\kappa)}=-\kappa\left(\partial_{s}^{2}\boldsymbol{\tau}-\kappa_{1}\psi\boldsymbol{D}^{(1)}-\kappa_{2}\psi\boldsymbol{D}^{(2)}\right),\qquad\boldsymbol{F}^{(\gamma)}=\gamma\left(\psi-\phi\right)\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right),$
(C.4)
$\displaystyle\boldsymbol{F}=\boldsymbol{F}^{(\kappa)}+\boldsymbol{F}^{(\gamma)}+T\boldsymbol{\tau}.$

Just like in the case of intrinsically straight and untwisted fibers,
substituting this representation into (27) gives the force and torque density
applied to the fluid as

$\displaystyle\boldsymbol{f}=\partial_{s}\boldsymbol{F},\qquad\boldsymbol{n}=n^{\parallel}\boldsymbol{\tau},\qquad
n^{\parallel}=\kappa\left(\kappa_{1}\left(\boldsymbol{D}^{(1)}\cdot\partial_{s}\boldsymbol{\tau}\right)+\kappa_{2}\left(\boldsymbol{D}^{(2)}\cdot\partial_{s}\boldsymbol{\tau}\right)\right)+\gamma\partial_{s}\psi.$
(C.5)

Free fiber boundary conditions are modified in the case of intrinsic curvature
and twist. As before, we require $\boldsymbol{N}=\boldsymbol{0}$ and
$\boldsymbol{F}=\boldsymbol{0}$ at the free end, which implies that

$\displaystyle\psi=\phi,\qquad\partial_{s}\boldsymbol{\tau}=\kappa_{2}\boldsymbol{D}^{(1)}-\kappa_{1}\boldsymbol{D}^{(2)},\qquad\partial_{s}^{2}{\boldsymbol{\tau}}=\kappa_{1}\phi\boldsymbol{D}^{(1)}+\kappa_{2}\phi\boldsymbol{D}^{(2)},\quad\text{and}\quad
T=0$ (C.6)

at a free end. The equation for $\partial_{s}\boldsymbol{\tau}$ is derived by
using (C.1) and subsituting
$\boldsymbol{D}^{(1)}=\boldsymbol{D}^{(2)}\times\boldsymbol{\tau}$ and
$\boldsymbol{D}^{(2)}=\boldsymbol{\tau}\times\boldsymbol{D}^{(1)}$.

###  C.1 Solving the Bishop ODE

In the case of intrinsic curvature and twist, the material frame is involved
in the forces and torques, and so it must be computed from the Bishop frame
and $\theta$ using (18). Computing the Bishop frame requires the solution of
the Bishop ODE (17), which is subject to the condition (67),

$\displaystyle\partial_{s}\boldsymbol{b}^{(1)}=\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)\times{\boldsymbol{b}^{(1)}}$
(C.7)
$\displaystyle\boldsymbol{b}^{(1)}\left(L/2\right)=\boldsymbol{D}^{(1)}\left(L/2\right).$

Figure 9: The $L^{2}$ function errors in the Bishop vector
$\boldsymbol{b}^{(1)}(s)$. We use the fiber (87) with $q=7$ and generate a
reference solution using $100$ collocation points (these vectors have unit
length to 13 digits). We then solve (C.9) using $N$ points and measure the
$L^{2}$ function error in $\boldsymbol{b}^{(1)}(s)$ on a common grid. We
separate the total error (blue circles) into a component parallel to
$\boldsymbol{b}^{(1)}(s)$ (red squares) and a component perpendicular to
$\boldsymbol{b}^{(1)}(s)$ (yellow diamonds).

To discretize the Bishop ODE (C.7), we first integrate both sides of (C.7) to
obtain an integral equation for $\boldsymbol{b}^{(1)}$. Since each fiber has
at most $40$ points on its centerline, we will use dense linear algebra to
solve for $\boldsymbol{b}^{(1)}$ at all of the collocation points. The
discrete integral form of (C.7) is

$\boldsymbol{b}^{(1)}+\boldsymbol{c}=\boldsymbol{D}^{\dagger}\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)\times\boldsymbol{b}^{(1)}.$
(C.8)

where $\boldsymbol{c}$ is an unknown integration constant and
$\boldsymbol{D}^{\dagger}$ is the pseudo-inverse of the Chebyshev
differentiation matrix. Combining with the boundary condition in (C.7), and
denoting with $\boldsymbol{I}(s)$ evaluation at $s$, we get the system

$\displaystyle\begin{pmatrix}\boldsymbol{I}-\boldsymbol{D}^{\dagger}\left(\left(\boldsymbol{\tau}\times\partial_{s}\boldsymbol{\tau}\right)\times\right)&\boldsymbol{I}\\\
\boldsymbol{I}(L/2)&\boldsymbol{0}\end{pmatrix}\begin{pmatrix}\boldsymbol{b}^{(1)}\\\
\boldsymbol{c}\end{pmatrix}=\begin{pmatrix}\boldsymbol{0}\\\
\boldsymbol{D}^{(1)}(L/2)\end{pmatrix},$ (C.9)

which we solve for $\boldsymbol{b}^{(1)}$.

In Fig. 9, we plot the errors in $\boldsymbol{b}^{(1)}(s)$ as a function of
the grid size for the fiber (87) with $q=7$ that satisfies the free fiber BCs
(the refined solution is generated with $N=100$, in which all vectors have
unit length to tolerance $10^{-13}$). We see that the errors in the parallel
and perpendicular direction are the same, which means that choosing to
renormalize $\boldsymbol{b}^{(1)}$ will not increase the overall error (in
Fig. 9 we do not renormalize $\boldsymbol{b}^{(1)}$).

##  Appendix D Principle of virtual work

In this appendix, we show that the Lagrange multiplier forces
$\boldsymbol{\Lambda}(s)$ in both the Kirchhoff and Euler models do no work on
the fluid. When we substitute (33) and (35) into the corresponding Stokes
equation for the pointwise fluid velocity $\boldsymbol{u}(\boldsymbol{x})$, we
obtain

$\displaystyle\nabla\pi(\boldsymbol{x})=\mu\nabla^{2}\boldsymbol{u}(\boldsymbol{x})+\int_{0}^{L}\left(\partial_{s}\left(\boldsymbol{F}^{(\kappa)}(s)+\boldsymbol{F}^{(\gamma)}\right)(s)\right)\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{X}(s)\right)\,ds$
$\displaystyle\
+\frac{\boldsymbol{\nabla}}{2}\times\int_{0}^{L}n^{\parallel}(s)\boldsymbol{\tau}(s)\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{X}(s)\right)\,ds+\left(\boldsymbol{\mathcal{A}}\boldsymbol{\Lambda}\right)\left(\boldsymbol{x}\right),$
(D.1)

where we have defined a linear operator

$\displaystyle\left(\boldsymbol{\mathcal{A}}\boldsymbol{\Lambda}\right)(\boldsymbol{x})=\int_{0}^{L}\left(\partial_{s}\boldsymbol{\Lambda}(s)\right)\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{X}(s)\right)\,ds+\frac{\boldsymbol{\nabla}}{2}\times\int_{0}^{L}\left(\boldsymbol{\tau}(s)\times\boldsymbol{\Lambda}(s)\right)\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{X}(s)\right)\,ds.$
(D.2)

Starting with the Kirchhoff model, let us consider the constraint (26), which
can be written as

$\partial_{s}\boldsymbol{U}+\boldsymbol{\tau}\times\boldsymbol{\Psi}=\boldsymbol{0}.$
(D.3)

Recalling that $\boldsymbol{U}$ and $\boldsymbol{\Psi}$ are also convolutions
with the regularized function as defined in (5) and (8), we write this
constraint as

$\displaystyle\partial_{s}\left(\int\boldsymbol{u}(\boldsymbol{x})\delta_{\hat{a}}\left(\boldsymbol{X}(s)-\boldsymbol{x}\right)d\boldsymbol{x}\right)+\boldsymbol{\tau}(s)\times\int\left(\frac{\boldsymbol{\nabla}}{2}\times\boldsymbol{u}(\boldsymbol{x})\right)\delta_{\hat{a}}\left(\boldsymbol{X}(s)-\boldsymbol{x}\right)d\boldsymbol{x}=0.$
(D.4)

It is not difficult to show that (D.4) can be written as
$\boldsymbol{\mathcal{A}}^{*}\boldsymbol{u}=\boldsymbol{0}$, where $*$ denotes
an $L^{2}$ adjoint. It follows that

$\displaystyle\left\langle\boldsymbol{\mathcal{A}}\boldsymbol{\Lambda},\boldsymbol{u}\right\rangle=\left\langle\boldsymbol{\Lambda},\boldsymbol{\mathcal{A}}^{*}\boldsymbol{u}\right\rangle=\boldsymbol{0},$
(D.5)

meaning the force $\boldsymbol{\Lambda}$ does no work on the fluid when the
velocity $\boldsymbol{u}$ is constrained to satisfy (26).

The same principle continues to hold in the Euler model as well. Specifically,
since in the Euler model (44) gives $\boldsymbol{\Lambda}=T\boldsymbol{\tau}$,
the forcing term
$\boldsymbol{\mathcal{A}}\boldsymbol{\Lambda}=\boldsymbol{\mathcal{B}}T$,
where

$\boldsymbol{\mathcal{B}}T=\int_{0}^{L}\partial_{s}\left(T(s)\boldsymbol{\tau}(s)\right)\delta_{\hat{a}}\left(\boldsymbol{x}-\boldsymbol{X}(s)\right)\,ds.$
(D.6)

It is straightforward to show that the principle of virtual work
$\langle\boldsymbol{\mathcal{B}}T,\boldsymbol{u}\rangle=0$ is (weakly)
equivalent to the inextensibility constraint (21),

$\boldsymbol{\mathcal{B}}^{*}\boldsymbol{u}=\boldsymbol{\tau}\cdot\partial_{s}\left(\int\boldsymbol{u}(\boldsymbol{x})\delta_{\hat{a}}\left(\boldsymbol{x}-\boldsymbol{X}(s)\right)\,d\boldsymbol{x}\right)=\boldsymbol{\tau}\cdot\partial_{s}\boldsymbol{U}=0.$
(D.7)

Note that both (D.5) and (D.7) assume the free fiber boundary conditions (45).
For clamped fibers, there may be work performed by the forcing at the
endpoints (e.g., the motor spinning the fiber in Section 5.3). The saddle-
point structure of our formulation shows that the dynamical Euler equations we
solve are a generalized gradient descent dynamics on the manifold of
inextensible fibers.

By combining (D.4) and (D.5) we can express the (vanishing) work performed by
the constraint forces in terms of the fiber motion

$\displaystyle\left\langle\boldsymbol{\mathcal{A}}\boldsymbol{\Lambda},\boldsymbol{u}\right\rangle$
$\displaystyle=\int_{0}^{L}\boldsymbol{\Lambda}(s)\cdot\left(\partial_{s}\boldsymbol{U}(s)+\boldsymbol{\tau}(s)\times\boldsymbol{\Psi}(s)\right)\,ds$
(D.8)
$\displaystyle=-\int_{0}^{L}\boldsymbol{U}(s)\cdot\boldsymbol{\lambda}(s)\,ds-\int_{0}^{L}\boldsymbol{n}_{\Lambda}(s)\cdot\boldsymbol{\Psi}(s)\,ds=0,$
(D.9)

where
$\boldsymbol{n}_{\Lambda}=\boldsymbol{\tau}\times\boldsymbol{\Lambda}=\boldsymbol{\tau}\times\boldsymbol{\Lambda}^{\perp}$
is the contribution to the torque density (36) from $\boldsymbol{\Lambda}$. In
the Euler model, $\boldsymbol{\Lambda}^{\perp}=\boldsymbol{0}$ and (D.9)
becomes (58), which we used to enforce $\boldsymbol{\Lambda}(s)$ be
tangential. By contrast, for the Kirchhoff model (D.9) vanished for any
$\boldsymbol{\Psi}$ since (D.3) holds, and therefore there is no additional
constraint on $\boldsymbol{\Lambda}$ from the principle of virtual work.

##  Appendix E Comparison of Kirchhoff and Euler models

In this appendix, we derive and confirm the asymptotic results (55) and (56)
that were given in Section 2.5. After using the estimates (54) to estimate the
difference between the two formulations, we develop a second order method for
the Kirchhoff and Euler models which can robustly handle jumps in the forces
at the endpoints. We then apply the method to a bent, twisted filament to show
the equivalence of the two models (away from the endpoints) as
$\epsilon\rightarrow 0$.

###  E.1 Asymptotic estimates

To estimate the strength of $\Delta\boldsymbol{\Lambda}$, we use the scalings
(54), repeated here for clarity

$\displaystyle\boldsymbol{\mathcal{M}}_{\text{tt}}\sim\log{\epsilon},\qquad\boldsymbol{\mathcal{M}}_{\text{tr}}\sim\log{\epsilon},\qquad\boldsymbol{\mathcal{M}}_{\text{rt}}\sim\log{\epsilon},\qquad\boldsymbol{\mathcal{M}}_{\text{rr}}\sim\epsilon^{-2},$
(E.1)

Beginning with the perpendicular Lagrange multipliers
$\Delta\boldsymbol{\Lambda}^{\perp}=\boldsymbol{\Lambda}^{\perp}$, we first
observe that the mismatch $\boldsymbol{m}$ in (52) is
$\mathcal{O}(\log{\epsilon})$ at most, since we don’t expect
$\boldsymbol{F}_{E}$ as defined in (53) to change with $\epsilon$. So
$\Delta\boldsymbol{\Lambda}^{\perp}$ cannot exceed
$\mathcal{O}(\epsilon^{2}\log{\epsilon})$, since otherwise the angular
velocity from $\boldsymbol{\mathcal{M}}_{\text{rr}}$ would be too large to be
consistent with the $\mathcal{O}(\log{\epsilon})$ mismatch. If
$\Delta\boldsymbol{\Lambda}^{\perp}\sim\epsilon^{2}\log{\epsilon}$, then the
rotational velocity $\Delta\boldsymbol{\Psi}$ is
$\mathcal{O}(\log{\epsilon})$, and therefore of the correct order to cancel
the mismatch $\boldsymbol{m}$, while the translational velocity
$\Delta\boldsymbol{U}=\boldsymbol{\mathcal{M}}_{\text{tt}}\left(\partial_{s}\Delta\boldsymbol{\Lambda}\right)+\boldsymbol{\mathcal{M}}_{\text{tr}}\left(\boldsymbol{\tau}\times\Delta\boldsymbol{\Lambda}\right)$,
is of the order
$\Delta\boldsymbol{U}\sim\epsilon^{2}\left(\log{\epsilon}\right)^{2}$ (if we
assume that $\Delta\boldsymbol{\Lambda}$ and its derivative are of the same
order, which should be true away from the endpoints).777Recall that
$\Delta\boldsymbol{\Lambda}^{\parallel}$ is a Lagrange multiplier that makes
the perturbation velocity $\Delta\boldsymbol{U}$ inextensible. Since the
interior velocity is at worst
$\Delta\boldsymbol{U}\sim\epsilon^{2}\left(\log{\epsilon}\right)^{2}$ and the
translational mobility is $\mathcal{O}(\log{\epsilon})$, the Lagrange
multiplier
$\Delta\boldsymbol{\Lambda}^{\parallel}\sim\epsilon^{2}\log{\epsilon}$ at
worst, just like $\Delta\boldsymbol{\Lambda}^{\perp}$.  This implies that in
the fiber interior
$\partial_{s}\Delta\boldsymbol{U}\sim\epsilon^{2}\left(\log{\epsilon}\right)^{2}$,
which gives (55).

The difference in parallel angular velocity $\Delta\Omega^{\parallel}$ is
given asymptotically by

$\displaystyle\Delta\Omega^{\parallel}=\boldsymbol{\tau}\cdot\Delta\boldsymbol{\Psi}$
$\displaystyle=\boldsymbol{\tau}\cdot\left(\boldsymbol{\mathcal{M}}_{\text{rt}}\left(\partial_{s}\Delta\boldsymbol{\Lambda}\right)+\boldsymbol{\mathcal{M}}_{\text{rr}}\left(\boldsymbol{\tau}\times\Delta\boldsymbol{\Lambda}\right)\right)\approx\boldsymbol{\tau}\cdot\left(\boldsymbol{\mathcal{M}}_{\text{rt}}\left(\partial_{s}\Delta\boldsymbol{\Lambda}\right)\right),$
(E.2)

where the last approximate equality is because
$\boldsymbol{\mathcal{M}}_{\text{rr}}\left(\boldsymbol{\tau}\times\Delta\boldsymbol{\Lambda}\right)$
is orthogonal to $\boldsymbol{\tau}$ to leading order (see Section 2.1.3). In
fact, in the fiber interior
$\boldsymbol{\mathcal{M}}_{\text{rr}}\left(\boldsymbol{\tau}\times\Delta\boldsymbol{\Lambda}\right)\cdot\boldsymbol{\tau}\sim\log{\epsilon}\left\|\Delta\boldsymbol{\Lambda}\right\|\sim\epsilon^{2}\left(\log{\epsilon}\right)^{2}$.
Adding the rot-trans contribution, since
$\Delta\boldsymbol{\Lambda}\sim\epsilon^{2}\log{\epsilon}$ and
$\boldsymbol{\mathcal{M}}_{\text{rt}}\sim\log{\epsilon}$, we get
$\Delta\Omega^{\parallel}\sim\epsilon^{2}\left(\log{\epsilon}\right)^{2}$. To
non-dimensionalize this, we need to estimate the parallel rotational velocity
$\Omega^{\parallel}_{E}$ defined in (47). When there is no twist, there is no
parallel torque, so
$\Omega^{\parallel}_{E}=\boldsymbol{\tau}\cdot\left(\boldsymbol{\mathcal{M}}_{\text{rt}}\partial_{s}\boldsymbol{F}_{E}\right)\sim\log{\epsilon}$,
and the relative difference should be
$\left\|\Delta\Omega^{\parallel}\right\|/\left\|\Omega_{E}^{\parallel}\right\|\sim\epsilon^{2}\log{\epsilon}$.
When there is twist, the disparity is even larger, since
$\Omega^{\parallel}_{E}\sim\epsilon^{-2}$, which explains the two different
cases in (56).

Summarizing, the mismatch $\boldsymbol{m}$ in (51) is compensated for by
$\Delta\boldsymbol{\Psi}^{\perp}$, which grows as $\epsilon\rightarrow 0$,
even as $\Delta\boldsymbol{\Lambda}$ decays. These are the components of the
velocity that are discarded in the Euler model, so it does not matter if they
differ from those of the Kirchhoff model. That said, because
$\Delta\boldsymbol{\Lambda}$ must be zero at the endpoints by the boundary
conditions (31), the perpendicular torques cannot generate angular velocity
there, and the term $\Delta\boldsymbol{\Psi}^{\perp}$ cannot compensate for
the mismatch $\boldsymbol{m}$ at the endpoints. Instead, the mismatch there
has to be compensated for by terms involving
$\partial_{s}\Delta\boldsymbol{\Lambda}$, meaning that the derivative of
$\Delta\boldsymbol{\Lambda}$ must be large near the boundaries, as we see in
Fig 10. This makes a spectral method a poor choice for the Kirchhoff rod
formulation, and is the reason why we develop a second-order method to solve
the Kirchhoff equations next.

###  E.2 Second-order numerical method

Following the method outlined in [62], we introduce a number of links $N$ and
set $\Delta s=L/N$. The discretization of the fiber is _staggered_ : the
“nodes” are located at $s_{j}=j\Delta s$, where $j=0,\dots N$, so that there
are actually $N+1$ nodes. The “links” are centered on $s_{j+1/2}=(j+1/2)\Delta
s$, where $j=0,\dots,N-1$, so that there are $N$ links.

####  E.2.1 Euler model

In our second-order discretization of the Euler problem (62), we define
$\boldsymbol{\lambda}$, $\boldsymbol{f}^{(\kappa)}$, and
$\boldsymbol{f}^{(\gamma)}$ at link centers $(j+1/2)\Delta s$, for $j=0,\dots
N-1$. All mobility operators in (7) and (10) (including the rot-rot mobility)
are then replaced with the obvious discretizations of the integrals which sum
over the links with weight $\Delta s$. For example,

$\left(\boldsymbol{\mathcal{M}}_{\text{tt}}\boldsymbol{f}\right)(s_{i+1/2})\approx\sum_{j=0}^{N-1}\mathbb{M}_{\text{tt}}\left(\boldsymbol{X}_{i+1/2},\boldsymbol{X}_{j+1/2}\right)\boldsymbol{f}_{j+1/2}\Delta
s:=\left(\boldsymbol{M}_{\text{tt}}\boldsymbol{f}\right)_{i+1/2}.$ (E.3)

Note that this discretization does not correspond to a more common
discretization of an inextensible bead-link polymer chain, which would place
forces on $N+1$ blobs connected by $N$ inextensible links, e.g., [52]. We use
the staggered discretization instead for a more direct comparison to the
Kirchhoff model.

To complete the second-order discretization of (62), we have to define the
matrices $\boldsymbol{K}$ and $\boldsymbol{K}^{*}$. In the second-order case,
these matrices are exact transposes of each other, so we only have to define
one of them. The matrix $\boldsymbol{K}$ is a midpoint discretization of the
continuum operator (57), with action given by

$\displaystyle\left(\boldsymbol{K}\boldsymbol{\alpha}\right)_{j+1/2}=\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\boldsymbol{U}$\kern-1.00006pt}}}+\sum_{k=0}^{j-1}$
$\displaystyle\Delta
s\left(\alpha_{1,k+1/2}\boldsymbol{n}_{1,k+1/2}+\alpha_{2,k+1/2}\boldsymbol{n}_{2,k+1/2}\right)$
(E.4) $\displaystyle+$ $\displaystyle\frac{\Delta
s}{2}\left(\alpha_{1,j+1/2}\boldsymbol{n}_{1,j+1/2}+\alpha_{2,j+1/2}\boldsymbol{n}_{2,j+1/2}\right).$

The midpoint quadrature assigns a weight of $\Delta s$ to all prior nodes
$1/2,3/2,\dots j-1/2$, and then a weight of $\Delta s/2$ to the current node
$j+1/2$. The matrix $\boldsymbol{K}^{*}$ is the transpose of $\boldsymbol{K}$.

Finally, the parallel angular velocity (64) is obtained at the link centers by
dotting the full angular velocity,
$\boldsymbol{\mathcal{M}}_{\text{rt}}\boldsymbol{f}+\boldsymbol{\mathcal{M}}_{\text{rr}}\left(n^{\parallel}\boldsymbol{\tau}\right)$
(computed at the link centers using the forces computing at the link centers),
with the tangent vectors $\boldsymbol{\tau}$, which are also defined at the
link centers.

####  E.2.2 Kirchhoff model

The discretization of the Kirchhoff model, which we take from [62], is more
complicated because the forces $\boldsymbol{\Lambda}$ are defined at the
_nodes_ , and not the link centers. To communicate between the nodes and
links, we let $\boldsymbol{D}_{\text{lb}}$ be the $N\times(N+1)$
differentiation matrix, which gives derivates at the links from values at the
nodes. We also define $\boldsymbol{D}_{\text{bl}}$ as the $(N-1)\times N$
differentiation matrix, which gives derivates at interior nodes from values at
the link centers. Likewise the $N\times(N+1)$ averaging matrix
$\boldsymbol{A}_{\text{lb}}$ gives averages at links from nodes, and the
$(N-1)\times N$ averaging matrix $\boldsymbol{A}_{\text{bl}}$ gives averages
at interior nodes from links. This gives the $(N-1)\times(N+1)$ system of
equations discretizing (51)

$\displaystyle\left[\boldsymbol{D}_{\text{bl}}\left(\boldsymbol{M}_{\text{tt}}\boldsymbol{D}_{\text{lb}}+\boldsymbol{M}_{\text{tr}}\boldsymbol{C}\boldsymbol{A}_{\text{lb}}\right)+\boldsymbol{A}_{\text{bl}}\boldsymbol{C}\left(\boldsymbol{M}_{\text{rt}}\boldsymbol{D}_{\text{lb}}+\boldsymbol{M}_{\text{rr}}\boldsymbol{C}\boldsymbol{A}_{\text{lb}}\right)\right]\Delta\boldsymbol{\Lambda}=-\boldsymbol{m}$
(E.5)
$\displaystyle=-\left(\boldsymbol{D}_{\text{bl}}\boldsymbol{M}_{\text{tt}}+\boldsymbol{A}_{\text{bl}}\boldsymbol{C}\boldsymbol{M}_{\text{rt}}\right)\left(\boldsymbol{f}^{(\kappa)}+\boldsymbol{f}^{(\gamma)}+\boldsymbol{\lambda}_{E}\right)-\left(\boldsymbol{D}_{\text{bl}}\boldsymbol{M}_{\text{tr}}+\boldsymbol{A}_{\text{bl}}\boldsymbol{C}\boldsymbol{M}_{\text{rr}}\right)\left(n^{\parallel}\boldsymbol{\tau}\right)$
(E.6)

which we solve for $\Delta\boldsymbol{\Lambda}$ at the $N-1$ _interior blobs_
away from the fiber boundaries. Because $\boldsymbol{\Lambda}$ is defined at
$N+1$ points, the other two equations are the BCs
$\Delta\boldsymbol{\Lambda}_{j}=\boldsymbol{0}$ for $j=0,N$, see (31). The
matrix $\boldsymbol{C}$ represents cross products with
$\boldsymbol{\tau}_{j+1/2}=\boldsymbol{\tau}\left((j+1/2)\Delta s\right)$,
where in the second-order method
$\boldsymbol{\tau}=\partial_{s}\boldsymbol{X}$ is evaluated analytically at
the links and
$\left(\boldsymbol{C}\boldsymbol{X}\right)_{j+1/2}=\boldsymbol{\tau}_{j+1/2}\times\boldsymbol{X}_{j+1/2}$.
Notice that the calculation of the mismatch (E.6) requires calculation of the
Euler quantities $\boldsymbol{f}^{(\kappa)}$ and $\boldsymbol{f}^{(\gamma)}$
and $\boldsymbol{\lambda}_{E}$ at the $N$ links, which we first obtain by
solving the Euler problem to second order accuracy using the method in Section
E.2.1.

###  E.3 Numerical results

To numerically compare the Euler and Kirchhoff constitutive models, we will
apply the second-order methods to the fiber (87) with $q=1$ and $L=2$ (see
Fig. 1 for a visual). This fiber is chosen to satisfy the free fiber boundary
conditions (34). We will also assume an instantaneous twist

$\psi=\sin{(2\pi s)},$ (E.7)

so that all of the free fiber boundary conditions (34) are satisfied. Because
$\boldsymbol{X}$ satisfies the free-fiber BCs, we can compute the elastic
force to spectral accuracy simply by differentiating $\boldsymbol{\tau}$ three
times analytically, and then setting
$\boldsymbol{f}^{(\kappa)}=-\kappa\partial_{s}^{3}{\boldsymbol{\tau}}$. For
this section, we will compute all geometric quantities, for example
$\partial_{s}\boldsymbol{\tau}$, $\partial_{s}^{2}{\boldsymbol{\tau}}$ and
$\partial_{s}\psi$, analytically. We will use a viscosity of $\mu=1$ and
aspect ratios $\hat{\epsilon}={\hat{a}}/L$ that are a function of the
regularized singularity radius ${\hat{a}}$; the corresponding value of the
true fiber radius $a\approx{\hat{a}}/1.12$ (see (1)).

####  E.3.1 Bending but no twisting

Figure 10:  Solving the mismatch problem (51) with bending stiffness
$\kappa=1$, twist modulus $\gamma=0$, and slenderness
$\hat{\epsilon}={\hat{a}}/L=0.005$. All quantities are shown in the basis of
the local triad, with the $\boldsymbol{\tau}$ component in blue and the
$\boldsymbol{D}^{(1/2)}$ component in red/yellow. In the top row, from left to
right, we show $\partial_{s}\boldsymbol{U}_{E}$, (negative) mismatch
$-\boldsymbol{m}$, and the corresponding $\Delta\boldsymbol{\Lambda}$. Notice
how the shapes of $\Delta\boldsymbol{\Lambda}$ and $\boldsymbol{m}$ match,
which implies that $\Delta\boldsymbol{\Lambda}(s)\approx
c(\hat{\epsilon})\boldsymbol{m}(s)$ away from the fiber boundaries. In the
bottom row, from left to right, we show the resulting velocity perturbations
$\Delta\Omega^{\parallel}$, $\partial_{s}\Delta\boldsymbol{U}$, and
$\boldsymbol{\tau}\times\Delta\boldsymbol{\Psi}$. In the fiber interior,
$\boldsymbol{\tau}\times\Delta\boldsymbol{\Psi}\approx-\boldsymbol{m}$, and
$\partial_{s}\Delta\boldsymbol{U}(s)\approx\boldsymbol{0}$. At the endpoints,
$\partial_{s}\Delta\boldsymbol{U}$ reaches values as negative as $-10$, so we
truncate the $y$ axis to show the behavior in the fiber interior.

We first solve the mismatch problem (51) with $\kappa=1$ and $\gamma=0$ using
our second-order method with $N=1/\hat{\epsilon}$ points. In Fig. 10, we plot
$\boldsymbol{m}$, $\Delta\boldsymbol{\Lambda}$, $\Delta\boldsymbol{U}$,
$\boldsymbol{\tau}\times\Delta\boldsymbol{\Psi}$, and
$\partial_{s}\Delta\boldsymbol{U}$ in the basis of the local triad for
$\hat{\epsilon}={\hat{a}}/L=0.005$. We observe several features: first, the
mismatch $\boldsymbol{m}$ is smooth until it reaches the endpoints, after
which it jumps to a large value at the boundaries. In the fiber interior, the
function $\Delta\boldsymbol{\Lambda}$ is approximately a constant multiple of
$\boldsymbol{m}$, but at the endpoints it transitions to zero (as required by
the boundary condition). Away from the fiber boundaries, the mismatch
$\boldsymbol{m}$ is almost entirely compensated for by the angular velocity
perturbation $\boldsymbol{\tau}\times\Delta\boldsymbol{\Psi}$, and
$\partial_{s}\Delta\boldsymbol{U}$ is small in the fiber interior, validating
(55).

Figure 10 also shows that $\boldsymbol{\tau}\times\Delta\boldsymbol{\Psi}$ by
itself does not compensate for the mismatch $\boldsymbol{m}$ near the
boundaries, since the dominant contribution from torque
$\Delta\boldsymbol{\Lambda}$ must be zero there. The mismatch $\boldsymbol{m}$
near the boundaries must therefore be compensated for by
$\partial_{s}\Delta\boldsymbol{U}$, which means that $\Delta\boldsymbol{U}$ is
highly nonsmooth at the boundaries. Since $\Delta\boldsymbol{U}$ comes from
$\Delta\boldsymbol{\Lambda}$, it is not surprising that
$\Delta\boldsymbol{\Lambda}$ is also not smooth at the boundaries, as it takes
on its largest values near the endpoints, only to have to fall back to zero as
required by the boundary conditions.

(a) $\kappa=1,\gamma=0$

(b) $\kappa=0,\gamma=1$

Figure 11: Decay of (properly normalized) magnitude of
$\Delta\boldsymbol{\Lambda}$ (yellow), $\partial_{s}\Delta\boldsymbol{U}$
(blue), and $\Delta\Omega^{\parallel}$ (red) as slenderness
$\hat{\epsilon}\rightarrow 0$ for: (a) $\gamma=0$ and $\kappa=1$ with
$N=1/\hat{\epsilon}+1$ grid points; and, (b) $\gamma=1$ and $\kappa=0$ and
$N=4/\hat{\epsilon}+1$ grid points. We use an interior “norm” ($L^{2}$ norm on
$0.05\leq s/L\leq 0.95$). All quantities scale like $\hat{\epsilon}^{2}$ when
properly normalized, with the exception of $\Delta\Omega^{\parallel}$ when
there is twist, which scales faster than $\hat{\epsilon}^{3}$. This
establishes (55) and (56) numerically.

We will focus on the equivalence of the Euler and Kirchhoff models in the
fiber _interior_ , or points $\mathcal{O}(1)$ from the fiber endpoints. We
therefore define an interior “norm” for all functions which is simply the
$L^{2}$ norm on $0.05\leq s/L\leq 0.95$. In Fig. 11(a), we plot the relative
difference in the velocity and forces using this interior norm. We observe
convergence of order $\hat{\epsilon}^{2}$ for all quantities. This verifies,
and in fact exceeds, the predictions of (55) and (56) of
$\hat{\epsilon}^{2}\log{\hat{\epsilon}}$. The reason is that the mismatch is
actually dominated by $\mathcal{O}(1)$ terms for the values of
$\hat{\epsilon}$ we are interested in, and so our estimates are off by a
factor of $\log{\hat{\epsilon}}$.

####  E.3.2 Twisting but no bending

Figure 12:  Solutions of the mismatch problem (51) with twist modulus
$\gamma=1$, bending modulus $\kappa=0$, and slenderness ratio
$\hat{\epsilon}={\hat{a}}/L=0.005$. The behavior is qualitatively the same as
the mismatch problem without twist, with the perturbations in
$\Delta\Omega^{\parallel}$ and $\partial_{s}\Delta\boldsymbol{U}$ being small
relative to the values obtained from the Euler method. Because the twist force
$\boldsymbol{f}^{(\gamma)}$ is zero at the endpoints, the mismatch there is
small, and $\partial_{s}\Delta\boldsymbol{U}$ is smoother than seen in Fig.
10.

We now proceed to the case of finite twist modulus, $\gamma=1$, with zero bend
modulus. We first plot the quantities related to the mismatch problem (51) in
Fig. 12. We see similar behavior to the case without twist and with bending
force. The velocity difference $\partial_{s}\Delta\boldsymbol{U}$ is again
small, since most of the mismatch is compensated for by
$\boldsymbol{\tau}\times\Delta\boldsymbol{\Psi}$. A major difference from the
non-twisted case is that $\Delta\boldsymbol{\Lambda}$ and
$\Delta\boldsymbol{U}$ are smoother at the boundaries, since the force when
$\kappa=0$, given in (38), is zero at the fiber boundaries because of the free
fiber boundary conditions (34). Because the force decays to zero there, the
mismatch is much smaller at the boundaries (but is still not zero). This means
that $\boldsymbol{\tau}\times\Delta\boldsymbol{\Psi}$ can almost compensate
for the entire mismatch everywhere. At the boundaries, whatever small mismatch
there is must be compensated for by $\partial_{s}\Delta\boldsymbol{U}$, and so
as $\hat{\epsilon}$ becomes smaller we do see some localized peaks in
$\partial_{s}\Delta\boldsymbol{U}$ near the endpoints.

In Figure 11(b), we plot the convergence of the velocities and forces for
$\gamma=1$ and $\kappa=0$. For $\Delta\boldsymbol{\Lambda}$ and
$\Delta\boldsymbol{U}$, we observe a similar rate of convergence
(approximately $\hat{\epsilon}^{2}$) to Fig. 11(a). In accordance with (56),
the parallel angular velocity $\Omega^{\parallel}_{E}$ is larger when there is
twist, and so as a result the relative rate at which
$\Delta\Omega^{\parallel}$ approaches zero is faster. In (56) we estimated
$\hat{\epsilon}^{4}\left(\log{\hat{\epsilon}}\right)^{2}$; in Fig. 11(b) we
simply show it is faster than $\hat{\epsilon}^{3}$, but the line isn’t
straight, which indicates some kind of log correction. This verifies (55) and
(56) numerically for finite twist modulus.

##  Appendix F Discretization of kinematic operators

In this appendix, we discuss the spectral discretization of the kinematic
operators (57) and (60) on a type 1 Chebyshev grid of size $N$. For simplicity
of notation, we will consider $s=0$ to be a clamped end so that
$\boldsymbol{U}_{0}=\boldsymbol{0}$ and
$\boldsymbol{\alpha}\equiv\boldsymbol{\Omega}$. The continuum operators
$\boldsymbol{\mathcal{K}}$ and $\boldsymbol{\mathcal{K}}^{*}$ (c.f. (57) and
(60)) in this case are defined as

$\displaystyle\boldsymbol{U}(s)=\left(\boldsymbol{\mathcal{K}}\boldsymbol{\alpha}\right)(s)=\int_{0}^{s}\boldsymbol{\Omega}(s^{\prime})\times\boldsymbol{\tau}(s^{\prime})\,ds^{\prime},$
(F.1)
$\displaystyle\boldsymbol{\mathcal{K}}^{*}\boldsymbol{\lambda}:=\boldsymbol{\tau}(s)\times\int_{s}^{L}\boldsymbol{\lambda}(s^{\prime})\,ds^{\prime}=\boldsymbol{0}.$
(F.2)

These operators must be discretized with care because of nonlinearities, which
lead to aliasing on the spectral collocation grid. Our approach in this paper,
which is somewhat different from what we have done previously [48, Sec.
4.1.2], is to discretize these operators on a type 1 Chebyshev grid of size
$2N$, and then appropriately upsample and downsample to a grid of size $N$.
This is in contrast to [48], where we set $\Omega^{\parallel}(s)=0$ on the $N$
point grid.

As in [48, Sec. 4.1.2], we let $\boldsymbol{J}$ be a matrix such that
$\boldsymbol{J}\boldsymbol{\Omega}$ is an inextensible motion
$\int_{0}^{s}\boldsymbol{\Omega}(s^{\prime})\times\boldsymbol{\tau}(s^{\prime})ds^{\prime}$
on a grid of size $2N$. Unlike [48, Sec. 4.1.2], however, we will use
$\boldsymbol{J}$ to project off extensible motions on the $2N$ grid, rather
than trying to write an explicit basis for inextensible motions on the $N$
point grid. To define $\boldsymbol{J}$, we let
$\boldsymbol{C}(\boldsymbol{\tau})$ be the matrix that computes cross products
with $\boldsymbol{\tau}$ on the $2N$ grid, i.e.,
$\boldsymbol{C}\boldsymbol{x}=\boldsymbol{\tau}\times\boldsymbol{x}$, with
$\boldsymbol{C}^{T}=-\boldsymbol{C}$. We use the pseudoinverse of the
Chebyshev differentiation matrix [68, c. 6] for integration, with
$\boldsymbol{D}^{\dagger}_{0}$ denoting integration from 0 to $s$, which has
$L^{2}$ adjoint
$\left(\boldsymbol{D}^{\dagger}_{0}\right)^{*}=-\boldsymbol{D}^{\dagger}_{L}$
(integration from $s$ to $L$). The discretization of $\boldsymbol{J}$ and its
$L^{2}$ adjoint $\boldsymbol{J}^{*}$ on a grid of size $2N$ can then be
written as

$\displaystyle\boldsymbol{J}=\boldsymbol{D}_{0}^{\dagger}\boldsymbol{C},\qquad\qquad\boldsymbol{J}^{*}=\boldsymbol{C}\boldsymbol{D}_{L}^{\dagger}.$
(F.3)

Here $\boldsymbol{J}$ is a $6N\times 6N$ matrix that discretizes (F.1), and
$\boldsymbol{J}^{*}$ is a $6N\times 6N$ matrix that discretizes (F.2).

To write a discretization on the $N$ point grid, we introduce an extension
matrix $\boldsymbol{E}$ and restriction matrix $\boldsymbol{R}$. The extension
matrix is well defined: to evaluate $\boldsymbol{E}f_{N}$, where $f_{N}$ is a
polynomial of degree $N-1$ on the $N$ point grid, we simply evaluate $f$ on
the grid of size $2N$. Following [48, Eq. (62)], we let the restriction matrix
$\boldsymbol{R}$ be the matrix that gives the best possible (in the $L^{2}$
sense) polynomial approximation of degree $N-1$ on the $N$ grid to a
polynomial of degree $2N-1$ on the $2N$ point grid, i.e., if $f_{2N}$ is a
function on the $2N$ point grid, we want to choose $\boldsymbol{R}$ such that
$\left\|\boldsymbol{E}\boldsymbol{R}f_{2N}-f_{2N}\right\|_{L^{2}}$ is
minimized. This can be accomplished by setting

$\boldsymbol{R}=\left(\boldsymbol{E}^{T}\boldsymbol{W}_{2N}\boldsymbol{E}\right)^{-1}\boldsymbol{E}^{T}\boldsymbol{W}_{2N},$
(F.4)

where $\boldsymbol{W}_{2N}$ is a diagonal matrix of Clenshaw-Curtis weights on
the $2N$ point grid. We do this for both $\boldsymbol{K}$ and
$\boldsymbol{K}^{*}$ to obtain the discretizations of
$\boldsymbol{\mathcal{K}}$ and $\boldsymbol{\mathcal{K}}^{*}$,

$\boldsymbol{K}=\boldsymbol{R}\boldsymbol{J}\boldsymbol{E},\qquad\text{and}\qquad\boldsymbol{K}^{*}=\boldsymbol{R}\boldsymbol{J}^{*}\boldsymbol{E}$
(F.5)

where $\boldsymbol{K}$ and $\boldsymbol{K}^{*}$ are both $3N\times 3N$
matrices. It is straightforward to modify $\boldsymbol{K}$ and
$\boldsymbol{K}^{*}$ to allow a nonzero $\boldsymbol{U}_{0}$ for a free end.

To update the fiber tangent vectors from the velocity
$\boldsymbol{U}=\boldsymbol{K}\boldsymbol{\alpha}$, we follow (63). This
update to the tangent vectors projects off any parallel components in
$\boldsymbol{U}_{s}$, which come from numerical discretization artifacts.
Therefore, the parallel component of $\boldsymbol{\Omega}$ on the $N$ point
grid does not affect the dynamics; it is included merely to lessen aliasing
artifacts.

As discussed in Section 4.1, the main difference between our discretization
here and that of [48, Sec. 4.1.2] is the treatment of the null space of the
continuum operator $\boldsymbol{\mathcal{K}}$. Because we project off the null
space of $\boldsymbol{\mathcal{K}}$ on a grid of size $2N$ and downsample the
result, the rank of the matrix $\boldsymbol{K}$ (and $\boldsymbol{K}^{*}$) is
a function of the fiber curvature. For a rigid fiber, the rank of
$\boldsymbol{K}$ is $2N+1$, which is the same as when we directly write a
basis for inextensible motions on an $N$ point grid. For curved fibers, the
rank is larger because inextensible motions are aliased to extensible ones on
the $N$ point grid. For the fiber (87) with $q=7$ for instance, for large $N$
we observe $2N+23$ nonzero singular values of $\boldsymbol{K}$, so there are
22 additional positive singular values on top of the $2N+1$ nonzero singular
values for a rigid fiber. Hence, typically the matrix $\boldsymbol{K}$ does
not have full rank, so the Schur complement
$\boldsymbol{K}^{*}\boldsymbol{M}^{-1}\boldsymbol{K}$ that we use to solve the
saddle point system (62) is singular. We therefore use the pseudo-inverse of
$\boldsymbol{K}^{*}\boldsymbol{M}^{-1}\boldsymbol{K}$ to solve the saddle
point system (pinv function in Matlab with default tolerance).

##  Appendix G Nearly singular quadrature

In this appendix, we discuss the nearly singular quadrature scheme we use to
evaluate the integrals of the RPY kernels in Section 4.2 to spectral accuracy.
The idea behind this scheme has already been published in [1, 66, 48], but in
those works the nearly-singular integral is actually a singular finite part
integral because the domain is $[0,L]$ instead of $D(s)$ as defined in (69).
This section presents a review of the quadrature scheme and how we modify it
when the domain changes. We do this sequentially for the Stokeslet integral,
doublet integral, and rotlet integral.

###  G.1 Stokeslet integral

For the Stokeslet integral (70), the nearly-singular integral that we must
evaluate is the Stokeslet minus the leading order singularity,

$\boldsymbol{U}_{\text{f}}^{\text{(int,
S)}}(s)=\int_{D(s)}\left(\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}\left(s^{\prime}\right)-\frac{1}{8\pi\mu}\left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s-s^{\prime}|}\right)\boldsymbol{f}(s)\right)\,ds^{\prime},$
(G.1)

which has the same integrand as the finite part integral in slender body
theory, but with the different domain of integration

$D(s)=\begin{cases}\left(0,s-2{\hat{a}}\right)\cup\left(s+2{\hat{a}},L\right)&2{\hat{a}}\leq
s\leq L-2{\hat{a}}\\\ \left(s+2{\hat{a}},L\right)&s<2{\hat{a}}\\\
\left(0,s-2{\hat{a}}\right)&s>L-2{\hat{a}}\end{cases}.$ (G.2)

We employ a singularity subtraction scheme that ensures that the second term
in the integral (G.1) cancels the leading order $1/|s^{\prime}-s|$ singularity
in the first. The next singularity is $\text{sign}(s^{\prime}-s)$, which means
that the near singular integral (G.1) can be written as

$\displaystyle\boldsymbol{U}_{\text{f}}^{\text{(int,
S)}}(s)=\int_{D(s)}\boldsymbol{g}_{\text{tt}}(s,s^{\prime})\frac{s^{\prime}-s}{|s^{\prime}-s|}\,ds^{\prime}=\frac{L}{2}\int_{D^{\prime}(\eta)}\boldsymbol{g}_{\text{tt}}(\eta,\eta^{\prime})\frac{\eta^{\prime}-\eta}{|\eta^{\prime}-\eta|}\,d\eta^{\prime},\qquad\text{where}$
(G.3)
$\displaystyle\boldsymbol{g}_{\text{tt}}(s,s^{\prime})=\left(\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}\left(s^{\prime}\right)|s^{\prime}-s|-\frac{1}{8\pi\mu}\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)\right)\frac{1}{s^{\prime}-s}$

is a smooth function, $\eta=-1+2s/L$ is a rescaled arclength coordinate on
$[-1,1]$, and $D^{\prime}(\eta)=[0,\eta_{\ell}]\cup[\eta_{h},L]$ is defined
from $D(s)$. The function $\boldsymbol{g}_{\text{tt}}$ is nonsingular at
$s=s^{\prime}$ with the finite limit

$\lim_{s^{\prime}\to
s}\boldsymbol{g}_{\text{tt}}(s,s^{\prime})=\frac{1}{8\pi\mu}\left(\frac{1}{2}\left(\boldsymbol{\tau}(s)\partial_{s}\boldsymbol{\tau}(s)+\partial_{s}\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)+\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\partial_{s}\boldsymbol{f}(s)\right).$
(G.4)

Furthermore, $\boldsymbol{g}_{\text{tt}}$ is smooth, so we can express it in a
truncated Chebyshev series on $[-1,1]$,

$\frac{L}{2}\boldsymbol{g}_{\text{tt}}(\eta,\eta^{\prime})\approx\sum_{k=0}^{N-1}\boldsymbol{c}_{k}(\eta)T_{k}(\eta^{\prime}),$
(G.5)

where $\boldsymbol{c}_{k}$ is a vector of 3 coefficients for each $\eta$.
Substituting the Chebyshev expansion (G.5) into the integrand (G.3), we obtain

$\displaystyle\boldsymbol{U}_{\text{f}}^{\text{(int, S)}}(\eta)$
$\displaystyle=\sum_{k=0}^{N-1}\boldsymbol{c}_{k}(\eta)\int_{D^{\prime}(\eta)}T_{k}\left(\eta^{\prime}\right)\frac{\eta^{\prime}-\eta}{|\eta^{\prime}-\eta|}\,d\eta^{\prime}=\sum_{k=0}^{N-1}\boldsymbol{c}_{k}(\eta)q^{(S)}_{k}(\eta)=\boldsymbol{c}^{T}(\eta)\boldsymbol{q}^{(S)}(\eta),$
(G.6) $\displaystyle\text{where}\qquad q_{k}^{(S)}(\eta)$
$\displaystyle=\int_{D^{\prime}(\eta)}T_{k}(\eta^{\prime})\frac{\eta^{\prime}-\eta}{|\eta^{\prime}-\eta|}\,d\eta^{\prime}=-\int_{-1}^{\eta_{\ell}}T_{k}(\eta^{\prime})\,d\eta^{\prime}+\int_{\eta_{h}}^{L}T_{k}(\eta^{\prime})\,d\eta^{\prime},$
(G.7)

are integrals that can be precomputed to high accuracy for each $\eta$ on the
Chebyshev collocation grid. The precomputed integrals in (G.7) can then be
used in an adjoint method to accelerate the repeated evaluation of (G.6).
Specifically, the coefficients of $\boldsymbol{g}_{\text{tt}}(\eta)$ for each
$\eta$ can be written as
$\boldsymbol{c}(\eta)=\boldsymbol{V}^{-1}\boldsymbol{g}_{\text{tt}}(\eta,\eta^{\prime})$,
where $\boldsymbol{V}$ is the matrix that maps coefficients of the basis
functions to values on the collocation grid. Then
$\boldsymbol{c}^{T}(\eta)=\boldsymbol{g}^{T}(\eta,\eta^{\prime})\boldsymbol{V}^{-T}$,
and $\boldsymbol{U}_{\text{f}}^{\text{(int,
S)}}(\eta)=\boldsymbol{g}^{T}(\eta,\eta^{\prime})\boldsymbol{V}^{-T}\boldsymbol{q}^{(S)}(\eta)$.
The matrix vector product $\boldsymbol{V}^{-T}\boldsymbol{q}^{(S)}(\eta)$ can
then be precomputed for each $\eta$, and so at each time step all that is
needed to compute the integral (G.1) for each $\eta$ are the values of
$\boldsymbol{g}_{\text{tt}}(\eta,\eta^{\prime})$ on the collocation grid for
$\eta^{\prime}$ being another collocation point. Evaluating
$\boldsymbol{U}_{\text{f}}^{\text{(int, S)}}$ therefore has a cost of
$\mathcal{O}\left(N^{2}\right)$.

We note that our use of Chebyshev polynomials differs from the approach of
[66, 48], which uses a monomial basis for the expansion (G.5). Although the
integrals (G.7) with monomials can be computed analytically, using the adjoint
method with monomial coefficients requires inverting the Vandermonde matrix
$\boldsymbol{V}$, which is too ill-conditioned with $N\gtrsim 40$ points in
double precision. The Chebyshev expansion is an improvement in the sense that
computing coefficients from values on the Chebyshev grid is a well-conditioned
problem, and so the technique we use to evaluate finite part integrals does
not limit the number of points on the fiber or require the use of panels
(composite quadrature).

###  G.2 Doublet integral

The next part of the translational mobility (68) is the integral of the
doublet kernel. We begin by applying the same singularity subtraction
technique we used in (70) for the Stokeslet integral

$\displaystyle\boldsymbol{U}_{\text{f}}^{(D)}(s)=$
$\displaystyle\int_{D(s)}\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}\left(s^{\prime}\right)\,ds^{\prime}$
(G.8) $\displaystyle=\frac{1}{8\pi\mu}$
$\displaystyle\int_{D(s)}\left(\frac{\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s^{\prime}-s|^{3}}\right)\boldsymbol{f}(s)\,ds^{\prime}$
$\displaystyle+$
$\displaystyle\int_{D(s)}\left(\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}\left(s^{\prime}\right)-\frac{1}{8\pi\mu}\left(\frac{\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s^{\prime}-s|^{3}}\right)\boldsymbol{f}(s)\right)\,ds^{\prime}$
$\displaystyle=$ $\displaystyle\boldsymbol{U}_{\text{f}}^{(\text{inner,
D})}(s)+\boldsymbol{U}_{\text{f}}^{\text{(int, D)}}(s).$

The term $\boldsymbol{U}_{\text{f}}^{(\text{inner, D})}(s)$ is given for all
$s$ by

$\displaystyle
8\pi\mu\boldsymbol{U}_{\text{f}}^{(\text{inner,D})}(s)=\frac{2}{3}\left(\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)\begin{cases}\dfrac{1}{4}-\dfrac{{\hat{a}}^{2}}{2s^{2}}-\dfrac{{\hat{a}}^{2}}{2(L-s)^{2}}&2{\hat{a}}<s<L-2{\hat{a}}\\\
\dfrac{1}{8}-\dfrac{{\hat{a}}^{2}}{2(L-s)^{2}}&s\leq 2{\hat{a}}\\\
\dfrac{1}{8}-\dfrac{{\hat{a}}^{2}}{2s^{2}}&s\geq L-2{\hat{a}}\end{cases}$
(G.9)

The nearly-singular doublet integral is handled in exactly the same way as the
Stokeslet one, but the singular behavior is different, namely,

$\displaystyle\boldsymbol{U}_{\text{f}}^{\text{(int,
D)}}(s)=\int_{D(s)}\frac{(s^{\prime}-s)}{|s^{\prime}-s|^{3}}\,\boldsymbol{g}_{D}(s^{\prime},s)\,ds^{\prime},$
$\displaystyle\text{where}\qquad\boldsymbol{g}_{D}(s^{\prime},s)=\left(\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{f}\left(s^{\prime}\right)-\frac{1}{8\pi\mu}\left(\frac{\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s^{\prime}-s|^{3}}\right)\boldsymbol{f}(s)\right)\frac{|s^{\prime}-s|^{3}}{(s^{\prime}-s)}$
(G.10)

is a smooth function with the finite limit

$\lim_{s^{\prime}\rightarrow
s}\boldsymbol{g}_{D}(s^{\prime},s)=\frac{1}{8\pi\mu}\left(\left(\boldsymbol{I}-3\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}^{\prime}(s)+\frac{1}{2}\left(\boldsymbol{\tau}(s)\partial_{s}\boldsymbol{\tau}(s)+\partial_{s}\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)\right).$
(G.11)

If we expand $\boldsymbol{g}_{D}$ in terms of Chebyshev polynomials, the
efficient evaluation of (G.10) requires precomputing integrals of the form

$q_{k}^{\text{(D)}}(\eta)=\int_{D^{\prime}(\eta)}\frac{(\eta^{\prime}-\eta)}{|\eta^{\prime}-\eta|^{3}}T_{k}\left(\eta^{\prime}\right)\,d\eta^{\prime},$
(G.12)

for all $\eta$ on the Chebyshev grid. Notice that the integrals (G.12) are not
defined if $\eta^{\prime}=\eta$ is included in the integration domain, but
$D^{\prime}(\eta)$ does _not_ contain $\eta$.

###  G.3 Rot-trans mobility

The long-range integral for the rot-trans mobility (72) is given by an
integral of the rotlet kernel over the domain $D(s)$. Using singularity
subtraction, we compute this integral as

$\displaystyle\Omega_{\text{f}}^{\text{(R)}}(s)=$
$\displaystyle\int_{D(s)}\frac{\left(\boldsymbol{f}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)\cdot\boldsymbol{\tau}(s)}{R(s^{\prime})^{3}}ds^{\prime}$
$\displaystyle=$
$\displaystyle\int_{D(s)}\left(\frac{\left(\partial_{s}\boldsymbol{\tau}(s)\times\boldsymbol{f}(s)\right)\cdot\boldsymbol{\tau}(s)}{2|s-s^{\prime}|}\right)\,ds^{\prime}$
$\displaystyle+$
$\displaystyle\int_{D(s)}\left(\frac{\left(\boldsymbol{f}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)\cdot\boldsymbol{\tau}(s)}{R(s^{\prime})^{3}}-\frac{\left(\partial_{s}\boldsymbol{\tau}(s)\times\boldsymbol{f}(s)\right)\cdot\boldsymbol{\tau}(s)}{2|s-s^{\prime}|}\right)\,ds^{\prime}$
$\displaystyle:=$ $\displaystyle\Omega_{\text{f}}^{(\text{inner,
R})}(s)+\Omega_{\text{f}}^{(\text{int, R})}(s),$ (G.13)

where
$\boldsymbol{R}(s^{\prime})=\boldsymbol{X}(s)-\boldsymbol{X}(s^{\prime})$ and
$R=\left\|\boldsymbol{R}\right\|$. Since the leading order singularity in
$\Omega_{\text{f}}^{(\text{inner, R})}$ of $1/|s-s^{\prime}|$ is the same as
that of the Stokeslet, the expression for $\Omega_{\text{f}}^{(\text{inner,
R})}$ is also given by (71), but with
$\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)$
replaced by
$1/2\left(\partial_{s}\boldsymbol{\tau}(s)\times\boldsymbol{f}(s)\right)\cdot\boldsymbol{\tau}(s)=1/2\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)\cdot\boldsymbol{f}(s)$.

The nearly singular integral $\Omega_{\text{f}}^{(\text{int, R})}(s)$ in
(G.13) can then be put in the form (G.3) (with scalar
$g_{\text{rt}}(s,s^{\prime})$) by setting

$g_{\text{rt}}(s,s^{\prime})=\left(\frac{|s-s^{\prime}|}{R(s^{\prime})^{3}}\left(\boldsymbol{R}(s^{\prime})\times\boldsymbol{\tau}(s)\right)\cdot\boldsymbol{f}(s^{\prime})-\frac{\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)\cdot\boldsymbol{f}(s)}{2}\right)\frac{1}{s^{\prime}-s},$
(G.14)

which has the finite limit

$\lim_{s^{\prime}\rightarrow
s}g_{\text{rt}}(s,s^{\prime})=-\frac{1}{2}\partial_{s}\boldsymbol{f}(s)\cdot\left(\partial_{s}\boldsymbol{\tau}(s)\times\boldsymbol{\tau}(s)\right)-\frac{1}{6}\boldsymbol{f}(s)\cdot\left(\partial_{s}^{2}{\boldsymbol{\tau}}(s)\times\boldsymbol{\tau}(s)\right).$
(G.15)

Now the scalar function $g_{\text{rt}}(s,s^{\prime})$ can be expanded in a
Chebyshev series as in (G.6), and the precomputed integrals
$q^{(S)}_{k}(\eta)$ in (G.7) used once again to evaluate the second integral
in (G.13).

###  G.4 Trans-rot mobility

For the trans-rot mobility (74), the singularity subtraction approach is to
set

$\displaystyle\boldsymbol{U}_{\text{n}}^{(\text{R})}(s)=$
$\displaystyle\int_{D(s)}\frac{\left(\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)n^{\parallel}(s^{\prime})}{R(s^{\prime})^{3}}ds^{\prime}$
$\displaystyle=$
$\displaystyle\int_{D(s)}\left(\frac{\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)n^{\parallel}(s)}{2|s-s^{\prime}|}\right)\,ds^{\prime}$
$\displaystyle+$
$\displaystyle\int_{D(s)}\left(\frac{\left(\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)n^{\parallel}(s^{\prime})}{R(s^{\prime})^{3}}-\frac{\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)n^{\parallel}(s)}{2|s-s^{\prime}|}\right)\,ds^{\prime}$
$\displaystyle:=$ $\displaystyle\boldsymbol{U}_{\text{n}}^{(\text{inner,
R})}(s)+\boldsymbol{U}_{\text{n}}^{(\text{int, R})}(s).$ (G.16)

where $\boldsymbol{U}_{\text{n}}^{(\text{inner, R})}$ is given by (71), but
with
$\left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{f}(s)$
replaced by
$1/2\left(\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)\right)n^{\parallel}(s)$.
The nearly singular integral $\boldsymbol{U}_{\text{n}}^{(\text{int, R})}(s)$
can then be put in the form (G.3) by setting

$\displaystyle\boldsymbol{g}_{\text{tr}}(s,s^{\prime})$
$\displaystyle=\left(\frac{|s-s^{\prime}|}{R(s^{\prime})^{3}}\left(\boldsymbol{\tau}(s^{\prime})\times\boldsymbol{R}(s^{\prime})\right)n^{\parallel}(s^{\prime})-\left(\frac{\boldsymbol{\tau}(s)\times\partial_{s}\boldsymbol{\tau}(s)}{2}\right)n^{\parallel}(s)\right)\frac{1}{s^{\prime}-s},$
(G.17) where $\displaystyle\lim_{s^{\prime}\rightarrow
s}\boldsymbol{g}_{\text{tr}}(s,s^{\prime})=-\frac{1}{2}\partial_{s}n^{\parallel}(s)\left(\partial_{s}\boldsymbol{\tau}(s)\times\boldsymbol{\tau}(s)\right)-\frac{1}{3}n^{\parallel}(s)\left(\partial_{s}^{2}{\boldsymbol{\tau}}(s)\times\boldsymbol{\tau}(s)\right).$
(G.18)

###  G.5 Rot-rot mobility

Although we use asymptotic evaluation for the rot-rot mobility as discussed in
Section 2.1.3, it is straightforward to derive a nearly singular quadrature
scheme to evaluate the RPY integral for
$\boldsymbol{\mathcal{M}}_{\text{rr}}$. Substituting the tensor
$\mathbb{M}_{\text{rr}}$ from (A.4) into the integral (10) and using the
approximation
$\left\|\boldsymbol{X}(s)-\boldsymbol{X}(s^{\prime})\right\|\approx
s-s^{\prime}$ when $s-s^{\prime}=\mathcal{O}({\hat{a}})$, the mobility is
given by

$\displaystyle
8\pi\mu\boldsymbol{\Psi}_{\text{n}}(s)=-\frac{1}{2}\int_{D(s)}\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^{\prime})\right)\boldsymbol{n}(s^{\prime})\,ds^{\prime}$
(G.19)
$\displaystyle+\frac{1}{{\hat{a}}^{3}}\int_{D^{c}(s)}\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 of these integrals can be evaluated to high accuracy using the
doublet scheme of Section G.2, while the second integral can be computed using
direct Gauss-Legendre quadrature on $D^{c}(s)$ (split into two pieces). The
parallel angular velocity $\Psi_{\text{n}}^{\parallel}$ can then be obtained
by taking the inner product of $\boldsymbol{\Psi}_{\text{n}}$ with
$\boldsymbol{\tau}$.

###  G.6 Convergence of quadratures

Figure 13:  Convergence of the nearly singular quadrature schemes for the
(left) trans-trans mobility (68) and trans-rot mobility (74), and (right) rot-
trans mobility (72) and rot-rot mobility (76) (here we are comparing the the
full RPY integral of doublets as in (10)). We use the fiber (87) with $q=7$
and compute force and torque on the collocation grid of $N$ points using (82),
(85), and (86) with $\kappa=\gamma=1$. We then compute the velocities
$\boldsymbol{U}_{\text{f}}$, $\boldsymbol{U}_{\text{n}}$,
$\Omega_{\text{f}}^{\parallel}$, and $\Omega_{\text{n}}^{\parallel}$ on the
collocation grid by applying the discrete mobility matrices given in Section
4.2. In all cases, we compute a reference solution by upsampled direct
quadrature and report the _relative_ $L^{2}$ errors in the velocity.

We next test the accuracy of the (action of the) mobility matrices
$\boldsymbol{M}_{\text{tt}}$, $\boldsymbol{M}_{\text{tr}}$,
$\boldsymbol{M}_{\text{rt}}$, and $\boldsymbol{M}_{\text{rr}}$ by comparing
against upsampled direct quadrature. We use the fiber (87) with $q=7$, set up
a collocation grid of $N$ points, and compute the force
$\boldsymbol{f}=\boldsymbol{f}^{(\kappa)}+\boldsymbol{f}^{(\gamma)}$ and
torque density $n^{\parallel}$ using (82), (85), and (86) with
$\kappa=\gamma=1$. We apply the forward mobility operators as discretized in
Section 4.2 to obtain the (translational and rotational) velocity at each of
the collocation points, and compare the result to a reference solution which
uses direct oversampled quadrature. The relative error is computed as the
difference between the two, normalized by the $L^{2}$ norm of the reference
solution.

Figure 13 shows the results for the trans-trans mobility (68) and trans-rot
mobility (74) (left), and rot-trans mobility (72) and rot-rot mobility (76)
(right) for two different values of $\hat{\epsilon}$, and various values of
$N$ (see Section 4.2 for discretization details). In the case of the trans-
trans, rot-trans, and trans-rot mobilities, we obtain spectral convergence and
about 3 digits of accuracy when $N=32$. For $\hat{\epsilon}=10^{-2}$, $N_{2}$
must be at least 6 to see this accuracy, but for $\hat{\epsilon}=10^{-3}$ a
smaller $N_{2}\approx 4$ is acceptable since the relative size of the domain
$|s-s^{\prime}|\leq 2{\hat{a}}$ shrinks by a factor of 10. For
$\hat{\epsilon}=10^{-2}$, we previously found $N_{2}\geq 8$ to give positive
eigenvalues in the mobility, so accuracy is a weaker restriction than the
eigenvalues of $\boldsymbol{M}_{\text{tt}}$ being strictly positive.

In Fig. 13, we also show the accuracy of the rot-rot mobility as defined in
(76) relative to the actual RPY integral of doublets as defined in (10) and
(A.4). This “convergence” plot is fundamentally different from the others,
since in this case we are measuring the error of the asymptotic formula _and_
its approximation by a finite number of collocation points, rather than just
the latter. Therefore, the relative error saturates around $N=16$ at a value
that strongly depends on $\hat{\epsilon}$. For $\hat{\epsilon}=10^{-2}$, the
relative error saturates around $2\times 10^{-2}$, while for
$\hat{\epsilon}=10^{-3}$ the saturated relative error is about $4\times
10^{-4}$. This reflects the fact that the full asymptotic formula for the rot-
rot mobility [46, Appendix C.3] is $\mathcal{O}(\hat{\epsilon}^{2})$ in the
fiber interior ($2{\hat{a}}\leq s\leq L-2{\hat{a}}$) and
$\mathcal{O}(\hat{\epsilon})$ at the fiber endpoints.

## References

  * [1] Ludvig af Klinteberg and Alex H Barnett.  Accurate quadrature of nearly singular line integrals in two and three dimensions by singularity swapping.  BIT Numerical Mathematics, pages 1–36, 2020.
  * [2] Ludvig af Klinteberg, Davoud Saffar Shamshirgar, and Anna-Karin Tornberg.  Fast ewald summation for free-space stokes potentials.  Research in the Mathematical Sciences, 4(1):1–32, 2017.
  * [3] Bhagwan D Agarwal, Lawrence J Broutman, and K Chandrashekhara.  Analysis and performance of fiber composites.  John Wiley & Sons, 2017.
  * [4] Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, and Peter Walter.  Molecular biology of the cell.  Garland Science, 2002.
  * [5] Helge I Andersson, Elena Celledoni, Laurel Ohm, Brynjulf Owren, and Benjamin K Tapley.  An integral model based on slender body theory, with applications to curved rigid fibers.  Physics of Fluids, 33(4):041904, 2021.
  * [6] Jared L Aurentz and Lloyd N Trefethen.  Block operators and spectral discretizations.  SIAM Review, 59(2):423–446, 2017.
  * [7] Florencio Balboa Usabiaga, Bakytzhan Kallemov, Blaise Delmotte, Amneet Bhalla, Boyce Griffith, and Aleksandar Donev.  Hydrodynamics of suspensions of passive and active rigid particles: a rigid multiblob approach.  Communications in Applied Mathematics and Computational Science, 11(2):217–296, 2017.
  * [8] Miklós Bergou, Max Wardetzky, Stephen Robinson, Basile Audoly, and Eitan Grinspun.  Discrete elastic rods.  In ACM SIGGRAPH 2008 papers, pages 1–12. 2008.
  * [9] Jacob J Blum and Michael Hines.  Biophysics of flagellar motility.  Quarterly reviews of biophysics, 12(2):103–180, 1979.
  * [10] John F Brady and Georges Bossis.  Stokesian dynamics.  Annual review of fluid mechanics, 20(1):111–157, 1988.
  * [11] Christopher Brennen and Howard Winet.  Fluid mechanics of propulsion by cilia and flagella.  Annual Review of Fluid Mechanics, 9(1):339–398, 1977.
  * [12] 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.
  * [13] Isaac R Bruss, Heena K Mutha, Katherine Stoll, Brent Collins, Vinh Nguyen, David JD Carter, Michael P Brenner, and Kasey J Russell.  Twirling, whirling, and tensioning: Plectoneme formation and suppression in flexible filaments.  Physical Review Research, 1(3):032020, 2019.
  * [14] Mireille Maria Anna Elisabeth Claessens, C Semmrich, L Ramos, and AR Bausch.  Helical twist controls the thickness of f-actin bundles.  Proceedings of the National Academy of Sciences, 105(26):8819–8822, 2008.
  * [15] Ricardo Cortez.  The method of regularized stokeslets.  SIAM Journal on Scientific Computing, 23(4):1204–1225, 2001.
  * [16] 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.
  * [17] 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.
  * [18] NJ De Mestre and WB Russel.  Low-reynolds-number translation of a slender cylinder near a plane wall.  Journal of Engineering Mathematics, 9(2):81–91, 1975.
  * [19] Steven Delong, Boyce E Griffith, Eric Vanden-Eijnden, and Aleksandar Donev.  Temporal integrators for fluctuating hydrodynamics.  Physical Review E, 87(3):033302, 2013.
  * [20] Ellis Harold Dill.  Kirchhoff’s theory of rods.  Archive for History of Exact Sciences, pages 1–23, 1992.
  * [21] Tobin A Driscoll and Nicholas Hale.  Rectangular spectral collocation.  IMA Journal of Numerical Analysis, 36(1):108–132, 2015.
  * [22] Olivia Du Roure, Anke Lindner, Ehssan N Nazockdast, and Michael J Shelley.  Dynamics of flexible fibers in viscous flows and fluids.  Annu. Rev. Fluid. Mech., 51:539–572, 2019.
  * [23] Andrew M Fiore, Florencio Balboa Usabiaga, Aleksandar Donev, and James W Swan.  Rapid sampling of stochastic displacements in brownian dynamics simulations.  The Journal of chemical physics, 146(12):124116, 2017.
  * [24] 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.
  * [25] Carlos Floyd, Haoran Ni, Ravinda Gunaratne, Radek Erban, and Garegin A Papoian.  On stretching, bending, shearing and twisting of actin filaments i: Variational models.  arXiv preprint arXiv:2112.01480, 2021.
  * [26] Mohit Garg and Ajeet Kumar.  A slender body theory for the motion of special cosserat filaments in stokes flow.
  * [27] Zydrunas Gimbutas and Leslie Greengard.  Computational software: Simple fmm libraries for electrostatics, slow viscous flow, and frequency-domain wave propagation.  Communications in Computational Physics, 18(2):516–528, 2015.
  * [28] Thomas Götz.  Interactions of fibers and flow: asymptotics, theory and numerics.  dissertation. de, 2001.
  * [29] Weicheng Huang and M Khalid Jawed.  Numerical simulation of bundling of helical elastic rods in a viscous fluid.  Computers & Fluids, page 105038, 2021.
  * [30] Kenta Ishimoto and Eamonn A Gaffney.  An elastohydrodynamical simulation study of filament and spermatozoan swimming driven by internal couples.  IMA Journal of Applied Mathematics, 83(4):655–679, 2018.
  * [31] Mehdi Jabbarzadeh and Henry C Fu.  A numerical method for inextensible elastic filaments in viscous fluids.  Journal of Computational Physics, page 109643, 2020.
  * [32] Robert E Johnson.  An improved slender-body theory for stokes flow.  Journal of Fluid Mechanics, 99(2):411–431, 1980.
  * [33] 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.
  * [34] 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.
  * [35] Joseph B Keller and Sol I Rubinow.  Slender-body theory for slow viscous flow.  Journal of Fluid Mechanics, 75(4):705–714, 1976.
  * [36] Lyndon Koens.  Tubular-body theory for viscous flows.  Physical Review Fluids, 7(3):034101, 2022.
  * [37] Cornelius Lanczos.  The variational principles of mechanics.  Courier Corporation, 2012.
  * [38] Eric Lauga and Thomas R Powers.  The hydrodynamics of swimming microorganisms.  Reports on Progress in Physics, 72(9):096601, 2009.
  * [39] 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.
  * [40] Sookkyung Lim, Anca Ferent, X Sheldon Wang, and Charles S Peskin.  Dynamics of a closed rod with twist and bend in fluid.  SIAM Journal on Scientific Computing, 31(1):273–302, 2008.
  * [41] 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.
  * [42] Dag Lindbo and Anna-Karin Tornberg.  Fast and spectrally accurate summation of 2-periodic stokes potentials.  arXiv preprint arXiv:1111.1815, 2011.
  * [43] Sune Lomholt and Martin R Maxey.  Force-coupling method for particulate two-phase flow: Stokes flow.  J. Comput. Phys., 184(2):381–405, 2003.
  * [44] Rui Ma and Julien Berro.  Structural organization and energy storage in crosslinked actin assemblies.  PLoS computational biology, 14(5):e1006150, 2018.
  * [45] Yi Man, William Page, Robert J Poole, and Eric Lauga.  Bundling of elastic filaments induced by hydrodynamic interactions.  Physical Review Fluids, 2(12):123101, 2017.
  * [46] Ondrej Maxian and Aleksandar Donev.  Slender body theories for rotating filaments.  arXiv preprint arXiv:2203.12059, 2022.
  * [47] Ondrej Maxian, Aleksandar Donev, and Alex Mogilner.  Interplay between brownian motion and cross-linking kinetics controls bundling dynamics in actin networks.  Biophysical Journal, 2022.
  * [48] 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.
  * [49] Ondrej Maxian, Raul P Peláez, Alex Mogilner, and Aleksandar Donev.  Simulations of dynamically cross-linked actin networks: morphology, rheology, and hydrodynamic interactions.  PLOS Computational Biology, 17(12):e1009240, 2021.
  * [50] 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.
  * [51] Ehssan Nazockdast, Abtin Rahimian, Denis Zorin, and Michael Shelley.  A fast platform for simulating semi-flexible fiber suspensions applied to cell mechanics.  J. Comput. Phys., 329:173–209, 2017.
  * [52] Frank TM Nguyen and Michael D Graham.  Impacts of multiflagellarity on stability and speed of bacterial locomotion.  Physical Review E, 98(4):042419, 2018.
  * [53] Anton Oberbeck.  Uber stationare flussigkeitsbewegungen mit berucksichtigung der inner reibung.  J. reine angew. Math., 81:62–80, 1876.
  * [54] Sarah D Olson, Sookkyung Lim, and Ricardo Cortez.  Modeling the dynamics of an elastic rod with intrinsic curvature and twist using a regularized stokes formulation.  Journal of Computational Physics, 238:169–187, 2013.
  * [55] Lorenzo Pareschi and Giovanni Russo.  Implicit–explicit runge–kutta schemes and applications to hyperbolic systems with relaxation.  Journal of Scientific computing, 25(1):129–155, 2005.
  * [56] Charles S Peskin.  Flow patterns around heart valves: a numerical method.  J. Comput. Phys., 10(2):252–271, 1972.
  * [57] Charles S Peskin.  The immersed boundary method.  Acta Numer., 11:479–517, 2002.
  * [58] Thomas R Powers.  Dynamics of filaments and membranes in a viscous fluid.  Reviews of Modern Physics, 82(2):1607, 2010.
  * [59] Michael Reichert.  Hydrodynamic interactions in colloidal and biological systems.  PhD thesis, 2006.
  * [60] Jens Rotne and Stephen Prager.  Variational treatment of hydrodynamic interaction in polymers.  J. Chem. Phys., 50(11):4831–4837, 1969.
  * [61] Mordecai B Rubin.  Cosserat theories: shells, rods and points, volume 79.  Springer Science & Business Media, 2013.
  * [62] 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.
  * [63] Michael J Shelley and Tetsuji Ueda.  The nonlocal dynamics of stretching, buckling filaments.  Advances in multi-fluid flows, pages 415–425, 1996.
  * [64] Michael J Shelley and Tetsuji Ueda.  The stokesian hydrodynamics of flexing, stretching filaments.  Physica D: Nonlinear Phenomena, 146(1-4):221–245, 2000.
  * [65] James W Swan and John F Brady.  Simulation of hydrodynamically interacting particles near a no-slip boundary.  Physics of Fluids, 19(11):113306, 2007.
  * [66] Anna-Karin Tornberg.  Accurate evaluation of integrals in slender-body formulations for fibers in viscous flow.  arXiv preprint arXiv:2012.12585, 2020.
  * [67] 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.
  * [68] Lloyd N Trefethen.  Spectral methods in MATLAB, volume 10.  Siam, 2000.
  * [69] Reinhard Vogel.  The bacterial flagellum: Modeling the dynamics of the elastic filament and its transition between polymorphic helical forms.  2012\.
  * [70] Hirofumi Wada and Roland R Netz.  Non-equilibrium hydrodynamics of a rotating filament.  EPL (Europhysics Letters), 75(4):645, 2006.
  * [71] Hirofumi Wada and Roland R Netz.  Model for self-propulsive helical filaments: kink-pair propagation.  Physical review letters, 99(10):108102, 2007.
  * [72] 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.
  * [73] Benjamin J Walker and Eamonn A Gaffney.  Regularised non-uniform segments and efficient no-slip elastohydrodynamics.  Journal of Fluid Mechanics, 915, 2021.
  * [74] Timothy A Westwood and Eric E Keaveny.  Coordinated motion of active filaments on spherical surfaces.  arXiv preprint arXiv:2109.13578, 2021.
  * [75] 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.
  * [76] Wen Yan and Michael Shelley.  Universal image systems for non-periodic and periodic stokes flows above a no-slip wall.  Journal of Computational Physics, 375:263–270, 2018.
  * [77] PJ Zuk, E Wajnryb, KA Mizerski, and P Szymczak.  Rotne–prager–yamakawa approximation for different-sized particles in application to macromolecular bead models.  Journal of Fluid Mechanics, 741, 2014.