Regularized Stokeslet rings - an efficient method for axisymmetric Stokes flow, with application to the growing pollen tube
James Tyrrell David J Smith Rosemary J Dyson





C>c<=.08em
=.05em
=.03em
=.65ex
=0pt
=.4ex
=0pt
==.5em
=.5em


"


	style              = Matlab-editor,
	basicstyle         = ,
	escapechar         = ",
	mlshowsectionrules = true,
automata,topaths
mdetect-weight=true, detect-family=true


arrows.metaOMXMnSymbolEMnLargeSymbolsOMXMnSymbolEmnMnLargeSymbolsboldOMXMnSymbolEbnOMXMnSymbolEmn
	<-6>  MnSymbolE5
	<6-7>  MnSymbolE6
	<7-8>  MnSymbolE7
	<8-9>  MnSymbolE8
	<9-10> MnSymbolE9
	<10-12> MnSymbolE10
	<12->   MnSymbolE12
OMXMnSymbolEbn
	<-6>  MnSymbolE-Bold5
	<6-7>  MnSymbolE-Bold6
	<7-8>  MnSymbolE-Bold7
	<8-9>  MnSymbolE-Bold8
	<9-10> MnSymbolE-Bold9
	<10-12> MnSymbolE-Bold10
	<12->   MnSymbolE-Bold12
undefined
undefined

MnLargeSymbols'164MnLargeSymbols'164
MnLargeSymbols'171MnLargeSymbols'171-
[Corresponding author:	]r.j.dyson@bham.ac.ukUniversity of Birmingham, School of Mathematics, Edgbaston, Birmingham, B15 2TT, United Kingdom
The method of regularized Stokeslets, based on the divergence-free exact solution to the equations of highly viscous flow due to a spatially-smoothed concentrated force, is widely employed in biological fluid mechanics. Many problems of interest are axisymmetric, motivating the study of the azimuthally-integrated form of the Stokeslet which physically corresponds to a ring of smoothed forces. The regularized fundamental solution for the velocity (single layer potential) and stress (double layer potential) due to an axisymmetric ring of smoothed point forces, the `regularized ringlet', is derived in terms of complete elliptic integrals of the first and second kind. The relative errors in the total drag and surrounding fluid velocity for the resistance problem on the translating, rotating unit sphere, as well as the condition number of the underlying resistance matrix, are calculated; the regularized method is also compared to 3D regularized Stokeslets, and the singular method of fundamental solutions. The velocity of Purcell's toroidal swimmer is calculated; regularized ringlets enable accurate evaluation of surface forces and propulsion speeds for non-slender tori. The benefits of regularization are illustrated by a model of the internal cytosolic fluid velocity profile in the rapidly-growing pollen tube. Actomyosin transport of vesicles in the tube is modelled using forces immersed in the fluid, from which it is found that transport along the central actin bundle is essential for experimentally-observed flow speeds to be attained. The effect of tube growth speed on the internal cytosolic velocity is also considered. For axisymmetric problems, the regularized ringlet method exhibits a comparable accuracy to the method of fundamental solutions whilst also allowing for the placement of forces inside of the fluid domain and having more satisfactory convergence properties. 
Regularized Stokeslet rings – an efficient method for axisymmetric Stokes flow, with application to the growing pollen tube
    R.J. Dyson
    December 30, 2023
===========================================================================================================================



§ INTRODUCTION

The Stokes equations for incompressible flow at zero Reynolds number are used extensively to model the viscous-dominated regime of microscale flows, particularly biological flows associated with cilia-driven transport, and the motility and feeding of flagellated cells such as bacteria, spermatozoa, algae and choanoflagellates. For an overview, see Lauga and Powers <cit.>. The fundamental solution of the Stokes flow equation, which corresponds to the flow driven by a single spatially-concentrated force is often referred to as the Oseen tensor or Stokeslet. The linearity of the Stokes flow equations enables the construction of solutions to problems involving moving boundaries with complex geometry through integral sums of Stokeslets, forming the basis for the method of fundamental solutions, slender body theory and boundary integral methods. The latter numerical method has the principle major advantage of avoiding the need to mesh the fluid volume, which has enabled highly accurate and efficient simulation of biological flow systems for several decades <cit.>. Indeed more approximate methods based on line distributions of Stokeslets and higher order singularities also enabled major progress in this area before the present era of computationally-intensive research. For review see the earlier work of Chwang, Wu and co-authors
<cit.>, who also explore a wide range of applications as part of a series of papers on low-Reynolds number flow <cit.>.

Nevertheless, two implementational issues arise with methods based on singular solutions. The first is that boundary integrals of solutions with a 1/r-type singularity can be technically complex to evaluate on or near the boundary. Line integrals associated with models of slender bodies such as cilia and flagellar are `more singular', and can require careful distinction between the inside and outside of the body. Moreover, there are cases in which immersed forces due to, e.g. many suspended moving particles, are desired to be modelled by an immersed volumetric force. Cortez et al. developed the method of regularized Stokeslets <cit.> based on the exact, divergence-free solution to the Stokes flow equations due to a concentrated but spatially-smoothed (regularized) force. This approach has enabled the use of Stokeslet methods in a wider range of applications, such as those in which an inducing force is present in  the interior of the fluid domain (as either a point in ℝ^2 or a point/curve in ℝ^3). 

While conceptually elegant, the standard implementation of the method of regularized Stokeslets is computationally expensive, motivating the development of boundary element discretization <cit.>, line integration <cit.>, and meshless interpolation <cit.> among other approaches. Many diverse biological flow problems of interest exhibit rotational symmetry, examples including spherical `squirmer' swimmers <cit.>, the conceptual toroidal swimmer of Purcell <cit.>, and cytosolic flow in elongating pollen tubes <cit.>. Thus motivated, in this paper we study axisymmetric Stokes flows in which the singular and regularized Stokeslets can be integrated azimuthally to yield an axisymmetric ring of point forces. The singular solution to this problem is already known <cit.>; the regularized solution, which we term the `regularized ringlet', is newly derived. This solution forms the basis for an efficient axisymmetric method of regularized Stokeslets.

We begin in Section <ref> by introducing the singular and regularized Stokeslets, and review their application in solving the resistance problem for a rigid body translating in a viscous fluid. In Section <ref> we derive the regularized ringlet, whose analytical solution is given in Section <ref> and is applied to the resistance problem for the translating and rotating sphere in Sections <ref>–<ref>. The double layer potential (DLP), relevant to bodies undergoing volume-changing deformation, is considered in Section <ref>, with analytical evaluation of the azimuthal integral in the DLP given in Appendix <ref>. In Section <ref> we consider the case of Purcell's toroidal swimmer <cit.>, in which the method of regularized ringlets enables the calculation of propulsion speeds which are in excellent agreement with analytical results for both slender and non-slender tori. We then present a study of fluid flow in the angiosperm pollen tube in Section <ref>, an illustrative example highlighting the benefits of being able to apply smoothed forces inside the fluid (two prominent features of regularization). In Section <ref> we conclude with a summary of results and a discussion of future work.



§ SINGULAR AND REGULARIZED STOKESLET SOLUTIONS

For the viscous–dominated very low Reynolds number flow associated with microscopic length scales and slow velocities, incompressible Newtonian flow is well–approximated by the steady Stokes flow equations

    μ∇^2𝐮   = ∇ p - 𝐅, 
    ∇·𝐮   = 0,

where μ is dynamic viscosity, p the pressure, 𝐮 the velocity, and 𝐅 the applied force per unit volume. In the case of a singular force of the form 𝐅(𝐱_0)=𝐠^pδ(𝐱_0-𝐱) for arbitrary point force 𝐠^p, arbitrary point 𝐱 at which the singularity is located, and where δ is the Dirac delta function, the fundamental solution <cit.> for 𝐮 is given by

    u_i(𝐱_0)=1/8πμS_ij(𝐱_0,𝐱)g^p_j,

where

    S_ij(𝐱_0,𝐱)=δ_ij/|𝐱_0-𝐱| + (x_0,i-x_i)(x_0,j-x_j)/|𝐱_0-𝐱|^3,

is known as the Stokeslet. 

The singularity in the Stokeslet solution can be eliminated without loss of incompressibility by regularization of the force 𝐅, as described by Cortez and colleagues <cit.>. The Dirac delta function is replaced with 𝐅(𝐱_0)=𝐠^pϕ_ε(𝐱_0-𝐱) where ϕ_ε is a radially symmetric, smooth `cutoff' function with the property ∫_ℝ^3ϕ_ε(𝐱) d𝐱=1. This is in essence applying the force over a small ball, varying smoothly from a maximum at its centre to ≈ 0 sufficiently far away, instead of using an infinite point force as in the classical Stokeslet solution. The numerical parameter ε dictates the radius of support of the force, and as ε→ 0 the classical solution is recovered. Solutions for 𝐮 using regularized Stokeslets differ from those found using the singular Stokeslet only near the point where the force is applied. Following Cortez et al.<cit.> we take

    ϕ_ε(𝐱_0-𝐱)=15ε^4/8π(|𝐱_0-𝐱|^2+ε^2)^7/2,

which yields

    S_ij^ε(𝐱_0,𝐱)=δ_ij|𝐱_0-𝐱|^2+2ε^2/(|𝐱_0-𝐱|^2+ε^2)^3/2+(x_0,i-x_i)(x_0,j-x_j)/(|𝐱_0-𝐱|^2+ε^2)^3/2.


By considering a solid body D moving through the fluid, it can be shown that

    ∫_ℝ^3 u_i(𝐱)ϕ_ε(𝐱_0-𝐱)   dV(𝐱) =  1/8πμ∫_∂ D  S_ij^ε (𝐱_0,𝐱) g_j^a   dS(𝐱),
 
where 𝐠^a is the force per unit area exerted by the body surface (denoted ∂ D) on the surrounding fluid <cit.>. The above equation is exact; replacing the left hand side with the velocity u_i(𝐱_0) introduces an error 𝒪(ϵ^p) where p=1 on or near the body surface, and p=2 sufficiently far away.

Discretising Equation (<ref>) using N Stokeslets on the surface of the solid body D enables the approximation of the fluid velocity at any point 𝐱_0 via a numerical quadrature formula

    u_i(𝐱_0) = 1/8πμ∑_n=1^N ∑_i=1^3 S_ij^ε (𝐱_0,𝐱_n)g^a_n,jW_n,
 where g^a_n,j denotes the jth component of the force per unit area applied at the point 𝐱_n (a Stokeslet location) and W_n is the quadrature weight associated with the nth particle. The value of W_n is dependent on the geometry of the body surface ∂ D and in Cortez and colleagues' work has units of area.



 §.§ Derivation of the regularized ringlet

Consider a specific case of Equation (<ref>) using a cylindrical (r,θ,z) coordinate system in which the body D exhibits rotational symmetry about the z axis. This symmetry enables analytical integration azimuthally, reducing the surface discretisation to a line discretisation and increasing accuracy. In doing so, we are effectively placing `rings' of regularized Stokeslets at positions  𝐱_n = (r_n,θ,z_n) for n=1,...,N, and θ∈[0,2π) (see Figure <ref> for a diagram of a single ring). This is analogous to covering the surface of the body in `strips' instead of the patches used in a standard 3D Cartesian discretisation. With a surface parameterization 𝐱(s,θ) where 0≤ s≤ℓ denotes arclength and 0≤θ≤ 2π, the boundary integral equation (<ref>) reads

    u_i(𝐱_0) = 1/8πμ∫_∂ D  S_ij^ε (𝐱_0,𝐱) g^a_j(𝐱)   dS(𝐱) = 1/8πμ∫_s=0^ℓ(∫_θ=0^2π S_ij^ε (𝐱_0,𝐱(s,θ)) g^a_j(s,θ)   r(s)   dθ)  ds.


Converting to cylindrical polar coordinates, we introduce the transformation matrix 

    Θ(θ) = ([  cosθ -sinθ     0;  sinθ  cosθ     0;     0     0     1 ]).

Letting indices i,j,k and α,β,γ correspond to Cartesian and cylindrical polar bases respectively (with Einstein summation convention employed for both sets), it follows that g_j^a = Θ_jα(θ)g_α^a. Assuming that velocity is evaluated at fluid point 𝐱_0=(r_0,θ_0,z_0) in cylindrical polars, it further follows that u_i = Θ_iα(θ_0)u_α. Recognizing that Θ^-1=Θ^T, substitution of the cylindrical forms of the velocity and force vectors into Equation (<ref>) thus yields

    u_α(r_0,θ_0,z_0)     = Θ_α i(θ_0)u_i(x_0,y_0,z_0) 
        = 1/8πμΘ_α i(θ_0)∫_0^ℓ∫_0^2π S_ij^ε (𝐱_0,𝐱(s,θ)) Θ_j β(θ)g^a_β(s)   r(s)   dθ  ds.

Under the assumption of axisymmetry, it is sufficient to only consider cases y_0=θ_0=0. This results in the Cartesian x-z and the cylindrical polar r-z planes coinciding, such that the transformation matrix with θ_0=0 simply yields the identity matrix and Equation (<ref>) reduces to

    u_α(r_0,z_0)    = 1/8πμ∫_0^ℓ  g_β^a(s)  [δ_α ir(s)∫_0^2π S_ij^ε (𝐱_0,𝐱(s,θ)) Θ_j β(θ)   dθ]ds, 
       = 1/8πμ∫_0^ℓ  g_β^a(s) R^ε_αβ(𝐱_0,𝐱(s))   ds,

in which the `ringlet' kernel

    R^ε_αβ (𝐱_0,𝐱(s)) := δ_α ir(s)∫_0^2π S^ε_ij(𝐱_0,𝐱(s,θ))Θ_jβ(θ)   dθ,
 is implicitly defined. Here, δ_α i is the Kronecker delta function in which δ_α i=1 if α = i and 0 otherwise. Unlike the Stokeslet which has the symmetric property S_ij^ε(𝐱_0,𝐱) ≡ S_ij^ε(𝐱,𝐱_0), ordering of arguments in the ringlet is important; the first and second arguments in R_αβ^ε(𝐱_0,𝐱) denote the fluid point and ring location respectively, with the ring radius r(s) being the crucial non-symmetric term.

The various terms in the Stokeslet S_ij^ε can be evaluated in cylindrical polar coordinates via

    (x_0,1-x_1)     = r_0-rcosθ, 
     (x_0,2-x_2)    = -rsinθ, 
     (x_0,3-x_3)    = z_0-z, 
     |𝐱_0-𝐱|^2    = (r_0-rcosθ)^2+(rsinθ)^2 + (z_0-z)^2.
 
The resulting form of each S^ε_ij is given in Appendix <ref> and used in the evaluation of R_αβ^ε in Equation (<ref>). These integrals R_αβ^ε yield the regularized fundamental solution for an axisymmetric ring of concentrated force (the regularized ringlet). Analytical evaluation reveals that R_rθ^ε=R_zθ^ε=R_θ r^ε=R_θ z^ε=0 such that the rotational problem for g_θ^a decouples from g_r^a and g_z^a. The remaining nonzero R_αβ^ε yield the equations

    ([ u_r(𝐱_0); u_z(𝐱_0) ]) =1/8πμ∫_0^ℓ[[ R_rr^ε(𝐱_0,𝐱(s)) R_rz^ε(𝐱_0,𝐱(s)); R_zr^ε(𝐱_0,𝐱(s)) R_zz^ε(𝐱_0,𝐱(s)) ]]([ g^a_r(s); g^a_z(s) ])   ds,
 and 
    u_θ(𝐱_0) = 1/8πμ∫_0^ℓ R^ε_θθ(𝐱_0,𝐱(s)) g^a_θ(s)   ds.

Utilising Equations (<ref>)–(<ref>) in tandem models axisymmetric problems with or without azimuthal rotation, in which the fluid experiences a constant force in each principal direction (𝐫̂,θ̂,𝐳̂) at points along which the ringlet is located. This could be used, for example, to model the flow around a mobile axisymmetric body rotating about an axis defined by its direction of translation. 

Approximating the integrals in Equations (<ref>)–(<ref>) numerically using a series of N rings yields a system of equations of the form

    ([ u_r(𝐱_0); u_z(𝐱_0) ])    = 1/8πμ∑_n=1^N [[ R_rr^ε(𝐱_0,𝐱_n) R_rz^ε(𝐱_0,𝐱_n); R_zr^ε(𝐱_0,𝐱_n) R_zz^ε(𝐱_0,𝐱_n) ]]([ g^a_r (𝐱_n); g^a_z (𝐱_n) ])w_n, 
     u_θ(𝐱_0)    = 1/8πμ∑_n=1^N R^ε_θθ(𝐱_0,𝐱_n) g^a_θ(𝐱_n)w_n,
 where g^a_r(𝐱_n),g^a_θ(𝐱_n),g^a_z(𝐱_n) are the radial, azimuthal, and axial components of the forces per unit area applied at ringlet location 𝐱_n respectively, and w_n is the quadrature weight associated with 𝐱_n for numerical integration over s. The quantity w_n has units of length unlike its counterpart W_n in Cortez and colleagues' work (units of length squared). It is also possible to recompile the force per unit area 𝐠^a and quadrature weight w_n into a force per unit length 𝐠^l, such that Equation (<ref>) can alternatively represent the fluid velocity induced by a series of rings. 

By considering the fluid velocity at each individual ringlet location, an invertible system can be produced. In the irrotational case (u_θ≡ 0), this takes the block matrix form

    𝖦^l = 8πμ[ [ 𝖱_rr^ε 𝖱_rz^ε; 𝖱_zr^ε 𝖱_zz^ε ]]_𝖱^ε^-1𝖴,
 where

    𝖦^l = [[ g^l_r (𝐱_1);           ⋮; g^l_r (𝐱_N); g^l_z (𝐱_1);           ⋮; g^l_z (𝐱_N) ]],    𝖴 =  [ [ u_r (𝐱_1);         ⋮; u_r (𝐱_N); u_z (𝐱_1);         ⋮; u_z (𝐱_N) ]],
 and

    𝖱_αβ^ε=[[ R_αβ^ε(𝐱_1,𝐱_1) R_αβ^ε(𝐱_1,𝐱_2)               … R_αβ^ε(𝐱_1,𝐱_N); R_αβ^ε(𝐱_2,𝐱_1) R_αβ^ε(𝐱_2,𝐱_2)               … R_αβ^ε(𝐱_2,𝐱_N);               ⋮               ⋮               ⋱               ⋮; R_αβ^ε(𝐱_N,𝐱_1) R_αβ^ε(𝐱_N,𝐱_2)               … R_αβ^ε(𝐱_N,𝐱_N) ]].
 
Hence the forces needed to induce a given prescribed velocity in the fluid may be found (the resistance problem). This works for both a series of translating rings and a translating axisymmetric body, with the force per unit length 𝐠^l essentially `absorbing' both the force per unit area 𝐠^a and the quadrature weight w_n in the latter case. Inclusion of rotation involves the formulation of a similar invertible system for Equation (<ref>) which can be solved separately.



 §.§ Analytical evaluation of the regularized ringlet

The nonzero elements of the regularized ringlet can be expressed in the form

    R_rr^ε(𝐱_0,𝐱_n)    = r_n(-r_0r_nI_0+(2τ-(z_0-z_n)^2)I_1 - 3r_0r_nI_2), 
     R_rz^ε(𝐱_0,𝐱_n)    = r_n(z_0-z_n)(r_0I_0-r_nI_1), 
     R_zr^ε(𝐱_0,𝐱_n)    = r_n(z_0-z_n)(-r_nI_0+r_0I_1), 
     R_zz^ε(𝐱_0,𝐱_n)    = r_n((τ + (z_0-z_n)^2 + ε^2)I_0 - 2r_0r_nI_1), 
     R_θθ^ε(𝐱_0,𝐱_n)    = r_n(r_0r_nI_0 + (τ+ε^2)I_1 - 3r_0r_nI_2),
 in which τ:=r_0^2+r_n^2+(z_0-z_n)^2+ε^2, and 
    I_n    := ∫_0^2πcos^nθ/(τ-2r_0r_ncosθ)^3/2  dθ, 
        = 4k^3/(4r_0r_n)^3/2∫_0^π/2(2cos^2θ-1)^n/(1-k^2cos^2θ)^3/2 dθ,

with k^2 := 4r_0r_n/(τ+2r_0r_n). Equation (<ref>) is found by using the double angle formula for cosθ as well as symmetry arguments about π/2. Following the example of Pozrikidis <cit.>, the integrals I_n can be computed by first expanding the numerator of the integrand in Equation (<ref>) to obtain a series of polynomial integrals with respect to cosθ. Letting

    I_n' := 2^n(4k^3/(4r_0r_n)^3/2) ∫_0^π/2cos^2nθ/(1-k^2cos^2θ)^3/2 dθ,
 
it follows that 
    I_0 = I_0',    I_1 = I_1' - I_0,    I_2 = I_2'-2I_1'+I_0.
 The individual integrals I_n' can be expressed in terms of complete elliptic integrals of the first and second kind, which are respectively defined

    F = F(k) := ∫_0^π/2dθ/(1-k^2sin^2θ)^1/2  and   E = E(k) := ∫_0^π/2(1-k^2sin^2θ)^1/2  dθ.
 The solutions for each I_n' (as can be found in Section 2.58 of Gradshteyn and Ryzhik <cit.>) are given by

    I_0'    = 4k^3/(4r_0r_n)^3/2(1/1-k^2E) , 
    
      I_1'    = 8k^3/(4r_0r_n)^3/2(1/k^2(1-k^2)E - 1/k^2F), 
    
      I_2'    = 16k^3/(4r_0r_n)^3/2(2-k^2/k^4(1-k^2)E-2/k^4F),

from which it follows that

    I_0    = 4k^3/(4r_0r_n)^3/2(1/1-k^2E) , 
    
      I_1    = 4k^3/(4r_0r_n)^3/2(2-k^2/k^2(1-k^2)E - 2/k^2F), 
    
      I_2    = 4k^3/(4r_0r_n)^3/2(k^4-8k^2+8/k^4(1-k^2)E-4(2-k^2)/k^4F).

Substitution of Equations (<ref>)–(<ref>) back into Equations (<ref>)–(<ref>) yields the complete solution for the regularized ringlet

    R_rr^ε(𝐱_0,𝐱_n)     = k/r_0r_n(r_n/r_0)^1/2[(τ+(z_0-z_n)^2)F + 4r_0^2r_n^2 - τ(τ+(z_0-z_n)^2)/τ-2r_0r_nE], 
    
    R_rz^ε(𝐱_0,𝐱_n)     = k(z_0-z_n)/r_0(r_n/r_0)^1/2[F + 2r_0^2-τ/τ-2r_0r_nE], 
     R_zr^ε(𝐱_0,𝐱_n)     = -k(z_0-z_n)/(r_0r_n)^1/2[F + 2r_n^2-τ/τ-2r_0r_nE], 
     R_zz^ε(𝐱_0,𝐱_n)     = 2k(r_n/r_0)^1/2[F + (z_0-z_n)^2+ε^2/τ-2r_0r_nE], 
     R_θθ^ε(𝐱_0,𝐱_n)     = k/r_0r_n(r_n/r_0)^1/2[(2τ-ε^2)F + 8r_0^2r_n^2+τ(ε^2-2τ)/τ-2r_0r_nE].

The solutions given by Equations (<ref>)–(<ref>) can be readily evaluated except when r_0=0 or r_n=0. In the limit as r_n→ 0, all R_αβ^ε→ 0. In the limit as r_0 → 0, both R_zr^ε and R_zz^ε tend to finite values while R_rr^ε,R_rz^ε,R_θθ^ε→ 0. This behaviour is described in detail in Appendix <ref>. Streamlines for the flow induced by the regularized ringlet with associated unit forces in both 𝐫̂ and 𝐳̂ directions are given in Appendix <ref>.

The form of the ringlet solutions R_αβ^ε is similar to those for the ring of singular Stokeslets (as detailed by Pozrikidis in <cit.> with the exception of the newly derived R_θθ^0)

    R_rr^0(𝐱_0,𝐱_n)      =  k/r_0r_n(r_n/r_0)^1/2[(r_0^2+r_n^2+2(z_0-z_n)^2)F - 2(z_0-z_n)^4 +3(z_0-z_n)^2(r_0^2+r_n^2) + (r_0^2-r_n^2)^2/(z_0-z_n)^2+(r_0-r_n)^2E], 
     R_rz^0(𝐱_0,𝐱_n)    =  k(z_0-z_n)/r_0(r_n/r_0)^1/2[F + r_0^2 - r_n^2 - (z_0-z_n)^2/(z_0-z_n)^2 + (r_0-r_n)^2E], 
     R_zr^0(𝐱_0,𝐱_n)    =  -k(z_0-z_n)/(r_0r_n)^1/2[F - r_0^2 - r_n^2 + (z_0-z_n)^2/(z_0-z_n)^2 + (r_0-r_n)^2E], 
     R_zz^0(𝐱_0,𝐱_n)    = 2k(r_n/r_0)^1/2[F+(z_0-z_n)^2/(z_0-z_n)^2 + (r_0-r_n)^2E], 
     R_θθ^0(𝐱_0,𝐱_n)     = k/r_0r_n(r_n/r_0)^1/2[2(r_0^2+r_n^2+(z_0-z_n)^2)F  - 4(z_0-z_n)^4 + 4(z_0-z_n)^2(r_0^2+r_n^2) + 2(r_0^2-r_n^2)^2/(z_0-z_n)^2 + (r_0-r_n)^2E],
 and in the limit as ε→ 0 our solutions are equivalent to their singular counterparts. This can be verified by substitution of ε=0 into R_αβ^ε and is a result of the cutoff function ϕ_ε approaching a delta distribution as ε→ 0. 

Equations (<ref>)–(<ref>) and (<ref>)–(<ref>) provide the exact analytical solution for the fluid velocity at any point due to the drag force per unit area on the surface of a generalized axisymmetric body. Using Equations (<ref>)–(<ref>) in place of (<ref>)–(<ref>) yields the numerical solution based on discretization over the arclength s. In the case of a single ring, removing the integral over s and replacing the force per unit area 𝐠^a with a force per unit length 𝐠^ℓ yields the exact analytical solution for the fluid velocity induced by the force acting along the ring in 3D space.


 §.§ On the double layer potential

A more complete formulation of Equation (<ref>) for the fluid velocity induced by a translating body D is given by

    ∫_ℝ^3 u_j(𝐱)ϕ_ε(𝐱_0-𝐱)   dV(𝐱) = 1/8πμ∫_∂ D  S_ij^ε (𝐱_0,𝐱) g_i^a(𝐱)   dS(𝐱) + 1/8π∫_∂ D u_i(𝐱)T_ijk^ε(𝐱_0,𝐱)n_k(𝐱)   dS(𝐱),
 
where the first and second integrals on the right hand side are known as the single layer potential (SLP) and the double layer potential (DLP) respectively. The stress tensor T_ijk^ε present in the DLP is given by

    T_ijk^ε(𝐱_0,𝐱) = -6(x_0,i-x_i)(x_0,j-x_j)(x_0,k-x_k)/(|𝐱_0-𝐱|^2 + ε^2)^5/2 - 3ε^2 (x_0,i-x_i)δ_jk + (x_0,j-x_j)δ_ik + (x_0,k-x_k)δ_ij/(|𝐱_0-𝐱|^2 + ε^2)^5/2.


The DLP can be neglected for problems in which the condition ∫_∂ D𝐮·𝐧̂  dS =0 is satisfied <cit.>, which is the case throughout this manuscript. However we provide the expressions for the DLP in Appendix <ref>, which may be of value for future studies of systems in which the condition ∫_∂ D𝐮·𝐧̂  dS =0 is violated (e.g. bubbles).



§ SIMPLE EXAMPLES AND TEST CASES

In Section <ref> the expression for the regularized fundamental solution for an axisymmetric ring of concentrated forces, the regularized ringlet, was derived. In the following, we demonstrate the validity of the method through application to simple cases of motion. 

The first two cases concern the translation and rotation of the unit sphere in a Stokesian fluid, treated independently in Sections <ref> and <ref> respectively. In Section <ref> a more complicated example is considered: the propulsion of `Purcell's toroidal swimmer' <cit.>, powered by tank treading of the torus surface. In considering these different cases, it is shown that the method of regularized ringlets can be used to model the surface motion of axisymmetric bodies in each principal direction 𝐫̂, θ̂, 𝐳̂ in a cylindrical coordinate system.



 §.§ Translating unit sphere

The validity of the regularized ringlet method is illustrated by solving the resistance problem for the translating unit sphere. Given a prescribed surface velocity (-𝐳̂), Equation (<ref>) yields a Fredholm first kind integral equation for the unknown force distribution <cit.>. The method of regularized ringlets (implemented here via Matlab) can be used to solve this problem.

The sphere is parametrized in the rz plane by 𝐩 = cosφ 𝐫̂ + sinφ 𝐳̂ for φ∈ [-π/2,π/2], then discretized as 
    φ_n = πn-1/2/N-π/2  for  n=1,…,N.
 The velocity boundary condition 𝐮=-𝐳̂ is prescribed at each 𝐱_n:=𝐩(φ_n), and the resulting linear system is solved to yield the required force densities 𝐠^ℓ at each of these locations. The drag is then calculated as 
    ∑_n=1^N ∫_θ=0^2π g_z^ℓ(𝐱_n)   r(𝐱_n)dθ = 2π∑_n=1^N r(𝐱_n)g_z^ℓ(𝐱_n),
 which is compared with the Stokes law value of -6π. The relative errors are given for various values of N and regularization parameter ε in Table <ref> alongside the condition number of the resistance matrix 𝖱^ε in Table <ref>.



For given N, excessively small ε results in the drag error becoming non–monotonic. For given ε, increasing N eventually ceases to result in a further reduction in the relative error. This is often the case with regularized Stokeslet methods (eg. see Figure <ref>, Figure <ref>, and references Cortez et al.<cit.>, Gallagher et al.<cit.>). 

A thorough comparison of our results for the translating unit sphere with those of Cortez et al.<cit.> can be found in Appendix <ref>, in which it is found that using just N ringlets yields consistently more accurate solutions than using 3N^2 regularized Stokeslets in a standard 3D discretization of the sphere surface. A comparison to results obtainable using the axisymmetric method of fundamental solutions (a singular Stokeslet method) is given in Appendix <ref>, in which the relative error in the fluid velocity is also discussed. In general, although the singular method can be tuned to give smaller relative errors in either the fluid velocity or total drag separately, it cannot do so simultaneously; regularized ringlets display more satisfactory convergence properties and are the more effective method to minimize errors in both fluid velocity and total drag. 


 §.§ Rotating unit sphere

The solution for the steady motion of a Stokesian fluid surrounding a solid sphere rotating uniformly about a central axis is well known and can be found in e.g. Hydrodynamics by Lamb <cit.>. If the sphere rotates around its z axis in an (r,θ,z) cylindrical coordinate system with an angular velocity Ω = ω_0 𝐞̂_z, the resulting angular velocity of the fluid is given by ω = (a/γ)^3ω_0𝐞̂_θ where γ = √(r^2 + z^2). This can be written in terms of the linear velocity (more readily usable in the Stokeslet formulae) over the entire domain as, 
    𝐮 =  r(a/γ)^3ω_0𝐞̂_θ   ∀γ≥ a, 
     rω_0𝐞̂_θ   ∀γ < a,
 where γ≥ a corresponds to the surrounding fluid velocity and γ < a to the solid body rotation of the sphere respectively. The zero Reynolds number torque on this sphere is given by, 
    𝐓=-8πμ a^3Ω,
 the derivation of which can be found in <cit.>. This torque is associated with a drag force per unit area on the surface of the sphere given by -3μω_0(r/a)𝐞̂_θ, as detailed in Appendix <ref>. The force per unit length used in the method of regularized ringlets thus takes form 
    g^l_θ = (3πμ/N)(r/a)ω_0,

although for the resistance problem this is not prescribed. 

The sphere surface is again parametrized in the rz plane by 𝐩=cosφ 𝐫̂+sinφ 𝐳̂ for φ∈[-π/2,π/2], discretized using N ringlets at locations 𝐱_n. 
Letting ω_0=-1, the velocity 𝐮(𝐱_n)=-r(𝐱_n)θ̂ is prescribed at each 𝐱_n and the resistance matrix is constructed to yield the required force densities 𝐠^ℓ at each of these locations. The torque is then calculated as

    ∑_n=1^N ∫_θ=0^2π g_θ^ℓ(𝐱_n)   r^2(𝐱_n)dθ = 2π∑_n=1^N r^2(𝐱_n)g_θ^ℓ(𝐱_n),
 in which r(𝐱_n) is squared since the torque is a moment, the product of distance and force. Comparing to the value -8π, relative errors are given for various values of N and regularization parameter ε in Table <ref> alongside the condition number of the resistance matrix 𝖱_θ^ε in Table <ref>.


These results are similar to those of the translating unit sphere, although for given ε, N the relative error in the drag is generally slightly larger and the condition number of the resistance matrix slightly smaller than for the results using the same ε,N on the translating unit sphere. 


 §.§ Purcell's toroidal swimmer

The torus is the simplest geometry capable of describing self–propelled organisms <cit.>. Purcell's toroidal swimmer <cit.> describes one such organism, the geometry of which can be seen in Figure <ref>. Inward rotation of the torus surface produces a net force in the direction of motion of the outermost surface (against which the torus is propelled). The magnitude of this net force (and resultant propulsion speed of the torus) is dependent on the speed with which the surface of the torus rotates as well as the slenderness of the torus. It has been suggested that this mechanism could describe how a DNA miniplasmid could be turned into a self–propelled nanomachine <cit.>. Three modes of locomotion are considered by Leshansky and Kenneth <cit.>, corresponding to tank treading of (a) an incompressible surface, in which the tangential surface velocity is largest on the inner surface, (b) a weakly compressible surface, in which the tangential surface velocity is constant, (c) a highly compressible surface, in which the tangential surface velocity is largest on the outer surface. In the following, we restrict ourselves to looking at the case of constant tangential surface velocity. 



The torus geometry is reduced to a single slenderness parameter s_0 = b/a, where b and a refer to the major and minor radii of the torus respectively. The torus surface is parametrized in terms of angle η∈ [0,2π) in the rz plane such that ds = a  dη and η=0 corresponds to the outermost radial point on the surface of the torus, traversed in an anticlockwise direction. For the free swimming torus, the rigid body translation 𝐔 and rotation 𝐮^(s) of the torus surface ∂ D are related to the force per unit area 𝐠^a exerted by the torus on the surrounding fluid by, 
    U_α(𝐱_0) + u^(s)_α(𝐱_0) =  1/8πμ∫_s=0^ℓ  R_αβ^ε(𝐱_0,𝐱) g^a_β(𝐬)   ds   ∀𝐱_0 ∈∂ D
 subject to the condition of zero net force in the 𝐳̂ direction, 
    ∫_∂ D𝐠_z^a(𝐱)   dS(𝐱) = 0,
 where s is the arclength parameterization of the cross section of the torus surface ∂ D in the rz plane. We note that the additional free swimming conditions of zero net force in 𝐫̂ and zero total moment (as outlined by Phan–Thien et al.<cit.>) are automatically satisfied by axisymmetry and g_θ^a ≡ 0 respectively. 

The propulsion speed U of the rotating torus for any given value of s_0 and rotation speed u^(s) can be determined by considering two separate situations: one in which motion is purely translational in the 𝐳̂ direction (the toroidal glider with 𝐔=𝐳̂, 𝐮^(s)≡0) and one in which motion is purely rotational (the anchored toroidal pump with 𝐔≡0, 𝐮^(s) = η̂ where η̂ is the unit vector whose direction varies over s, pointing tangential to the surface in the anticlockwise direction at all points). The glider and pump have associated force distributions 𝐠^gld and 𝐠^pmp respectively. Using the regularized ringlet, surface motions can be prescribed (Figures <ref>&<ref>) in order to compute the associated force distributions (Figures <ref>&<ref>) responsible for producing each motion. For the toroidal glider the radial force is assumed to be zero (g_r^gld≡ 0), hence 
    U_r(𝐱_0)     = 1/8πμ∫_η=0^2π  R_rz^ε(𝐱,𝐱_0) g_z^gld(𝐱)   a  dη = 0   ∀ 𝐱_0 ∈∂ D, 
     U_z(𝐱_0)     = 1/8πμ∫_η=0^2π  R_zz^ε(𝐱,𝐱_0) g_z^gld(𝐱)   a  dη = 1   ∀ 𝐱_0 ∈∂ D,
 which, after solving for the unknown force distribution 𝐠_z^gld, yields the net axial force in the 𝐳̂ direction 
    G^gld = ∫_∂ D  g_z^gld(𝐱)   dS(𝐱),
 which is nonzero. For the anchored toroidal pump the surface velocity is given by 𝐮^(s)=η̂=(-sinη,cosη) in (r,z) coordinates. It follows that 
    u_r^(s)(𝐱_0)   =1/8πμ∫_η=0^2π R_rβ^ε(𝐱,𝐱_0)g_β^pmp(η) a dη = -sinη,   ∀ 𝐱_0 ∈∂ D, 
     u_z^(s)(𝐱_0)   =1/8πμ∫_η=0^2π R_zβ^ε(𝐱,𝐱_0) g_β^pmp(η) a dη =cosη,   ∀ 𝐱_0 ∈∂ D,
 from which it can be determined that the net axial 𝐳̂ force is given by 
    G^pmp = ∫_∂ D g_z^pmp(η)   dS(𝐱).
 By linearity of the Stokes flow equations we may subtract the gliding solution from the pump solution to rewrite the system in the form 
    u_r^(s)(𝐱_0)   =1/8πμ∫_η=0^2π R_rβ^ε(𝐱,𝐱_0)g_β^a(η) a dη,   ∀ 𝐱_0 ∈∂ D, 
     -G^pmp/G^gld + u_z^(s)(𝐱_0)   =1/8πμ∫_η=0^2π R_zβ^ε(𝐱,𝐱_0)g_β^a(η) a dη,   ∀ 𝐱_0 ∈∂ D,
 where 𝐠^a = 𝐠^pmp - (G^pmp/G^gld)g^gld_z𝐳̂ and the net force is equal to 
    ∫_∂ Dg^a_z(η)   dS(𝐱) = 0,
 as required for the free–swimmer. The propulsion velocity of the swimming torus is thus given by 𝐔 = - (G^pmp/G^gld)𝐳̂, opposing the direction of outer surface motion. The speed U=|𝐔| is dependent on both the rotation speed u^(s):=|𝐮^(s)| and the slenderness ratio s_0, so both G^gld and G^pmp must be recomputed whenever one of these parameters is changed. However, Leshansky and Kenneth <cit.> were able to show that the propulsion speed depends linearly on the rotation speed, such that by considering the scaled propulsion speed U/u^(s) it is only necessary vary s_0 to be able to consider all possible propulsion speeds resulting from a given constant rotational surface velocity. 




Our results for the scaled propulsion velocity found using a discretization of the torus surface using N=100 regularized ringlets are compared with those obtained by Leshansky and Kenneth <cit.>, who tackled the same problem using a line distribution of rotlets at the torus centreline (inaccurate as s_0 → 1) and an exact series solution via expansion in toroidal harmonics (Figure <ref>). Of particular note, in the limit as s_0 → 1 it is found that using N=100 rings with ε=0.01 in the regularized ringlet method yields a scaled propulsion velocity of 0.6684, representing just a 0.513% error when compared to the series solution value of 0.665. This error can be reduced to <0.1% by using N=1000 rings, at which point the value of the scaled propulsion velocity as calculated by the method of regularized ringlets is 0.6656. This is a significant improvement over the solution found using a line distribution of rotlets, in which the error is  >1% for all s_0 ≤ 6. 



In addition to giving  values for the scaled propulsion speed that are in excellent agreement with the series solution of Leshansky and Kenneth <cit.>, the regularized ringlet solution also provides the force required at each point on the torus surface to produce the swimming motion. The series and rotlet solutions do not yield this information, with the propulsion speed instead being calculated according to the net drag force on the toroidal glider in these methods. This, in combination with Figures <ref>&<ref>, highlights why the centreline rotlet solution is inaccurate for small values of s_0; as slenderness decreases, the difference between the magnitude of the force required at the inner- and outer-most surfaces to produce rigid body translation of the torus grows large. In Figure <ref> (in which s_0=2), the drag force on the outer surface is approximately 5.7 times larger than that on the inner surface. The constant centreline force associated with the rotlet cannot account for this discrepancy, whereas the full discretization of the torus surface using regularized ringlets can. 

Figures <ref>&<ref> show the streamlines and magnitude of the fluid velocity in a region near the force–free toroidal swimmer with slenderness parameter s_0=2, undergoing uniform counterclockwise surface rotation with unit angular velocity. This results in propagation of the torus in the direction -𝐳̂. Fluid passing through the central hole of the torus is caught in closed streamlines, in agreement with the results of Leshansky and Kenneth <cit.>.




	
	
	
	

	
	
	
	


§ CYTOSOLIC FLOW IN THE POLLEN TUBE - AN ILLUSTRATIVE EXAMPLE OF THE BENEFITS OF REGULARIZATION

In Section <ref>, the method of regularized ringlets was applied to simple cases of motion to illustrate validity and applicability. These canonical examples demonstrate the use of ringlets in situations with moving boundaries. However, perhaps the most important feature of the regularized method is the ability to place source points directly in the fluid, in which case the regularized ring represents a curve in 3D space. The regularization parameter ε can then be used to control the spreading of the force (the size of the region over which it is applied). This is useful in situations in which the exact location at which a force is applied in a fluid is unknown or when a force is applied over a large area. In the following, these features of regularization are explored in more detail by considering the example of fluid flow in the angiosperm pollen tube. A more extensive study of the biomechanics of pollen tube growth (with a particular focus on the transport and distribution of elements of the cytoplasm) will follow in a future paper.



 §.§ A brief overview of the pollen tube

The pollen tube, a cellular protuberance originating from the pollen grain, is a vital component of the fertilisation process in plants. Responsible for the delivery of sperm cells from the pollen grain to the ovule, the tube may have to grow over 30cm in length at speeds in excess of 1cm.h^-1<cit.>. This rapid elongation of the tube occurs via tip growth, with expansion localised close to the apical region of the cell <cit.>. To prevent rupture, new cell membrane and cell wall material must be targeted to the sites of highest expansion. Spherical vesicles in the cytosolic region act as the delivery vectors for this material <cit.> and are carried by motor myosin proteins along actin filaments in a process referred to as actomyosin transport. These actin filaments, arranged in bundles in the periphery and centre of the tube, are oriented in such a way that `full' peripheral vesicles travel towards the growing apex where they secrete their contents into the wall. In order to maintain the correct ratio of structural components, `empty' vesicles (comprised mostly of membrane) are also secreted by the wall back into the cytoplasm, where it is hypothesised that they are picked up by myosin on the central actin bundle and travel away from the apex. The combined movement of these vesicles induces a flow in the cytosol, known as `cytoplasmic streaming' or `cyclosis', which further aids in the cycling of vesicles towards and away from the apical region. To our knowledge, no comprehensive attempt has been made to produce a complete model of cytosolic flow in the pollen tube in a manner deriving directly from physical principles. In this section, we show how this is easily achievable using the method of regularized ringlets. 



 §.§ Mathematical model


  §.§.§ Geometry and boundary conditions

The typical geometry of the tube (in 2D cylindrical coordinates) along with the two vesicle populations and actin bundles can be seen in Figure <ref>, in which the static cylindrical shank is joined to a growing apical hemispherical cap. The boundary of the domain is split into four sections: the symmetry boundary Ω_sym running down the centreline of the tube, the static, impermeable peripheral wall Ω_imp, the growing apical hemisphere Ω_grw and the artificial basal boundary Ω_bas  (where the computational domain is truncated).



The Reynolds number for fluid flow in the pollen tube can be evaluated using the typical flow speed (≈1 m.s^-1), tube radius (≈8.13 m) and the kinematic viscosity of water (≈e6 m^2.s^-1) to find Re≈e-5, firmly in the regime of Stokes flow. Although pollen tube growth is typically oscillatory, acceleration and deceleration are small compared to growth speed itself. This enables the use of the steady Stokes equations.

Since the cytosolic flow is induced by the actomyosin transport of vesicles along cytoskeletal actin bundles, forces must be applied inside the computational domain. Singular stokeslet methods typically do not allow for this unless the surfaces of individual vesicles are discretized, but their small size (≈100nm radius) compared to the typical length scale of the problem and their large number density make this impractical. One solution is to use the method of regularized ringlets instead, placing rings in series along the centreline of the peripheral actin bundle with the parameter ε being used to control the bundle thickness. 

The boundary conditions for the fluid velocity 𝐮 are given by,

    [              𝐮=𝐮_g  on Ω_grw,                 𝐮=0 on Ω_imp, u_r=∂ u_z / ∂ r =0  on Ω_sym, ]

in accordance with the assumptions of tip growth and axisymmetry. No restriction is placed on 𝐮 on the artificial boundary Ω_bas. For the growth velocity 𝐮_g in the hemispherical apex, the normal displacement growth assumption of Dumais <cit.> is employed to define 
    𝐮_g = 𝐮_g(φ)  v_gsinφ(cosφ,sinφ),
 where φ is the angle between the outward-pointing surface normal and the positive r axis (varying from 0 at the point where the hemisphere joins the shank to π/2 at the extreme apex) and v_g is the growth speed of the tube. An example of how the growth of the boundary varies over the apical hemisphere can be seen in Figure <ref>, where the maximum growth speed v_g is equal to 0.1  μm s^-1 and the wall velocity in the adjacent shank (not pictured) is equal to zero. 



  §.§.§ Actin bundle model

A section of the pollen tube is modelled ranging from the extreme apex to a point equal to six tube radii distal where the domain is artificially truncated. In dimensionless values, the extreme apex is thus given by (r,z) = (0,6) with z=0 being the distal truncation line. The central line of axisymmetry is given by r=0, with r=1 denoting the peripheral boundary in the shank. The hemispherical apex is the upper-right quarter circle of radius 1, centred at (0,5). The peripheral actin bundle (PAB) is considered adjacent to the pollen tube wall, with its width equal to one fifth of the pollen tube radius (0.2) and extending at an angle φ = π/5 into the apical hemisphere (see Figure <ref>) in accordance with the confocal microscopy imaging of Lovy-Wheeler et al.<cit.>. The central actin bundle (CAB) is located further away from the apex than the PAB, and is modelled with a reduced thickness (half that of the PAB) in cylindrical coordinates to account for the central axisymmetry of the tube.



By placing rings of regularized Stokeslets in series along the centre line of the PAB, it is possible to carefully select an appropriate value for the regularization parameter ε such that the region over which the majority of the force distribution is applied is roughly the same as the bundle thickness. The same procedure is employed for the CAB, only with individual regularized Stokeslets rather than ringlets since the centreline coincides with r=0. Using ε = 0.05 it is found that ≈93% of the total applied force is contained within a region of radius 0.1 which corresponds well to the approximate thickness of the actin bundles. This is shown in Figure <ref> where the maximal value of ϕ_ε has been scaled to 1 for the sake of clarity. Smaller values of ε result in an even larger percentage of the total force being contained within 0 ≤ρ≤ 0.1 but are increasingly skewed towards ρ=0. 

We note that our modeling ignores the presence of a third f-actin structure, the short actin bundles observed in the extreme apex <cit.>; these bundles are transient and significantly less dense than the thick peripheral and central actin bundles. Further, due to their proximity to the growing tube wall, any effect that these short actin bundles may have on the fluid velocity is likely insignificant compared to the effect of the growth of the wall.



  §.§.§ Frames of reference

In each of the velocity profiles to be produced, cytosolic velocities are prescribed along the actin bundles in accordance with values taken from the STICS analysis of Bove et al.<cit.>. This is done under the assumption that the large-scale directed movement of vesicles must be a consequence of cytosolic flow (induced by the smaller-scale movement of individual vesicles along actin). The fluid velocity is defined in terms of two different coordinate systems: the static lab frame where 𝐮_L=(u_r,u_z) and the moving tip frame with 𝐮_T=(u_r,u_z-v_g). The lab frame represents the fluid velocity relative to a stationary observation point (in which 𝐮_L·𝐧̂≠ 0 on Ω_grw due to boundary growth), whereas the tip frame represents the fluid velocity relative to the growing tip (resulting in a static domain with 𝐮_T·𝐧̂ at all boundary points except on Ω_bas). The steady geometry of the domain for the advancing tip frame means that 𝐮_T satisfies the steady Stokes equations, from which it follows that 𝐮_L is also a solution (since the two differ by the uniform flow field (0,v_g) only).



  §.§.§ Two model scenarios for comparison

In Section <ref>, we investigate the role of the CAB in vesicle transport. Fluid velocity is prescribed on the PAB only (see Figure <ref>) in order to compare the resultant central fluid velocity to experimental results. The centreline of the PAB is modelled using the union of the straight line extending from (r,z) = (0.9,-2) to (0.9,5) and the curve (r,z) = (0.9cosφ,5+0.9sinφ) for φ∈ [0,π/5]. The regularized ringlet placement is extended to z=-2 to ensure the velocity profile at z=0 is consistent with the rest of the tube. The fluid velocity on the straight line is given by (u_r,u_z) = (0,0.5) with the fluid velocity on the curve being tangential and of constant magnitude, that is, (u_r,u_z) = 0.5(sinφ,cosφ) for φ∈[0,π/5]. On the peripheral wall (again extended to z=-2), the velocity is 0 in the shank and prescribed according to the growth speed v_g of the tube in the apical hemisphere using the velocity function 𝐮_g(φ) =v_gsinφ(cosφ,sinφ) for φ∈ [0,π/2]. 

In Section <ref>, the velocity of the growing boundary is varied in order to determine the effect of the growth speed of the tube on the cytosolic velocity. Fluid velocity is prescribed on each of the PAB, the CAB, and the growing boundary accordingly (see Figure <ref>). Fluid velocity along the PAB is kept constant at 0.5 m.s^-1 based on the observation that cytoplasmic streaming rates are typically independent of pollen tube growth rates <cit.>. Fluid velocity on the centreline of the central actin bundle, given by r=0 for -2 ≤ z ≤ 4, is also prescribed. It is assumed that vesicles are sufficiently closely packed on this bundle that the fluid velocity at its centre can be approximated by the speed of vesicles on actin (1 m.s^-1), giving (u_r,u_z) = (0,-1). Since these node locations are at r=0, standard regularized Stokeslets must be used here.

A regularisation parameter of ε=0.05 is chosen in all cases and the ringlets are linearly spaced a distance approximately 0.025 apart, resulting in a smooth velocity profile. The flow velocity for z < 3 (not pictured) always matches the flow velocity at z=3 almost exactly, with no further change occurring in the 𝐳̂ direction.



 §.§ Results


  §.§.§ Role of the central actin bundle

The CAB has long been hypothesized to aid in the removal of vesicles from the apical region, but direct observation of vesicle transport along actin is often hindered by the small size of vesicles (typically below the resolution limit of conventional confocal microscopes). Higher resolution imaging methods such as evanescent wave microscopy have been used to observe long-range vesicle movement (presumably a result of actomyosin transport) in the periphery of the tube <cit.>, but the limited penetration depths available in these methods (≤400nm) do not allow for imaging of the central region. STICS analysis suggests a way of determining whether the central bundle participates in vesicle transport. Directed vesicle movement in the periphery of the tube is not seen to exceed speeds of ≈0.5 m.s^-1, smaller than the ≈0.8 m.s^-1 observed in the centre <cit.>. By using the method of regularized ringlets and prescribing fluid velocity only along, 

  *  the PAB, where |𝐮_T| = 0.5 m.s^-1,

  *  the tube boundary with apical growth speed v_g = 0.1 m.s^-1,

the resulting fluid velocity in the centre should provide further insight insight into whether the central bundle participates in vesicle transport. 

The results of this investigation can be seen in Figure <ref>, with both lab (Figure <ref>) and tip (Figure <ref>) frames shown. The shape of the velocity profile in Figure <ref> is in excellent agreement with the STICS analysis of Bove et al.<cit.>, with a wider band of basal flow through the center than apical flow in the periphery. However, despite an arguably exaggerated prescribed peripheral fluid velocity of 0.5 m.s^-1, fluid velocity in the centre does not achieve speeds of 0.8 m.s^-1. This is a strong indication that the CAB must also participate in the transport of vesicles, particularly considering that that our current implementation of the Stokes equations does not account for variations in local fluid viscosity (known to be larger in the presence of filamentous actin networks <cit.>, reducing the fluid velocity induced by any given force). The increased flow speed through the centre of the tube can be easily accounted for by inclusion of additional drag from actomyosin vesicle transport (known to reach speeds of up to 2 m.s^-1<cit.>) along the CAB, with the largest velocities being observed at the very centre as a result of the reduced cytosolic volume in this region. 



  §.§.§ The influence of growth speed

We now turn our attention to the effect of tube growth speed on the cytosolic velocity profile. A further three new velocity profiles are produced, based on three different growth speeds for the tube (0,0.1,  and  0.2   m.s^-1). These speeds (approximately) correspond to that of a static tube, the typical growth rate cited in Bove et al.<cit.>, and the average growth rate measured by Vidali et al.<cit.> for the Lilium longiflorum species. Fluid velocity is prescribed along 

  *  the PAB (|𝐮_T|=.5 m.s^-1), 

  *  the CAB (u_z=-1 m.s^-1), 

  *  and the tube boundary (v_g = 0,0.1,0.2 m.s^-1),

in each case.

Figures <ref>–<ref> show the velocity profiles for each of these three tube growth speeds in the lab frame. Significant differences can be seen between the three profiles, both in the magnitude and direction of the fluid velocity. In particular, at larger growth speeds there is a wider band of cytosolic flow in the positive 𝐳̂ direction in the peripheral region, and the central band of basal flow is both narrower and of a reduced magnitude. This is an expected consequence of mass conservation, since in a tube with a faster growth speed more fluid must flow towards the apical region to fill the increasing space.

Figures <ref>–<ref> show the velocity profiles for each of these three tube growth speeds in the tip frame. Here, the growth speed of the tube has been subtracted from the 𝐳̂ component of the fluid velocity for each corresponding velocity profile in the lab frame. Interestingly, the differences between the velocity profiles are significantly less pronounced in this frame of reference. Differences in the magnitude of the fluid velocity still persist, but the overall shape of the profiles bear a striking similarity. The persistent shape of the velocity profiles seen in Figures <ref>–<ref> could help explain the observed similarities in the distribution of the apical vesicle population and other elements of the cytoplasm across multiple pollen tube species and throughout different phases of oscillatory growth. 



§ CONCLUSIONS AND FUTURE WORK

In this paper, the regularized fundamental solution for the velocity (single layer potential) and stress (double layer potential) due to an axisymmetric ring of smoothed point forces (the `regularized ringlet') was derived, expanding on the work of Cortez and colleagues <cit.>. The velocity solution, written in the form of complete elliptic integrals of the first and second kind, tends to the singular solution of Pozrikidis <cit.> in the limit as the regularization parameter ε tends to 0. 

For the resistance problem on the translating and rotating unit sphere, the method of regularized ringlets was shown to produce small relative errors in drag calculations and small condition numbers for the underlying resistance matrices. Further testing in the case of the translating unit sphere shows that the regularized ringlet performs favourably compared to the traditional method of regularized Stokeslets both in terms of accuracy (Appendix <ref>) and speed (Appendix <ref>) and that it possesses more satisfactory convergence properties than the singular ringlet (Appendix <ref>). 

The applicability of the regularized ringlet to fluid flow problems involving body surface motion in the 𝐫̂ direction (perpendicular to the line of axisymmetry) was established in the example of Purcell's toroidal swimmer. Using the regularized ringlet method, we were able to reproduce Leshansky and Kenneth's <cit.> results for the scaled propulsion velocity of the toroidal swimmer, propelled by surface tank treading against the direction of motion of its outer surface. Our results show a significant improvement over the asymptotic solution found by integrating a centre-line distribution of rotlets in the limit as slenderness decreases (s_0 → 0). The use of regularized ringlets also yields the drag force at all points on the torus surface for the toroidal swimmer, information that is not readily available using the series or rotlet solutions.

In order to elucidate the further benefits of regularization, we studied the case of cytosolic flow in the growing pollen tube. Here the inducing drag force is the result of the directed actomyosin transport of vesicles inside the fluid, making application of singular Stokeslet methods particularly challenging. Since the area over which the force is applied is also relatively large (and somewhat indeterminate), being able to use the regularization parameter ε to control the spreading of the force is also vital to the solution. 

Using the regularized ringlet solution for the steady Stokes equations, we were able to show that it is highly likely that the central actin bundle plays a role in the removal of vesicles from the apical pool. This has long been hypothesised, with vesicle movement through the central region seen to reach speeds of up to 2 m.s^-1<cit.> (similar to those of actomyosin transport), but direct evidence has been lacking. Fluid velocity profiles based on drag induced by actomyosin transport of vesicles along the peripheral bundle alone shows that cytosolic flow in the central region does not reach these speeds, strongly suggesting that actomyosin transport takes place on the central actin bundle. 

During further investigation, we were able to show that the shape of the cytosolic velocity profile relative to the moving tip is largely independent of the growth speed of the tube (under the assumption that the actin profile is able to keep up with the advancing tip). This provides some insight into how the tube is able to consistently orient its internal architecture across multiple pollen tube species and in spite of changes in growth speed during different phases of oscillatory growth. In a future paper, we will continue our work on the pollen tube by using the cytosolic velocity profiles we have produced in tandem with an advection-diffusion-reaction equation for the spatio-temporal distribution of vesicles in the pollen tube. Investigating the conditions necessary to produce the `inverted vesicle cone' almost universally observed in the apex of the angiosperm pollen tube yields new results concerning the appropriate values for different parameters pertaining to tube growth. We believe that the regularized Stokeslet ring method will be of benefit to further studies in the fluid dynamics of axisymmetric growth and beyond.


The authors acknowledge funding from the Engineering and Physical Sciences Research Council (EPSRC) in the form of a Doctoral Training Award (EP/M508202/1 to JT) and grants (EP/N021096/1 to DJS;  EP/M00015X/1 to RJD).

The authors would also like to thank Prof Anja Geitmann and Dr Youssef Chebli (both of the Department of Plant Science, Faculty of Agriculture and Environmental Sciences, McGill University, Sainte-Anne-de-Bellevue, Canada) for their valuable insights on pollen tube physiology and growth. 








43
fxundefined [1]
 ifx#1
fnum [1]
  #1firstoftwo
 secondoftwo
 
fx [1]
  #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0] `
12 `$12
  `&12 `#12 `1̂2 `_12 `%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

 [Lauga and Powers(2009)]lauga2009hydrodynamics
  authorauthorE. Lauga and authorT. R. Powers, titletitleThe hydrodynamics of
  swimming microorganisms, @noop journaljournalReports on Progress in Physics volume72, pages096601 (year2009)NoStop
 [Phan-Thien et al.(1987)Phan-Thien, Tran-Cong, and Ramia]phan1987boundary
  authorauthorN. Phan-Thien, authorT. Tran-Cong, and authorM. Ramia, titletitleA boundary-element analysis
  of flagellar propulsion, @noop journaljournalJournal of Fluid Mechanics volume184,pages533 (year1987)NoStop
 [Ramia et al.(1993)Ramia,
  Tullock, and Phan-Thien]ramia1993role
  authorauthorM. Ramia, authorD. Tullock,  and authorN. Phan-Thien,titletitleThe role of hydrodynamic interaction
  in the locomotion of microorganisms, @noop journaljournalBiophysical journal volume65, pages755 (year1993)NoStop
 [Shum et al.(2010)Shum,
  Gaffney, and Smith]shum2010modelling
  authorauthorH. Shum, authorE. A. Gaffney,  and authorD. J. Smith,titletitleModelling bacterial behaviour close to
  a no-slip plane boundary: the influence of bacterial geometry, @noop
  journaljournalProceedings of the Royal
  Society of London A volume466, pages1725 (year2010)NoStop
 [Chwang and Wu(1975)]chwang1975hydromechanicsP2
  authorauthorA. T. Chwang and authorT. Y.-T. Wu, titletitleHydromechanics of
  low-Reynolds-number flow. Part 2. Singularity method for Stokes flows,@noop journaljournalJournal of Fluid
  Mechanics volume67, pages787
  (year1975)NoStop
 [Chwang and Wu(1974)]chwang1974hydromechanicsP1
  authorauthorA. T. Chwang and authorT. Y.-T. Wu, titletitleHydromechanics of
  low-Reynolds-number flow. Part 1. Rotation of axisymmetric prolate bodies,@noop journaljournalJournal of Fluid
  Mechanics volume63, pages607
  (year1974)NoStop
 [Chwang(1975)]chwang1975hydromechanicsP3
  authorauthorA. T. Chwang, titletitleHydromechanics of
  low-Reynolds-number flow. Part 3. Motion of a spheroidal particle in
  quadratic flows, @noop journaljournalJournal of Fluid Mechanics volume72, pages17 (year1975)NoStop
 [Chwang and Wu(1976)]chwang1976hydromechanicsP4
  authorauthorA. T. Chwang and authorT. Y. Wu, titletitleHydromechanics of
  low-Reynolds-number flow. Part 4. Translation of spheroids, @noop journaljournalJournal of Fluid Mechanicsvolume75, pages677 (year1976)NoStop
 [Johnson and Wu(1979)]johnson1979hydromechanicsP5
  authorauthorR. E. Johnson and authorT. Y. Wu, titletitleHydromechanics of
  low-Reynolds-number flow. Part 5. Motion of a slender torus, @noop
  journaljournalJournal of Fluid Mechanicsvolume95, pages263 (year1979)NoStop
 [Huang and Chwang(1986)]huang1986hydromechanicsP6
  authorauthorL. H. Huang and authorA. T. Chwang, titletitleHydromechanics of
  low-Reynolds-number flow. Part 6. Rotation of oblate bodies, @noop
  journaljournalJournal of Engineering
  Mathematics volume20, pages307
  (year1986)NoStop
 [Cortez(2001)]cortez2001method
  authorauthorR. Cortez, titletitleThe method of regularized
  Stokeslets, @noop journaljournalSIAM Journal on Scientific Computing volume23,pages1204 (year2001)NoStop
 [Cortez et al.(2005)Cortez,
  Fauci, and Medovikov]cortez2005method
  authorauthorR. Cortez, authorL. Fauci, andauthorA. Medovikov, titletitleThe method of regularized Stokeslets in three
  dimensions: analysis, validation, and application to helical swimming,@noop journaljournalPhysics of
  Fluids volume17, pages031504
  (year2005)NoStop
 [Smith(2009)]smith2009boundary
  authorauthorD. J. Smith, titletitleA boundary element
  regularized Stokeslet method applied to cilia-and flagella-driven flow,@noop journaljournalProceedings of the
  Royal Society of London A volume465, pages3605 (year2009)NoStop
 [Cortez(2018)]cortez2018regularized
  authorauthorR. Cortez, titletitleRegularized Stokeslet
  segments, @noop journaljournalJournal of Computational Physics volume375,pages783 (year2018)NoStop
 [Smith(2018)]smith2018nearest
  authorauthorD. J. Smith, titletitleA nearest-neighbour
  discretisation of the regularized stokeslet boundary integral equation,@noop journaljournalJournal of
  Computational Physics volume358, pages88 (year2018)NoStop
 [Blake(1971)]blake1971spherical
  authorauthorJ. R. Blake, titletitleA spherical envelope
  approach to ciliary propulsion, @noop journaljournalJournal of Fluid Mechanics volume46, pages199 (year1971)NoStop
 [Purcell(1977)]purcell1977life
  authorauthorE. M. Purcell, titletitleLife at low Reynolds
  number, @noop journaljournalAmerican Journal of Physics volume45, pages3 (year1977)NoStop
 [Chebli et al.(2013)Chebli,
  Kroeger, and Geitmann]chebli2013transport
  authorauthorY. Chebli, authorJ. Kroeger,  and authorA. Geitmann,titletitleTransport logistics in pollen tubes,@noop journaljournalMolecular Plantvolume6, pages1037 (year2013)NoStop
 [Pozrikidis(1992)]pozrikidis1992boundary
  authorauthorC. Pozrikidis, @noop titleBoundary integral
  and singularity methods for linearized viscous flow (publisherCambridge University Press, year1992)NoStop
 [Leshansky and Kenneth(2008)]leshansky2008surface
  authorauthorA. M. Leshansky and authorO. Kenneth, titletitleSurface tank treading:
  Propulsion of Purcell’s toroidal swimmer, @noop journaljournalPhysics of Fluids volume20, pages063104 (year2008)NoStop
 [Fairweather and Karageorghis(1998)]fairweather1998method
  authorauthorG. Fairweather and authorA. Karageorghis, titletitleThe method of
  fundamental solutions for elliptic boundary value problems, @noop journaljournalAdvances in Computational
  Mathematics volume9, pages69
  (year1998)NoStop
 [Young et al.(2006)Young,
  Jane, Fan, Murugesan, andTsai]young2006method
  authorauthorD. L. Young, authorS. J. Jane,
  authorC. M. Fan, authorK. Murugesan, and authorC. C. Tsai, titletitleThe method of fundamental solutions for 2D and 3D Stokes
  problems, @noop journaljournalJournal of Computational Physics volume211,pages1 (year2006)NoStop
 [Gradshteyn and Ryzhik(2014)]gradshteyn2014table
  authorauthorI. S. Gradshteyn and authorI. M. Ryzhik, @noop titleTable of integrals,
  series, and products (publisherAcademic Press, year2014)NoStop
 [Gallagher et al.(2018)Gallagher, Choudhuri, and Smith]gallagher2018sharp
  authorauthorM. T. Gallagher, authorD. Choudhuri, and authorD. J. Smith, titletitleSharp quadrature error
  bounds for the nearest-neighbor discretization of the regularized stokeslet
  boundary integral equation, @noop journaljournalSIAM Journal on Scientific Computing volume(in press) (year2018)NoStop
 [Lamb(1945)]lamb1945hydrodynamics
  authorauthorH. Lamb, @noop titleHydrodynamics  (publisherCourier Corporation, year1945)NoStop
 [Rubinow and Keller(1961)]rubinow1961transverse
  authorauthorS. I. Rubinow and authorJ. B. Keller, titletitleThe transverse force on a
  spinning sphere moving in a viscous fluid, @noop journaljournalJournal of Fluid Mechanics volume11, pages447 (year1961)NoStop
 [Thaokar et al.(2007)Thaokar, Schiessel, and Kulic]thaokar2007hydrodynamics
  authorauthorR. M. Thaokar, authorH. Schiessel,  and authorI. M. Kulic,titletitleHydrodynamics of a rotating torus,@noop journaljournalEuropean Physical
  Journal B volume60, pages325
  (year2007)NoStop
 [Taylor(1952)]taylor1952analysis
  authorauthorG. I. Taylor, titletitleAnalysis of the swimming
  of long and narrow animals, @noop journaljournalProc. R. Soc. Lond. A volume214,pages158 (year1952)NoStop
 [Kulić et al.(2005)Kulić, Thaokar, and Schiessel]kulic2005twirling
  authorauthorI. Kulić, authorR. Thaokar,  and authorH. Schiessel,titletitleTwirling DNA ringsâ€”Swimming
  nanomotors ready for a kickstart, @noop journaljournalEPL volume72, pages527 (year2005)NoStop
 [jiro (MathWorks staff)()]grabit
  authorauthorjiro
  (MathWorks staff), @noop titlegrabit function
  (File Exchange), howpublished<https://uk.mathworks.com/matlabcentral/fileexchange/7173-grabit>,noteAccessed: 10/12/2018NoStop
 [Bove et al.(2008)Bove,
  Vaillancourt, Kroeger, Hepler, Wiseman, and Geitmann]bove2008magnitude
  authorauthorJ. Bove, authorB. Vaillancourt,
  authorJ. Kroeger, authorP. K. Hepler, authorP. W. Wiseman, and authorA. Geitmann, titletitleMagnitude and direction of vesicle dynamics in growing pollen tubes
  using spatiotemporal image correlation spectroscopy and fluorescence recovery
  after photobleaching, @noop journaljournalPlant Physiology volume147, pages1646 (year2008)NoStop
 [Booy et al.(1992)Booy,
  Krens, and Bino]booy1992analysis
  authorauthorG. Booy, authorF. A. Krens,  and authorR. J. Bino,titletitleAnalysis of pollen-tube growth in
  cultured maize silks, @noop journaljournalSexual Plant Reproduction volume5,pages227 (year1992)NoStop
 [Fayant et al.(2010)Fayant,
  Girlanda, Chebli, Aubin,
  Villemure, and Geitmann]fayant2010finite
  authorauthorP. Fayant, authorO. Girlanda,
  authorY. Chebli, authorC.-É. Aubin, authorI. Villemure, and authorA. Geitmann, titletitleFinite element model of polar growth in pollen tubes, @noop
  journaljournalPlant Cell volume22, pages2579 (year2010)NoStop
 [Dumais et al.(2006)Dumais,
  Shaw, Steele, Long, andRay]dumais2006anisotropic
  authorauthorJ. Dumais, authorS. L. Shaw,
  authorC. R. Steele, authorS. R. Long, and authorP. M. Ray, titletitleAn anisotropic-viscoplastic model of plant cell
  morphogenesis by tip growth, @noop journaljournalInternational Journal of Developmental Biology volume50, pages209 (year2006)NoStop
 [Lovy-Wheeler et al.(2005)Lovy-Wheeler, Wilsen, Baskin, andHepler]lovy2005enhanced
  authorauthorA. Lovy-Wheeler, authorK. L. Wilsen, authorT. I. Baskin,  and authorP. K. Hepler,titletitleEnhanced fixation reveals the apical
  cortical fringe of actin filaments as a consistent feature of the pollen
  tube, @noop journaljournalPlantavolume221, pages95 (year2005)NoStop
 [Fu et al.(2001)Fu,
  Wu, and Yang]fu2001rop
  authorauthorY. Fu, authorG. Wu, andauthorZ. Yang, titletitleRop GTPase–dependent dynamics of tip-localized
  F-actin controls tip growth in pollen tubes, @noop journaljournalJournal of Cell Biology volume152, pages1019 (year2001)NoStop
 [Vidali et al.(2001)Vidali,
  McKenna, and Hepler]vidali2001actin
  authorauthorL. Vidali, authorS. T. McKenna, and authorP. K. Hepler, titletitleActin polymerization is
  essential for pollen tube growth, @noop journaljournalMolecular Biology of the Cell volume12, pages2534 (year2001)NoStop
 [Wang et al.(2006)Wang,
  Teng, Wang, Li, Sheng, Zheng, Šamaj, Baluška, and Lin]wang2006imaging
  authorauthorX. Wang, authorY. Teng, authorQ. Wang, authorX. Li, authorX. Sheng, authorM. Zheng, authorJ. Šamaj, authorF. Baluška, and authorJ. Lin, titletitleImaging of
  dynamic secretory vesicles in living pollen tubes of Picea meyeri using
  evanescent wave microscopy, @noop journaljournalPlant Physiology volume141, pages1591 (year2006)NoStop
 [Wong et al.(2004)Wong,
  Gardel, Reichman, Weeks,
  Valentine, Bausch, and Weitz]wong2004anomalous
  authorauthorI. Y. Wong, authorM. L. Gardel,
  authorD. R. Reichman, authorE. R. Weeks, authorM. T. Valentine, authorA. R. Bausch, and authorD. A. Weitz, titletitleAnomalous diffusion probes microstructure dynamics of entangled
  F-actin networks, @noop journaljournalPhysical Review Letters volume92,pages178101 (year2004)NoStop
 [De Win(1997)]de1997quantitative
  authorauthorA. H. N.De Win, titleQuantitative analysis of
  organelle movements in pollen tubes, @noop Ph.D. thesis,schoolRadboud University (year1997)NoStop
 [Alves and Silvestre(2004)]alves2004density
  authorauthorC. J. S.Alves and authorA. L.Silvestre, titletitleDensity
  results using Stokeslets and a method of fundamental solutions for the Stokes
  equations, @noop journaljournalEngineering Analysis with Boundary Elements volume28, pages1245 (year2004)NoStop
 [Barrero-Gil(2012)]barrero2012method
  authorauthorA. Barrero-Gil, titletitleThe method of
  fundamental solutions without fictitious boundary for solving Stokes
  problems, @noop journaljournalComputers & Fluids volume62, pages86 (year2012)NoStop
 [Chen et al.(2005)Chen,
  Young, Tsai, and Murugesan]chen2005method
  authorauthorC. W. Chen, authorD. L. Young,
  authorC. C. Tsai, andauthorK. Murugesan, titletitleThe method of fundamental solutions for inverse
  2D Stokes problems, @noop journaljournalComputational Mechanics volume37,pages2 (year2005)NoStop


§ REGULARIZED STOKESLETS IN CYLINDRICAL COORDINATES FOR RINGLET EVALUATION

The expressions S_ij^ε(𝐱_0,𝐱) for fluid point 𝐱_0=(r_0,0,z_0) and ring point 𝐱_n=(r_n,θ,z_n) in cylindrical coordinates are

    S_11^ε   = 1/r̂_ε^3[2(r_0-r_ncosθ)^2 + (r_nsinθ)^2 + (z_0-z_n)^2 + 2ε^2], 
    
    S_12^ε   = 1/r̂_ε^3[-(r_0-r_ncosθ)(r_nsinθ)],
    
    
    S_13^ε   = 1/r̂_ε^3[(r_0-r_ncosθ)(z_0-z_n)],
    
    
    S_21^ε   = 1/r̂_ε^3[-(r_0-r_ncosθ)(r_nsinθ)],
    
    
    S_22^ε   = 1/r̂_ε^3[(r_0-r_ncosθ)^2 + 2(r_nsinθ)^2 + (z_0-z_n)^2 + 2ε^2 ],
    
    
    S_23^ε   = 1/r̂_ε^3[-(r_nsinθ)(z_0-z_n)],
    
    
    S_31^ε   = 1/r̂_ε^3[(r_0-r_ncosθ)(z_0-z_n)],
    
    
    S_32^ε   = 1/r̂_ε^3[-(r_nsinθ)(z_0-z_n)],
    
    
    S_33^ε   = 1/r̂_ε^3[(r_0-r_ncosθ)^2 + (r_nsinθ)^2 + 2(z_0-z_n)^2 + 2ε^2],
 in which r̂_ε = ((r_0-r_ncosθ)^2 + (r_nsinθ)^2 + (z_0-z_n)^2 + ε^2)^1/2.

§ EVALUATING THE LIMITING BEHAVIOUR OF R_ΑΒ^Ε

To understand the behaviour of R_αβ^ε as r_n,r_0→ 0, it is easiest to consider R_αβ^ε in the form

    R_rr^ε(𝐱_0,𝐱_n)    = r_n(-r_0r_nI_0+(2τ-(z_0-z_n)^2)I_1 - 3r_0r_nI_2), 
     R_rz^ε(𝐱_0,𝐱_n)    = r_n(z_0-z_n)(r_0I_0-r_nI_1), 
     R_zr^ε(𝐱_0,𝐱_n)    = r_n(z_0-z_n)(-r_nI_0+r_0I_1), 
     R_zz^ε(𝐱_0,𝐱_n)    = r_n((τ + (z_0-z_n)^2 + ε^2)I_0 - 2r_0r_nI_1), 
     R_θθ^ε(𝐱_0,𝐱_n)    = r_n(r_0r_nI_0 + (τ+ε^2)I_1 - 3r_0r_nI_2),

in which

    I_0    = 4k^3/(4r_0r_n)^3/2(1/1-k^2E) , 
    
    I_1    = 4k^3/(4r_0r_n)^3/2(2-k^2/k^2(1-k^2)E - 2/k^2F), 
    
    I_2    = 4k^3/(4r_0r_n)^3/2(k^4-8k^2+8/k^4(1-k^2)E-4(2-k^2)/k^4F),

with τ = r_0^2 + r_n^2 + (z_0-z_n)^2 + ε^2 and k^2 := 4r_0r_n/(τ+2r_0r_n). Evaluating lim_k→ 0I_n for each n∈{0,1,2} and substituting these back into Equations (<ref>)–(<ref>) yields the desired results for lim_r_n→ 0R_αβ^ε and lim_r_0→ 0R_αβ^ε.


The first step in evaluating these limits is to observe that from the definition of k, it follows that k→ 0 as either r_0→ 0 or r_n→ 0 (or both). Further, noting that

    lim_k→ 04k^3/(4r_0r_n)^3/2 = lim_k→ 04(4r_0r_n/τ+2r_0r_n)^3/2/(4r_0r_n)^3/2= 4/τ^3/2,
 in which the value of τ as k→ 0 depends on whether r_0→ 0 or r_n→ 0 (or both), it is observed that I_0,I_1,I_2 contain common finite term outside of the larger brackets. Evaluating the remaining parts of I_0,I_1,I_2 in the limit as k→ 0 requires employing the power series expansions of the complete elliptic integrals <cit.>, such that 

    F(k)    = π/2(1 + 1/4k^2 + 9/64k^4 + …), 
    
    E(k)    = π/2(1 - 1/4k^2 - 9/64k^4 + …).


Letting E ∼π/2 it follows that

    lim_k→ 0(1/1-k^2E) = π/2.
 Similarly, letting E ∼π/2(1-1/4k^2) and F ∼π/2(1+1/4k^2) it can be found that 

    lim_k→ 0(2-k^2/k^2(1-k^2)E - 2/k^2F) = 0.

Finally, using E ∼π/2(1-1/4k^2 - 9/64k^4) and F ∼π/2(1+1/4k^2 + 9/64k^4) yields

    lim_k→ 0(k^4-8k^2+8/k^4(1-k^2)E-4(2-k^2)/k^4F) = -π/8.

Compiling all of the above gives

    [  I_0→ 2π/τ^3/2,         I_1→ 0, I_2→ -π/2τ^3/2,       as  k→ 0, ]

which upon substitution into the expressions for R_αβ^ε yields 

    lim_r_n→ 0R_αβ^ε   ≡ 0, 
    lim_r_0→ 0 R_rr^ε    = lim_r_0→ 0 R_rz^ε = lim_r_0→ 0 R_θθ^ε  = 0, 
    lim_r_0→ 0 R_zr^ε    = -2π r_n^2(z_0-z_n)/τ^3/2,
    lim_r_0→ 0 R_zz^ε    = 2π r_n(τ+(z_0-z_n)^2+ε^2)/τ^3/2.


§ FLUID FLOW INDUCED BY UNIT RINGLET FORCES

Figure <ref> shows the fluid velocity induced by a ringlet located at (r,z)=(0.5,0.5) with associated unit force in 𝐫̂ and 𝐳̂ directions, both under free–space conditions ((a) & (c)) and in the presence of a cylindrical wall at r=1 ((b) & (d), with the wall also represented by Stokes ringlets). In the presence of the bounding wall, closed streamlines and toroidal eddies are observed.



§ EVALUATING THE DOUBLE LAYER POTENTIAL

Recall the form of the double layer potential

    (DLP)_i = 1/8π∫_∂ D u_j(𝐱)T_ijk^ε(𝐱_0,𝐱)n_k(𝐱)   dS(𝐱),
 
in which the stress tensor T_ijk^ε is given by

    T_ijk^ε(𝐱_0,𝐱) = -6(x_0,i-x_i)(x_0,j-x_j)(x_0,k-x_k)/(|𝐱_0-𝐱|^2 + ε^2)^5/2 - 3ε^2 (x_0,i-x_i)δ_jk + (x_0,j-x_j)δ_ik + (x_0,k-x_k)δ_ij/(|𝐱_0-𝐱|^2 + ε^2)^5/2.


As with the Stokeslet S_ij^ε in the SLP, the stress tensor T_ijk^ε can be expressed in cylindrical coordinates via

    (x_0,1-x_1)     = r_0-r_ncosθ, 
     (x_0,2-x_2)    = -r_nsinθ, 
     (x_0,3-x_3)    = z_0-z_n, 
     |𝐱_0-𝐱|^2    = (r_0-r_ncosθ)^2+(r_nsinθ)^2 + (z_0-z_n)^2.
 
and the Cartesian and cylindrical forms of the flow vector 𝐮 are related via u_j = Θ_jβu_β in which

    Θ(θ) = ([  cosθ -sinθ     0;  sinθ  cosθ     0;     0     0     1 ]),

Recalling that the azimuthal component of the normal to an axisymmetric body is zero, the transformation from polar to Cartesian coordinates is n_i = Ψ_iαn_α where 

    Ψ(θ) = ([ cosθ    0    0; sinθ    0    0;    0    0    1 ]).
 Let 
    Q_αβγ = r δ_α i∫_0^2πΘ_jβT_ijk^εΨ_k γ  dθ,
 such that 
    (DLP)_α = 1/8π∫_0^L Q_αβγu_βn_γ  ds.

For fixed α, nonzero elements of Q_αβγ are given by

    [                                                                                          Q_α 1 1 = r_n δ_α i∫_0^2π(T^ε_i 1 1cos^2θ + T^ε_i 1 2sinθcosθ + T^ε_i 2 1sinθcosθ + T^ε_i 2 2sin^2θ) dθ,;                                                                                          Q_α 1 3                                               = r_nδ_α i∫_0^2π(T^ε_i 1 3cosθ + T^ε_i 2 3sinθ)dθ,;                                                                                          Q_α 2 1 = r_nδ_α i∫_0^2π(- T^ε_i 1 1sinθcosθ - T^ε_i 1 2sin^2θ + T^ε_i 2 1cos^2θ + T^ε_i 2 2sinθcosθ)dθ,;                                                                                          Q_α 2 3                                              = r_nδ_α i∫_0^2π(-T^ε_i 1 3sinθ + T^ε_i 2 3cosθ)dθ,;                                                                                          Q_α 3 1                                               = r_nδ_α i∫_0^2π(T^ε_i 3 1cosθ + T^ε_i 3 2sinθ)dθ,;                                                                                          Q_α 3 3                                                                    = r_nδ_α i∫_0^2π(T^ε_i 33)dθ, ]

such that each Q_αβγ is a linear sum of terms of the form

    ∙_ijk:=r_n∫_0^2πT_ijk^ε∙  dθ.
 
Letting

    J_m,n := r_n∫_0^2πsin^mθcos^nθ/(τ-2r_0r_ncosθ)^5/2 dθ,
 in which J_m,n = 0 for m odd, each of the necessary ∙_ijk can be written as: 

    [                                                                              cos^2θ_111 = -(6r_0^3+9ε^2r_0)J_0,2 + (18r_0^2r_n + 9ε^2r_n)J_0,3 - 18r_0r_n^2J_0,4 + 6r_n^3J_0,5,;                                                                            sinθcosθ_112                          = (3ε^2 r_n + 6r_0^2r_n)J_2,1 - 12r_0r_n^2J_2,2 + 6r_n^3J_2,3,;                                                                              sin^2θ_122                          = -3ε^2r_0J_2,0 + 3ε^2r_nJ_2,1 - 6r_0r_n^2J_4,0 + 6r_n^3J_4,1,;                                                                                cosθ_113                              = -3(z_0-z_n)((2r_0^2+ε^2)J_0,1-4r_0r_nJ_0,2+2r_n^2J_0,3),;                                                                                sinθ_123                                                     = 6r_n(z_0-z_n)(r_0J_2,0-r_nJ_2,1),;                                                                                   1_133                                            = (6(z_0-z_n)^2 + 3ε^2)(-r_0J_0,0+r_nJ_0,1),;                                                                                        ;                                                                            sinθcosθ_211                             = (3ε^2r_n+6r_0^2r_n)J_2,1 - 12r_0r_n^2J_2,2 + 6r_n^3J_2,3,;                                                                              sin^2θ_212                         = -3ε^2r_0J_2,0 + 3ε^2 r_nJ_2,1 - 6r_0r_n^2J_4,0 + 6r_n^3J_4,1,;                                                                              cos^2θ_221                        = -3ε^2 r_0J_0,2 + 3ε^2 r_nJ_0,3 + 6r_n^3J_2,3 - 6r_0r_n^2J_2,2,;                                                                            sinθcosθ_222                                                           = 9ε^2r_nJ_2,1 + 6r_n^3J_4,1,;                                                                                sinθ_213                                                     = 6r_n(z_0-z_n)(r_0J_2,0-r_nJ_2,1),;                                                                                cosθ_223                                           = - 3ε^2(z_0-z_n)J_0,1 -6r_n^2(z_0-z_n)J_2,1,;                                                                                        ;                                                                              cos^2θ_311                              = -3(z_0-z_n)((2r_0^2+ε^2)J_0,2-4r_0r_nJ_0,3+2r_n^2J_0,4),;                                                                            sinθcosθ_312                                           = 6r_0r_n(z_0-z_n)J_2,1-6r_n^2(z_0-z_n)J_2,2,;                                                                              sin^2θ_322                                             = -3ε^2(z_0-z_n)J_2,0-6r_n^2(z_0-z_n)J_4,0,;                                                                                cosθ_313                                               =(6(z_0-z_n)^2+3ε^2)(-r_0J_0,1+r_nJ_0,2),;                                                                                sinθ_323                                                       = (6r_n(z_0-z_n)^2+3r_nε^2)J_2,0,;                                                                                   1_333                                                 = -(6(z_0-z_n)^3 + 9ε^2(z_0-z_n))J_0,0. ]

Using the double angle formulae for sin and cos we are able to express J_m,n purely in terms of even powers of cos, such that

    J_m,n = λ· 2^m ∫_0^π/2(2cos^2θ-1)^n(1-cos^2θ)^m/2cos^mθ/(1-k^2cos^2θ)^5/2 dθ  for  m  even,
 
where λ = r_n(k/√(r_0r_n))^5/8 = 4r_n(√(τ+2r_0r_n))^-5. Note that the upper limit of integration was first reduced by application of the double angle formulae, followed by the even parity of the resulting integrand about π/2. If we further define

    C_m = λ∫_0^π/2cos^2mθ/(1-k^2cos^2θ)^5/2  dθ,
 then expanding Equation (<ref>) for the relevant values of m,n yields


    [ J_0,0 =   +C_0,        ; J_0,1 =    -C_0  +2C_1,        ; J_0,2 =    +C_0   -4C_1  +4C_2,        ; J_0,3 =    -C_0   +6C_1  -12C_2  +8C_3,        ; J_0,4 =    +C_0   -8C_1  +24C_2  -32C_3 +16C_4,        ; J_0,5 =    -C_0  +10C_1  -40C_2  +80C_3  -80C_4 +32C_5,;        ; J_2,0 =  4(+C_1  -C_2),        ; J_2,1 =  4(-C_1   +3C_2 -2C_3),        ; J_2,2 =  4(+C_1   -5C_2   +8C_3 -4C_4),        ; J_2,3 =  4(-C_1   +7C_2  -18C_3  +20C_4 -8C_5),;        ; J_4,0 = 16(+C_2   -2C_3  +C_4),        ; J_4,1 = 16(-C_2   +4C_3   -5C_4 +2C_5), ]


in which the integrals C_m can be expressed in terms of complete elliptic integrals of the first and second kind (F and E respecitvely) with elliptic modulus k as

    C_0     = λ/3(1-k^2)(-F + 2(2-k^2)/1-k^2E), 
    
    C_1     = λ/3k^2(1-k^2)(-F + 1+k^2/1-k^2E), 
    
    C_2     = λ/3k^4(1-k^2)((2-3k^2)F + 2(2k^2-1)/1-k^2E), 
    
    C_3     = λ/3k^6(1-k^2)((8-9k^2)F - 3k^4-13k^2+8/1-k^2E), 
    
    C_4     = λ/3k^8(1-k^2)((16-16k^2-k^4)F - 2k^6 + 4k^4 - 24k^2 + 16/1-k^2E), 
     
    C_5     = λ/15k^10(1-k^2)(128-120k^2-9k^4-4k^6)F - 8k^8+11k^6+27k^4-184k^2+128/1-k^2E).


§ TRANSLATING UNIT SPHERE

Following the example of Cortez et al.<cit.>, consider the translating unit sphere with velocity 𝐮=-𝐳̂. The sphere is parametrized and discretized in the same manner as in Section <ref>. Assuming the sphere experiences zero azimuthal spin, the fluid velocity at any point 𝐱_0 can be approximated using the N ringlets via Equation (<ref>). Using intervals of equal size, the quadrature weight associated with numerical integration over 𝐩 is simply w_n=π/N ∀ n. Now considering the velocity evaluated at the location of each ring in the rz plane yields a system of equations which can be written in matrix form

    𝖴 = 1/8Nμ[ [ 𝖱_rr^ε 𝖱_rz^ε; 𝖱_zr^ε 𝖱_zz^ε ]]𝖦^a,

as outlined in Equations (<ref>) - (<ref>). By setting 
    g^a_r(𝐱_i)=0,    g^a_z(𝐱_i)=3μ/2aU,   ∀ i=1,... ,N,
 where μ,a,U are the fluid viscosity, sphere radius, and sphere speed (in z) respectively, the classical solution for Stokes flow 
    u_r(𝐱_i)=0,    u_z(𝐱_i)=U,   ∀ i=1,... ,N,
 should follow. Note that this calculation is independent of the value of μ by cancellation (besides the implicit requirement that Re=UL/μ≪ 1 for the Stokes equations to be valid), and that we use a=1, U=-1 as described in the outline of the problem. 

The first test involves using N=50 regularized rings so that the grid size is given by π/50≈0.065. This is chosen such that the minimum distance between adjacent rings in the axisymmetric discretization of the sphere surface is approximately the same as the distance between adjacent points in Cortez and colleagues'  discretization using 3N^2 regularized Stokeslets. The regularization parameter ε is varied between 0.005 and 0.1, and the error in the ℓ^2 norm for the z component of the flow field is recorded in each case. This error is defined as 
    |u_z+1|_2 :=√(∑_i=1^N (u_z(𝐱_i)+1)^2/N).
 Division by N is necessary for the sake of comparison of errors with later tests where the value of N will change in order to alter the grid size. Initial results using regularized ringlets are shown in Figure <ref> 
and are favourable compared to those using regularized Stokeslets. The regularized ringlet method appears to be optimal for a lower value of ε than the regularized Stokeslet method, with a minimal error found at ε≈0.015 with our method and ε≈0.025 with that of the regularized Stokeslet, as well as being slightly more accurate for almost all values of ε tested. Interestingly, the magnitude of the errors using the two methods briefly appear to coincide near the point at which the regularized Stokeslet error is minimized.



The second test involves varying the grid size for a fixed value of ε=0.01. Regularized ringlet results are shown in Figure <ref> and again compare very favourably to results using regularized Stokeslets. For larger grid sizes our errors are significantly reduced compared to those found using regularized Stokeslets. This is at least in part a result of the ringlet method being better suited to handling the small value of ε=0.01. With both methods, the error eventually stops decreasing as the grid size tends towards 0 since in this regime the regularization error dominates.



The final test again looks at the effect of varying the value of ε on the magnitude of the numerical error, this time using N=124 ringlets for a grid size approximately equal to 0.026 (note that we do not use N=125 ringlets to avoid placing a ring at the point (r,z)=(1,0), which would result in a singular velocity when ε=0). The velocity error is compared at two distinct points: (r,z)=(1,0) lying on the surface of the sphere, and (r,z)=(1.5,0) lying a distance of half the sphere radius away. Ringlet results are shown in Figures <ref>&<ref>, and once more match regularized Stokeslet results very closely. The magnitude of the error is linear with respect to ε on the surface of the sphere and quadratic a sufficient distance away.



§ COMPARISON TO SINGULAR SOLUTIONS

Regularized solutions for Stokes flow have the advantage of being simple to implement and readily usable without needing to worry about the presence of singularities in the computational domain. This does not mean that singular solutions cannot be used, however; so long as the appropriate care is taken to deal with the singularities in some way, singular solutions can also yield excellent results. The method of fundamental solutions <cit.> is a popular choice for implementing singular Stokeslet solutions, in which a fictitious boundary is placed outside of the computational domain and adjacent to the physical boundary of the problem considered. Stokeslets (or `source points') are placed on this fictitious boundary and are associated with collocation points (typically of an equal number) on the physical boundary, with the force density for each Stokeslet being calculated using the resistance matrix such that the physical boundary conditions are satisfied. What the appropriate distance between the fictitious and physical boundaries should be is difficult to determine a priori, and in some sense this distance can be regarded as a regularization parameter for the singular problem <cit.>. If the separation distance is too small, the proximity between the Stokeslet singularities and the physical boundary may lead to inaccurate solutions, whereas if the distance is too large the resistance matrix may become ill-conditioned <cit.>. In some cases, placement of the fictitious boundary may also be constrained by the geometry of the problem itself, leading to solutions that are far from optimal. The difficulty in balancing all of these factors is one of the reasons for the popularity of regularized methods. 



In the case of the axisymmetric ring of singular Stokeslets, some interesting behaviour occurs in the limit as the source and collocation points coincide. In this limit, k tends to unity from which it follows that F →∞ and E→ 1 in Equations (<ref>)–(<ref>). By employing the asymptotic expansion F ≈ -lnr̂ + … in which r̂ = |𝐱_0 - 𝐱_n|, we find that R_θθ^0 ≈ 2R_rr^0≈ 2R_zz^0 ≈ - 4lnr̂ + … which all tend to infinity as r̂→ 0 but at a significantly slower rate than the individual Stokeslet (≈r̂^-1). The matrix elements R_rz^0 and R_zr^0 are similarly divergent but typically take values in the range [-1,1]. Figure <ref> shows a comparision of the magnitude of the singularities as r̂→ 0 for R_rr^0,R_zz^0 ≈ -2lnr̂ (red circles) versus S_ij^0≈2r̂^-1 for i=j (blue squares). A log-log plot must be employed due to the speed with which the Stokeslet singularity increases for small r̂. 



As a result of the slow rate with which the singularity in the axisymmetric Stokeslet ring approaches infinity, the method of fundamental solutions can be employed with a separation distance between the fictious and physical boundaries that is several orders of magnitude smaller than is possible using individual Stokeslets. This allows for the discovery of accurate solutions with a well-conditioned underlying resistance matrix. We illustrate this by once again conducting Matlab simulations for the resistance problem on the translating unit sphere (velocity -𝐳̂) using both regularized and singular Stokeslet rings.

As before, the sphere is parametrized in 2D using the arc segment 𝐩 = cosφ 𝐫̂ + sinφ 𝐳̂ for φ∈ [-π/2,π/2] with N=400 points placed along its length at locations 𝐱_n. For the method regularized ringlets, these 𝐱_n denote ringlet locations. The velocity 𝐮(𝐱_n)=-𝐳̂ is prescribed at all ringlet locations and the resistance matrix is constructed to yield the required force densities 𝐠^ℓ(𝐱_n). The condition number of the resistance matrix 𝖱^ε is determined using the "cond" function in Matlab. The analytical solution for the drag on the unit sphere with unit velocity -𝐳̂ is known to be equal to (3/2)𝐳̂, so that the mean relative error for the calculation of the force term resulting from the resistance matrix is given by

    e(g^a_z) = 1/N∑_n=1^N | g^a_z(𝐱_n) + 3/2/-3/2|,
 in which we recall that 𝐠^a is the force per unit area exerted by the sphere on the fluid and so has opposite sign to the drag. The fluid velocity at the near– and far–field locations (r,z) = (0,1.1), (0,1.5) is evaluated using the force densities 𝐠^a according to Equation (<ref>). The relative error in these velocities is given by

    e(1.1) = |u_z(0,1.1) + 1315/1331/-1315/1331|,    e(1.5) = |u_z(0,1.5) + 23/27/-23/27|,

based on the analytical solution given by Cortez et al.<cit.>. Table <ref> shows a summary of the results obtained for each of these four metrics for various values of regularisation parameter ε.



For the method of fundamental solutions, let 𝐱_n denote collocation points on the sphere surface. N source points are placed at positions 𝐱_s = (1-10^-b)·𝐱_c where b is a parameter that represents the separation distance between the physical and fictitious boundaries (Figure <ref>). The velocity 𝐮(𝐱_n)=-𝐳̂ is prescribed at all collocation points and the resistance matrix is constructed to yield the required force densities 𝐠^ℓ(𝐱_s) at each source point. The accuracy and applicability of the method is measured using the same four metrics as for the method of regularized ringlets. Table <ref> shows a summary of results for various values of separation distance 10^-b, analogous to the results in the regularized case from Table <ref>.



Similar results are achievable using both methods, although the singular method generally yields more well-conditioned resistance matrices and can produce relative errors of a smaller magnitude for the drag calculation and fluid velocities. It should however be noted that it is not possible to minimize each error in the singular calculation concurrently (as the smallest values for e(1.1) and e(1.5) are generally associated with larger values for the condition number and e(g_z^a)). Despite axisymmetry enabling a drastic reduction of the separation distance for the singular problem, the question of what distance is considered `optimal' still persists. This is not an issue for the regularized ringlet, in which case each error achieves a minimal value for similar values of ε (≈ 0.0025 to 0.005). For excessively small ε (< 0.05), the error becomes non–monotonic as is often the case with regularized Stokeslet methods (eg. see Figure <ref>, Cortez et al.<cit.>, Gallagher et al.<cit.>). The same is also true for small b in the singular case. For both regularized and singular methods, reduction of regularization parameter ε or b always appears to result in a reduction of the condition number of the resistance matrix.

Although the singular method can be tuned to give smaller relative errors in either the fluid velocity or total drag separately, it cannot do so simultaneously; regularized ringlets display more satisfactory convergence properties and are the more effective method to minimize errors in both fluid velocity and total drag. 



§ COMPUTATIONAL SPEED

Table <ref> shows a comparison of the computational times t_R and t_S (measured in seconds),  associated with constructing the 2N × 2N ringlet matrix R^ε (in the irrotational fluid case) and the 3N × 3N Stokeslet matrix S^ε respectively. The increase in computational time for computing R^ε is a result of needing to compute the complete elliptic integrals F(k) and E(k) for all combinations of ring locations (𝐱_m,𝐱_n)  ∀  m,n ∈1,...,N, which requires the construction of two further 2N × 2N matrices E and F. The computational time needed for this isolated operation, t_e is also listed in Table <ref>. The total additional time needed to construct R^ε is modest, typically between 10% and 20% of the time needed for S^ε. 

The computational time associated with F and E (and by extension, R^ε) can be reduced by evaluating F(k) and E(k) to a lower degree of accuracy; the Matlab function "ellipke(k,TOL)" calculates F(k) and E(k) to the accuracy defined by TOL, which has a default value of 2^-52≈ 2.2 · 10^-16 (double-precision accuracy). This is a far greater accuracy than we are typically able to achieve using regularized Stokeslet methods, suggesting a larger value of TOL will suffice. 



We note that by consideration of the size of the matrices involved and the typically small value of t_e, we should hypothetically be able to achieve t_R = (4/9)t_S + t_e < t_S for any given value of N. In practice, this is not the case. The elements S_ij^ε for i,j ∈{1,2,3} share a common form that enables them to be encoded in matrix form 𝖲^ε very efficiently. The same is not true of the elements R_αβ^ε for α,β∈{r,z}, hence why t_R > t_S  ∀  N in Table <ref>.

Inclusion of fluid rotation (such that an additional N × N matrix must be constructed for 𝖱_θ^ε incurs an additional cost of small value, approximately equal to (t_R - t_e)/4 for any given value of N. We divide by 4 under the assumption that the cost associated with constructing the additional N × N matrix 𝖱_θ^ε is one quarter of the cost associated with constructing the irrotational 2N × 2N matrix 𝖱^ε. 

The true value of working with ringlets in 2D can be seen by by evaluating the cost of inverting the matrices 𝖱_θ^ε and 𝖱^ε versus 𝖲^ε (as is necessary in the resistance problem for evaluating the force associated with a given boundary velocity). These matrices are of size N× N, 2N × 2N and 3N × 3N respectively. We investigate the cost of solving a linear system 𝐗 = 𝐀∖𝐛 in which 𝐀 and 𝐛 are a M × M matrix and a M × 1 vector of normally distributed random data respectively, with M  ∈ {N,2N,3N}. We denote the time taken to solve this system by t_N,t_2N, and t_3N in each case. The results are given in Table <ref>, in which it is clear to see that t_3N is significantly larger than the sum of t_N and t_2N for all values of N tested. In practice, this means that using a 3D Stokeslet implementation with 𝖲^ε for solving an axisymmetric resistance problem will always be significantly more costly than our 2D ringlet implementation with 𝖱^ε, 𝖱_θ^ε.

As well as producing smaller resistance matrices with a reduced associated computational cost for a given number of nodes, our axisymmetric ringlet method also requires far fewer nodes in order to achieve the same level of accuracy as the traditional regularized Stokeslet method. In Appendix <ref> it was shown that for the case of the translating unit sphere the axisymmetric discretization of the sphere surface with N rings produces results that are consistently more accurate than the traditional 3D patch discretization using 3N^2 nodes. To put this into context, the result using 1000 ringlets in 2D (computational time ≈0.1s) thus corresponds to using 1,500,000 regularized Stokeslets in 3D. A lower bound estimate for the associated computational time (using t_3N∝ 3N with data from Table <ref>) for method of regularized Stokeslets using this number of nodes yields a time in excess of 500 days. Assuming each element of the resulting matrix 𝖲^ε is stored as a double-precision floating-point number, this would also require ≈1PB of storage.



§ FORCE CALCULATION FOR THE ROTATING SPHERE

The leading order expression (in terms of Re) for the torque on the rotating sphere is given by, 
    𝐓=-8πμ a^3Ω,
 the derivation of which can be found in <cit.>. This torque is associated with a drag force per unit area on the surface of the sphere given by 𝐟=-3μω_0(r/a)𝐞̂_θ, which can be verified by considering the identity, 
    𝐓 = ∬_S 𝐱×𝐟  dS,
 where S denotes the sphere surface. Multiplication by 𝐞̂_z yields, 
    -8πμ a^3ω_0 = 𝐞̂_z ·∬_S 𝐱×𝐟  dS.
 In cylindrical coordinates, we have 
    𝐞̂_z·(𝐱×𝐟)     = 𝐞̂_z·((r𝐞̂_r+θ𝐞̂_θ+z𝐞̂_z)×(-3μω_0(r/a)𝐞̂_θ)) 
       = -3μω_0(r/a)𝐞̂_z·(r𝐞̂_z-z𝐞̂_r) 
       = -3μω_0(r^2/a).
 Converting to a spherical system (r,θ,φ) in which θ denotes the azimuthal angle and φ the polar angle, we substitute r^2=a^2sin^2φ and dS = a^2sinφ dθ dφ such that, 
    𝐞̂_z ·∬_S 𝐱×𝐟  dS    = -3μ a^3ω_0 ∫_θ=0^2π(∫_φ=0^πsin^3φ  dφ) dθ, 
       = -6πμ a^3 ω_0 ∫_φ=0^πsin^3φ  dφ, 
       = -6πμ a^3 ω_0 ·4/3 = -8πμ a^3 ω_0,
 as required.