Unsteady feeding and optimal strokes of model ciliates
Abstract
The flow field created by swimming microorganisms not only enables their locomotion but also leads to advective transport of nutrients. In this paper we address analytically and computationally the link between unsteady feeding and unsteady swimming on a model microorganism, the spherical squirmer, actuating the fluid in a timeperiodic manner. We start by performing asymptotic calculations at low Péclet number (Pe) on the advectiondiffusion problem for the nutrients. We show that the mean rate of feeding as well as its fluctuations in time depend only on the swimming modes of the squirmer up to order , even when no swimming occurs on average, while the influence of nonswimming modes comes in only at order . We also show that generically we expect a phase delay between feeding and swimming of 1/8th of a period. Numerical computations for illustrative strokes at finite Pe confirm quantitatively our analytical results linking swimming and feeding. We finally derive, and use, an adjointbased optimization algorithm to determine the optimal unsteady strokes maximizing feeding rate for a fixed energy budget. The overall optimal feeder is always the optimal steady swimmer. Within the set of timeperiodic strokes, the optimal feeding strokes are found to be equivalent to those optimizing periodic swimming for all values of the Péclet number, and correspond to a regularization of the overall steady optimal.
I Introduction
In order to be able to swim in viscous fluids, microorganisms must undergo nontimereversible sequences of shape changes referred to as swimming strokes (Lighthill, 1975; Purcell, 1977; Lauga & Powers, 2009). Through the noslip boundary condition, these strokes induce a net flow field around the organism and a distribution of viscous stresses which lead to locomotion. This swimminginduced flow also impacts hydrodynamic interactions with neighboring organisms (Drescher et al., 2009; Michelin & Lauga, 2010b) or material boundaries (Berke et al., 2008; Lauga et al., 2006), the overall dynamics of suspensions of cells (Kessler, 1986; Pedley & Kessler, 1992; Sokolov et al., 2007; Saintillan & Shelley, 2008a; Evans et al., 2011) and the feeding ability of organisms (Childress et al., 1987; Short et al., 2006).
Cellular motility is essential to many biological functions, from reproduction (Suarez & Pacey, 2006) to escaping agressions (Crawford & Purdie, 1992; Hamel et al., 2011). It also allows organisms to travel toward better local environments for example to seek (or escape) light, nutrient, or heat. The performance of the particular stroke displayed by a single microorganism, or that of a suspension of such swimmers, also results in the modification of the bulk stress and effective viscosity of a flow (Batchelor, 1970; Ishikawa, Simmonds & Pedley, 2007), or of its mixing properties (Saintillan & Shelley, 2008b; Leptos et al., 2009; Kurtuldu et al., 2011), an effect that is suspected to play an important role on largescale biomixing in the ocean for example (Doostmohammadi, Stocker & Ardekani, 2012).
The metabolism of many microorganisms relies on the absorption of diffusing nutrients present in their vicinity, ranging from dissolved gases and lowweight proteins, to more complex molecular compounds and, in the case of large organisms such as the protozoon Paramecium, smaller bacteria whose runandtumble motion is equivalent to a diffusive process at the scale of Paramecium (Berg, 1993). For a particular microorganism, the impact of the stroke on its feeding ability can be thought of as twofold: (i) through the motility resulting from the stroke, the organism can travel toward nutrientrich regions; (ii) by stirring nutrients in its immediate vicinity, the strokeinduced flow modifies, and possibly enhances, local concentration gradients.
The competition of advective and diffusive effects on the dynamics of a particular nutrient is quantified in the Péclet number, , where and are the characteristic diffusive and advective timescales respectively, where , and are the typical size of the organism, the characteristic flow velocity, and the nutrient diffusivity, respectively. Depending on the nutrient considered, Pe can vary by several orders of magnitude, even for a given microorganism.
Performing its stroke represents an energetic cost for the organism, as it must work against the fluid to overcome viscous dissipation. How far it can swim or how much nutrient it can absorb is therefore, in theory, limited by the finite amount of energy it has available. Considering that energy losses other than hydrodynamic can be accounted for by a fixed metabolic efficiency, the optimization of the swimming stroke to maximize either motility or feeding can therefore be formulated as follows: for a fixed amount of energy available to deform its shape, what is the optimal stroke of a particular microorganism maximizing either (i) the net displacement (optimal swimming problem) or (ii) the amount of a particular nutrient absorbed at the surface of the organism (optimal feeding problem)? In the latter case, the optimal stroke does not necessarily require a net displacement of the cell, as the organism can potentially just sit in a given location and stir the fluid around it. The optimal feeding stroke may also depend on the particular nutrient considered and the relative importance of advection and diffusion through the value of Pe.
The optimization problems described above are closely linked to the question of optimality with respect to a specific biological function, which can take two different forms: optimal shape or optimal gait. In the former, one is interested in the optimal morphology of the swimmer (e.g. its aspect ratio, the use of flagella vs. cilia,…) and compares different species of microorganisms. In the latter, the focus is placed on a given organism, and the goal is to determine the sequence of body deformations that performs best (Tam & Hosoi, 2007; Spagnolie & Lauga, 2010; Michelin & Lauga, 2010a, 2011; Tam & Hosoi, 2011a, b).
In this work, we focus on the optimal gait of a particular swimmer model, the socalled squirmer. This canonical model, consisting of a spherical microorganism imposing a tangential velocity at its surface, was introduced as a socalled enveloppe model for ciliated microorganisms (Lighthill, 1952; Blake, 1971). Ciliates, such as Paramecium, swim in viscous flows using the coordinated beating of a large number of small cilia distributed over their surface (Blake & Sleigh, 1974; Brennen & Winnet, 1977). In the squirmer model, the flow field can be determined analytically through the projection of the stroke on orthogonal squirming modes. Because of its simplicity, this model has been used to study a large variety of problems related to swimming microorganisms, including hydrodynamic interactions (Ishikawa, Simmonds & Pedley, 2006), mixing (Lin, Thiffeault & Childress, 2011), suspension rheology (Ishikawa & Pedley, 2007), collective dynamics and instabilities (Ishikawa et al., 2007; Evans et al., 2011), and feeding (Magar et al., 2003; Magar & Pedley, 2005; Doostmohammadi et al., 2012).
Recently, Michelin & Lauga (2010a) determined the optimal timeperiodic swimming strokes (i.e. those maximizing the swimming velocity for fixed energetic cost) of such a model microorganism, and identified their main properties. In a subsequent contribution, Michelin & Lauga (2011) considered the optimization of the stroke for feeding in the particular case of a steady surface velocity. Although such strokes correspond to nonperiodic displacements of the surface, the results shed some light on the link between swimming and feeding, and in particular it was shown that optimal swimming strokes and optimal feeding strokes were essentially identical regardless of Pe, a result that is not a priori intuitive due to the fundamental differences in the impact of swimming on feeding at low or high Pe: at low Pe, swimming only impacts marginally the nutrient distribution, but enables the organism to travel toward regions with richer nutrient content, while at high Pe, swimming also impacts feeding through stirring and strong advection of the nutrient in the vicinity of the organism surface.
The validity of these conclusions, and in particular the intimate relationship between optimal swimming and optimal feeding, remains however to be addressed in the general case of unsteady strokes. Magar & Pedley (2005) showed that in the particular limit of large Pe and small surface displacement, an equivalent steady problem could be defined. However, the unsteady effects of advection and diffusion in the general case of both finite swimmer displacement and finite Pe number remain unclear. In this paper, we specifically focus on the unsteady swimming problem. We first address analytically and computationally the link between unsteady feeding and unsteady swimming. We then derive, and use, an adjointbased optimization algorithm to determine the optimal unsteady strokes maximizing feeding rate for a fixed energy budget.
The paper is organized as follows. In §II, the squirmer model is briefly presented, and the swimming and feeding problems are posed mathematically. In §III, the unsteady feeding rate is determined in the asymptotic limit of small Pe. The impact of the swimming stroke and of the Péclet number on the feeding rate is further analyzed in §IV using numerical simulations, providing an important insight on the link between swimming and feeding. Section V presents the result of the stroke optimization with respect to feeding and conclusions and perspectives are finally presented in §VI.
Ii Swimming and feeding of a model ciliate
ii.1 The squirmer model
The present work focuses on a particular model microorganism, the squirmer, illustrated in figure 1. It is a spherical organism of radius which prescribes periodic tangential deformations of its surface with a frequency , in order to swim in a viscous fluid of dynamic viscosity and density . The present analysis is restricted to purely axisymmetric deformations of so that the swimming velocity is parallel to the axis of symmetry , with no rotation. In this paper, we will seek optimal strokes maximizing the feeding rate of the organism for a given amount of energy available during each period to perform its surface deformation (and possibly its swimming). This average rate of energy consumption, , is identified with the rate of work applied on the fluid by the swimmer at its surface, or, equivalently, the total mechanical energy dissipated in the fluid through viscous effects during one period. It is related to the typical surface velocity scale by
(1) 
The squirmer is swimming in a continuous suspension of a given nutrient (e.g. bacteria, large proteins/molecules, heat…) characterized by a farfield concentration and a diffusivity , and advected by the flow created by the surface stroke. On the swimmer boundary, the nutrient is instantaneously absorbed and processed at the surface so that , with the equilibrium concentration at the surface determined by the processing mechanism. A more realistic, but more complex, boundary condition was proposed by Magar et al. (2003) and Magar & Pedley (2005), taking into account such effects as the resistance of the membrane to nutrient absorption, and the finite diffusion and processing time of the nutrient within the cell. The instantaneous nutrient uptake by the organism through diffusion at its boundary, , is given by
(2) 
In the case of a purely rigid sphere, with no advection, a steady nutrient flux is achieved through diffusion . In the following, we focus on the modification of the concentration field by the organism and define the rescaled concentration field .
Three distinct timescales are present in the problem: (i) a diffusive timescale , (ii) an advective timescale and (iii) the stroke period , while only the latter two were present in the purely swimming problem (Michelin & Lauga, 2010a) and only the first two in the steady feeding problem (Michelin & Lauga, 2011). The Péclet number, , is a measure of the relative importance of advective and diffusive effects near the surface of the squirmer, and is equal to
(3) 
A second independent timescale ratio can be defined either as a characteristic of the stroke, for example the relative velocity , or as a periodbased Péclet number . In the following, all equations and quantities are nondimensionalized using , , , and as reference quantities.
Swimming problem
Due to the small size of the organisms considered, the Reynolds number, , a relative measure of inertia and viscous effects in the flow, is always much smaller than one, and the velocity and pressure fields satisfy Stokes equations. The swimming problem in the reference frame attached to the organism is therefore
(4) 
with the boundary conditions on the swimmer surface and at infinity given by
(5)  
(6) 
Note that the prescribed surface field, , is the stroke imposed by the organism and at the origin of both locomotion and stirring. The stroke is assumed to be axisymmetric, therefore the surface velocity only depends on and , with the polar angle measured from the swimming direction (figure 1). In Stokes flow, the swimmer can not sustain any net hydrodynamic force, therefore we have
(7) 
where is the unit normal vector pointing into the fluid ( here). Note that we have assumed the swimmer to be neutrally buoyant. The solution to the swimming problem in (4)–(7) is obtained by decomposing the surface velocity onto the squirming modes (Blake, 1971; Michelin & Lauga, 2010a)
(8) 
with
(9) 
where is the th order Legendre polynomial. The values of the pressure field and streamfunction are then obtained as
(10)  
(11) 
with
(12)  
(13)  
(14) 
In the decomposition above, the first mode is the only one that contributes to the swimming motion (we have for all times) and is referred to as the swimming mode, or “treadmill”. All remaining modes (including the socalled stresslet, , characterizing the modification of the bulk stress by the swimmer) correspond to higher order singularities in the farfield flow and do not contribute to the swimming motion.
The dimensionless energetic cost, , is computed as (Michelin & Lauga, 2010a)
(15) 
with
(16) 
and is equal to the rate of working of the squirmer on the fluid through its boundary actuation or, equivalently, to the total energy loss through viscous dissipation in the fluid domain. In the following, we define as the timeaveraging operator over one stroke period. With this definition, is the typical nondimensional surface velocity of the swimmer. Following Lighthill (1975), the stroke swimming efficiency,  or scaled energy cost , is defined as the ratio of the energetic cost of pulling a rigid sphere with constant velocity and the energetic cost of swimming at the same average velocity, obtained here as (Michelin & Lauga, 2010a):
(17) 
Feeding problem
To evaluate the amount of nutrient absorbed at the surface of the organism, the nondimensional advectiondiffusion problem must be solved
(18) 
together with the farfield behavior and purely absorbing boundary conditions on the swimmer surface (figure 1)
(19)  
(20) 
In equation (18) the parameter can also be understood as the periodbased Péclet number. The flow field, , originates from the organism stroke and is obtained from the squirming mode amplitudes, , using (11), (13) and (14). The feeding performance of the stroke is evaluated using the ratio quantifying the net gain in nutrient uptake in comparison with the purely diffusive case (). The relative nutrient flux, , is therefore nondimensional and given by
(21) 
Eulerian vs. Lagrangian description
A given periodic stroke, be it swimming or nonswimming, can be mathematically described following two different approaches:

By prescribing at each instant, a periodic surface velocity on each point fixed in the swimmer frame, , or equivalently a set of functions . We will refer to this description in the following as the Eulerian periodic stroke.

By prescribing periodic trajectories, , of material surface points labeled by their reference position on the sphere . We will refer to this description in the following as the Lagrangian periodic stroke. The surface velocity and mode amplitudes, , can then be obtained from as (Michelin & Lauga, 2010a)
(22) (23)
In both descriptions, the flow velocity is periodic and completely determined by the periodic functions . However, in the Eulerian formulation, material surface points do not necessarily have periodic trajectories. Indeed, periodic Lagrangian strokes only represent a subset of periodic Eulerian strokes, namely the ones guaranteeing that every surface point comes back to its original position at the end of a full stroke period. Despite its shortcomings regarding the description of material point trajectories, the Eulerian approach has been the most popular for models of swimmers because of its simplicity, and in particular the possibility to consider steady strokes corresponding to steady surface and flow velocities (Ishikawa et al., 2006; Short et al., 2006; Doostmohammadi et al., 2012; Evans et al., 2011, to cite only a few).
ii.2 Optimal swimming and optimal feeding
For a given amount of energy available to perform a periodic stroke, an organism might have different optimal surface motions depending on the biological function of interest: migration (swimming problem) or nutrient uptake (feeding problem). A priori, those two objectives should lead to different optimal strokes, if anything because the optimal feeding stroke may depend on nutrient diffusivity through the value of Pe while the swimming problem does not depend on it.
As emphasized earlier, a periodic stroke can be defined in two different ways, either from an Eulerian point of view (periodic flow field) or from a Lagrangian point of view (periodic material displacement). In our recent contributions, we presented the result of the optimal swimming problem (for both Eulerian and Lagrangian strokes) (Michelin & Lauga, 2010a) and of the optimal feeding problem in the Eulerian steady framework only (Michelin & Lauga, 2011). A brief summary of these results is first presented here.
We start by remarking that, for the swimming problem, Eulerian optimal strokes are necessarily steady and each mode, , is independent of time. This is a direct consequence of the absence of history effect in the swimming problem: the swimming velocity and the energetic cost only depend on the instantaneous surface velocity. The optimal Eulerian stroke is then obtained by choosing the surface velocity distribution maximizing instantaneously the efficiency . From (17) we see that the Eulerian optimal swimming stroke is simply obtained by putting all the energy into the swimming mode, namely . The resulting treadmill swimmer, with an efficiency , is therefore the overall optimal for locomotion (Leshansky et al., 2007; Michelin & Lauga, 2010a).
In the case of the feeding problem, the presence of a timederivative in the advectiondiffusion equation introduces history effects, and the optimal Eulerian feeding stroke is therefore not necessarily steady. Focusing on the simplified problem of steady strokes, Michelin & Lauga (2011) showed using adjointbased optimization that the optimal steady feeding stroke is essentially the same as the optimal steady swimming stroke, a result which, surprisingly, remains true for all Péclet number.
That result was not obvious a priori. The value of the mean feeding rate of the organism for a given stroke is a strong function of the diffusivity of the nutrient whose distribution around the organism is qualitatively different in the diffusive and advective regimes (Magar et al., 2003; Michelin & Lauga, 2011). The optimal feeding rate, , depends strongly on Pe, but the stroke to achieve this optimal value does not. This result is important biologically as it implies that, for a given organism, a unique optimal stroke maximizes the nutrient uptake regardless of the details of its diffusive transport. For all Pe, and in the steady Eulerian framework, maximizing feeding and maximizing swimming are therefore equivalent problems.
Although simpler conceptually and mathematically, the Eulerian framework is not appropriate to describe periodic deformations of a material surface, such as, for example, the strokes of ciliated cells. To impose periodicity of the surface motion, it is necessary to turn to the Lagrangian approach and to consider the unsteady swimming and feeding problems. Michelin & Lauga (2010a) showed numerically that the optimal Lagrangian swimming stroke could be decomposed into two different parts: an effective stroke, dominated by the swimming mode, , and producing a forward velocity, and a recovery stroke during which material points (e.g. cilia tips) are brought back to their original position with frontlike dynamics to minimize their (negative) impact on the swimming velocity. This front, or wave, is reminiscent of metachronal waves observed on ciliated organisms (Brennen & Winnet, 1977) and results from a small phaseshift in the motion of neighboring surface points leading to global symmetrybreaking at the wholeorganism level. When the squirmer model is used to represent a ciliate, the cilia length constrains the maximum displacement of the surface and therefore limits the ability of the swimmer to not only approach the optimal Eulerian stroke (treadmill) during the effective stroke but also to reduce the impact of the recovery stroke on the swimming motion. Using a constrained optimization algorithm, the direct relationship between swimming efficiency and surface displacement amplitude was obtained, and Michelin & Lauga (2010a) showed that the optimal efficiency of could be reached asymptotically.
The optimization of the Lagrangian feeding stroke however remains at this point an open question; it is the focus of the present paper. The analysis of the nutrient uptake is first addressed analytically at small Pe. The general unsteady feeding problem is then considered numerically before turning to its optimization.
Iii Unsteady feeding at low Pe: Asymptotics, scalings, and optimum
In this section we focus on the feeding problem in the asymptotic limit of dominant diffusion (). For a given stroke, this is equivalent to the asymptotic analysis of the advectiondiffusion problem in the limit .
iii.1 Steady and unsteady boundary layers
For a steady velocity field, finding the asymptotic expansion of the scalar concentration, , and surface flux, , in the limit corresponds to a variation on the classical mass transfer problem near a sedimenting sphere (Acrivos & Taylor, 1962; Magar et al., 2003; Michelin & Lauga, 2011). It is based on matching two different solutions for the scalar field : near the surface of the sphere, diffusive effects are dominant, and advection only appears as higher order corrections, while in the farfield, a balance of both advection and diffusion leads to the proper decay of .
In the case of an unsteady velocity field, both terms on the left handside of (18) do not have the same scaling in the farfield. As a result the decay of the concentration field at infinity is not the same whether one considers the timeaverage of or its fluctuations around the mean, and a double boundary layer problem must be considered:

in the near field, , diffusion dominates and the absorbing boundary condition () at the surface of the swimmer is satisfied;

in the unsteady boundary layer (UBL), , a balance between diffusive effects and rate of change of the local concentration ensures the proper farfield decay for the timedependent fluctuations of the concentration field ;

in the steady boundary layer (SBL), , a balance between advection by the steady velocity field and diffusion ensures the farfield decay of the timeaverage concentration .
iii.2 Asymptotic problem formulation
Decomposing the mode amplitudes, , as well as the concentration field, , and feeding rate, , into their Fourier components, we write
(24) 
The advectiondiffusion equation becomes then
in the near field,  (25)  
in the UBL,  (26)  
in the SBL,  (27) 
In (25)–(27), the following linear operators have been defined
(28)  
(29)  
(30)  
(31) 
and (resp. ) is identical to in (28) after replacing by (resp. ), and is defined as after replacing by . The following boundary conditions must also be satisfied:
(32)  
(33)  
(34) 
iii.3 Matched Asymptotic Expansion
A regular series expansion in of , and is then performed up to . We write
(35) 
At each order, and are to be matched for and , while and are to be matched in the limit and .
Here, the nonhomogeneous forcing (32) only acts on the steadystate component of the concentration field, and is transmitted to the timedependent components by advection. Therefore, from the scalings of the different terms in equations (25)(26),
(36) 
Order
At this order, advection is neglected and the solution is simply the steady diffusive solution , which satisfies both nearfield and farfield boundary conditions. Therefore, for all . The resulting feeding rate is
(37) 
Order
Using (36), and for all . The steady components , , and satisfy
(38)  
(39)  
(40) 
Solving these equations and matching , and up to leads to
(41) 
and the resulting feeding rate remains unmodified at this order.
Order
Next, the advection diffusion equation is expanded up to in each region.

In the near field, :
(42) whose general solution satisfying the nearfield boundary condition, , is obtained as
(43) where are constants to be determined after matching with the UBL solution.

In the unsteady boundary layer, :
(44) whose general solution compatible with the boundary condition at infinity (for ) is
(45) (46)
where the are constants to be determined in the matching process. Matching , , and , up to leads to
(48)  
(49)  
(50) 
and the resulting feeding rate expansion is
(51) 
Up to this order, we see that the results of the classical lowPe asymptotic expansion for a steady velocity field are recovered and the mean feeding rate only depends on the average swimming velocity. In order to capture the leading order unsteady contribution to the feeding problem, the expansion must be carried out to the next order.
Order
From (21), we see that only the computation of the azimuthal average, , of the th Fourier component of the concentration field
(52) 
is necessary in order to compute the correction to the nutrient uptake.

This equation can be solved explicitly for using the farfield boundary condition for the nonconstant Fourier components and we get
(55) (56) with defined for using .

In the steady boundary layer, the equation for is identical to that at the previous order and the general solution takes the same form, see (47).
iii.4 Discussion
The asymptotic analysis obtained in (57)–(59) provides some important physical insight into the relationship between the swimming motion and the nutrient uptake on the surface of the swimmer. As for the steady case, the leading order advective correction to the feeding rate is linear in Pe and only depends on the average velocity of the organism (Acrivos & Taylor, 1962; Magar et al., 2003). At this order in Pe, there is a direct correlation between swimming and feeding and only the mean feeding rate is modified, fluctuations in time being negligible (higher order).
The next order correction marks a fundamental difference between the steady and unsteady problems: in the steady case, all squirming modes contribute to the next correction at order (Michelin & Lauga, 2011). Instead, in the unsteady feeding problem, a new correction to (both its mean value in time and fluctuations) appears at order , which depends solely on the swimming velocity of the organism (through all the Fourier components, , of the swimming velocity, , with no other squirming modes), and dominates the contribution of nonswimming modes that will only enter at order . For all timeperiodic strokes, the instantaneous feeding rate is therefore completely determined up to by the characteristics of the swimming velocity of the organism.
This result has a major consequence for strokes that swim instantaneously () but do not swim on average (). In this case, the leadingorder improvement to the feeding rate is solely governed by the zeromean fluctuations of . Nonswimming modes only contribute to higher order corrections, even if they have nonzero time averages. Consequently, for an organism that does not have a net swimming motion (e.g. a timereversible swimmer), an instantaneous zeromean swimming motion still presents a feeding advantage over stirring strokes where the cell stays in the same position at each instant ().
Our asymptotic expansion also provides some information on the relative phase of swimming and feeding. For an unsteady swimming velocity, , with a single dominant Fourier component, the instantaneous feeding rate has a delay on the swimming velocity (since in 58). A maximum in the feeding rate is therefore expected to take place after the peak swimming velocity, with a delay of 1/8 of a period.
Note that the total nutrient flux is fully determined by the body velocity up to . Whether the organism is swimming (forcefree) or is an actuated rigid sphere (forced motion) does not actually come into play here. All the conclusions above are therefore valid for nonbuoyant swimmers, but also for oscillating rigid spheres in Stokes flow, for which the present results represent a generalization of classical steady mass transfer results (Acrivos & Taylor, 1962) to unsteady motions (see Appendix B for more details).
In summary, our analytical results show that for low Pe, feeding is completely determined by swimming for any periodic stroke. Optimization of the feeding rate for a fixed amount of available energy is therefore equivalent in this limit to maximizing the swimming velocity under the same constraint, namely the swimming efficiency optimization problem. At low Péclet number, the Lagrangian optimal swimming and optimal feeding strokes are therefore identical, which confirms the result obtained in the steady framework by Michelin & Lauga (2011). In addition, similarly to the result for swimming, we get the result that at low Péclet number the optimal unsteady feeding problem is actually steady. This can be seen from (57) where the steady Fourier mode, , carries a higher weight than the other Fourier components compared to their relative importance in the rate of working.
Iv Unsteady feeding at finite Pe: simulations
To confirm the lowPe results obtained analytically, we now turn to characterizing the feeding performance of different strokes for intermediate and large Pe. Eulerian periodic strokes are determined by prescribing for all , while Lagrangian periodic strokes are described by giving the trajectories of material points where is the current position of the material point and its mean position. Alternatively, those strokes will be defined by , with . For illustration we consider three particular swimming and nonswimming Lagrangian periodic strokes:

Stroke A is the numerical optimal swimmer identified in Michelin & Lauga (2010a) which has swimming efficiency ;

Stroke B is a less efficient swimmer obtained using surface deformations in the form of a simple progressive wave:
(60) with and ;

Stroke C takes the same form as stroke B but with and . Stroke C represents a timereversible (or “reciprocal”) deformation, and therefore has no net swimming motion, .
All three strokes display nonzero instantaneous swimming, but only strokes A and B show swimming on average. Stroke C differs thus from purely stirring strokes for which the organism is strictly still at each instant. The trajectories of material surface points are shown for strokes A, B and C in figure 2. Mathematically, from the knowledge of , the mode amplitudes are obtained using (23).
iv.1 Numerical solution of the advection diffusion problem
For a given set of mode amplitudes, , the advectiondiffusion equation in (18) is solved spectrally in time for each azimuthal component of the concentration field
(61) 
The functions satisfy therefore the following systems of ordinary differential equations for and :
(62) 
with boundary conditions
(63)  
(64) 
In (62), and are thirdorder scalar tensors defined in Appendix A. Equations (62)–(64) are discretized on an exponentiallystretched grid in to concentrate points near the surface of the swimmer (see Michelin & Lauga, 2011, for more details), and the solution is then found iteratively. In typical simulations, the resolution used was points for the grid, – Legendre polynomials for the azimuthal dependence, – points in time, and – squirming modes to describe the swimming stroke.
Alternatively, the advectiondiffusion equation can be marched in time for each azimuthal component, , using an explicit timestepping scheme for the advective terms and Crank Nicholson for the diffusion term. In the following, the advectiondiffusion equation is solved spectrally in time except for strokes that do not swim on average (e.g. stroke C) for which the iterative algorithm does not converge properly or fast enough, and the timemarching approach is used in that case.
Computationally, it is observed that the instantaneous nutrient flux converges rapidly with the number of squirming modes used to represent the swimming stroke, as shown in figure 3. The convergence is even faster for the average nutrient flux: describing stroke A with only the first two squirming modes significantly speeds up the computations while introducing an error smaller than on the average feeding rate. Similar numerical tests performed on less efficient swimmers than stroke A (that is, swimming strokes for which mode 1 is not dominant) did not modify this observation significantly, and restricting the computation to only 2 or 3 squirming modes typically introduces an error smaller than . This rapid convergence of the mean and fluctuating feeding rate is yet another indication that the swimming motion controls the feeding ability of the organism and higherorder modes only act as a small correction to the average feeding rate.
iv.2 Impact of the swimming stroke on the feeding performance
Figures 4 and 5 show the concentration field around the squirmer for five successive and equispaced instants of a full period, for (figure 4) and (figure 5), and for the three different strokes. For strokes A and B, at lower Péclet number, the nutrient concentration field only shows a weak front–back anisotropy as diffusion dominate over advection, confirming the observations on steady strokes of Magar et al. (2003) and Michelin & Lauga (2011). As Pe is increased, sharper concentration gradients can be seen on the front of the squirmer. This results in an increased average feeding rate for increasing Pe as was observed for steady strokes (Michelin & Lauga, 2011). The main difference with the steady results is that in the unsteady scenario, the velocity of the squirmer changes (and possibly reverses sign) inducing a fluctuation in this frontback anisotropy and in the boundary layer thickness. For stroke C, which does not swim on average, the nutrient concentration field shows a strong isotropy, even at larger Pe, with much weaker concentration gradients resulting in a very weak modification of the nutrient uptake .
Stroke  

A  
B  
C 
Comparing the results obtained for the different strokes in Table 1, we see that stroke A is clearly more efficient than strokes B and C from a feeding point of view, and stroke A also corresponds to a “better” swimmer. This is consistent with the increase of the feeding rate with the instantaneous swimming velocity that enables the formation of sharp concentration gradients in front of the squirmer. For stroke C, the periodic reversal of the swimming velocity over the period, and the absence of net displacement, results in the impossibility to maintain sharp concentration gradients at the front of the body and to swim toward regions with richer nutrient content, reducing its feeding ability significantly.
Looking at the temporal variations of the swimming velocity and feeding rate throughout the stroke period (bottom frames of figures 4 and 5), a phase delay between the former and the latter is clearly identified for stroke A and B, and for all Pe considered. For stroke C, a similar delay is observed between the peaks in velocity magnitude (positive or negative) and the peaks in feeding rate: for this stroke, the feeding rate frequency is twice that of the swimming velocity because of the exact symmetry between the two half stroke periods. The presence of this time delay in all strokes is consistent with the results of the lowPe asymptotic analysis in §III and can be interpreted as the time necessary for the concentration gradient (and possibly boundary layer) to reestablish at the front of the cell when its velocity starts increasing again.
iv.3 Impact of the Péclet number on the feeding performance
It was observed previously that the value of Pe plays an important role in the feeding ability of the cell. This is investigated further here by looking at the impact of Pe on the instantaneous feeding rate for strokes A, B and C. The instantaneous feeding rate, , is decomposed into its mean value, , the amplitude of its fluctuations in time, , and its normalized profile, , so we write
(65) 
where and . Similar quantities are also defined for the swimming velocity: , , and . For a given stroke (A, B or C), the variation of these three quantities with Pe is displayed in figure 6.
For swimming strokes, it is observed that, for low Pe, the modification in the mean feeding rate, , scales linearly with Pe (strokes A and B). This is consistent with the asymptotic analysis of Section III and with the steady results in Michelin & Lauga (2011). In such a diffusiondominated regime, swimming enables the cell to sweep a region of fresher nutrients with an effective crosssection radius that is independent of the swimming velocity (because of the predominance of diffusion) and of the order of the size of the cell. At higher Pe, the reduced importance of diffusion over advection reduces the effective crosssection radius and increases at a lower rate with Pe. For strokes with no net swimming motion (stroke C), the modification in the mean feeding rate scales as a higher power, , for , consistently with the results of the asymptotic analysis.
For both swimming and nonswimming strokes, the amplitude of the feeding rate fluctuations, , varies as for , consistently with our asymptotic results. On figure 6 the fluctuations profile, , is also represented and compared to the leading order prediction of the asymptotic analysis. We see a very good agreement at low Pe which persists even at high Pe for efficient swimming strokes such as stroke A. This confirms that the feeding rate (both its mean value and its fluctuations) is determined at leading order by the swimming mode and corrections from the other modes only play marginal roles. Again, a clear phase delay between the swimming velocity and feeding rate is observed for all Pe, and for the least efficient swimmers considered (B and C), this delay seems to increase with Pe.
When Pe becomes large, another significant difference appears between strokes with zero (stroke C) or nonzero (strokes A and B) mean swimming velocity. For strokes A and B, the average feeding rate continues to increase with Pe, albeit more slowly. From the largePe steady results by Michelin & Lauga (2011), we expect to scale as , when the increase in feeding rate with swimming is driven by the concentration boundary layer thickness around the cell. In contrast, for nonswimming strokes, reaches a maximum for a finite value of Pe () beyond which an increase in Pe actually results in a decrease of the feeding rate. Moreover, beyond a second critical value ( for this particular stroke), the mean feeding rate falls below , and for large Pe, swimming actually penalizes feeding as it reduces the net feeding rate below the level of the purely diffusive regime (). This somehow surprising result can be understood as follows. In stroke C, the sphere swims forward during half of a period leaving behind it a nutrientdepleted wake. In the second half of the stroke, the cell swims backward into this region of poor nutrient concentration, resulting in a reduced flux at the boundary.
iv.4 The optimal unsteady stroke is steady
As we discussed above, the optimal Eulerian swimming stroke is necessarily steady. The same conclusion can not be drawn a priori for the feeding problem due to the timedependence of the advection diffusion equation (see §II.2). We saw however that it was true analytically at low Péclet number. Numerically, it also seems to hold as illustrated in figure 7. We performed numerical simulations on a large collection of unsteady Eulerian periodic and Lagrangian periodic strokes (8500 in total), ranging from very efficient to poor swimmers. For all values of Pe, the feeding rate is seen to be always less than that obtained with the optimal steady feeding stroke (treadmill). As for the optimal swimming stroke, the optimal Eulerian unsteady feeding stroke must therefore also be steady. Furthermore, figure 7 demonstrates that the more efficient the unsteady stroke is for swimming, the closer it can get to the optimal feeding rate.
iv.5 Feeding and swimming
In the previous sections, a relationship between the swimming velocity and the feeding rate was clearly identified suggesting that at leading order, the mean feeding rate is determined by the swimming velocity and Pe. More precisely, and in the light of the steady results of Michelin & Lauga (2011), one expects the feeding rate to be determined by the swimming Péclet number, , defined as
(66) 
which measures the relative importance of advection of nutrients by the net displacement of the cell and diffusion. This is clearly the case at leading order for low Pe, as seen in (57).
In order to test the validity of this conjecture at higher Pe, we plot in figure 8 the mean feeding rate as a function of the “swimming Péclet number” , for the same large collection of unsteady strokes as in the previous section. All data points collapse rather well on a single curve, that corresponds exactly to the results for the steady treadmill swimmer (Michelin & Lauga, 2011). The agreement is particularly good for larger , corresponding to more efficient swimming strokes where the swimming motion dominates. The collapse of all the data points on that curve indicates that at leading order, for all strokes and all Pe, the mean feeding rate is determined by the mean swimming velocity.
Figure 8 shows however that a significant number of points do not follow that leading order trend and are located above the grey treadmill curve. Indeed, for swimming strokes with poor efficiency (including those with ), the contribution from the mean swimming velocity to the mean feeding process is no longer dominant and the influence of other squirming modes, or from timevariations of the swimming velocity, cannot be neglected, so remains strictly greater than one.
V Optimal unsteady feeding
The results presented in the previous sections and in Michelin & Lauga (2011) suggest that (i) swimming determines feeding, at least at leading order, and as a result (ii) optimal swimming and optimal feeding strokes are essentially identical. In this section, result (ii) is confirmed directly by performing an optimization of the swimming stroke maximizing the average nutrient uptake for a fixed energetic cost. The approach and methods presented below are based on the frameworks presented in Michelin & Lauga (2010a, 2011) and generalized here to the unsteady feeding problem for periodic Lagrangian strokes.
v.1 Adjoint optimization framework
The rescaled nutrient concentration satisfies the advectiondiffusion problem, (18)–(20), and the mean feeding rate, , is given by
(67) 
where is the outward normal unit vector. Considering a small perturbation, , in the velocity field, at leading order and for fixed Pe (or equivalently, fixed energetic cost) the resulting modification in mean feeding rate, , is obtained at leading order as
(68) 
where is the resulting linear perturbation in the nutrient concentration field satisfying
(69) 
with Dirichlet boundary conditions, , both on the surface of the swimmer and in the farfield. The last term in (69) guarantees that is constant and is obtained from and using (15) as
(70) 
From (69), the change in mean feeding rate for constant Pe can be computed as
(71) 
where