Axisymmetric spheroidal squirmers and self-diffusiophoretic particles

We study, by means of an exact analytical solution, the motion of a spheroidal, axisymmetric squirmer in an unbounded fluid, as well as the low Reynolds number hydrodynamic flow associated to it. In contrast to the case of a spherical squirmer --- for which, e.g., the velocity of the squirmer and the magnitude of the stresslet associated with the flow induced by the squirmer are respectively determined by the amplitudes of the first two slip ("squirming") modes --- for the spheroidal squirmer each squirming mode either contributes to the velocity, or contributes to the stresslet. The results are straightforwardly extended to the self-phoresis of axisymmetric, spheroidal, chemically active particles in the case when the phoretic slip approximation holds.


INTRODUCTION
Self-propulsion of micro-organisms, such as bacteria, algae and protozoa, plays an important role in many aspects of nature. Whether a bacteria tries to reach a nutrient rich area or a sperm cell an unfertilized egg, motility often yields a substantial advantage over competitors. Due to their small size and velocity, the viscous friction experienced by the micro-organisms swimming in water is very strong compared to inertial forces [1,2]; consequently, they have developed swimming mechanisms adapted to these circumstances.
The motion of motile bacteria and other small organisms is typically induced by the beating of thin, thread-like appendages, the so called flagella; however, these species exhibit a broad morphology, in that they may posses either a single flagellum (e.g., monotrichous bacteria), several flagella (e.g., E. coli.), or a carpet of many small flagella, called cilia (e.g. Opalina). Focusing for the moment on the example of cilia covered microorganisms, it was proposed by Lighthill that the emergence of motility can be understood without a detailed modeling of the complex, synchronous beating of the cilia [3]. Instead, the effect of this beating pattern, i.e., the time-averaged surface flow induced by the envelope of cilia tips, has been modeled as providing a prescribed active flow velocity v (actuation) at the surface of the particle (both within the tangent plane of the surface, v s , and in the direction normal to the surface, v n ). The model, known by now as the "squirmer" model, was subsequently corrected and extended by his student Blake [4].
In its simplest and most used form, the tangential slip velocity of a spherical squirmer is taken to be the superposition of the fore-aft asymmetric and fore- * poehnl@hawaii.edu † popescu@is.mpg.de ‡ uspal@hawaii.edu aft symmetric modes with the slowest decaying contributions to flow in the surrounding liquid, i.e., the leading order modes. Combined, they make for a squirmer endowed with motility (from the leading fore-aft asymmetric mode) and with a hydrodynamic stresslet (from the leading fore-aft symmetric mode) [5][6][7][8][9]. Considerable progress has already been made in understanding the behavior of spherical squirmers, and many interesting questions have been answered, such as: what does the flow field around a squirmer look like and how it compares with the ones produced by simple microorganisms [4,10,11]; what happens if the swimmer is not in free space, but rather disturbed by boundaries [5] or external flows [8]; and how do pairs or even swarms of these particles interact [6]. Although this model can also be used to better understand the motion of microorganisms, e.g. Volvox [12], observed in experiments, its restriction to spherical swimmers limits a wider application. For example, Paramecia, one of the most studied ciliates, has an elongated body [13][14][15], which prevents a straightforward application of the traditional squirmer model. However, the aforementioned questions are also of interest for these and other elongated micro-organisms. A generalization of the squirmer model to non-spherical shapes is thus a natural, useful development.
Driven by advances in technology that allow increasingly sophisticated manufacturing capabilities, in the last decade significant efforts have been made towards the development of artificial swimmers [16][17][18][19][20]. The envisioned lab-on-a-chip devices and micro-cargo carriers [21,22], e.g., for targeted drug deliveries and nanomachines focused on monitoring and dissolving harmful chemicals [23,24], all need to perform precise motions on microscopic length scales. A better understanding of the general framework of microscale locomotion is required in order to optimally design and control such devices, and theoretical models facilitate new steps along this path. As one example, chemically active colloids achieve selfpropulsion by harvesting local free energy. They catalyze a chemical reaction in the surrounding fluid and propel due to the ensuing chemical gradient. Spherical "Janus" particles belong to this group, and it has been shown that the squirmer model can, in various circumstances, capture essential features of their motion [25]. However, there are many chemically active colloids with non-spherical shape, and rod-like particles are especially prevalent in experimental studies (see, e.g., [19,[26][27][28]).
It is well known that passive, non-spherical colloids exposed to an ambient flow exhibit significant qualitative differences -such as alignment with respect to the direction of the ambient flow [29], Jeffery orbits by ellipsoids in shear flow [30], nematic ordering arising from steric repulsion [31], or noise-induced migration away from confining surfaces [32] -in comparison to their spherical counterparts. It is thus reasonable to expect that endowing such objects with a self-propulsion mechanism will lead to rich, qualitatively novel dynamical behaviors, some of which may be advantageous for, while others may hinder, applications. Accordingly, it is important to develop an in-depth understanding of the shape-dependent behavior, and significant experimental and theoretical efforts have been made in this direction (see, e.g., Refs. [33][34][35][36][37][38][39][40][41][42] as well as the insightful reviews provided in Refs. [1,[43][44][45][46][47][48][49][50]).
Similar to the case of spherical swimmers, physical insight into the phenomenology exhibited by elongated swimmers can be gained from a corresponding squirmer model. For a model spheroidal swimmer moving by small deformations of its surface, Ref. [51] derived an analytic solution to the corresponding Stokes flow, from which the velocity of the swimmer could be calculated. In Ref. [52] it has been pointed out that an exact squirmer model for a spheroidal particle with a prescribed axi-symmetric, tangential slip velocity (active actuation of the fluid) on its surface can be written down by employing an available analytical solution for the axi-symmetric Stokes flow around a spheroidal object [53]. However, the approach has been used in the context of a somewhat restricted model particle, involving an additional fore-aft asymmetry of the surface slip velocity, because of the particular interest, for that work, in the question of "hydrodynamically stealthy" microswimmers. While capturing the self-propulsion velocity of the swimmer, this removes certain characteristics of the flow, inter alia those that contribute to the corresponding stresslet [54] and would allow distinguishing between, e.g., "puller" and "pusher" type squirmers.
We note that also a strongly truncated model, based on an ansatz for the slip velocity in the form of two terms which, in the corresponding limit of a sphere, reproduce the first two modes of the spherical squirmer model of Lighthill and Blake, has been discussed by Ref. [55]. This approach, however, is significantly affected by the fact that -as noticed in Ref. [51] and also discussed here (see Sec. IV) -for spheroidal squirmers both their velocity and their stresslet involve significant contributions from the higher order slip modes, in contrast to the case of a spherical squirmer for which only the first two modes contribute to those observables.
In this paper we employ the available analytical solution for axi-symmetric Stokes flow around a spheroidal object [53] to study the velocity and the induced hydrodynamic flow field around a spheroidal squirmer with a tangential slip velocity possessing axial symmetry, but otherwise unconstrained. The model squirmer is introduced in Sec. II. The series representation of the incompressible, axi-symmetric, creeping flow field around a spheroid [53] is succinctly summarized in Sec. III. In Section IV we discuss the velocity and the flow field corresponding to the spheroidal squirmer, with particular emphasis on illustrating the contributions from the modes of various order. Additionally, in Sec. IV A we discuss the straightforward extension of these results to deriving the flow field around a chemically active self-phoretic spheroid (for a similar mapping in the case of spherical particles see Ref. [56]). The final Sect. V is devoted to a summary of the results and to the conclusions of the study.

II. MODEL
The model system we consider is that of a spheroidal, rigid and impermeable particle immersed in an incompressible, unbounded, quiescent Newtonian liquid through which it moves due to a prescribed "slip velocity" (active actuation) at its surface (see Fig. 1). The slip velocity v s , which is tangential to the surface Σ of the particle (i.e., n · v s ≡ 0 on Σ, with n denoting the outer (into the fluid) normal to Σ), is assumed to preserve the axial symmetry and to be constant in time, but it is otherwise arbitrary. The surface slip v s is part of the model and thus it is a given function (or, alternatively as in, c.f., Sec. IV A, it is determined as the solution of a separate problem). There are no external forces or torques acting on either the particle or the liquid.
Due to the surface actuation, a hydrodynamic flow around the particle is induced and the particle sets in motion; we assume that the linear size L of the particle, the viscosity µ and the density of the liquid, and the magnitude |v s | of the slip velocity are such that the Reynolds number Re := |v s |L/µ is very small. Accordingly, after a short transient time a steady state hydrodynamic flow is induced around the particle and the particle translates steadily with velocity U (with respect to a fixed system of coordinates, the "laboratory frame"). (Owing to the axial symmetry, in the absence of thermal fluctuations, which are neglected in this work, there is no rigid-body rotation of the particle in this model.) The analysis is more conveniently carried out in a system of coordinates attached to the particle (co-moving system). This is chosen with the origin at the center of the particle and such that the z-axis is along the axis of symmetry of the particle (thus U = U e z ). The semiaxis of the particle are accordingly denoted by b z and b x (see Fig. 1); their ratio r e = bx bz is called the aspect ratio and determines the slenderness of the particle. The 1. (a) Two-dimensional (2d) cut of a prolate spheroidal particle with semi-major axis bz and semi-minor axis bx; also shown are the x and z axes of the co-moving system of coordinates and the unit vectors eτ and e ζ of the prolate spheroidal coordinate system; (b) 2d cut of a chemically-active selfphoretic particle; the chemically active part is the cap-like region shown in white, while the rest of the surface (black) is chemically inert. The quantity χ ("coverage") denotes a certain measure for the extent of the active part (see Sec. IV A in the main text).
In the co-moving system of coordinates, the flow v and the velocity U are determined as the solution of the Stokes equations where is the Newtonian stress tensor, with p denoting the pressure (enforcing incompressibility), I denoting the unit tensor, and () T denoting the transpose, subject to: -the boundary conditions (BC) and -the force balance condition (overdamped motion of the particle in the absence of external forces) The hydrodynamic flow v lab in the laboratory frame, if desired, is then obtained as v lab (P) = v(P) + U, with P denoting the observation point.
The boundary-value problem above can be straightforwardly solved numerically by using, e.g., the Boundary Element Method (BEM) (see, e.g., Refs. [8,[57][58][59]) as well as analytically. In the following, we shall focus on the analytical solution of this problem; the corresponding numerical solutions obtained by the BEM are presented in Appendix E and used as a means of testing the convergence of the series representation of the analytical solution. Before proceeding, we note that we will focus the discussion on the case of prolate shapes, i.e., r e < 1 (see Fig. 1(a)); the case of an oblate squirmer (r e > 1) can be obtained from the results for the prolate shapes via a certain mapping (see, e.g., Ref. [52] and Appendix A), while the case of a sphere is obtained from the results corresponding to a prolate spheroid by taking the limit r e → 1 − .

III. HYDRODYNAMIC FLOW AND VELOCITY OF A PROLATE SQUIRMER
The calculation of the hydrodynamic flow v induced by the squirmer and of the velocity U of the squirmer are most conveniently carried out by employing the modified prolate spheroidal coordinates ( [53]. These are defined in the co-moving system of reference, which has the origin in the center of the particle, and are related to the Cartesian coordinates via [53] x is a purely geometric quantity which ensures smooth convergence into spherical coordinates in the limit b x → b z . The unit vectors e τ and e ζ (see Fig.  1) are related to the ones of the Cartesian coordinates via and the corresponding Lamé metric coefficients are given by The iso-surfaces τ = const are confocal prolate spheroids with common center O; the surface of the particle corresponds to and the values τ > τ 0 (1 ≤ τ < τ 0 ) correspond to the exterior (interior) of the prolate ellipsoid, respectively. The coordinate ζ takes the values ζ = ±1 at the z = ±b z apexes, respectively, and ζ = 0 on the equatorial (x − y plane cut) circle. Note that, as shown in Fig. 1(a), e τ coincides with the normal n while e ζ and e φ span the tangent plane.

A. Velocity of the squirmer
The velocity of the squirmer can be determined from Eqs. (1) -(4) as a linear functional of the axisymmetric slip velocity v s = v s (ζ)e ζ without explicitly solving for the flow v. This follows via a straightforward application of the Lorentz reciprocal theorem [60,. For the case of the prolate spheroid with U = U e z , this renders (for details see, e.g., Refs. [52,54,61]) Therefore, U can be considered as known; accordingly, the boundary value problem defined by Eqs. (1) -(4) is specified and the solution v can be determined.

B. Stokes stream function and hydrodynamic flow
Since the boundary value problem defined by Eqs.
The general "semiseparable" solution for the stream function in prolate coordinates has been derived by Ref. [53] in the form of a series representation where G n and H n are the Gegenbauer functions of the first and second kind, respectively (see Ref. [62]). The functions g n (τ ) and h n (τ ) are given by certain linear combinations of G k (τ ) and H k (τ ), the coefficients of which are fixed by the corresponding boundary conditions (see below). The general solution above is applied to our particular system as follows. Noting that for our system the solution v(r) should be bounded (i.e., |v(r)| < ∞ everywhere), none of the terms involving the Gegenbauer functions of the second kind H n (ζ), which are divergent along the z-axis (i.e., ζ = ±1) [60,, can be present in the solution; accordingly, h n (τ ) ≡ 0 for all n ≥ 2. Furthermore, due to the same requirement of bounded magnitude of the flow, the terms involving the Gegenbauer functions of the first kind G 0 (ζ) = 1 and G 1 (ζ) = −ζ also cannot be present in the solution because they lead to divergences of v ζ at ζ = ±1 (see Eqs. (11) and (7)); accordingly, g 0 (τ ) ≡ 0 and g 1 (τ ) ≡ 0. Therefore, for our system only the functions g n≥2 will be of interest; these are given by [53] where the constants {C n , D n } n≥2 , {E n } n≥4 , and {F n } n≥2 , are fixed, as noted above, by requiring that the solution satisfies the BCs, Eqs. (3), and the force balance condition, Eq. (4).

IV. SQUIRMER AND SQUIRMING MODES
The form of Eq. (14e) suggests (for the expansion of the slip velocity v s ) the functions as a suitable basis over the space of square integrable functions f (ζ) satisfying f (ζ = ±1) = 0 (note that V n has a parametric dependence on τ 0 > 1). Defining, in this space, the weighted scalar product and choosing the weight as one infers (form the known properties of the associated Legendre polynomials) that indeed the set {V n } n≥1 is an orthogonal basis, Accordingly, the slip velocity function v s = v s (ζ)e ζ has a unique representation in terms of a series in the functions The coefficients of the expansion are written in the form above to ensure that in the limit of a sphere (b x → b z , i.e, τ 0 → ∞) one arrives at the usual form employed for a spherical squirmer, i.e., that of an expansion in associated Legendre polynomials P 1 n (cos ϑ) (see, e.g., Ref. [4,63]). This can be seen by noticing that in the limit τ 0 → ∞ one has τ 0 (τ 2 0 − ζ 2 ) −1/2 P 1 n (ζ) → P 1 n (ζ) and that b x → b z implies ζ → cos(ϑ); accordingly, it follows that, in the limit of the shape approaching that of a sphere, τ 0 V n (ζ)e ζ → −P 1 n (cos ϑ)e ϑ and the expansion in Ref. [4] is matched identically upon changing B n → For a given slip velocity, thus a given set of amplitudes B n of the slip modes V n (ζ), Eqs. (14d) and (20) evaluated for n ≥ 1 provide a system of coupled linear equations for the last unknown coefficients in the stream function, C n≥3 and D n≥2 . Inspection of this system reveals that it splits into two subsystems of coupled equations, one involving only the coefficients of even index, n = 2k, and the other one involving only the coefficients of odd index n = 2k + 1 (see Appendix C). Furthermore, from Eq. (20) it can be inferred that in the case that the slip velocity is given by a pure slip mode V n0 (ζ), i.e., B n = δ n,n0 for n ≥ 1, then the parity of n 0 selects one of the two subsystems. If n 0 is even, then all the coefficients C k and D k of even index k vanish and the stream function, Eq. (12), involves only the functions g k of odd index k, while for odd n 0 all the coefficients C k and D k of odd index k vanish and the stream function, Eq. (12), involves only the functions g k of even index k. Finally, we note that a pure slip mode V n0 (ζ) selects, as discussed above, either the odd or even number terms in the series representation of the the stream function, Eq. (12), but not only a single squirmer mode, as it is the case for the spherical squirmer. Accordingly, even simple distributions of active slip velocities on the surface of the particle can give rise to quite complex hydrodynamic flows around the squirmer.
With these general results, we can proceed to the discussion of prolate squirmers. We will focus on the usual quantities employed to characterize an axisymmetric squirmer [3,4,51], i.e., the velocity U of the squirmer, the magnitude S of the stresslet associated with the squirmer (which determines the far-field hydrodynamic flow of the squirmer), and the general characteristics of the hydrodynamic flow v(r) around the squirmer.
By combining Eqs. (14e), (19), (18), and (15), and noting that the polynomials P 1 n (x) are even (odd) func- tions of x for n odd (even), one arrives at Accordingly, it follows that (a) a squirmer may exhibit self-motility (i.e., U = 0) only if the slip velocity involves at least an odd index n slip mode V n ; (b) in contrast to the case of a spherical squirmer, for which the velocity is determined solely by the slip mode n = 1 irrespective of the details of the slip velocity v s , for a spheroidal squirmer all the slip modes of odd index contribute to the velocity (see also Fig. 2); consequently, (c) spheroidal squirmers with B 1 = 0 can be self-motile (due to contributions from other odd index slip modes), and spheroidal squirmers with B 1 = 0 can yet be non-motile if the contributions from other slip modes of odd index n precisely balance the contribution of the mode B 1 (which clearly pinpoints the shortcomings of a model with only two slip modes as in Ref. [55]).
In terms of the dependence on the slenderness parameter r e (which determines the value of τ 0 , see Eq. (8)), there are two findings. (d) for every slip mode n, the contribution |U n | (in absolute value) is a decreasing function of r e (see Fig. 2); second, (e) while at low values of the aspect ratio the contributions |U n | of the n > 1 slip modes are significant, the contributions from the modes n ≥ 3 decay steeply with increasing r e and become negligible, compared to |U 1 | (which remains non-zero), as the aspect ratio r e of the spheroid approaches that of a sphere (r e → 1 − ). This ensures a smooth transition into the spherical case, where, as previously mentioned, higher mode squirmers are not motile and U ∝ B 1 . Fi-nally, we note that, by comparison with the expression in Eq. (14c), one infers that the series in the last line of Eq. In what concerns the stresslet S(τ ), which (similarly to the case of a spherical squirmer) allows classification into pullers (positive stresslet, S > 0), pushers (negative stresslet, S < 0), and neutral swimmers (vanishing stresslet, S = 0), we note that it can also be expressed in terms of the slip velocity u s (ζ) as [54]: = n≥1,n even µB n S n (τ 0 ) where A denotes the surface area of the spheroid and As in the case of the velocity, there are certain significant differences from the case of a spherical squirmer (see Fig.  3): (a) a necessary condition for a prolate squirmer to exhibit a nonvanishing streslet is that at least one even (k = 2n) mode V k contributes to the slip velocity; (b) as for the velocity U (where more than a single mode contributes), the stresslet depends on all even squirmer modes; hence, (c) the stresslet contribution of B 2 = 0 can be offset or even inverted by other even, active modes B k = 0; (d) At a given aspect ratio r e , the contribution (in absolute value) |S 2n | is a decreasing function of n; and (e) while for elongated shapes (small values r e ) the contributions |S n | from the slip modes n ≥ 4 are significant, these contributions are steeply decreasing towards zero with increasing r e towards the value 1 − . In contrast, |S 2 | is increasing with r e . As in the case of the velocity, this behavior ensures the smooth transition to the case of a spherical shape, where S ∝ B 2 .
Since the stresslet is, by definition, the amplitude of the r −2 far-field term in the flow field (in the lab frame) of the squirmer, the series in the last line of Eq. (22) can be connected with one of the coefficients in the expansion of the stream function as follows. In the laboratory frame, which is related to the one (ψ particle ) in the co-moving frame via the slowest decaying term with the distance τ from the squirmer is − C3 90 G 0 (τ )G 3 (ζ); by Eq. (10), this term leads to a contribution ∼ C 3 /τ 2 to the flow. Accordingly, S ∝ C 3 and the pusher or puller squirmers (S = 0) indeed exhibit the expected far-field hydrodynamics, while for the neutral squirmers (S = 0) the far-field flow necessarily decays at least as ∼ 1/τ 3 .
Turning now to the flow field around the prolate squirmer, we will discuss separately the flow fields generated by the first few pure slip modes, i.e., the cases B n = δ n,n0 with n 0 = 1, 2, 3, 4; these flows are shown, in the laboratory frame, in Fig. 4. From the discussion above, we know that only a subset (either the odd index ones, if n 0 is even, or vice versa) of the terms in the series representation, Eq. (12), of the stream functions contributes to the flow. Since the metric factors are even functions of ζ, while G n (ζ) is an odd (even) function of ζ when n is odd (even), the flow has the following foreaft symmetries. For n 0 odd, the stream function involves the functions G k of index k an even number, and thus ψ(τ, ζ) = ψ(τ, −ζ); this implies v τ (τ, ζ) = −v τ (τ, −ζ) and v ζ (τ, ζ) = v ζ (τ, −ζ) (see figures 5 and 6); i.e., the z− (x−) flow components are fore-aft (anti)symmetric (see Fig. 4), and, accordingly, it contributes to the motility because it provides a "fore to back" streaming. Vice versa, for n 0 an even number one has v τ (τ, ζ) = v τ (τ, −ζ) and v ζ (τ, ζ) = −v ζ (τ, −ζ), i.e., the x− (z−) flow component are fore-aft symmetric (see Fig. 4); consequently, they cannot be associated to a motile particle.
Finally, we note that the similar analysis and results for the case of an oblate squirmer can be obtained from the available ones for a prolate squirmer via a simple transformation of the coordinates system and of the stream function (i.e., a mapping), as discussed in the Appendix A. As an illustration of using this mapping, the results shown in Fig. 4, corresponding to a prolate squirmer, have been used to determine the flows around the corresponding (i.e., of slenderness parameters r e = 1/r e ) oblate squirmers induced by the pure slip modes B n = δ n,n0 with n 0 = 1, 2, 3, 4; these are shown in Fig. 7. The analysis of the velocity U and stresslet S (see Figs. 8 and 9) for oblate spheroids leads to conclusions that are similar with those drawn in the case of prolate shapes. The only significant difference is that now for oblate squirmers the contributions of from the higher orders n slip modes decay to zero with decreasing aspect ratio r e → 1 + ; again, this ensures a smooth transition to the case of a spherical shape, where only B 1 or B 2 are relevant.

A. Self-phoretic particle
Squirmers with a wide range of active slip modes occur naturally in the context of model self-phoretic particles. One of the often employed realizations of such systems consists of micrometer-sized silica or polystyrene spherical particles partially coated with a Pt layer and immersed in an aqueous peroxide solution [59,65,66]. The catalytic decomposition of the peroxide at the Pt side creates gradients in the chemical composition of the suspension; as in the case of classic phoresis [67], these gradients, in conjunction with the interaction between the colloid and the various molecular species in solution, give rise to self-phoretic motility [68]. The mechanism of steady-state motility can be intuitively understood in terms of the creation of a so-called phoretic slip velocity tangential to the surface of the particle [20,67,69]. For such chemically active, axi-symmetric, spherical particles in unbounded solutions, the approximation of phoretic slip velocity leads to a straightforward mapping [25,56] onto a squirmer model; accordingly, the translational velocity of the particle and the hydrodynamic flow around the particle in terms of the phoretic slip can be directly inferred from the corresponding results in Ref. [4].
A similar mapping can be developed from a spheroidal, self-phoretic colloid to a spheroidal squirmer model for spheroidal, chemically active colloids. Following Ref. [61] the slip velocity of such particles can be written as where b is the so-called phoretic mobility (for simplicity, here assumed to be a constant) over the surface of the particle, and Q l (τ ) denotes the Legendre polynomial of the second kind [62]. The coverage χ = h−1 is defined in terms of the height of the active cap (measured, from the bottom apex, in units of b z ), i.e., χ = −1 corresponds to a chemically inactive spheroid and χ = +1 corresponds to all the surface being active. The coverage dependent coefficients c l (χ) (see Ref. [61]) describe the decoration of the surface of the particle by the chemically active element (e.g., Pt in the example discussed above) as an expansion in Legendre polynomials P l (ζ). Knowledge of these parameters allows us to formulate a slip velocity boundary condition similar to (20), i.e.
by identifying the effective squirmer modesB l = b c·τ0 c l Q l (τ 0 ), the desired mapping is achieved. As an illustration of this mapping, we show in Fig. 10 (left column) the flow fields for particles with parameters r e = 0.3, 0.5 and χ = 0, −0.3, respectively, and compare with the corresponding results obtained by direct numerical solutions obtained using BEM [64]. (Note that, if necessary, the accuracy of the analytical estimate can be systematically improved simply by increasing the order of the truncation in the series expansion of the stream function (see the Appendix D)).

V. SUMMARY AND CONCLUSION
We have studied in detail the most general axisymmetric spheroidal squirmer and we have shown that, in analogy with the situation for spherical shapes [56], model chemically active, self-phoretic colloids can be mapped onto squirmers. By using the semiseparable ansatz derived in Ref. [53] for the stream function, and representing the active slip (squirming) velocity of the squirmer in a suitable basis (Eq. (19), chosen such that in the limit of a spherical shape it smoothly transforms into the usually employed expansion for the classical spherical squirmer [4], the velocity of the squirmer, the stresslet  The main conclusion emerging from the study is that for spheroidal squirmers (or self-phoretic particles) the squirming modes beyond the second are, in general, as important as the first two ones in what concerns the contributions to the velocity and stresslet of the particle (and, implicitly, to the flow field around the particle, even in the far-field). The velocity is contributed by all the odd-index components (but none of the even-index ones) of the slip velocity; accordingly, in contrast with the case of spherical squirmers, it is possible to have spheroidal squirmers with a non-zero first mode and yet not motile, as well as ones missing the first slip mode and yet motile. Similarly, the stresslet value is contributed by all the even-index slip modes; thus, distinctly from the case of spherical squirmers, one can have pushers/pullers even if the second slip mode is vanishing, as well as neutral squirmers in spite of a non-vanishing second slip mode. Finally, even a single slip mode leads to a large number of non-vanishing terms in the series expansion of the stream function, and thus spheroidal squirmers with simple distributions of slip on their surface can lead to very complex flows around them (see Figs. 4 and 7).
This raises the interesting speculative question as whether the spheroidal shape is providing an evolutionary advantage; i.e., with small modifications of the squirming pattern, e.g., switching from a sole B 1 mode to a sole B 3 mode, a microrganism could maintain its velocity unchanged but dramatically alter the topology of the flow around it (compare the first and third rows in Fig. 4). In other words, compared to a squirmer with spherical shape, for which multiple modes must be simultaneously activated in order to change the structure of the flow, a spheroidal squirmer possesses simple means for acting in hydrodynamic disguise, which can be advantageous as either predator or prey.

A. Oblate spheroids
The similarity between oblate and prolate spheroids allows us to obtain the flow field around an oblate microswimmer of aspect ratio r e > 1 as a series in the oblate coordinates by using a mapping from the results, in prolate coordinates, corresponding to a prolate microswimmer with aspect ratio r e = 1/r e < 1 (and vice versa).

Analytical results
Numerical results (BEM) FIG. 10. The flow field (streamlines and velocity magnitude (color coded background)) induced by a chemically active prolate particle moving by self-phoresis, calculated either analytically via the stream function (left column) or numerically, i.e., by directly solving the corresponding Laplace and Stokes equations using the BEM [64] (right). The results shown correspond to the cases of half (χ = 0, top two rows) and less than half (χ = −0.3, bottom two rows) coverage, and two values, re = 0.3 and re = 0.5, of the aspect ratio, respectively. The arrows on the particles indicate the directions of their motion. The red area depicts the chemically active region. The motion of the particle and the direction of the flow corresponds to the choice b < 0; the characteristic velocity U0 is defined as in Ref. [61] Noting that the oblate coordinates and the metric factors can be obtained from the expressions of the corresponding prolate coordinates via the transformations [53] τ → iλ c → −ic , one concludes that the equation and boundary conditions obeyed by the stream function in oblate coordinates in the domain outside the oblate of aspect ratio r e = 1/r e can be obtained, by using the same transformation, from the ones in prolate coordinates outside a prolate of aspect ratio r e . Accordingly, it follows that ψ obl,r e = 1/re (λ, ζ) = ψ pro,re (τ = iλ, ζ). The flow field then follows as the curl of the stream function, i.e., v obl,r e (λ, ξ, ϕ) = ∇ × ψ obl,r e (λ, ξ)

B. Gegenbauer functions
The Gegenbauer functions of first kind G n , are also known as the Gegenbauer polynomials C α n with parameter α = −1/2 [62]; they are defined in terms of the Legendre polynomials P n as and fulfill the orthogonality relation δ n,m n, m ≥ 2 .
(30) The Gegenbauer functions can be related to the associated Legendre polynomials P 1 n , which are defined in terms of the Legendre polynomials P n as Using the relations x 2 −1 n d dx P n (x) = xP n (x) − P n−1 (x) and (n + 1)P n+1 (x) = (2n + 1)xP n (x) − nP n−1 (x), one arrives at the following relation between the l-th associated Legendre polynomial and the l + 1-th Gegenbauer function.
C. Decoupling of the even-and odd-index modes in the stream function expansion The coefficients C n and D n appear in the functions g k (τ ), entering the series expansion of the stream function at various indexes k (e.g., C n≥4 appears at both k = n − 2 and k = n); thus the infinite system of linear equations is strongly coupled. However, since the even and odd terms are not entering the same equations, the system splits into two decoupled subsystems, which are solved by using different methods.
The first subsystem consists of the even numbered terms in the series expansion of ψ(τ, ζ) and the set of conditions (Eqs. (14a), (14c), (14d) and (20)). Because each of the Eqs. (14a) and (14c) fix one of the even index coefficients (i.e., C 2 and F 2 ), the Eqs. (14d) and (20) evaluated at n = 2 involve only two unknowns, C 4 and D 2 , and thus can be solved as a sub-subsystem. With C 4 and D 2 known, the rest of the coefficients, up to the order N max at which the system is truncated (i.e., B k is set to zero for k > N max ), are solved iteratively by noting that Eqs. (14d) and (20) evaluated at k = n involve only two unknowns, the rest of the coefficients being already determined in the previous iterations up to k = n − 2.
The second subsystem includes all odd-index terms in the series expansion of ψ(τ, ζ) and Eqs. (14d) and (20). Complementing the previous case, only the even squirmer modes contribute. However, unlike in the case of the first subsystem, there are no lower level decouplings; accordingly, this subsystem is solved in the standard manner by truncation at the cut-off N max (above which all the coefficients are set to zero) and inversion of the resulting finite system of linear equations.
The choice of N max is done by varying the value of the cut-off and testing the changes in the coefficients. For the cases we analyzed in this work, a value N max = 16 was found to be sufficient.

D. Example of effects of too strong truncations for the case self-phoretic particles
To show the importance of the higher orders in the case of a self-phoretic swimmer, we compare the analytic results for keeping different numbers of squirmer modes B n in figure 11. Even for an aspect ratio r e = 0.8 (i.e., close to a spherical shape), the higher orders have significant influence on the flow field around the particle: e.g., compare the occurrence of a region of high magnitude flow near the point (0, −1.5) (behind the particle) instead of the correct location at near the point (0, 1.5) (in front of the particle).

E. Quantitative comparison to BEM
In figures 12 -14 we present a more detailed comparison of the results obtained by using the series representation (top row) with numerical results obtained by using the BEM to directly solve the governing equations (second row). In addition, the relative error, defined as ∆v = (v x,ana − v x,num ) 2 + (v z,ana − v z,num ) 2 is shown in the bottom row. The comparison confirms the expected quantitative agreement, with regions of sig-   In both cases the analytical results are shown in the top row, the numerical results (BEM [64]) in the middle row, and the relative error between the two (Eq. (31)) in the bottom row. The red area depicts the chemically active region.