# An integral model based on slender body theory, with applications to curved
rigid fibers

Helge I. Andersson  Department of Energy and Process Engineering, The
Norwegian University of Science and Technology, 7491 Trondheim, Norway  Elena
Celledoni  Department of Mathematical Sciences, The Norwegian University of
Science and Technology, 7491 Trondheim, Norway  Laurel Ohm  Courant Institute
of Mathematical Sciences, New York University, New York, New York 10012, USA
Brynjulf Owren  Department of Mathematical Sciences, The Norwegian University
of Science and Technology, 7491 Trondheim, Norway  Benjamin K. Tapley
Department of Mathematical Sciences, The Norwegian University of Science and
Technology, 7491 Trondheim, Norway

###### Abstract

We propose a novel integral model describing the motion of curved slender
fibers in viscous flow, and develop a numerical method for simulating dynamics
of rigid fibers. The model is derived from nonlocal slender body theory (SBT),
which approximates flow near the fiber using singular solutions of the Stokes
equations integrated along the fiber centerline. In contrast to other models
based on (singular) SBT, our model yields a smooth integral kernel which
incorporates the (possibly varying) fiber radius naturally. The integral
operator is provably negative definite in a non-physical idealized geometry,
as expected from PDE theory. This is numerically verified in physically
relevant geometries. We propose a convergent numerical method for solving the
integral equation and discuss its convergence and stability. The accuracy of
the model and method is verified against known models for ellipsoids. Finally,
a fast algorithm for computing dynamics of rigid fibers with complex
geometries is developed.

##  1 Introduction

The dynamics of thin fibers immersed in fluid play an important role in many
biological and engineering processes, including microorganism propulsion [30,
5, 18, 33], rheological properties of fiber suspensions used to create
composite materials [29, 12, 8], and deposition of microplastics in the ocean
[20]. Here the term ‘fiber’ is used to refer to a particle with a very large
aspect ratio. In many of the applications mentioned, the cross sectional
radius of the fiber is small compared to the length scales of the surrounding
fluid, which can be well approximated locally by Stokes flow. This allows for
the development of computationally tractable mathematical models describing
the interaction between the fiber and the surrounding fluid.

Slender body theory [10, 15, 16, 19, 37] is a popular tool for reducing
computational costs of simulating these thin fibers by approximating the
fibers as one dimensional curves. Here we propose an integral model based on
slender body theory which involves a smooth kernel and incorporates the
(possibly varying) fiber radius in a natural way. Since the integral kernels
are smooth, the model resembles the method of regularized Stokeslets [7] with
an arclength-dependent regularization similar to [40]; however, we derive our
model from usual (singular) Stokeslets and doublets. As such, we avoid
introducing additional parameters into the basic model. The model relies on
the asymptotic cancellation of angular-dependent terms along the fiber surface
(see Section 3 for details), leaving an expression that retains a dependence
on the fiber radius in a natural way.

Furthermore, we calculate the spectrum of our integral operator in the toy
scenario of a straight-but-periodic fiber with constant radius. In this toy
geometry, our integral operator is negative definite, as is the well-posed PDE
operator of [24, 25] which it is designed to approximate (see [22]). This is
in contrast to other models based on (non-regularized) slender body theory
which rely on further asymptotic expansion with respect to the fiber radius
[15, 16, 19]. These models exhibit an instability as the eigenvalues of the
operator cross zero at a high but finite wavenumber.

The model we derive initially yields a first-kind Fredholm integral equation
for the force density along the fiber centerline. Such integral equations are
known for being ill-posed [17, Chapt. 15.1], as they do not necessarily have a
bounded inverse at the continuous level. Numerical discretization alone can
provide sufficient regularization to invert first-kind integral equations at
the discrete level, but to make our model more suitable for inversion, we use
an integral identity to regularize the expression into a second-kind equation.
The second-kind regularization preserves the asymptotic accuracy of the model
while improving the conditioning and invertibility of the corresponding
numerical method. The regularization also serves to ensure that the
discretized operator is negative definite, even in the presence of numerical
errors, by bounding the spectrum away from zero. We distinguish this type of
regularization from the method of regularized Stokeslets, since our
regularization is not a key component of the model derivation. In particular,
we can directly compare our model _with_ regularization to our model _without_
, which we will do repeatedly throughout the paper. We also distinguish this
regularization from the procedure used by Tornberg and Shelley [37], since we
are not correcting for a high wavenumber instability. This allows us to
compare the numerical behavior of our regularized and unregularized models at
the discrete level even for very fine discretization. Moreover, the
regularization used here affects all directions (both normal and tangent to
the slender body centerline) in the same way.

The solution of the resulting second-kind Fredholm integral equation is a
force density along the slender body centerline which we integrate to find the
total force and torque on the rigid fiber. We present a numerical method based
on the Nyström method for solving second-kind Fredholm integral equations [1,
Chapt. 12.4]. Based on the results of [1], we show that the method is
convergent and that the accuracy of the method depends on the quadrature
error. Numerical tests confirm the convergence results. Not surprisingly, we
note significant improvements in the conditioning of the second-kind versus
first-kind formulation of the model. We also numerically verify the spectral
properties of the model in different geometries.

We develop an algorithm for dynamic simulations of a rigid fiber. The rigidity
of the fiber can be exploited such that only matrix-vector products need to be
performed within the time loop, resulting in a fast algorithm for computing
dynamics. We compare the dynamics of our model to the well-studied dynamics of
a slender prolate spheroid [14, 3, 6]. We then apply our model to compare the
dynamics of curved fibers whose centerlines deviate randomly from straight
lines by varying magnitudes.

The structure of the paper is as follows. Section 2 presents the slender body
model, which is derived in greater detail and justified via spectral
comparisons with other slender body theories in Section 3. In Section 4 we
discuss a method for numerically solving Fredholm integral equations and
integrating the result, and demonstrate the convergence of the method for our
model. Section 5 outlines a fast algorithm for computing the dynamics of a
rigid slender fiber in viscous flow. We apply the dynamical algorithm to
simulate the dynamics of fibers with complex shapes. Finally, we comment on
conclusions and outlook for the model in Section 6.

###  1.1 Fiber geometry

We begin by introducing some notation for the slender geometries considered
throughout the paper. Fix $\epsilon$, $L$ with $0<\epsilon\ll L$ and let
$\bm{X}_{\text{ext}}:[-\sqrt{L^{2}+\epsilon^{2}},\sqrt{L^{2}+\epsilon^{2}}]\to\mathbb{R}^{3}$
denote the coordinates of a $C^{2}$ curve in $\mathbb{R}^{3}$, parameterized
by arclength $s$. Defining $\bm{e}_{\rm
s}(s)=\frac{d\bm{X}_{\text{ext}}}{ds}/\left\lvert\frac{d\bm{X}_{\text{ext}}}{ds}\right\rvert$,
the unit tangent vector to $\bm{X}_{\text{ext}}(s)$, we parameterize points
near $\bm{X}_{\text{ext}}(s)$ with respect to the orthonormal frame
$(\bm{e}_{\rm s}(s),\bm{e}_{n_{1}}(s),\bm{e}_{n_{2}}(s))$ defined in [25].
Letting

$\bm{e}_{r}(s,\theta):=\cos\theta\bm{e}_{n_{1}}(s)+\sin\theta\bm{e}_{n_{2}}(s),$

we define the slender body $\Sigma_{\epsilon}$ as

$\Sigma_{\epsilon}:=\big{\\{}\bm{x}\in\mathbb{R}^{3}\;:\;\bm{x}=\bm{X}_{\text{ext}}(s)+\rho\bm{e}_{r}(s,\theta),\;\rho<\epsilon
r(s),\;s\in[-\sqrt{L^{2}+\epsilon^{2}},\sqrt{L^{2}+\epsilon^{2}}]\big{\\}}.$
(1)

Here the radius function $r\in
C^{2}(-\sqrt{L^{2}+\epsilon^{2}},\sqrt{L^{2}+\epsilon^{2}})$ is required to
satisfy $0<r(s)\leq 1$ for each
$s\in(-\sqrt{L^{2}+\epsilon^{2}},\sqrt{L^{2}+\epsilon^{2}})$, and $r(s)$ must
decay smoothly to zero at the fiber endpoints $\pm\sqrt{L^{2}+\epsilon^{2}}$.
There are many admissible radius functions $r$ which can be considered. For
the simulations in this paper, we will use a thin prolate spheroid as our
geometrical model for a slender fiber. In this case, the radius function
$r(s)$ is given by

$r(s)=\frac{1}{\sqrt{L^{2}+\epsilon^{2}}}\sqrt{L^{2}+\epsilon^{2}-s^{2}}.$ (2)

We consider the subset

$\bm{X}:=\\{\bm{X}_{\text{ext}}(s)\;:\;-L\leq s\leq L\\}$ (3)

extending from focus to focus of the prolate spheroid (2), and define
$\bm{X}(s)$ to be the effective centerline of the slender body so that
$r=O(\epsilon)$ at the effective endpoints $s=\pm L$.

The slender body model described in Section 2 may also be used in the case of
a closed curve, in which case we take $\bm{X}(L)=\bm{X}(-L)$ and consider
$s\in\mathbb{R}/2L$. We may take the radius function $r\equiv 1$ in this case.

##  2 Slender body model

To describe the motion of the thin fiber $\Sigma_{\epsilon}$ (1) in Stokes
flow, we will use an expression derived from nonlocal slender body theory [10,
15, 16, 37]. Letting $\bm{f}(s,t)$ denote the force per unit length exerted by
the fiber on the surrounding fluid at time $t$, we approximate the velocity
$\frac{\partial\bm{X}}{\partial t}$ of the fiber relative to a given
background flow $\bm{u}_{0}$ by

$\displaystyle 8\pi\mu\bigg{(}\frac{\partial\bm{X}}{\partial
t}-\bm{u}_{0}(\bm{X}(s,t),t)\bigg{)}$
$\displaystyle=-2\log\eta\,\bm{f}(s,t)-\int_{-L}^{L}\bigg{(}\bm{S}_{\epsilon,\eta}+\frac{\epsilon^{2}r^{2}(s^{\prime})}{2}\bm{D}_{\epsilon}\bigg{)}\bm{f}(s^{\prime},t)\,ds^{\prime},$
(4) $\displaystyle\bm{S}_{\epsilon,\eta}(s,s^{\prime},t)$
$\displaystyle=\frac{{\bf
I}}{(|\overline{\bm{X}}|^{2}+\eta^{2}\epsilon^{2}r^{2}(s))^{1/2}}+\frac{\overline{\bm{X}}\overline{\bm{X}}^{\rm
T}}{(|\overline{\bm{X}}|^{2}+\epsilon^{2}r^{2}(s))^{3/2}}$ (5)
$\displaystyle\bm{D}_{\epsilon}(s,s^{\prime},t)$ $\displaystyle=\frac{{\bf
I}}{(|\overline{\bm{X}}|^{2}+\epsilon^{2}r^{2}(s))^{3/2}}-\frac{3\overline{\bm{X}}\overline{\bm{X}}^{\rm
T}}{(|\overline{\bm{X}}|^{2}+\epsilon^{2}r^{2}(s))^{5/2}}$ (6)

where $\overline{\bm{X}}(s,s^{\prime},t)=\bm{X}(s,t)-\bm{X}(s^{\prime},t)$.
Here $\eta>1$ is a regularization parameter chosen so that (4) is a second-
kind Fredholm equation for $\bm{f}$. Notice that $\eta$ must also appear in
the first term of $\bm{S}_{\epsilon,\eta}$ in order to retain the asymptotic
consistency of the model (4). This is due to an integral identity (13) used to
convert the integral model from a first-kind equation for for $\bm{f}$. The
model accounts for a varying radius $r(s)$ through the denominators of each
term as well as the coefficient of $\bm{D}_{\epsilon}$. Note that since $r(s)$
is nonzero for $-L\leq s\leq L$, the integral kernel is smooth for each
$s\in[-L,L]$. We provide a more detailed derivation of (4)–(6) in Section 3.

The model given by equations (4)–(6) and the analysis in Section 3 can be used
to describe both flexible and rigid fibers. In Section 5 we apply our model to
the dynamics of a rigid fiber, since the invertibility properties of (4)–(6)
make the model especially suitable for simulating rigid filaments.

In the case of a rigid fiber, at each time $t$ we additionally impose the
constraint

$\frac{\partial\bm{X}}{\partial t}=\bm{v}+\bm{\omega}\times\bm{X}(s),$ (7)

where $\bm{v}$, $\bm{\omega}\in\mathbb{R}^{3}$ are the given linear and
angular velocity of the fiber (see [11, 23, 36]). We then use (4) to solve for
the total force $\bm{F}(t)$ and torque $\bm{T}(t)$ exerted on the slender body
at time $t$ via

$\int_{-L}^{L}\bm{f}(s,t)\,ds=\bm{F}(t),\qquad\int_{-L}^{L}\bm{X}(s,t)\times\bm{f}(s,t)=\bm{T}(t).$
(8)

Note that solving for $\bm{F}$ and $\bm{T}$ involves inverting the expression
(4), so we are particularly concerned with the invertibility of the integral
equation.

##  3 Derivation and justification of the slender body model

Our model for the motion of the fiber is based on slender body theory [10, 15,
16, 37]. According to slender body theory, the fluid velocity $\bm{u}^{\rm
SB}(\bm{x},t)$ at any point $\bm{x}$ away from the fiber centerline
$\bm{X}(s,t)$ is approximated by the integral expression

$\displaystyle 8\pi\mu\big{(}\bm{u}^{{\rm
SB}}(\bm{x},t)-\bm{u}_{0}(\bm{x},t)\big{)}$
$\displaystyle=-\int_{-L}^{L}\bigg{(}\mathcal{S}\big{(}\bm{x}-\bm{X}(s^{\prime},t)\big{)}+\frac{\epsilon^{2}r^{2}(s^{\prime})}{2}\mathcal{D}\big{(}\bm{x}-\bm{X}(s^{\prime},t)\big{)}\bigg{)}\bm{f}(s^{\prime},t)\,ds^{\prime}$
(9) $\displaystyle\mathcal{S}(\bm{x})$ $\displaystyle=\frac{{\bf
I}}{\left\lvert\bm{x}\right\rvert}+\frac{\bm{x}\bm{x}^{\rm
T}}{\left\lvert\bm{x}\right\rvert^{3}},\;\mathcal{D}(\bm{x})=\frac{{\bf
I}}{\left\lvert\bm{x}\right\rvert^{3}}-\frac{3\bm{x}\bm{x}^{\rm
T}}{\left\lvert\bm{x}\right\rvert^{5}},$

where $\bm{u}_{0}(\bm{x},t)$ is the fluid velocity in the absence of the fiber
and $\mu$ is the fluid viscosity. The expression
$\frac{1}{8\pi\mu}\mathcal{S}(\bm{x})$ is the free space Green’s function for
the Stokes equations in $\mathbb{R}^{3}$, commonly known as the Stokeslet,
while
$\frac{1}{8\pi\mu}\mathcal{D}(\bm{x})=\frac{1}{16\pi\mu}\Delta\mathcal{S}(\bm{x})$
is a higher order correction to the velocity approximation, often known as a
doublet. The force-per-unit-length $\bm{f}(s,t)$ exerted by the fluid on the
body is distributed between the generalized foci of the slender body at $s=\pm
L$.

In the stationary setting, the velocity field given by (9) is an
asymptotically accurate approximation to the velocity field around a three-
dimensional semi-flexible rod satisfying a well-posed slender body PDE,
defined in [24, 25] as the following boundary value problem for the Stokes
equations:

$\displaystyle-\mu\Delta\bm{u}+\nabla p$
$\displaystyle=0,\quad{\rm{div}\,}\,\bm{u}=0\qquad\text{ in
}\mathbb{R}^{3}\backslash\overline{\Sigma_{\epsilon}}$ (10)
$\displaystyle\int_{0}^{2\pi}(\bm{\sigma}\bm{n})\big{|}_{(\varphi(s),\theta)}\mathcal{J}_{\epsilon}(\varphi(s),\theta)\varphi^{\prime}(s)\,d\theta$
$\displaystyle=-\bm{f}(s)\hskip 59.75095pt\text{ on
}\partial\Sigma_{\epsilon}$
$\displaystyle\bm{u}\big{|}_{\partial\Sigma_{\epsilon}}$
$\displaystyle=\bm{u}(s),\hskip 62.59596pt\text{ unknown but independent of
}\theta$ $\displaystyle\left\lvert\bm{u}\right\rvert\to 0$
$\displaystyle\text{ as }\left\lvert\bm{x}\right\rvert\to\infty.$

Here $\bm{\sigma}=\mu\big{(}\nabla\bm{u}+(\nabla\bm{u})^{\rm T}\big{)}-p{\bf
I}$ is the fluid stress tensor, $\bm{n}(\bm{x})$ denotes the unit normal
vector pointing into $\Sigma_{\epsilon}$ at
$\bm{x}\in\partial\Sigma_{\epsilon}$, $\mathcal{J}_{\epsilon}(s,\theta)$ is
the Jacobian factor on $\partial\Sigma_{\epsilon}$, and
$\varphi(s):=\frac{s\sqrt{L^{2}+\epsilon^{2}}}{L}$ is a stretch function to
address the discrepancy between the extent of $\bm{f}$ and the extent of the
actual slender body surface. Given a force density $\bm{f}\in C^{1}(-L,L)$
which decays like $r(s)$ at the fiber endpoints ($\bm{f}(s)\sim r(\varphi(s))$
as $s\to\pm L$), the difference between the slender body approximation
$\bm{u}^{{\rm SB}}$ and the solution of (10) is bounded by an expression
proportional to $\epsilon\left\lvert\log\epsilon\right\rvert$. Note that
$r(s)$ need not be spheroidal (2) for this error analysis to hold, but $r(s)$
must decay smoothly to zero at the physical endpoints of the fiber at
$s=\pm\sqrt{L^{2}+\epsilon^{2}}$.

A key component of the well-posedness theory for the slender body PDE to which
(9) is an approximation is the fiber integrity condition on
$\bm{u}\big{|}_{\partial\Sigma_{\epsilon}}$. The fiber integrity condition
requires the velocity across each cross section $s$ of the slender body to be
constant; i.e. the velocity $\bm{u}(\bm{x})$ at any point
$\bm{x}(s,\theta)=\bm{X}(s)+\epsilon
r(s)\bm{e}_{r}(s,\theta)\in\partial\Sigma_{\epsilon}$ satisfies
$\partial_{\theta}\bm{u}(\bm{x}(s,\theta))=0$. This is to ensure that the
cross sectional shape of the fiber does not deform over time. An important
aspect of the accuracy of slender body theory is that the expression (9)
satisfies this fiber integrity condition to leading order in $\epsilon$.
Specifically, by Propositions 3.9 and 3.11 in [24, 25], respectively, the
angular dependence in $\bm{u}^{\rm SB}(\bm{x})$ over each cross section $s$ of
the slender body is only $\mathcal{O}(\epsilon\log\epsilon)$.

Another important general feature of the slender body PDE (10) is that the
operator mapping the force data $\bm{f}(s)$ to the $\theta$-independent fiber
velocity $\bm{u}|_{\partial\Sigma_{\epsilon}}(s)$ is negative definite (see
[22]; note that the sign convention for $\bm{f}$ is opposite).

Now, the velocity expression (9) is singular at $\bm{x}=\bm{X}(s,t)$ and can
be used only away from the fiber centerline; however, (9) presents a starting
point for approximating the velocity of the slender body itself. Various
methods can be used to obtain an expression for the relative velocity of the
fiber centerline $\frac{\partial\bm{X}(s,t)}{\partial t}$ which depends only
on the arclength parameter $s$ and time $t$ [7, 10, 15, 16, 19, 21, 28, 37].
Here we consider a different approach to deriving a limiting centerline
expression from (9) which evidently results in a negative definite integral
operator mapping $\bm{f}$ to $\bm{u}|_{\partial\Sigma_{\epsilon}}$. We then
regularize this first-kind integral equation in an asymptotically consistent
way to yield the second-kind integral equation (4). We outline our approach
here and provide a more detailed justification in Section 3.1.

The first step in approximating $\frac{\partial\bm{X}(s,t)}{\partial t}$ is to
evaluate (9) on the surface of the slender body at
$\bm{x}=\bm{X}(s,t)+\epsilon r(s)\bm{e}_{r}(s,\theta,t)$. Written out, the
velocity field along the fiber surface is given by

$\displaystyle 8\pi\mu\bigg{(}\bm{u}^{\rm
SB}(\bm{x}(s,\theta,t),t)-\bm{u}_{0}(\bm{X}(s,t),t)\bigg{)}=$ (11)
$\displaystyle\hskip 28.45274pt-\int_{-L}^{L}\bigg{(}\frac{{\bf
I}}{\left\lvert\bm{R}\right\rvert}+\frac{\overline{\bm{X}}\overline{\bm{X}}^{\rm
T}+\epsilon r(\overline{\bm{X}}\bm{e}_{r}^{\rm
T}+\bm{e}_{r}\overline{\bm{X}}^{\rm
T})+\epsilon^{2}r^{2}\bm{e}_{r}\bm{e}_{r}^{\rm
T}}{\left\lvert\bm{R}\right\rvert^{3}}$ $\displaystyle\hskip
56.9055pt+\frac{\epsilon^{2}r^{2}(s^{\prime})}{2}\bigg{(}\frac{{\bf
I}}{\left\lvert\bm{R}\right\rvert^{3}}-3\frac{\overline{\bm{X}}\overline{\bm{X}}^{\rm
T}+\epsilon r(\overline{\bm{X}}\bm{e}_{r}^{\rm
T}+\bm{e}_{r}\overline{\bm{X}}^{\rm
T})+\epsilon^{2}r^{2}\bm{e}_{r}\bm{e}_{r}^{\rm
T}}{\left\lvert\bm{R}\right\rvert^{5}}\bigg{)}\bigg{)}\bm{f}(s^{\prime},t)\,ds^{\prime},$

where unless otherwise specified, we have $r=r(s)$,
$\overline{\bm{X}}=\overline{\bm{X}}(s,s^{\prime},t)=\bm{X}(s,t)-\bm{X}(s^{\prime},t)$
and $\bm{R}=\bm{R}(s,s^{\prime},\theta,t)=\overline{\bm{X}}+\epsilon
r(s)\bm{e}_{r}(s,\theta,t)$. Now, along the fiber surface, the expression (11)
satisfies the fiber integrity condition to leading order in $\epsilon$; i.e.
the terms containing $\bm{e}_{r}(s,\theta,t)$ in (11) vanish to
$\mathcal{O}(\epsilon\log\epsilon)$. In particular, both the Stokeslet and
doublet include a $\theta$-dependent term with
$\epsilon^{2}r^{2}\bm{e}_{r}\bm{e}_{r}^{\rm T}$ in the numerator. Due to the
form of $\bm{R}$ in the denominator, both of these terms are $\mathcal{O}(1)$
at $s=s^{\prime}$; however, upon integrating in $s^{\prime}$, these terms
cancel each other asymptotically to order $\epsilon\log\epsilon$ (see
estimates 3.62 and 3.65 in [24] and estimates 3.40 and 3.43 in [25]).
Furthermore, the terms $\epsilon r(\overline{\bm{X}}\bm{e}_{r}^{\rm
T}+\bm{e}_{r}\overline{\bm{X}}^{\rm T})$ in both the Stokeslet and doublet
approximately integrate to zero in $s^{\prime}$, while the $\bm{e}_{r}$ term
in each denominator from
$\left\lvert\bm{R}(s,\theta,t)\right\rvert^{2}=\left\lvert\overline{\bm{X}}\right\rvert^{2}+2\epsilon
r\bm{e}_{r}\cdot\overline{\bm{X}}+\epsilon^{2}r^{2}$ is also only
$O(\epsilon\log\epsilon)$ (see Propositions 3.9 and 3.11 in [24], [25],
respectively).

Due to these cancellations and the fact that dropping these terms still
approximates the slender body PDE solution of [24, 25] to at least
$O(\epsilon\log\epsilon)$, we may eliminate all terms containing
$\bm{e}_{r}(s,\theta,t)$ in (11) to obtain a $\theta$-independent expression
which approximates the velocity of the fiber itself:

$\displaystyle 8\pi\mu\bigg{(}\frac{\partial\bm{X}}{\partial
t}-\bm{u}_{0}(\bm{X}(s,t),t)\bigg{)}=-\int_{-L}^{L}\bigg{(}\frac{{\bf
I}}{(|\overline{\bm{X}}|^{2}+\epsilon^{2}r^{2}(s))^{1/2}}+\frac{\overline{\bm{X}}\overline{\bm{X}}^{\rm
T}}{(|\overline{\bm{X}}|^{2}+\epsilon^{2}r^{2}(s))^{3/2}}$ (12)
$\displaystyle\hskip
113.81102pt+\frac{\epsilon^{2}r^{2}(s^{\prime})}{2}\bigg{(}\frac{{\bf
I}}{(|\overline{\bm{X}}|^{2}+\epsilon^{2}r^{2}(s))^{3/2}}-\frac{3\overline{\bm{X}}\overline{\bm{X}}^{\rm
T}}{(|\overline{\bm{X}}|^{2}+\epsilon^{2}r^{2}(s))^{5/2}}\bigg{)}\bigg{)}\bm{f}(s^{\prime},t)\,ds^{\prime}.$

The expression (12) serves as the model underlying our final slender body
velocity expression (4). In Section 3.1, we show that in a simplified setting,
(12) results in a negative definite operator mapping the force density
$\bm{f}$ to the fiber velocity $\frac{\partial\bm{X}}{\partial t}$, whereas
other models which rely on further asymptotic expansion of (11) about
$\epsilon=0$ do not, and incur high wavenumber instabilities. This phenomenon
is well known for the Keller–Rubinow model [16, 10], but for other possible
centerline expressions, including models similar to Lighthill [19], this high
wavenumber instability has not been documented previously. It seems that our
model (12) may be the simplest that can be obtained by expanding from (11)
while still guaranteeing a negative definite operator.

Now, since the integral operator in (12) has a smooth kernel, the expression
(12) yields a first-kind Fredholm integral equation for $\bm{f}$ when the
fiber velocity $\frac{\partial\bm{X}}{\partial t}$ is supplied. Describing the
motion of a rigid fiber involves inverting this expression to solve for
$\bm{f}$, which in general is an ill-posed problem for a first-kind equation.
Thus we want to regularize the integral operator (12) to create a second-kind
integral equation while keeping the same order of accuracy in the map
$\bm{f}\mapsto\frac{\partial\bm{X}}{\partial t}$.

We first note that, for $\eta>1$, we have the following identity:

$\int_{-L}^{L}\bigg{(}\frac{1}{(|\overline{\bm{X}}|^{2}+\epsilon^{2}r^{2}(s))^{1/2}}-\frac{1}{(|\overline{\bm{X}}|^{2}+\eta^{2}\epsilon^{2}r^{2}(s))^{1/2}}\bigg{)}g(s^{\prime})\,ds^{\prime}=2\log\eta\,g(s)+\mathcal{O}(\eta\epsilon\log(\eta\epsilon)).$
(13)

###### Proof.

By Lemma 3.8 in [25], for $a>0$ sufficiently small, we have

$\displaystyle\int_{-L}^{L}$
$\displaystyle\bigg{(}\frac{g(s^{\prime})}{(|\overline{\bm{X}}|^{2}+a^{2}r^{2}(s))^{1/2}}-\frac{g(s^{\prime})}{|\overline{\bm{X}}|}+\frac{g(s)}{\left\lvert
s-s^{\prime}\right\rvert}\bigg{)}\,ds^{\prime}$ (14) $\displaystyle\hskip
113.81102pt=\log\bigg{(}\frac{2(L^{2}-s^{2})+2\sqrt{(L^{2}-s^{2})^{2}+a^{2}r^{2}(s)}}{a^{2}r^{2}(s)}\bigg{)}+\mathcal{O}(a\log
a).$

Subtracting (14) with $a=\eta\epsilon$ from (14) with $a=\epsilon$ and using
that

$\displaystyle\left\lvert\log\bigg{(}\frac{(L^{2}-s^{2})+\sqrt{L^{2}+\epsilon^{2}r^{2}}}{(L^{2}-s^{2})+\sqrt{L^{2}+\eta^{2}\epsilon^{2}r^{2}}}\bigg{)}\right\rvert=\left\lvert\log\bigg{(}\frac{(L^{2}-s^{2})+\sqrt{L^{2}+\epsilon^{2}r^{2}}}{(L^{2}-s^{2})+\sqrt{L^{2}+\eta^{2}\epsilon^{2}r^{2}}}\bigg{)}-\log(1)\right\rvert\leq
C\epsilon^{2},$

we obtain (13). ∎

Using (13), we replace the first term in the integrand of (12) to obtain (4).
We can compare the expression (4) to that of Tornberg and Shelley in [37],
where a regularization of the Keller–Rubinow model is used to obtain a second-
kind integral equation for $\bm{f}$. One thing to note is that, due to the
form of the local term in our model (4), the effect of the regularization
parameter $\eta$ is the same in all directions (both tangent and normal to the
fiber centerline). This is not necessarily the case for the Tornberg and
Shelley model (see Section 3.1.3 for a spectral comparison given a simplified
fiber geometry).

###  3.1 Spectral comparison of slender body integral operators

In this subsection we provide evidence that our model (4) is well suited for
approximating the map $\frac{\partial\bm{X}}{\partial t}\mapsto\bm{f}$ needed
to simulate the motion of a rigid fiber. Here we consider the spectrum of the
integral operator taking the force density $\bm{f}$ to the fiber velocity
$\frac{\partial\bm{X}}{\partial t}$ in the non-physical but nevertheless
instructive case of a straight, periodic fiber with constant radius
$\epsilon$. In this scenario we can explicitly calculate the eigenvalues of
both the slender body PDE operator (10) as well as the integral operator (12)
and related models. This allows us to directly compare the properties of
different models in the same simple setting and serves as a starting point for
understanding more complicated geometries. In particular, we expect this
analysis to roughly capture the high wavenumber behavior of these models in
different geometries – on length scales much smaller than the variation in
curvature and fiber radius. The high wavenumber behavior is of particular
interest for the invertibility and stability of the slender body theory
integral operator.

For comparison, we first recall the form of the eigenvalues of the slender
body PDE (10), calculated in [22]. In Section 3.1.2, we consider the model
(12), before regularization, and show that the integral operator is negative
definite. We compare the spectrum of (12) to three other possible models based
on slender body theory which do not result in negative definite operators.
Then in Section 3.1.3, we consider the regularized version of our model (4)
and compare its spectrum to the regularized model of Tornberg and Shelley
[37]. We note that in our model, a uniform regularization parameter appears to
give the best approximation of the slender body PDE spectrum in directions
both normal and tangent to the slender body centerline, whereas in the
Tornberg–Shelley model, the parameter required by the tangential direction may
not be optimal in the normal direction.

####  3.1.1 Spectrum of the slender body PDE

Here we consider a straight, periodic fiber with constant radius $\epsilon$.
We take the fiber centerline to be 2-periodic and lie along the $z$-axis,
$\bm{X}(z)=z\bm{e}_{z}$, $z\in\mathbb{R}/2\mathbb{Z}$, and for simplicity take
$\mu=1$ and zero background flow. We consider the stationary setting and omit
the time dependence in our notation; in particular, we denote the fiber
velocity by $\overline{\bm{u}}(z)$ to distinguish from the fluid velocity away
from the fiber.

We consider this scenario because we can explicitly calculate the eigenvalues
of the slender body PDE (10) as well as various possible integral expressions
for approximating the map $\bm{f}\mapsto\overline{\bm{u}}$. In particular, the
eigenvectors of this map can be decomposed into tangential ($\bm{e}_{z}$) and
normal ($\bm{e}_{x},\bm{e}_{y}$) directions and are given by
$\bm{f}_{m}(z)=e^{i\pi kz}\bm{e}_{m}$, $m=x,y,z$. We may then explicitly solve
for $\lambda^{m}_{k}$ satisfying

$\overline{\bm{u}}(z)=\lambda^{m}_{k}\bm{f}_{m}(z),\qquad m=x,y,z$ (15)

for both the slender body PDE operator and various approximations based on
slender body theory. To avoid logarithmic growth of the corresponding bulk
velocity field at spatial infinity, we will ignore translational modes ($k=0$)
in the following spectral analysis. Clearly these modes are important,
especially for a rigid body; however, we are mainly interested in the high
wavenumber behavior of these operators. High wavenumber instabilities are a
known issue for nonlocal slender body theory [10, 31, 37], and the following
analysis likely captures the behavior of these models at high wavenumbers
(small length scales) even in curved geometries.

To begin, the eigenvalues of the slender body PDE operator (10) mapping
$\bm{f}$ to $\overline{\bm{u}}$ were calculated in [22], Proposition 1.4. Note
that the sign convention in this paper is opposite, as we are considering
$\bm{f}$ to be the hydrodynamic force exerted _by_ rather than _on_ the
slender body. For the slender body PDE, the eigenvalues satisfying (15) in the
tangential and normal directions, respectively, are given by

$\lambda^{m}_{k}=\begin{cases}-\frac{2K_{0}K_{1}+\pi\epsilon\left\lvert
k\right\rvert\big{(}K_{0}^{2}-K_{1}^{2}\big{)}}{4\pi^{2}\epsilon\left\lvert
k\right\rvert K_{1}^{2}},&m=z\\\
-\frac{2K_{0}K_{1}K_{2}+\pi\epsilon\left\lvert
k\right\rvert\big{(}K_{1}^{2}(K_{0}+K_{2})-2K_{0}^{2}K_{2}\big{)}}{2\pi^{2}\epsilon\left\lvert
k\right\rvert\big{(}4K_{1}^{2}K_{2}+\pi\epsilon\left\lvert k\right\rvert
K_{1}(K_{1}^{2}-K_{0}K_{2})\big{)}},&m=x,y\end{cases}$ (16)

where each $K_{j}=K_{j}(\pi\epsilon\left\lvert k\right\rvert)$, $j=0,1,2$, is
a $j^{\text{th}}$ order modified Bessel function of the second kind. Note that
both sets of eigenvalues $\lambda^{z}_{k}$ and
$\lambda^{x}_{k},\lambda^{y}_{k}$ are strictly negative and decay to 0 at a
rate proportional to $1/\left\lvert k\right\rvert$ as $\left\lvert
k\right\rvert\to\infty$. We will compare our approximation and various other
slender body approximations to (16).

####  3.1.2 Pre-regularization comparison

Before we consider the regularized version (4) of our model, we consider the
base model (12) and compare its spectrum to other existing models based on
slender body theory, before regularization. In the straight-but-periodic
scenario, our model (12) becomes the periodization of the expression

$\overline{\bm{u}}(z)=-\frac{1}{8\pi}\int_{-1}^{1}\bigg{(}\frac{\bf
I}{(\overline{z}^{2}+\epsilon^{2})^{1/2}}+\frac{\overline{z}^{2}\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{3/2}}+\frac{\epsilon^{2}}{2}\bigg{(}\frac{\bf
I}{(\overline{z}^{2}+\epsilon^{2})^{3/2}}-3\frac{\overline{z}^{2}\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{5/2}}\bigg{)}\bigg{)}\bm{f}(z-\overline{z})\,d\overline{z}.$
(17)

For this geometry, we may calculate the eigenvalues $\lambda^{m}_{k}$
satisfying (15), which are given by

$\displaystyle\lambda^{m}_{k}$
$\displaystyle=\begin{cases}\displaystyle-\frac{1}{8\pi}\int_{-1}^{1}\frac{2\overline{z}^{4}+2\epsilon^{2}\overline{z}^{2}+\frac{3}{2}\epsilon^{4}}{(\overline{z}^{2}+\epsilon^{2})^{5/2}}e^{-i\pi
k\overline{z}}\,d\overline{z},&m=z\\\
\displaystyle-\frac{1}{8\pi}\int_{-1}^{1}\frac{\overline{z}^{2}+\frac{3}{2}\epsilon^{2}}{(\overline{z}^{2}+\epsilon^{2})^{3/2}}e^{-i\pi
k\overline{z}}\,d\overline{z},&m=x,y.\end{cases}$ (18)

These integrals may be computed explicitly to obtain

$\lambda^{m}_{k}=\begin{cases}\displaystyle-\frac{1}{8\pi}\bigg{(}(4+\pi^{2}\epsilon^{2}k^{2})K_{0}(\pi\epsilon\left\lvert
k\right\rvert)-2\pi\epsilon\left\lvert k\right\rvert
K_{1}(\pi\epsilon\left\lvert k\right\rvert)\bigg{)},&m=z\\\
\displaystyle-\frac{1}{8\pi}\bigg{(}2K_{0}(\pi\epsilon\left\lvert
k\right\rvert)+\pi\epsilon\left\lvert k\right\rvert
K_{1}(\pi\epsilon\left\lvert k\right\rvert)\bigg{)},&m=x,y.\end{cases}$ (19)

Here $K_{0}$ and $K_{1}$ are zero and first order modified Bessel functions of
the second kind, respectively. The eigenvalues $\lambda^{m}_{k}$ lie along the
curves plotted in Figure 1. Importantly, these eigenvalues satisfy the
following lemma.

######  Lemma 3.1.

For all $\left\lvert k\right\rvert\geq 1$ and $m=x,y,z$, the eigenvalues
$\lambda^{m}_{k}$ given by (19) satisfy $\lambda^{m}_{k}<0$.

###### Proof.

The case $m=x,y$ is immediate, since $K_{0}(t)>0$ and $K_{1}(t)>0$ for any
$t>0$.

For the tangential direction $m=z$, we first note that, by Lemma 1.16 in [22],
we have

$1\leq\frac{K_{1}(t)}{K_{0}(t)}\leq 1+\frac{1}{2t}$

for all $t>0$. Letting $g(t)=(4+t^{2})K_{0}(t)-2tK_{1}(t)$, it suffices to
show that $g(t)/K_{0}(t)>0$. But

$\frac{g(t)}{K_{0}(t)}=4+t^{2}-2t\frac{K_{1}(t)}{K_{0}(t)}\geq
3+t^{2}-2t>(t-\sqrt{3})^{2}\geq 0.$

∎

Now, at a continuous level, regularization is necessary to make sense of
inverting the integral operator (17), since $K_{0}$ and $K_{1}$ decay
exponentially as $\left\lvert k\right\rvert\to\infty$. However, at a discrete
level, numerical approximation of (17) will be invertible, albeit with a large
condition number, due to Lemma 3.1. This negativity does not hold for other
popular slender body approximations which rely on further asymptotic expansion
of (12) with respect to $\epsilon$ to obtain a limiting centerline velocity
expression. In particular, we consider the models of Keller and Rubinow [16]
and of Lighthill [19].

The Keller–Rubinow model, proposed in [16] and further studied by [10, 15, 31,
37], is equivalent to a full matched asymptotic expansion of (11) about
$\epsilon=0$. In the straight-but-periodic setting, the Keller–Rubinow
expression for the slender body velocity is given by

$8\pi\overline{\bm{u}}(z)=-\bigg{(}({\bf I}-3\bm{e}_{z}\bm{e}_{z}^{\rm
T})-2\log(\pi\epsilon/8)({\bf I}+\bm{e}_{z}\bm{e}_{z}^{\rm
T})\bigg{)}\bm{f}(z)-({\bf I}+\bm{e}_{z}\bm{e}_{z}^{\rm
T})\frac{\pi}{2}\int_{-1}^{1}\frac{\bm{f}(z-\overline{z})-\bm{f}(z)}{\left\lvert\sin(\pi\overline{z}/2)\right\rvert}\,d\overline{z}.$
(20)

The eigenvalues of the periodic Keller–Rubinow operator taking $\bm{f}$ to
$\overline{\bm{u}}$ have been calculated in [10, 31, 37] and are given by

$\lambda^{m}_{k}=\begin{cases}\displaystyle\frac{1}{4\pi}\big{(}1+2\log(\pi\epsilon\left\lvert
k\right\rvert/2)+2\gamma\big{)},&m=z\\\
\displaystyle-\frac{1}{8\pi}\big{(}1-2\log(\pi\epsilon\left\lvert
k\right\rvert/2)-2\gamma\big{)},&m=x,y.\end{cases}$ (21)

Here $\gamma\approx 0.5772$ is the Euler gamma.

In both the tangent and normal directions, however, the Keller–Rubinow
approximation runs into stability issues at moderately high wavenumbers,
apparent in Figure 1 at $\left\lvert
k\right\rvert=\frac{2e^{-\gamma-1/2}}{\pi\epsilon}\approx 0.217/\epsilon$
(tangent) and $\left\lvert
k\right\rvert=\frac{2e^{-\gamma+1/2}}{\pi\epsilon}\approx 0.589/\epsilon$
(normal). In particular, the curve containing the eigenvalues
$\lambda^{m}_{k}$ crosses zero and becomes negative. This is an issue both
because the slender body PDE eigenvalues (16) are strictly negative, and
because, for arbitrary $\epsilon$, there is no clear way to guarantee that
$\lambda^{m}_{k}\neq 0$, especially for more complicated fiber geometries.
Thus some sort of regularization of (20) is necessary before approximating the
inverse map $\overline{\bm{u}}\mapsto\bm{f}$.

Figure 1: The eigenvalues $\lambda^{m}_{k}$ of the operator mapping
$\bm{f}\mapsto\overline{\bm{u}}$ in various slender body models lie along the
curves plotted here in the case of a straight-but-periodic fiber. Our model
(blue) results in strictly negative eigenvalues in both the tangential and
normal directions, as does the slender body PDE (dotted). The Keller–Rubinow
approximation (green) exhibits instabilities at wavenumbers $\left\lvert
k\right\rvert\approx 0.2/\epsilon$ (tangential direction) and $\left\lvert
k\right\rvert\approx 0.6/\epsilon$ (normal direction) as the eigenvalues of
the operator mapping $\bm{f}\mapsto\overline{\bm{u}}$ become negative. For the
modified Lighthill models, the normal direction eigenvalues $\lambda^{x}_{k}$
and $\lambda^{y}_{k}$ (red) remain negative at high wavenumber, but in the
tangential direction, the eigenvalues of Modified Lighthill 1 (red) become
negative when $\left\lvert k\right\rvert>0.5/\epsilon$. Furthermore, the
tangential eigenvalues of Modified Lighthill 2 (magenta) do not agree with the
slender body PDE at low wavenumber.

In addition to the Keller–Rubinow model, we consider what we will term the
_modified Lighthill_ approach to deriving a fiber velocity approximation. This
approach takes advantage of the fact that the doublet term of (12) only has an
$O(1)$ contribution to the fiber velocity very close to $s^{\prime}=s$, and
thus can be integrated asymptotically to leave only a local term. This results
in a model similar to that of Lighthill [19], which was derived via different
reasoning but also includes a local doublet term and a nonlocal Stokeslet
contribution (see Remark 3.2).

There are two ways to consider the nonlocal Stokeslet contribution. The first
expression, which we will term Modified Lighthill 1, is given by the
periodization of

$\overline{\bm{u}}(z)=-\frac{1}{8\pi}\bigg{(}({\bf
I}-\bm{e}_{z}\bm{e}_{z}^{\rm T})\bm{f}(s)+\int_{-1}^{1}\bigg{(}\frac{\bf
I}{(\overline{z}^{2}+\epsilon^{2})^{1/2}}+\frac{\overline{z}^{2}\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{3/2}}\bigg{)}\bm{f}(z-\overline{z})\,d\overline{z}\bigg{)}.$
(22)

Here the local term $({\bf I}-\bm{e}_{z}\bm{e}_{z}^{\rm T})$ comes from
asymptotically integrating the doublet term of (11) (see estimate 3.65 of [24]
for more detail). Note that in (22), the Stokeslet term inside the integral is
equal to $\bm{f}/\epsilon$ when $\overline{z}=0$.

For the second expression, which we will call Modified Lighthill 2, the
$\bm{e}_{z}\bm{e}_{z}^{\rm T}$ component of the Stokeslet term is normalized
to give the same order contribution at $\overline{z}=0$ as in (11); namely,
$({\bf I}+\bm{e}_{z}\bm{e}_{z}^{\rm T})\bm{f}/\epsilon$. This yields the
periodization of the expression

$\overline{\bm{u}}(z)=-\frac{1}{8\pi}\bigg{(}({\bf
I}-\bm{e}_{z}\bm{e}_{z}^{\rm T})\bm{f}(s)+\int_{-1}^{1}\frac{{\bf
I}+\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{1/2}}\bm{f}(z-\overline{z})\,d\overline{z}\bigg{)}.$
(23)

######  Remark 3.2.

The actual model proposed by Lighthill in [19], written in the periodic,
straight setting, has the form

$\overline{\bm{u}}(z)=-\frac{1}{8\pi}\bigg{(}2({\bf
I}-\bm{e}_{z}\bm{e}_{z}^{\rm
T})\bm{f}(z)+\int_{\left\lvert\overline{z}\right\rvert>q}\frac{{\bf
I}+\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{\left\lvert\overline{z}\right\rvert}\bm{f}(z-\overline{z})\,d\overline{z}\bigg{)};\quad
q=\epsilon\sqrt{e}/2.$ (24)

At first glance, this looks like a slightly different model from (22) and
(23), due to the 2 in front of the $({\bf I}-\bm{e}_{z}\bm{e}_{z}^{\rm
T})\bm{f}(z)$ term. However, the extra factor here is precisely due to the
removal of the section $\left\lvert\overline{z}\right\rvert\leq q$ from the
integral term. Indeed, if we consider the integrand of (22), we note that

$\displaystyle\int_{-q}^{q}\bigg{(}\frac{\bf
I}{(\overline{z}^{2}+\epsilon^{2})^{1/2}}+\frac{\overline{z}^{2}\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{3/2}}\bigg{)}\bm{f}(z-\overline{z})\,d\overline{z}$
$\displaystyle=\big{(}2\log(2q/\epsilon)({\bf I}+\bm{e}_{z}\bm{e}_{z}^{\rm
T})-2\bm{e}_{z}\bm{e}_{z}^{\rm T}\big{)}\bm{f}(z)+O(\epsilon^{2}/q^{2})$
$\displaystyle=({\bf I}-\bm{e}_{z}\bm{e}_{z}^{\rm
T})\bm{f}(z)+O(\epsilon^{2}/q^{2})$

for $q$ as in (24). Now, this particular choice of $q$ is not large relative
to $\epsilon$, so the $O(\epsilon^{2}/q^{2})$ error term is not small
asymptotically. However, this is merely a heuristic and we will not be
considering the expression (24) in greater depth here. Furthermore, the
expressions (22) and (23) are more amenable to calculating eigenvalues.

The eigenvalues of (22) are given by

$\displaystyle\lambda^{m}_{k}$
$\displaystyle=\begin{cases}\displaystyle-\frac{1}{4\pi}\bigg{(}2K_{0}(\pi\epsilon\left\lvert
k\right\rvert)-\pi\epsilon\left\lvert k\right\rvert
K_{1}(\pi\epsilon\left\lvert k\right\rvert)\bigg{)},&m=z\\\
\displaystyle-\frac{1}{8\pi}\bigg{(}1+2K_{0}(\pi\epsilon\left\lvert
k\right\rvert)\bigg{)},&m=x,y.\end{cases}$ (25)

Now the normal eigenvalues $\lambda^{x}_{k}$ and $\lambda^{y}_{k}$ are always
negative. However, there is still a high wavenumber instability in the tangent
direction. In particular, $\lambda^{z}_{k}=0$ when $\pi\epsilon\left\lvert
k\right\rvert\approx 1.55265$, and becomes positive at higher wavenumbers (see
Figure 1). Thus the instability issue is not fully resolved by expanding only
the doublet term of (11).

For Modified Lighthill 2, the eigenvalues of (23) are given by

$\displaystyle\lambda^{m}_{k}$
$\displaystyle=\begin{cases}\displaystyle-\frac{1}{2\pi}K_{0}(\pi\epsilon\left\lvert
k\right\rvert),&m=z\\\
\displaystyle-\frac{1}{8\pi}\bigg{(}1+2K_{0}(\pi\epsilon\left\lvert
k\right\rvert)\bigg{)},&m=x,y.\end{cases}$ (26)

Here the eigenvalues $\lambda^{x}_{k}$ and $\lambda^{y}_{k}$ in the normal
directions are identical to (25), but the tangential eigenvalues
$\lambda^{z}_{k}$ are very different. In fact, they are too different: Recall
that near $t=0$, the modified Bessel functions $K_{0}(t)$ and $K_{1}(t)$
satisfy

$K_{0}(t)=-\log(t/2)-\gamma+O(t^{2});\quad tK_{1}(t)=1+O(t^{2}).$ (27)

Therefore, at low wavenumber ($k=O(1)$), the tangential eigenvalues of
Modified Lighthill 2 (23) look like

$\lambda^{z}_{k}=\frac{1}{2\pi}(\log(\pi\epsilon\left\lvert
k\right\rvert/2)+\gamma)+O(\epsilon^{2}k^{2}).$

This does not agree with the low wavenumber behavior of the slender body PDE
(16) (see Figure 1). It appears that the normalization in Modified Lighthill 2
(23) results in the wrong model.

For the sake of completeness, we also consider a modification of our model
(12) in which the $\overline{\bm{X}}\overline{\bm{X}}^{\rm T}$ terms are
normalized as in Modified Lighthill 2 (23) to yield a nonzero contribution to
the fiber velocity when $s=s^{\prime}$. In the case of the periodic straight
centerline, the modified version of our model becomes the periodization of

$\overline{\bm{u}}(z)=-\frac{1}{8\pi}\int_{-1}^{1}\bigg{(}\frac{{\bf
I}+\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{1/2}}+\frac{\epsilon^{2}}{2}\frac{{\bf
I}-3\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{3/2}}\bigg{)}\bm{f}(z-\overline{z})\,d\overline{z}.$
(28)

The eigenvalues of (28) are given by

$\displaystyle\lambda^{m}_{k}$
$\displaystyle=\begin{cases}\displaystyle-\frac{1}{4\pi}\bigg{(}2K_{0}(\pi\epsilon\left\lvert
k\right\rvert)-\pi\epsilon\left\lvert k\right\rvert
K_{1}(\pi\epsilon\left\lvert k\right\rvert)\bigg{)},&m=z\\\
\displaystyle-\frac{1}{8\pi}\bigg{(}2K_{0}(\pi\epsilon\left\lvert
k\right\rvert)+\pi\epsilon\left\lvert k\right\rvert
K_{1}(\pi\epsilon\left\lvert k\right\rvert)\bigg{)},&m=x,y.\end{cases}$ (29)

Now, the eigenvalues $\lambda^{x}_{k}$ and $\lambda^{y}_{k}$ in the directions
normal to the fiber are unchanged from our original expression (19). However,
the tangent eigenvalues $\lambda^{z}_{k}$ are now given by the same expression
as Modified Lighthill 1 (25), which we recall exhibits a high wavenumber
instability (Figure 1).

The takeaway here is that, at least in the case of a straight, periodic fiber,
our model (12), before regularization, captures the negative-definiteness of
the the slender body PDE and provides a better approximation than other
slender body models (20), (22), (23), (28).

####  3.1.3 Regularized comparison

To make our model truly suitable for inversion, we need to regularize the
integral kernel as in (4). In the straight-but-periodic setting, the operator
in (4) becomes the periodization of

$\displaystyle 8\pi\overline{\bm{u}}(z)$
$\displaystyle=-2\log\eta\,\bm{f}(z)-\int_{-1}^{1}\bigg{(}\frac{\bf
I}{(\overline{z}^{2}+\eta^{2}\epsilon^{2})^{1/2}}+\frac{\overline{z}^{2}\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{3/2}}$ (30) $\displaystyle\hskip
113.81102pt+\frac{\epsilon^{2}}{2}\bigg{(}\frac{\bf
I}{(\overline{z}^{2}+\epsilon^{2})^{3/2}}-3\frac{\overline{z}^{2}\bm{e}_{z}\bm{e}_{z}^{\rm
T}}{(\overline{z}^{2}+\epsilon^{2})^{5/2}}\bigg{)}\bigg{)}\bm{f}(z-\overline{z})\,d\overline{z}.$

The eigenvalues of (30) are then given by

$\lambda^{m}_{k}=\begin{cases}\displaystyle-\frac{1}{8\pi}\bigg{(}2\log\eta+2K_{0}(\eta\pi\epsilon\left\lvert
k\right\rvert)+(2+\pi^{2}\epsilon^{2}k^{2})K_{0}(\pi\epsilon\left\lvert
k\right\rvert)-2\pi\epsilon\left\lvert k\right\rvert
K_{1}(\pi\epsilon\left\lvert k\right\rvert)\bigg{)},&m=z\\\
\displaystyle-\frac{1}{8\pi}\bigg{(}2\log\eta+2K_{0}(\eta\pi\epsilon\left\lvert
k\right\rvert)+\pi\epsilon\left\lvert k\right\rvert
K_{1}(\pi\epsilon\left\lvert k\right\rvert)\bigg{)},&m=x,y.\end{cases}$ (31)

For $\eta>1$, the spectrum of our operator is bounded away from 0 and (30) is
a second-kind integral equation for $\bm{f}$.

We can compare the behavior of (30) with the Tornberg–Shelley regularization
of the Keller–Rubinow model. In [37, 31], the high wavenumber instability in
(20) is removed by replacing the denominator of the integral term, which
vanishes at $\overline{z}=0$, with an expression proportional to $\epsilon$ at
$\overline{z}=0$. Using the relation

$\int_{-1}^{1}\bigg{(}\frac{\pi}{\left\lvert 2\sin(\pi
z/2)\right\rvert}-\frac{1}{\left\lvert
z\right\rvert}\bigg{)}\,dz=-2\log(\pi/4),$

to rewrite (20), a regularization $\delta\epsilon$, $\delta>0$, is added to
the denominator to obtain

$8\pi\overline{\bm{u}}(z)=-\bigg{(}({\bf I}-3\bm{e}_{z}\bm{e}_{z}^{\rm
T})+2\log(\delta)({\bf I}+\bm{e}_{z}\bm{e}_{z}^{\rm T})\bigg{)}\bm{f}(z)-({\bf
I}+\bm{e}_{z}\bm{e}_{z}^{\rm
T})\int_{-1}^{1}\frac{\bm{f}(z-\overline{z})}{(\overline{z}^{2}+\delta^{2}\epsilon^{2})^{1/2}}\,d\overline{z}.$
(32)

Here we have also used that the second term in the original Keller–Rubinow
integral expression can now be integrated up to $O(\epsilon^{2})$ errors to
nearly cancel the logarithmic term in (20), leaving only $\log(\delta)$. The
idea is to then choose $\delta$ such that all eigenvalues of the operator
taking $\bm{f}\mapsto\overline{\bm{u}}$ are negative. Since the integral
kernel is now smooth, (32) is now a second-kind integral equation for
$\bm{f}$.

The eigenvalues of this $\delta$-regularized Keller–Rubinow operator are given
by

$\displaystyle\lambda^{m}_{k}$
$\displaystyle=\begin{cases}\displaystyle-\frac{1}{4\pi}\bigg{(}-1+2\log\delta+2K_{0}(\delta\pi\epsilon\left\lvert
k\right\rvert)\bigg{)},&m=z\\\
\displaystyle-\frac{1}{8\pi}\bigg{(}1+2\log\delta+2K_{0}(\delta\pi\epsilon\left\lvert
k\right\rvert)\bigg{)},&m=x,y.\end{cases}$ (33)

Since $K_{0}$ is positive, $\lambda^{z}_{k}$ is guaranteed to be negative and
bounded away from 0 as long as $\delta>\sqrt{e}$ (see Figure 2).

Figure 2: The eigenvalues of our regularized model (30) with $\eta=1.1$ and
the Tornberg–Shelley $\delta$-regularized model (32) with
$\delta=\sqrt{e}+0.1$ lie along the blue and red curves, respectively. Note
that the regularization parameter $\eta$ in our model affects the tangential
and normal eigenvalues in a similar way; in particular, $\eta>1$ is required
in both cases to ensure that (30) is a second-kind integral equation. In the
$\delta$-regularized model, the tangential direction requires
$\delta>\sqrt{e}$, but the normal direction does not, resulting at least
visually in a greater disparity between the $\lambda^{x}_{k},\lambda^{y}_{k}$
for the PDE (dotted) and the $\delta$-regularized approximation.

Note that in our model (30), the regularization parameter $\eta$ affects the
spectrum of the operator mapping $\bm{f}$ to $\overline{\bm{u}}$ in the same
way in both the tangential and normal directions. In particular, in both
directions, $\eta>1$ is required to obtain the desired second-kind integral
equation. In the Tornberg–Shelley model, the bound $\delta>\sqrt{e}\approx
1.649$ is required to ensure negativity of the tangential eigenvalues, but
this lower bound does not apply to the normal direction; in fact,
$\delta>e^{-1}\approx 0.368$ is sufficient for ensuring strictly negative
normal eigenvalues. This may mean that our model can achieve better agreement
with the slender body PDE in both the tangent and normal directions at the
same time.

In [22], it is shown that using the $\delta$-regularized model (32) to
approximate the map $\overline{\bm{u}}\mapsto\bm{f}$ yields $\epsilon^{2}$
convergence to the slender body PDE for sufficiently smooth
$\overline{\bm{u}}$. It is also shown that the constant in the resulting error
estimate has the form $C_{1}\delta^{2}(1+\log\delta)+C_{2}/(-1+\log\delta)$
for constants $C_{1}$ and $C_{2}$. We expect that a similar error estimate and
analogous $\eta$ dependence hold for our model (30); i.e. the constant should
look like $C_{1}\eta^{2}+C_{2}/\log\eta$. If $C_{1}\approx C_{2}$, this yields
an optimal $\eta$ of approximately $1.5$. This should give a rough guideline
for a good choice of $\eta$ for more general curved geometries, at least in
the periodic setting.

##  4 Numerical discretization of the slender body model

We turn now to a numerical method for simulating thin rigid fibers in flows.
We begin by discussing a general method for numerically solving Fredholm
integral equations where the solution must be integrated (i.e. to find the
total force and torque on a rigid fiber). We apply these general methods to
the slender body model (4) and perform convergence tests. We note improvements
in conditioning and stability for the second kind ($\eta>1$) versus first kind
($\eta=1$) integral equation. Finally, we look at the spectrum of the
discretized integral operator in different geometries to verify the negative
definite nature of the operator.

###  4.1 Solving the second-kind Fredholm integral equation

Denote by $\mathbf{K}:L^{2}([-L,L],\mathbb{R}^{3})\rightarrow
L^{2}([-L,L],\mathbb{R}^{3})$ the integral operator

$\mathbf{K}[\bm{f}](s):=\int_{-L}^{L}\\!K(s,s^{\prime})\bm{f}(s^{\prime})ds^{\prime}.$
(34)

Then a Fredholm integral equation of the first kind reads

$\bm{y}(s)=\mathbf{K}[\bm{f}](s).$ (35)

It is well known that the inversion of such an integral operator is an ill-
posed problem, meaning that the solution may not be unique or not even exist
[1, 13, 17]. Furthermore, small perturbations to the left hand side of (35)
can lead to relatively large perturbations of the solution $\bm{f}(s)$. The
ill-posedness of this problem can be circumvented by regularizing the integral
operator into a second-kind Fredholm integral equation, which takes the form

$\bm{y}(s)=(\alpha\mathbf{I}+\mathbf{K})[\bm{f}](s)$ (36)

for some parameter $\alpha$. Discretization of (36) yields a linear system
with a far better condition number. The connection between equation (36) and
our model is illustrated in Section 4.2.

Numerical methods for solving Fredholm integral equations are well documented
[13, 39] and the approach we adopt is based on the Nyström method [1, Chapt.
12.4]. The main additional consideration for rigid fibers is that after
numerically inverting a second-kind Fredholm integral equation, linear
functionals (8) need to be applied to the resulting $\bm{f}(s)$ to find the
total force and torque.

We consider now the numerical approximation of a general linear functional of
$\bm{f}(s)$, given by

${\phi}_{M}(\bm{f})=\int_{-L}^{L}\\!M(s)\bm{f}(s)\mathrm{d}s.$ (37)

Here $M(s)\in\mathbb{R}^{3\times 3}$ is a bounded, smooth operator and
$\bm{f}(s)$ is found by numerically inverting a second-kind Fredholm integral
equation of the form (36). The numerical method is obtained discretizing the
equation (36) by replacing the integral with a convergent quadrature formula
with nodes $-L=s_{1}<s_{2}<...<s_{n}=L$ and weights
$\bm{w}=(w_{1},w_{2},...,w_{n})^{T}\in\mathbb{R}^{n}$, and requiring the
numerical approximation $\bm{f}_{i}^{[n]}\approx\bm{f}(s_{i})$ to satisfy

$\bm{y}(s_{i})=\alpha\,\bm{f}_{i}^{[n]}+\sum_{j=1}^{n}w_{j}K(s_{i},s_{j})\bm{f}_{j}^{[n]}\quad\text{for}\quad
i=1,...,n.$ (38)

Introducing the vectors
$\underline{\bm{f}}^{[n]}=((\bm{f}_{1}^{[n]})^{T},...,(\bm{f}_{n}^{[n]})^{T})^{T}$
and $\underline{\bm{y}}=(\bm{y}(s_{1})^{T},...,\bm{y}(s_{n})^{T})^{T}$,
equation (38) can be written compactly as

$\underline{\bm{y}}=\left(\alpha\,I+\underline{K}\,\underline{W}\right)\,\underline{{\bm{f}}}^{[n]}.$
(39)

Here $I$ denotes the $3n\times 3n$ identity matrix, and

$\displaystyle\underline{W}=\mathrm{diag}({\bm{w}})\otimes{\bf
I},\quad\text{and}\quad\underline{K}=\left(\begin{array}[]{ccc}K(s_{1},s_{1})&\dots&K(s_{1},s_{n})\\\
\vdots&\ddots&\vdots\\\ K(s_{n},s_{1})&\dots&K(s_{n},s_{n})\\\
\end{array}\right)\in\mathbb{R}^{3n\times 3n}$ (43)

with $\otimes:\mathbb{R}^{n_{1}\times m_{1}}\times\mathbb{R}^{n_{2}\times
m_{2}}\rightarrow\mathbb{R}^{(n_{1}n_{2})\times(m_{1}m_{2})}$ the Kronecker
product of matrices and ${\bf I}$ the $3\times 3$ identity matrix. We then
approximate (37) by the same quadrature formula

$\displaystyle{\phi}_{M}(\bm{f})\approx\sum_{i=1}^{n}w_{i}M(s_{i}){\bm{f}}_{i}^{[n]}=(\mathbb{1}^{T}\otimes{\bf
I})\underline{M}\,\underline{W}\,\underline{{\bm{f}}}^{[n]}:={\phi}^{[n]}_{M},$
(44)

where

$\underline{M}=\left(\begin{array}[]{ccc}M(s_{1})&&0\\\ &\ddots&\\\
0&&M(s_{n})\\\ \end{array}\right)\in\mathbb{R}^{3n\times 3n}$ (45)

and $\mathbb{1}=(1,\dots,1)^{T}\in\mathbb{R}^{n}.$ Here we have used
${\phi}^{[n]}_{M}$ to denote the approximation of ${\phi}_{M}({\bm{f}})$
obtained by quadrature. After inserting the solution of (39), we obtain

${\phi}^{[n]}_{M}=(\mathbb{1}^{T}\otimes{\bf
I})\underline{M}\,\underline{W}\left(\alpha
I+\underline{K}\,\underline{W}\right)^{-1}\underline{\bm{y}}.$ (46)

######  Remark 4.1.

The linear functional (37) has a regularizing effect on numerical solutions to
first-kind integral equation (35) in the case of a constant fiber radius. This
is illustrated via singular value expansion for a thin ring using $M={\bf I}$
in Appendix B. This effect applies whenever such a functional is the final
quantity of interest, such as in the case of a rigid slender body.

####  4.1.1 Convergence and error bounds

We are interested in obtaining an estimate for the error when approximating
(37) by its discrete approximation (46), which we denote by

$\mathbf{d}^{[n]}={\phi_{M}}({\bm{f}})-\phi^{[n]}_{M}=\int_{-L}^{L}\\!M(s)\bm{f}(s)ds-\sum_{j=1}^{n}w_{j}M(s_{j})\bm{f}_{j}^{[n]}.$
(47)

This error will depend on the error committed in the numerical approximation
of (36) by the solution $\underline{\bm{f}}^{[n]}$ of (39). For this reason,
we first analyze the convergence of Nyström’s method in using (39) to
approximate the solution of (36) [1, Chapt. 12.4]. At each quadrature node, we
define the error of this approximation as

$\mathbf{e}^{[n]}_{i}:=\bm{f}(s_{i})-\bm{f}_{i}^{[n]},\quad\text{for}\quad
i=1,\dots,n,$ (48)

and let
$\mathbf{\underline{e}}^{[n]}:=((\mathbf{e}^{[n]}_{1})^{T},\dots,(\mathbf{e}^{[n]}_{n})^{T})^{T}$
denote the error vector. We want to show that
$\|\mathbf{\underline{e}}^{[n]}\|_{\infty}\rightarrow 0$ as
$n\rightarrow\infty$. Let
$\underline{{\bm{f}}}:=(\bm{f}(s_{1})^{T},\dots,\bm{f}(s_{n})^{T})^{T}$ and
define
$\mathbf{\underline{\tau}}^{[n]}:=(\mathbf{\tau}_{1}^{T},\dots,\mathbf{\tau}_{n}^{T})^{T}$
with components

$\mathbf{\tau}_{i}:=\bm{y}(s_{i})-\alpha\bm{f}(s_{i})-\sum_{j}^{n}K_{i,j}w_{j}\bm{f}(s_{j}),$
(49)

the truncation error for the discrete second kind equation (39) – i.e. the
residual obtained replacing $\underline{\bm{f}}^{[n]}$ by
$\underline{{\bm{f}}}$ in (39). We obtain

$\qquad\left(\alpha\,I+\underline{K}\,\underline{W}\right)\,\underline{{\bm{f}}}=\underline{\bm{y}}-\mathbf{\underline{\tau}}^{[n]}.$
(50)

It is easily seen using (36) that

$\mathbf{\tau}_{i}=\int_{-L}^{L}K(s_{i},s^{\prime})\bm{f}(s^{\prime})ds^{\prime}-\sum_{j}^{n}K_{i,j}w_{j}\bm{f}(s_{j}),$
(51)

which is simply quadrature error, and for any convergent quadrature formula we
have

$\lim_{n\rightarrow\infty}\|\mathbf{\tau}^{[n]}\|_{\infty}=0.$ (52)

We next bound the norm of the error $\mathbf{\underline{e}}^{[n]}$ by the norm
of $\mathbf{\underline{\tau}}^{[n]}$ to prove the convergence of the method.
Subtracting (39) from (50) we obtain a linear system satisfied by
$\mathbf{\underline{e}}^{[n]}$:

$\left(\alpha\,I+\underline{K}\,\underline{W}\right)\,\mathbf{\underline{e}}^{[n]}=-\mathbf{\underline{\tau}}^{[n]}.$
(53)

From [1, Chapt. 12.4] Theorem 12.4.4 and equation (12.4.51) we have that for
sufficiently large $n$, say $n\geq n^{*}$, the matrix
$\left(\alpha\,I+\underline{K}\,\underline{W}\right)$ is invertible and

$\|\left(\alpha\,I+\underline{K}\,\underline{W}\right)^{-1}\|_{\infty}\leq
C_{1}\qquad\forall n\geq n^{*}.$ (54)

Thus we can conclude that

$\|\mathbf{\underline{e}}^{[n]}\|_{\infty}\leq\|\left(\alpha\,I+\underline{K}\,\underline{W}\right)^{-1}\|_{\infty}\,\|\mathbf{\underline{\tau}}^{[n]}\|_{\infty}\leq
C_{1}\,\|\mathbf{\underline{\tau}}^{[n]}\|_{\infty}.$ (55)

Since $C_{1}$ is independent of $n$ for $n\geq n^{*}$ and
$\|\mathbf{\underline{\tau}}^{[n]}\|_{\infty}\rightarrow 0$ as
$n\rightarrow\infty$, this implies that

$\lim_{n\rightarrow\infty}\|\mathbf{\underline{e}}^{[n]}\|_{\infty}=0.$

Consider now the quadrature error

$\displaystyle\delta^{[n]}:=$
$\displaystyle\int_{-L}^{L}M(s)\bm{f}(s)ds-\sum_{j=1}^{n}w_{j}M(s_{j})\bm{f}(s_{j}).$
(56)

From (47) we obtain

$\mathbf{d}^{[n]}=\delta^{[n]}-\sum_{j=1}^{n}w_{j}M(s_{j})\mathbf{e}_{j},$
(57)

and using (53) the total discretization error for our methods is given by

$\mathbf{d}^{[n]}=(\mathbb{1}^{T}\otimes\mathbf{I})\underline{W}\,\underline{M}(\alpha
I+\underline{K}\,\underline{W})^{-1}\underline{\tau}^{[n]}+\delta^{[n]}$ (58)

Since both $\delta^{[n]}$ and $\mathbf{\tau}^{[n]}$ are quadrature errors,
$\|(\alpha I+\underline{K}\,\underline{W})^{-1}\|\leq C_{1}$ for all $n\geq
n^{*}$, and $M$ is bounded, the method converges at the same rate as the
underlying quadrature.

###  4.2 Application to the slender body model and convergence tests

We next apply our numerical method from Section 4.1 to approximate the force
and torque on a slender body. Note that the equations (8) are given by setting
$M(s)=\mathbf{I}$ and $M(s)=\widehat{\bm{X}}(s)$ in the functional (37). That
is,

$\bm{F}={\phi}_{\mathbf{I}}(\bm{f})\quad\text{and}\quad\bm{T}={\phi}_{\mathbf{\widehat{\bm{X}}}}(\bm{f}).$
(59)

Letting $\alpha=2\log(\eta)$ and

$\displaystyle
K(s,s^{\prime})=\bm{S}_{\epsilon,\eta}(s,s^{\prime})+\frac{\epsilon^{2}r^{2}(s^{\prime})}{2}\bm{D}_{\epsilon}(s,s^{\prime}),$
(60)
$\bm{y}(s)=-8\pi\mu(\bm{v}-\widehat{\bm{X}}(s)\boldsymbol{\omega}-\bm{u}_{0}(\bm{X}(s,t),t)),$
(61)

our model (4) is of the form (36), and we may write the discretization of (4)
in the form (39). Here we have introduced the hat operator
$\widehat{\cdot}:\mathbb{R}^{3}\rightarrow\mathfrak{so}(3)$ which maps vectors
in $\mathbb{R}^{3}$ to $3\times 3$ skew symmetric matrices by

$\bm{\omega}=\left(\begin{array}[]{c}\omega_{1}\\\ \omega_{2}\\\ \omega_{3}\\\
\end{array}\right)\mapsto\widehat{\boldsymbol{\omega}}=\left(\begin{array}[]{ccc}0&-\omega_{3}&\omega_{2}\\\
\omega_{3}&0&-\omega_{1}\\\ -\omega_{2}&\omega_{1}&0\\\ \end{array}\right).$
(62)

Here $\mathfrak{so}(3)$ is the Lie algebra of $SO(3)$, and such that
$\mathbf{\boldsymbol{\omega}}\times\bm{v}=\hat{\boldsymbol{\omega}}\bm{v}$ for
$\boldsymbol{\omega},\bm{v}\in\mathbb{R}^{3}$.

Denote the numerical approximations to (59) by

$\bm{F}^{[n]}={\phi}^{[n]}_{\mathbf{I}}\quad\text{and}\quad\bm{T}^{[n]}={\phi}^{[n]}_{\mathbf{\widehat{\bm{X}}}}.$
(63)

Defining the matrices $\Phi$ and $\Psi\in\mathbb{R}^{3\times 3n}$ as

$\displaystyle\Phi=$
$\displaystyle(\mathbb{1}^{T}\otimes\mathbf{I})\underline{W}\left(\alpha
I+\underline{K}\,\underline{W}\right)^{-1},$ (64) $\displaystyle\Psi=$
$\displaystyle(\mathbb{1}^{T}\otimes\mathbf{I})\underline{W}\,\underline{X}\left(\alpha
I+\underline{K}\,\underline{W}\right)^{-1},$ (65)

we may then write equations (63) as

$\displaystyle\bm{F}^{[n]}$
$\displaystyle=\Phi\underline{\bm{y}}\quad\text{and}\quad\bm{T}^{[n]}=\Psi\underline{\bm{y}}.$
(66)

In the next section we perform convergence tests for our discrete model (66)
for both a thin ring and a prolate spheroid. With these geometries we are able
to calculate accurate reference solutions against which we can compare the
accuracy of our numerical solution. We provide an error bound (see (75) below)
depending on the speed of convergence of the quadrature rule, this error bound
is confirmed by our numerical experiments. Furthermore, we will look at how
the conditioning of the linear system associated with the discretized integral
operator improves as the regularization parameter $\eta$ is increased from
$\eta=1$ to $\eta>1$.

####  4.2.1 Convergence of numerical method for closed loop geometry

By applying the formula (58), we now show how one can achieve spectral
convergence in the case of a closed fiber geometry with constant radius
$\epsilon$ and periodic integration domain. In this setting, we will use
trapezoidal quadrature. We begin by bounding the norms of the integration
kernels to which we apply the trapezoidal quadrature rules to, namely the
integrals (34) and (37). Using this, and some smoothness assumptions, we are
able bound the quadrature errors $\tau_{i}^{[n]}$ and $\delta^{[n]}$ using
classical error estimates. This leads to a bound on the total error
$\mathbf{d}^{[n]}$ for both the force and torque calculation.

Let $C_{2}$ be a constant such that

$\|\bm{f}(s^{\prime})\|_{\infty}\leq C_{2}\quad\text{ for }\quad s\in[-L,L].$
(67)

From the definition of $K(s,s^{\prime})$ (equations (5), (6), and (60)) in the
constant radius case, we observe that

$\|K(s,s^{\prime})\|_{\infty}\leq\frac{3}{2\epsilon}$ (68)

with equality when $s=s^{\prime}$. From equation (59) we have
$\|M(s)\|_{\infty}=1$ for the force calculation, while for the torque
calculation, $M(s)=\widehat{\bm{X}}(s)$ and therefore

$\|M(s)\|_{\infty}\leq\max_{s\in[-L,L]}\|\bm{X}(s)\|_{1}.$ (69)

Therefore we can bound the integration kernels of (34) and (37) by

$\|K(s,s^{\prime})\bm{f}(s^{\prime})\|_{\infty}\leq\frac{3}{2\epsilon}C_{2}$
(70)

and

$\|M(s)\bm{f}(s)\|_{\infty}\leq\|M(s)\|_{\infty}C_{2}.$ (71)

Note that in the constant radius case, $K(s,s^{\prime})$ has the same
regularity as $\bm{X}(s)$. If we assume that $\bm{X}(s)$, $\bm{f}(s)$ and
$M(s)$ are analytic, then using [38, Theorem 3.2] we can bound the trapezoidal
rule quadrature error from equation (49) by

$\|\tau_{i}^{[n]}\|_{\infty}\leq\frac{6LC_{2}}{\epsilon(e^{an}-1)}\quad\text{for}\quad
i=1,...,n.$ (72)

Similarly, we can bound equation (56) by

$\|\delta^{[n]}\|_{\infty}\leq\frac{4L\|M(s)\|_{\infty}C_{2}}{e^{an}-1}.$ (73)

Here $a$ is some constant. Using equation (58), the total discretization error
is therefore bounded as

$\|\mathbf{d}^{[n]}\|_{\infty}\leq\left(\|(\mathbb{1}^{T}\otimes\mathbf{I})\underline{W}\,\underline{M}(\alpha
I+\underline{K}\,\underline{W})^{-1}\|_{\infty}\frac{3}{2\epsilon}+\|M(s)\|_{\infty}\right)\frac{4LC_{2}}{e^{an}-1}.$
(74)

Using that $\|\underline{M}\|_{\infty}\leq\|{M(s)}\|_{\infty}$,
$\|\underline{W}\|_{\infty}=\frac{2L}{n}$ and $C_{1}$ is given by equation
(54), this simplifies to

$\|\mathbf{d}^{[n]}\|_{\infty}\leq\left(\frac{6C_{1}L}{2\epsilon}+1\right)\frac{4L\|{M(s)}\|_{\infty}C_{2}}{e^{an}-1}.$
(75)

Hence, the method shares the same exponential convergence as the underlying
trapezoidal rule. We remark that one could perform an analogous analysis for
open ended fiber geometries with, e.g., Gauss-Lobatto quadrature, and derive
similar results. Furthermore, we also remark that one could require less
stringent regularity assumptions on the integration on the kernels or the
fiber centreline $\bm{X}(s)$, e.g., $M(s)\bm{f}(s)\in C^{2m+2}[-L,L]$. Then
[2, Thm. 5.5] can be used to derive asymptotic error estimates for
$\tau_{i}^{[n]}$ and $\delta^{[n]}$ of order $O(h^{2m+2})$. Nonetheless, we do
observe spectral convergence in numerical experiments in the following
sections, as predicted by the bound (75).

####  4.2.2 Thin ring translating with unit velocity

As a convergence test, we use (66) to calculate the force on a thin ring of
unit length in the $xy$-plane translating in the $z$ direction with unit
velocity in zero background flow. We will consider both the first- and second-
kind formulations of the model. In this setting, the force on the ring can be
calculated to arbitrarily high precision by evaluating elliptic integrals,
which can be used as a reference solution. For a circular centerline
parametrized by

$\bm{X}(s)=\left(\frac{\cos(\pi s)}{2\pi},\frac{\sin(\pi
s)}{2\pi},0\right)^{T},$

the $z$-component of our unregularized ($\eta=1$) model becomes

$8\pi\mu=-\int_{-\frac{1}{2}}^{\frac{1}{2}}{\frac{\sqrt{2}\,\pi\left(3\,{\epsilon}^{2}{\pi}^{2}-\cos\left(2\pi\,(s-s^{\prime})\right)+1\right)}{\left(2\,{\epsilon}^{2}{\pi}^{2}-\cos\left(2\pi\,(s-s^{\prime})\right)+1\right)^{3/2}}}f^{z}(s^{\prime})\,ds^{\prime}.$
(76)

As in the straight-but-periodic geometry of Section 3.1, the eigenfunctions of
this operator are the Fourier modes $f_{k}^{z}(s)=\exp(i2\pi ks)$. The force
$\bm{F}=(F,0,0)^{T}$ is therefore given by

$F=\int_{-\frac{1}{2}}^{\frac{1}{2}}\\!f^{z}(s){\rm
d}s=\frac{8\pi\mu}{\lambda^{z}_{0}}$ (77)

where $\lambda^{z}_{0}$ is the $k=0$ eigenvalue. This can be found by
evaluating the integral in equation (76) with $f^{z}(s)=f^{z}_{0}(s)=1$, which
gives

$\lambda^{z}_{0}=-c_{\epsilon}\left(2\,\phi_{K}\left(c_{\epsilon}\right)+\phi_{E}\left(c_{\epsilon}\right)\right).$
(78)

Here $c_{\epsilon}=\sqrt{\left({\epsilon}^{2}{\pi}^{2}+1\right)^{-1}}$, and

$\phi_{K}(x)=\int_{0}^{1}\\!{\frac{1}{\sqrt{1-{\theta}^{2}}\sqrt{1-{x}^{2}{\theta}^{2}}}}\,{\rm
d}\theta\quad\text{and}\quad\phi_{E}(x)=\int_{0}^{1}\\!{\frac{\sqrt{1-{x}^{2}{\theta}^{2}}}{\sqrt{1-{\theta}^{2}}}}\,{\rm
d}\theta$ (79)

are the complete elliptic integrals of the first and second kind,
respectively.

For $\epsilon=0.05,0.025,0.01$ and $0.005$, equation (76) is discretized using
trapezoidal quadrature, and our numerical method is used to approximate $F$ by
equation (66). Figure 3 plots the error as a function of $n$ for four
different values of $\epsilon$. We observe spectral convergence of the error
to machine precision, which is consistent with our error estimates from
expression (75). We note that the condition number of the unregularized
discrete integral operator grows exponentially as $n$ increases, as shown in
Figure 4(a). However, because we are considering a rigid fiber with constant
radius, the regularizing effect of computing $F$ comes into effect (see Remark
4.1) and thus the effect of this ill-conditioning is not apparent in the final
force calculation. This may be contrasted with the case of the prolate
spheroid, where, as we will see in Section 4.2.3, the conditioning does have a
noticeable effect on the error. Nevertheless, we note that by setting $\eta>1$
we can improve the condition number of the linear system (see Figure 4(b)).

Figure 3: The approximate drag force $F^{[n]}$ on a thin ring translating
broadwise with unit velocity converges with spectral accuracy to the true
force $F$.

(a) Unregularized ($\eta=1$).

(b) Regularized ($\eta=1.5$).

Figure 4: The condition numbers associated with the discretized versions of
the unregularized ($\eta=1$) and regularized ($\eta=1.5$) slender body models
for calculating the force on a thin ring. Note the change in scale between the
two figures.

####  4.2.3 Prolate spheroid with artificial fluid velocity field

We next use (66) to compute the drag force for a stationary prolate spheroid
immersed in an artificial fluid velocity field. The particle centerline is
aligned in the $z$-direction, parameterized by $\bm{X}(s)=(0,0,s)^{T}$,
$s\in[-1,1]$. The fluid velocity field $\bm{u}(s)=(u(s),0,0)^{T}$ is designed
such that $\bm{f}(s)=(f^{x}(s),0,0)^{T}$ is a known analytic function. We
choose this function to be a Gaussian
$f^{x}(s)=\exp\left(-\frac{s^{2}}{\epsilon^{2}}\right)$ such that the force
decays to zero at the fiber endpoints and use high order Gauss-Lobatto
quadrature for the discretization of the integral operator. Denote the set of
$n$ quadrature nodes by $\\{s_{i}\\}^{n}_{i=1}$. Inserting the above
expression for $f^{x}(s)$ into our model (17), the fluid velocity at $s_{i}$
is found by solving the integral

$u(s_{i})=\frac{-1}{8\pi}\left(2\log(\eta)\,\exp\left(-\frac{s_{i}^{2}}{\epsilon^{2}}\right)+\int_{-1}^{1}{\frac{{\epsilon}^{2}r\left(s_{i}\right)^{2}+\frac{1}{2}{\epsilon}^{2}r\left(s^{\prime}\right)^{2}+\left(s_{i}-s^{\prime}\right)^{2}}{\left({\epsilon}^{2}r\left(s_{i}\right)^{2}+\left(s_{i}-s^{\prime}\right)^{2}\right)^{3/2}}}\exp\left(-\frac{s^{\prime
2}}{\epsilon^{2}}\right)ds^{\prime}\right)$ (80)

where the ellipsoidal radius function is given by equation (2). We also take
the viscosity $\mu=1$. To solve for $u(s_{i})$ for $i=1,...,n$, the integral
in equation (80) is evaluated to machine precision using MATLAB’s built-in
integral function, which uses adaptive quadrature. For this fluid velocity
field, the total force $\bm{F}=(F,0,0)^{T}$ on the ellipsoid is found by

$F=\int_{-1}^{1}\exp\left(-\frac{s^{2}}{\epsilon^{2}}\right)ds=\sqrt{\pi}\,\epsilon\,\mathrm{erf}(\frac{1}{\epsilon}).$
(81)

We compute numerical approximations to $F$ using equation (66) for four
choices of $\epsilon$. We initially set $\eta=1$ and compute these numerical
approximations for the non-regularized, first-kind equation. The errors are
presented in Figure 5(a). We see that the error converges spectrally up to a
certain point where the method begins to diverge due to numerical
instabilities and poor conditioning of the discrete integral operator, which
is plotted in Figure 5(b).

However, by choosing $\eta>1$, we can amend the condition number and therefore
improve the accuracy of the numerical solution. In Figure 6, we fix
$\epsilon=0.025$ and calculate the errors for four choices of $\eta$. We see
from Figure 6(a) that the error converges spectrally to machine precision for
all such choices of $\eta$. This level of accuracy was unattainable for the
non-regularized problem. Furthermore, we observe from Figure 6(b) that the
condition number of the discrete integral operator is bounded by a value that
becomes smaller for larger $\eta$.

(a)

(b)

Figure 5: The errors (a) and condition numbers (b) associated with the
unregularized ($\eta=1$) numerical method for the calculation of the force on
a prolate spheroid for different values of $\epsilon$.

(a)

(b)

Figure 6: The errors (a) and condition numbers (b) associated with the
regularized numerical method for the calculation of the force on a prolate
spheroid for $\epsilon=0.025$. Similar results are observed for other values
of $\epsilon$.

###  4.3 Spectrum of the slender body operator in different geometries

One important unresolved question about the slender body model (4) is the
effect of different geometries, including curvature, endpoints, and non-
uniform fiber radius, on the spectrum of the integral operator. Although we
cannot analytically determine the spectrum of the continuous operator in
general, we _can_ determine the eigenvalues of the discrete operator
$(2\log(\eta)I+\underline{K}\,\underline{W})$ (39). We consider first the
unregularized version $\eta=1$, recalling that in the straight-but-periodic
geometry of Section 3.1, the continuous operator was provably negative
definite. Ideally we would like to see evidence that this negative
definiteness persists in general geometries, as this would be the physically
correct behavior and also would agree with the underlying slender body PDE
operator (10).

We begin by calculating the eigenvalues $\\{\lambda_{i}\\}_{i=1}^{3n}$ of
$\underline{K}\,\underline{W}$ for the thin ring. Letting
$\lambda_{\max}=\max_{i}(\lambda_{i})$, in Figure 7(a) we plot
$\lambda_{\max}$ versus $n$ for five different values of $\epsilon$. Note that
for very large $n$ relative to $\epsilon^{-1}$ (roughly $n=O(\epsilon^{-2})$),
we begin to see numerical error resulting in very small positive eigenvalues
of $\underline{K}\,\underline{W}$ (denoted by red markers). However, the
magnitude of these positive eigenvalues are on the order of machine precision
and therefore can be attributed to round-off errors.

We next consider the effects of endpoints and a non-uniform radius by
calculating the eigenvalues of $\underline{K}\,\underline{W}$ for a slender
prolate spheroid (2), keeping in mind the above level of numerical error. In
Figure 7(b) we again plot $\lambda_{\max}$ versus $n$ for four different
values of $\epsilon$. Again for $n=O(\epsilon^{-2})$ we begin to see small
positive eigenvalues which are significantly larger than for the thin ring
(around $O(10^{-10})$). However, the magnitude of the positive eigenvalues is
still very small and bounded as $n$ increases. It is not clear whether this is
a numerical artifact or an actual eigenvalue crossing 0 for the continuous
operator. At any rate, the non-regularized operator would never actually be
used for simulations with such large $n$ because the condition number of
$\underline{K}\,\underline{W}$ is prohibitive (see Figure 5(b)). It appears
that a very reasonable choice of regularization parameter $\eta$ will ensure
that none of these near-zero eigenvalues actually cross zero.

(a) Thin ring

(b) Spheroid

Figure 7: Magnitude of the maximum eigenvalue of the non-regularized discrete
slender body operator $\underline{K}\,\underline{W}$. Blue markers mean that
$\lambda_{\max}<0$ implying $\underline{K}\,\underline{W}$ is negative
definite, while red markers mean that $\lambda_{\max}>0$.

As a final test, we calculate the spectrum of $\underline{K}\,\underline{W}$
for randomly but systematically generated curvy fibers with complicated shapes
(Figure 8). Here the magnitude of the fiber’s deviation from a straight line
is controlled by a small parameter $\delta\geq 0$. The fiber shapes are
generated by interpolating $m$ points $(x_{i},y_{i},z_{i})\in\mathbb{R}^{3}$,
$i=1,...,m$, with cubic splines. Here $z_{i}=(i-1)\frac{2L}{m}$ while
$x_{i},y_{i}\in[-\delta,\delta]$ are given by a random number generator and
are of size at most $\delta$. Setting $\delta=0$ corresponds to a straight
fiber. Examples of the fiber centerline for $m=10$ and four different values
of $\delta$ are given in Figure 8.

(a)

(b)

(c)

(d)

Figure 8: The centerlines of four curved fiber shapes.

We fix $\epsilon=0.1$ and use the spheroidal radius function (2). Taking
$m=10$, we generate 6 different curvy fibers for different magnitudes
$\delta\in[0,\frac{1}{10}]$. For each fiber we compute the spectrum
$\\{\lambda^{\delta}_{i}\\}_{i=1}^{n}$ of its corresponding (non-regularized)
integral operator $\underline{K}\,\underline{W}$. We plot the most positive
eigenvalue $\lambda^{\delta}_{\max}=\max_{i}(\lambda^{\delta}_{i})$ for each
fiber in Figure 9(a). For each value of $\delta$ we note that there is an
eigenvalue crossing zero when $n=O(\epsilon^{-2})$. As $\delta$ increases and
the magnitude of the curviness of the fiber increases, we can note a slight
increase in the magnitude of the largest positive eigenvalue, but
$\lambda^{\delta}_{\max}$ is still small – roughly $O(10^{-8})$. Again, we can
be assured to have a negative spectrum bounded away from 0 by a reasonable
choice of regularization $\eta>1$. This effect is displayed in Figure 9(b),
which shows the maximum eigenvalue $\lambda^{\delta,\eta}_{\mathrm{max}}$ of
the now regularized discrete integral operator
$(2\log(\eta)I+\underline{K}\,\underline{W})$ for a fixed value of $\epsilon$
and $\delta$ and varying values of $\eta$. We see here that for all choices of
$\eta>1$ in this range, the spectrum of
$(2\log(\eta)I+\underline{K}\,\underline{W})$ remains negative definite.

(a) $\lambda^{\delta}_{\max}$ for different values of $\delta$ and $\eta=1$ .

(b) $\lambda^{\delta,\eta}_{\mathrm{max}}$ for different values of $\eta$ with
$\epsilon=0.1$ and $\delta=0.001$.

Figure 9: The maximum eigenvalues for the unregularized (a) and regularized
(b) discrete integral operators for the curved fibers. The color blue denotes
a negative maximum eigenvalue and red denotes a positive maximum eigenvalue.

##  5 Dynamics of curved rigid fibers

We next use the slender body model (4) and the discretization procedure of
Section 4 to simulate the dynamics of curved rigid fibers in Stokes flow.
After outlining the dynamical equations, we validate the model against known
dynamics for a slender prolate spheroid. Finally, we compare the rotational
dynamics of randomly curved fibers as in Figure 8 to straight fibers.

###  5.1 Dynamical equations

The dynamics of the slender body are governed by the rigid body equations. The
angular momentum $\bm{m}$ of a rigid particle with torque $\bm{T}(t)$ is found
by solving

$\dot{\bm{m}}=\bm{m}\times\boldsymbol{\omega}+\bm{T},$ (82)

where $\boldsymbol{\omega}=J^{-1}\bm{m}$ for moment of inertia tensor $J$.
Each of these quantities are given in the fiber frame of reference. The fiber
orientation (with respect to a fixed inertial frame of reference) is specified
using Euler parameters $q\in\mathbb{R}^{4}$ which satisfy the constraint
$||q||_{2}=1$ and are determined by solving the ODE

$\dot{q}=\frac{1}{2}q\cdot w,$ (83)

where $w=(0,\boldsymbol{\omega}^{\mathrm{T}})^{\mathrm{T}}\in\mathbb{R}^{4}$
and $\cdot$ denotes the Hamilton product of two quaternions [9]. A vector
$\bm{x}_{p}$ in the fiber reference frame can be rotated to a vector
$\bm{x}_{T}$ in an inertial co-translating reference frame via
$\bm{x}_{T}=Q\bm{x}_{p}$ where the rotation matrix $Q$ is the image of $q$
under the Euler-Rodriguez map. We refer the reader to [9] for details on
quaternion algebra and rigid body mechanics.

The translational dynamics are given by Newton’s second law

$\dot{\bm{p}}=\bm{F},$ (84)

where $\bm{p}=\bm{v}/m$ is the inertial frame linear momentum for a fiber of
mass $m$. The position of the fiber center of mass is found by solving

$\dot{\bm{x}}=\bm{v}.$ (85)

The ODEs (82) - (85) are integrated using the second order Strang splitting
method of [34].

Recall the equations (66) for $\bm{F}^{[n]}$ and $\bm{T}^{[n]}$. Since
$\bm{F}^{[n]}$ and $\bm{T}^{[n]}$ depend linearly on the linear and angular
momenta $\bm{p}$ and $\bm{m}$, we may update them according to the linear
equation

$\left(\begin{array}[]{c}\bm{F}^{[n]}\\\ \bm{T}^{[n]}\\\
\end{array}\right)=A\left(\begin{array}[]{c}\bm{p}\\\ \bm{m}\\\
\end{array}\right)+\bm{b},$ (86)

where $A$ is a negative definite dissipation matrix and $\bm{b}$ is due to the
background fluid velocity and is independent of $\bm{p}$ and $\bm{m}$. We have
that

$A=\left(\begin{array}[]{cc}\Phi\,\left(\mathbb{1}\otimes(I/m)\right),&\Phi\,\left(-\underline{X}(\mathbb{1}\otimes
J^{-1})\right)\\\
\Psi\left(\mathbb{1}\otimes(I/m)\right),&\Psi\left(-\underline{X}(\mathbb{1}\otimes
J^{-1})\right)\end{array}\right)\quad\mathrm{and}\quad\bm{b}=-\left(\begin{array}[]{c}\Phi\underline{\bm{u}}\\\
\Psi\underline{\bm{u}}\\\ \end{array}\right),$ (87)

where $m$ and $J$ are the filament mass and moment of inertia tensor,
respectively. We have also introduced the vector
$\underline{\bm{u}}=(\bm{u}_{0}(\bm{X}(s_{1}))^{T},...,\bm{u}_{0}(\bm{X}(s_{n}))^{T})^{T}$
containing the background fluid velocities at the location of the quadrature
nodes along the centerline.

####  5.1.1 Overview and cost of algorithm

The algorithm used to compute the dynamics of a slender fiber is as follows:

  1. 1.

Define particle geometry $\bm{X}(s)$, $\epsilon$, regularization parameter
$\eta$ and discretization parameter $n$.

  2. 2.

Choose a quadrature rule and compute the matrices $\underline{W}$ and
$\underline{K}$.

  3. 3.

Compute the matrices $\Phi$, $\Psi$ and $A$ from equations (64), (65) and
(87).

  4. 4.

Time loop: for $t=0,\Delta t,...,m\Delta t$

a) Compute $\bm{F}^{[n]}$ and $\bm{T}^{[n]}$ using equation (86)

b) Numerically integrate the ODEs (82) - (85)

For step (2), we use the trapezoidal quadrature rule for closed fibers (i.e.,
a periodic integration interval) or Gauss-Lobatto quadrature rule for fibers
with open ends. For step (4b), we use a splitting method [34]. We note that
for simulations where the fluid velocity field is calculated from a direct
numerical simulation of the Navier-Stokes equations, the fluid field needs to
be approximated onto the centerline of the particle using an interpolation
method [35].

The above algorithm exploits the rigidity of the fiber by using the fact that
$A$, $\Phi$ and $\Psi$ are constant in time and therefore can be computed
outside of the time loop. The calculation of these matrices, which involves
solving a linear system, is the most costly operation in the algorithm but
only needs to be done once. If, for example, Gaussian elimination is used,
this step has complexity of $O(n^{3})$. Within the time loop, however, the
most costly operation is the calculation of $\bm{F}^{[n]}$ and $\bm{T}^{[n]}$,
which involves only $3\times 3n$ by $3n\times 1$ matrix-vector products, which
has $O(n)$ complexity. We assume that the cost of numerically integrating the
ODEs is negligible compared to this. For a single fiber, the total complexity
of the algorithm is therefore $O(n^{3}+nm)$, where $m$ is the total number of
time steps used in the simulation. Hence, for simulations where many time
steps are needed, the algorithm scales by $O(n)$. We remark that for problems
where the background flow is zero, the cost of computing $\bm{F}^{[n]}$ and
$\bm{T}^{[n]}$ is independent of $n$ (after $A$ has been computed) and
therefore is $O(1)$. This is relevant, for example, when simulating fibers
sedimenting in a still fluid under the influence of gravity [26].

###  5.2 Numerical validation of model dynamics

####  5.2.1 Dissipation matrix of a prolate spheroid

Here we compare our model and numerical method with accurate closed form
expressions for the force and torque given by Brenner [3] and Jeffery [14].
These expressions are valid for an ellipsoid when the fluid Jacobian is
approximately constant throughout the volume of the particle. When the flow is
linear, these terms are essentially exact and therefore serve as a good
reference model against which to validate our model.

The purpose of this numerical experiment is therefore twofold. Firstly, we aim
to show that our model converges to the reference model as
$\epsilon\rightarrow 0$. This is primarily to validate the accuracy of the
model. However, the numerical approximation of the force and torques also
introduces a numerical error that is related to the discretization parameter
$n$. Clearly, taking $n$ too small means that we will not exploit the accuracy
of the model to its entirety. On the other hand, it is unwise to take $n$ as
large as possible as this will incur unnecessary computational costs that go
to minimizing numerical error beyond the accuracy of the model. So the second
question we address here is what is an ideal choice of discretization
parameter to use such that the numerical error is roughly the same as the
modeling error.

Using $\eta=1+\epsilon^{2}$, the dissipation matrix for our slender body model
$A$ is numerically approximated by equation (87). The reference dissipation
matrix $A_{sph}$ is found using the closed form expressions from Jeffery and
Brenner, which are given in Appendix A. Denote the six eigenvalues of $A$ and
$A_{sph}$, by $\lambda_{i}$ and $\lambda_{i}^{sph}$, respectively. Note that
due to symmetry of the spheroid, $\lambda_{1}=\lambda_{2}$ and
$\lambda_{4}=\lambda_{5}$ and similarly for the eigenvalues of $A_{sph}$.
Furthermore, the slender body model is essentially a one dimensional filament
and therefore $\lambda_{6}=0$ meaning that spinning motion about the
centerline doesn’t dissipate. This is in contrast to the Jeffrey term, which
does dissipate spinning motion. We remark that this phenomenon only occurs in
the case where the centerline is perfectly straight. Hence for curved fiber
geometries where the application of the slender body is most useful, this
nonphysical phenomenon is not observed. Note that for this geometry the
dissipation matrices are diagonal and therefore the eigenvalues are directly
proportional to the calculation of $\bm{F}^{[n]}$ and $\bm{T}^{[n]}$ in zero
background flow.

The eigenvalues of $A$ are calculated using equation (87) after discretizing
equation (36) on the Gauss-Lobatto nodes. The values
$|\lambda_{i}-\lambda_{i}^{sph}|$ for $i=1,3,4$ are plotted in Figure 10 as a
function of the discretization parameter $n$. We see that $\lambda_{i}$
converges exponentially to a point near $\lambda^{sph}_{i}$, which is likely
due to the slender body modelling error. As $\epsilon$ decreases, we make two
observations. First, for large $n$ the rate at which $\lambda_{i}$ converges
to $\lambda_{i}^{sph}$ is approximately
$-\epsilon^{2}\eta^{2}\log(\epsilon\eta)$, as seen by the horizontal dash-dot
lines. Second, as $\epsilon$ decreases, the convergence rate slows down and
one must use a larger value of $n$ to reach the most accurate solution. This
means that one must pay careful attention to the choice of $n$ when taking
$\epsilon$ to be very small. In fact, we observe empirically that the
convergence rate is approximately bounded by $e^{-4\epsilon n}$. Motivated by
this, we will take $n$ in future experiments to be approximately the
intersection of these two lines, that is

$n\approx-\frac{\log(-\epsilon^{2}\log(\epsilon))}{4\epsilon}.$ (88)

(a)

(b)

(c)

(d)

Figure 10: The difference in the dissipation matrix eigenvalues
$|\lambda_{i}-\lambda_{i}^{sph}|$, $i=1,3,4$ as a function of $n$ for three
different values of $\epsilon$. The black dashed lines are $e^{-4\epsilon n}$
and the horizontal dash-dot lines are
$-\epsilon^{2}\eta^{2}\log(\epsilon\eta)$

####  5.2.2 Prolate spheroids rotating in shear flow

Now we calculate the dynamics of a prolate spheroid in shear flow
$\bm{u}=(z,0,0)^{T}$ using our model and compare it with that of the accurate
Jeffrey model. The fiber is initially aligned at rest in the $z$-direction and
its rotational dynamics are calculated by integrating equation (82) on the
interval $t\in[0,100]$ using the splitting method of [34] with a small step
size of $h=0.01$. The simulation was repeated with $h=0.05$ with no
significant changes to the results and it is therefore concluded that time
integration errors are negligible. We repeat the experiment for 20 values of
$\epsilon$ logarithmically spaced in the interval $[0.1,0.001]$ and choose $n$
using equation (88) and $\eta=1+\epsilon^{2}$. As the spheroids are
axisymmetric, they only experience a torque about their $y$ axis, hence all of
other angular momentum components are zero (to machine precision). Three
examples of the rotational dynamics are shown in Figure 11. It is seen here
that as $\epsilon$ becomes smaller, the dynamics more closely resemble the
Jeffery model.

The relative difference between the angular momenta of the Jeffery and slender
body solutions are calculated and averaged over the simulation. This average
relative error is then plotted against the corresponding value of $\epsilon$
in Figure 12. We see that the average relative error decreases with
$\epsilon$. It is observed that in the region $0.01<\epsilon<0.1$ the error
converges at a faster rate than in the region $0.001<\epsilon<0.01$. This
could be partially explained by the fact that wider particles (larger
$\epsilon$) experience a greater resistive force as seen by the regions where
$m_{y}$ nearly reaches zero. This means that the particle spends more time in
the shear plane where the fluid velocity is zero and hence the slender body
model does not experience a large torque. However, the fluid gradient is non-
zero in this orientation and therefore the Jeffery model, which depends only
on the fluid gradient, still experiences a constant torque. This means that
compared to the Jeffery model, thicker fibers will see a greater difference in
the torque term when the fiber is aligned in the shear plane than thinner
fibers.

(a)

(b)

(c)

Figure 11: The $y$ component of a spheroid rotating in shear flow for three
different values of $\epsilon$. The solid line is the our slender body
expression and the dashed line is due to Jeffery. Figure 12: The relative
difference in $m_{y}$ between the slender body and Jeffery solutions averaged
over the interval $[0,100]$.

###  5.3 Dynamics of randomly curvy fibers

Here we simulate the dynamics of the randomly curvy fibers of Figure 8 as they
rotate in shear flow. We show how the rotational variables deviate from a
straight fiber as $\delta$ becomes larger.

We generate 100 different fiber shapes with $m=10$ using 10 different values
of $\delta$ logarithmically spaced in the interval $[5\times 10^{-5},5\times
10^{-2}]$. The 100 fibers are placed in shear flow $\bm{u}=(z,0,0)^{T}$ and
their rotational dynamics are calculated on the interval $t\in[0,100]$. The
moment of inertia tensor is approximated by placing point masses along the
centerline and using the formula

$J_{i,i}=\sum_{j=1}^{k}m_{j}(X_{i}(s_{j})-c_{i})^{2},\quad\text{for}\quad
i=1,...,3$ (89)

where $X_{i}(s_{j})$ is the $i$th component of the centerline function at the
point $s_{j}$ on the centerline and $c_{i}$ is the $i$th component of the
fiber center of mass. We weight $m_{j}$ by the cross sectional radius and use
a very large value for $k$, e.g., $k=10^{4}$. Here we take $\epsilon=0.1$ and
use the spheroidal radius function (2) along with $\eta=1+\epsilon^{2}$.

Figure 13(a) shows the angular momentum $\bm{m}$ of three fibers compared to
the $\delta=0$ case. As the $\delta=0$ fiber is perfectly straight, it does
not exhibit spinning motion and its angular momentum is purely in the $m_{y}$
component. This is in contrast to the fibers with a non-zero value of
$\delta$, in which case some of the momentum is transferred to $m_{x}$. We
therefore compare the value $\sqrt{m_{x}^{2}+m_{y}^{2}}$ between the fibers to
account for this. We see here that the $\delta=0.017783$ solution is visually
very similar to the $\delta=0$ solution. We notice a significant difference
between the other two solutions. Figure 13(b) shows the angle $\theta$ between
the $z$-axis of the particle reference frame (that is, a frame that is
rotating with the fiber) and the $x$-axis of a fixed inertial reference frame.
As the $\delta\neq 0$ fibers are not symmetric, they slowly rotate out of the
$xz$-plane and therefore after a long time, we see much more significant
discrepancies in $\theta$.

To quantify the effect that $\delta$ has on the angular momentum, we calculate
the difference in the angular momentum $\Delta m$ by subtracting off the
$\delta=0$ solution and averaging over the time interval $t\in[92,100]$, which
corresponds to roughly one period of rotation. This value is averaged over all
the fibers with similar values of $\delta$ and is expressed as a percentage of
the $\delta=0$ solution, which we denote by $\%\Delta m$. The results are
plotted in Figure 14(a). We notice that the $\%\Delta m$ is linearly
proportional to $\delta$. We observe that at the end of the simulation the
$\delta=0.0003$ fibers correspond to roughly 1% discrepancy in angular
momentum and $\delta=0.0015$ corresponds to roughly 7.5% discrepancy.

The difference in $\theta$ after one rotation as a function of $\delta$ is
displayed in Figure 14(b). The $\delta=0.0003$ solution corresponds to about a
$3^{\circ}$ difference in $\theta$ and the $\delta=0.0015$ solution
corresponds to about an $8^{\circ}$ difference.

(a)

(b)

Figure 13: The rotational variables of four fibers with different values of
$\delta$. Figure (a) shows the angular momentum and Figure (b) is the angle
between the fiber’s long axis and the $x$-axis of the inertial frame.

(a)

(b)

Figure 14: Figure (a) shows the difference in angular momentum $\Delta m$
between the curved fibers and the $\delta=0$ solution after 100 time units.
The black dashed line is $O(\delta)$. Figure (b) shows the discrepancy
$\Delta\theta$ in the angle between the centerline and the $x$-axis after
roughly one rotation.

##  6 Conclusions

We have developed an integral model for the motion of a thin filament in a
viscous fluid based on nonlocal slender body theory. The model relies on
standard singular Stokeslets and doublets but makes use of the fiber integrity
condition – the near-cancellation of angular-dependent terms along the fiber
surface – to derive an asymptotically accurate fiber velocity expression
depending only on arclength. The kernel of the resulting integral operator is
smooth and retains dependence on the (possibly varying) fiber radius in a
natural way. We can show that this integral operator is negative definite in
the simplified geometry of a straight-but-periodic filament, and we expect
similar high wavenumber behavior for curved filaments with constant radius. It
is less clear how a non-constant radius affects the spectrum; however,
numerical tests indicate that the discretized integral operator is very close
to negative definite. Nevertheless, to ensure invertibility, we develop an
asymptotically consistent regularization to convert the first-kind Fredholm
integral equation for the force density along the fiber into a second-kind
equation and show that this second-kind regularization improves the stability
and conditioning of the discretized equation. We develop a numerical method
for solving the integral equation based on the the Nyström method [1] and show
how constraining the fiber motion to be rigid can be exploited for fast
computation of fiber dynamics. We validate the method and model against the
prolate spheroid model of Jeffery [14], and apply the method to study the
rotational deviation of randomly curved rigid fibers from straight fibers.

While the fibers considered here are rigid, the model can also be used to
simulate the dynamics of semiflexible filaments. The invertibility properties
of the integral equation make it particularly well suited for handling
simulations involving inextensible fibers, where an additional line tension
equation must be solved at each time step [37, 21]. We may also consider the
effects of different choices of radius functions on the model properties,
similar to what is done in [40], although we note the necessity of smooth
decay in our radius function near the fiber endpoints.

To build on the dynamic simulations for rigid fibers, we aim to consider the
effects of fiber shape on particle deposition and aggregation. We are
especially interested in more complicated background flows, including
suspensions of rigid fibers in turbulence. The novel modelling approach
advocated herein will enable earlier explorations based on the point-particle
approach [4] to be extended to curved fibers particles.

##  7 Acknowledgments

This work has received funding from the European Unions Horizon 2020 research
and innovation programme under the Marie Sklodowska-Curie grant agreement (No.
691070) as well as the SPIRIT project (No. 231632) under the Research Council
of Norway FRIPRO funding scheme. E. Celledoni, B. Owren and B. K. Tapley would
like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge,
for support and hospitality during the programme Geometry, compatibility and
structure preservation in computational differential equations (2019) where
part of the work on this paper was undertaken. E. Celledoni and B. Owren also
acknowledge funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skodowska-Curie grant agreement No
860124. L. Ohm was supported by a University of Minnesota Doctoral
Dissertation Fellowship and NSF Postdoctoral Research Fellowship Grant No.
2001959. The authors are listed in alphabetical order.

## References

  * [1] K. Atkinson and W. Han.  Theoretical numerical analysis, volume 39.  Springer, 2005.
  * [2] K. E. Atkinson.  An introduction to numerical analysis.  1978\.
  * [3] H. Brenner.  The stokes resistance of an arbitrary particle—iv arbitrary fields of flow.  Chemical Engineering Science, 19(10):703–727, 1964.
  * [4] N. R. Challabotla, L. Zhao, and H. I. Andersson.  On fiber behavior in turbulent vertical channel flow.  Chemical Engineering Science, 153:75–86, 2016.
  * [5] S. Chattopadhyay and X.-L. Wu.  The effect of long-range hydrodynamic interaction on the swimming of a single bacterium.  Biophys. J., 96(5):2023–2028, 2009.
  * [6] A. T. Chwang and T. Y.-T. Wu.  Hydromechanics of low-Reynolds-number flow. Part 2: Singularity method for Stokes flows.  J. Fluid Mech., 67(4):787–815, 1975.
  * [7] R. Cortez and M. Nicholas.  Slender body theory for Stokes flows with regularized forces.  Commun. Appl. Math. Comput. Sci., 7(1):33–62, 2012.
  * [8] X. Fan, N. Phan-Thien, and R. Zheng.  A direct simulation of fibre suspensions.  J. Non-Newton. Fluid Mech., 74(1):113–135, 1998.
  * [9] H. Goldstein, C. Poole, and J. Safko.  Classical mechanics, 2002.
  * [10] T. Götz.  Interactions of fibers and flow: asymptotics, theory and numerics.  Doctoral dissertation, University of Kaiserslautern, 2000.
  * [11] K. Gustavsson and A.-K. Tornberg.  Gravity induced sedimentation of slender fibers.  Phys. Fluids, 21(12):123301, 2009.
  * [12] J. Hämäläinen, S. B. Lindström, T. Hämäläinen, and H. Niskanen.  Papermaking fibre-suspension flow simulations at multiple scales.  J. Engrg. Math., 71(1):55–79, 2011.
  * [13] P. C. Hansen.  Numerical tools for analysis and solution of fredholm integral equations of the first kind.  Inverse problems, 8(6):849, 1992.
  * [14] G. B. Jeffery.  The motion of ellipsoidal particles immersed in a viscous fluid.  Proc. R. Soc. Lond. A, 102(715):161–179, 1922.
  * [15] R. E. Johnson.  An improved slender-body theory for Stokes flow.  J. Fluid Mech., 99(02):411–431, 1980.
  * [16] J. B. Keller and S. I. Rubinow.  Slender-body theory for slow viscous flow.  J. Fluid Mech., 75(4):705–714, 1976.
  * [17] R. Kress, V. Maz’ya, and V. Kozlov.  Linear integral equations, volume 82.  Springer, 1989.
  * [18] E. Lauga and T. R. Powers.  The hydrodynamics of swimming microorganisms.  Rep. Progr. Phys., 72(9):096601, 2009.
  * [19] J. Lighthill.  Flagellar hydrodynamics.  SIAM review, 18(2):161–230, 1976.
  * [20] J. Martin, A. Lusher, R. C. Thompson, and A. Morley.  The deposition and accumulation of microplastics in marine sediments and bottom water from the irish continental shelf.  Sci. Rep, 7(1):10772, 2017.
  * [21] O. Maxian, A. Mogilner, and A. Donev.  An integral-based spectral method for inextensible slender fibers in stokes flow.  arXiv preprint arXiv:2007.11728, 2020.
  * [22] Y. Mori and L. Ohm.  Accuracy of slender body theory in approximating force exerted by thin fiber on viscous fluid.  arXiv preprint arXiv:2008.06829, 2020.
  * [23] Y. Mori and L. Ohm.  An error bound for the slender body approximation of a thin, rigid fiber sedimenting in Stokes flow.  Res. Math. Sci., 7(8), 2020.
  * [24] Y. Mori, L. Ohm, and D. Spirn.  Theoretical justification and error analysis for slender body theory.  Comm. Pure Appl. Math., 73(6):1245–1314, 2020.
  * [25] Y. Mori, L. Ohm, and D. Spirn.  Theoretical justification and error analysis for slender body theory with free ends.  Arch. Ration. Mech. Anal., 235:1905–1978, 2020.
  * [26] R. Newsom and C. Bruce.  The dynamics of fibrous aerosols in a quiescent atmosphere.  Physics of Fluids, 6(2):521–530, 1994.
  * [27] A. Oberbeck.  Uber stationare flussigkeitsbewegungen mit berucksichtigung der inner reibung.  J. reine angew. Math., 81:62–80, 1876.
  * [28] L. Ohm, B. K. Tapley, H. I. Andersson, E. Celledoni, and B. Owren.  A slender body model for thin rigid fibers: validation and comparisons.  Proc. of MEKiT’19, 10th Nat. Conf. on Comp. Mech., 2019.
  * [29] C. J. Petrie.  The rheology of fibre suspensions.  J. Non-Newton. Fluid Mech., 87(2):369–402, 1999.
  * [30] B. Rodenborn, C.-H. Chen, H. L. Swinney, B. Liu, and H. Zhang.  Propulsion of microorganisms by a helical flagellum.  Proc. Natl. Acad. Sci., 110(5):E338–E347, 2013.
  * [31] M. J. Shelley and T. Ueda.  The Stokesian hydrodynamics of flexing, stretching filaments.  Phys. D, 146(1):221–245, 2000.
  * [32] C. Siewert, R. Kunnen, M. Meinke, and W. Schröder.  Orientation statistics and settling velocity of ellipsoids in decaying turbulence.  Atmospheric research, 142:45–56, 2014.
  * [33] S. E. Spagnolie and E. Lauga.  Comparative hydrodynamics of bacterial polymorphism.  Phys. Rev. Lett., 106(5):058103, 2011.
  * [34] B. Tapley, E. Celledoni, B. Owren, and H. I. Andersson.  A novel approach to rigid spheroid models in viscous flows using operator splitting methods.  Numer. Algorithms, pages 1–19, 2019.
  * [35] B. K. Tapley, H. I. Andersson, E. Celledoni, and B. Owren.  Computational methods for tracking inertial particles in discrete incompressible flows.  arXiv preprint arXiv:1907.11936, 2019.
  * [36] A.-K. Tornberg and K. Gustavsson.  A numerical method for simulations of rigid fiber suspensions.  J. Comput. Phys., 215(1):172–196, 2006.
  * [37] A.-K. Tornberg and M. J. Shelley.  Simulating the dynamics and interactions of flexible fibers in Stokes flows.  J. Comput. Phys., 196(1):8–40, 2004.
  * [38] L. N. Trefethen and J. Weideman.  The exponentially convergent trapezoidal rule.  siam REVIEW, 56(3):385–458, 2014.
  * [39] S. Twomey.  On the numerical solution of fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature.  Journal of the ACM (JACM), 10(1):97–101, 1963.
  * [40] B. J. Walker, M. P. Curtis, K. Ishimoto, and E. A. Gaffney.  A regularised slender-body theory of non-uniform filaments.  Journal of Fluid Mechanics, 899:A3, 2020.

##  Appendix A Dissipation matrix of a prolate spheroid

The non-dimensionalized body frame resistance tensor $R_{1}$ for a spheroid
with aspect ratio $\lambda$ was derived by Oberbeck [27] and is given by

$R_{1}=16\pi\lambda~{}\mathrm{diag}\left(\frac{1}{\chi_{0}+\alpha_{0}},\frac{1}{\chi_{0}+\beta_{0}},\frac{1}{\chi_{0}+\lambda^{2}\gamma_{0}}\right).$
(90)

The constants $\chi_{0}$, $\alpha_{0}$, $\beta_{0}$ and $\gamma_{0}$ were
calculated by Siewert [32] and are presented for a prolate ($\lambda>1$)
spheroid

$\displaystyle\chi_{0}=\frac{-\kappa_{0}\lambda}{\sqrt{\lambda^{2}-1}},$ (91)
$\displaystyle\alpha_{0}=\beta_{0}=\frac{\lambda^{2}}{\lambda^{2}-1}+\frac{\lambda\kappa_{0}}{2(\lambda^{2}-1)^{3/2}},$
(92)
$\displaystyle\gamma_{0}=\frac{-2}{\lambda^{2}-1}-\frac{\lambda\kappa_{0}}{(\lambda^{2}-1)^{3/2}},$
(93)
$\displaystyle\kappa_{0}=\ln\left(\frac{\lambda-\sqrt{\lambda^{2}-1}}{\lambda+\sqrt{\lambda^{2}-1}}\right).$
(94)

The torques $\mathbf{N}=(N_{x},N_{y},N_{z})^{\mathrm{T}}$ that describe the
rotational forces acting on an ellipsoid in creeping Stokes flow in the body
frame were calculated by Jeffery [14] and are presented in their non-
dimensional form with zero background flow

$\displaystyle N_{x}=$
$\displaystyle-\frac{16\pi\lambda}{3(\beta_{0}+\lambda^{2}\gamma_{0})}\left[(1+\lambda^{2})\omega_{x}\right],$
(95) $\displaystyle N_{y}=$
$\displaystyle-\frac{16\pi\lambda}{3(\alpha_{0}+\lambda^{2}\gamma_{0})}\left[(1+\lambda^{2})\omega_{y}\right],$
(96) $\displaystyle N_{z}=$
$\displaystyle-\frac{32\pi\lambda}{3(\alpha_{0}+\beta_{0})}\omega_{z}.$ (97)

Here $\boldsymbol{\omega}=(\omega_{x},\omega_{y},\omega_{z})^{\mathrm{T}}$ is
the body frame angular velocity, which is related to body frame angular
momentum by $\bm{m}=J\boldsymbol{\omega}$. Taking derivatives of $\mathbf{N}$
with respect to $\bm{m}$ gives for the rotational dissipation matrix

$R_{2}=-\frac{16\lambda}{3}\mathrm{diag}\left(\frac{(1+\lambda^{2})}{(\beta_{0}+\lambda^{2}\gamma_{0})},\frac{(1+\lambda^{2})}{(\alpha_{0}+\lambda^{2}\gamma_{0})},\frac{2}{(\alpha_{0}+\beta_{0})}\right)J^{-1}.$
(98)

The full dissipation matrix used for the calculation in Figure 10 is given by

$A_{sph}=\left(\begin{array}[]{cc}R_{1}&0\\\ 0&R_{2}\\\ \end{array}\right).$
(99)

##  Appendix B Regularizing effect of rigid body integration

Here we consider the regularizing effect of a linear functional of the form
(37) on numerical solutions to the first-kind equation (35) when the fiber
radius is constant. We consider specifically $M={\bf I}$.

First note that when the fiber radius is constant, the integral kernel of (35)
is symmetric. Therefore, if $\bm{f}(s)$ is a solution to equation (35) then so
is $\bm{f}(s)+\bm{f}_{o}(s)$, where $\bm{f}_{o}(s)=-\bm{f}_{o}(-s)$ is any odd
function. However, this nonuniqueness is not an issue if we are only
interested in the integral of $\bm{f}(s)+\bm{f}_{o}(s)$ over the same interval
as such an odd function would vanish.

This can further be illustrated using the singular value expansion. Let
$\bm{u}_{i}$ and $\bm{v}_{i}$ be the (orthogonal) left and right singular
vectors and $\sigma_{i}$ be the singular values of $\mathbf{K}$. It is
possible to express the integral operator as

$\mathbf{K}[f]=\sum_{i=1}^{\infty}\sigma_{i}(\bm{f},\bm{v}_{i})\bm{u}_{i},$
(100)

where $(\cdot,\cdot)$ is the inner product on $L^{2}([-L,L],\mathbb{R}^{3})$.
Similarly, we can expand the right hand side of the integral equation (35) in
terms of the basis $\bm{u}_{i}$:

$\bm{y}(s)=\sum_{i=1}^{\infty}y_{i}\bm{u}_{i}.$ (101)

From equation (35), we now obtain the relations

$\sigma_{i}f_{i}=y_{i},\quad\forall i.$ (102)

We see here that we can only solve for $f_{i}$ only if $\sigma_{i}\neq 0$;
otherwise the operator is singular. Our objective is to compute the integral
of the solution $\bm{f}(s)$, which is a linear functional applied to
$\bm{f}(s)$ that we denote by $\phi$. By the Reisz lemma this can be
represented by a function $\bar{\phi}\in L^{2}([-L,L],\mathbb{R}^{3})$ through
${\phi}[\bm{f}]=(\bar{\phi},\bm{f})$. Expanding this function in the basis
$\bm{v}_{i}$ gives

$\bar{\phi}=\sum_{i=1}^{\infty}\bar{\phi}_{i}\bm{v}_{i}$ (103)

where $\bar{\phi}_{i}=(\bar{\phi},\bm{v}_{i})$. We then obtain

${\phi}(\bm{f})=\sum_{i=1}^{\infty}\frac{\bar{\phi}_{i}y_{i}}{\sigma_{i}}.$
(104)

Since $\sum_{i=1}^{\infty}y_{i}<\infty$ and
$\sum_{i=1}^{\infty}\bar{\phi}_{i}<\infty$, these coefficients tend to zero
and can therefore be approximated by some finite truncation as long as
$\sigma_{i}$ doesn’t decay too rapidly. Therefore the integration step has a
regularizing effect on the ill-conditioned first-kind equation (35); however,
we are still not guaranteed that a solution exists in the first place.

Note that we can perform a similar analysis for the second-kind equation (36)
to obtain

${\phi}(\bm{f})=\sum_{i=1}^{\infty}\frac{\bar{\phi}_{i}y_{i}}{\alpha+\sigma_{i}}.$
(105)

In this case, even if $\sigma_{i}$ rapidly decays to zero, ${\phi}(\bm{f})$
exists and can be approximated by a finite sum.