THE ROLE OF THE DOUBLE LAYER POTENTIAL IN REGULARIZED STOKESLET MODELS OF SELF-PROPULSION

: The method of regularized stokeslets is widely-used to model microscale biological propulsion. The method is usually implemented with only the single layer potential, with the double layer potential being neglected, despite this formulation often not being justiﬁed a priori due to non-rigid surface deformation. We describe a meshless approach enabling inclusion of the double layer which is applied to several Stokes ﬂow problems in which neglect of the double layer is not strictly valid: the drag on a spherical droplet with partial slip boundary condition, swimming velocity and rate of working of a force-free spherical squirmer, and trajectory, swimmer-generated ﬂow and rate of working of undulatory swimmers of varying slenderness. The resistance problem is solved accurately with modest discretization on a notebook computer with the inclusion of the double layer ranging from no-slip to free slip limits; neglect of the double layer potential results in up to 24% error, conﬁrming the importance of the double layer in applications such as nanoﬂuidics, in which partial slip may occur. The squirming swimmer problem is also solved for both velocity and rate of working to within a few percent error when the double layer potential is included, but the error in the rate of working is above 250% when the double layer is neglected. The undulating swimmer problem by contrast produces a very similar value of the velocity and rate of working for both slender and non-slender swimmers, whether or not the double layer is included, which may be due to the deformation’s ‘locally rigid body’ nature, providing empirical evidence that its neglect may be reasonable in many problems of interest. Inclusion of the double layer enables us to conﬁrm robustly that slenderness provides major advantages in efﬁcient motility despite minimal qualitative changes to the ﬂow ﬁeld and force distribution.


Introduction
In our contribution to this Special Issue, we discuss the double layer potential in the regularized stokeslet boundary integral equation, focusing on the circumstances in which it can be formally eliminated, before describing a method for its 'meshless' implementation in the manner of the almost ubiquitous single layer regularized stokeslet method. This approach is then applied to three representative resistance and swimming problems in microscale biological fluid mechanics in which the double layer potential may be important. In the Introduction below we give a brief overview of the context of the work, before describing the mathematical background to the method of regularized stokeslets and its computational implementation, also discussing the nearest-neighbor approach that will be used in this manuscript. In Materials and Methods and supporting Appendices we describe our implementation of the double layer potential with a meshless point cloud and its application to the problems considered herein. In Results we discuss convergence and efficiency of the method, and demonstrate the calculation of flow fields, swimming velocity/trajectory and rate of working, focusing on the effect of including or neglecting the double layer term in situations in which this is not formally justified. Finally in the Discussion, findings are summarized and potential future applications and areas of further methodological development are considered.

Literature review
Microscale biological flow is typically dominated by viscous forces relative to inertia; in the absence of non-Newtonian effects, the fluid dynamics may be approximated by the Stokes flow equations, − ∇p + µ∇ 2 u = 0, ∇ · u = 0, where u = u(x, t) is fluid velocity, p = p(x, t) is pressure, and µ is dynamic viscosity. In systems associated with microscale propulsion and pumping, e.g. by cilia and flagella, these equations are typically accompanied by the no-slip, no-penetration boundary condition u(X, t) =Ẋ(t), for points X(t) on the boundary. Due to the Stokes flow equations being linear and having no explicit time dependence, mathematical complexity in microscale biological flow arises from the shape and movement of the boundaries, for example the movement of the cell body through the fluid, beating of cilia and flagella, and peristaltic contractions of bounding walls. Prior to 2001, the majority of modeling studies either used local force-drag approximations [1], slender body theory [2,3], perhaps accompanied by image systems [4], or the standard boundary element method [5][6][7]. The latter approach is distinguished by both accuracy and efficiency, however the learning curve in their application is steep due to both the need to compute singular integrals, and requirement to generate a 'true surface mesh', i.e. a set of surface points which are assigned to the elements of a surface partitioning. Methods based on volumetric discretizations (e.g. finite difference/element/volume methods) have typically been less popular in the field due to the relatively high computational cost associated with constructing and moving a body-fitted volumetric mesh. A significant exception is the immersed boundary method [8][9][10], which exploits a regular discretization and moreover enables the solution of the nonlinear equations arising from finite Reynolds number and viscoelasticity -although much of the application of this method in microscale flow has been for 2D flow problems rather than fully 3D.
Since 2001, the Method of Regularized Stokeslets, proposed by R. Cortez [11] and further developed in three dimensions [12,13] has become a popular approach for these systems due to its ease of implementation, particularly in its original 'Nyström' form, which discretizes the boundary integral equation using a single quadrature rule, then applies collocation at each quadrature point. The method removes the need for a volumetric mesh, and also the requirements to generate a true surface mesh or evaluate weakly-singular integrals, as would be needed by the standard boundary element method. In practical terms, it is necessary only to generate a list of points covering the surfaces in the flow, for example the bodies and flagella of swimming cells, beating cilia and epithelial surfaces. The problem of computing drag on a body undergoing prescribed motion (resistance problem) is then reduced to inverting a linear system to find a vector of unknown stokeslet strengths (forces) at each discretization point. The problem of computing translation and rotation of a body due to movements prescribed only in a frame moving with the body (swimming problem) can be expressed similarly as inverting a linear system for the force augmented with additional unknowns for the velocity and angular velocity, and additional constraints on the total force and moment [14]. Recent examples of the success of the method of regularized stokeslets implemented via a Nyström discretization include: elucidating the role of bacteria flagellar polymorphism in bundling/unbundling and hence run-and-tumble behaviour [15], estimation of hydrodynamic forces triggering protist contraction [16], simulating the dynamics of elongated fibres in shear [17,18], modeling of flows due to sperm-templated microrobots [19], and modelling eukaryotic flagella synchronisation [20].
Other key examples over the 4 of 26 boundary integral representation of any problem in which the boundary motion cannot be decomposed into rigid body motions. Elimination of the double layer is possible in certain situations, for example flows around non-rigid volume-conserving bodies, however a modification to the single layer density is required, which complicates the calculation of the total force and moment (and solution of swimming problems which have zero total force and moment), and may affect the calculation of the rate of working. Notable exceptions include the Spagnolie & Lauga's application of the method of regularized stokeslets with both single and double layer potential to model slip-velocity squirmers interacting with a plane boundary [61], and Montenegro-Johnson, Lauga and colleagues' work on both phoretic [37] and undulating [36] swimmers.
Given that many problems of interest, particularly those involving bending flagella, undulating cell bodies, effective slip flows produced by cilia [62], and nanofluidics applications [63] cannot strictly be represented by combinations of rigid body motions, we will continue this line of research in the present manuscript. We will describe how to implement the double layer potential in the context of the nearest-neighbor regularized stokeslet method, before testing the method for two problems with known analytical solutions, the calculation of the drag on a translating sphere with free slip or partial slip boundary condition, and the calculation of the translational velocity and rate of working of a slip-velocity spherical squirmer in an infinite fluid. Finally we apply the method to undulatory slender and non-slender swimmers (to assess whether slenderness affects the surface deformation in a way which influences the importance of the double layer potential) both in infinite fluid and between plane boundaries, and assess the effect of the neglect of the double layer potential on the calculation of rate of working.
In the remainder of the Introduction, we briefly recapitulate the regularized stokeslet boundary integral equation and the circumstances in which the elimination of the double layer potential is justified, before describing the methodology for meshless discretization of the double layer potential, and some example applications.

Mathematical background of regularized stokeslets and the double layer potential
The linearity of equations (1) enables boundary conditions to be satisfied by taking discrete sums, line or surface integrals of fundamental solutions. The classical (singular) stokeslet or Oseen tensor is the solution to, where δ is the three dimensional Dirac delta function andê k is the unit basis vector pointing in the k-direction. Defining, the pair of tensors (S ij , P j ) then provide the solutions u = (8πµ) −1 (S 1j , S 2j , S 3j ) and p = (8π) −1 P j to equations (3), with σ ij = (8πµ) −1 T ijk the corresponding stress. The regularized stokeslet [11,12] is the exact solution to the Stokes flow equations with spatially-smoothed point force, ∇ · u = 0.
where φ (x) is a family of "blob" functions which approximates δ(x) as → 0 and k = 1, 2, 3. With the most frequently-used choice of the regularized stokeslet pressure, velocity and stress fields are defined by Hence define respectively the pressure, velocity and stress due to the spatially-smoothed point force.
The method of regularized stokeslets is then based on formulating a boundary integral equation formulated in terms of the near-singular fundamental solutions (10)- (13). Following [12], a 'regularized stokeslet version' of the Lorentz reciprocal theorem can be derived for any Stokes flow (u, p) and point y, In equation (14) and throughout, repeated i, j, k indices imply summation over {1, 2, 3}. Equipped with equation (14), we now consider Stokes flow in the unbounded threedimensional space external to a volume D with smooth, orientable surface ∂D. This volume could represent, for example, a sedimenting rigid body, a droplet, or a swimming cell. We will take D to be topologically open, so that its boundary ∂D ⊂ R 3 \ D. Define B R to be a ball of radius R sufficiently large to contain D, and consider the region consisting of the boundary ∂D and exterior volume of fluid which is given by Ω R = B R \ D. Integrating over the volume Ω R , we have for any point y either exterior to D or on its surface, 1 8πµ Denoting the outward pointing unit normal to Ω R by n (inward pointing to D on ∂D), applying the Divergence Theorem, and finally taking the limit as R → ∞ then yields the (exact) boundary integral equation, where f i := −σ ik n k is the traction, i.e. the force per unit area exerted on the body by the fluid. For the standard singular boundary integral equation, the right hand side of equation (16) would involve a Dirac delta function, thereby evaluting precisely to c(y)u j (y), where c(y) = 1 if y is interior to the fluid region R 3 \ D. In the case that y ∈ D then this coefficient is adjusted based on the solid angle exterior pointing into the fluid at y; in the usual case that the surface is smooth then c(y) = 1/2. In the regularized case, the right hand side of equation (16) can only be written approximately in terms of u j (y) and an error which is linear in vicinity of ∂D and quadratic otherwise. Given that calculations typically involve boundary collocation we will use the worst case of p = 1 throughout, yielding the regularised approximate boundary integral equation, valid for y ∈ R 3 \ D, i.e. on the boundary ∂D or exterior to D, whereκ(y) is the mean curvature of ∂D at y. The dependency on curvature is specific to the non-elimiated double layer equation and associated approximation of R 3 \D u j φ dV as opposed to R 3 u j φ dV (see ref. [38], which extends the analysis of Cortez et al. [12]). The left hand sides of equations (16) and (17) consist of two terms: a surface integral of regularized stokeslets S ij multiplying the traction, referred to (by analogy with electric potential theory) as the single layer potential (SLP), and a surface integral of the regularised stokeslet stress T ijk n k multiplying the surface velocity, referred to as the double layer potential (DLP). A similar equation can be derived for situations involving completely or partially bounded Stokes flow (Appendix A).
The 'more nearly-singular' (O( −2 ) as x → y) behaviour of the stress tensor T ijk has motivated most studies using the method of regularized stokeslets to eliminate this term, leaving a single layer equation. There are (at least) two situations in which this elimination can be mathematically justified through identical arguments to singular boundary integral theory -one of which involves an important caveat.
Suppose that D is undergoing rigid body motion, equivalent to the flow throughout the interior and surface D ∪ ∂D being given by u i (x) = U i + ε ijk Ω j x k where U is velocity and Ω is angular velocity. Then the double layer potential can be simplified as (Appendix B), Substituting into equation (16) we then have, for y on the surface of, or exterior to the body, The right hand side can then be approximated [12] by u j (y) + O( ) (note the lack of dependence onκ in this case). By similar arguments, equation (18) can also be approximated, up to regularization error, DLP j (y; u; ∂D) = c(y)u j (y) + O(κ(y) ), for y ∈ ∂D.
This identity, which is an approximate form of a well-established result for the standard boundary integral equation [64], will be useful for solving the swimming problem with the double layer potential. The second situation in which the double layer potential may be eliminated refers to bodies which are not necessarily rigid, but which do not change volume -for example flexible swimming cells or certain models of ciliates involving a surface tangential slip velocity, i.e.
This condition implies the existence of a (unique) solution to the Stokes flow equations in the region interior to D, with velocity u * i (x) and stress tensor σ * ij (x) satisfying u * i (y) = u i (y) for all y ∈ D (a 'fictitious solution'). Applying equation (14) and (16) leads to (Appendix C) the single layer approximate boundary integral equation, where for brevity here and below the regularization error associated with the approximation of the right hand side by u j (y) will be omitted but implied.
The key issue for practical application is that equation (22) is not formulated in terms of the physical traction f , which may present difficulties if the total force and moment must be determined, or are part of the specification of the problem. For example, the resistance problem involves specifying u j (y) for y ∈ ∂D and determining the force and moment, If the adjusted density q but not the physical traction f is known, then neither F nor M can be determined directly, therefore the resistance problem cannot in general be solved from equation (22) accurately if the body is not undergoing a rigid body motion. This issue will be exhibited in section 3.1.
The swimming problem involves specifying the force and moment F and M, along with the surface velocity relative to a body frame moving with an unknown rigid body motion U + Ω ∧ x, and then determining the surface traction f along with U and Ω. Again, the fact that equation (A6) involves only q is problematic. In the particular case F = 0 = M, characteristic of inertialess and neutrally-buoyant swimming, the fictitious flow stress can be shown to result in zero total force and moment (Appendix (D)), and hence the conditions in terms of the adjusted single layer density, are equivalent to the conditions in terms of the physical density f , The zero-force and moment, volume conserving swimming problem can therefore be formulated in terms of the single layer potential only, giving (up to numerical error) equivalent results for the swimming velocity U and angular velocity Ω -although quantities derived from f such as the rate of working ∂D u · f dS may not in general be given accurately by the equivalent calculation from q. This issue will be evident in the model of the squirming swimmer, (Results section 3.2).

Numerical discretization
The final aspect to introduce is the spatial discretization of the force and quadrature as typically implemented for the Nyström or nearest-neighbor discretizations. First consider the resistance problem, i.e. calculation of the total force and moment due to rigid body motion u j (y) = U j + jk Ω k y . From equation (19) we have (subject to regularization error), . The simplicity of equation (27) and absence of the need for a true mesh with local geometry, or indeed singular quadratures, constitute significant practical advantages over the standard boundary integral method.
Nevertheless, this simplicity comes at the cost of the matrix system size depending on the regularization parameter. Recall that equation (19) possesses an error which is linear in the regularization parameter. Moreover it can be shown that the quadrature error of equation (27) is O(h 2 / ) where h is the quadrature spacing [29]. Therefore, reducing with the aim of reducing the regularization error then demands a finer spatial discretization in order to control the quadrature error. The size of the resulting linear system then grows, rapidly escalating the computational cost. As discussed above, strategies to address this issue include improving the order of the regularization error [53,55] and boundary elements/regularized stokeslet segments [30,33].
The approach we will focus on is nearest-neighbor interpolation, which preserves the simple meshless nature of the method of regularized stokeslets, while significantly improving efficiency and accuracy. This method [41,42] utilizes coarse force and fine quadrature discretizations, with nearest-neighbor interpolation being used to project from the fine to coarse discretization. This approach removes the coupling between the system size and regularization parameter and for non-overlapping force/quadrature sets it also removes the dependence of the quadrature error on the regularization parameter. Moreover the method can be implemented in terms of basic linear algebra operations to exploit associated hardware and software optimizations [44].
The basis for the nearest-neighbor method is the use of two point cloud discretizations, one coarse set {x [1], . . . , x[N]} which is sufficient to capture the variation of the traction, and which dictates the size of the linear system, and another finer quadrature discretization set {X [1], . . . , X[Q]} which has sufficient resolution to capture the rapidly-varying kernel S ij (x, y) for y in the vicinity of x. The slowly-varying force per unit area at a quadrature point f (X[q]) is approximated by its value at the nearest point on the coarse force set f (x[n]); defining ν[q,n] = 1 and ν[q, n] = 0 for n =n, and a similar approximation is made for the surface metric dS. Hence, the discrete problem takes the slightly modified form, can be interpreted as the force on each quadrature point neighbouring x[n]. The error analysis of this method is given in ref. [29] and discussed further in ref. [55]; the most important aspects are that with quadrature spacing h q , the quadrature error is O(h 2 q / ) at collocation points which are contained in the quadrature set, and O(h 3 q /δ 2 ) at collocation points which have minimum distance δ from the quadrature set (hence no asymptotic dependence as → 0); the force discretization error and hence linear system size has no intrinsic dependence on . In the following section we focus on the discretization of the double layer potential within this framework.

Materials and Methods
The additional complications to including the double layer potential in the method of regularized stokeslets are: 1.
Evaluation of the higher order near-singularity ( The need to calculate the surface normal n(y) without a true mesh/local geometry. 3.
The fact that the surface metric dS(y) must be calculated explicitly for the implementation of the double layer term, as opposed to being absorbed into the force density, as is the case for the single layer term.
Issue 1. is straightforward to solve, using (the regularized version of) a well-known technique from the boundary element literature. Equation (18) with the specific choice of rigid body motion U j = δ jk , Ω j = 0 and y ∈ ∂D yields (up to regularization error), Recall that if ∂D is smooth in the neighborhood of y then c(y) = 1/2. This identity is a mathematical fact about the double layer potential that can be exploited to construct an approximation for the inner quadrature even when the surface motion is not rigid, as follows: splitting the surface integral into an inner component on ∂D δ (y) = {x ∈ ∂D : |x − y| < δ} sufficiently small that the variation in the velocity is negligible, and an outer component on ∂D \ ∂D δ (y), we may approximate the double layer potential by, Consider the application of the Nyström discretization, with collocation at y = x[m]. Taking δ as the radius of the part of the surface belonging to x[m], the double layer operator evaluated at a collocation point can then be approximated by, where A[n] is the surface metric at x[n] (this term may also be multiplied by a quadrature weight if a higher-order quadrature rule is used). Via the nearest-neighbor discretization, the double layer operator can be approximated by, If evaluating the double layer potential at any point y ∈ ∂D then the simpler approximation may be used. We briefly note that in the MATLAB R implementation accompanying this paper, the identity T ijk (x, y) = −T ijk (y, x) is used, so that matrix rows correspond to collocation points and matrix columns to quadrature points. Issues 2. and 3. require an approximate reconstruction of the local geometry from the quadrature discretization. Below we describe a method to accomplish this using principal component analysis, implemented via the accompanying MATLAB R code that can be applied to any point cloud discretization.
Given . The surface metric and normal can then be approximated using a similar formula to the boundary element method [65], The accuracy of this method in approximating the surface area of a sphere for various values of point spacing h is shown in figure 1, showing approximately linear convergence for each of J = 4, 6, 9.
which approaches the limits for [N] and [F] as λ → 0, ∞ respectively. The resistance problem can be formulated as a second kind Fredholm integral equation by moving the double layer potential in equation (17) across to the right hand side, −SLP j (y; f ; ∂D) = DLP j (y; u; ∂D) + c(y)u j (y).
For problem [N], the velocity u i (y) = δ j1 on ∂D and so we formulate the following problem for the the traction f , for all y ∈ ∂D, the last line following by equation (20). For problem [F], only the normal component of surface velocity is known, i.e. u i (x)n i (x) = δ i1 n i (x) = n 1 (x); the system is completed by the requirement that the unknown traction satisfies the free slip condition f i t 1 i = 0 = f i t 2 i . Hence we have the augmented system (with all unknowns on the left hand side), for all y ∈ ∂D. Equation (39) involves solving for both f and u on the surface; the single and double layer potentials are discretized as in equation (32).
Equation (39) can be generalized to provide the following formulation of problem P: for all y ∈ ∂D. A key step in the numerical implementation is to utilize the metric dS[q] and nearest-neighbor matrix ν[q, n] in order to relate the F j [n] (force) unknowns of the discrete system to the f j (y) (force per unit area) variable of the continuous problem. For the purposes of comparison, we will also evaluate the effect of (without formal justification) eliminating the double layer potential, replacing the double layer term of equation (40) (39), the total drag on the sphere is calculated numerically from the weighted sum, the velocity field at any point y in the fluid exterior to D is calculated from the discrete surface force and velocity solutions via,

Slip-velocity squirmer
The slip-velocity spherical squirmer is one of the canonical models of zero Reynolds number propulsion [62,67,68]; this model provides a coarse-grained representation of microorganisms which produce a surface flow due to the action of many beating cilia. Scaling coordinates so that the radius of the squirmer is 1, we denote the angle relative to the north pole of the squirmer as Θ(y) (so that 0 Θ(y) π) for |y| = 1; we also define t Θ to be the basis vector in the tangent space pointing towards the north pole. Taking a single squirming mode, the velocity boundary condition in a frame of reference moving with the squirmer is, (u j (y)) bf n j (y) = 0, (δ ij − n i (y)n j (y))(u j (y)) bf = sin Θ(y)t Θ i (y), for all |y| = 1, where the superscript 'bf' denotes a frame of reference in which the swimmer is stationary, referred to as the body frame. If the squirmer frame of reference has translational velocity U and angular velocity Ω about the origin, then u j (y) = U j + ε jk Ω k y + (u j (y)) bf .
Because ∂D u j n j dS = 0, the double layer potential can be eliminated in a similar manner to equation (26) to provide a single layer boundary integral equation in terms of the adjusted density q, − SLP j (y; q; ∂D) − U j − ε jk Ω k y = (u j (y)) bf .
The left hand side of the equation is a linear operator on the unknowns q(x), U and Ω.
As discussed in Appendix D, to determine the swimming velocity U and angular velocity about the origin Ω, the force and moment-free conditions can be expressed in terms of q as given in equation (24). Equations (43) - (45) and equation (24) therefore provide a complete model from which the swimming velocity can be determined, although the physical density f and hence derived quantities such as rate of working are not available. Alternatively, equipped with the double layer implementation, we may avoid the elimination of the double layer and work directly with the physical traction f . Again making use of (20), the integral equation is replaced by, − SLP j (y; q; ∂D) − U j + jk Ω k y = (u j (y)) bf 2 + DLP j y; u bf ; ∂D all y ∈ ∂D.
Equations (43), (44), (46) and (25) form a complete model for the squirming swimmer in terms of the physical traction f . Once the traction is computed, the rate of working is given by the integral, An analogous quantity W q can be computed from the adjusted density for the purpose of comparison.
The numerical model can be compared directly to the analytical results of Blake [62]. For a single squirming mode, the exact swimming velocity and rate of working are given by,

Undulating swimmer
The undulating swimmer, resembling a worm or headless eukaryotic flagellum, is described by an origin X 0 (t) (in this case the front tip of the swimmer) and body frame B(t) = (b 1 (t), b 2 (t), b 3 (t)), where b j = (b 1j , b 2j , b 3j ) T are orthogonal unit vectors. Then a point y bf in body frame coordinates has laboratory frame coordinates, Denoting the translational velocity U =Ẋ 0 and angular velocity Ω such thatḃ j = Ω ∧ b j , the kinematics of the swimmer can be expressed by, The prescribed waveform swimming problem involves specifying a time-varying function y bf i (t) at each material point on the surface of the swimmer, hence prescribing the body frame velocityẏ bf i . At any instant in time, given the current position X 0 and orientation B, the aim is to solve for (U, Ω, f ).
To model an undulatory swimmer in an unbounded fluid, the single layer boundary integral equation in terms of the adjusted density q is, Alternatively, working with the full boundary integral model in terms of the physical traction f , we have, Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 23 September 2021 doi:10.20944/preprints202109.0400.v1 The notationẋ bf in the final term is chosen intentionally to denote a variable of integration (as opposed to the free variableẏ bf ). Equation (53) can be adapted to include the effect of stationary boundaries; details are given in Appendix E. The waveform of the centreline of the cell is prescribed in the body frame; because the present study is proof-of-principle rather than a biologically-focused analysis of any specific species, we make use of the activated beat of a sperm flagellum by Dresdner & Katz [69]ỹ (x, t) = (0.1087x + 0.0543) sin(2πx − t).
To construct the centreline of the waveform in the body frame, the waveform (x,ỹ) is translated and rotated so that, the first point of the cell is at y c = 0 and the flagellum is aligned with the y bf 1 -axis. The centreline waveform is parameterised as y c (s, t) where s ∈ [0, 1] is arclength; the normaln(s, t) and binormalb(s, t) are then constructed from the centreline tangent t = ∂ s y c . The surface of the swimmer is then defined by, where the thickness function is defined to be Note that due to the use of the arclength parameterisation, the length of the swimmer is ensured to be constant throughout. The trajectory of the swimmer is formulated as a pair of ordinary differential equations Equations (57) combined with the initial position and orientation are solved with the 4th order Runge-Kutta method with timestep ∆t = 0.02, the inner problem of finding the instantaneous rates of change U and Ω through solution of equations (53) and (25).
All computations were carried out on a mid-specification notebook computer (Lenovo Yoga 9 with 8GB RAM, 11th Gen Intel R Core TM i5-1135G7, 2.40GHz, 4 core) running Windows 10.

Rate of working
The instantaneous rate of working W * (t) is calculated as equation (47) using the physical traction * = f for the full boundary integral model and the adjusted density * = q for the single layer model. This quantity is then averaged over one beat cycle to yield,

Results
In this Section, we present computational experiments with and without the inclusion of the double layer potential for: (1) the resistance problems for a sphere with no-slip, free-slip and partial slip boundary conditions, comparing results against known analytical results; (2) the swimming problem for a squirming sphere, calculating the flow field and comparing the velocity and rate of working against analytical results; (3) undulatory swimming, examining the effect of slenderness on velocity and rate of working, and examining the effect of nearby plane boundaries on the flow field.  Figure  A1.

Resistance problem with no-slip, free-slip and partial slip boundary conditions
Results for the partial slip resistance problem as a function of slip length λ are given in Figure 3, showing very close correspondence to the exact solution over 5 orders of magnitude (red dots compared with solid blue line). These results also confirm that unjustified neglect of the double layer potential for the resistance problem for non-rigid boundary conditions results in a significant error, which grows as the slip length increases (black circles compared with solid blue line) to around 24%.

Slip velocity squirmer 1
The flow field for the slip velocity squirmer is shown in Figure 4a; the source dipole the swimming velocity U and rate of working W with quadrature spacing is shown in 5 Appendix F, Figure A2. With 7200 DoF, a relative error in U below 1% is found with 6 all values of quadrature spacing used, whereas the relative error in W is rather larger, 7 with values of 1-3% being achievable with h q 0.033. Taking these results together, 8 results accurate to within a few percent error for both velocity and rate of working can be 9 computed within around a minute without specialist computing hardware. 10 With the single layer only formulation in terms of the adjusted density q, the velocity 11 was accurate to within 1.2% of the exact value with the most refined discretization (7200 12 DoF and h q = 0.0167). The calculation of rate of working in terms of q yielded a value of 13 W q = 95.8 as compared with the exact value W = 12π = 37.7, a 254% relative error.  15 Results in this section utilize a discretization with 324 DoF, and quadrature spacing 16 h q = 0.00245 for the most slender swimmer and h q = 0.00516 for the least slender; com-17 puted instantaneous velocity and angular velocity were found to be within 3% of results 18 with the force and quadrature discretizations four times finer (Appendix G, Table A1). A and below the swimmer, with vortical behaviour being evident both above and below the 34 swimmer. 35 Finally, we consider the calculation of mean rate of working over a full beat cycle, 36 as shown in Figure 12. From a methodological perspective, comparing the calculation 37 with the physical potentialW f to the single layer only calculation with the adjusted 38 densityW q yields a remarkably small discrepancy of only 1.6-2.0%. Focusing on physical 39 interpretation, as the slenderness ratio varies from 0.02 to 0.06, the mean rate of working 40 varies almost linearly, increasing by 72% from 0.61 to 1.05. Combined with the dramatic 41 reduction in swimming velocity, the slender swimmer is much more efficient than the 42 non-slender swimmer.

44
This manuscript reviewed the derivation of the method of regularized stokeslets, 45 focusing on the formal justification, or otherwise, of eliminating the double layer poten- 46 tial, before describing a method based on principal components analysis to enable its the classic boundary element study of Phan-Thien et al. [6], which showed that for rotating 69 helical flagella attached to a spheroidal cell body, efficiency increased as the flagellum is 70 made thinner relative to the cell body. A similar effect has also been reported in the context 71 of self-motile ribbons, with wider swimmers moving more slowly than narrow [71]. In 72 stark contrast to the squirming swimmer case, neglect of the double layer potential resulted 73 in a rather small discrepancy in mean rate of working, on the order of the numerical error. 74 The latter finding provides evidence that, while not formally justified, neglect of the double 75 layer potential may be acceptable for undulating swimmers such as worms and sperm. We 76 postulate that the character of undulating motion -i.e. the relative rotation of adjacent 77 segments, may mean that the swimmer is behaving in a similar way to a collection of 78 rigid bodies, and hence arguments for the formal neglect of the double layer potential 79 may transfer approximately; the latter may form a topic for future work. This manuscript 80 also did not assess the role of the regularization parameter and regularization error; the 81 effect of reducing , utilizing higher order blobs [53] or Richardson extrapolation [55] in 82 problems involving the double layer potential with the aim of reducing the regularization 83 error may also form a topic for future investigation.

84
In addition to applications in microscale biological flow, the ability to accurately and 85 rapidly model flows and bodies with partial slip boundary conditions may be valuable in 86 the field of nanofluidics [72]. On the scales of tens of nanometers, object surface roughness 87 properties, material hydrophobicity, and a range of other phenomena entail that the no-88 slip condition is not appropriate, and partial slip provides a more realistic description of 89 the boundary-fluid interface behaviour [73]. Generally the simulation tools used in the 90 nanofluidics community are relatively costly particle-based methods such as molecular 91 dynamics [74] and dissipative particle dynamics [75]. While inherently deterministic,

92
Stokes flow is generally considered to give a valid mean description of fluid flow down to 93 the 10 nm scale [76] and thereby provides potential advantages in efficiency.

94
In summary, the method of regularized stokeslets has proved a highly versatile,  115 We discuss below the modification of equation (17) for two situations in which the 116 flow is partially or completely bounded: Flow which is exterior to a volume D but completely enclosed within a bounded The boundary integral equations can then be derived as follows. In case 1 (flow completely bounded), equation (15) is replaced by the same pair of volume integrals over C \ D instead of Ω R . Applying the divergence theorem and again defining n to be the unit normal pointing out of C \ D and traction f i = −σ ik n k we find that, In case 2 (the flow is partially bounded), equation (15) is modified by integrating over Ω R ∩ P, leading to Again denoting f i = −σ ik n k and taking the limit as R → ∞ yields an identical equation to the first case (with C replaced by P), In the commonly-occurring case that the boundary ∂P is stationary, then u j = 0 on the Under the circumstances that D is undergoing rigid body motion with velocity U and angular velocity Ω about the origin so that u = U + Ω ∧ x, the double layer potential can be decomposed as, for all y ∈ R 3 \ D, (A4) the penultimate line using symmetry of T ijk and antisymmetry of ε i k in the indices i and k, 129 and the final line following from the definition of the regularised stokeslet stress T ijk . The velocity and stress field pair (u * i , σ * ij ) satisfying u * = u on ∂D is referred to as a fictitious solution because the interior of the physical object represented by D does not in general consist of fluid. Nevertheless, this mathematical solution satisfies equation (14), which can then be integrated over the volume D to yield, 1 8πµ for all y on the surface of or exterior to D. Applying the divergence theorem, noting u * i = u i on ∂D, and denoting f * i := −σ * ij n j then yields, (A6) The change in signs relative to equation (16) is due to the fact that n points out of D and into R 3 \ D. Adding equations (16) and (A6) gives,