On the limitations of some popular numerical models of flagellated microswimmers: importance of long-range forces and flagellum waveform

For a sperm-cell-like flagellated swimmer in an unbounded domain, several numerical models of different fidelity are considered based on the Stokes flow approximation. The models include a regularized Stokeslet method and a three-dimensional finite-element method, which serve as the benchmark solutions for several approximate models considered. The latter include the resistive force theory versions of Lighthill, and Gray and Hancock, as well as a simplified approximation based on computing the hydrodynamic forces exerted on the head and the flagellum separately. It is shown how none of the simplified models is robust enough with regards to predicting the effect of the swimmer head shape change on the swimmer dynamics. For a range of swimmer motions considered, the resulting solutions for the swimmer force and velocities are analysed and the applicability of the Stokes model for the swimmers in question is probed.

For a sperm-cell-like flagellated swimmer in an unbounded domain, several numerical models of different fidelity are considered based on the Stokes flow approximation. The models include a regularized Stokeslet method and a three-dimensional finite-element method, which serve as the benchmark solutions for several approximate models considered. The latter include the resistive force theory versions of Lighthill, and Gray and Hancock, as well as a simplified approximation based on computing the hydrodynamic forces exerted on the head and the flagellum separately. It is shown how none of the simplified models is robust enough with regards to predicting the effect of the swimmer head shape change on the swimmer dynamics. For a range of swimmer motions considered, the resulting solutions for the swimmer force and velocities are analysed and the applicability of the Stokes model for the swimmers in question is probed.

Introduction
Flagellated microswimmers are cells or micrometre-sized robots that swim by moving appendages called flagella. Bacteria flagella appear as helical filaments rigidly rotated by a motor complex attached to the cell wall, eukaryotic flagella, instead, move by propagating sinusoidal waves in a whip-like fashion. The reason for this difference lies in the specific structure of eukaryotic cilia and flagella: the axoneme [1]. The ability to model flagellated microswimmers mathematically and numerically is relevant to a variety of applications in the fields of & 2019 The Authors. Published by the Royal Society under the terms of the Creative biology, medicine, medical diagnostics and engineering [2]. Beside improving our understanding of the physical phenomenon, accurate models can inform the design of effective microfluidic devices to sort microswimmers by motility [3][4][5] or suggest the design of efficient artificial microswimmers [6].
In this study, we are specifically concerned with the hydrodynamical modelling of sperm-cell swimmers. The dimensionless ratio between the inertia and viscous forces acting on these cells, namely, the Reynolds number, is of the order of Re ¼ UL/n % 10 22 , while the ratio between the characteristic viscous time scale and the time scale representing the rate of deformation of the swimmer body, i.e. the frequency Reynolds number, is of the order of Re v ¼ L 2 v/n % 10 21 . Here L % 50 mm and U % 2 Â 10 À4 m s À1 are the characteristic length and velocity of the swimmer, v % 32 s 21 is the beating frequency of the flagellum and n ¼ 10 26 m 2 s 21 is the kinematic viscosity of water. For Re and Re v ( 1, the flow field is typically computed by integrating the Stokes equations in a time-independent zero-Reynolds number framework (e.g. [7,8]). It can be noted, however, that neither of the two standard Reynolds number definitions include any scale associated with a change of the swimmer shape such as a characteristic wavelength of the flagellum motion that is not uniquely defined by the beating frequency. Thus, in the latter case, the applicability of the common criteria of ignoring the unsteady and inertial effects based on the two Reynolds numbers being of o(1) can be debated. It can be noted that classical studies [9] avoid this controversy by considering simplified flagellated swimmer models, which, for example, cannot capture the hydrodynamically important details of the flagellum waveform near the open ends, and assume that the wavelength of the swimmer's motion is equivalent to its linear size. Although for a simple sperm cell swimming along a straight trajectory the wavelength is more or less equal to the length of the flagellum, in general, for more complex trajectories of the sperm cell that includes sharp turns or other organisms and artificial swimmers this is not the case.
The resistive force theory (RFT) [10], whose foundations were laid in [11], is a further simplified model applicable to the motion of slender bodies of which beating flagella is a good example. The RFT neglects the long-range hydrodynamical interactions and evaluates the viscous forces exerted on the immersed body as a function of the local velocities only. This model presents many advantages: it has a low computational cost when compared with the numerical integration of the Stokes equations allowing for proof of concept calculations [12,13], and it is simple enough to serve as a starting point for further analytical derivations [14,15], yet its accuracy is debated. Early works discussed the best choice for the model parameters, namely the normal and tangential hydrodynamical friction coefficients; different proposals were put forward by Lighthill [16] and Gray & Hancock [17]. Recent experimental tests showed that either choices badly capture the behaviour of helical flagella for the range of shapes present in nature [18,19]. Other studies calibrated the parameters to match experimental observations [20] or the results obtained by integrating Stokes [21]. The necessity of calibrating the model versus experimental observations or the results of more sophisticated models highlights the unsuitability of the RFT approximation for applications outside the range of calibration.
The RFT is applied to the flagellum, the contribution of the approximately ellipsoidal cell-body or 'head' is evaluated through analytical expressions [22] and added to the flagellum contribution to compute the entire cell dynamics. On this same line, one may represent the flagellum by a model of choice and still approach the problem by separately studying the cell body and flagellum dynamics before simply summing their contributions [23]. This procedure is naturally embodied in the RFT and is, for example, implied by studies that look for the optimal flagellum shape, neglecting to include the head in the calculations (e.g. [10,14]).
In this paper, we study the motion of a single sperm cell in an infinite domain to address three issues: (i) the accuracy of approximating swimming as the linear superposition of the dynamics of separate body-parts (head plus flagellum) versus a full-body description, (ii) the accuracy of modelling swimming with the RFT, which is one of the possible approximations based on the above superposition that ignores long-range hydrodynamical interactions, and (iii) the validity of the quasisteady and inertia-less assumption for swimming at the micro-scale, where the quasi-steady hypothesis entails assuming that the flow instantaneously adapts to the body deformations.
To address the first two issues, we contrast the results obtained by applying the simplified approaches and the full hydrodynamical model, § §3.1 and 3.2. In particular, we compare the swimming velocities, trajectory and force distribution on the swimmer body. We then study locomotion of swimmers with different head shapes, §3.3, and show that the simplified approaches fail to identify the most hydrodynamically efficient swimmer. To tackle the third point, we calculate the propulsive matrix, that is the matrix that relates the forces generated by the moving flagellum to the rigid-body velocities of the swimmers, and we analyse its eigenvalues and eigenvectors to identify a criterion that establishes when the inertia-less quasi-steady hypothesis is valid, §3.5.
It should be pointed out that investigation of the applicability limits of classical semi-analytical microswimmer models has been a popular topic since the original seminal papers by Hancock [11,17], who developed slender body theory (SBT) and its algebraic approximation known as the RFT. For example, Higdon [24] developed a version of the finite amplitude SBT to study the importance of taking into account head-flagellum interaction by considering an analytical Green's function solution of the spherical head. The latter was superimposed on the flagellum solution without taking the effect of flagellum on Green's function of the head into account, thus, ignoring some of the non-local hydrodynamic interactions. In the same work, the solutions of RFT with Lighthill's coefficients were compared with the SBT results and large errors (25-50%) were reported in the case of sufficiently large-cell bodies relative to the flagellum length.
Johnson [25] developed a further advanced version of SBT through matched asymptotic expansions with the error term being quadratic in slenderness. This theory was used in the study of Johnson & Brokaw [9] to show that RFT is satisfactory for use in analysis of mechanisms for the control of flagellar bending. Along a similar line of research, Brokaw [26] simulated the behaviour of a spermatozoa flagellum by an active shear system controlled by the curvature of the flagellum. In that work, a good agreement of RFT and SBT was reported.
In accordance with the derivation assumptions [11], SBT models not only require a small thickness of the flagellum relative to its length, but also the amplitude of the displacement being small compared with the wavelength. Within the applicability range, however, SBT models showed excellent agreement with more sophisticated models such as regularized Stokeslet method (RSM) [18], boundary element method [27], finite-element method (FEM) [28], finite volume method solutions and the experimental data [29]. Notably, the test cases considered in these studies typically involve idealized models such as a small rigid helix spiral moving in water at a distance of at least one slender-body length from boundaries, which falls under remit of the SBT approximations. Later studies that focused on the dynamics of such swimmers close to no-slip and/or free-slip boundaries generally adopted a boundary element method approach [30][31][32].
Despite this extensive prior research on the zero-Reynolds number propulsion of flagellated microorganisms, a systematic comparison of RFT, which can be seen as an abridged version of SBT, and a direct solution of the Stokes equation in the case of a realistic, flexible flagellum of a sperm cell [12] with and without fully accounting for non-local hydrodynamic effects of the head/flagellum interaction, has not been performed yet. Neither has a systematic study of the waveform shape on the sperm swimmer's dynamics been conducted with an examination of the common assumptions such as neglecting the flow unsteadiness in accordance with the Stokes' flow assumption. Both of these are a distinct focus and novelty of the present paper.
The paper is organized as follows: in §2, we introduce the numerical methods used to simulate the flagellated swimmer motion. The RSM is presented in §2.1, the swimming problem details in the context of the RSM are discussed in §2.2, while in §2.3 the geometry and beating movement of the flagellum are defined. The application of a finite-element method for the same microswimmer problem is introduced in §2.4. In §3, we report the numerical results including a validation of the RSM versus the FEM code, a modal analysis of the swimming velocities and a visualization of the flow field induced by the swimmer ( §3.4). We summarize the main results in §4.

Regularized Stokeslet method
The fundamental solution of the incompressible forced Stokes equation À rp þ mr 2 u þ fd(x À y) ¼ 0 ( 2 :1) and r Á u ¼ 0, for a point force f acting on y in an unbounded domain is the Stokeslet J (r) royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 180745 where r ¼ jx 2 yj, r ¼ (x 2 y), I is the identity matrix, m is the dynamic viscosity of the fluid and d is the d-function.
Since the Stokes equation is linear, the flow field generated by an immersed body with a deforming boundary S(t) can be represented through a continuous distribution of Stokeslets [33] u(x, t) ¼ ð S(t) f(y) Á J (r) dS y : (2:4) For complex geometries, as is the case of flagellated microswimmers, the integral is computed numerically by discretizing the immersed surface: for n, the points on which the velocity is evaluated, m, the grid points on the surface of the immersed body and A m the quadrature weights where r nm ¼ x n 2 x m is the radius vector from point n to m. In our calculations, the RSM grid is built on the surface of the swimmer head and on the cylindrical surface of the flagellum of radius F r , see §2.3 for details on the swimmer geometry. The Stokeslets are singular kernels, their singularity can be dealt with by replacing the point force with an approximate point force with local support. In practice, as proposed by Cortez et al. [34], fd(r) can be replaced by ff e (r) yielding the regularized Stokeslet for which 97% of the force is within a radius e [34], where e is the regularization parameter that requires calibration. By running some tests for the flow past an ellipsoid (see §2.2) we have found, consistently with [18], that our numerical results minimize the error when e is between one-third and one-half of the grid spacing.

The swimming problem
The velocity of a swimmer in Stokes flow can be decomposed into a rigid-body translation v j , a rigid rotation v j and the body deformation (head and flagellum) u BC j (x). For convenience, we express these velocities in the frame of reference of the swimmer. The head is non-motile (u BC j ¼ 0), while the flagellum moves with the beating motion introduced in §2.3. We consider the case of a flagellum beating on the z ¼ 0 plane.
In the context of the RSM, the swimming problem is solved by inverting the system and F n,i x n,j 1 ijk ¼ 0, (2:10) to find the forces F n,i and the swimming velocities v j and v j . Here F n,i ¼ f n,i A n , for n ¼ 1, 2 . . . N and in equations (2.9) and (2.10) the summation over the repeated index convention is adopted, with i ¼ 1, 2, 3 as well as j and k. The conditions (2.9) and (2.10) derive from the fact that forces and torques need to balance exactly since inertia is absent.
An alternative but equivalent approach consists in computing the propulsive matrix coefficients C x , C y , m z C vx (2:11) and solving the system for v x , v y and v. The coefficients C and m are computed by solving the Stokes equation separately for an arbitrary (unitary for convenience) rigid translation of the swimmer body and an arbitrary solid body rotation. They correspond to the surface integral on the swimmer body of the x, equation (2.12), and y, equation (2.13), component of the force density f i , and the z, equation (2.14), component of the torque density for, respectively, a unitary v x , v y and v. The known terms F B x , F B y , T B z are the integrals of the forces and torque due to the flagellum beating only. This approach requires solving equation (2.8) four times (for an arbitrary v x , v y , v and for the flagellum beating) but has the advantage of producing better conditioned matrices. We recall that a direct consequence of the reciprocal theorem is that the resistance matrix (2.11) is symmetric.
If we approximate the swimming problem by treating the flagellum and the head separately, the propulsive matrix for the frame of reference located on the head centroid becomes where the superscript 't' and 'h' stand for the tail (flagellum) contribution and the head contribution. We will next refer to this approach as the head þ tail (H þ T) model. If the head is spherical C vx,h 3 . Expressions for C vx x C vy y and m v z for a prolate ellipsoid with b ¼ c as the minor semi-axes and a as the major semi-axis are derived in [22], and is the eccentricity, and m v z is calculated for a rotation about a minor axis. To test the accuracy of the RSM, we study the flow past a prolate ellipsoid for N e different stretching ratios, b/a, in the range 0.125 , b/a , 1. A comparison between the analytical result, i.e. expressions (2.16) -(2.18), and the numerical result obtained with the RSM is shown in figure 1. After calibration, the regularization parameter e is chosen to be e ¼ 0: , where S is the surface of the ellipsoid and N is the number of grid points the ellipsoid is represented by. The numerical error, plotted in the inset of figure 1, is computed as where C stands for the coefficient C vx x , C vy y or m v z for a prolate ellipsoid divided by the same coefficient for a sphere of equal volume. Note that for a given volume, the surface area increases as the aspect ratio decreases.
We have compared the results for two different distributions of the grid points: the case of equally spaced points on the surface of the ellipsoid and an the case of an uneven grid (a discretization based on spherical coordinates with points located at equal azimuthal and polar angle intervals). We finally 5 royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 180745 chose to adopt the latter since the results are not very sensitive to the type of grid once the resolution is large enough.

Swimming parameters
We denote the mean flagellar curvature as K 0 , the flagellum frequency as n, l the wavelength and A 0 the amplitude of the wave [20]. Following [21], we first use the values: which give a similar flagellum waveshape compared with [12]. Additionally, we consider the flagellum radius to be F r ¼ 0.25 mm, and the spherical swimmer head to have diameter of L head ¼ 5 mm (figure 2). We then study the case K 0 ¼ 0 which corresponds to a swimmer following a rectilinear rather than circular trajectory. Furthermore, to investigate the effect of the characteristic wavelength on the swimmer dynamics we perform some simulations with different l: l 0 /4, l 0 /2, 2l 0 , l ! 1 (figure 2b).
Following [12,20], for an angle C measured along the flagellum arc length s equal to the flagellum coordinates are These values, together with the null-velocity distribution on the swimmer's head, form the known term of the system of equations (2.8).
As a remark, we stress that the swimming problem is typically solved in the frame of reference of the swimmer (e.g. in [12]) with the origin of the axes on the centre of the head as shown in figure 2. Although arbitrary, this choice is convenient since this is the natural frame to express the velocities of the flagellum. Note, however, that the centre of mass of the head is not the centre of mass of the entire body, the latter moves as the flagellum itself changes shapes and mostly falls outside the swimmer body. A different choice of the frame of reference leads to different values of forces, torques and velocities, but once the results are recasted in a common frame, e.g. the laboratory frame, velocities and trajectories coincide. This is a consequence of the fact that the torque for systems with null-force resultants is independent of the location of the frame of reference. In this case, the torque resultant is zero too.
When l ! 1 the expressions for the coordinates simplify and the integrals in du can be easily computed where cos(sk) and sin(sk) can be expanded in a Taylor series to yield series of, respectively, even or odd powers of sk. In the special case of the rectilinear swimmer, i.e. K 0 ¼ 0, the expansion simplifies into even, for u BC x , and odd, for u BC y , power series of 2A 0 s cos(nt) that correspond, in the frequency domain, to spectra with only even or odd modes different from zero.
After drawing the flagellum centreline, we use the local Frenet -Serret frame to build the cylindrical surface of radius F r , which we discretize by approximately evenly spaced points. The distance between the points is chosen in such a way that the corresponding surface area approximately equals the surface area A relative to the points on the head (see end of §2.2). In conclusion, the regularization parameter for the entire swimmer surface is e ¼ 0:5 ffiffiffiffi A p .

Finite-element method
To cross-verify solutions obtained with the RSM, the system of the governing three-dimensional Stokes equations is solved numerically in the reference frame fixed with the centre of the spherical head of the swimmer. The time period of the flagellum motion is discretized into 100 uniform time steps, which amount was found sufficient for accuracy. Each time moment corresponds to a particular configuration of the flagellum wave shape (2.20) and (2.21). For each shape of the swimmer, the same open-boundary computational box domain around the swimmer is specified. The box size is large enough to simplify the specification of numerical boundary conditions at the external boundaries. The domain is discretized by tetrahedral grid elements with applying a sufficient refinement near the swimmer boundary to resolve both the head and the flagellum surface. The total number of grid cells in the model is about 440 000 and the grid details are shown in figure 3. The grid is generated by using the software magentablack'gmsh'. For each waveform configuration, a converged flow solution is obtained with applying non-slip condition on the swimmer surface with the velocity of the fluid equal to the velocity of the swimmer boundary, which consists of the rigid head and the flexible flagellum parts, and the full slip condition at all external boundaries. The solution obtained is found to be virtually insensitive to any further increase of the computational domain size or a further grid refinement. By integrating the forces on the flagellum surface, the drag force components and the torque specified on the right-hand side of equations (2.12) -(2.14) is calculated. In a similar way, the coefficients for the left-hand side of the same equations, which correspond to the two elementary rectilinear motions in the plane of the swimmer and the elementary rotation of the swimmer around its head centre as of a rigid body, are computed. These amount to four boundary value problems, which correspond to the same governing equations (2.24) and (2.25), the same computational domain, but different boundary conditions. These four problems are solved numerically with a FEM for each time moment in accordance with a particular phase of the swimming cycle ( figure 4). Details of the FEM methods for numerical solution are summarized below. Following the standard approach [35], the FEM with second-order base functions is implemented in the framework of the penalty method, which requires minimization of the following functional: with penalty parameter l , where e xx , e yy , e zz , e xy , e xz , e yz are components of the strain rate tensor (e xx , e yy , e zz , e xy , e xz , e yz ) ¼ @u @x , @v @y , @w @z , 1 2 @u @y þ @v @x , 1 2 @u @z þ @w @x , 1 2 @v @z þ @w @y ; D ¼ @u @x þ @v @y þ @w @z ; and f x , f y , f z are internal forces. This results in a sparse system of linear algebraic equations that is solved using a direct method based on lower -upper (LU) decomposition. The magentablackIntel Math Kernel Library solver is used for solution of the linear system of equations.

Numerical results and data analysis
In figure 5, we plot the x-and y-components of the velocity and the z-component of the angular velocity in the frame of reference of the swimmer for different numerical methods and head shapes. Different points on the swimmer body draw different trajectories on the x-y plane; in figure 6 we display the trajectories of the head centroid. An inspection of the frequency spectra of v x , v y and v reveals that the signal can be reconstructed within a 0.6% error by retaining the first five terms of the Fourier series expansion: P 4 k¼0 a k cos (knt þ f k ). The error is computed as: maxjv 2 v reconstructed j/maxjvj, where v reconstructed is the reconstructed signal from the truncated Fourier series and v the original signal. When the parameter K 0 is set to zero (the rectilinear swimmer) the curves v y and v have zero mean and are described within a 0.6% error by the first and third mode; keeping only the mode k ¼ 1 guarantees an 8% error on v y and a 6% error on v. Differently, v x is described within a 0.3% error by an expansion in the even modes k ¼ 0, 2, 4. We have additionally verified for the rectilinear swimmer that the temporal variation of the angular frequency of the swimmer, v, obtained numerically can be reasonably well approximated (within 6%) by a single harmonic function of the beating frequency, n, regardless of the numerical discretization applied (e.g. 100, 200 and 400 points per the flagellum length).
We stress that since the system is linear no mechanism is in place to allow for the creation of non-zero modes from the forcing/boundary condition, in fact, the non-zero modes detected in v x , v y and v reflect the complex spectrum of u BC x and u BC y (equation (2.22)), which produces a coupling of the x-and y-coordinates of the local reference system of the flagellum. The fact that for the rectilinear swimmer only the even modes are excited in v x and the odd ones in v y matches the u BC x and u BC y spectra for the special case K 0 ¼ 0 as discussed at the end of §2.3. Similar considerations hold when comparing the spectra of F B x , F B y T B z , that is the known term of the resistive matrix system and the unknowns v x , v y , v.   Figure 5. v x , v y components of the velocity in the swimmer frame of reference (a,b) and angular velocity (c) for a cell with a spherical head (solid-blue, circle magenta and dashed-green) or an ellipsoidal head with aspect ratio 0.5 (dotted-red) and 0.25 (dashed-dot black). The circle-magenta points have been computed by an FEM code for comparison, the dashed-green curve were computed by using the simplified approach discussed in §3.2. The frame of reference for these calculations is located on the head centroid.  We have verified that for the calculations presented in this paper the motility matrix coefficients are symmetric within numerical precision as dictated by the reciprocal theorem.
As a further remark, note that the curvature parameter K 0 . 0 of the swimmer waveform in (2.20) corresponds to a circular trajectory in the absolute frame of reference as shown in figure 6. Accordingly, this should give rise to apparent accelerations in the swimmer's frame, which are not accounted for in the Stokes model. To justify the neglect of these accelerations, we want to evaluate the order of magnitude of these terms first. The difference between the accelerations in the non-inertial frame a r and those in the inertial frame a f are where the first term on the right-hand side represents the centrifugal acceleration, the second the Coriolis acceleration and the third the Euler acceleration. The force associated to a r in our calculations is at most of the order O(1 Â 10 À15 ), that is three orders of magnitude smaller than the smallest coefficients in the motility matrix, thus negligible as initially hypothesized. Still, it can be argued that even a small unbalanced force can build up into a non-negligible effect for the swimmer trajectory over a time period long enough compared with the swimmer cycle. Therefore, to confirm that the effect of the non-inertial forces on the trajectory of the swimmer is small, we compared the swimmer's trajectories with and without taking the apparent accelerations into account in accordance with the 'instantaneous' coordinate and velocity of the swimmer calculated numerically. Over a few circular trajectory periods, the swimmer trajectories with and without taking the apparent accelerations into account virtually coincided, which finally justifies the neglect of these terms.

Resistive force theory
According to the RFT, the viscous forces applied to the flagellum centreline depend on the flagellum velocities through where K T and K N are the tangential and normal friction coefficients, t ¼ (cosC, sinC) the tangent to the flagellum centreline and I the identity matrix. Equation (3.2) can alternatively be expressed as [21] f ¼ RKR À1 u BC (3:3) Table 2. Values of the time averaged and rescaled root mean square error for the swimming velocities computed with the RFT-L, the RFT-GH, the head þ tail (H þ T) and the FEM model. The results are for the rectilinear (K 0 ¼ 0 rad s 21 ) and curved (K 0 ¼ 7735.5 rad s 21 ) swimmer with l ¼ l 0 . and K is a diagonal matrix with K T and K N on the diagonal. We proceed by (i) substituting the expressions for C and u BC provided in §2. Following [18] we contrast the results for two possible choices of the friction coefficients: those derived by Lighthill (RFT-L) [16] and those suggested by Gray & Hancock (RFT-GH) [17] K T,L ¼ 2pm ln (0:18l=F r ) , K N,L ¼ 4pm ln (0:18l=F r ) þ 1=2 (3:5) and : We also compare the results obtained with the RFT with those obtained with the RSM. In figure 7 and table 2, we quantify the error of the former as where 'h i' is the average over one beating period, u Stk refers to the value computed with the regularized Stokeslet model, and the denominator rescales the root mean square error by the range of variability of the quantity under consideration, being it v x , v y or v. For the reference case of the swimmer with l ¼ l 0 (solid blue line in figure 5), we find that Lighthill coefficients outperform Gray and Hancock's both for the curved and rectilinear swimmer ( figure 7 and table 2). However, the model performances depend on the geometry, as clearly seen in figure 7. For l , l 0 , we find that the RFT-GH model gives better answers than RFT-L and for l ¼ l 0 /4 both RFT models are deemed unreliable. We have also checked that the simplified model referred to as the head þ tail model (H þ T), which will be discussed in the next section, performs better than the RFT (table 2).  Note that, as stressed in [20], what really matters to reproduce the kinematics of the motion correctly is the ratio between the tangential and the normal friction coefficients, K N /K T . The choice of the friction coefficients for the head is also only important in relative terms: the absolute values of C vx x and C vy y are irrelevant as far as they maintain the right proportion with one another and the flagellum coefficients, this is because the kinematics results from a force balance. The choice of using Perrin's formulae for the Brownian motion of an ellipsoid provides coefficients whose ratio is comparable to the ratio for an ellipsoid in a Stokes flow as given by (2.16) and (2.17) and whose absolute values (at least for the b/a ¼ 0.5 aspect ratio) are not very dissimilar (table 3). However, given the size and speeds involved, resorting to the Stokes law seems more physically based. Table 3 reveals that only the RFT-L model has both K T and K N /K T within the range indicated in [20]. While in [20] it is found that the RFT reproduces the trajectories reliably, in [18] it is reported that the drag values computed with RFT are inaccurate. This may not necessarily be in contrast since for the last calculations it is the absolute value of the coefficients that matters. However, it is most likely the case that while the customarily chosen RFT coefficients well fit the data for the swimming of a spermatozoon with reference parameters, they do not match the results for different geometries, e.g. larger or smaller ls or helical flagella as in [18].

Inaccuracy of treating the head and flagellum separately and further comparisons with the RFT solutions
Although being accurate compared with the RFT, our results indicate that the simplified model, consisting in evaluating the friction coefficients as the sum of the head and the tail contribution calculated independently, leads to notable errors. This approximation is convenient for first estimates since it reduces the computational cost when analytical solutions are available, e.g. for spherical or ellipsoidal heads, but neglects the interaction between the tail and the head. The system (2.8), in short, in this form each component of the vector of local forces f, being it located on the head or on the tail, can be split up into two contributions: one due to the points located on the head and one due to the points located on the tail. With superscript 'h' and 't' denoting the head and tail, respectively, and greek letters indicating the points on the surface of the swimmer we have the approximate approach neglects the head -tail interaction terms: J À1 a h b t u tot b t and J À1 a t b h u tot b h . For the case of a spherical head, we have verified that the approximate method overestimates the leading coefficients of equations (2.12) and (2.13), C vx x and C vy y , by approximately 23.5% and 17% and incorrectly captures many others (table 4). The difference between the resistive matrix coefficients calculated with the two methods is expressed by   [12,20] are divided by 0.7 mPa s) and ratio between the normal and longitudinal friction coefficients for the head and flagellum. Following [12,20] the head friction coefficient is computed for an aspect ratio b/a ¼ 0.5 and L head ¼ 10 mm. Perrin's formula and K N and K T as in [12] 5.7571 Â 10 25 1.14392 0.5429 1.89 Table 4. Difference between the resistive matrix coefficients and known terms of system (2.12) -(2.14) computed with the simplified assumption of representing the swimmer as the superposition of the head and tail flow field (resistive matrix (2.15)) and the full-swimmer representation. The error is computed according to formula (3.7). The second row in this table indicates the order of magnitude of the corresponding coefficient. where C refers to any coefficient and known term of equation (2.12) -(2.14). A positive value signifies that on average the simplified approach overestimates the coefficient. Remarkably, the overestimate of C vx x and C vy y leads to swimming trajectories with different radii as shown in figure 6 (compare the solid blue and dashed-green curve).
The simplified approach results in a different distribution of the forces on the surface of the flagellum as shown in figure 8; note that the values differ in particular on the left end where the flagellum is attached to the head. In the full-body model the values are zero because of the presence of the boundary, while in the simplified model the values are relatively large since this is a free end; however, they are not as high as on the tip of the tail given the lower velocities. Note again that even if inaccurate, the head þ tail model clearly outperforms the RFT (figure 8c), for which the force distribution only qualitatively resembles the first two models.
In figure 9, we plot in dashed-red the total forces and torque exerted by the head to the tail versus the swimmer linear and angular velocities. These results were obtained for the full-body regularized Stokeslet model calculation. We find that the points do not lie on a straight line as would be expected for the flow past an isolated sphere or ellipsoid, instead, they trace closed curves. The deviation from a straight line quantifies the contribution due to the head-tail interaction and further demonstrates the limitations of the approximate approach. Lines fitted to the dashed-red curves of figure 9 have slopes smaller than the theoretical coefficient 6pma by 8.5% in x-direction, larger by 9.8% in y-direction and smaller than 8pma 3 by 6.9% for the torque in z. Similar results hold when comparing the theoretical friction coefficients for prolate ellipsoids with the slopes of lines fitted to analogous curves for swimmers with ellipsoidal heads.

Sensitivity to the head shape
We study how swimming is affected by the head shape; in particular, we consider the case of swimmers with a prolate spheroidal head. We keep the head volume constant while varying the minor to major axis ratio b/a. In figures 5 and 6, we report results for ellipsoidal heads with an aspect ratio of 0.5 and 0.25 (dotted red and dashed-dotted black curves). Despite the identical beating movement, the tail integral of forces and torque varies as the head shape varies, given the different swimming velocities and head-tail interaction.
We have already quantified in figure 7 and table 2 the error of the RFT in predicting the absolute value of the velocity, here we report, in addition, a slight discrepancy between the RFT and the regularized Stokeslet model in identifying the fastest swimmer. In figure 10, we plot the net displacement and rotation angle Q for one period as a function of the minor to major axis ratio b/a. Note that for both values of K 0 the swimmer that swims the farthest is the one with b/a ¼ [0.375 0.5] for the Stokeslet model, b/a ¼ 0.25 for the RFT models. We have estimated the error bar associated with the chosen numerical resolution of RSM to be 0.42%. This value corresponds to the x-velocity component error of the RSM solution for solving the analytical problem reported in figure 1 at the same numerical resolution as the swimmer problem. The x-velocity has been selected as the most sensitive solution component since it corresponds to the maximum discrepancy between the RSM solution and the reference FEM solution (table 2). The swimmers that possess the largest v x average velocity according to the RSM are those with b/a ¼ [0.625 1], while for the RFT models the maximum v x average velocity is achieved for b/a ¼ 0.625.
The results in figure 10a,c reiterate that for the wavenumber of choice the RFT-L model is more accurate than the RFT-GH model. Note that both the RFT models badly capture the behaviour for the most stretched shapes. Observe also that the swimmer with a non-null curvature swims farther than the rectilinear swimmer (covers larger distances in one period).
For the K 0 ¼ 3375.5 rad s 21 swimmer, the net displacement and angle Q determine the curved trajectories reported on figure 10b. For a given displacement, larger angles correspond to trajectories with smaller radii, so that swimmers with more elongated heads display trajectories with larger radii. The behaviour is monotonic with b/a across all models ( figure 10d), surprisingly, the RFT-L model predicts the trajectory radius even better than the head þ tail model despite the lower accuracy in the velocities.

Flow field around the swimmer
In figure 11, we show the average flow field u f ,v f and the root mean square velocity fluctuations around the swimmer in the frame of reference of the swimmer, here the brackets denote the time average.  Figure 9. In dashed-red: surface integral on the swimmer flagellum of minus the x-and y-components of the force ÀF t x , ÀF t y (a,b) and z-component of the torque ÀT t z (c), versus, respectively, the swimmer velocity in x, y, and the angular velocity. In dotteddashed blue and dotted black, we show similar curves for the integral of the forces on the flagellum due to the beating motion only (right-hand side of system (2.12) -(2.14)), the dotted black case corresponds to the values in the coordinates of the diagonal system. All the variables in these plots are normalized by their maximum absolute value and for this reason are marked by a hat.
We compare the case of swimmers that draw circular and straight trajectories, and for the latter we report results for the swimmer with b/a ¼ 0.375. The swimmer flow field resembles the pattern produced by two counterrotating vortex dipoles, one centred at the head, the other centred at about x ¼ 4 Â 10 26 . Note that in figure 11a, which corresponds to the case of K 0 ¼ 7735.5 rad s 21 , the second dipole structure is offset with respect to the first one, this reflects the asymmetric beating movement visualized in figure 2 and responsible for the curved trajectory.
In figure 12, we display the average flow field profile for a cross section in x and y that spans the interval [0, 5 Â 10 24 ] m, we observe the r 22 scaling emerging at large enough distances for respectively the u f and v f component of the velocity. Close to the swimmer the velocity profile does not follow a clear power law and the u f component dominates.

Eigenvalues of the propulsive matrix system
The system of equations (2.12) -(2.14) can be diagonalized to remove the effect of the coordinate coupling. When doing so, the curves that on the v j À F B j , v j À T B j planes draw close loops collapse almost perfectly, and somewhat surprisingly, to a single line as expected for bodies of fixed shape (e.g. a sphere whose line slope on the v j 2 F j plane would be 6pma). See figure 9 and compare the dashed-dotted blue curve versus the black dots curve. The eigenvalues L x (t), L y (t), L z (t) and eigenvectors of the system are in general a function of time; however, the dotted-black curves of figure 9 can be well fitted by lines ÀF B x % L x v x , ÀF B y % L y v y , ÀT B z % L z v of slope L x ¼ 0.00011 N s m 21 , L y ¼ 0.00013 N s m 21 , L z ¼ 3.7 Â 10 214 N s m. These are not the slopes that appear in figure 9, where quantities are normalized by their maximum absolute values to allow for comparison. As for the eigenvalues, only the x and y axes appear to change their orientation in time while the z axis along which v and T B z are directed is fixed. In conclusion, despite the fact that the swimmer body goes through cyclical deformations, we are able to identify single time-independent friction coefficients L x , L y , L z that in opportunely rotated frames (the diagonal ones) relate F B j to v j , and v to T B z . This fact is consistent with the assumption of a quasi-steady flow described by the time-independent Stokes equations.   Figure 12. (a) Log -log plot of the u f and v f average velocity profile respectively along the y ¼ 0 and x ¼ 0 cross section on the z ¼ 0 plane for the K 0 ¼ 0 swimmer with a spherical head. The u f profile differs on the front and rear of the swimmer as marked by the arrows, while the v f profile is symmetric with respect to the y ¼ 0 line. Note the length over which the profile has been computed, one order of magnitude larger than the plots in figure 11, and the r 22 far field scaling. (b) Torque exerted by the flagellum due to the beating motion only versus the swimmer angular velocity for the diagonalized system and the cases: l ¼ l 0 (black), l ¼ l 0 /2 (blue), l ¼ 2l 0 (magenta), l ¼ 100 m (red). All the variables in these plots are normalized by their maximum absolute value.
However, for larger l: l ¼ 2l 0 and l ¼ 100 m, the latter representing the case l ! 1, the friction coefficients of the propulsion matrix show a more marked dependence on time since the points corresponding to the diagonalized system do not lie on a line but draw loops or S-shaped curves in the v j À T B j plane (figure 12b). Therefore, the flow is sensitive to the change of shape of the beating flagellum and this suggests that the quasi-steady and inertia-less assumptions may break down. It is important to point out that while the frequency Reynolds number is unchanged, the Reynolds number increases as l increases given the larger swimming velocities. The highest Reynolds number for the l ¼ 100 m case is Re ¼ 0.035.

Conclusion
We study numerically the motion of a flagellated microswimmer, specifically a sperm-cell swimmer, in an infinite domain. We first simulate locomotion by using two different numerical methods: the RSM and the FEM, finding a very good agreement between the two.
We find that the RFT performs reasonably well for swimming parameters close to laboratory observations: the normalized root mean square error for the swimming velocities is within 5-15% depending on the choice of the model parameters. However, the model is unreliable for smaller values of the wavelength, specifically, the case of wavelengths of about 1 4 and 1 2 of the total flagellum length, while the reference case has a wavelength of about one flagellum length. These results are consistent with the findings of previous studies that were focused on other types of flagella such as prokaryotic or bacteria-like type [18,19,36]. For example, in [18] it was reported that RFT fails to provide an accurate description of helical shapes relevant to bacteria, while in [36] it was concluded that for the broad range of geometry parameters of helical flagellated swimmers considered, RFT never gives accurate results. It was also pointed out that despite being unable to capture the full dynamics, RFT sometimes provides accurate solutions for single quantities [36]. This is an observation that we also made in reference to the RFT model with the choice of parameters suggested by Lighthill which is able to accurately predict the swimming trajectory radius. Finally, in [19] it was noted that when studying the complex geometry of superhelices, the experimental results are in excellent agreement with the calculations performed with the RSM, which is not the case of the RFT model. We also find that the RFT approach fails to correctly predict the optimum shape of the swimmer's head for fastest swimming in the rectilinear case.
We find that the simplified approach that consists in studying swimming as the linear superposition of the head and the tail contribution separately, referred to as the head þ tail model, leads to inaccurate results and we precisely quantify the error for the velocities, trajectory and force distribution. The inaccuracy of the method originates from the fact that the interaction terms between the head and the flagellum are neglected, or, else, from a physical perspective, the model fails to account for the frontrear symmetry breaking of the flow past the ellipsoidal head due to the presence of the flagellum. For a circular swimmer the radius difference between the trajectory of the head þ tail and full-body model is comparable to the radius difference between the trajectory of a swimmer with a spherical head and a swimmer with an ellipsoidal head of minor to major axis ratio 0.25. This suggests that this simplified approach as well as the RFT is unsuitable to perform optimization studies on hydrodynamically efficient body and beating shapes.
Finally, we have revealed by diagonalizing the propulsion matrix that the flow is rather insensitive to the cyclic deformation of the swimmer body for our choice of parameters K 0 , A 0 and l 0 . In fact, during one swimming cycle the points that correspond to different instant of time and therefore different shape configurations, lie on an almost perfect straight line on the v x -F B x , v y -F B y and v z -T B z plane, a behaviour typical of fixed-shape objects for which the friction coefficients are given. Consider, however, that the friction coefficient is approximately constant provided that the axes are opportunely rotated as the flagellum moves, this is where the analogy with fixed-shape objects ends. We have also observed that if we fix the frequency Reynolds number and change the wavenumber l of the flagellum travelling wave by choosing larger values, the curves on the v z -T z plane start displaying a markedly nonlinear behaviour. This suggests that for these cases the flow 'sees' the object deforming and hence could be prone to time-and inertia-dependent behaviours. This also calls for a more accurate definition of the relevant Reynolds number in the case of the flagellated swimmers compared with those two commonly used in the literature, which are based on the swimmer length or its beating frequency.
Data accessibility. Additional data is provided as Electronic Supplementary Material. Authors' contributions. C.R. designed the research, performed the calculations for the RSM and the RFT, analysed the data and wrote the manuscript. M.Z. performed the FEM simulations, analysed the data and contributed to the writing of the corresponding text. S.K. designed the research, supervised the project and co-wrote the manuscript. All authors have read and approved the manuscript before submission.
Competing interests. The authors declare no competing interests. Funding. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 703526.