Rheotaxis of spheroidal squirmers in microchannel ﬂow: Interplay of shape, hydrodynamics, active stress, and thermal ﬂuctuations

Microswimmers exposed to microchannel ﬂows exhibit an intriguing coupling between propulsion, shape, hydrodynamics, and ﬂow which gives rise to distinct swimming behaviors. We employ a generic coarse-grained model of prolate spheroidal microswimmers, denoted as squirmers, exposed to channel ﬂow to shed light onto their transport properties. The embedding ﬂuid is implemented by the multiparticle collision dynamics approach (MPC), a particle-based mesoscale simulation method, which includes thermal ﬂuctuations. Speciﬁcally, the inﬂuence of swimmer shape—spherical vs spheroidal—, active stress—pusher, ciliate, puller—, and thermal ﬂuctuations on their rheotactic behavior is analyzed. The microswimmers accumulate at the conﬁning walls at very low ﬂow rates. With increasing ﬂow strength, squirmers are depleted from the walls, and at high ﬂow rates are also depleted from the channel center. The squirmers show pronounced cross-channel swimming between the conﬁning walls with mixed oscillating and rotational motions due to thermal ﬂuctuations. This strongly affects their rheotactic behavior. In particular, spherical pullers and ciliates swim upstream, whereas spherical pushers essentially swim downstream. The anisotropic shape of spheroidal squirmers enhances wall and center depletion and the alignment of the propulsion direction parallel to the ﬂow, which leads to preferred downstream swimming for all active stresses. This emphasizes the importance of swimmer shape and hydrodynamic wall interactions on the transport properties of a microswimmer such as Volvox and Opalina , for example.


I. INTRODUCTION
The swimming behavior of motile microorganisms is determined by the response to external fields such as gravity, chemical or thermal gradients, geometrical restrictions, or flow fields.The reaction to specific fields, e.g., chemical or thermal gradients, requires sensing [1][2][3] and adaptation, whereas responses to other fields depend upon mechanical forces and torques [1,3,4], e.g., external flows.A broad spectrum of microorganisms is habitually exposed to flows.Biological microswimmers, such as spermatozoa, swim in the oviduct, bacteria populate the human gastrointestinal tract, pathogens move along blood vessels [1,[5][6][7][8], and phytoplankton in oceans and lakes are exposed to turbulent fluid motion [9].From a technological perspective, microfluidic devices open a route to control microbial spreading with far-reaching impact in medical and technical applications [4,8,10,11].In addition, synthetic microswimmers are being designed with the goal to transport cargo, e.g., to serve as drug carriers in human blood vessels.In all these cases-specifically for the rational design of microfluidic devices and synthetic swimmers-a detailed understanding of the behavior of such microswimmers in channel flows is desirable or even indispensable.
These studies provide valuable insight into various aspects of swimming of biological microorganisms in the vicinity of walls and under flow.In general, the swimming behavior depends on the characteristics of the microswimmer itself, namely, its shape and generated flow field.Typically, self-propelled particles are hydrodynamically classified as pushers-where thrust is generated at the rear, e.g., E. coli-, pullers-with thrust generated in the front part of the cell, e.g., Chlamydomonas Reinhardtii-, or ciliates, e.g., Volvox and Opalina [36,37].Most studies on swimming behavior in microchannel flows have been performed for pushers, specifically, E. coli [5][6][7]12,13,16,17,21].Here their anisotropic shape and the chirality of their flagellar bundle is considered to be responsible for positive rheotaxis and migration in the vorticity direction of the flow (vorticity FIG. 1. Sketch of a spheroidal squirmer in a channel exposed to pressure-driven flow.migration) [6,13,16,33].In comparison, little is known about the nonequilibrium swimming behavior of pullers and ciliates, aside from some studies on phytoplankton cells [9] and on Dunaliella primolecta, a puller with internal rotation, which exhibits swimming along the vorticity direction in shear flow [38].
The interplay of flow, confinement, and propulsion mechanism is even more subtle for synthetic microswimmers.Phoretic microswimmers display rheotaxis and upstream swimming [39,40], and their rheotactic behavior can be controlled by hydrodynamic forces, chemical catalysis, and acoustic propulsion [40].Moreover, catalytic Janus particles in channel flow exhibit vorticity migration [41], underlining the complex behavior emerging due to confinement, chemical gradients, and external flow.
Insight into the intriguing properties of microswimmers in flows can be obtained by computer simulations [42].Here we apply a hybrid simulation approach, combining the multiparticle collision dynamics (MPC) method for the fluid, with a coarse-grained representation of the microswimmera squirmer [43][44][45][46][47][48][49][50][51][52][53]-to study the properties of elongated microswimmers in channel flow.MPC is a particle-based mesoscale hydrodynamics simulation approach for fluids [54][55][56][57][58], which can efficiently be implemented in simulations on graphics processing units (GPUs) [59].A squirmer is a spheroidal colloid with a prescribed slip velocity on its surface.It has originally been introduced to describe ciliated microswimmers such as Paramecia and Volvox [43,60].It is currently employed to studies of a broad class of microswimmers, both biological and synthetic.A squirmer is typically characterized by two modes, accounting for its swimming velocity and its active stress.The latter distinguishes between pushers, pullers, and ciliates (often denoted as neutral squirmers) [45,49,52,53,61].Fluid and squirmers are confined in a channel of two parallel no-slip walls and periodic boundary conditions in the other spatial dimensions (Fig. 1).An applied constant force on the fluid particles gives rise to a parabolic flow profile (Poiseuille flow) [62,63].We consider various fluid velocities, different active stresses, and two particle aspect ratios.
At low flow velocities, our simulations show significant wall accumulation of the squirmers independent of their active stress but with a stress-specific preferred alignment, as is well established [5,19,50,61,[64][65][66][67][68].Increasing flow leads to swimmer depletion at walls and at high shear rates additionally to depletion in the channel center, in agreement with previous observations in simulations, theory, and experiments [14,69,70].The squirmers continuously traverse the channel from one wall to the other and back, with the propulsion direction showing oscillatory and rotational motion, which is reminiscent of the trajectories of athermal self-propelled particles (SPP) [7].However, thermal fluctuations destroy the theoretically predicted stable trajectories, fixed points, and limit cycles.This severely changes the squirmers' rheotactic behavior, and spherical pushers predominantly swim downstream rather than upstream, in contrast to analytical predictions at low flow rates [7].In contrast, spherical ciliates and pullers exhibit upstream swimming next to walls for a wide range of flow rates.The anisotropic shape of spheroidal squirmers affects the rheotactic behavior, and squirmers swim predominately downstream, independent of the active stress.Hence, wall hydrodynamic interactions and shape determine the rheotactic and swimming behavior of squirmer-type microswimmers.
The paper is organized as follows.Section II briefly summarizes theoretical predictions for simple athermal selfpropelled point particles in channel flow.The microswimmer model and the simulation approach are introduced in Sec.III.Section IV presents the results for spherical squirmers and Sec.V those of spheroidal squirmers.The results are discussed and conclusions are provided in Sec.VI.

II. MICROSWIMMER IN POISEUILLE FLOW-THEORETICAL BACKGROUND
To illustrate certain basic features of microswimmers in Poiseuille flow, we briefly review the behavior of a confined, athermal SPP in a parabolic flow field in two dimensions [7].The translational equation of motion of the particle position r = r(cos φ, sin φ) T is where U 0 is the swimming velocity, e the propulsion direction (|e| = 1), and v f = 4u m (H/2 − y)(H/2 + y)e x /L 2 the unperturbed position-dependent fluid flow velocity in the x direction (|e x | = 1) with no-slip lines located at y = ±H/2.
The orientation e changes in response to the local vorticity, = ∇ × v f , as with = | | and e φ the unit vector along φ in polar coordinates.In particular, the equations of motion for the angle φ and the y coordinate are or, written as a second-order differential equation, This is the equation of a mathematical pendulum with a stable fixed point at φ = π .The integral of motion is the (normalized) energy Figure 2 illustrates the phase space (y-φ) for various values of the conserved quantity E , with its two parts-oscillating and rotational trajectories-separated by the separatrix (E = 2) and the unstable fixed point at X = (φ, y) = (0, 0), equivalently, X = (2π, 0).For E < 2, a SPP exhibits an oscillatory motion, denoted as swinging in Ref. [7], where 0 < φ < 2π .The covered range of angles and the displacement along y depends on the value of E .The value E = 0 corresponds to the stable fixed point X = (π, 0), with straight upstream motion in the center between the walls for u m /U 0 < 1 and downstream motion for u m /U 0 > 1.In case of E > 2, the rotational regime, φ covers the full range of possible values 0 < φ 2π ; this is often denoted as tumbling regime [7].More generally, a SPP swims upstream, ẋ/U 0 < 0, for u m /U 0 < 1 − E , and downstream for u m /U 0 > 1 − E .
The simple SPP model yields deterministic trajectories due to the absence of thermal fluctuations.Moreover, there is no wall accumulation due to lack of hydrodynamic and steric wall interactions, whereas wall accumulation is a hallmark of active matter and microswimmers [1,5,17,27,66,[71][72][73][74][75].Taking hydrodynamic wall interactions into account within the force-dipole approximation, the fixed point X is found to be unstable for pushers and stable for pullers [7].
Moreover, a confined slender-body swimmer without hydrodynamic wall interactions has been considered in presence of thermal fluctuations [19].The slender-particle rotation includes contributions from the imposed flow via Jeffery's equation [76] and rotational diffusion.The impenetrability of the walls is taken into account by the condition of a zero normal component of the translation flux at a wall [19].This approach yields wall accumulation, swimmer depletion in the channel center, and near-wall upstream swimming by shear alignment of the particles within the accumulation layer adjacent to a wall.This led to the conclusion that particle-wall hydrodynamic interactions are not necessary to explain these phenomena [69].
In addition, the rheotactic properties of athermal spherical diffusiophoretic Janus-type active particles have been analyzed [35].Accounting for the propulsion mechanism via an effective slip velocity, the transport of squirmers and catalytic Janus particles is considered, wherein the latter case, the solute (fuel) concentration is taken into account explicitly.Rheotactic behavior may emerge via self-trapping near a hard wall by auto-chemoattraction of the produced chemicals [35,77].Swimmer rotation and alignment by flow, combined with the solute concentration field, then leads to locking of the propulsion direction in the shear-gradient plane and upstream or downstream swimming at a steady height above the wall.Positive rheotaxis can be achieved by tailoring the surface chemistry of a catalytically active Janus particle [35].

III. MODELS FOR MICROSWIMMER AND FLUID A. Squirmer
A squirmer is modeled as a prolate spheroidal rigid body of mass M with the prescribed surface velocity [49,52], in terms of the spheroidal coordinates (ζ , τ, ϕ).They are related with the Cartesian coordinates according to The coefficient β accounts for the active stress [45,49], where β < 0 corresponds to a pusher, β > 0 to a puller, and β = 0 to a ciliate.A sphere of radius R follows for b x = b z ≡ R.
A more general representation of the flow field, including infinitely many squirming modes, is presented in Ref. [78].As a consequence, the swim velocity and active stress of a spheroidal squirmer depend not only on B 1 and B 2 , respectively, but include contributions from further modes.
In general, the propulsion direction e, expressed in spherical coordinates, is e = (cos φ sin θ, sin φ sin θ, cos θ ) T , (9) with the azimuthal angle 0 < φ 2π and the polar angle 0 θ π .The solution of the rigid-body equations of motion of a squirmer is described in Refs.[49,52], as well as in the Supplemental Material [79].

B. Multiparticle collision dynamics
In MPC, the fluid is represented by point particles of equal mass m undergoing subsequent streaming and collision steps [54][55][56].Flow is induced by a gravitational force acting on every fluid particle along the x axis of the Cartesian reference frame (Fig. 1) [56,80].In the streaming step, the fluid particles move ballistically transverse to the flow direction and experience a constant force along the flow direction during a time interval h, denoted as collision time [63].For the collision step, we employ the stochastic rotation dynamics (SRD) version of MPC with angular momentum conservation (MPC-SRD+a) [57,58,81], where particles are sorted into the cells of a cubic lattice of lattice constant a, and their relative velocities, with respect to the center-of-mass velocity of each cell, are rotated around a randomly oriented axis by a fixed angle α.For every cell, mass, momentum, and angular momentum are conserved.Since energy is not conserved in the collision step, we apply the Maxwell-Boltzmann scaling approach, a cell-level canonical thermostat at temperature T [82,83].The latter ensures Maxwellian distributed velocities.Partition of the system into collision cells leads to violation of Galilean invariance, which is reestablished by a random shift of the collision-cell lattice at every collision step [84,85].The MPC algorithm is highly parallel, hence we implement it on a graphics processing unit (GPU) for a high-performance gain [59].Further details are described in Ref. [79].

C. Wall interactions
No-slip boundary conditions at the walls are implemented for the MPC fluid by the bounce-back rule and phantom wall particles [62,63].The repulsive wall interactions of a squirmer are captured by a truncated Lennard-Jones potential for the closest distance of a point on the squirmer surface and the wall.This is straightforward for a sphere but somewhat more demanding for a spheroid.Details are presented in Refs.[49,52,79].The wall potential allows us to control the squirmer-wall distance.If not otherwise stated, we set the closest approach to approximately 1.4a, which is roughly 40% of the sphere radius and 20% of the spheroid major semiaxis.This insures sufficient fluid between squirmer and wall, and the emergence of a suitable flow field within the MPC approach.Simulation studies of passive spherical colloids indicate that for such close wall approaches lubrication forces are already significant, specifically for the dynamics normal to the wall [86].However, the transport parallel to the wall is determined by far-field wall hydrodynamic interactions [67].

D. Parameters
The mean number of MPC particles per collision cell is N c = 10-corresponding to the mass density ρ = m N c = 10m/a 3 -, the rotation angle α = 130 • , and the time step h = 0.05 ma 2 /(k B T ), which yields a fluid viscosity of η = 7.2 mk B T /a 4 ; k B is the Boltzmann constant and T is the temperature.
Spherical squirmers of radius R = 3.6a = σ /2, where σ is the diameter, and spheroidal squirmers with b z = 6a and b x = 3a (aspect ratio b z /b x = 2) are considered, which are neutrally buoyant, i.e., M = (4π/3)b 2 x b z ρ.We set B 1 = 0.05 √ k B T /m, which corresponds to the swimming velocity )/m for a sphere and U 0 = 0.04 √ (k B T )/m for a spheroid.The rotational diffusion coefficients of a sphere and a spheroid around the minor axis are D R = 1.2 × 10 −4 / mk B Ta 2 and D ⊥ R = 6.8 × 10 −5 / mk B Ta 2 , respectively.With the definition of the Péclet number Pe = U 0 /(σ D R ), we obtain Pe ≈ 38.5 for a sphere (σ = 2R) and Pe ≈ 49 for a spheroid (σ = 2b z ).These values are in the range of synthetic [4] and biological [12] microswimmers.The Reynolds number is Re = ρU 0 R/η = 0.16 and 0.21 for a sphere and spheroid, respectively, when we approximate the spheroid by a sphere with the same volume.
The distance between the confining walls is H = 5σ .Parallel to the walls, we set L x = 5H and L z = H.Hence, the ratio of persistence length in the bulk, l p = U 0 /(2D R ), and the wall separation is l p /H = Pe/10 > 1.The strength of the flow field can be characterized by the Weissenberg number Wi = γ /(2D R ), where γ = gρH/(2η) = 4u m /H is the shear rate at a wall, with mg the applied (gravitational) driving force and u m the maximum flow velocity in the channel center.In terms of the Péclet number and for the applied geometry, the Weissenberg number can be expressed as Wi = (u m /U 0 )(2Pe/5).We characterize the flow strength by the ratio u m /U 0 in the following.
Dilute systems of squirmers are considered with a packing fraction of approximately 0.03.Flow simulations of passive spherical and spheroidal colloids at this density yield a parabolic flow profile of the form U x /u m = 4y(H − y)/H 2 , with the maximum velocity u m /u f m = 0.9, where u f m is the theoretical maximum velocity of the MPC fluid.The value u m is somewhat smaller than the value of the MPC fluid, which is attributed to an increased viscosity partly due to the finite colloid density.In the following, velocities are normalized with respect to the maximum flow velocity u m of the passive colloids in the flow.Explicitly, we consider the ratios u m /U 0 = 0.28, 0.57, 1.14, 2.25, 4.5, 6.75, and 11.3.This covers the range 4 Wi 190 for spheres and 5 Wi 250 for spheroids.

IV. RESULTS FOR SPHERICAL SQUIRMERS
The rheotactic properties of spherical squirmers adjacent to walls and under flow have been studied theoretically and numerically [7,35,50,66], typically without thermal fluctuations.Here we demonstrate that such fluctuations substantially affect their swimming behavior in Poiseuille flow.9)] for various flow rates and active stresses, corresponding to the phase-space plots of Fig. 2. The flow lines indicate the (average) swimming behavior, which would be followed by a squirmer in absence of thermal fluctuations and are reminiscent of the deterministic flow lines in Fig. 2. Evidently, the distribution function P(y, φ) is rather inhomogeneous and no deterministic trajectories exist anymore as a consequence of thermal fluctuations.Instead, distinct areas of high probability are present which depend on the applied flow velocity.
At low flow rates, u m /U 0 0.3, P(y, φ) exhibits strong maxima adjacent to the walls with preferred angles.This accumulation implies squirmer depletion from the channel center, but the finite bulk density indicates a pronounced cross-channel migration of squirmers.An example trajectory is displayed in Fig. 4(a).For the low flow rate, the squirmer is mainly aligned parallel to the flow direction.Starting on the upper wall, the squirmer swims upstream toward the lower wall with a propulsion orientation slightly larger than π .At the lower wall, the orientation changes to a value slightly smaller than π , while the squirmer moves along the wall, and it swims toward the upper wall again (cf.movie M4 [79]).This corresponds to the oscillatory motion as described in Sec.II.However, fluctuations perturb the predicted deterministic motion and other features appear, such as local reorientation near the walls, as reflected in Fig. 4(a), e.g., at time t/h = 1.5 × 10 3 and yσ ≈ 1.6.This dynamics is typical and independent of β.
The phase-space plots of Figs.3(d)-3(f) for the larger flow velocity u m /U 0 = 2.25 closely resemble the structure of the SPP phase space of Fig. 2, since the probability for positions off the walls increases.Moreover, the range of angles φ broadens.As for the smaller flow rate, the squirmer oscillate (swing) between the walls with high probability for angles in the interval π/2 φ 3π/2.Athermal theoretical models, i.e., Pe = ∞, taking only far-field hydrodynamics into account, predict a fixed point at X = (φ, y/σ ) = (π, 2.5), where the fixed point is stable for pullers and unstable for pushers [7,87].For ciliates, X is marginally stable (center) [88].This is reflected in Figs.3(d)-3(f).P(y, φ) of the pusher is small at X , but a "band" of high probability exists reminiscent of the limit cycle of an athermal system [7].In the case of pullers, the trajectories inside of the (athermal) limit cycle are attracted toward X [Fig.3(f)].Ciliates still exhibit a pronounced "band" of high probability in the vicinity of the limit cycle.Noteworthy, the probability at the walls is still very prominent.In fact, there are (athermal) fixed points close to a wall with preferred orientations [50].In the presence of noise, these points of highest probability correspond to maxima in the density and the angular distribution function (Fig. 5).
At even higher flow rates, u m /U 0 > 6.7, rotational trajectories appear, with a high probability for trajectories covering the whole range of angles 0 As illustrated in Fig. 4(b), pushers swim downstream for u m /U 0 = 6.75, with a large variation of the propulsion angle.Rotational parts appear typically close to walls and in the center (cf.movie M10 [79]).Rotational trajectories are frequently interrupted by oscillations, and vice versa, where for oscillating parts squirmers are able to traverse the channel persistently from one wall and back, and the azimuthal angle varies between 0 < φ < 2π only.
Similar distribution functions are obtained for β = ±1, as shown in Fig. S1 [79].There is a somewhat larger probability P(y, φ) in the channel center for low flow velocities u m /U 0 < 0.3 and a smaller one for higher flow rates compared to the cases β = ±3.
The differences in the trajectories-oscillating and rotational motion-at the various flow rates are reflected in the autocorrelation function e(t ) • e(0) of the propulsion direction and its component along the flow direction (Fig. S2 [79]).The oscillatory motion for flow velocities u m /U 0 < 5 implies a power-law, or near power-law, decay of e(t ) • e(0) as a function of time, with an active-stress-dependent exponent.Rotation (tumbling) is revealed by oscillations, where the correlation e x (t )e x (0) is both positive and negative at short times for u m /U 0 > 6.For the lower flow rates, the correlation functions decay exponentially on average at longer times.
Evidently, hydrodynamic fluctuation perturbs the pure oscillatory or rotational scenario.The probability densities in Fig. 3 indicate the most probable values for the squirmer position and orientation.However, an individual squirmer can assume any angle and position.Moreover, there are no more strict stable or unstable fixed points or limit cycles.Yet, the high probabilities adjacent to walls in Fig. 3, specifically for u m /U 0 = 0.28, are remnants of a stable fixed point at the wall for nonfluctuating systems [18,50].

B. Density distribution
A generic effect of active particles is their accumulation at walls, which leads to an increased wall density [5,8,19,22,50,[64][65][66][67]74,75,89] and is also apparent in Figs.5(a)-5(c) for our systems (see also Figs.S3(a)-S3(c) [79]).However, wall accumulation depends significantly on the applied flow strength.For low active stresses and weak flows, the wall density is rather similar for pushers, ciliates, and pullers.The squirmers swim toward the wall, and only after rotation by slow diffusional processes are they able to leave again.Accumulation reduces with increasing flow, which implies an increasing number of squirmers in the channel center.A similar drop in wall density under flow has been observed for a dumbbell-type microswimmer [70,90] and slender-body active particles [69], and is attributed to a preferred parallel alignment of the propulsion direction by shear flow, which implies fewer wall impacts and a faster detachment.The wall depletion is rather similar for ciliates and pullers but is less pronounced for pushers.Pushers (β < 0) prefer to stay at a wall, and significantly higher flow rates are required to achieve a comparable depletion as for β 0. This can be explained by preferences in the distribution of propulsion directions adjacent to a wall (cf.discussion in Sec.IV C).
For large flow velocities u m /U 0 10, a maximum of the squirmer density appears roughly in the middle between the wall and the channel center.This is reminiscent of flowinduced concentration dips in the channel center observed in experiments [14] and simulations [19,70,90].The maximum is due to a shear-induced flow trapping of the squirmers, which depends on the strength of the active stress and is most pronounced for β 0. Qualitatively, we explain the effect by an increased residence time of the squirmer in the sheartrapping region, where a density peak emerges and shifts toward the channel center with increasing flow strength (cf.Sec.VI for more details).

C. Sphere azimuthal orientation (Flow-gradient plane)
Microswimmers are known to exhibit a particular orientation at a wall due to hydrodynamic interactions, which depends on their active stress: pushers preferentially align parallel and pullers perpendicular to a wall [1,5,12,16,18,65,67].This behavior is reflected in the probability distribution function P w (φ) of the azimuthal angle in the layer 0 < y/σ < 1 next to the wall at y = 0 displayed in Figs.5(d)-5(f) at weak flows [see also Figs.S3(d)-S3(f)].The distribution function P w (φ) exhibits pronounced maxima at φ = 0 (2π ) and π for pushers, corresponding to down-and upstream orientation, respectively.For ciliates, an upstream maximum in the vicinity of φ π is present [Fig.5(e)], and the maximum for pullers at φ 3π/2 indicates alignment toward the wall [Fig.5(f)].With increasing flow velocity, the distribution P w of pushers develops a maximum for angles π/2 < φ π , i.e., the preferred orientation is upstream and away from the wall.The maximum of the distribution function for ciliates and pullers is in the range π φ < 3π/2; hence they are preferentially aligned upstream and toward the wall.At high flow rates, P w (φ) becomes broad and rather symmetric with a maximum at approximately φ = π , corresponding to parallel alignment with the flow.
The results in terms of squirmer alignment adjacent to a wall are qualitatively consistent with theoretical and simulation studies [50,67].These studies predict stable fixed points or oscillations next to walls with preferred angles, where the orientation for pushers and pullers points away and toward the wall, respectively.However, the theoretical considerations predict downstream swimming for pushers and upstream swimming for pullers in the stationary state [50].Thus noise changes the swimmer behavior qualitatively.
The alignment in the channel center 2 y/σ 3 is rather different but is also influenced by confinement, since the scaled persistence length l p /H for the swimming motion is larger than unity, as is reflected by the probability distribution function P c (φ) in Figs.6(a) and S6 (see also Fig. S3 for β = ±1 [79]).All squirmer types exhibit a preferred upstream orientation at low flow rates, u m /U 0 < 1, with a rather broad and flat maximum.However, for u m /U 0 5, the distribution function shows two maxima symmetric with respect to φ = π which are shifted to smaller and larger angles, respectively, with increasing flow rate.Naturally, the latter is related to the symmetry of the setup but dynamically appears as a consequence of the squirmers' movement from the wall at y = 0 to that at y = H, with a preferred orientation of the squirmers, and vice versa.Figure 6(b) depicts the dependence of the maximum φ m of P c (φ) for φ < π on the flow velocity for various β.The continuous broadening of the maxima at φ = π implies an abrupt emergence of peaks at φ m < π for all active stresses, reminiscent to a pitchfork bifurcation.For larger flow velocities, φ m decreases nearly linearly.Interestingly, the data for β = 0, ±1 are rather similar, suggesting a weak dependence on active stress only.The values for |β| = 3 deviate slightly, specifically, for β = 3 the transition appears at large u m /U 0 only (Fig. S6(c) [79]).The "critical" flow velocity for the transition depends on the magnitude and sign of the active stress.For pullers, larger β values imply a shift of the critical flow velocity to larger values, whereas for pushers smaller β lead to a splitting at small flow strengths.
the polar angle θ [Eq.( 9)] in the accumulation layer (y/σ < 1) reveals a substantial migration in the vorticity direction and a strong dependence on the active stress (Fig. S7).Specifically, pushers exhibit strong migration nearly along the vorticity direction (θ = 0, π), with a weak dependence on flow strength only.Ciliates show (weak) vorticity migration for low flow strengths, whereas a uniform distribution is obtained at flow velocities u m /U 0 > 1.2.In qualitative difference, pullers preferentially align along the flow direction (θ = π/2), with a distribution function only weakly depending on flow velocities.For the smaller active stresses |β| = 1, a similar behavior is obtained, however, with lower alignment probabilities and a stronger preference for a uniform distribution, as shown in Fig. S4.
In the channel center 2 < y/σ < 3, the distribution of polar angles is rather uniform for all active stresses and flow velocities.But for weak flows u m /U 0 1, the squirmers exhibit a small preference for swimming in the vorticity direction (Figs.S4 and S7).
Swimming in the vorticity direction of synthetic spherical Janus particles has been observed experimentally [41], and the preferred orientation is explained by a balance of contributions to the angular velocity of the particle from shear flow, bottom heaviness, and swimming near a planar substrate.In contrast, migration in the vorticity direction of our squirmer pushers and ciliates is solely determined by hydrodynamics, shear flow, and active stress.

E. Local squirmer velocity
Figure 8 illustrates the impact of flow on the local velocity U x of squirmers along the flow direction.For all β and u m /U 0 < 5, the streaming velocities are smaller than those of passive colloids.At high flow strengths, u m /U 0 6.75, the velocity profile approaches that of passive colloids in the channel center, the better the larger β.However, in the vicinity of the wall the squirmer velocity is always smaller than that of the passive colloid.This difference in the active and passive profiles indicates the presence of a large fraction of squirmers swimming upstream, particularly at low flow strengths and large β.
Steric and hydrodynamic wall interactions affect the squirmer velocity, especially for low flow strengths.In the accumulation layer adjacent to the wall, the squirmer velocity is high as a consequence of the strong alignment of the propulsion direction [Figs.5(d)-5(f)].The strongly varying velocity curves for the various active stresses reflect the significant interference of the squirmer flow field with the no-slip wall.In particular, negative velocities appear for β 0 as long as u m /U 0 0.6, which indicates upstream swimming (positive rheotaxis).In contrast, the velocities U x for pushers with β = −3 are always positive, with a weak u m /U 0 dependence for u m /U 0 0.65.This is in qualitative agreement with previous squirmer simulations [50].
The rheotactic behavior of pushers depends on the magnitude of the active stress, and they exhibit (weak) positive rheotaxis next to a wall for β = −1, whereas the velocity profile of pullers is far less affected, as shown in Fig. S5.The average squirmer velocity in the layer y/σ < 1.5 next to a wall is presented in Fig. 9, which emphasizes the significant upstream velocity of ciliates and pullers for all β.In comparison, the upstream velocity of pushers is rather weak, even for β = −1.A weaker upstream velocity for pushers also has to be found for phoretic rods [91].
Downstream swimming of pushers is in stark contrast to the behavior of bacteria, such as E. coli, or sperm, which exhibit positive rheotaxis while swimming at walls and channels [6,13,16,33,92].Evidently, the observed positive rheotaxis of flagella-propelled bacteria cannot simply be explained by the propulsion mechanism (pusher).The obtained features may be related to the shape of the microswimmer, as bacteria are typically rather elongated than spherical; however, a spheroidal pusher squirmer does not show positive rheotaxis either, as discussed in Sec.VI.

V. RESULTS FOR SPHEROIDAL SQUIRMERS
The shape of microswimmers plays a major role while they interact sterically or hydrodynamically with confining walls [18,52,93].To elucidate the influence of asphericity on the rheotactic behavior, we consider spheroidal squirmers with the aspect ratio b z /b x = 2.

A. Phase space: Joint probability distribution function
Figure 10 displays joint probability distribution functions P(y, φ) for spheroidal pushers with β = −1.The distribution functions qualitatively differ from those of spherical squirmers, Fig. 3, in particular, the probability for squirmers adjacent to walls and the preference for parallel alignment with respect to the walls is higher.Moreover, the distribution functions are rather similar for pushers, ciliates, and pullers, as shown in Fig. S8; the active-stress differences are less evident compared to those for spherical squirmers.Nevertheless, the characteristics of the trajectories are similar.As indicated by the flow lines and is illustrated in Fig. 11, the squirmers steadily migrate from one wall to the other while changing their orientation continuously.
At the low flow strength u m /U 0 = 0.28, both downstream and upstream sequences appear, and the angle φ varies mainly in narrow ranges in the vicinity of φ = 0, π, and 2π , respectively.Hence, steric wall and shear alignment imply  oscillation (swinging) motion of the athermal swimmer of Sec.II [7], as depicted in Fig. 11.However, fluctuations imply a transition between upstream and downstream swimming sequences, and no stable oscillatory motion appears.At high flow rates, u m /U 0 = 6.75, the angle φ changes by 2π while the squirmer traverses several times the space between the walls.Here, oscillations are interrupted by rotational motions, and vice versa, in a random sequence.As for spherical squirmers, fluctuations destroy the dynamical patterns predicted by far-field and athermal theories, and squirmers typically cover the whole distance between the walls rather than being confined in certain spatial areas.In any case, there is no stable fixed point at X = [π, H/(2σ )], and limit cycles are unstable.However, there are (athermal) fixed points adjacent to the walls, corresponding to the maxima in the density, Fig. 12(a), and angular distribution functions, Fig. 12(b) (see also Fig. S10).
Active stress substantially affects P(y, φ) for pushers, as indicated in Fig. S12 for β = −3.Specifically, the probability for wall accumulation and parallel orientation of the major axis with the wall increases (Fig. S13).For β = 3, the squirmers self-organize into flowerlike structures at a wall; hence we refrain from a discussion of their flow behavior here.
As for spherical squirmers, the propulsion autocorrelation function e(t ) • e(0) of Fig. S9 decays exponentially or by a power law for the respective flow rate.Again, oscillations of the various curves indicate the tumbling motion of the squirmers when u m /U 0 > 2.2.

B. Density distribution
The density distribution function, displayed in Figs. 12 and  S10, shows approximately equally large density maxima at y/σ ≈ 0.37, the closest approach along the minor squirmer axis, for pushers and ciliates, but somewhat higher values for pullers at vanishing flow rates.The density decreases with increasing flow strength and ρ increases in the channel center.For flow velocities u m /U 0 5, squirmers are depleted next to a wall, similar to spheres; however, depletion of spheroids is more pronounced, which is related to the rotation of the extended spheroid next to a wall.As for spherical squirmers, a peak in the density distribution appears in the shear-trapping region at y/σ ≈ 1.3 for large flow velocities u m /U 0 > 6.Note that the density in the channel center is nearly equal to the bulk density.Again, we explain the maximum by an enhanced residence time of the microswimmer in the shear-trapping region (cf.discussion in Sec.VI).The maximum moves toward the channel center with increasing flow strength.The shape of the minimum in the channel center is well fitted by a parabola with a curvature, which increases with increasing flow velocity, consistent with previous simulation and theoretical studies [90].Remarkably, the differences between pushers, ciliates, and pullers is less pronounced compared to spherical squirmers.

C. Spheroid azimuthal orientation (Flow-gradient plane)
Figures 12(b) and 12(c) show the distribution function of the azimuthal angle φ for ciliates in the vicinity of a wall, P w , and in the channel center, P c , for various flow rates.We display results for ciliates only, because the behavior for pushers and pullers is very similar, as shown in Fig. S10, in contrast to the results of Fig. 5 for spheres, which show qualitative differences for the various β values.Adjacent to a wall [Fig.12(b)], steric and hydrodynamic interactions lead to a preferred parallel wall alignment, independent of the active stress and the strength of the flow.However, pushers point preferentially slightly away from the wall, with maxima of the distribution function at φ 0 and φ π , while pullers point toward the wall with maxima at φ 2π and φ π (Fig. S10).With increasing flow strength, the maximum at φ ≈ 0 gradually vanishes and only that at φ ≈ π remains.Hence in strong flows the spheroidal squirmers point upstream adjacent to walls, but, on average, move with the flow downstream.The situation is rather different in the channel center.For flow strengths u m /U 0 1, upstream (φ ≈ π ) and downstream (φ ≈ 0) orientations occur with approximately equal probability.With increasing u m /U 0 5, downstream orientations become dominant.The upstream-oriented squirmers are rotated by the flow in the downstream direction when they leave the accumulation layer.This is the rotational motion caused by the flow (cf.movie M23 [79]).
Consistent with the distribution P(y, φ) (Fig. S13), the probability P w (φ) for downstream orientation of pushers with β = −3 (Fig. S13) is significantly higher than that for β = −1 for all considered flow rates (Fig. S10).This emphasizes the relevance of active stress for microswimmer migration at walls.

D. Spheroid polar orientation (Flow-vorticity plane)
Spheroidal squirmers exhibit substantial vorticity migration adjacent to a wall at smaller flow velocities, which is qualitatively similar for all considered active stresses but is most pronounced for pushers (Figs. 13, S11, S14).Evidently, spheroidal pullers behave differently from spherical pullers, which we attribute to specific squirmer-wall hydrodynamic and steric interactions by the spheroidal shape.With increasing u m /U 0 , the distribution functions for ciliates and pushers flatten and develop a small maxima at θ = π/2, suggesting a preferred alignment along the flow direction.Combined with the high transport velocity, such squirmers quickly traverse the channel from one wall to the other, hardly leaving time for a cross-stream (vorticity) orientation next to walls.Pushers preferentially swim in the flow-vorticity direction even at strong flows, only for u m /U 0 = 11.3 and β = −1, nearly equal probabilities for all angles are obtained.However, for β = −3, swimming in the vorticity direction is even more pronounced at high flow strengths (Fig. S14).Here the spheroidal squirmers exhibit a rotation along the major axis, induced by the local shear flow.
In the channel center, the distribution of polar angles depends only weakly on the flow strength for all active stresses.Yet for weak flows u m /U 0 < 2.3, the squirmers exhibit a small preference for swimming in the vorticity direction which changes to a preferred alignment along the flow direction for larger flow velocities (Figs. 13, S11, and S14).

E. Velocity distribution
Simulations of passive spheroidal colloids yield the same parabolic flow profile as for spheres, i.e., the profile is independent of colloid shape.Velocity profiles for spheroidal squirmers are displayed in Fig. 14 for β = 0, ±1.The impact of flow is less pronounced compared to that of spherical squirmers [Fig.8(d)-8(f)].In particular, there are essentially only positive U x values, and the profiles U x (y) vary far less with flow strength.This applies specifically to pushers for β = −3 with their nearly flow-strength-independent profiles (Fig. S15).Remarkably, there is essentially no positive rheotaxis for spheroidal squirmers, in contrast to analytical calculations in absence of fluctuations [94].
As for spheres, the velocities U x of squirmers are always smaller than those of passive colloids for flow velocities below a "critical" value, which depends on the active stress.For example, U x /u m < 0.9 for u m /U 0 < 2.4 and β = −1, whereas u m /U 0 < 6 and β = 1.Again, this reflects a large fraction of upstream swimming squirmers.Interestingly, for large u m /U 0 , in contrast to spherical squirmers, U x /u m can exceed the value of passive colloids in the channel center due to a preferred downstream orientation of the spheroids.

VI. DISCUSSION AND CONCLUSIONS
Flow, wall hydrodynamic interactions, active stress, thermal fluctuations, and shape determine the transport of squirmers in channels.Various aspects of these have been addressed previously, but not all of them, and their interplay, have been taken into account simultaneously before.
Our squirmer simulations exhibit various qualitative features of microswimmers at walls and in microchannel flow which are consistent with theoretical and experimental observations for Janus particles, bacteria, and spherical ABPs.At weak flows, spherical as well as spheroidal squirmers accumulate at walls, and specifically, spheres exhibit a preferred active-stress-dependent orientation-pushers parallel and pullers perpendicular to walls.An increasing flow rate leads to squirmer depletion at walls, and for high flow rates, a density maximum appears in a shear-trapping region between the wall and channel center.This is reminiscent of a density dip in the channel center found before in experiments on bacteria [14] and simulations of active dumbbells and rods [69,70,90].Wall depletion is attributed to flow-induced changes in the distribution of the azimuthal angle, from a preference of the propulsion direction pointing toward a wall to pointing away from it.In contrast to previous interpretations of the center density dip in terms of a fast crossing of the channel center by the microswimmer [70], we explain the increased squirmer density by an increased residence time in the shear-trapping region, because there is no obvious reason why squirmers well aligned with the flow direction should exhibit an increased transport velocity in the flow-gradient direction.Shear flow leads to squirmer tumbling, which appears above a critical shear rate in the trapping region in the vicinity of a wall and moves with increasing flow velocity closer to the channel center.Specifically, spheroidal squirmers approaching a wall are preferentially oriented upstream by the flow.While leaving the wall, rotation by shear flow leads to a preferred downstream orientation in the channel center.This tumbling motion extends the relative time spent by a squirmer in the trapping region.In fact, several rotations are possible.Hence, the dynamics corresponds to the rotational motion of an athermal SPP of Sec.II at large flow rates.
Although the depletion effect and density modulation can all be qualitatively explained by steric interactions and flow only [69,70,90], hydrodynamic interactions, most prominently the active stress, imply qualitative and quantitative differences.These are reflected, e.g., in the phase-space plots of Figs. 3 and 10 and the velocity profiles of Figs. 8 and 14.In particular, the effect of thermal fluctuations severely affects the squirmer transport, specifically (i) fixed-point and limiting cycles of athermal systems [7,35] are blurred, and (ii) for any activity and active stress trajectories exhibit random combinations of oscillatory and rotational parts.Both aspects allow squirmers to explore the whole channel volume, although with flow-velocity and active-stress-dependent probability.This is reflected in the phase-space distribution function P(y, φ), which shows areas of preferred positions and azimuthal orientations of the propulsion direction.In particular, extremal values of P(y, φ) indicate surface accumulation, and an enhanced transport normal to the walls appears with increasing flow rate.Fixed points and limit cycles of athermal systems are replaced by extremal values of P(y, φ).The squirmer trajectories reflect the dynamical behavior and exhibit oscillatory and rotational features with a random switch between them, which is indicated by the flow lines of Fig. 3 and 10.
Rheotaxis is determined by the preferred orientation of the propulsion direction.Our simulations of spherical pushers exhibit weak net upstream swimming only at weak active stresses, β = −1, while downstream swimming occurs for β = −3.However, ciliates and pullers exhibit strong positive rheotaxis for β = 1 and β = 3 as long as u m /U 0 1.The hydrodynamic wall interactions of ciliates and even more of pullers imply an upstream orientation of the propulsion direction inclined toward a wall, which favors upstream swimming.In contrast, pushers preferentially point away from a wall and upstream swimming is absent.
Theoretical and simulation studies of athermal spherical squirmers predict a stable fixed point next to a wall for pullers (β > 0) and upstream swimming [35].Positive rheotaxis is predicted for the wall shear rate γ R/U 0 < 0.003 at β = 3, corresponding to u m /U 0 < 0.0015H/σ = 0.075 for our setup.We find positive rheotaxis for a much wider range of flow velocities at β = 1 and 3, specifically at walls.Moreover, theoretical calculations predict "locking" of the polar angle θ in the shear-gradient plane.As shown in Fig. S7, pullers exhibit a preference for such an alignment, but the distribution of the angle θ can be broad, in particular for small β (Fig. S4).Thus, confinement of the propulsion direction in the shear-gradient plane is not a prerequisite for upstream swimming.Remarkably, pushers display a pronounced vorticity migration.Whether this influences the stability of upstream swimming needs to be studied in more detail.
For spheroidal squirmers, the most noticeable difference to spherical squirmers is their distinct alignment at a wall, with a maximum in the probability distribution function of the azimuthal angle close to φ = 0, 2π and φ = π .This is due to steric wall interactions on the one hand and shear flow on the other.Hence the orientation toward a wall is far less pronounced compared to spheres due to the elongated shape.As spheroids leave a wall, they are rotated by the shear flow in the shear-trapping region and the propulsion direction points preferentially downstream; as a consequence, downstream swimming is favored, independent of the active stress.The observed selection of preferred downstream orientation of spheroids in the channel center seems to be a hallmark of elongated microswimmers in the presence of thermal fluctuations.Studies of dumbbell microswimmers suggest that the effect is dominated by flow and is independent of hydrodynamic interactions [70].
We observe a very strong vorticity migration at low flow velocities for any active stress, which weakly depends on the squirmer shape.Again, the effect on rheotaxis needs to be analyzed analytically.Swimming in the vorticity direction of catalytic spherical Janus particles has been observed experimentally [41], and the preferred orientation is explained by a balance of contributions to the angular velocity of the particle from shear flow, bottom heaviness, and swimming near a planar substrate.In contrast, migration in the vorticity direction of our squirmer pullers and ciliates is solely determined by hydrodynamics, shear flow, and active stress.
Experiments and theoretical calculations of gold-platinum phoretic rods find a pronounced upstream swimming for ciliates and pullers [91], in contrast to our spheroidal squirmer results.This discrepancy with the squirmer results could be related to additional effects by the concentration profiles of the reactant and products.The effect of the latter has been shown to contribute to upstream alignment of spherical Janus particles [35].However, simulation studies of (A-B-A)-type triblock rods with blocks of different slip velocities suggest that particular features of the hydrodynamic flow field suffice to create a preferred alignment of the propulsion direction toward a wall and hence to achieve positive rheotaxis [91].
Our approach predicts downstream swimming of spheroidal squirmers, independent of the active stress.This is in contrast to experimental results on flagellated bacteria, which swim upstream.Evidently, a pusher-type active stress is insufficient to describe upstream swimming of bacteria in channel flows.This is in line with extensive experimental and simulation studies which suggest a strong influence of the bacteria architecture and hydrodynamics on their rheotactic behavior, specifically, the presence of a helical flagellar bundle [6,13,16,33].Bacteria swimming in circles at walls are oriented toward the vorticity direction by shear flow due to the chirality of their flagella bundle and thus exhibit positive rheotaxis.For even larger flow rates, cells in the bulk swim in the vorticity direction with an upstream preference.Wobbling of the whole cell (Jeffrey orbits) enhances the effect [16].
Phenomenological nonhydrodynamic models neglecting, e.g., active stresses and (partly) wall steric interactions qualitatively describe some of the experimentally observed behaviors, such as positive rheotaxis or formation of a density dip in the channel center at sufficiently high flow strengths [19,90].Evidently, such models capture part of the phenomenology but provide only limited insight into the detailed underlying physical mechanisms.Nevertheless, such models can be rather useful, since they are able to predict various aspects of microswimmers in flow but may miss others completely.
We focused on a single Péclet number (Pe ≈ 50) that is representative for catalytic Janus particles [4] and tumbling bacteria [12] and varied the flow rate.Further studies are desirable to resolve the influence of the Péclet number on the rheotactic behavior of squirmers.At higher densities, microswimmers exhibit an intriguing collective behavior at walls and in thin films.Here, squirmer simulations can provide valuable insight into collective effects and their influence on rheotaxis, such as modifications of viscosity [95][96][97].Our simulations of spheroidal pushers with the active stress β = 3 show hydrodynamically organized flower-type structures at walls.The puller flow field of the spheroids, inclined with respect to the wall, leads to attraction and cluster formation for β = 3 but not for β = 1.This once more emphasizes the significance of active stress for the properties of microswimmers adjacent to walls.

FIG. 2 .
FIG. 2. Phase space (y-φ) of a confined self-propelled particle (SPP) in Poiseuille flow for (a) u m /U 0 = 5/4 and (b) u m /U 0 = 5.The color code designates the various values of E of Eq. (5).White lines are "flow" lines with arrows indicating the flow direction.The black line is the separatrix, separating rotational from oscillating trajectories.
and 0 ϕ 2π ; b x and b z denote the length of the minor and major semiaxis of the prolate spheroid.All points with τ = τ 0 = b z /c lie on the spheroid's surface.The propulsion velocity U 0 is related with the coefficient B 1 via[49]

FIG. 4 .
FIG. 4. Trajectory of a spherical squirmer (pusher, β = −3) for (a) u m /U 0 = 0.28 and (b) u m /U 0 = 6.75.The color code of the trajectory represents time and the arrows indicate the propulsion direction with the color code of the color wheel for the azimuthal angle φ.See also movies M4 and M10 for illustration.

FIG. 6 .
FIG. 6.(a) Probability distribution function P c (φ) of the azimuthal angle φ [Eq.(9)] for spherical squirmers (ciliates, β = 0) in the channel center 2 y/σ 3. The color code is the same as in Fig. 5. Results for β = ±3 are presented in Fig. S6 and for β = ±1 in Fig. S3.(b) Dependence of the maximum φ m (φ m < π) of the distribution function P(φ) in (a) on the flow strength u m /U 0 for the displayed active stresses.The vertical lines indicate the abrupt transition in φ m .

FIG. 8 .
FIG. 8. Velocity U x /u m along the flow direction of spherical squirmers as a function of their position between the walls.The flow velocities are u m /U 0 = 0.28, 0.57, 1.14, 2.25, 4.5, 6.75, and 11.3 (bottom to top), and the active stresses are (a) β = −3 (pusher), (b) β = 0 (ciliate), and (c) β = 3 (puller).The upper yellow curve is the flow profile of the respective system of passive colloids.The gray area indicates upstream velocities.

FIG. 9 .
FIG. 9. Average squirmer velocity in the layer 0 y/σ < 1.5 near a wall as a function of the flow velocity u m /U 0 for the active stresses β = 0, ±1, and ±3.

FIG. 10 .
FIG. 10.Probability distribution function P(y, φ) of the position y of spheroidal pushers (β = −1) within the channel normal to the walls and the azimuthal angle φ.The flow velocities are (a) u m /U 0 = 0.28, (b) 2.25, and (c) 6.75.The white lines are "flow" lines with arrows indicating the flow direction.The distribution function is presented logarithmically.For an illustration of the squirmer dynamics, see movies M13-M24 [79].

FIG. 12 .
FIG. 12. (a) Density distribution normal to a wall for spheroidal squirmers (ciliates, β=0) and various flow velocities.The bulk density ρ bulk is the number density of squirmers in the channel.(b) Probability distribution function P w (φ) of the azimuthal angle φ [Eq.(9)] in the layer 0 < y/σ < 1. (c) Probability distribution function P c (φ) in the channel center 2 y/σ 3. The flow velocities are u m /U 0 = 0.28, 1.14, 2.25, 6.75, and 11.3 [bottom to top at y/σ = 1.5 in (a)].The distribution function is normalized such that the integral of P(φ) is unity where P(φ) is the distribution function for the whole channel.

FIG. 14 .
FIG. 14. Velocity U x /u m along the flow direction of spheroidal squirmers as a function of their position between the walls.The flow velocities are u m /U 0 = 0.28, 1.14, 2.25, 6.75, and 11.3 [bottom to top at y/σ = 0.5 in (c)] and the active stresses (a) β = −1 (pusher), (b) β = 0 (ciliate), and (c) β = 1 (puller).The upper yellow curve is the flow profile of the respective system of passive spheroidal colloids.The gray area indicates upstream velocities.