Slender body theories for rotating filaments

Abstract

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

1. Introduction

Whether an organism has one or millions of cells, slender filaments play a vital role in its motion and maintenance. For instance, these types of filaments make up the flagellar appendages that bacteria use to propel themselves in run and tumble motions (Berg 2004; Lauga & Powers 2009). In animal cells, slender filaments are a key component in the cell cytoskeleton, which provides a scaffolding for the cell (Alberts et al. 2002), and adapts to promote cellular division, migration and resistance to deformation (Eghiaian, Rigato & Scheuring 2015).

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

Recent work has focused on placing these theories back in the context of three-dimensional well-posed partial differential equations (PDEs) (Koens & Lauga 2018; Mori, Ohm & Spirn 2020a,b; Mori & Ohm 2021). In particular, it has been shown that SBT can be derived from a three-dimensional boundary integral equation (Koens & Lauga 2018; Garg & Kumar 2022), and that its solution is an $O(\epsilon )$ approximation to a Stokes PDE with mixed Dirichlet–Neumann boundary data and a boundary condition that the fibre maintain the integrity of its cross-section (Mori et al. 2020a,b; Mori & Ohm 2021). A second class of recent work has developed SBTs for different geometries, including ribbons (Koens & Lauga 2016) and bodies with non-circular cross-section (Borker & Koch 2019).

Despite the breadth of literature on SBT and the novel theories for complex geometries, there remains an open problem which is simple to formulate, but difficult to solve, having to do with the rotation of a filament. Consider a filament with one clamped end and one free end, where a motor spins the clamped end with frequency $\omega$. At some critical frequency $\omega _c$, the applied torque overcomes the bending moment, the steady twirling state becomes unstable, and large-scale filament deformations are introduced. This instability, which has been studied extensively in Wolgemuth, Powers & Goldstein (2000), Lim & Peskin (2004), Wada & Netz (2006), Lee et al. (2014) and Maxian et al. (2022), has applications in the motion of bacterial flagella (Chen & Berg 2000; Yuan et al. 2010) and supercoiling during DNA transcription (Liu & Wang 1987). It has been shown that the critical frequency $\omega _c$ scales like $1/\epsilon ^2$, and thus that the torque required to induce it is approximately independent of the fibre aspect ratio (Powers 2010, equation (62)). Because of this, fibres being driven at the critical frequency have angular velocities which are sufficiently large as to also affect the translational dynamics. But as of yet, we are not aware of any SBT that can consistently account for the full hydrodynamics of a slender, translating and rotating filament, and in particular for the coupling between the two. In fact, we recently showed that the failure to account for rotation–translation coupling causes the SBT-based analysis of Wolgemuth et al. (2000) to overestimate $\omega _c$ by as much as 20 % (Maxian et al. 2022, § 5.3).

To understand the reason for the failures of previous SBTs, in § 2 we present a survey of an old and new approach to SBT, and an explanation of the challenges involved when we try to extend these approaches to account for rotation. Fundamentally, the issues that we point out can be traced to the difference between the resistance problem, which seeks the force leading to a given velocity, and the mobility problem, which is the inverse of the resistance problem and has different asymptotic behaviour. The problem that we consider here is a mobility problem, since the torque driving the filament is constant while the angular velocity changes with $\epsilon$. It is therefore not necessarily surprising that prior SBTs, which were designed with the resistance problem in mind, cannot be applied to the spinning filament mobility problem.

One potential workaround to this is the recent approach of Koens (2022), which is designed to exactly solve the resistance (and therefore mobility) problem by solving a series of Fredholm integral equations of the second kind to find the surface traction. These integral equations are similar to that of SBT, but are instead forced by a function of the previous order traction. The approach, however, offers little insight into the local resistance behaviour under high rotation, and it is unclear how many iterations would be required to accurately capture the non-local translation–rotation coupling that is not correctly treated by SBT. For this reason, we seek an alternative approximation which takes advantage of the fibre slenderness to minimize the number of non-local integral evaluations and allow for dynamic simulations with many fibres.

In § 3 of this paper, we offer a line integral of regularized singularities as a potential alternative to traditional SBT. Mathematically, the approach of the original slender body pioneers (Batchelor 1970a; Keller & Rubinow 1976; Johnson 1980) was to distribute singularities along the centreline, then evaluate the flow from these singularities on a fibre cross-section. The flow is non-singular since it is evaluated a fibre radius away from the location of the singularities. In regularized singularities, we instead distribute the force and torque over a ‘cloud’ or ‘blob’ region of characteristic size ${\hat {a}} \sim a$, which is defined by some regularized delta function $\delta _{\hat {a}}$. Throughout this paper, $\delta _{\hat {a}}$ will be a surface $\delta$ function on a sphere of radius ${\hat {a}}$, but other approaches are also possible. In this case, the velocity field is non-singular because the force and torque are distributed over a region of finite size approximately equal to the fibre radius. We can then convolve the velocity field (or half the vorticity field) with the regularized delta function again to obtain the velocity (or angular velocity) of the fibre centreline. In this way, we derive fundamental regularized singularities for all components of the mobility operator, then perform asymptotics on line integrals of these regularized singularities. The sum of these components gives a model of the filament dynamics with the correct physical behaviour that can be used in practice without approximating precisely a specific geometry; after all, the geometry of an actin filament (Alberts et al. 2002) or flagella (Berg 2004, figure 5.3) is not perfectly spheroidal or even tubular, and so one could argue that the regularized singularity model is just as much of an approximation as SBT.

Representing filaments using regularized singularities is not a novel idea, as it was first employed for numerical purposes in the immersed boundary method of the early 1970s to model the flow patterns inside the heart (Peskin 1972). Later variations, including the force coupling method (Lomholt & Maxey 2003), method of regularized Stokeslets (Cortez 2001; Cortez, Fauci & Medovikov 2005) and Rotne–Prager–Yamakawa (RPY) tensor (Rotne & Prager 1969; Wajnryb et al. 2013), focused on choosing the regularized delta function such that the Stokes equations are solvable analytically. While not a unique or necessarily superior choice, we prefer to use the RPY kernel because of its simple far-field analytical form (as a Stokeslet plus a source doublet), its symmetric positive definite property, and because of its long history of use in polymer suspensions (e.g. Butler & Shaqfeh 2005; Beck & Shaqfeh 2006; Wang et al. 2012; Schoeller et al. 2021). The analysis here could be conducted for any regularization function where the mobility between two points can be determined analytically.

Recently, the accuracy of using regularized singularities to approximate slender translating filaments has come under scrutiny both from a numerical (Bringley & Peskin 2008) and analytical (Cortez & Nicholas 2012; Maxian, Mogilner & Donev 2021; Ohm 2021; Zhao & Koens 2021) perspective. Generally speaking, the recent analysis has affirmed (Cortez & Nicholas 2012; Maxian et al. 2021; Ohm 2021) that regularization methods can effectively match the SBT centreline and far-field fluid velocity with a judicious choice of regularization radius ${\hat {a}}$, while others (Ohm 2021) have shown that the global flow field near the slender body cannot be matched unless very specific conditions on the regularization function are met (Zhao & Koens 2021, § 4). While the latter restriction is somewhat disappointing, it has little effect in practice: since regularized singularities can match the SBT filament centreline velocity, the only impact of the flow near the filament being incorrect is in computing disturbance flows on other nearly touching filaments. These disturbance flows change sharply over distances of $O(\epsilon )$, which makes a one-dimensional SBT approach for nearly touching fibres questionable if not implausible. One of the goals of this paper is therefore to determine how the centreline motion induced by RPY singularities compares with that of SBT when rotation is considered. For instance, does it reproduce the relationship between torque and rotational velocity for a straight cylinder? And what about the feedback of rotational velocity on translational motion for curved filaments?

These are the main questions we try to answer in this paper. To do so, we first look at how far we can get with traditional SBT in § 2. We show that the singularity methods of Keller and Rubinow and Johnson yield well-defined mobilities for translation from force (translation–translation) and rotation from torque (rotation–rotation) terms, but that both leave an unknown $O(1)$ constant in the rotation–translation coupling term which is relevant in the mobility problem. Then, we perform asymptotic expansions on line integrals of regularized RPY singularities to derive a consistent SBT for filaments rotating and translating with speeds of arbitrary order. We show that the regularized singularity radius ${\hat {a}}$ can be chosen to match either the three-dimensional translation–translation or rotation–rotation mobility, but not both, which implies that the set of RPY singularities we consider cannot accurately model a three-dimensional cylinder. We also show that, unlike SBTs that are more faithful to the filament geometry, the regularized RPY singularity approach gives symmetric rotation–translation coupling terms, and naturally extends to the fibre endpoints.

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

2. Three-dimensional SBTs for twist

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

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

(2.1)\begin{equation} \hat{\boldsymbol{X}}(s,\theta)=\boldsymbol{X}(s)+a \rho(s) \left(\boldsymbol{e}_n(s) \cos{\theta}+\boldsymbol{e}_b(s) \sin{\theta}\right):=\boldsymbol{X}(s)+a \rho(s) \hat{\boldsymbol{r}}(s,\theta), \end{equation}

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

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

(2.2)\begin{equation} \hat{\boldsymbol{U}}(s,\theta) = \boldsymbol{U}(s) + a \varPsi^\parallel(s) (\boldsymbol{\tau}(s) \times \hat{\boldsymbol{r}}(s,\theta)). \end{equation}

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

(2.3)\begin{equation} \boldsymbol{u}(\hat{\boldsymbol{X}}(s,\theta))=\hat{\boldsymbol{U}}(s,\theta). \end{equation}

2.1. Singularity representations

We first try to apply an older SBT approach to the spinning filament problem: the singularity approach of Johnson (1980), and later Götz (2001). The approach of Johnson is to look for a series of singularities along the fibre centreline to represent the flow field globally (everywhere outside the fibre). The singularities are chosen so that the boundary condition (2.3) is satisfied with an error $O(\epsilon )$.

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

(2.4)\begin{equation} \boldsymbol{u}^{(u)}(\boldsymbol{x}) = \frac{1}{8{\rm \pi}\mu} \int_0^L \left(\mathbb{S}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)+\frac{\left(a \rho(s)\right)^2}{2}\mathbb{D} \left(\boldsymbol{x},\boldsymbol{X}(s)\right)\right)\,\boldsymbol{f}(s) \, {\rm d}s \end{equation}

where

(2.5a,b)\begin{equation} \mathbb{S}\left(\boldsymbol{x},\boldsymbol{X}(s)\right) = \left(\frac{\boldsymbol{I}}{r}+\frac{\boldsymbol{r}\boldsymbol{r}^T}{r^3}\right), \quad \mathbb{D}\left(\boldsymbol{x},\boldsymbol{X}(s)\right) = \left(\frac{\boldsymbol{I}}{r^3}-3\frac{\boldsymbol{r}\boldsymbol{r}^T}{r^5}\right), \end{equation}

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

(2.6)$$\begin{gather} 8 {\rm \pi}\mu \boldsymbol{u}^{(u)}(\hat{\boldsymbol{X}}(s,\theta)) =\boldsymbol{U}_{f}^{(LD)}(s)+\boldsymbol{U}_{f}^{{(FP)}} + O(\epsilon), \end{gather}$$

(2.7)$$\begin{gather}\boldsymbol{U}_{f}^{(LD)}(s)=\ln{\left(\frac{4s(L-s)}{a^2}\right)} \left(\boldsymbol{I}+\boldsymbol{\tau}(s) \boldsymbol{\tau}(s)\right) \,\boldsymbol{f}(s) + \left(\boldsymbol{I}-3\boldsymbol{\tau}(s) \boldsymbol{\tau}(s)\right) \,\boldsymbol{f}(s), \end{gather}$$

(2.8)$$\begin{gather}\boldsymbol{U}_{f}^{{(FP)}}(s) = \int_0^L \left(\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(t)\right)\,\boldsymbol{f}(t)-\left( \frac{\boldsymbol{I}+\boldsymbol{\tau}(s) \boldsymbol{\tau}(s)}{|t-s|}\right)\,\boldsymbol{f}(s)\right), \end{gather}$$

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

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

(2.9)\begin{equation} n^\parallel = 4 {\rm \pi}\mu a^2 \varPsi^\parallel, \end{equation}

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

(2.10a,b)\begin{equation} 8 {\rm \pi}\mu \boldsymbol{u}^{(\varPsi)}\left(\boldsymbol{x}\right)=\int_0^L \mathbb{R}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)n^\parallel(s) \, {\rm d}s, \quad \mathbb{R}\left(\boldsymbol{x},\boldsymbol{X}(s)\right)= \frac{\boldsymbol{\tau}(s) \times \boldsymbol{r}}{r^3}. \end{equation}

Asymptotically, this velocity is equal to

(2.11)\begin{align} 8{\rm \pi} \mu \boldsymbol{u}^{(\varPsi)}(\hat{\boldsymbol{X}}\left(s,\theta\right)) &\approx \frac{2 n^\parallel}{a}\left(\boldsymbol{\tau} \times \hat{\boldsymbol{r}}\right) +\frac{n^\parallel}{2} \left(\ln{\left(\frac{4s(L-s)}{ a^2}\right)}-2\right) \left(\boldsymbol{\tau} \times {\partial s}{\boldsymbol{\tau}}\right)+{\boldsymbol{U}}^{{(FP)}} \nonumber\\ &\quad +n^\parallel \kappa \cos{\theta}\left(\boldsymbol{\tau} \times \hat{\boldsymbol{r}}\right) + O\left(\epsilon \log{\epsilon}\right), \end{align}

where

(2.12)\begin{equation} {\boldsymbol{U}}^{{(FP)}}(s)=\int_0^L \left(\mathbb{R}\left(\boldsymbol{X}(s),\boldsymbol{X}(t)\right)n^\parallel(t) - \frac{1}{2}\left(\frac{\boldsymbol{\tau}(s)\times{\partial s}{\boldsymbol{\tau}}(s) }{|t-s|}\right)n^\parallel(s)\right)\end{equation}

on the fibre surface. Substituting (2.9) into (2.11), we see that the first term is exactly the angular velocity of (2.2), while the second and third terms on the first line of (2.11) give an additional constant velocity on the cross-section (the rotation–translation coupling term). However, we still have an angular-dependent term on the second line of (2.11): since $\boldsymbol {\tau } \times \hat {\boldsymbol {r}}=\boldsymbol {e}_b \cos {\theta } - \boldsymbol {e}_n \sin {\theta }$, the term in the second line has the angle dependence

(2.13)\begin{equation} n^\parallel \kappa \cos{\theta}(\boldsymbol{\tau} \times \hat{\boldsymbol{r}})= \frac{n^\parallel \kappa}{2} \left( \boldsymbol{e}_b +\boldsymbol{e}_b \cos{2\theta} -\boldsymbol{e}_n \sin{2\theta}\right). \end{equation}

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

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

2.2. Boundary integral representations

Recently, effort has been made to place the asymptotic theories of Keller and Rubinow and Johnson on a more rigorous three-dimensional footing. Foremost among these are Mori et al. (2020a,b) and Mori & Ohm (2021), who show that SBT is an $O(\epsilon )$ approximation of a three-dimensional PDE with non-standard boundary conditions, and Koens & Lauga (2018), who show that the translational SBT of Johnson and Keller and Rubinow can be obtained from asymptotic expansion of the single layer boundary integral equation (Pozrikidis et al. 1992, § 4.1)

(2.14) \begin{equation} 8 {\rm \pi}\mu \, \hat{\boldsymbol{U}}(s,\theta) = \int_0^{2{\rm \pi}} \int_0^L \mathbb{S}{\hat{\boldsymbol{X}}(s,\theta),\hat{\boldsymbol{X}} (s',\theta')} \hat{\boldsymbol{f}}(s',\theta'){\rm d}s' \,{\rm d}\theta'. \end{equation}

Note that Koens and Lauga begin with the full boundary integral representation, but the single layer representation is sufficient as long as there is no change in the fibre volume. However, it is important to point out that $\hat {\boldsymbol {f}}$ in (2.14) is not in general the force density on the fibre surface. Notable exceptions are a filament undergoing rigid body motion or a ‘filament’ whose interior is a fluid of viscosity equal to that of the fibre exterior.

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

(2.15)\begin{equation} \boldsymbol{f}(s) = \int_0^{2 {\rm \pi}}\hat{\boldsymbol{f}}\left(s,\theta\right) {\rm d}\theta \end{equation}

by the classical SBT (2.6), and that, for a straight cylinder, the torque-angular velocity relationship (2.9) is recovered (Koens & Lauga 2018, equation (4.11)) for the total parallel torque on a cross-section,

(2.16)\begin{equation} n^\parallel(s) = \boldsymbol{\tau}(s) \boldsymbol{\cdot} \int_0^{2{\rm \pi}} a\rho(s) (\hat{\boldsymbol{r}}(s,\theta) \times \hat{\boldsymbol{f}}(s,\theta)) {\rm d}\theta. \end{equation}

Here we briefly discuss how the approach of Koens and Lauga might be extended to account for rotation in curved filaments. The main issue is that, if the boundary integral equation (2.14) is only expanded to leading order (i.e. with an error of $O(\epsilon )$), it is impossible to make observations about the rotational velocity, since the velocity on the cross-section induced by an angular velocity $\varPsi ^\parallel$ is $\sim \epsilon \varPsi ^\parallel$. Thus, more terms are necessary in the asymptotics. Unfortunately, because the expansions must be conducted around both $\theta$ and $s$ in (2.14), even expanding a single layer formulation to the next order is complex, as can be seen in the results of slender phoretic theory (Katsamba, Michelin & Montenegro-Johnson 2020), which completely extends the approach of Koens and Lauga to next order for the case of a body which swims via a fluid slip along a self-generated surface chemical concentration gradient.

To get a flavour for the challenges involved, we consider an asymptotic expansion of the isotropic part of the single-layer potential (cf. Koens & Lauga 2018, (5.1))

(2.17)\begin{equation} \boldsymbol{U}_{i1}(s,\theta):=\int_{0}^{2{\rm \pi}} \int_{0}^L \frac{\boldsymbol{f}(s',\theta')}{R(\hat{\boldsymbol{X}}(s,\theta),\hat{\boldsymbol{X}}(s',\theta'))} {\rm d}\theta' \, {\rm d}s', \end{equation}

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

(2.18)\begin{equation} \hat{\boldsymbol{f}}(s,\theta) = \frac{1}{2{\rm \pi}}\,\boldsymbol{f}(s) + \frac{1}{2{\rm \pi}}\left(\,\boldsymbol{f}_{c,1}(s) \cos{\theta } + \boldsymbol{f}_{s,1}(s) \sin{\theta }\right), \end{equation}

which we substitute into (A5). After integrating over $\theta '$, the final inner expansion for the surface velocity is

(2.19) \begin{align} \boldsymbol{U}_{i1}(s,\theta) &=\ln{\left(\frac{4s(L-s)}{a^2}\right)}\,\boldsymbol{f}(s)+(L-2s){\partial s}{\boldsymbol{f}}(s)+\boldsymbol{f}_{c,1}(s) \cos{\theta}\nonumber\\ &\quad -\frac{a \kappa}{2} \left(\,\boldsymbol{f}(s)\cos{\theta} + \left(\frac{1}{4}\cos{2\theta} + \frac{1}{2}\left(\ln{\left(\frac{4s(L-s)}{a^2 }\right)}-2\right)\right)\,\boldsymbol{f}_{c,1}(s)\right.\nonumber\\ &\quad \left.-\,\frac{1}{2}\,\boldsymbol{f}_{s,1}(s)\cos{\theta}\sin{\theta}\right)+O(\epsilon^2). \end{align}

We use the variable $\theta$ here for notational convenience; more precisely, $\theta$ should be measured relative to the curve torsion angle $\theta _i$(s), which is defined in Appendix A. That is, $\theta$ in (2.18) and (2.19) stands for $\theta -\theta _i(s)$.

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

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

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

3. Slender body theory from regularized singularities

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

An alternative approach, which we consider in this section, is to relax the fidelity to a true cylindrical geometry and instead represent the cylinder using regularized singularities distributed along the centreline. The idea of the regularized singularity approach is to replace the delta function in the force and torque with a regularized delta function $\delta _{\hat {a}}$, where ${\hat {a}}$ is the ‘regularization radius’ that is chosen to model the three-dimensional problem. Our choice of kernel is the RPY tensor (Rotne & Prager 1969; Wajnryb et al. 2013), for which $\delta _{\hat {a}}(\boldsymbol {r})=\delta (r-{\hat {a}})/(4{\rm \pi} {\hat {a}}^2)$, i.e. the regularized delta function is a surface delta function on a sphere of radius ${\hat {a}}$. This choice is more faithful to the original cylindrical geometry than delta functions with infinite support (Cortez 2001; Maxey & Patel 2001; Cortez et al. 2005) or those whose support is tied to a numerical fluid grid (Peskin 2002). In particular, since each regularized singularity is a sphere, it is not difficult to imagine this resulting in a cylinder of constant radius ${\hat {a}}$ along a length $L$, with hemispherical caps at each end. Still, since the regularized delta function is the same in the axial and radial directions, the spheres are not rings as we would need to make a proper cylinder. Also, the RPY tensor itself is an approximation to the dynamics of spheres, since it does not include stresslet terms required to keep the spheres rigid in flow (Fiore & Swan 2018), even though stresslets are of the same multipole order as torques, which are included. Because of these approximations, the intuitive picture of a cylinder composed of surface spheres is not quite the correct one, and we will require ${\hat {a}} > a$ for our RPY formulae to match those of SBT.

To obtain the centreline velocity of a filament made of regularized singularities, we first solve for the fluid velocity $\boldsymbol {u}(x)$ and corresponding pressure ${\rm \pi} ( \boldsymbol {x} )$ by (It is important to emphasize again that $\boldsymbol {u}(\boldsymbol {x})$ and ${\rm \pi} (\boldsymbol {x})$ are expected to be good approximations of the true fluid velocity and pressure only sufficiently far from the fibre centreline, but not $O(\epsilon )$ away from it.) solving the Stokes equations (Maxian et al. 2022, § 2.1)

(3.1)\begin{equation} \boldsymbol{\nabla}{\rm \pi}(\boldsymbol{x})=\mu \nabla^{2}\boldsymbol{u}(\boldsymbol{x})+ \int_0^L \left(\left(\,\boldsymbol{f}(s) + \frac{\boldsymbol{\nabla}}{2} \times n^\parallel(s)\boldsymbol{\tau}(s)\right)\delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{X}(s)\right)\right) {\rm d}s, \end{equation}

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

(3.2)$$\begin{gather} \boldsymbol{U} \left(s\right)= \int \boldsymbol{u}(\boldsymbol{x}) \, \delta_{{\hat{a}}}\left(\boldsymbol{X}(s)-\boldsymbol{x}\right) {\rm d}\boldsymbol{x} , \end{gather}$$

(3.3)$$\begin{gather}\varPsi^\parallel \left(s\right)=\boldsymbol{\tau}(s) \boldsymbol{\cdot} \int \delta_{{\hat{a}}}\left(\boldsymbol{X}(s)-\boldsymbol{x}\right)\frac{\boldsymbol{\nabla}}{2}\times\boldsymbol{u}\left(\boldsymbol{x}\right) {\rm d}\boldsymbol{x}. \end{gather}$$

This yields a mobility that is guaranteed to be symmetric positive definite (Kallemov et al. 2016, footnote 3). Here we only consider the parallel component of rotational velocity, since, for unshearable filaments, the other components can be derived from the translational velocity $\boldsymbol {U}(s)$ (Maxian et al. 2022, § 2).

Using superposition, the translational and rotational velocity of the fibre centreline can be defined as integrals of the fundamental regularized singularities

(3.4)\begin{gather} \left.\begin{gathered} \boldsymbol{U}(s) = \int_0^L \left(\mathbb{M}_{tt}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^\prime)\right)\,\boldsymbol{f}(s^\prime)+\mathbb{M}_{tr}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^\prime)\right)n^ \parallel(s')\boldsymbol{\tau}(s')\right) {\rm d}s^\prime \\ \varPsi^\parallel(s) = \boldsymbol{\tau}(s) \boldsymbol{\cdot} \int_0^L \left(\mathbb{M}_{rt}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^\prime)\right)\,\boldsymbol{f}(s^\prime)+\mathbb{M}_{rr}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^\prime)\right)n^ \parallel(s')\boldsymbol{\tau}(s')\right) {\rm d}s^\prime \end{gathered}\right\} \end{gather}

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

(3.5)\begin{equation} \left.\begin{gathered} \mathbb{M}_{tt}\left(\boldsymbol{x},\boldsymbol{y}\right) = \int \delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{z}\right) \int \mathbb{S}{\boldsymbol{z},\boldsymbol{w}}\delta_{{\hat{a}}}\left(\boldsymbol{w}-\boldsymbol{y}\right) {\rm d}\boldsymbol{w} \, {\rm d}\boldsymbol{z},\\ \mathbb{M}_{tr}\left(\boldsymbol{x},\boldsymbol{y}\right) = \int \delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{z}\right) \int \mathbb{S}{\boldsymbol{z},\boldsymbol{w}} \frac{\boldsymbol{\nabla}_{\boldsymbol{w}}}{2} \times\delta_{{\hat{a}}}\left(\boldsymbol{w}-\boldsymbol{y}\right) {\rm d}\boldsymbol{w} \, {\rm d}\boldsymbol{z},\\ \mathbb{M}_{rt}\left(\boldsymbol{x},\boldsymbol{y}\right) = \int \delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{z}\right) \frac{\boldsymbol{\nabla}_{\boldsymbol{z}}}{2} \times \int \mathbb{S}{\boldsymbol{z},\boldsymbol{w}}\delta_{{\hat{a}}}\left(\boldsymbol{w}-\boldsymbol{y}\right) {\rm d}\boldsymbol{w} \,{\rm d}\boldsymbol{z}= \mathbb{M}_{tr}^T\left(\boldsymbol{y},\boldsymbol{x}\right),\\ \mathbb{M}_{rr}\left(\boldsymbol{x},\boldsymbol{y}\right) = \int \delta_{{\hat{a}}}\left(\boldsymbol{x}-\boldsymbol{z}\right) \frac{\boldsymbol{\nabla}_{\boldsymbol{z}}}{2} \times \int \mathbb{S}{\boldsymbol{z},\boldsymbol{w}} \frac{\boldsymbol{\nabla}_{\boldsymbol{w}}}{2} \times\delta_{{\hat{a}}}\left(\boldsymbol{w}-\boldsymbol{y}\right) {\rm d}\boldsymbol{w} \, {\rm d}\boldsymbol{z}. \end{gathered}\right\} \end{equation}

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

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

In this section, we therefore perform asymptotic expansions on the integrals (3.4) to obtain an SBT for a twisting filament. Our strategy is standard matched asymptotics and similar to the approaches of Johnson (1980) and Götz (2001) for SBT. Because we use regularized singularities, however, all of the integrals take place on the fibre centreline, and so the asymptotics are in one dimension, rather than on the fibre cross-section (two dimensions). It therefore becomes much simpler to obtain expressions for the velocity at the endpoints (see Appendix B), where we only have to redefine the domain of integration in the inner expansion. Thus, a physical endpoint regularization, which is an ingredient of any SBT-based numerical method (Maxian et al. 2021, § 2.1), comes with our choice to use regularized singularities as an approximation to the three-dimensional fibre geometry.

For each integral in (3.4), we first compute an outer expansion by considering the region where $|s-s'|$ is $O(1)$. Then in the inner expansion, we consider the part of the integrals where $|s-s'|$ is $O({\hat {a}})$. In this case, we follow (Götz 2001) and introduce the rescaled variable

(3.6)\begin{equation} \xi = \frac{s'-s}{{\hat{a}}}, \end{equation}

so that $\xi$ is $O(1)$ in the inner expansion, and the domain of $\xi$ is $[-s/{\hat {a}},(L-s)/{\hat {a}}]$. With this definition, it is straightforward to integrate over $\xi$ and obtain an inner expansion, which is done in two pieces because the RPY tensor changes at $s=2{\hat {a}}$ (the case of overlapping spheres). The breakup of the integral into two pieces is reminiscent of the classic approach of Lighthill (1976) (see also Lauga & Powers (2009), figure 6), who represented the flow due to a translating filament as the sum of an integral of Stokeslets (at distances greater than some intermediate length $q \gg a$) and another integral of Stokeslets and doublets (at distances smaller than $q$). The RPY approach is fundamentally different, however, because here the RPY regularization imposes the splitting of the integral at the distance $2{\hat {a}}$. For translation–translation, the doublet term can still make a contribution to the flow field on this length scale; its precise $O(1)$ contribution is computed in (B4). Adding the inner and outer solutions together and subtracting the common part then gives the final matched asymptotic solution. Unless otherwise specified, all expansions will be carried out to $O(\hat {\epsilon })$, i.e. in this section $\approx$ denotes equality up to $O(\hat {\epsilon })$.

3.1. Translation from force

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

(3.7)\begin{equation} \boldsymbol{U}_{f}=\int_0^L \mathbb{M}_{tt}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^\prime)\right)\,\boldsymbol{f}(s^\prime) \, {\rm d}s'. \end{equation}

The following derivation is a repeat of that given in Maxian et al. (2021, Appendix A), but it is included here for completeness, and because the set of steps for the other three blocks of the mobility operator is identical. When we substitute the definition of $\mathbb {M}_{tt}$ for RPY singularities from Wajnryb et al. (2013, (3.12)) into (3.7), we have the integral

(3.8)\begin{align} 8{\rm \pi} \mu\boldsymbol{U}_{f}\left(s\right) &= \int_{R> 2{\hat{a}}} \left(\mathbb{S}\left(\boldsymbol{X}(s),\boldsymbol{X}(s')\right)+\frac{2{\hat{a}}^2}{3}\mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s')\right)\right)\,\boldsymbol{f}(s')\, \,{\rm d}s' \nonumber\\ &\quad +\int_{R \leq 2{\hat{a}}} \left(\left(\frac{4}{3{\hat{a}}}-\frac{3R\left(s'\right)}{8{\hat{a}}^2}\right)\boldsymbol{I}+\frac{1}{8{\hat{a}}^2R\left(s'\right)} \left(\boldsymbol{R}\boldsymbol{R}\right)\left(s'\right)\right)\,\boldsymbol{f}\left(s'\right) {\rm d}s', \end{align}

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

In the outer expansion, we consider the part of the integral (3.8) where $|s-s'|$ is $O(1)$. In this case, the doublet term in (3.8) is insignificant and we obtain the outer velocity by integrating the Stokeslet over the fibre centreline,

(3.9)\begin{equation} 8{\rm \pi} \mu \boldsymbol{U}_{f}^{({outer})}(s) = \int_{|s-s'| > 2{\hat{a}}} \mathbb{S}{\boldsymbol{X}(s),\boldsymbol{X}(s')}\,\boldsymbol{f}\left(s'\right) {\rm d}s'. \end{equation}

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

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

(3.10)\begin{align} 8 {\rm \pi}\mu \boldsymbol{U}_{f}^{({inner})}(s) &= \left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\,\boldsymbol{f}(s) \ln{\left(\frac{(L-s)s}{4{\hat{a}}^2}\right)} \nonumber\\ &\quad + \left(\left(4-\frac{{\hat{a}}^2}{3s^2}-\frac{{\hat{a}}^2}{3(L-s)^2}\right)\boldsymbol{I} + \left(\frac{{\hat{a}}^2}{s^2}+\frac{{\hat{a}}^2}{(L-s)^2}\right) \boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\,\boldsymbol{f}(s), \end{align}

in the fibre interior. The complete formulae including endpoint modifications are given in (B7); the $O({\hat {a}}^2)$ terms are included in (3.10) because they are necessary to give $\mathcal {C}^2$ continuity with the endpoint formulae.

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

(3.11)\begin{equation} 8{\rm \pi} \mu \boldsymbol{U}_{f}^{({common})}(s) = \int_{|s-s'| > 2{\hat{a}}} \left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s-s'|}\right)\,\boldsymbol{f}(s) \, {\rm d}s'. \end{equation}

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

(3.12)\begin{equation} \boldsymbol{U}_{f}(s) = \boldsymbol{U}_{f}^{({inner})}(s) + \boldsymbol{U}_{f}^{({outer})}(s)-\boldsymbol{U}_{f}^{({common})}(s). \end{equation}

This can be written as

(3.13)\begin{align} \boldsymbol{U}_{f}(s) = \boldsymbol{U}_{f}^{({inner})}(s) + \frac{1}{8{\rm \pi}\mu}\int_{|s-s'|> 2{\hat{a}}} \left(\mathbb{S}{\boldsymbol{X}(s),\boldsymbol{X}(s')} \,\boldsymbol{f}\left(s'\right) - \left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s-s'|}\right)\,\boldsymbol{f}(s)\right) {\rm d}s', \end{align}

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

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

(3.14)\begin{align} 8 {\rm \pi}\mu \boldsymbol{U}_{f}(s) &= \left(\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)\right)\,\boldsymbol{f}(s) \ln{\left(\frac{(L-s)s}{4{\hat{a}}^2}\right)} +4\,\boldsymbol{f}(s) \nonumber\\ &\quad+ \int_{0}^L \left(\mathbb{S}{\boldsymbol{X}(s),\boldsymbol{X}(s')} \,\boldsymbol{f}\left(s'\right) - \left(\frac{\boldsymbol{I}+\boldsymbol{\tau}(s)\boldsymbol{\tau}(s)}{|s-s'|}\right)\,\boldsymbol{f}(s)\right) {\rm d}s', \end{align}

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

3.2. Translation from torque

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

(3.15)\begin{equation} {\boldsymbol{U}}(s)=\int_0^L \mathbb{M}_{tr}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^\prime)\right)n^\parallel(s') \boldsymbol{\tau}(s') \, {\rm d}s'. \end{equation}

When we substitute the definition of $\mathbb {M}_{tr}$ for RPY singularities from Wajnryb et al. (2013, (3.16)) into (3.15), we have the integral

(3.16)\begin{align} 8{\rm \pi} \mu {\boldsymbol{U}}(s) &= \int_{R > 2{\hat{a}}} \frac{\boldsymbol{\tau}(s') \times \boldsymbol{R}(s')}{R(s')^3}n^\parallel(s') \, {\rm d}s'\nonumber\\ &\quad + \frac{1}{2{\hat{a}}^2}\int_{R \leq 2{\hat{a}}} \left(\frac{1}{{\hat{a}}}-\frac{3R(s')}{8{\hat{a}}^2}\right)\left(\boldsymbol{\tau}(s') \times \boldsymbol{R}(s')\right)n^\parallel(s') \,{\rm d}s', \end{align}

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

In the outer expansion, we consider the part of the integral (3.16) where $|s'-s|$ is $O(1)$. In this case the correction term is insignificant and we get the outer velocity

(3.17)\begin{equation} 8{\rm \pi} \mu {\boldsymbol{U}}^{{(outer)}}(s) = \int_{|s-s'| > 2{\hat{a}}} \frac{\boldsymbol{\tau}(s') \times \boldsymbol{R}(s')}{R(s')^3}n^\parallel(s') \, {\rm d}s'. \end{equation}

In the inner expansion, we use the inner asymptotics (B1) to write the inner expansion as

(3.18)\begin{equation} 8{\rm \pi}\mu {\boldsymbol{U}}^{{(inner)}}(s) = \frac{1}{2}\left(\boldsymbol{\tau}(s) \times {\partial s}{\boldsymbol{\tau}}(s)\right)\left(\ln{\left(\frac{s(L-s)}{4{\hat{a}}^2}\right)}+\frac{7}{6}\right) \end{equation}

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

The common part for translation–rotation coupling is the inner expansion written in terms of the outer variables,

(3.19)\begin{equation} 8{\rm \pi}\mu {\boldsymbol{U}}^{{(common)}}(s) =\int_{|s-s'| > 2{\hat{a}}} \left(\frac{\boldsymbol{\tau}(s) \times {\partial s}{\boldsymbol{\tau}}(s)}{2|s-s'|}\right)n^\parallel(s) \, {\rm d}s'. \end{equation}

To $O(\hat {\epsilon })$, the velocity due to torque $\boldsymbol {n}(s)=n^\parallel(s)\boldsymbol {\tau }(s)$ on the fibre centreline is given by

(3.20)$$\begin{gather} {\boldsymbol{U}}(s)= {\boldsymbol{U}}^{{(inner)}}(s)+{\boldsymbol{U}}^{{(outer)}}(s)-{\boldsymbol{U}}^{{(common)}}(s) \end{gather}$$

(3.21)$$\begin{gather}{\boldsymbol{U}}(s) = {\boldsymbol{U}}^{{(inner)}}(s)+ \frac{1}{8{\rm \pi}\mu} \int_{R > 2{\hat{a}}} \left(\frac{\boldsymbol{\tau}(s') \times \boldsymbol{R}(s')}{R(s')^3}n^\parallel(s')- \frac{\boldsymbol{\tau}(s) \times {\partial s}{\boldsymbol{\tau}}(s)}{2|s-s'|} n^\parallel(s)\right) {\rm d}s' . \end{gather}$$

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

(3.22)\begin{align} 8{\rm \pi} \mu {\boldsymbol{U}}(s) &=\frac{1}{2}\left(\boldsymbol{\tau}(s) \times {\partial s}{\boldsymbol{\tau}}(s)\right)\left(\ln{\left(\frac{s(L-s)}{4{\hat{a}}^2}\right)}+\frac{7}{6}\right)\nonumber\\ &\quad + \int_{0}^L \left(\left(\frac{\boldsymbol{\tau}(s') \times \boldsymbol{R}(s')}{R(s')^3}\right)n^\parallel(s')-\left(\frac{\boldsymbol{\tau}(s) \times {\partial s}{\boldsymbol{\tau}}(s)}{2|s-s'|}\right) n^\parallel(s) \right) {\rm d}s' , \end{align}

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

3.3. Rotation from force

We next calculate the parallel rotational velocity $ {\varPsi }^\parallel (s)$ due to the force density $\boldsymbol {f}(s)$ on the fibre centreline,

(3.23)\begin{equation} {\varPsi}^\parallel(s)=\boldsymbol{\tau}(s) \boldsymbol{\cdot} \int_0^L \mathbb{M}_{rt}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^\prime)\right)\,\boldsymbol{f}(s') \, {\rm d}s'. \end{equation}

When we substitute the definition of $\mathbb {M}_{tr}$ from (Wajnryb et al. 2013, (3.16)) into (3.23), we have

(3.24)\begin{align} 8{\rm \pi} \mu {\varPsi}^\parallel(s)&=\int_{R > 2{\hat{a}}} \frac{\left(\,\boldsymbol{f}(s') \times \boldsymbol{R}(s')\right)}{R(s')^3}\boldsymbol{\cdot} \boldsymbol{\tau}(s) \, {\rm d}s'\nonumber\\ &\quad + \frac{1}{2{\hat{a}}^2}\int_{R < 2{\hat{a}}} \left(\frac{1}{{\hat{a}}}-\frac{3R(s')}{8{\hat{a}}^2}\right)\left(\,\boldsymbol{f}(s') \times \boldsymbol{R}(s')\right)\boldsymbol{\cdot} \boldsymbol{\tau}(s)\, {\rm d}s'. \end{align}

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

(3.25)$$\begin{gather} 8{\rm \pi} \mu {\varPsi}^\parallel(s) = \frac{1}{2}\left(\ln{\left(\dfrac{s(L-s)}{4{\hat{a}}^2}\right)}+ \dfrac{7}{6}\right)\left(\boldsymbol{\tau}(s) \times {\partial s}{\boldsymbol{\tau}}(s)\right) \boldsymbol{\cdot} \,\boldsymbol{f}(s)+{\varPsi}^{(FP)}(s) \end{gather}$$

(3.26)$$\begin{gather}{\varPsi}^{(FP)}(s) = \int_0^L \left(\frac{\boldsymbol{R}(s') \times \boldsymbol{\tau}(s)}{R(s')^3} \boldsymbol{\cdot} \,\boldsymbol{f}(s')-\frac{\boldsymbol{\tau}(s) \times {\partial s}{\boldsymbol{\tau}}(s)}{2|s-s'|}\boldsymbol{\cdot} \,\boldsymbol{f}(s)\right) {\rm d}s', \end{gather}$$

in the fibre interior. To obtain a formula accurate to $O(\hat {\epsilon })$ up to and including the fibre endpoints, we replace the local part in (3.25) with $ {\varPsi }^{\parallel, {inner}}$ defined in (B12).

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

(3.27)\begin{equation} \int_0^L {\varPsi}^\parallel(s) n^\parallel(s) \, {\rm d}s = \int_0^L {\boldsymbol{U}}(s) \boldsymbol{\cdot} \,\boldsymbol{f}(s) \, {\rm d}s. \end{equation}

The only non-trivial part of this calculation is the rotlet term in the finite part integrals, in which the order of integration in $s$ and $s'$ must be swapped.

3.4. Rotation from torque

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

(3.28)\begin{equation}{\boldsymbol{\varPsi}}(s)=\int_0^L \mathbb{M}_{rr}\left(\boldsymbol{X}(s),\boldsymbol{X}(s^\prime)\right)\boldsymbol{n}(s') \, {\rm d}s', \end{equation}

and specialize to the case of parallel velocity and torque later. We substitute the definition of $\mathbb {M}_{rr}$ from Wajnryb et al. (2013, (3.14)) into (3.28) to obtain the vector rotational velocity,

(3.29)\begin{gather} 8{\rm \pi} \mu {\boldsymbol{\varPsi}}(s) ={-}\frac{1}{2}\int_{R > 2{\hat{a}}} \mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s')\right) \boldsymbol{n}(s') \, {\rm d}s'\nonumber\\ \quad +\frac{1}{{\hat{a}}^3}\int_{R \leq 2{\hat{a}}} \left(\left(1-\frac{27R(s')}{32{\hat{a}}}\!+\frac{5 R(s')^3}{64 {\hat{a}}^3}\right)\boldsymbol{I} \!+\!\left(\frac{9}{32{\hat{a}} R(s')}-\frac{3R(s')}{64{\hat{a}}^3}\right)\left(\boldsymbol{R}\boldsymbol{R}\right)(s')\right)\boldsymbol{n}(s') \, {\rm d}s'. \end{gather}

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

In the outer expansion, we consider the part of the integral (3.29) where $|s'-s|$ is $O(1)$,

(3.30)\begin{equation} 8{\rm \pi} \mu {\boldsymbol{\varPsi}}^{{(outer)}}(s) ={-}\tfrac{1}{2}\int_{|s-s'| > 2{\hat{a}}} \mathbb{D}\left(\boldsymbol{X}(s),\boldsymbol{X}(s')\right) \boldsymbol{n}(s') \, {\rm d}s'. \end{equation}

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

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

(3.31)\begin{equation} 8{\rm \pi} \mu{\boldsymbol{\varPsi}}^{{(inner)}}(s)=\left(p_I(s) \boldsymbol{I}+p_\tau(s) \boldsymbol{\boldsymbol{\tau}}(s)\boldsymbol{\tau}(s)\right)\boldsymbol{n}(s), \end{equation}

where

(3.32)\begin{equation} \left.\begin{gathered} p_I (s)= \dfrac{1}{{\hat{a}}^2} \left(\dfrac{9}{8}+\dfrac{{\hat{a}}^2}{4}\left(\dfrac{1}{s^2}+\dfrac{1}{(L-s)^2}\right)\right)\\ p_\tau(s) = \dfrac{1}{{\hat{a}}^2}\left(\dfrac{9}{8}-\dfrac{3{\hat{a}}^2}{4}\left(\dfrac{1}{s^2}+\dfrac{1}{(L-s)^2}\right)\right) \end{gathered}\right\} \end{equation}

in the fibre interior. The functions $p_I$ and $p_\tau$ are different at the fibre endpoints, and are given in (B15). As in the translation case, the $O({\hat {a}}^2)$ terms in (3.31) are negligible in the fibre interior, but are necessary to give $\mathcal {C}^2$ continuity with the endpoint formulae in (B15).

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

(3.33)\begin{equation}{\boldsymbol{\varPsi}}= {\boldsymbol{\varPsi}}^{{(inner)}}. \end{equation}

Specializing now to the case of parallel angular velocity from parallel torque, $\boldsymbol {n}=n^\parallel \boldsymbol {\tau }$, and $ {\varPsi }^\parallel = {\boldsymbol {\varPsi }} \boldsymbol {\cdot } \boldsymbol {\tau }$, and dropping the negligible terms in the fibre interior from (3.31), we obtain the simple formula

(3.34)\begin{equation} 8 {\rm \pi}\mu {\varPsi}^\parallel(s) = \frac{9n^\parallel(s)}{4{\hat{a}}^2} \end{equation}

in the fibre interior.

3.5. Comparison with Keller–Rubinow–Johnson

We recall that the goal of our regularized singularity approach is to find ${\hat {a}}$ such that the RPY-based filament centreline velocity matches that of SBT. With this in mind, let us compare our RPY-based SBT with the more rigorous three-dimensional SBT results we have from § 2. We first compare the theories when we have a well-defined SBT, i.e. for translation–translation and rotation–rotation. Then, we use the form of the translation–rotation RPY coupling and our results in § 2 to postulate an SBT equation for translation–rotation coupling with a single unknown parameter. The parameter will be determined in the next section.

3.5.1. Translation–translation and rotation–rotation

All three-dimensional SBTs agree that the translational velocity in the case when $\varPsi ^\parallel =0$ is given by

(3.35)\begin{equation} 8 {\rm \pi}\mu \boldsymbol{U}_{f}=\ln{\left(\frac{4s(L-s)}{a ^2}\right)} \left(\boldsymbol{I}+\boldsymbol{\tau} \boldsymbol{\tau}\right) \,\boldsymbol{f} + \left(\boldsymbol{I}-3\boldsymbol{\tau} \boldsymbol{\tau}\right) \,\boldsymbol{f}+\boldsymbol{U}_{f}^{{(FP)}}+O(\epsilon), \end{equation}

while our RPY-based theory (3.14) gives

(3.36)\begin{equation} 8 {\rm \pi}\mu \boldsymbol{U}_{f}=\ln{\left(\frac{s(L-s)}{4{\hat{a}} ^2}\right)} \left(\boldsymbol{I}+\boldsymbol{\tau} \boldsymbol{\tau}\right) \,\boldsymbol{f} + 4\,\boldsymbol{f}+\boldsymbol{U}_{f}^{{(FP)}}+O(\epsilon). \end{equation}

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

(3.37)\begin{equation} 8 {\rm \pi}\mu {\varPsi}^\parallel{=}\frac{2 n^\parallel }{a^2}, \end{equation}

while our RPY-based theory gives (3.34).

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

(3.38)\begin{equation} {\hat{a}}={\hat{a}}_{tt} = a \frac{e^{3/2}}{4} \approx 1.12 a, \end{equation}

then the two translational theories match (Maxian et al. 2021, Appendix A), but for rotation we obtain

(3.39)\begin{equation} 8{\rm \pi}\mu {\varPsi}^\parallel{=}\frac{9n^\parallel(s)}{4{\hat{a}}_{tt}^2}=\frac{36n^\parallel(s)}{e^3 a^2} \approx 1.79 \frac{n^\parallel}{a^2}, \end{equation}

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

On the other hand, if we take

(3.40)\begin{equation} {\hat{a}}={\hat{a}}_{rr} = \sqrt{9/8} a \approx 1.06 a, \end{equation}

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

(3.41)\begin{align} 8 {\rm \pi}\mu \boldsymbol{U}&=\ln{\left(\frac{s(L-s)}{4{\hat{a}}_{rr} ^2}\right)} \left(\boldsymbol{I}+\boldsymbol{\tau} \boldsymbol{\tau}\right) \,\boldsymbol{f} + 4\,\boldsymbol{f}+\boldsymbol{U}_{f}^{{(FP)}}\nonumber\\ & \approx \ln{\left(\frac{4s(L-s)}{a ^2}\right)}\left(\boldsymbol{I}+\boldsymbol{\tau} \boldsymbol{\tau}\right) \,\boldsymbol{f} + \left(1.1\boldsymbol{I}-2.9\boldsymbol{\tau} \boldsymbol{\tau}\right) \,\boldsymbol{f} + \boldsymbol{U}_{f}^{{(FP)}}, \end{align}

which gives a small $O(1)$ error relative to translational SBT (3.35). The need to choose a different regularized singularity radius for the two different mobility components shows that, despite the intuitive picture of making a cylinder out of spheres, our choice of spherically symmetric regularized singularities is insufficient to exactly model a cylinder. It is clear, however, that the RPY approximation captures the correct form of all the relevant terms, but with $O(1)$ error in the coefficients.

3.5.2. Translation–rotation coupling

The more compelling question is whether we can use our findings to gain insights into the unknown $O(1)$ coefficient in translation–rotation SBT. Our RPY-based result is (rearranging (3.22) using properties of logs and rounding the $O(1)$ coefficient to two decimal places)

(3.42)\begin{equation} {8{\rm \pi} \mu {\boldsymbol{U}} \approx \left(\ln{\left(\frac{4 s(L-s)}{{\hat{a}}^2}\right)}-1.61\right)\left(\boldsymbol{\tau} \times {\partial s}{\boldsymbol{\tau}}\right)\frac{n^\parallel}{2}+{\boldsymbol{U}}^{{(FP)}},} \end{equation}

with the corresponding symmetric result for $ {\varPsi }^\parallel$ in (3.25). Since § 3.5.1 shows that the RPY models capture the essential physics of the interaction of the fibre with the fluid, but with different coefficients, we can postulate that the rotation–translation coupling in a true tubular geometry is also of the form (3.42), but with different coefficients. Furthermore, since there is an $O(1)$ error in the boundary condition (2.3) in the result (2.11) of Johnson, we can expect that any modifications to (2.11) will be $O(1)$ (not $\log {\epsilon }$). Combining these two assumptions, we postulate a rotation–translation coupling term of the form

(3.43)\begin{equation} 8{\rm \pi} \mu {\boldsymbol{U}} = \left(\ln{\left(\frac{4s(L-s)}{ a^2}\right)}-k\right) \left(\boldsymbol{\tau} \times {\partial s}{\boldsymbol{\tau}}\right)\frac{n^\parallel}{2}+{\boldsymbol{U}}^{{(FP)}}+O(\epsilon) \end{equation}

for some unknown constant $k$, with the corresponding symmetric form for $ {\varPsi }^\parallel$. In the next section, we fit the constant $k$ to three-dimensional boundary integral results. This gives us an empirical three-dimensional SBT for rotation–translation coupling to which we can compare the RPY theory. It is important to note that our work does not rigorously justify (3.43) for a specific three-dimensional fibre geometry.

4. Comparing with three-dimensional boundary integral

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

(4.1)\begin{equation} \rho(s) = \sqrt{1-s^2}. \end{equation}

The radius at each point $s$ is $a \rho (s)$, where we will vary $a$ to study the limit $\epsilon =a/L \rightarrow 0$. While the result (3.42) was derived by assuming a constant regularization radius ${\hat {a}}$, it is unchanged to leading order when considering a spatially varying radius function ${\hat {a}}(s)$, for which the polydisperse RPY kernel must be used (Zuk et al. 2014). The choice of an ellipsoidally decaying radius function removes any singular behaviour of (3.43) near the endpoints and speeds up the convergence of the boundary integral method (relative to a cylinder with spherical caps, for example).

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

(4.2)\begin{equation} 8 {\rm \pi}\mu \hat{\boldsymbol{U}}(s,\theta) = \int_{{-}1}^1 \int_{0}^{2{\rm \pi}} \mathbb{S}(\hat{\boldsymbol{X}}(s,\theta),\hat{\boldsymbol{X}} (s',\theta'))\, \hat{\boldsymbol{f}} (s',\theta') {\rm d}\theta' \,{\rm d}s' , \end{equation}

for $\hat {\boldsymbol {f}}(s,\theta )$, where $\hat {\boldsymbol {X}}$ and $\hat {\boldsymbol {U}}$ are defined in (2.1) and (2.2).

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

1. (i) Given the surface velocity $\hat {\boldsymbol {U}}(s,\theta )= \boldsymbol {U}(s) + a \varPsi ^\parallel (s) (\boldsymbol {\tau }(s) \times \hat {\boldsymbol {r}}(s,\theta ))$, solve (4.2) using the numerical method of § 4.1 to obtain the surface force density $\hat {\boldsymbol {f}}(s,\theta )$.

2. (ii) Compute the centreline force $\boldsymbol {f}(s)$ and parallel torque $n^\parallel (s)$ from $\hat {\boldsymbol {f}}(s,\theta )$ via (2.15) and (2.16).

3. (iii) Compute the SBT velocity induced from the BI force as

   (4.3)\begin{align} 8 {\rm \pi}\mu\boldsymbol{U}_{SB}&= \ln{\left(\epsilon^{{-}2}\right)} \left(\boldsymbol{I}+\boldsymbol{\tau} \boldsymbol{\tau}\right) \,\boldsymbol{f} + \left(\boldsymbol{I}-3\boldsymbol{\tau} \boldsymbol{\tau}\right) \,\boldsymbol{f}+\boldsymbol{U}_{f}^{{(FP)}}\nonumber\\ &\quad + \left(\ln{\left(\epsilon^{{-}2}\right)}-k\right) \left(\boldsymbol{\tau} \times {\partial s}{\boldsymbol{\tau}}\right) \frac{n^\parallel}{2}+{\boldsymbol{U}}^{{(FP)}}, \end{align}
   (4.4)\begin{equation} 8 {\rm \pi}\mu\varPsi^\parallel_{SB}= \frac{2n^\parallel}{a^2 \rho^2}+\frac{1}{2}\left(\ln{\left(\epsilon^{{-}2}\right)}-k\right)\left(\boldsymbol{\tau} \times {\partial s}{\boldsymbol{\tau}}\right) \boldsymbol{\cdot} \,\boldsymbol{f}+{\varPsi}^{(FP)}, \end{equation}
   where the finite part integrals are defined in (2.8), (2.12) and (3.26), and we have used the local drag coefficients for ellipsoidally decaying radius (4.1).

4. (iv) Calculate the $L^2$ norm of the error $\boldsymbol {U}_{SB}-\boldsymbol {U}$ and $\varPsi ^\parallel _{SB} - \varPsi ^\parallel$, and plot it as a function of $\epsilon$.

Note that we apply the forward SBT operator, which gives velocity from force/torque, in (4.3) and (4.4). This avoids inverting the slender body mobility, which is an ill-posed problem (Götz 2001; Tornberg & Shelley 2004). To evaluate the finite part integrals, we use the singular quadrature scheme developed in Tornberg (2020) and Maxian et al. (2021, 2022).

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

4.1. Numerical method

To solve (4.2), we need a way to integrate the singularity at $(s',\theta ') = (s,\theta )$. We follow the idea of Batchelor (1970b) and Koens (2022) in adding and subtracting the flow from a translating, rigid prolate ellipsoid, which has constant $\hat {\boldsymbol {f}}$ along its surface. The geometry of the prolate ellipsoid is chosen to exactly cancel the singularity in (4.2). Specifically, we introduce an ellipsoid with geometry

(4.5)\begin{equation} \hat{\boldsymbol{X}}_e\left(s',\theta',s,\theta\right) = \boldsymbol{X}(s) + \ell(s,\theta) \left(s'-s_c(s)\right) \boldsymbol{\tau}(s) + a_e(s) \hat{\boldsymbol{r}}(s,\theta') \sqrt{1-\left(s'\right)^2}. \end{equation}

This ellipsoid has half-axis length $\ell$ and radius $a_e$. The point $s_c(s)$ is the coordinate on the ellipsoid at which we want to cancel the singularity at $\hat {\boldsymbol {X}}(s,\theta )$. Following Koens (2022) again, we cancel the singularity by imposing the conditions

(4.6)\begin{equation} \left.\begin{gathered} \hat{\boldsymbol{X}}_e \left(s_c,\theta',s,\theta\right) =\hat{\boldsymbol{X}} \left(s,\theta\right) \\ \partial_{s'} \hat{\boldsymbol{X}}_e \left(s_c,\theta,s,\theta\right)= {\partial s}{\hat{\boldsymbol{X}}} (s,\theta) \quad \partial_{\theta'} \hat{\boldsymbol{X}}_e \left(s_c,\theta,s,\theta\right)=\partial_\theta \hat{\boldsymbol{X}} (s,\theta) \end{gathered}\right\}. \end{equation}

Substituting the definition of $\hat {\boldsymbol {X}}$ from (2.1) and $\hat {\boldsymbol {X}}_e$ from (4.5), the matching conditions (4.6) give three equations for the unknown parameters $\ell, a_e$ and $s_c$ in (4.5),

(4.7)$$\begin{gather} a \rho = a_e \sqrt{1-s_c^2}, \end{gather}$$

(4.8)$$\begin{gather}1+a \rho \left(\boldsymbol{\tau} \cdot {\partial s}{\hat{\boldsymbol{r}}} \right) = \ell , \end{gather}$$

(4.9)$$\begin{gather}a {\partial s}{\rho} ={-}\frac{a_e s_c}{\sqrt{1-s_c^2}}. \end{gather}$$

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

(4.10)\begin{equation} s_c = \frac{\rho - \sqrt{\rho^2 + 4 \left({\partial s}{\rho}\right)^2}}{2\rho}, \end{equation}

with $a_e$ determined from $s_c$ using (4.7). When the fibre is a straight ellipsoid, $\rho (s) = \sqrt {1-s^2}$ and ${\partial s} {\hat {\boldsymbol {r}}}=0$, so $a_e=a$ and $s_c=s$, and therefore the effective ellipsoid reduces to the ellipsoid itself.

Because the traction multiplied by the area element on a translating prolate ellipsoid is constant, the singularity subtraction scheme for (4.2) can now be written as (Koens 2022, (5))

(4.11)\begin{align} \hat{\boldsymbol{U}}(s,\theta) &= \boldsymbol{M}_e(s,\theta) \,\boldsymbol{f}(s,\theta)+ \frac{1}{8{\rm \pi}\mu} \int_{{-}1}^1 \int_{0}^{2{\rm \pi}} \left(\mathbb{S}\left(\hat{\boldsymbol{X}}(s,\theta),\hat{\boldsymbol{X}} \left(s',\theta'\right)\right) \,\boldsymbol{f}\left(s',\theta'\right)\right.\nonumber\\ &\quad\left. -\,\mathbb{S}\left(\hat{\boldsymbol{X}}_e(s_c,\theta,s,\theta),\hat{\boldsymbol{X}}_e\left(s',\theta',s,\theta \right)\right) \,\boldsymbol{f}\left(s,\theta\right)\right) {\rm d}\theta' \,{\rm d}s' , \end{align}

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

(4.12a,b)\begin{equation} e = \sqrt{1-\frac{a_e^2}{\ell^2}}, \quad q = \log{\left(\frac{1+e}{1-e}\right)}, \end{equation}

the resistance on a translating prolate ellipsoid is given in Keaveny & Shelley (2011, (44)–(45)) and in Padrino, Sprittles & Lockerby 2020 as

(4.13a,b)\begin{equation} \xi_\parallel= \frac{16 {\rm \pi}\ell e^3 \mu}{(1+e^2)q-2e} \qquad \xi_\perp= \frac{32 {\rm \pi}\ell e^3 \mu}{(3e^2-1)q+2e}. \end{equation}

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

(4.14)\begin{equation} \boldsymbol{M}_e = \frac{4 {\rm \pi}}{\xi_\parallel}\boldsymbol{\tau}(s)\boldsymbol{\tau}(s) + \frac{4{\rm \pi}}{\xi_\perp}\left(\boldsymbol{I}-\boldsymbol{\tau}(s) \boldsymbol{\tau}(s)\right). \end{equation}

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

4.2. Rotating curved fibre

To extract the unknown coefficient $k$ in the translation–rotation coupling term in (4.3), we consider the half-turn helix

(4.15)\begin{equation} \boldsymbol{X}(s) = \frac{1}{\sqrt{2}}\left(\frac{1}{\rm \pi}\cos{\left({\rm \pi} s\right)},\frac{1}{\rm \pi}\sin{\left({\rm \pi} s\right)}, s\right), \end{equation}

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

In Appendix C, we verify first-order accuracy in obtaining the cross-sectional force $\boldsymbol {f}(s)$ in (2.15) from the traction jump $\hat {\boldsymbol {f}}(s,\theta )$ in (4.11) as we refine $N_\theta$ and $N$ together. We also show how the errors $||{\boldsymbol {U}_{SB}-\boldsymbol {U}}||_{L^2}$, which are obtained from the traction jump via (2.15) and then (4.3), change as we refine $N_\theta$ and $N$ for three different values of $k$ ranging from $k= 2$ (which corresponds to throwing out the second line of the rotlet line integral (2.11)), to $k=3$ (which corresponds to simplifying (2.11) to separate the angle-dependent motions from pure translational ones, see (5.2)). For translational velocity, the finest feasible discretization is not sufficiently refined to see complete saturation of $||{\boldsymbol {U}_{SB}-\boldsymbol {U}}||$ (see Appendix C), so we report error bars in figure 1 that give the difference in $||{\boldsymbol {U}_{SB}-\boldsymbol {U}}||$ between our last two levels of refinement based on the results in Appendix C, the numerical error in the BI calculations likely causes an underestimate of the SBT error for $k=2$ and $k=2.5$ (since the error increases from the second-finest to finest discretization), while for $k=2.85$ and $k=3$ we likely overestimate the SBT error (since the error decreases from the second-finest to finest discretization).

Figure 1. Error in the rotation–translation SBT relative to the boundary integral calculation for various $\epsilon$ and $k$ and two different fibre shapes. The fibre shapes we consider, shifted so that their midpoints coincide, are shown in (a). For each fibre, we prescribe a rotational velocity $\varPsi ^\parallel \equiv 1/a^2$ with translational velocity $\boldsymbol {U}=\boldsymbol {0}$. We then solve the boundary integral equation (4.2) for the surface traction, use (2.15) and (2.16) to obtain the centreline force density and parallel torque density, and then (4.3) to obtain an SBT velocity $\boldsymbol {U}_{SB}$ from these densities. We plot the $L^2$ norm of $\boldsymbol {U}_{SB}$ (which should be zero) for several different values of $k$. For $k \neq 2.85$, results are shown only for the half-helix (4.15), and the error is clearly $O(1)$ with respect to $\epsilon$. For $k=2.85$, the error decreases with $\epsilon$ at a rate of roughly $\epsilon ^1$, independent of the fibre shape considered. The data points here come from the most refined discretization in figure 2, while the error bars are the difference between the most refined and second-most refined calculation.

Figure 2. Convergence plot for the curved rotating filament (4.15) studied in § 4.2. We show the self-convergence of our numerical method for (a) the force $\boldsymbol {f}$ and (b) the parallel torque $n^\parallel$ on each cross-section. The ‘error’ is the $L^2$ difference of each quantity relative to the next level of refinement, normalized by the $L^2$ norm of the most refined solution.

To find the approximately optimal value for $k$, we systematically vary $k$ between 2 and 3 in increments of 0.05, with the step size being relatively large because we are constrained by the accuracy of our boundary integral calculations. In figure 1, we show the SBT errors for $k=2, 2.5$ and 3, which represent evenly spaced values between 2 and 3, and compare them with our optimal value of $k=2.85$, which has the smallest errors for $\epsilon \leq 0.01$. For $k=2$ (which is our first guess from the first line of the rotlet integral expansion (2.11)), and $k=2.5$, we observe an error which saturates to a constant as we decrease $\epsilon$, meaning that the error in the SBT (4.3) when $k=2$ or $2.5$ is asymptotically $O(1)$. Increasing to $k=2.85$, we observe saturated errors that are roughly proportional to $\epsilon$, which is exactly what we want to achieve in SBT. When we increase again to $k=3$, we see errors which initially decrease with $\epsilon$, but appear to plateau to an $O(1)$ value. It appears then, that there is an optimal $k$ around $k=2.85$ which makes the SBT formula (4.3) accurate to $O(\epsilon )$.

To verify that $k=2.85$ is independent of the fibre geometry, in figure 1 we show as a light blue line the results when we perform the exact same test on a fibre with tangent vector

(4.16) \begin{equation} \boldsymbol{\tau}(s) = \frac{1}{\sqrt{2}}\left(\cos{\left(s^3(s-L)^3\right)}, \sin{\left(s^3 (s-L)^3\right)}, 1\right), \end{equation}

with the position obtained by integration. This fibre forms an ‘S’ shape, which is different from the half-turn helix of (4.15) (see figure 1a). While the shape is different, figure 1 shows that using $k=2.85$ is once again sufficient to match the SBT theory (4.3) with the three-dimensional boundary integral calculation, which suggests that the optimal $k$ is independent of the fibre shape.

While this is by no means a proof that the correct rotation–translation coupling coefficient in (4.3) is $k=2.85$ (indeed, there is no a priori reason to believe (4.3) or a similarly simple equation holds rigorously from (2.2) and (2.14)–(2.16)) it does provide us with a good educated guess that can be used to cheaply approximate boundary integral simulations with small error, and get an idea of how close our RPY formulae are to the true three-dimensional dynamics. In particular, using an RPY radius of ${\hat {a}} \approx 1.86 a$ in (3.42) gives the formula (3.43) with $k=2.85$. This value of the RPY radius is comparable, but fairly larger than, what we obtained by matching the translation–translation (${\hat {a}}_{tt} \approx 1.12 a$) and rotation–rotation (${\hat {a}}_{rr} \approx 1.06 a$) mobilities in § 3.5.

5. Conclusion

We examined the deficiencies of existing SBTs for rotating filaments, and then proposed a semiempirical SBT for tubular filaments derived from line integrals of RPY singularities. We showed that the singularity-based SBTs of Keller & Rubinow (1976) and Johnson (1980), which were designed with the resistance problem in mind, both suffer from the same issues for rotation–translation coupling in the mobility problem: when the fibre is driven by an applied torque at its base, their theories imply that torque induces an $O(1)$ angular-dependent translational velocity on the fibre cross-section, which violates the SBT boundary condition that the fibre only translates and rotates. In addition, these previous SBTs exhibit a lack of symmetry in the $2 \times 2$ mobility operator because of assumptions about the relative order of magnitude of the fibre translational and rotational velocities. As we showed previously (Maxian et al. 2022, § 5.3), neglecting rotation–translation coupling can lead to overestimation of the torque required to induce large-scale filament deformations. It is therefore important to derive an SBT which properly treats rotation–translation coupling.

As a workaround to the true problem of obtaining an SBT for a three-dimensional cylindrical (or prolate-spherioidal) filament, we used a different definition of a filament as a line of infinitely many RPY regularized singularities, similar to what has been done previously by Cortez and Nicholas for regularized Stokeslets (Cortez & Nicholas 2012). We then defined each component of the $2 \times 2$ mobility operator (translation–translation, rotation–translation, etc.) as a line integral of the corresponding fundamental regularized singularity given in (3.5). Because these regularized singularities, which define the mobility for a pair of ‘blobs’, give a symmetric $2 \times 2$ mobility, our mobility for fibres is automatically symmetric, which alleviates one of the two major problems with previous SBTs. The other issue, the $O(1)$ angular dependence of the translational velocity on a cross-section, is sidestepped because we use an isotropic regularized singularity function (a delta function on a sphere). This comes at a price: by using a series of spheres to represent the fibre centreline, we lose fidelity to the true three-dimensional geometry and the fluid flows near the filament surface. While this is far from ideal, it gives a mobility which is physical and simple to compute, since all that results is a local drag term plus a finite part integral for the translation–translation, rotation–translation and translation–rotation terms, for which special quadrature schemes have already been developed (Tornberg 2020).

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

Using the RPY kernels simplified a number of potential issues that arise for more rigorous models of the three-dimensional fibre geometry. First, we were able to easily derive mobility formulae near the fibre endpoints simply by modifying the domain of integration in the inner expansion (see Appendix B). This is in contrast to the more rigorous three-dimensional case, where a non-singular endpoint flow is obtained by assuming an ellipsoidally decaying radius function and distributing the centreline singularities only between the generalized foci of the body (Johnson 1980). Furthermore, the RPY mobility immediately generalizes to multiple filaments, since the flow induced on filament $j$ can also be computed by an integral of the RPY kernel over filament $i$. In SBT, while the pointwise velocity induced by singularities might be constant on a cross-section of fibre $i$, the same might not be true on a nearby fibre $j$, which leads to an ambiguity in how the induced flow on the centreline of fibre $j$ is assigned. This is especially relevant when fibres are in near contact, and shows why only three-dimensional boundary integral approaches are precise when fibres are $O(\epsilon )$ apart; such calculations are, however, very complex and expensive for suspensions of many fibres (Yan et al. 2022).

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

(5.1)\begin{equation} 8 {\rm \pi}\mu {\boldsymbol{U}} \approx \left(\ln{\left(\epsilon^{{-}2}\right)}-2.85\right) \left(\boldsymbol{\tau} \times {\partial s}{\boldsymbol{\tau}}\right)\frac{n^\parallel}{2}+{\boldsymbol{U}}^{{(FP)}}+O(\epsilon) \end{equation}

for the translational velocity due to parallel torque on an ellipsoidal filament, with the finite part integral ${\boldsymbol {U}}^{{(FP)}}$ defined in (2.12). Of course, this result is an estimate for a particular filament, and there is no way to show it is universal, although tests on multiple geometries gave the same $O(1)$ constant of $2.85$. Still, our guess for the rotation–translation coupling tensor is informed by the asymptotics of the RPY tensor, and so, if it is incorrect or non-universal, the RPY mobility (3.42) must miss some of the essential physics in the modelling of slender filaments. Given the number of studies that have used the RPY tensor to model filament hydrodynamics (e.g. Wang et al. 2012; Kallemov et al. 2016; Schoeller et al. 2021; Westwood & Keaveny 2021), this would be concerning.

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

(5.2)\begin{align} 8{\rm \pi} \mu \boldsymbol{u}^{(\varPsi)}\left(\hat{\boldsymbol{X}}\left(s,\theta\right)\right) &= \frac{2 n^\parallel}{a}\left(\boldsymbol{\tau} \times \hat{\boldsymbol{r}}\right) +\frac{n^\parallel}{2} \left(\ln{\left(\epsilon^{{-}2}\right)}-3\right) \left(\boldsymbol{\tau} \times {\partial s}{\boldsymbol{\tau}}\right)+{\boldsymbol{U}}^{{(FP)}}\nonumber\\ &\quad +\frac{n^\parallel \kappa}{2}\left(\boldsymbol{e}_b \cos{2\theta} -\boldsymbol{e}_n \sin{2\theta}\right)+ O\left(\epsilon\right). \end{align}

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

Given that it is difficult to find such a singularity representation, if it exists at all, our empirical approach in § 4 offers some promise, if only we could reach even smaller $\epsilon$. In the case when $\varPsi ^\parallel \sim 1/\epsilon ^2$, a more accurate numerical approach is necessary to do so, since in that case the surface tractions scale like $1/\epsilon$, and it becomes difficult to confirm $O(\epsilon )$ scaling of asymptotic expressions. And finally, there is an opportunity to analyse what Stokes-related PDE such an SBT is the solution of (this is an extension of the work of Mori, Ohm and others for translation (Mori et al. 2020a,b; Mori & Ohm 2021)).