title: Regularised non-uniform segments and efficient no-slip elastohydrodynamics
authors: Benjamin J. Walker Eamonn A. Gaffney
text: 
  	
  
  Forecasting with Multiple Seasonality
      Tianyang Xie
  School of Statistics 
  
  University of Minnesota Twin Cities
  
  Minneapolis, US 
  
  tianyang.xie.xty@gmail.com
  Jie Ding
  School of Statistics 
  
  University of Minnesota Twin Cities
  
  Minneapolis, US 
  
  dingj@umn.edu
      December 30, 2023
  ==============================================================================================================================================================================================================================
  
  
  
  
  
  The elastohydrodynamics of slender bodies in a viscous fluid have long been
  the source of theoretical investigation, being pertinent to the microscale
  world of ciliates and flagellates as well as to biological and engineered
  active matter more generally. Though recent works have overcome the severe
  numerical stiffness typically associated with slender elastohydrodynamics,
  employing both local and non-local couplings to the surrounding fluid, there
  is no framework of comparable efficiency that rigorously justifies its
  hydrodynamic accuracy. In this study, we combine developments in filament
  elastohydrodynamics with a recent slender-body theory, affording algebraic
  asymptotic accuracy to the commonly imposed no-slip condition on the surface
  of a slender filament of potentially non-uniform cross-sectional radius.
  Further, we do this whilst retaining the remarkable practical efficiency of
  contemporary elastohydrodynamic approaches, having drawn inspiration from the
  method of regularised Stokeslet segments to yield an efficient and flexible
  slender-body theory of regularised non-uniform segments.
  
  
  
  
  
  § INTRODUCTION
  
  
  
  
  The coupled elastohydrodynamics of flexible slender filaments are of intense
  interest to a breadth of active research communities, ranging from theoretical
  to experimental studies of filaments from the perspectives of synthetic
  sensors to those rooted in the biology and mechanics of cilia and flagella
  <cit.>.
  A comprehensive summary of the field is given in the recent review of
  <cit.>, which notes a particular need for further theoretical
  development in this area. Indeed, up until recently, problems involving
  filament elastohydrodynamics have been largely out of reach due to severe
  numerical stiffness associated with the dynamics of a slender body in a
  viscous fluid, with few studies being able to utilise large computing
  resources to combat this issue <cit.>.
  However, the work of <cit.> sought to address such problems,
  integrating the governing equations of elasticity in space in order to
  generate a coarse-grained framework with greatly reduced numerical stiffness.
  Despite being a recent development in the field, this approach has already
  been extended by <cit.> and <cit.> to include
  improved non-local hydrodynamics, applied to the model biological problem of
  flagellar efficiency <cit.>, and extended to motion in three
  dimensions <cit.>.
  
  
  Common to these recent models, as well as to other treatments of slender
  filaments at zero Reynolds number, are simplified representations of
  slender-body hydrodynamics. The aforementioned work of <cit.>
  utilises resistive force theory, a local relation between motion and drag that
  has seen widespread use since its advent in the 1950s
  <cit.>. More refined and complex are slender-body
  theories, which capture the non-local coupling of kinematics and associated
  forces via an integral relation, as considered in the early studies of
  <cit.> and later refined by
  <cit.>. Use of these slender theories in numerical applications
  often necessitates the use of many-point quadrature rules or specialised
  techniques to evaluate the integral of a rapidly varying or singular kernel,
  issues also found in methods derived from the boundary integral formulation of
  Stokes equations, as summarised by <cit.>. In the early 2000s,
  <cit.> circumvented such issues of numerical complexity by instead
  considering solutions of the regularly forced Stokes equations, leading to a
  regularised Green's function and an associated regularised theory. In turn,
  drawing from significant earlier study of singular slender-body theories, this
  led to commonplace use of a regularised slender-body theory ansatz for flow
  around a slender filament in terms of a force density f⃗, typically an
  integral over the centreline of the filament of the form
  
      u⃗(x⃗) = ∫K^ϵ(x⃗,s')f⃗(s')s' ,
  
  where u⃗(x⃗) is the fluid velocity at a point x⃗ and
  K^ϵ is a regular integral kernel. The parameter ϵ
  represents a lengthscale of the regularisation, which in studies of filament
  dynamics has often been taken to be the filament radius without rigorous
  justification
  <cit.>, with
  circular cross sections invariably assumed. The general ansatz of
  <ref> is also commonly used in conjunction with the
  hydrodynamic no-slip condition, though is evaluated not on the surface of the
  body, but on the filament centreline. With many approaches taking the integral
  kernel K^ϵ to simply be the regularised point force
  Green's function in the appropriate domain, application of this approximate
  relation does not guarantee that the no-slip boundary condition is satisfied
  on the surface of the body, with particular issues arising at the endpoints of
  the flagellum, where more than a velocity Green's function can be required
  <cit.>.
  
  Building upon the singular work of <cit.> and the classical
  solution of <cit.> for a prolate ellipsoid, the recent theory of
  <cit.> surpasses these general shortfalls and leverages a
  particular choice of kernel K^χ, along with a systematically
  justified and spatially dependent regularisation parameter χ, to satisfy
  the no-slip boundary condition on the surface of a slender body up to errors
  algebraic in the body aspect ratio. This theory retains the non-singular
  nature and accompanying numerical simplicity of the general regularised
  ansatz, whilst affording systematically justified accuracy and
  parameterisation. With such features having been absent from the recent
  efficient frameworks of <cit.>, the
  primary aim of this study is to incorporate the theory of <cit.>
  into the coarse-grained elastohydrodynamic framework of <cit.>,
  enabling the efficient simulation of slender bodies with asymptotically
  justified hydrodynamic accuracy in the no-slip condition. In doing so, we will
  additionally attempt to address concerning oscillations present in the force
  density solutions of these frameworks, which reportedly persist even with
  improved filament discretisations <cit.>.
  
  However, whilst the incorporation of the simple ansatz of <cit.>
  may be achieved with relative ease, integration of the regular but rapidly
  varying kernels may limit the speed of computation if performed with
  quadrature, as implemented in the original work of <cit.>. Having
  built upon the works of <cit.> and <cit.>,
  respectively, <cit.> and <cit.> avoid such
  expensive computation by analytically integrating the kernel over the straight
  line segments that form the discretised centreline of the slender body, which
  we will refer to as the regularised Stokeslet segment (RSS) approach. Though
  complicated here by a non-constant regularisation parameter χ, we will
  aim to proceed in a similar fashion and remove the reliance on quadrature
  rules in order to realise a highly efficient numerical framework for the study
  of slender-body elastohydrodynamics.
  
  
  Hence, we will proceed by first defining the non-uniform filament problem,
  adopting and unifying the notation of <cit.> and
  <cit.> for slender-body kinematics. We then describe a
  modification of the coarse-grained framework of <cit.>, similar in
  form to that of <cit.>, and present the slender-body theory of
  <cit.> cast in dimensionless quantities. Having adopted a
  piecewise-constant discretisation of viscous force density, we then seek to
  perform the slender-body integrals analytically, Taylor expanding the
  regularisation parameter χ to yield symbolic tractability. We will then
  numerically evidence the improved satisfaction of the no-slip boundary
  condition on the surface of the filament attained with the presented
  methodology, in turn considering the computed profiles of force density along
  the centreline of the filament and their behaviour near the endpoints of the
  slender body.
  
  
  § THE NON-UNIFORM FILAMENT PROBLEM
  
  
  
  In this work we will consider the planar motions of a thin inextensible,
  unshearable, untwistable filament in a viscous fluid, with the filament
  centreline denoted x⃗(s,t)= x(s,t)+y(s,t), without loss of
  generality, where , are constant orthogonal unit vectors in a fixed
  inertial reference frame and span the plane of motion. Here, s∈[0,L] is an
  arclength parameter and time is denoted by t, where L is the length of the
  slender object. Distinct from the notation of the Introduction, this
  slenderness is captured by the dimensionless parameter ϵ, defined
  explicitly as
  
      ϵ = 2max_s∈[0,L]{η(s)}/L≪1 ,
  
  where η(s) is the non-negative radius of the filament at arclength s,
  having assumed local axisymmetry about the centreline. With the shape
  therefore entirely defined by the centreline and radius function, we may
  describe points on the surface of the filament as
  
      x⃗^S(s,ϕ) = x⃗(s) + η(s)e⃗_r(s,ϕ) ,
  
  where ϕ is a cross-sectional angle. Here, e⃗_r is a radial unit
  vector embedded in a transverse cross section to the centreline. For unit
  tangent, normal, and binormal unit vectors defined by the Frenet-Serret
  relations
  
      e⃗_t(s) = ∂x⃗/∂ s,     ∂e⃗_t/∂ s = θ_s e⃗_n(s),     e⃗_b(s) = e⃗_t(s) e⃗_n(s) ,
  
  where θ(s,t) defines the filament tangent angle relative to , we
  define
  
      e⃗_r(s,ϕ) = e⃗_n(s) cosϕ + e⃗_b(s) sinϕ .
  
  Here and throughout, subscripts of s denote derivatives with respect to
  arclength and we have omitted writing the inherent time dependence of the
  filament centreline and all derived quantities. These definitions are
  illustrated in <ref>.
  
  
  We discretise the filament centreline into N linear segments, with the
  endpoints of these segments denoted by x⃗(s_i) for uniformly spaced
  arclengths s_i=(i-1)L/N ∈[0,L], where i=1,…,N+1. We write
  t⃗_i for the unit tangent to each linear segment, noting that this is
  an approximation of e⃗_t(s) on the ith segment, and parameterise
  these discrete tangents by θ(s), itself discretised as
  θ(s)≈θ_i on the ith segment such that t⃗_i =
  cosθ_i + sinθ_i. With this piecewise linear
  discretisation of x⃗ in arclength, or equivalently a piecewise constant
  discretisation of θ, we may describe the position of the filament with
  only the N+2 quantities x_1,y_1,θ_1,…,θ_N, where x⃗_1
  = x_1+y_1. Explicitly, for j=1,…,N+1 we have
  
      x⃗_j = x⃗_1 + ∑_i=1^j-1(cosθ_i + sinθ_i) ,
  
  where  is the constant segment length, equivalently defined as =
  L/N. Differentiating with respect to time, denoting time derivatives with a
  dot, this gives the linear velocity of the material point x⃗_j as
  
      ẋ⃗̇_j = ẋ⃗̇_1 + ∑_i=1^j-1(-sinθ_i + cosθ_i)θ̇_i .
  
  We may concisely write this latter linear relation as
  
      Qθ̇⃗̇ = Ẋ⃗̇ ,
  
  where θ⃗ = [x_1,y_1,θ_1,…,θ_N]^T, X⃗ =
  [x_1,y_1,…,x_N+1,y_N+1]^T and Q is the linear
  operator encoding <ref>, the latter having dimension
  (2N+2)× (N+2) and given explicitly in the work of <cit.>.
  Hence, we may readily cast expressions involving Ẋ⃗̇ in terms of
  the reduced variables θ⃗ and their time derivatives.
  
  The equations governing the surrounding fluid medium will be the familiar
  Newtonian Stokes equations, valid in the inertia-free limit of zero Reynolds
  number, which we will assume throughout. This limit is relevant to a broad
  range of biological and physical circumstances, for example the small-scale
  beating of spermatozoan flagella or the bending of cilia in flow. The Stokes
  equations may be briefly stated as
  
      μ∇^2u⃗ = ∇ p ,   ∇·u⃗=0 ,
  
  where u⃗ is the fluid velocity, μ is the associated viscosity and
  p is the pressure. Here we will also assume that the flow is in an unbounded
  domain in the exterior of the filament, and decays to zero in the far field.
  
  
  § NO-SLIP ELASTOHYDRODYNAMICS
  
  
  
  
  
   §.§ Coarse-grained mechanics
  
  Following <cit.>, we state the governing equations of
  elasticity for this slender inextensible unshearable filament in pointwise
  form as
  
      n⃗_s - f⃗   = 0⃗ ,
      m⃗_s + x⃗_s n⃗   = 0⃗ ,
  
  for contact force and couple denoted n⃗,m⃗ respectively and where
  a subscript of s denotes differentiation with respect to arclength. Here and
  throughout, the filament is passive, with no driving internal couple and
  f⃗ denotes the force per unit length applied on the surrounding fluid
  by the filament. Note that the external couple exerted by the fluid on the
  filament is O(ϵ^2), which will be negligible at the level of
  asymptotic approximation that we will consider in this work. To proceed, we
  integrate these equations with respect to arclength s, yielding
  
      - ∑_j=1^N∫_s_j^s_j+1f⃗(s)s   = n⃗(0) , 
      
      	- ∑_j=i^N∫_s_j^s_j+1 (x⃗(s)-x⃗_i)f⃗(s)s   = m⃗(s_i) ,    i = 1,…,N ,
  
  where we have decomposed the integrals into those over discrete segments and
  integrated the pointwise moment balance from s=s_i to s=s_N+1=L for
  i=1,…,N. In writing
  <ref> we have
  assumed that the filament is force and moment free at s=L, equivalent to
  imposing n⃗(L)=m⃗(L)=0⃗. We additionally assume that these
  conditions hold at the base, so that n⃗(0) = m⃗(0) =
  0⃗, though each of these boundary conditions may be readily replaced
  with those appropriate for particular problem settings, for example the
  clamping of one end of the filament. Recalling that the considered filament
  motion is purely planar, each term of <ref> is
  proportional to =, with m⃗(s_i)=m(s_i), so that
  <ref> collapses onto N scalar equations. We
  adopt a simple constitutive law, writing m⃗(s_i)=EIθ_s(s_i)≈
  EI(θ_i-θ_i-1)/ for bending stiffness EI, valid for
  i=2,…,N.
  
  Illustrated in <ref>, we discretise the force
  density f⃗, adopting a piecewise constant representation that is
  distinct from that of θ. Denoting the value taken by f⃗ at the
  segment endpoints x⃗_i by f⃗_i, for i=1,…,N+1, we
  discretise f⃗ as
  
      f⃗(s) = {[             f⃗_i ,     s∈[s_i,s_i+/2);           f⃗_i+1 , s∈[s_i+/2,s_i+1) , ].
  
  where i∈{2,…,N-1} is such that s∈[s_i,s_i+1). This is
  equivalent to stating that, on segments i=2,…,N-1, the value taken by
  f⃗ is equal to that at the closest segment endpoint, with the ith
  segment effectively split into two halves. The definition on the first and
  last segments is similar, though the segment is not precisely split into two
  equal parts, which will enable a concise description of the slender body
  theory in <ref>. Defining e=√(1-ϵ^2) to be the
  effective filament eccentricity, on the first segment we take
  
      f⃗(s) = {[        f⃗_1 , s∈[s_1,s_L^⋆);        f⃗_2 , s∈[s_L^⋆,s_2) ].   for  s_L^⋆ = 1/2(L(1-e)/2+) ,
  
  whilst on the last segment we analogously have
  
      f⃗(s) = {[          f⃗_N ,   s∈[s_N,s_R^⋆);        f⃗_N+1 , s∈[s_R^⋆,s_N+1) ].   for  s_R^⋆ = 1/2(L(1+e)/2+) .
  
  Whilst this is somewhat cumbersome, with the first and last segments being
  treated differently to the others, we have found that it yields significant
  advantages over simpler piecewise constant and linear schemes found in the
  literature. In particular, attempts at a piecewise linear approximation, as in
  <cit.>, result in large endpoint oscillations in the computed
  values of f⃗, akin to those found in the regularised Stokeslet segment
  methodology of <cit.> and are examined further in
  <ref>, where we evidence a lack of such
  oscillations in the approach presented in this study. A natural alternative,
  in which f⃗ is constant on each segment, yields equivalently
  undesirable results, with the methodology becoming numerically intractable due
  to stiffness when considering nearly straight filaments. Indeed, the same
  issue is present in the scheme proposed by <cit.>, which
  utilises this intuitive discretisation. Though these issues are circumvented
  by the approach presented in this work, the source of this sensitive numerical
  dependence of the filament problem on discretisation remains unclear, and
  warrants future investigation.
  
  
  Returning to the now-discretised filament problem, the force-density
  dependence of <ref>
  may be cast as a simple linear operator, denoted by B, allowing us to write
  
      -BF⃗ = R⃗ ,
  
  where R⃗=[0,0,m(s_1),…,m(s_N)]^T encodes the bending moments and
  total force acting on the filament, whilst
  F⃗=[f⃗_1·,f⃗_1·,…,f⃗_N+1·,
  f⃗_N+1·]^T is the vector of discretised force densities. The
  first two rows B_1,B_2 of B represent total
  force balance over the filament, and are given explicitly by
  
      B_1    = /2[1+d,0,2-d,0,2,0,…,2,0,2-d,0,1+d,0] ,
      B_2    = /2[0,1+d,0,2-d,0,2,0,…,2,0,2-d,0,1+d] ,
  
  where d=L(1-e)/(2). The remaining rows B_i+2 encode the
  integrated moment balance equations for i=1,…,N, the expressions for
  which are given in <ref>.
  
  We now suppose that an invertible linear operator A may be
  constructed such that
  
      Ẋ⃗̇ = AF⃗(1 + O(ϵ)) ,
  
  which we will find explicitly in <ref>. Upon substitution of
  <ref> into <ref>, also making use of <ref>, we
  obtain the leading-order coarse-grained linear system
  
      -BA^-1Qθ̇⃗̇ = R⃗ ,
  
  where B encodes the integrated equations of elasticity,
  A represents the hydrodynamic relation between velocity and force
  density, Q links the kinematic descriptions of the filament, and
  R⃗ is the elastic response of the filament to bending.
  
  Finally, we non-dimensionalise lengths by filament half-length L/2, forces with 4EI/L^2, and time with some characteristic time scale T. This yields the dimensionless system
  
      -E_h B̂Â^-1Q̂θ̇̂̇⃗̇̂̇ = R̂⃗̂ ,    E_h = πμ L^4/2EI T ,
  
  where the notation ·̂ denotes dimensionless quantities, and we note
  that the rescaled arclength parameter is ŝ = 2s/L ∈[0,2]. The
  elastohydrodynamic number E_h here is analogous to that of
  <cit.>, though differs by a factor of 16 due to differing choices
  of lengthscale. We have the explicit relations
  
      B = L^2/4B̂ ,   A = 1/8πμÂ ,   Qθ̇⃗̇ = L/2TQ̂θ̇̂̇⃗̇̂̇ ,   R⃗ = 2EI/LR̂⃗̂ ,
  
  between dimensional and dimensionless quantities, having multiplied the force
  balance equations by /2 and absorbed the dimensional scalings of
  x_1,y_1 in to Q for covenience, writing θ̂⃗̂ =
  (x̂_1,ŷ_1,θ_1,…,θ_N)^T. In what follows, we will
  drop the ·̂ notation for dimensionless variables, though for later
  convenience we first write
  
      η(s) = L/2η̂(ŝ) =
      	ϵ L/2η̃(ŝ)
  
  and immediately drop the tilde on η̃(ŝ)∼ O(1). For clarity, the
  points on the surface of the filament may now be written in terms of
  dimensionless quantities as
  
      x⃗^S(s,ϕ) = x⃗(s) + ϵη(s)e⃗_r(s,ϕ) .
  
  
  
  
   §.§ Non-uniform hydrodynamics
  
  Before describing the slender-body theory of <cit.> that we will
  use to relate forces and flow, we first recapitulate the well-known
  regularised singularities of <cit.> and <cit.> on
  which it is built. Following <cit.>, for points
  α⃗,β⃗ and particular choices of mollifier the regularised
  Stokeslet S^χ and potential dipole D^χ are
  given by
  
      S^χ(α⃗,β⃗) = (α⃗-β⃗^2 
      	+ 2χ)I/(α⃗-β⃗^2 + χ)^3/2 +
      	  Q(α⃗,β⃗)/(α⃗-β⃗^2+χ)^3/2 ,
      	  
      D^χ(α⃗,β⃗) = -(α⃗-β⃗^2 - 2χ)I/(α⃗-β⃗^2
      	+ χ)^5/2 + 3Q(α⃗,β⃗)/(α⃗-β⃗^2+χ)^5/2 ,
  
  where Q(α⃗,β⃗) =
  (α⃗-β⃗)⊗(α⃗-β⃗), I is the 3×3
  identity tensor, and χ is the regularisation parameter.
  
  Throughout this section, it will be convenient to consider functions of
  filament arclength instead as functions of a shifted arclength parameter
  s'∈[-1,1], which will greatly simplify the notation associated with the
  slender-body theory of <cit.>. We will consistently abuse
  notation and write x⃗(s) ≡x⃗(s'), where s=s'+1 and other
  functions of arclength are treated analogously. In particular, this enables us
  to concisely define the arclength-dependent regularisation parameter
  χ=χ(s'), which may be written as
  
      χ(s') = ϵ^2[(1-s'^2) - η^2(s')] .
  
  Recalling the effective filament eccentricity as e=√(1-ϵ^2), the
  dimensionless ansatz of <cit.> for the fluid velocity at a point
  y⃗ in terms of the force per unit length f⃗(s') may now be
  written as
  
      u⃗(y⃗) =
      	∫_-e^e[S^χ(s')(y⃗,x⃗(s'))
      	-
      	1-e^2/2e^2(e^2-s'^2)D^χ(s')(y⃗,x⃗(s'))]f⃗(s')s' ,
  
  noting that the dimensional factor of 8πμ has been absorbed by the
  scalings of <ref>. Note that the limits in the integral are
  between -e and e, rather than -1 and 1, as inherited ultimately from
  the  solution for a translating prolate ellipsoid, as
  detailed in <cit.>. Taking y⃗=x⃗^S(s_i',ϕ), where
  s_i' = s_i - 1 are the shifted dimensionless arclengths corresponding to the
  discrete points s_i, we may apply this ansatz at the filament surface to
  generate the N+1 vector equations
  
      u⃗(x⃗^S(s_i',ϕ)) =
      	∫_-e^e[S^χ(s')(x⃗^S(s_i',ϕ),x⃗(s')).
      .
      	-
      	1-e^2/2e^2(e^2-s'^2)D^χ(s')(x⃗^S(s_i',ϕ),x⃗(s'))]f⃗(s')s' .
  
  We impose the no-slip condition u⃗(x⃗^S(s_i',ϕ)) =
  ẋ⃗̇^S(s',ϕ) on the surface of the filament, and may decompose
  
      ẋ⃗̇^S(s',ϕ) = ẋ⃗̇(s') + ϵω⃗(s')η(s')e⃗_r(s',ϕ) ,
  
  for centreline velocity ẋ⃗̇(s') and angular velocity
  ω⃗(s') measured about x⃗(s'), recalling that the filament
  is assumed to be unshearable. Supposing that ω⃗ is O(1) as
  ϵ→0, consistent with the filament being assumed untwistable
  and planar, at leading order we simply have
  
      ẋ⃗̇^S(s',ϕ) = ẋ⃗̇(s') + O(ϵ) ,
  
  independent of ϕ. Finally, we arrive at the leading-order relation
  
      ẋ⃗̇(s_i') ≈∫_-e^e[S^χ(s')(x⃗^S(s_i',ϕ),x⃗(s')).
      .
      	-
      	1-e^2/2e^2(e^2-s'^2)D^χ(s')(x⃗^S(s_i',ϕ),x⃗(s'))]f⃗(s')s' .
  
  For comparison, the equivalent expression used in the method of regularised
  Stokeslet segments, and indeed many regularised slender body theories
  <cit.>, may be written as
  
      ẋ⃗̇(s_i') ≈∫_-1^1S^ϵ/2(x⃗(s_i'),x⃗(s'))f⃗(s')
      		s' ,
  
  where we note in particular that the evaluations of the regularised Stokeslet
  kernel are on the filament centreline, not on the surface of the slender body.
  
  Though the left-hand side of <ref> is trivially
  independent of the cross-sectional angle ϕ, it is not clear if the
  integral is similarly independent. However, with the particular choice of
  regularisation parameter χ(s') given in <ref>,
  <cit.> showed that the integral of <ref> is in
  fact independent of ϕ at leading order in ϵ, with errors linear
  in the aspect ratio. Thus, <ref> satisfies this
  necessary condition. Moreover, its solution enables the no-slip condition on
  the surface of the filament to be satisfied to O(ϵ). With the force
  density f⃗ discretised as described in <ref> and
  incurring errors proportional to df⃗/ds,
  which are assumed small and may be verified a posteriori, yields the
  leading-order linear system
  
      Ẋ⃗̇ = AF⃗
  
  relating force densities on the filament to the centreline velocities.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Regularised non-uniform segments
  
  The entries of A may be readily computed with quadrature, as was
  performed in the original work of <cit.>. However, with the
  integral kernels rapidly varying in some regions, this can be prohibitively
  expensive in elastohydrodynamic simulations, where numerous evaluations of
  A are typically required. Inspired by the recent method of
  regularised Stokeslet segments, as detailed by <cit.> and used in
  the work of <cit.>, we seek to evaluate these integrals
  analytically on each linear segment of the filament, in general incurring
  discretisation errors as a result of the arclength-dependent regularisation
  χ(s'). We will refer to this approach as the method of regularised
  non-uniform segments (RNS).
  
  With f⃗ approximated as piecewise constant, computation of A
  reduces to performing a number of integrals over the discrete straight
  segments of the filament. We consider such an integral over the
  jth segment, parameterised by α∈[0,1], on which
  we write χ=χ(α) and the filament centreline is given by
  x⃗(α)=x⃗_j -
  αv⃗, where v⃗=x⃗_j-x⃗_j+1. This formulation is
  applicable even to the first and last segments, subject to the substitution of
  x⃗_1 by x⃗(-e) and of x⃗_N+1 by x⃗(e), owing to
  the chosen discretisation of f⃗. With this discretised f⃗ taking
  the values f⃗_j and f⃗_j+1 on different halves of this
  segment, we require expressions for the integrals evaluated on two subdomains,
  α∈[0,1/2] and α∈[1/2,1]. Further, analogous to
  <cit.> and <cit.> and given explicitly in
  <ref>, we note that this discretisation of f⃗
  has rendered each of these integrals as linear combinations of
  
      T_m,q^L = ∫_0^1/2α^mR^qα ,    T_m,q^R = ∫_1/2^1α^mR^qα ,
  
  where we define R=(x⃗^S_i-x⃗_j+αv⃗^2 +
  χ(α))^1/2 for a surface point
  x⃗^S_i=x⃗^S(s_i',ϕ), with ϕ arbitrary as above. These
  integrals may be readily performed in the case that R^2 is a quadratic
  function of α, as in the original method of regularised Stokeslet
  segments though here prohibited in general by χ(α).
  
  In order to recover this desirable property, we Taylor expand χ(α)
  about an endpoint of the segment, either α=0 or α=1, assuming
  sufficient smoothness of χ. The expansion point is chosen in order to
  minimise the error in the resulting integral, noting in particular that
  R(α) can become O(ϵ) if x⃗_j-αv⃗ nears
  x⃗_i^S, for example in the trivial case of i=j. With χ therefore
  plausibly the dominant term in this O(ϵ) neighbourhood, we choose to
  expand χ(α) about the segment endpoint that is closest to
  x⃗_i^S, denoting the value of α at this endpoint as
  α^⋆. Collecting powers of α, in each case this yields an
  expansion of the form
  
      R(α)^2 = A + B α + Cα^2 + E ,   E≤{[     1/6α^3v⃗^3sup^3χs'^3 ,                    α^⋆=0 ,; 1/6(1-α)^3v⃗^3sup^3χs'^3 ,                    α^⋆=1 . ].
  
  In the error term E we have bounded the third derivative of χ over the
  segment, and have cast the derivative in terms of the normalised arclength
  s' in order to unify our phrasing of model assumptions. With
  R(α)=O(ϵ) when x⃗_i^S
  -x⃗_j+αv⃗=O(ϵ), and R(α) strictly order unity
  otherwise, when α^⋆=0 this error term is subdominant if
  
      1/6α^3^3sup^3χs'^3 = {[                     O(ϵ^3) , where x⃗_i^S-x⃗_j+αv⃗=O(ϵ) ,;                       O(ϵ) ,                   otherwise, ].
  
  noting that v⃗≤, with a similar expression required for
  α^⋆=1. This imposes a weak restriction on the derivatives of
  χ and the discretisation length , recalling that
  χ=O(ϵ^2) everywhere. Assuming that such a restriction holds, we
  drop the error term E in what follows, approximating R(α) as a
  quadratic function on each segment. The segment-dependent coefficients A,B,C
  may be readily computed when expanding with α^⋆=0 or
  α^⋆=1, and for α^⋆=0 are given explicitly by
  
      A = x⃗_i^S - x⃗_j^2 + χ ,   
      	B = 2v⃗·(x⃗_i^S - x⃗_j) + v⃗χs' ,   
      	C = v⃗^2(1 + ^2χs'^2) ,
  
  where evaluations of χ and its derivatives for A,B,C are at α=0
  and we henceforth write R(α)^2 = A + Bα + Cα^2 for brevity.
  Omitted here for brevity, analogous expressions hold for A,B,C when
  α^⋆=1. As noted above and written explicitly in
  <ref>, the integral kernel may be decomposed into a
  linear combination of terms α^mR^q for
  (m,q)∈{(0,-1),(0,-3),(0,-5),(1,-3),(1,-5),(2,-3),(2,-5),(3,-5),(4,-5)}.
  For m>0, computation of these quantities may be performed simply via the
  recurrence relations
  
      T_m+1,q-2^L    = .α^mR^q/qC|_0^1/2 -m/qCT_m-1,q^L - B/2CT_m,q-2^L ,
      
      	T_m+1,q-2^R    = .α^mR^q/qC|_1/2^1 -m/qCT_m-1,q^R - B/2CT_m,q-2^R ,
  
  where q,C≠0. These are analogous to the recurrence of <cit.>
  and are similarly derived via integration by parts. Thus, explicit calculation
  of T_m,q^L and T_m,q^R is required only for m=0, with the relevant antiderivatives given in <ref>.
  
  Hence, the construction of the operator A proceeds simply and
  efficiently: the coefficients A,B,C are evaluated from precomputed values of
  χ and its derivatives, the integrals T_m,q^L,T_m,q^R are computed
  for m=0 using the given antiderivatives, further integrals for m>0 are
  computed via the recurrences of <ref>, and the
  entries of A are formed as linear combinations of these terms
  following <ref>. We additionally note that this
  process may be readily generalised to evaluation points that do not lie on the
  surface of the filament, in this case Taylor expanding about the segment
  endpoint that is closest to the evaluation point.
  
  
  § VERIFICATION AND EXAMPLES
  
  
  
  
  
   §.§ Efficiency and accuracy against quadrature
  
  Construction of the operator A via the method of regularised
  non-uniform segments introduces local approximations of the regularisation
  parameter χ wherever it is not simply a quadratic function of arclength,
  enabling analytic integration. We now compare this approach with quadrature in
  terms of both accuracy and efficiency in a practical parameter regime,
  considering three dimensionless radius functions η(s') of varying
  complexity:
  
      [                        (a)                √(1-s'^2) ,;                        (b) √(1-s'^2)(1-0.1cos2π s') ,;                        (c)  √(1-s'^2)(1.1+sin9π s') , ]
  
  each subject to normalisation and shown in <ref>. Considering a
  filament with a curved centreline, corresponding to the initial condition of
  <ref>a, with N=100 and ϵ=0.02 we compute A
  using both the RNS methodology and the inbuilt  routine in
  MATLAB®, with the numerical quadrature set to
  a tolerance of e-12 and denoting the results of these computations by
  A_RNS and A_Q, respectively. We write
  ℰ for the relative matrix infinity norm error between these two
  results, defined explicitly as
  
      ℰ = A_RNS - A_Q_∞/A_Q_∞ .
  
  These relative errors are shown in <ref>, each of which can be seen
  to be several orders of magnitude lower than the asymptotic slenderness
  parameter. The rapidly varying curvature of case (c) gives rise to the largest
  error, consistent with the restrictions imposed on the derivatives of χ
  in <ref>. Computations were performed on modest hardware
  (Intel® Core™ i7-6920HQ CPU),
  with the walltime for the RNS method being over two orders of magnitude less
  than that of the quadrature implementation, representing a significant
  improvement in computational efficiency for minimal reduction in accuracy.
  These observations of efficiency and accuracy hold for a range of considered
  body centrelines and radius functions, and are robust to variations in the
  slenderness parameter ϵ.
  
  
  
  
  
   §.§ Invariants of free-filament motion
  
  The coarse-grained framework for filament elasticity is similar to that
  presented and derived in the recent work of <cit.>, where it was
  extensively verified and benchmarked, utilising the stiff solver
   provided in MATLAB® with
  relative and absolute tolerances of e-6 <cit.>. However,
  due to the modification of considering a piecewise constant discretisation of
  the force density f⃗, akin to the study of <cit.>, we
  additionally verify the presented methodology in the case of a relaxing
  symmetric filament. Having taken N=40 and ϵ=0.01, in
  <ref> we showcase the simulated dynamics of an initially symmetric
  filament relaxing to a straight configuration, during which we see that
  symmetry is preserved. Owing to the filament having no net force or torque act
  upon it, the centre of mass should not deviate from its initial position.
  Computing the translation of the centre of mass over the motion, a
  quantitative measure of framework accuracy, in <ref>b we see that
  this approximate constancy is preserved numerically with errors on the order
  of e-3L, improved by an order of magnitude when compared to the
  previous methodology of <cit.>.
  
  
  
  
  
   §.§ Comparison against existing theories
  
  We now more thoroughly compare and contrast the presented elastohydrodynamic
  framework against two existing approaches, in particular the published RSS
  methodology of <cit.> and a resistive force theory (RFT)
  formulation based on that of <cit.>. The latter RFT method is as
  described in the work of <cit.>, though we make use of the
  resistive coefficients of <cit.>, with the normal
  resistive coefficient twice that of the tangential coefficient.
  
  
  
    §.§.§ A relaxing filament
  
  We simulate the free relaxation of a bent filament, with the θ_i
  initially equally spaced and increasing between -π/4 and π/4 to
  correspond to a filament of constant curvature, via each of the three
  methodologies, picking a common but arbitrary elastohydrodynamic number of
  E_h=9600 and setting N=40. The filament has aspect ratio 1:100,
  corresponding to ϵ=0.02 in the RNS framework and ϵ=0.01 in
  the RFT and RSS approaches. Simulating until a dimensional time of
  100, at which point the RNS solution is nearing complete
  relaxation to a straight configuration, we display snapshots of the computed
  solutions and some associated metrics in <ref>. Immediately
  evident is a qualitative similarity between the computations, though there is
  some pairwise disagreement throughout the motion. Most prominent are
  differences in the timescale of relaxation, as can be seen in the maximum
  curvature plot of <ref>b, with the RFT solution relaxing more
  slowly than the predictions by non-local theories. We more concretely quantify
  the overall differences between methodologies at a given time t via the
  measure D, defined for a computed solution x⃗(s,t) by
  
      D^2(t) = 1/L∫_0^L
      	x⃗(s,t) -
      	x⃗_RNS(s,t)_2^2s ,
  
  relative to the RNS solution x⃗_RNS. The evolution of this
  distance measure for the RFT and RSS approaches is shown in
  <ref>c, and demonstrates that, whilst differences between
  solutions are indeed small, being on the scale of ϵ in this
  particular case, these distinctions persist throughout the motion.
  
  With elastohydrodynamic simulations appearing broadly similar at the level of
  detail considered thus far, we also note a common computational efficiency of
  the frameworks, with even the more complex regularised non-uniform segments
  approach computing the relaxation dynamics in a number of seconds. Indeed,
  this is replicated throughout further testing for each of a wide array of
  initial conditions, and is robust to variations in the radius function
  η(s') and the filament aspect ratio. Thus, despite employing a more
  sophisticated slender-body ansatz, we see retained in the RNS methodology the
  desirable efficiency associated with the existing coarse-grained frameworks.
  
  
  
  
  
    §.§.§ A simple filament in flow
  
  From the agreement seen above in the case of a relaxing filament, one might
  expect that the theoretical refinement offered by the RNS approach over the
  simpler and cruder RSS methodology is minimal in practice. However, more
  significant differences are indeed present, as we now highlight via a simple
  example.
  
  We consider perhaps the most simple possible filament simulation: the dynamics
  of an initially straight filament in a uniform background flow, with a
  background flow u⃗_b incorporated into the current framework via the
  mapping u⃗↦u⃗-u⃗_b as in the work of
  <cit.>. The simulated filament should exhibit trivial motion and
  deformation, merely translating with the background flow and retaining its
  straight configuration. Both the RNS and RSS methodologies successfully
  replicate this behaviour, and solution time is negligible. However, a noted
  issue of methods based on regularised Stokeslet segments and similar
  approaches are endpoint oscillations in the computed force density f⃗,
  present in each of the works of
  <cit.>, which persist even with mesh
  refinement.
  
  Here, we explicitly compute the force density on a straight filament of aspect
  ratio 1:100 in a unit background flow u⃗_b=e⃗_y using both the RNS
  and RSS approaches, where e⃗_y is perpendicular to the filament
  tangent. In <ref>a-c we present the magnitude of the
  computed force density on the filament from s=0 to s=L for various body
  radius functions, appealing to symmetry and noting that the force density is
  identically zero in the direction of the filament tangent. In each case, we
  observe the oscillations of the RSS force density near the endpoints of the
  slender body, with the RSS solution being fundamentally independent of the
  radius function, whilst the piecewise-constant RNS solution essentially
  eliminates these oscillations. We have taken N=200 in
  <ref>c in order to capture the highly oscillatory radius
  function of <ref>c, consistent with the error analysis of
  <ref>, taking N=100 in the other cases.
  
  
  
  Perhaps more pertinent, and indeed the motivation behind the use of the ansatz
  of <cit.>, is the velocity boundary condition on the filament. We
  explicitly evaluate the flow velocity on the surface of the filament via both
  the RNS and RSS methods, sampling at 1000 uniformly spaced points on the
  surface, and show the infinity norm error in the computed velocity in
  <ref>d-f as a function of dimensionless shifted arclength
  s'. Notably, the RSS approach is consistently inaccurate along the length of
  the slender body, yielding approximately 5% errors over the entire surface,
  corresponding to five times the regularisation parameter of the RSS method.
  The RNS methodology significantly improves upon this, with limitingly small
  error along the majority of each of the slender bodies in both
  <ref>d and <ref>e, with errors of
  approximately 2ϵ near the endpoints of the slender body in
  <ref>e. In particular, these errors are on the same order
  as those found in the original evaluation of the slender body theory by
  <cit.>, with the impact of moving away from quadrature therefore
  minimal in all but <ref>f, which is improved by reducing
   to once again accommodate the oscillatory radius function. Thus, we
  observe that the use of the RNS methodology affords significant gains in the
  accuracy of the no-slip boundary condition over other approaches. Convergence
  of this velocity error as a function of N and ϵ is illustrated in
  <ref> for the case of <ref>b.
  
  
  § DISCUSSION
  
  
  
  Though the study of <cit.> vastly increased the computational
  efficiency of filament simulations, it did so whilst employing only resistive
  force theory, with this leading order hydrodynamic relation typically
  conferring errors logarithmic in the filament aspect ratio. Subsequent works
  have extended this framework to feature improved hydrodynamics
  <cit.>, each making use of a simple but non-local
  regularised ansatz. However, even these works neglect the boundary condition
  on the body surface, instead evaluating velocities along the centreline when
  linking fluid velocity to applied force density. Via the evaluations performed
  in <ref> of this work, we have evidenced
  the relative inaccuracy of such approaches, observing non-negligible errors in
  the computed surface velocity over the entire length of the slender body, with
  these hydrodynamic errors being a fundamental weakness of previous
  methodologies. Incorporating a refined hydrodynamic ansatz, the presented
  regularised non-uniform segment methodology significantly improves upon such
  errors, with discrepancies in the velocity boundary condition present only at
  the filament endpoints, given adequate discretisation to account for the level
  of variation in the cross-sectional radius function. In particular, the
  slender-body theory employed here inherently takes into account the complex
  shape of the filament, enabling the study of realistic slender-body geometries
  and replacing previous imprecise justifications with analytically derived
  quantifications of accuracy.
  
  However, a naive incorporation of the slender-body theory of
  <cit.> into a coarse-grained framework of filament elasticity
  yielded large computation times, sacrificing the efficiency typically
  associated with the underlying approach of <cit.>. Indeed, whilst
  the use of automated quadrature rules allows computation of the hydrodynamic
  operator to any desired degree of numerical accuracy, even the regular
  integral kernel of the ansatz of <cit.> was insufficient to
  enable rapid computation on par with the existing frameworks of
  <cit.>. Thus, exploiting a low-degree
  approximation of the unknown force density f⃗, we instead computed the
  necessary integrals analytically, mimicking the approach of <cit.>
  after Taylor expanding the generally non-quadratic regularisation parameter
  χ. Quantifying the errors associated with this approximation, we have
  evidenced a remarkable accuracy and efficiency of this approach, yielding a
  scheme for elastohydrodynamic simulation that is comparable in computational
  cost to existing methodologies, whilst simultaneously improving on their
  accuracy. Thus, the presented framework will enable rapid solution of the
  forward elastohydrodynamic problem, pertinent to modern Bayesian parameter
  inference techniques, for example, along with explorations of fluid-structure
  interactions in slender-body systems. Further, the method of regularised
  non-uniform segments will more generally enable rapid application of the
  slender-body theory of <cit.>, facilitating future investigative
  and explorative studies into filament dynamics.
  
  Whilst efficiency gains were made by adopting the general principle of the
  method of regularised Stokeslet segments, the regularised non-uniform segment
  approach avoids a pertinent issue associated with the principles of the former
  theory. Present in the works and published codes of
  <cit.> are severe variations in the
  computed force density f⃗ near the endpoints of the considered
  filaments, persisting or indeed worsening with increased refinement of
  approximation. With force density a fundamental component of such
  elastohydrodynamic frameworks, these apparent errors may contribute
  non-negligibly to simulated dynamics and applications, particularly given the
  reported significance of distal activity in recent model spermatozoa
  <cit.>. Thus, the absence of comparable oscillations in the RNS
  solutions represents a significant advantage over these existing
  methodologies. Curiously, the insertion of the slender body theory of
  <cit.> alone into the framework of <cit.> was not
  sufficient to achieve this, as discovered during the author's initial attempt
  at formulating the RNS methodology, which differs to the presented approach
  only by using a piecewise linear discretisation of force density f⃗.
  However, the combination of this improved ansatz and a lower order
  discretisation of f⃗ successfully removed the unphysical oscillations
  from the computed solutions, yielding the smooth profiles seen in
  <ref>, though detailed investigation of the Fredholm
  integral equation of <ref> is required in order to
  ascertain the source of such pervasive errors. Future work may also include
  trivial extensions to the study of active filaments and general background
  flows, affording justified accuracy to the wide range of elastohydrodynamic
  problems made tractable by the work of <cit.>.
  
  In summary, we have integrated the fundamental advance of <cit.>
  and the regularised slender-body theory of <cit.>, overcoming
  their respective shortfalls to yield a framework for the efficient and
  accurate simulation of slender-body elastohydrodynamics. The so-called
  regularised non-uniform segment approach retains the flexibility of its parent
  models, and hence may be applied to a wide variety of biological and
  biophysical problems to afford increased accuracy over earlier approaches.
  Further, complex axisymmetric geometries may now be reliably modelled using
  this framework, previously only realisable with reduced fidelity or
  drastically increased computational effort. Applicable even more generally,
  this study has markedly improved the efficiency of the slender-body theory of
  <cit.>, with this work overall facilitating both the accurate
  quantification and large scale no-slip simulation of slender elasticity and
  hydrodynamics.
  
  
  B.J.W. is supported by the UK Engineering and Physical Sciences Research
  Council (EPSRC), grant EP/N509711/1.
  
  Declaration of Interests. The authors report no conflict of interest.
  
  
  
  
  
  § MOMENT BALANCE AS A LINEAR SYSTEM
  
  
  For i=1,…,N, the rows B_i+2 of B encode the
  integrated moment balance in terms of the f⃗_j, resultant of
  integrating over segments i through N. For each i, the summation of
  <ref> may be written simply as
  
      ∑_j=i^N I⃗_j
  
  for integrals I⃗_j. For j=2,…,N-1, these are given by
  
      I⃗_j = /2[x⃗_j - x⃗_i - 1/4v⃗_j] f⃗_j + /2[x⃗_j - x⃗_i - 3/4v⃗_j] f⃗_j+1 ,
  
  where v⃗_j = x⃗_j - x⃗_j+1. For j=1, again writing
  d=L(1-e)/(2), we have the modified expression
  
      I⃗_1 = /2[(1+d)(x⃗_j-x⃗_i) - 1/4(1+d)^2v⃗_1]f⃗_1 
      + /2[(1-d)(x⃗_j-x⃗_i) + (1/4(1+d)^2 - 1)v⃗_1]f⃗_2 ,
  
  whilst for j=N we have
  
      I⃗_N = /2[(1-d)(x⃗_j-x⃗_i) - 1/4(1-d)^2v⃗_N]f⃗_N 
      + /2[(1+d)(x⃗_j-x⃗_i) + (1/4(1-d)^2 - 1)v⃗_N]f⃗_N+1 .
  
  These expressions are self-consistent, as taking d=0 in the latter two
  yields the expression for I⃗_j.
  
  
  
  § INTEGRALS AS A LINEAR COMBINATION
  
  
  We decompose the integral of <ref> over a straight
  segment with endpoints x⃗_j and x⃗_j+1, adopting a piecewise
  constant discretisation of the force density f⃗, such that it takes the
  value f⃗_j on the half of the segment nearest to x⃗_j, and
  f⃗_j+1 otherwise. The limits of integration are determined by
  requiring either the coefficient of f⃗_j or that of f⃗_j+1,
  and for brevity we omit such limits here and will refer instead to the
  placeholder T_m,q in lieu of T_m,q^L and T_m,q^R in what follows,
  which should be appropriately substituted. Parameterising the straight segment
  by α∈[0,1], with x⃗(α)=x⃗_j-αv⃗, where
  v⃗ = x⃗_j -
  x⃗_j+1, and taking K^ϵ to be the kernel of
  <ref>, we may write the integral over the part of
  the segment as
  
      ∫K^ϵ(x⃗,s')f⃗(s') s' =
      K^ϵ_I f⃗^⋆ ,
  
  where f⃗^⋆ is the constant force density over the domain of
  integration, which is either α∈[0,1/2] or α∈[1/2,1]. The
  operator K^ϵ_I is given explicitly by
  
      K^ϵ_I = v⃗(K_S - 1-e^2/2e^2[(e^2-s_j'^2)K_D_0 - 2s_j'v⃗K_D_1 - v⃗^2K_D_2] ) ,
  
  where the outermost v⃗ term arises due to the change of
  integration variable from s' to α. In turn, the terms
  K_S,K_D_0,K_D_1,K_D_2
  are given by
  
      K_S   = +C_0,1T_0,-1 + 1(C_0,3T_0,-3 + C_1,3T_1,-3 + C_2,3T_2,-3) ,
      K_D_0   = -C_0,1T_0,-3 + 3(C_0,3T_0,-5 + C_1,3T_1,-5 + C_2,3T_2,-5) ,
      K_D_1   = -C_0,1T_1,-3 + 3(C_0,3T_1,-5 + C_1,3T_2,-5 + C_2,3T_3,-5) ,
      K_D_2   = -C_0,1T_2,-3 + 3(C_0,3T_2,-5 + C_1,3T_3,-5 + C_2,3T_4,-5) .
  
  Finally, the coefficients C_0,1,C_0,3,C_1,3,C_2,3 are determined by
  the choice of Taylor expansion point, being either the left or right endpoint
  of the segment. When expanding about the left endpoint, where the shifted
  rescaled arclength parameter is s_j', we have
  
      C_0,1   = I ,
      C_0,3   = χ(s_j')I + w⃗w⃗^T ,
      C_1,3   = v⃗χs'(s_j')I + w⃗v⃗^T + v⃗w⃗^T ,
      C_2,3   = 1/2v⃗^2^2χs'^2(s_j')I + v⃗v⃗^T ,
  
  with v⃗ as defined previously. Here, w⃗ joins the evaluation
  point to the left endpoint of the segment, which, in the case of
  <ref>, is given as w⃗ = x⃗^S(s_i',ϕ)
  - x⃗_j but may be readily generalised to evaluation points off the
  surface of the filament. The corresponding expressions for expansion about the
  right endpoint are
  
      C_0,1   = I ,
      C_0,3   = (χ(s_j+1') - v⃗χs'(s_j+1') + 1/2v⃗^2^2χs'^2(s_j+1'))I + w⃗w⃗^T ,
      C_1,3   = (v⃗χs'(s_j+1') - v⃗^2^2χs'^2(s_j+1'))I + w⃗v⃗^T + v⃗w⃗^T ,
      C_2,3   = 1/2v⃗^2^2χs'^2(s_j+1')I + v⃗v⃗^T .
  
  
  
  
  § EXPLICIT ANTIDERIVATIVES
  
  
  Writing β=β(α)=B+2Cα for brevity, the antiderivatives of
  α^mR^q for m=0, q∈{-1,-3,-5} may be readily computed as
  
      ∫ R^-1α   = C^-1/2log(β + 2C^1/2R(α)) ,
      ∫ R^-3α   = -2/B^2-4ACβ/R(α) ,
      ∫ R^-5α   = -2/3(B^2-4AC)^2(B^2-8BCα-4C(3A+2Cα^2))β/R(α)^3 ,
  
  unless we are in the degenerate case, where B^2-4AC=0, which yields
  
      ∫ R^-1α   = C^-1/2(β)log(β) ,
      ∫ R^-3α   = -2C^1/2(β)/β^2 ,
      ∫ R^-5α   = -4C^3/2(β)/β^4 .
  
  Here we have assumed that C>0, consistent with our assumptions on the
  derivatives of χ and the definition of C in <ref>. The analysis
  of <cit.> and the assumptions of <ref> are
  sufficient to guarantee that the R(α) is nonzero on α∈[0,1],
  thus these integrals are indeed well defined.
  
  
  
  § CONVERGENCE OF SURFACE VELOCITY
  
  
  For the radius function in <ref>b, we compute the error in the
  surface velocity of a straight filament in unit background flow using the RNS
  methodology, as in <ref> though here
  sampling at 2000 points on the surface. The maximum error over the filament
  surface is reported in <ref>, showing substantial refinement as
  N increases for common values of slenderness parameter ϵ. Similar
  to the method of regularised Stokeslet segments
  <cit.>, for regimes with both large ϵ and N
  we see that the error increases dramatically, such that the method is highly
  inaccurate, occurring when the parameters are approximately past the threshold
  ϵ√(N) = 1, which is illustrated as a black dashed line in
  <ref>, though this relation is purely empirical. Notably, this
  typical breakdown occurs outside regimes of common relevance. Regions marked
  with crosses correspond to errors larger than the range of the colour axis.
  
  
  
  32
  natexlab#1#1#1#1 #1#1 #1#1#1#1#1#1
    #1#1#1#1 #1#1 #1#1
    #1#1#1#1#1#1#1#1
  
  [Ainley et al.(2008)Ainley, Durkin, Embid, Boindala &
    Cortez]Ainley2008
  Ainley, Josephine, Durkin, Sandra, Embid, Rafael,
    Boindala, Priya & Cortez, Ricardo 2008  The method of
    images for regularized Stokeslets.  Journal of Computational Physics
    227 (9),  4600–4616.
  
  [Chwang & Wu(1975)]Chwang1975
  Chwang, Allen T. & Wu, T. Yao-Tsu 1975
    Hydromechanics of low-Reynolds-number flow. Part 2. Singularity method
    for Stokes flows.  Journal of Fluid Mechanics  67 (4),
    787–815.
  
  [Cortez(2001)]Cortez2001
  Cortez, Ricardo 2001  The method of regularized
    Stokeslets.  SIAM Journal on Scientific Computing  23 (4),
    1204–1225.
  
  [Cortez(2018)]Cortez2018
  Cortez, Ricardo 2018  Regularized Stokeslet segments.
    Journal of Computational Physics  375,  783–796.
  
  [Cortez & Nicholas(2012)]Cortez2012
  Cortez, Ricardo & Nicholas, Michael 2012  Slender
    body theory for Stokes flows with regularized forces.  Communications
    in Applied Mathematics and Computational Science  7 (1),
    33–62.
  
  [Cox(1970)]Cox1970
  Cox, R. G. 1970  The motion of long slender bodies in a
    viscous fluid Part 1. General theory.  Journal of Fluid Mechanics
    44 (04),  791–810.
  
  [Curtis et al.(2012)Curtis, Kirkman-Brown, Connolly &
    Gaffney]Curtis2012
  Curtis, M. P., Kirkman-Brown, J. C., Connolly, T. J. &
    Gaffney, E. A. 2012  Modelling a tethered mammalian sperm
    cell undergoing hyperactivation.  Journal of Theoretical Biology
    309,  1–10.
  
  [Gillies et al.(2009)Gillies, Cannon, Green &
    Pacey]Gillies2009
  Gillies, Eric A., Cannon, Richard M., Green, Richard B. &
    Pacey, Allan A. 2009  Hydrodynamic propulsion of human
    sperm.  Journal of Fluid Mechanics  625,  445–474.
  
  [Gray(1928)]Gray1928
  Gray, James 1928 Ciliary movement.  Cambridge
    [England]: Cambridge University Press.
  
  [Gray & Hancock(1955)]Gray1955
  Gray, J. & Hancock, G. J. 1955  The Propulsion of
    Sea-Urchin Spermatozoa.  Journal of Experimental Biology
    32 (4),  802–814.
  
  [Guglielmini et al.(2012)Guglielmini, Kushwaha, Shaqfeh &
    Stone]Guglielmini2012
  Guglielmini, Laura, Kushwaha, Amit, Shaqfeh, Eric S. G. &
    Stone, Howard A. 2012  Buckling transitions of an elastic
    filament in a viscous stagnation point flow.  Physics of Fluids
    24 (12),  123601.
  
  [Hall-Mcnair et al.(2019)Hall-Mcnair, Montenegro-Johnson,
    Gadêlha, Smith & Gallagher]Hall-McNair2019
  Hall-Mcnair, Atticus L., Montenegro-Johnson, Thomas D.,
    Gadêlha, Hermes., Smith, David J. & Gallagher,
    Meurig T. 2019  Efficient implementation of elastohydrodynamics
    via integral operators.  Physical Review Fluids  4 (11),
    1–24.
  
  [Hancock(1953)]Hancock1953
  Hancock, G. J. 1953  The self-propulsion of microscopic
    organisms through liquids.  Proceedings of the Royal Society of London.
    Series A. Mathematical and Physical Sciences  217 (1128),
    96–121.
  
  [Ishimoto & Gaffney(2018)]Ishimoto2018a
  Ishimoto, Kenta & Gaffney, Eamonn A. 2018  An
    elastohydrodynamical simulation study of filament and spermatozoan swimming
    driven by internal couples.  IMA Journal of Applied Mathematics
    83 (4),  655–679.
  
  [Johnson(1980)]Johnson1980
  Johnson, Robert E. 1980  An improved slender-body theory
    for Stokes flow.  Journal of Fluid Mechanics  99 (2),
    411–431.
  
  [Keller & Rubinow(1976)]keller1976
  Keller, Joseph B & Rubinow, Sol I 1976  Slender-body
    theory for slow viscous flow.  Journal of Fluid Mechanics
    75 (4),  705–714.
  
  [Lighthill(1976)]lighthill1976
  Lighthill, James 1976  Flagellar hydrodynamics.
    SIAM review  18 (2),  161–230.
  
  [Moreau et al.(2018)Moreau, Giraldi &
    Gadêlha]Moreau2018
  Moreau, Clément, Giraldi, Laetitia & Gadêlha,
    Hermes 2018  The asymptotic coarse-graining formulation of
    slender-rods, bio-filaments and flagella.  Journal of The Royal Society
    Interface  15 (144),  20180235.
  
  [Neal et al.(2020)Neal, Hall-McNair, Kirkman-Brown, Smith &
    Gallagher]Neal2020
  Neal, Cara V., Hall-McNair, Atticus L., Kirkman-Brown,
    Jackson, Smith, David J. & Gallagher, Meurig T. 2020
    Doing more with less: The flagellar end piece enhances the propulsive
    effectiveness of human spermatozoa.  Physical Review Fluids
    5 (7),  073101.
  
  [Olson et al.(2013)Olson, Lim & Cortez]Olson2013
  Olson, Sarah D., Lim, Sookkyung & Cortez, Ricardo
    2013  Modeling the dynamics of an elastic rod with intrinsic
    curvature and twist using a regularized Stokes formulation.  Journal of
    Computational Physics  238,  169–187.
  
  [Pozrikidis(1992)]pozrikidis1992
  Pozrikidis, Constantine 1992 Boundary integral and
    singularity methods for linearized viscous flow.  Cambridge
    University Press.
  
  [Pozrikidis(2010)]Pozrikidis2010
  Pozrikidis, C. 2010  Shear flow over cylindrical rods
    attached to a substrate.  Journal of Fluids and Structures
    26 (3),  393–405.
  
  [Roper et al.(2006)Roper, Dreyfus, Baudry, Fermigier, Bibette
    & Stone]Roper2006
  Roper, Marcus, Dreyfus, Rémi, Baudry, Jean,
    Fermigier, M., Bibette, J. & Stone, H. A. 2006  On
    the dynamics of magnetically driven elastic filaments.  Journal of
    Fluid Mechanics  554,  167–190.
  
  [du Roure et al.(2019)du Roure, Lindner, Nazockdast &
    Shelley]DuRoure2019
  du Roure, Olivia, Lindner, Anke, Nazockdast, Ehssan N. &
    Shelley, Michael J. 2019  Dynamics of Flexible Fibers in
    Viscous Flows and Fluids.  Annual Review of Fluid Mechanics
    51 (1),  539–572.
  
  [Schoeller & Keaveny(2018)]Schoeller2018
  Schoeller, Simon F. & Keaveny, Eric E. 2018  From
    flagellar undulations to collective Motion: Predicting the dynamics of sperm
    suspensions.  Journal of the Royal Society Interface  15 (140),
     arXiv: 1801.08180.
  
  [Shampine & Reichelt(1997)]Shampine1997
  Shampine, Lawrence F. & Reichelt, Mark W. 1997  The
    MATLAB ODE Suite.  SIAM Journal on Scientific Computing
    18 (1),  1–22.
  
  [Simons et al.(2015)Simons, Fauci & Cortez]Simons2015
  Simons, Julie, Fauci, Lisa & Cortez, Ricardo 2015
    A fully three-dimensional model of the interaction of driven elastic
    filaments in a Stokes flow with applications to sperm motility.
    Journal of Biomechanics  48 (9),  1639–1651.
  
  [Smith(2009)]Smith2009d
  Smith, D. J. 2009  A boundary element regularized
    Stokeslet method applied to cilia- and flagella-driven flow.
    Proceedings of the Royal Society A: Mathematical, Physical and
    Engineering Sciences  465 (2112),  3605–3626,  arXiv:
    1008.0570.
  
  [Smith et al.(2019)Smith, Montenegro-Johnson &
    Lopes]Smith2019
  Smith, David J., Montenegro-Johnson, Thomas D. & Lopes,
    Susana S. 2019  Symmetry-Breaking Cilia-Driven Flow in
    Embryogenesis.  Annual Review of Fluid Mechanics  51 (1),
    105–128.
  
  [Walker et al.(2020)Walker, Curtis, Ishimoto &
    Gaffney]Walker2020b
  Walker, Benjamin J., Curtis, Mark P., Ishimoto, Kenta &
    Gaffney, Eamonn A. 2020  A regularised slender-body theory of
    non-uniform filaments.  Journal of Fluid Mechanics  899,
    A3.
  
  [Walker et al.(2019a)Walker, Ishimoto,
    Gadêlha & Gaffney]Walker2019d
  Walker, Benjamin J., Ishimoto, Kenta, Gadêlha, Hermes
    & Gaffney, Eamonn A. 2019a  Filament
    mechanics in a half-space via regularised Stokeslet segments.  Journal
    of Fluid Mechanics  879,  808–833.
  
  [Walker et al.(2019b)Walker, Ishimoto &
    Gaffney]Walker2019g
  Walker, Benjamin J., Ishimoto, Kenta & Gaffney, Eamonn A.
    2019b  A new basis for filament simulation in
    three dimensions ,  arXiv: 1907.04823.
  
  
  
  
  