Dynamics of squirmers in explicitly modeled polymeric fluids

Biological microswimmers such as bacteria and sperm cells often encounter complex biological fluid environments. Here we use the well-known squirmer microswimmer model to show the importance of the local fluid microstructure and non-continuum effects on their swimming speed in different polymeric and filamentous fluids. Surprisingly, we find that different squirmer types move at considerably different speed in filamentous fluids which cannot be explained by existing continuum models, but by considering the local fluid and polymer properties around the squirmers. Furthermore, direct squirmer-polymer interactions slow down in particular pushers by trapping large stiff filaments in a self-generated recirculation region in front of them.

Introduction.-Many microorganisms swim in viscous fluids, and their dynamics can often be described by low Reynolds number hydrodynamics [1,2].Furthermore, active colloids and droplets have been investigated extensivly as well-controllable model systems to study single and collective microswimmer behavior [3,4].A common minimal model to capture the hydrodynamics of biological and artificial microswimmers is the so-called squirmer [5][6][7].By applying different surface velocity fields to a rigid sphere, different types of microswimmers can be modeled: pushers such as bacteria, pullers such as algae, or source-dipole swimmers as leading order models for active colloids [3].For a squirmer moving in bulk analytic solutions for the swimming speed and flow field exist, and the squirmer has so far extensivly been employed to capture effects of self-generated hydrodynamic flow fields on microswimmer dynamics [3,7].
Often microswimmers move through more complex fluids which sometimes fail to be described as a simple continuous Newtonian fluid.Recently there has been great interest in understanding the behavior and dynamics of microswimmers in such fluids, which is less well understood compared to swimming in Newtonian fluids [8].So far the squirmer moving in complex, viscoelastic, and nonhomogeneous fluids has been studied theoretically and by using continuum fluid models.Examples for nonhomogeneous but Newtonian fluids include a squirmer moving in concentric viscous fluid layers [9], or in viscosity gradients [10,11].Furthermore viscoelastic [12][13][14] and shear-thinning effects [15][16][17] have been investigated, as well as porosity modeled by Brinkman theory [18][19][20].
Polymeric fluids are often inhomogeneous and consist of nano-and microstructures as a result of the specific polymer composition.In particular in complex biological fluids such as mucus or in collagen fiber solutions the size of these structures can be strongly inhomogeneous and be up to some micron [21,22].Furthermore, in out-of-equlibrium situations such as for driven or active particles moving in supramolecular solutions, the polymer relaxation time can be larger compared to the typical time a driven or active particle spends to move its own size, which can leave behind a polymer-free zone [23][24][25].To capture heterogeneity of biological fluids such as mucus, two-fluid models have been employed which capture density-and viscosityinhomogeneities around squirmers [9], as well as porosity [26].However, such models so far do not capture nonisotropic effects induced by the aforementioned wake depletion of polymers, which depends on the specific local flow field created by different squirmer types.
Here we investigate the dynamics of squirmers in explicitly modeled fluids including coarse-grained polymers consisting of a relatively small number of relatively large monomers (see Fig. 1) where entanglements and viscoelastic effects are expected to play a minor role, but nonhomogeneous effects are present, depending on polymer density and stiffness in different self-generated squirmer fluid flows.We use the local fluid and polymer properties measured in the vicinity of the squirmer to determine underlying physical mechanisms for the observed swimming behavior.We find considerable differences in swimming speeds between pushers and pullers, which can not be explained by existing continuum fluid models.We show that local viscosity gradients slow down squirmers, as well as trapping and steric hindrance of large stiff polymers in the recirculation region in front of pushers.
Model.-We consider a coupled multi-component system, consisting of (i) a rigid squirmer, (ii) bead-spring polymers and filaments, (iii) the MPCD background fluid, and (iv) two bounding walls located at x = ±S X /2, where V = S X S Y S Z is the total volume of the simulation box with periodic boundary conditions in y-and z-direction, see Fig. 1.The simulation box size is fixed by S X = 72a 0 , S Y = S Z = 48a 0 .The two walls are included to avoid a possible net fluid-flow and an associated violation of angular momentum conservation along the x direction, i.e. the average squirmer direction [24].In the following we desribe how we model the individual components and how they are coupled through hydrodynamic and steric interactions.

MPCD.
Squirmer and polymers are immersed in a Newtonian background fluid such as water, which is modeled by MPCD and solves the Navier Stokes equations on a coarse-grained level including thermal fluctuations [28,29,49].The fluid is represented by N f effective, pointlike fluid particles with mass m, positions x i and velocities v i with i = 1, . . ., N f .The basic dynamics of the fluid particles consists of alternating streaming and collision steps.In the streaming step they move ballistically for a time δt and their positions are updated according to After the streaming step particles are sorted into cubic cells of length a 0 (see Fig. 1(c)), and in the collision step all particles in a cell exchange momentum which updates their velocities where v ξ = (1/N ξ ) j v j is the instantaneous average velocity in the cell with N ξ the total number of fluid particles in a cell, and j = 1, . . ., N ξ considers all particles in cell ξ.The temperature T in the cell is kept constant using an Anderson thermostat, i.e. using random velocities v r drawn from a Maxwell-Boltzmann distribution with variance k B T /m.The terms v P and v L ensure local linear and angular momentum conservation, respectively.In order to restore Galilean invariance and to minimize correlation effects we perform a random shift of the cell grid [29].Details of the algorithm are given in the SI.
As basic units of length, mass and energy we choose a 0 , m and k B T , respectively, and times are measured in units of t 0 = ma 2 0 /k B T .The free fluid model parameters are then the streaming time step δt and the average number of particles in a cell n = ⟨N ξ ⟩.We use δt = 0.02t 0 and n = 10 which models viscous flows at low Reynolds number Re and sufficiently high Schmidt number Sc [29,50].In the absence of polymers the viscosity is then η 0 = 16.04 mk B T /a 4 0 [51] and Sc ≈ 120 [52].
Polymer model.We model polymers and filaments as N beads of diameter σ = a 0 and mass m p = 10m connected by stiff harmonic springs with rest length l 0 = σ and spring constant k bond = 10 5 k B T /a 2 0 , where r i , i = 1, . . ., N are the bead positions and ∆r i = r i −r i−1 .The total number N p of polymers is varied to set the desired polymer density ρ = N N p πσ 3 /(6(V − V Sq )), defined as the volume fraction of all monomers in the simulation box with V Sq the squirmer volume.To model semiflexible polymers and filaments we use a bending potential where k b is the bending stiffness.In total we consider six different polymer models, i.e. flexible polymers  [53], for distances between beads r < σ and zero otherwise.We use ϵ 0 = k B T and σ * = σ/2 1/6 .

Squirmer model.
We employ the simplest form of a squirmer, a rigid sphere of radius R = 4a 0 with orientation e and a static, axisymmetric and tangential surface velocity v s = B 1 (1+βe•r s )[(e•r s )r s −e] (see also Fig. 1(c)) which pushes fluid backwards and moves the squirmer forwards, and depends on parameters B 1 > 0 and β.In bulk continuum fluids of viscosity η at Re = 0 the squirmer swims along its direction e with speed V b = 2 3 B 1 , which is independent of η and β.However, β determines the type of the flow field, and the squirmer is called a pusher for β < 0, a puller for β > 0, and a neutral squirmer for β = 0 [3,35,54].The flow field of a squirmer is independent of viscosity η, but the power consumption for a fixed v s to maintain the flow field depends linearly on η [36].For a pusher the far-field flow is that of an extensile force dipole, for a puller a contractile force dipole, and for a neutral squirmer a source dipole [3].We set the squirmer parameters B 1 = 0.062a 0 /t 0 and β = {−3, 0, 3}.In the absence of polymers the Reynolds number is Re 0 = V b Rρ f /η ≈ 0.10 with ρ f = m 0 n/a 3 0 , and Re < Re 0 in the presence of polymers.

Simulation procedure.
Initially the squirmer is located at R 0 = {18a 0 , 0, 0} pointing in negative x direction (e 0 = {−1, 0, 0}).Before the actual simulation we initialize and equlibrate for each parameter set an ensemble of 16 independent polymer solutions: Fluid particles and polymers are randomly distributed in the simulation box, but are not allowed to overlap with the squirmer.We then pre-equilibrate polymers for at least 5 • 10 4 MPCD time steps; here we replace the squirmer by a hard non-moving sphere of radius R.After equilibration we perform the actual simulations for N t = 5 • 10 4 MPCD time steps.The typical squirmer persistence time [3] τ r = 4πR 3 η/(k B T ) ≈ 1.3 • 10 4 t 0 is much larger than the simulation time ∆T = N t δt = 1000t 0 , and the associated Péclet number is Pe = 3V b τ r /(2R) ≈ 200 [3].Thus the squirmer moves more or less along the −x direction for a distance ∆x ≲ V b ∆T ≈ 40a 0 and the final position X f ≳ X 0 − ∆x ≈ −22a 0 is thus always still at least 3R away from the bottom wall located at x = −36a 0 .We measure time-and ensemble-averaged quantities for the last 50% of the simulation time where the system has reached a quasi steady-state.
To simulate the dynamics of the coupled system we use a hydrid MPCD-MD scheme [29,35].In the streaming step fluid particles move according to Eq. ( 8) for a time δt.Positions and velocities of the polymer beads evolve through molecular dynamics (MD) using a Velocity Verlet algorithm [55] with time step δt P = δ t /10 (except for k b = 3000k B T where δt P = δ t /50), and the forces from the potentials [Eqs.(3) -( 5)].The squirmer position and orientation dynamics evolves using a Velocity Verlet algorithm with time step δt S = δt.When fluid particles overlap with the squirmer or the walls a bounceback rule is applied where the squirmer surface velocity v s has to be included and linear and angular momenta are exchanged accordingly [35,36].This ensures the desired no-slip boundary condition and correct angular dynamics of the squirmer [35,37].We consider non-adsorbing polymers, hence squirmers and polymers interact with each other through a purely repulsive WCA potential similar as in Eq. ( 5) but with σ * = (R + σ/2)/2 1/6 and ϵ 0 is replaced by ϵ S = 1000k B T to model almost inpenetrable squirmers.In the collision step not only fluid particles exchange momentum, see Eq. ( 2), but polymer beads are coupled to the fluid by including them in the collision step [28,29].Furthermore, virtual particles inside the squirmers contribute to the collision step to better resolve the local flow fields [35][36][37].For details of the simulation procedure see SI.

Squirmer velocity.
The steady-state squirmer velocities V = ⟨V • e⟩ [35][36][37] are calculated from time-and ensemble averages, and are compared to the respective velocities V 0 determined by simulations in the absence of polymers.We obtain V 0 = 0.0406a 0 /t 0 (β = 0), V 0 = 0.0390a 0 /t 0 (β = −3), and V 0 = 0.0409a 0 /t 0 (β = +3), which deviate by a few percent from the theoretical value V b = 2B 1 /3 = 0.0413a 0 /t 0 for a squirmer in an infinite domain at Re = 0 because of small hydrodynamic squirmer-wall interactions [56] and small effects from finite Re [57,58].Fig. 2(a) shows squirmer velocities in solutions of short flexible polymers (N = 12, k b = 0), where V decreases with polymer density ρ for all squirmer types.A similar trend can be observed for longer flexible polymers (N = 30 and N = 100) in Fig. 2(b,c), but the decrease becomes stronger with increasing N .Furthermore, with increasing polymer density ρ pushers slow down faster compared to neutral squirmers and pullers.The decrease of V with ρ and squirmer type becomes even stronger in semiflexible polymer solutions (Fig. 2(d,e)).The largest effect is observed in the stiffest filamentous solutions (N = 30, k b = 3000k B T , see Fig. 2(f)), where pushers become several times slower compared to pullers at high densities.Such strong deviations between pushers and pullers have so far not been predicted by continuum models.It is thus clear that the flow microstructure and local fluid-squirmer interaction plays an important role.
Notably, we also measured the orientational squirmer dynamics and did not identify enhanced rotational diffusion of squirmers, see SI Fig S1, in contrast to recent work discussed in Ref. [47].This we attribute to the relatively hard polymer-squirmer interaction and the lack of polymer adsorbtion in our work, in contrast to Ref. [47].3(g,h) (solid blue and green curves), and is indeed very similar to the polymerfree case decaying as ∼ r −2 and ∼ r −3 , respectively, similar, as expected in simple homogeneous viscous fluids [3,7].Indeed we have previously demonstrated that solutions of short polymers can be approximated by Newtonian fluids of bulk viscosity η > η 0 , even under non-equilibrium conditions [24,25].
To quantify the effect of fluid viscosity we measure the ρ-dependent viscosity η of all considered fluids (see SI for details), as shown in Fig. 2(g).Interestingly, when we plot the squirmer velocities V depending on η [Fig.2(h)] we identify a clear decrease with η, while in continuum Newtonian fluids squirmers swim at constant velocity, independent of η [7].However, as shown in Fig. 4  all squirmer types due to polymer depletion, with a reduced viscosity around the squirmer, but the general local polymer density depends on squirmer type.
Theoretical two-fluid models have been used around sedimenting colloids [59] and squirmers [9], which are surrounded by a thin layer of thickness δ with viscosity η 0 and an ambient bulk fluid of viscosity η > η 0 .Indeed it had been shown that squirmers always slow down, depending on η/η 0 and δ/R, independent of the squirmer type [9].The theoretical curve from Ref. [9] for a single value δ/R = 0.2 fits more or less all our data for neutral squirmers and pushers [Fig.2(h)], in particular at sufficiently small viscosities [inset of Fig. 2(h)], although the polymer distribution is highly non-isotropic and squirmertype dependent [Fig.4(a)].Furthermore, the theoretical two-fluid model velocity decay fits the simulation results for sufficiently small η, as shown for N = 12, k b = 0 in Fig. 3(g,h).The value of δ/R = 0.2 only gives an effective polymer depletion thickness, and fails to descibe pullers accurately which slow down much weaker.While for both pushers and pullers polymers are more depleted at the back compared to the front because of the finite diffusion time of the polymers (similar as observed around driven colloidal particles [23,24]), pullers show a significant polymer accumulation in front due to the contractile flow field [Fig.3(b)].We speculate that this locally enhanced polymer density in front is responsible for pullers to be faster compared to pushers at low η, motivated by the theoretical prediction of speed enhancement for squirmers surrounded by a high-viscosity layer for the case η 0 > η [9].
As can be seen in Fig. 2, at larger polymer density, stiffness (and hence larger viscosity), the velocities for different squirmer types can deviate strongly.In these cases, also the flow fields are considerably weaker compared to the low-viscosity polymer solutions, as we show in Fig. 3(c 3(g,h), can for both pushers and pullers fitted by the two-fluid model [9] (black dotted lines), again by assuming δ/R = 0.2.In addition we fit the curves by a different two-fluid model from Ref. [20], which again uses a layer of viscosity η 0 close to the squirmer, but the outer bulk fluid is described as a Brinkman fluid as a model for porous medium with screening length λ [18,19].These solutions are shown in Fig. 3(g,h) as orange dotted lines, for δ/R = 0.2 and λ = 0.1R.All in all, while these models capture the flow fields relatively well, neither the viscous two-fluid model [9], nor models of a squirmer in a Brinkman medium [18][19][20] predict deviations in swimming speed between different squirmer types.
Interestingly, we observe large velocity deviations between pushers and pullers in stiff filamentous solutions although the polymer distribution around a pusher and a puller are similar, as shown in Fig. 4(b).It shows an additional layer of high polymer density above the lowdensity layer in front of the squirmer due to frequent encounters between the moving squirmers and the relatively large filaments, which we quantify below by steric squirmer-polymer interactions.
Our explicit approach not only gives access to local polymer distribution around the squirmers, but also allows to measure steric forces F st = i ∇ i V W CA (|r i − R|) between polymers and squirmers, induced by the overlap potential Eq. ( 5), which contribute in addition to pure hydrodynamic effects to the squirmer velocities.We plot the time-and ensemble averaged forces F st = −⟨F st • e⟩ for different fluids depending on viscosity η in Fig. 5(a), normalized by the bulk Stokes force F 0 = 6πηRV b ≈ 50.0kB T /a 0 of a sphere at velocity V b in a polymer-free fluid of viscosity η 0 .Note that F st > 0, i.e. steric forces act against squirmer direction e, as expected.For sufficiently small viscosities η pushers experience only a very small steric hindrance, as a consequence of the extensile flow field [Fig.3] pushing polymeric material away in front of the squirmer, while pullers accumulate more polymers in front because of the contractile flow field leading to much stronger steric forces.See also the Supplementary Movies M1 (pusher) and M2 (puller) demonstrating the local polymer dynamics in the reference frame of the squirmer.Since pullers move faster compared to pushers at low viscosities [inset of Fig. 2(h)] we conclude that the effect on the swimming speed due to different local polymer density [Fig.4(a)] of pullers (high density in front) compared to pushers (high density at side) outplays the stronger steric hindrance of pullers compared to pushers at small and moderate viscosities.
To identify the effective steric friction γ st experienced by the squirmers due to direct interactions with polymers we calculate γ st = F st /V which we compare to the Stokes bulk friction coefficient in the respective fluid of viscosity η, γ b = 6πηR.Fig. 5(b) shows γ st /γ b for relatively small η/η 0 , reiterating the fact that pushers experience higher steric friction compared to pullers.The inset of Fig. 5(b) shows the steric friction for large η/η 0 , and we highlight results at two particular viscosities which correspond to squirmers moving in stiff filamentous solutions (N = 30, k b = 3000k B T ) at respective densities ρ = 0.2 and ρ = 0.3 where pushers move much slower compared to pullers [see Fig. 2(f)].This can be explained by the fact that here the steric friction for pushers is about twice as large as compared to pullers, in stark contrast as observed at lower densities and particularly when compared to (semi-)flexible polymers (see also SI Fig. S2).We visualize the dynamics of stiff filaments around swimming pushers and pullers in Supplementary Movies M3 (pusher) and M4 (puller).While short flexible polymers (M1 and M2) easily pass by squirmers and act as large tracers in the surrounding flow, the longer stiffer filaments are less mobile (M3 and M4).Due to their size, individual filaments occupy larger regions around the squirmer, and do not change their conformations much because of their stiffness, see also local end-to-end distance in SI Fig. S3.Stiff filaments accumulate much stronger and stay longer in front of the pusher (M3), compared to the puller (M4).This can be explained by the flow fields in the reference frame of the squirmer, as shown in Fig. 3(d-f).While filaments in front of pullers experience a clear sidewise flow, pushers develop a vortex in front, as already known for polymer-free squirmers [3,54].For filaments located in different parts of the vortex their center-of-mass velocity cancels out in contrast to the polymer velocities in the more unidirectonal flows in front of pullers.This leads to more frequent encounters, reorientation of filaments perpendicular to the squirmer (see also SI Fig. S3), strong steric hindrance and hence to a strongly reduced swimming speed for pushers compared to pullers and neutral squirmers, for which flow vortices in the front are absent.
In general, in a simplified way we can interpret the steric friction as an effective load of effective radius R st = γ st /(6πη) the squirmer pushes (mainly) in front of it while swimming.In most of the cases this load is small, γ st < γ b and hence R st < R, while for example for a pusher in stiff filamentous fluids the effective load can be larger than the squirmer (γ st > γ b → R st > R), see inset Fig. 5(b).
Finally, we note that we also measured stretching of flexible polymers when moving around the squirmer.However, this effect is small (see also SI Fig. S3) Also, while shear-thinning effects are known to create local viscosity gradients and slowing down of squirmers, this effect is independent of squirmer type [16].Furthermore, since our polymers only consist of a relatively small number of monomers (N ≤ 100), entanglement effects are expected to be negligible [25].Notably, even when viscoelastic effects are accounted for, the effect on squirmer speed is typically small [13].
Conclusion. -Our work demonstrates the importance of non-continuum effects on the dynamics of different squirmers in fluids consisting of large polymers and filaments.First, non-homogeneous local polymer densities lead to non-homogeneous and squirmer-type specific viscosity gradients, which slow down squirmers but can only partly be captured by effective two-fluid models.How the specific polymer density landscape affects the squirmer speed remains to be investigated in more detail by refined theoretical models in future work.Second, because of their large size, polymers at sufficiently large size and stiffness act as an effective load where the squirmer pushes against, quantified by an effective steric friction.In this sense our approach bridges the gap between available continuum hydrodynamic models and purely non-hydrodynamic models of active Brownian particles moving in polymer solutions [60,61], by capturing non-continuum effects of squirmer hydrodynamics in complex fluids.
Finally we want to note that our results are expected to be relevant for real biological microswimmers moving in solutions of polymers and filaments which are comparable in size with the microswimmer itself.Examples are bacteria or sperm cells moving in mucus [8,62], where both the hydrodynamic and steric effects are expected to be of relevance for their transport in these heterogeneous fluid environments.* * * AZ acknowledges funding from the Austrian Science Fund (FWF) through a Lise-Meitner Fellowship (Grant No M 2458-N36).The computational results presented have been achieved using the Vienna Scientific Cluster (VSC).

Supplementary Information: Dynamics of squirmers in explicitly modeled polymeric fluids
Here we present the details of the MPCD simulations of a squirmer moving in a solution of polymers.The detailed description of the motion of fluid particles and squirmers is based on Refs.[35,51].
Initial setup.-The system consists of (i) a fixed number of N f pointlike MPCD fluid particles of mass m at positions x i and velocities v i , i = 1, . . ., N f , of (ii) N p polymers consisting each of N monomers of diameter σ, mass m p , positions r i and velocities w i , i = 1, . . ., N N p , of (iii) a squirmer of radius R and surface velocity modes B 1 and β with position R, orientation e, velocity V and angular velocity Ω, and of (iv) two solid walls located at ±S X /2.We employ periodic boundary conditions in y and z direction.The toal simulation volume V = S X S Y S Z is given by S X = 72a 0 , S Y = 48a 0 and S Z = 48a 0 .The length a 0 is the length of the MPCD simulation box, and the basic length scale in our system.Furthermore we employ m as the basic mass scale in the system, and k B T of the energy scale of the system, where T is the temperature of fluid kept constant by using a thermostat (see below).Then the unit of time is Initially the squirmer is placed at position R 0 = {18a 0 , 0, 0}, orientation e 0 = {−1, 0, 0}, velocity V 0 = 0 and angular velocity Ω 0 = 0.The total number of fluid particles is chosen such to fulfill an average number n of particles per unit volume a 3 0 , and is given by N f = n(V − V Sq )/a 3 0 with V Sq = 4πR 3 /3.Two parameters define the transport properties such as the viscosity of the fluid, namely the time step δt between two fluid collisions (see below), and the value of n.
The squirmer is assumed to be neutrally buoyant, i.e. it has the same mass density as the fluid, ρ Sq = ρ f = nma −3 0 , and its mass M S is given by M S = ρ Sq V Sq , and its moment of inertia by Initially, the effective fluid particles are randomly distributed in the simulation domain but are not allowed to overlap with the squirmer.The fluid velocities v i are Gaussian distributed with zero mean and standard deviation The number of polymers N p in the system is given by fixing the total volume fraction ρ = N N p πσ 3 /(6(V − V Sq )) of monomers in the simulation domain.The mass of the monomers is set to m p = 10m to achieve good coupling between fluid and polymers [29].Polymers are initially randomly distributed in the simulation domain with velocities w i drawn from a Gaussian distribution with zero mean and standard deviation σ p = k B T /m p .
Streaming step.-In the streaming step squirmers and fluid particles move ballistically for a time δt, while the dynamics of the polymers is captured by a molecular dynamics (MD) scheme, as described below.

Motion of the squirmers.
Before we integrate the equations of motion for the squirmers we have to calculate the forces acting on the squirmer through the overlap with polymer beads.These forces are exactly the steric forces described in the main text, F st = i ∇ i V W CA (|r i − R|), obtained from Eq. ( 5) in the main text but with σ * = (R + σ/2)/2 1/6 and ϵ 0 is replaced by ϵ S = 1000k B T .Furthermore, the squirmer assumes a torque from the interaction with the polymers, We use a simple Velocity-Verlet integration for a time δt.It updates positions and orientations according to and translational and angular velocities according to Motion of the fluid particles.In the streaming step the fluid particles simply move ballistically by When a fluid particle hits a wall, it is moved half a time step back and its velocity is updated according to the bounce-back rule to fulfill the no-slip boundary condition at the wall, where the velocity is simply reversed, A. Zöttl Then, the fluid particle is moved forward for half a time step with the new velocity.A similar procedure is applied if a fluid particle hits the squirmer, but we have to account for the fact that the surface of the squirmer at position r S contains a local velocity v s (e, r s ) and translation and rotation of the squirmer [36,37], Here, r * i is the collision point of the particle on the surface of the squirmer and v s (e, r * i ) is the surface velocity of the squirmer at position r * i .When a fluid particle interacts with a squirmer or a wall, its momentum p i = mv i is modified.While a fixed wall absorbs this momentum, moving objects such as colloids or squirmers update their velocity and angular velocity such that the total momentum and angular momentum is conserved during the collision.The change of momentum for fluid particle i hitting the squirmer is simply which is then transferred to the squirmer.In general, if N * fluid particles hit the squirmer during time interval δt, its velocity V and angular velocity Ω is updated to where are the total linear and angular momentum transferred to the squirmer in the streaming step by the fluid particles.
Motion of the polymers.The polymer beads move by molecular dynamics with time step δt P = δ t /10 (except for k b = 3000k B T where δt P = δ t /50), and the monomers interact with each other and with the squirmer through the potentials Eq. ( 3)- (5) in the main text.A total force F i on each monomer is calculated in each time step as the gradient of the potentials.Then their positions and velocities are updated according to and Collision step.-In the collision step fluid particles interact with each other but also interact with the squirmer and the walls via virtual particles [29,63], and with polymer beads.In order to fulfill Galilean invariance but also to reduce memory effects and correlations in the collision step, the cell grid is first shifted randomly at each time step by a random vector s = {s 1 , s 2 , s 3 }, where the components s i are drawn from a uniform distribution s i ∈ [−a 0 /2, a 0 /2].Then, collision cells partly overlapping with squirmers and walls are filled with virtual particles, which are placed in the walls and in the squirmer.This increases the accuracy of the hydrodynamic flow fields significantly since otherwise these cells would have an average fluid particle number below the mean number n and hence locally a smaller viscosity [35].
The virtual particles are placed randomly distributed at the same density as the fluid, ρ f , in a layer of thickness a 0 inside a bounding wall and in layers of thickness √ 3a 0 inside the squirmers.It is thus guaranteed that grid cells overlapping with walls and squirmers are always completely filled.The positions and velocities of the virtual particles are denoted by xi and vi , respectively, with i = 1, . . ., N v and N v is the total number of virtual particles.Their velocities are drawn from a normal distribution with the usual standard deviation, σ 0 = k B T /m, and zero mean.In addition, virtual particles located inside the squirmer also assume the local surface velocity of the point r * i on the squirmer surface, which is closest to xi , plus the velocity of this point due to the translation and rotation of the squirmer [37]: where vr i are the random velocities.Then fluid particles, virtual particles and the monomers are sorted into the cells, which all take part in the collision step.We denote the number of fluid, virtual and monomer particles in cell ξ by N f ξ , N v ξ and N m ξ , respectively.The total number of particles in the cell is where the total number of monomers in a cell, N m ξ , is typically either zero or one since σ = a 0 .In order to perform the collision step, the mean velocity in a cell is computed, We use an Anderson thermostat [52,64], where random velocities for all fluid particles, virtual and monomers have to be computed, which are denoted by v r i , vr i and w r i , respectively, and which are again drawn from a Gaussian distribution with width σ 0 = k B T /m for fluid and virtual particles and σ p = k B T /m p for the monomers.To conserve linear momentum the change of total velocity due to the added random velocities, V ξ , has to be calculated, given by Next, the center of mass has to be calculated, and the relative positions x s i = x i − r s ξ , xs i = xi − r s ξ and r s i = r i − r s ξ .They are needed to calculate the inverse of the moment of inertia tensor I −1 ξ in each cell, where The random velocities added in the collision step change angular momentum by To compensate for ∆L ξ and conserve angular momentum in each cell, the following angular velocity is computed [64], and used to rotate the particle velocities in a cell.Thus, by adding the extra terms ω i to the new velocities, angular momentum is conserved without changing linear momentum.To summarize, in the collision step the particle velocities are updated according to [64] v Then momentum and angular momentum is transferred to the squirmers: The change of momentum for a virtual particle is p-9 A. Zöttl The momenta and angular momenta of all virtual particles are then added up and assigned to the squirmer.Thus, after a collision step the squirmer assumes an additional momentum ∆P c and angular momentum ∆L c , where the sum goes over all virtual particles located in the squirmer with total number N v S .So, after each collision step the squirmer velocity and angular velocity are updated according to Viscosity Measurements with MPCD Poiseuille flow.-We measure the viscosity η of all fluid similar as in our previous work [25].It follows the method proposed in Ref. [65] for Newtonian fluids in the absence of polymers.We use a system size of S X = 60a 0 , S Y = S Z = 30a 0 and use periodic boundary conditions in all three dimensions and fluid particles and polymers are initially placed randomly.
Then the fluid particles are subjected to a constant but small acceleration force f a = 10 −3 k B T /a 0 [66] in z direction in one half of the simulation box (x > 0) and to f a = −10 −3 k B T /a 0 in z direction in the other half (x < 0).Then the streaming step has to be modified and positions and velocities of the fluid particles are updated according to, where f a = f a ẑ.Even in the absence of any walls this results in a periodic Poiseuille flow profile for sufficiently small f a , i.e. shear-thinning is still negligible.The fitted maximum velocity v max of the quadratic flow profile, averaged over time and 16 simulation runs, is linearly related to the inverse viscosity η −1 [25,65], In the absence of polymers we measure v 0 max = 0.070a 0 /t 0 which corresponds to viscosity η 0 = 16.04 mk B T /a 4 0 , as obtained in previous work [38,51].In the presence of polymers v max < v 0 max in the presence of polymers which enables us to determine all the viscosities η > η 0 presented in Fig. 2

Fig. 1 :
Fig. 1: (a) Motion of a squirmer (orange) in the simulation box bounded by two rigid walls at x = ±SX /2 and filled with a MPCD fluid (not shown) and polymers (gray polygons).(b) Typical local environement around a squirmer.Individual bead-spring polymers are shown in different colors.(c) 2D sketch of the multi-component system: squirmer oriented along direction e with surface velocity vs at surface rs, polymers (yellow, purple, cyan), and MPCD fluid particles (black dots).The grid for the collision step is shown in black.Not shown are virtual particles inside the squirmer used in the collision step.

Fig. 2 :
Fig. 2: (a-f) Mean swimming speed V of neutral squirmers (β = 0), pushers (β = −3) and pullers (β = 3) in fluids of different length N and bending stiffness k b depending on polymer density ρ, normalized by the swimming speed V0 in the absence of polymers.Error bars indicate the standard deviations of time-averaged steady-state swimming speeds of 16 independent realizations.(g) Measured fluid viscosities η normalized by the polymer-free viscosity η0.(h) Swimming speed V depending on viscosity η.

Fig. 3 :
Fig. 3: (a-f) Flow field around the squirmer in the lab frame (a-c) and in the co-moving frame (d-f).(g,h) Flow-field decay around the squirmer equator in the tangential (g) and (h) radial direction, normalized by local tangential surface velocity V l ; blue: pushers, green pullers; solid lines: N = 12, k b = 0; dotted lines: N = 30, k b = 3000; black lines: from theoretical model of Ref. [9]; orange lines: from theoretical model of Ref. [20].
Flow field and polymer distribution.In order to better understand our results, we measure the local fluid flow and polymer properties around the squirmer.The flow field has been measured both in the laboratory frame of reference, and in the co-moving frame of the squirmer, and has been averaged over ensemble, time, and along the azimuthal direction.Examples of the flow fields u f (r) = u θ θ + u r rs in the lab frame are shown in Fig. 3(b,c), and are compared to the flow fields in the absence of polymers [Fig.3(a)].Pusher and puller flow fields in solutions of short flexible polymers (N = 12, k b = 0, ρ = 0.1) [Fig.3(b)] are very similar to polymer-free velocity fields.The decay of |u r |(r) and |u θ |(r) along the equator of the squirmer (θ = π/2) is shown in Fig.
(a) for N = 12, k b = 0, ρ = 0.1, the time-and ensemble-averaged local polymer density ρ l around squirmers is non-homogeneous, induced by static and dynamic depletion effects because of finite polymer size and finite polymer diffusion time, and which depends on the squirmer type.In a thin layer around the squirmer the polymer density is very low for
) for stiff filamentous solutions (N = 30, k b = 3000k B T ) at density ρ = 0.2, i.e. at viscosity η = 14.1η 0 .The decay of the flow field components |u θ |(r) and |u r |(r) along the equator of the squirmer, see blue and green dotted lines in Fig.

Fig. 5 :
Fig. 5: (a) Normalized forces Fst from steric interactions between squirmers and polymers in different polymer solutions of viscosity η.(b) Effective steric friction coefficients γst = Fst/V normalized by the respective Stokes bulk friction coefficients γ b = 6πηR for small η (and for large η in inset).Pink arrows indicate stiff filaments at ρ = 0.2 and 0.3, respectively.
(g) in the main text.Movie M2.M2 shows the dynamics of a puller (β = +3) and its local polymeric environment, characterized by polymer length N = 12, stiffness k b = 0 and density ρ = 0.1, in the reference frame of the moving squirmer.Movie M3.M3 shows the dynamics of a pusher (β = −3) and its local polymeric environment, characterized by polymer length N = 30, stiffness k b = 3000k B T and density ρ = 0.2, in the reference frame of the moving squirmer.Movie M4.M4 shows the dynamics of a puller (β = +3) and its local polymeric environment, characterized by polymer length N = 30, stiffness k b = 3000k B T and density ρ = 0.2, in the reference frame of the moving squirmer.

Fig. S1 :Fig. S2 :
Fig. S1: Measured rotational diffusion constant Dr compared to the theoretical value in the absence of polymers D 0 r , for different fluids at different densities ρ. blue: pushers, green: pullers, red: neutral squirmers.

Fig. S3 :
Fig. S3: (a,b) Local polymer orientation angle α compared to the average bulk value α0.(c,d) Local polymer end-to-end distance compared to the average value l0.Shown are results for differents quirmer types and different fluids.