11institutetext:  Faculty of Physics, University of Vienna, Kolingasse 14-16,
1090 Wien, Austria
Institute for Theoretical Physics, TU Wien, Wiedner Hauptstr. 8-10, 1040 Wien,
Austria

# Dynamics of squirmers in explicitly modeled polymeric fluids

A. Zöttl 11221122

###### Abstract

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.

##  1 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 non-
homogeneous fluids has been studied theoretically and by using continuum fluid
models. Examples for non-homogeneous 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 viscosity-inhomogeneities around
squirmers [9], as well as porosity [26]. However, such models so far do not
capture non-isotropic effects induced by the aforementioned wake depletion of
polymers, which depends on the specific local flow field created by different
squirmer types.

To capture non-continuum effects, we conduct in this work hydrodynamic
simulations of explicitely modeled polymers around moving squirmers. We use a
coarse-grained particle-based method called multiparticle collision dynamics
(MPCD) [27] where polymers can be coupled to a background fluid very
efficiently in the presence of thermal noise [28, 29]. MPCD has been used in
the last years extensively to model the hydrodynamics of microswimmers [30].
Important examples are bacteria [31], sperm [32], pathogens [33], Janus
particles [34] and squirmers [35]. For example, squirmers in MPCD fluids have
been succesfully modeled to study swimming and persistent random walks [36,
37], pairwise interactions [37], dynamics in channel flow [38], near surfaces
[39], under gravity [40], or in nematic fluids [41], as well as their
collective motion [42, 43, 44, 45, 46]. MPCD also allows to model out-of-
equilibrium multi-component systems, such as active or driven particles in
polymer solutions [25, 24, 47]. Recently the orientational dynamics of
squirmers in solutions of polymers has been studied [47], motivated by
experiments of active colloids in viscoelastic polymer solutions [48].

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 non-homogeneous 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.

##  2 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=\pm S_{X}/2$, where
$\mathcal{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.

Figure 1: (a) Motion of a squirmer (orange) in the simulation box bounded by
two rigid walls at $x=\pm S_{X}/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
$\mathbf{e}$ with surface velocity $\mathbf{v}_{s}$ at surface
$\mathbf{r}_{s}$, 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.

###  2.1 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, 49, 29]. The fluid is
represented by $N_{f}$ effective, pointlike fluid particles with mass $m$,
positions $\mathbf{x}_{i}$ and velocities $\mathbf{v}_{i}$ with
$i=1,\dots,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 $\delta t$ and their positions are updated according
to

$\mathbf{x}_{i}(t+\delta t)=\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t)\delta t\,.$
(1)

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

$\mathbf{v}_{i}(t+\delta
t)=\mathbf{v}_{\xi}(t)+\mathbf{v}_{r}(t)+\mathbf{v}_{P}(t)+\mathbf{v}_{L}(t),$
(2)

where $\mathbf{v}_{\xi}=(1/N_{\xi})\sum_{j}\mathbf{v}_{j}$ is the
instantaneous average velocity in the cell with $N_{\xi}$ the total number of
fluid particles in a cell, and $j=1,\dots,N_{\xi}$ considers all particles in
cell $\xi$. The temperature $T$ in the cell is kept constant using an Anderson
thermostat, i.e. using random velocities $\mathbf{v}_{r}$ drawn from a
Maxwell-Boltzmann distribution with variance $k_{B}T/m$. The terms
$\mathbf{v}_{P}$ and $\mathbf{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}=\sqrt{ma_{0}^{2}/k_{B}T}$. The free fluid model parameters are then the
streaming time step $\delta t$ and the average number of particles in a cell
$n=\langle N_{\xi}\rangle$. We use $\delta t=0.02t_{0}$ and $n=10$ which
models viscous flows at low Reynolds number Re and sufficiently high Schmidt
number Sc [50, 29]. In the absence of polymers the viscosity is then
$\eta_{0}=16.04\sqrt{mk_{B}T/a_{0}^{4}}$ [51] and $\textnormal{Sc}\approx 120$
[52].

###  2.2 Polymer model

We model polymers and filaments as $N$ beads of diameter $\sigma=a_{0}$ and
mass $m_{p}=10m$ connected by stiff harmonic springs with rest length
$l_{0}=\sigma$ and spring constant $k_{bond}=10^{5}k_{B}T/a_{0}^{2}$,

$V_{bond}=\frac{1}{2}k_{bond}\sum_{i=2}^{N}(|\Delta\mathbf{r}_{i}|-l_{0})^{2}$
(3)

where $\mathbf{r}_{i}$, $i=1,\dots,N$ are the bead positions and
$\Delta\mathbf{r}_{i}=\mathbf{r}_{i}-\mathbf{r}_{i-1}$. The total number
$N_{p}$ of polymers is varied to set the desired polymer density
$\rho=NN_{p}\pi\sigma^{3}/(6(\mathcal{V}-\mathcal{V}_{Sq}))$, defined as the
volume fraction of all monomers in the simulation box with $\mathcal{V}_{Sq}$
the squirmer volume. To model semiflexible polymers and filaments we use a
bending potential

$V_{b}=\frac{1}{2}k_{b}\sum_{i=3}^{N}\left(\frac{\Delta\mathbf{r}_{i}\cdot\Delta\mathbf{r}_{i-1}}{|\Delta\mathbf{r}_{i}||\Delta\mathbf{r}_{i-1}|}-1\right)^{2}$
(4)

where $k_{b}$ is the bending stiffness. In total we consider six different
polymer models, i.e. flexible polymers ($k_{b}=0$) of length
$N=\\{12,30,100\\}$, and semiflexible filaments of length $N=30$ with
stiffness $k_{b}=\\{30k_{B}T,300k_{B}T,3000k_{B}T\\}$. Polymers are considered
at different densities $\rho=\\{0.05,0.1,0.2,0.3\\}$.

All polymer beads interact with each other through a purely repulsive soft
Weeks-Chandler-Anderson (WCA) potential [53],

$V_{WCA}(r)=4\epsilon_{0}\left[{\left(\frac{\sigma^{\ast}}{r}\right)}^{12}-{\left(\frac{\sigma^{\ast}}{r}\right)}^{6}\right]+\epsilon_{0}$
(5)

for distances between beads $r<\sigma$ and zero otherwise. We use
$\epsilon_{0}=k_{B}T$ and $\sigma^{\ast}=\sigma/2^{1/6}$.

###  2.3 Squirmer model

We employ the simplest form of a squirmer, a rigid sphere of radius $R=4a_{0}$
with orientation $\mathbf{e}$ and a static, axisymmetric and tangential
surface velocity
$\mathbf{v}_{s}=B_{1}(1+\beta\mathbf{e}\cdot\hat{\mathbf{r}}_{s})[(\mathbf{e}\cdot\hat{\mathbf{r}}_{s})\hat{\mathbf{r}}_{s}-\mathbf{e}]$
(see also Fig. 1(c)) which pushes fluid backwards and moves the squirmer
forwards, and depends on parameters $B_{1}>0$ and $\beta$. In bulk continuum
fluids of viscosity $\eta$ at $\textnormal{Re}=0$ the squirmer swims along its
direction $\mathbf{e}$ with speed $V_{b}=\frac{2}{3}B_{1}$, which is
independent of $\eta$ and $\beta$. However, $\beta$ determines the type of the
flow field, and the squirmer is called a pusher for $\beta<0$, a puller for
$\beta>0$, and a neutral squirmer for $\beta=0$ [3, 35, 54]. The flow field of
a squirmer is independent of viscosity $\eta$, but the power consumption for a
fixed $\mathbf{v}_{s}$ to maintain the flow field depends linearly on $\eta$
[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
$\beta=\\{-3,0,3\\}$. In the absence of polymers the Reynolds number is
$\textnormal{Re}_{0}=V_{b}R\rho_{f}/\eta\approx 0.10$ with
$\rho_{f}=m_{0}n/a_{0}^{3}$, and $\textnormal{Re}<\textnormal{Re}_{0}$ in the
presence of polymers.

Figure 2: (a-f) Mean swimming speed $V$ of neutral squirmers ($\beta=0$),
pushers ($\beta=-3$) and pullers ($\beta=3$) in fluids of different length $N$
and bending stiffness $k_{b}$ depending on polymer density $\rho$, normalized
by the swimming speed $V_{0}$ 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 $\eta$ normalized by
the polymer-free viscosity $\eta_{0}$. (h) Swimming speed $V$ depending on
viscosity $\eta$.

###  2.4 Simulation procedure

Initially the squirmer is located at $\mathbf{R}_{0}=\\{18a_{0},0,0\\}$
pointing in negative x direction ($\mathbf{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\cdot
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\cdot 10^{4}$ MPCD time steps. The typical squirmer persistence
time [3] $\tau_{r}=4\pi R^{3}\eta/(k_{B}T)\approx 1.3\cdot 10^{4}t_{0}$ is
much larger than the simulation time $\Delta T=N_{t}\delta t=1000t_{0}$, and
the associated Péclet number is $\textnormal{Pe}=3V_{b}\tau_{r}/(2R)\approx
200$ [3]. Thus the squirmer moves more or less along the $-x$ direction for a
distance $\Delta x\lesssim V_{b}\Delta T\approx 40a_{0}$ and the final
position $X_{f}\gtrsim X_{0}-\Delta x\approx-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 $\delta t$. Positions and velocities of the polymer beads evolve
through molecular dynamics (MD) using a Velocity Verlet algorithm [55] with
time step $\delta t_{P}=\delta_{t}/10$ (except for $k_{b}=3000k_{B}T$ where
$\delta t_{P}=\delta_{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 $\delta t_{S}=\delta t$. When fluid particles
overlap with the squirmer or the walls a bounce-back rule is applied where the
squirmer surface velocity $\mathbf{v}_{s}$ has to be included and linear and
angular momenta are exchanged accordingly [36, 35]. This ensures the desired
no-slip boundary condition and correct angular dynamics of the squirmer [37,
35]. 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 $\sigma^{\ast}=(R+\sigma/2)/2^{1/6}$ and $\epsilon_{0}$ is replaced
by $\epsilon_{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 [36, 37, 35]. For
details of the simulation procedure see SI.

Figure 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].

##  3 Results and discussion

###  3.1 Squirmer velocity

The steady-state squirmer velocities
$V=\langle\mathbf{V}\cdot\mathbf{e}\rangle$ [36, 37, 35] 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}$ ($\beta=0$), $V_{0}=0.0390a_{0}/t_{0}$ ($\beta=-3$),
and $V_{0}=0.0409a_{0}/t_{0}$ ($\beta=+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 $\textnormal{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 $\rho$ 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 $\rho$ pushers
slow down faster compared to neutral squirmers and pullers. The decrease of
$V$ with $\rho$ 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.2 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
$\mathbf{u}_{f}(\mathbf{r})=u_{\theta}\hat{\boldsymbol{\theta}}+u_{r}\hat{\mathbf{r}_{s}}$
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$, $\rho=0.1$) [Fig.
3(b)] are very similar to polymer-free velocity fields. The decay of
$|u_{r}|(r)$ and $|u_{\theta}|(r)$ along the equator of the squirmer
($\theta=\pi/2$) is shown in Fig. 3(g,h) (solid blue and green curves), and is
indeed very similar to the polymer-free case decaying as $\sim r^{-2}$ and
$\sim r^{-3}$, respectively, similar, as expected in simple homogeneous
viscous fluids [7, 3]. Indeed we have previously demonstrated that solutions
of short polymers can be approximated by Newtonian fluids of bulk viscosity
$\eta>\eta_{0}$, even under non-equilibrium conditions [25, 24].

To quantify the effect of fluid viscosity we measure the $\rho$-dependent
viscosity $\eta$ of all considered fluids (see SI for details), as shown in
Fig. 2(g). Interestingly, when we plot the squirmer velocities $V$ depending
on $\eta$ [Fig. 2(h)] we identify a clear decrease with $\eta$, while in
continuum Newtonian fluids squirmers swim at constant velocity, independent of
$\eta$ [7]. However, as shown in Fig. 4(a) for $N=12$, $k_{b}=0$, $\rho=0.1$,
the time- and ensemble-averaged local polymer density $\rho_{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 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 $\delta$
with viscosity $\eta_{0}$ and an ambient bulk fluid of viscosity
$\eta>\eta_{0}$. Indeed it had been shown that squirmers always slow down,
depending on $\eta/\eta_{0}$ and $\delta/R$, independent of the squirmer type
[9]. The theoretical curve from Ref. [9] for a single value $\delta/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 squirmer-type dependent
[Fig. 4(a)]. Furthermore, the theoretical two-fluid model velocity decay fits
the simulation results for sufficiently small $\eta$, as shown for $N=12$,
$k_{b}=0$ in Fig. 3(g,h). The value of $\delta/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 $\eta$, motivated by the theoretical prediction of speed
enhancement for squirmers surrounded by a high-viscosity layer for the case
$\eta_{0}>\eta$ [9].

Figure 4: Local polymer density $\rho_{l}$ compared to bulk density $\rho$
around pushers ($\beta=-3$) and pullers ($\beta=+3$) in solutions of (a) short
polymers and (b) long stiff filaments.

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) for
stiff filamentous solutions ($N=30$, $k_{b}=3000k_{B}T$) at density
$\rho=0.2$, i.e. at viscosity $\eta=14.1\eta_{0}$. The decay of the flow field
components $|u_{\theta}|(r)$ and $|u_{r}|(r)$ along the equator of the
squirmer, see blue and green dotted lines in Fig. 3(g,h), can for both pushers
and pullers fitted by the two-fluid model [9] (black dotted lines), again by
assuming $\delta/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 $\eta_{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 $\lambda$ [18, 19].
These solutions are shown in Fig. 3(g,h) as orange dotted lines, for
$\delta/R=0.2$ and $\lambda=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 low-density 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.

###  3.3 Direct squirmer-polymer forces

Our explicit approach not only gives access to local polymer distribution
around the squirmers, but also allows to measure steric forces
$\mathbf{F}_{st}=\sum_{i}\boldsymbol{\nabla}_{i}V_{WCA}(|\mathbf{r}_{i}-\mathbf{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}=-\langle\mathbf{F}_{st}\cdot\mathbf{e}\rangle$ for different fluids
depending on viscosity $\eta$ in Fig. 5(a), normalized by the bulk Stokes
force $F_{0}=6\pi\eta RV_{b}\approx 50.0k_{B}T/a_{0}$ of a sphere at velocity
$V_{b}$ in a polymer-free fluid of viscosity $\eta_{0}$. Note that $F_{st}>0$,
i.e. steric forces act against squirmer direction $\mathbf{e}$, as expected.
For sufficiently small viscosities $\eta$ 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.

Figure 5:  (a) Normalized forces $F_{st}$ from steric interactions between
squirmers and polymers in different polymer solutions of viscosity $\eta$. (b)
Effective steric friction coefficients $\gamma_{st}=F_{st}/V$ normalized by
the respective Stokes bulk friction coefficients $\gamma_{b}=6\pi\eta R$ for
small $\eta$ (and for large $\eta$ in inset). Pink arrows indicate stiff
filaments at $\rho=0.2$ and $0.3$, respectively.

To identify the effective steric friction $\gamma_{st}$ experienced by the
squirmers due to direct interactions with polymers we calculate
$\gamma_{st}=F_{st}/V$ which we compare to the Stokes bulk friction
coefficient in the respective fluid of viscosity $\eta$, $\gamma_{b}=6\pi\eta
R$. Fig. 5(b) shows $\gamma_{st}/\gamma_{b}$ for relatively small
$\eta/\eta_{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 $\eta/\eta_{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 $\rho=0.2$ and
$\rho=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}=\gamma_{st}/(6\pi\eta)$ the
squirmer pushes (mainly) in front of it while swimming. In most of the cases
this load is small, $\gamma_{st}<\gamma_{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 ($\gamma_{st}>\gamma_{b}\rightarrow 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\leq 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].

##  4 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.

###### Acknowledgements.

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. [51, 35].

##  5 Initial setup

The system consists of (i) a fixed number of $N_{f}$ pointlike MPCD fluid
particles of mass $m$ at positions $\mathbf{x}_{i}$ and velocities
$\mathbf{v}_{i}$, $i=1,\dots,N_{f}$, of (ii) $N_{p}$ polymers consisting each
of $N$ monomers of diameter $\sigma$, mass $m_{p}$, positions $\mathbf{r}_{i}$
and velocities $\mathbf{w}_{i}$, $i=1,\dots,NN_{p}$, of (iii) a squirmer of
radius $R$ and surface velocity modes $B_{1}$ and $\beta$ with position
$\mathbf{R}$, orientation $\mathbf{e}$, velocity $\mathbf{V}$ and angular
velocity $\boldsymbol{\Omega}$, and of (iv) two solid walls located at $\pm
S_{X}/2$. We employ periodic boundary conditions in $y$ and $z$ direction. The
toal simulation volume $\mathcal{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 $t_{0}=a_{0}\sqrt{m/k_{B}T}$.

Initially the squirmer is placed at position
$\mathbf{R}_{0}=\\{18a_{0},0,0\\}$, orientation $\mathbf{e}_{0}=\\{-1,0,0\\}$,
velocity $\mathbf{V_{0}}=0$ and angular velocity $\boldsymbol{\Omega}_{0}=0$.
The total number of fluid particles is chosen such to fulfill an average
number $n$ of particles per unit volume $a_{0}^{3}$, and is given by
$N_{f}=n(\mathcal{V}-\mathcal{V}_{Sq})/a_{0}^{3}$ with $\mathcal{V}_{Sq}=4\pi
R^{3}/3$. Two parameters define the transport properties such as the viscosity
of the fluid, namely the time step $\delta 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, $\rho_{Sq}=\rho_{f}=nma_{0}^{-3}$, and its mass $M_{S}$
is given by $M_{S}=\rho_{Sq}\mathcal{V}_{Sq}$, and its moment of inertia by
$I_{S}=\frac{2M_{S}R^{2}}{5}=\frac{8\pi nR^{5}}{15}\frac{m}{a_{0}^{3}}$.

Initially, the effective fluid particles are randomly distributed in the
simulation domain but are not allowed to overlap with the squirmer. The fluid
velocities $\mathbf{v}_{i}$ are Gaussian distributed with zero mean and
standard deviation $\sigma_{0}=\sqrt{k_{B}T/m}$.

The number of polymers $N_{p}$ in the system is given by fixing the total
volume fraction $\rho=NN_{p}\pi\sigma^{3}/(6(\mathcal{V}-\mathcal{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
$\mathbf{w}_{i}$ drawn from a Gaussian distribution with zero mean and
standard deviation $\sigma_{p}=\sqrt{k_{B}T/m_{p}}$.

##  6 Streaming step

In the streaming step squirmers and fluid particles move ballistically for a
time $\delta t$, while the dynamics of the polymers is captured by a molecular
dynamics (MD) scheme, as described below.

###  6.1 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,
$\mathbf{F}_{st}=\sum_{i}\boldsymbol{\nabla}_{i}V_{WCA}(|\mathbf{r}_{i}-\mathbf{R}|)$,
obtained from Eq. (5) in the main text but with
$\sigma^{\ast}=(R+\sigma/2)/2^{1/6}$ and $\epsilon_{0}$ is replaced by
$\epsilon_{S}=1000k_{B}T$. Furthermore, the squirmer assumes a torque from the
interaction with the polymers,
$\mathbf{T}_{st}=\sum_{i}(\mathbf{R}-\mathbf{r}_{i})\times\boldsymbol{\nabla}_{i}V_{WCA}(|\mathbf{r}_{i}-\mathbf{R}|)$

We use a simple Velocity-Verlet integration for a time $\delta t$. It updates
positions and orientations according to

$\begin{split}\mathbf{R}(t+\delta
t)&=\mathbf{R}(t)+\left[\mathbf{V}(t)+\frac{1}{2M_{S}}\mathbf{F}_{st}(t)\delta
t\right]\delta t,\\\ \mathbf{e}(t+\delta
t)&=\mathbf{e}(t)+\left[\left(\boldsymbol{\Omega}(t)+\frac{1}{2I_{S}}\mathbf{T}_{st}(t)\delta
t\right)\times\mathbf{e}(t)\right]\delta t,\end{split}$ (6)

and translational and angular velocities according to

$\begin{split}\mathbf{V}(t+\delta
t)&=\mathbf{V}(t)+\frac{1}{2M_{S}}\left[\mathbf{F}_{st}(t)+\mathbf{F}_{st}(t+\delta
t)\right]\delta t,\\\ \boldsymbol{\Omega}(t+\delta
t)&=\boldsymbol{\Omega}(t)+\frac{1}{2I_{S}}\left[\mathbf{T}_{st}(t)+\mathbf{T}_{st}(t+\delta
t)\right]\delta t.\end{split}$ (7)

###  6.2 Motion of the fluid particles

In the streaming step the fluid particles simply move ballistically by

$\mathbf{x}_{i}(t+\delta t)=\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t)\delta t\,.$
(8)

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,

$\mathbf{v}_{i}^{\prime}=-\mathbf{v}_{i}\,,$ (9)

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 $\mathbf{r}_{S}$ contains a local velocity
$\mathbf{v}_{s}(\mathbf{e},\mathbf{r}_{s})$ and translation and rotation of
the squirmer [36, 37],

$\mathbf{v}_{i}^{\prime}=-\mathbf{v}_{i}+2\left[\mathbf{v}_{s}(\mathbf{e},\mathbf{r}_{i}^{\ast})+\boldsymbol{\Omega}\times(\mathbf{r}_{i}^{\ast}-\mathbf{R})+\mathbf{V}\right]\,.$
(10)

Here, $\mathbf{r}_{i}^{\ast}$ is the collision point of the particle on the
surface of the squirmer and $\mathbf{v}_{s}(\mathbf{e},\mathbf{r}_{i}^{\ast})$
is the surface velocity of the squirmer at position $\mathbf{r}_{i}^{\ast}$.

When a fluid particle interacts with a squirmer or a wall, its momentum
$\mathbf{p}_{i}=m\mathbf{v}_{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

$\Delta\mathbf{p}_{i}=m[\mathbf{v}_{i}^{\prime}-\mathbf{v}_{i}]\,,$ (11)

which is then transferred to the squirmer. In general, if $N^{\ast}$ fluid
particles hit the squirmer during time interval $\delta t$, its velocity
$\mathbf{V}$ and angular velocity $\boldsymbol{\Omega}$ is updated to

$\begin{split}\mathbf{V}^{\prime}&=\mathbf{V}+\frac{1}{M_{S}}\Delta\mathbf{P}^{s},\\\
\boldsymbol{\Omega}^{\prime}&=\boldsymbol{\Omega}+\frac{1}{I_{S}}\Delta\mathbf{L}^{s},\end{split}$
(12)

where

$\begin{split}\Delta\mathbf{P}^{s}&=-\sum_{i=1}^{N^{\ast}}\Delta\mathbf{p}_{i},\\\
\Delta\mathbf{L}^{s}&=-\sum_{i=1}^{N^{\ast}}\left[(\mathbf{r}_{i}^{\ast}-\mathbf{R})\times\Delta\mathbf{p}_{i}\right].\end{split}$
(13)

are the total linear and angular momentum transferred to the squirmer in the
streaming step by the fluid particles.

###  6.3 Motion of the polymers

The polymer beads move by molecular dynamics with time step $\delta
t_{P}=\delta_{t}/10$ (except for $k_{b}=3000k_{B}T$ where $\delta
t_{P}=\delta_{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
$\mathbf{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

$\mathbf{r}_{i}(t+\delta
t_{p})=\mathbf{r}_{i}(t)+\left[\mathbf{w}_{i}(t)+\frac{1}{2m}\mathbf{F}_{i}(t)\delta
t_{p}\right]\delta t_{p}$ (14)

and

$\begin{split}\mathbf{w}_{i}(t+\delta
t_{p})=\mathbf{w}_{i}(t)+\frac{1}{2m}\left[\mathbf{F}_{i}(t)+\mathbf{F}_{i}(t+\delta
t_{p})\right]\delta t_{p}\,.\end{split}$ (15)

##  7 Collision step

In the collision step fluid particles interact with each other but also
interact with the squirmer and the walls via virtual particles [63, 29], 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
$\mathbf{s}=\\{s_{1},s_{2},s_{3}\\}$, where the components $s_{i}$ are drawn
from a uniform distribution $s_{i}\in[-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, $\rho_{f}$, in a layer of thickness $a_{0}$ inside a bounding wall
and in layers of thickness $\sqrt{3}a_{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 $\bar{\mathbf{x}}_{i}$ and $\bar{\mathbf{v}}_{i}$, respectively,
with $i=1,\dots,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, $\sigma_{0}=\sqrt{k_{B}T/m}$, and zero mean. In addition, virtual
particles located inside the squirmer also assume the local surface velocity
of the point $\bar{\mathbf{r}}_{i}^{\ast}$ on the squirmer surface, which is
closest to $\bar{\mathbf{x}}_{i}$, plus the velocity of this point due to the
translation and rotation of the squirmer [37]:

$\bar{\mathbf{v}}_{i}=\bar{\mathbf{v}}_{i}^{\text{r}}+\mathbf{v}_{s}(\mathbf{e},\bar{\mathbf{r}}_{i}^{\ast})+\boldsymbol{\Omega}\times(\bar{\mathbf{r}}_{i}^{\ast}-\mathbf{R})+\mathbf{V},$
(16)

where $\bar{\mathbf{v}}_{i}^{\text{r}}$ 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 $\xi$ by $N_{\xi}^{f}$,
$N_{\xi}^{v}$ and $N_{\xi}^{m}$, respectively. The total number of particles
in the cell is $N_{\xi}=N_{\xi}^{f}+N_{\xi}^{v}+N_{\xi}^{m}$ where the total
number of monomers in a cell, $N_{\xi}^{m}$, is typically either zero or one
since $\sigma=a_{0}$. In order to perform the collision step, the mean
velocity in a cell is computed,

$\mathbf{u}_{\xi}=\frac{1}{mN_{\xi}^{f}+mN_{\xi}^{v}+m_{p}N_{\xi}^{m}}\left(m\sum_{i=1}^{N_{\xi}^{f}}\mathbf{v}_{i}+m\sum_{i=1}^{N_{\xi}^{v}}\bar{\mathbf{v}}_{i}+m_{p}\sum_{i=1}^{N_{\xi}^{m}}\mathbf{w}_{i}\right)\,.$
(17)

We use an Anderson thermostat [64, 52], where random velocities for all fluid
particles, virtual particles and monomers have to be computed, which are
denoted by $\mathbf{v}_{i}^{\text{r}}$, $\bar{\mathbf{v}}_{i}^{\text{r}}$ and
$\mathbf{w}_{i}^{\text{r}}$, respectively, and which are again drawn from a
Gaussian distribution with width $\sigma_{0}=\sqrt{k_{B}T/m}$ for fluid and
virtual particles and $\sigma_{p}=\sqrt{k_{B}T/m_{p}}$ for the monomers. To
conserve linear momentum the change of total velocity due to the added random
velocities, $\boldsymbol{\mathcal{V}}_{\xi}$, has to be calculated, given by

$\boldsymbol{\mathcal{V}}_{\xi}=\frac{1}{mN_{\xi}^{f}+mN_{\xi}^{v}+m_{p}N_{\xi}^{m}}\left(m\sum_{i=1}^{N_{\xi}^{f}}\mathbf{v}_{i}^{r}+m\sum_{i=1}^{N_{\xi}^{v}}\bar{\mathbf{v}}_{i}^{r}+m_{p}\sum_{i=1}^{N_{\xi}^{m}}\mathbf{w}_{i}^{r}\right)\,.$
(18)

Next, the center of mass has to be calculated,

$\mathbf{r}_{\xi}^{s}=\frac{1}{mN_{\xi}^{f}+mN_{\xi}^{v}+m_{p}N_{\xi}^{m}}\left(m\sum_{i=1}^{N_{\xi}^{f}}\mathbf{x}_{i}+m\sum_{i=1}^{N_{\xi}^{v}}\bar{\mathbf{x}}_{i}+m_{p}\sum_{i=1}^{N_{\xi}^{m}}\mathbf{r}_{i}\right)\,,$
(19)

and the relative positions
$\mathbf{x}_{i}^{s}=\mathbf{x}_{i}-\mathbf{r}_{\xi}^{s}$,
$\bar{\mathbf{x}}_{i}^{s}=\bar{\mathbf{x}}_{i}-\mathbf{r}_{\xi}^{s}$ and
$\mathbf{r}_{i}^{s}=\mathbf{r}_{i}-\mathbf{r}_{\xi}^{s}$. They are needed to
calculate the inverse of the moment of inertia tensor $\mathbf{I}^{-1}_{\xi}$
in each cell, where

$\mathbf{I}_{\xi}=m\sum_{i=1}^{N_{\xi}^{f}}(|\mathbf{x}_{i}^{s}|^{2}\,\textbf{{1}}-\mathbf{x}_{i}^{s}\otimes\mathbf{x}_{i}^{s})+m\sum_{i=1}^{N_{\xi}^{v}}(|\bar{\mathbf{x}}_{i}^{s}|^{2}\,\textbf{{1}}-\bar{\mathbf{x}}_{i}^{s}\otimes\bar{\mathbf{x}}_{i}^{s})+m_{p}\sum_{i=1}^{N_{\xi}^{m}}(|\mathbf{r}_{i}^{s}|^{2}\,\textbf{{1}}-\mathbf{r}_{i}^{s}\otimes\mathbf{r}_{i}^{s}).$
(20)

The random velocities added in the collision step change angular momentum by

$\Delta\boldsymbol{\mathcal{L}}_{\xi}=m\sum_{i=1}^{N_{\xi}^{f}}[\mathbf{x}_{i}^{s}\times(\mathbf{v}_{i}-\mathbf{v}_{i}^{\text{r}})]+m\sum_{i=1}^{N_{\xi}^{v}}[\bar{\mathbf{x}}_{i}^{s}\times(\bar{\mathbf{v}}_{i}-\bar{\mathbf{v}}_{i}^{\text{r}})]+m_{p}\sum_{i=1}^{N_{\xi}^{m}}[\mathbf{r}_{i}^{s}\times(\mathbf{w}_{i}-\mathbf{w}_{i}^{\text{r}})].$
(21)

To compensate for $\Delta\boldsymbol{\mathcal{L}}_{\xi}$ and conserve angular
momentum in each cell, the following angular velocity is computed [64],

$\boldsymbol{\omega}_{\xi}^{\text{AMC}}=\mathbf{I}_{\xi}^{-1}\Delta\boldsymbol{\mathcal{L}}_{\xi},$
(22)

and used to rotate the particle velocities in a cell. Thus, by adding the
extra terms $\boldsymbol{\omega}_{\xi}^{\text{AMC}}\times\mathbf{x}_{i}^{s}$,
$\boldsymbol{\omega}_{\xi}^{\text{AMC}}\times\bar{\mathbf{x}}_{i}^{s}$ or
$\boldsymbol{\omega}_{\xi}^{\text{AMC}}\times\mathbf{r}_{i}^{s}$ 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]

$\begin{split}\mathbf{v}_{i}^{\prime}&=\mathbf{u}_{\xi}+\mathbf{v}_{i}^{\text{r}}-\boldsymbol{\mathcal{V}}_{\xi}+\boldsymbol{\omega}_{\xi}^{\text{AMC}}\times\mathbf{x}_{i}^{s},\\\
\bar{\mathbf{v}}_{i}^{\prime}&=\mathbf{u}_{\xi}+\bar{\mathbf{v}}_{i}^{\text{r}}-\boldsymbol{\mathcal{V}}_{\xi}+\boldsymbol{\omega}_{\xi}^{\text{AMC}}\times\bar{\mathbf{x}}_{i}^{s}\,,\\\
\mathbf{w}_{i}^{\prime}&=\mathbf{u}_{\xi}+\mathbf{w}_{i}^{\text{r}}-\boldsymbol{\mathcal{V}}_{\xi}+\boldsymbol{\omega}_{\xi}^{\text{AMC}}\times\mathbf{r}_{i}^{s}.\\\
\end{split}$ (23)

Then momentum and angular momentum is transferred to the squirmers: The change
of momentum for a virtual particle is

$\Delta\bar{\mathbf{p}}_{i}=m(\bar{\mathbf{v}}_{i}^{\prime}-\bar{\mathbf{v}}_{i}).$
(24)

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 $\Delta\mathbf{P}^{c}$ and angular momentum
$\Delta\mathbf{L}^{c}$,

$\begin{split}\Delta\mathbf{P}^{c}&=\sum_{i=1}^{N_{S}^{v}}\Delta\bar{\mathbf{p}}_{i},\\\
\Delta\mathbf{L}^{c}&=\sum_{i=1}^{N_{S}^{v}}(\bar{\mathbf{r}}_{i}-\mathbf{R})\times\Delta\bar{\mathbf{p}}_{i},\\\
\end{split}$ (25)

where the sum goes over all virtual particles located in the squirmer with
total number $N_{S}^{v}$. So, after each collision step the squirmer velocity
and angular velocity are updated according to

$\begin{split}\mathbf{V}^{\prime}&=\mathbf{V}+\frac{1}{M_{S}}\Delta\mathbf{P}^{c},\\\
\boldsymbol{\Omega}^{\prime}&=\boldsymbol{\Omega}+\frac{1}{I_{S}}\Delta\mathbf{L}^{c}\,.\end{split}$
(26)

##  8 Viscosity Measurements with MPCD Poiseuille flow

We measure the viscosity $\eta$ 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,

$\begin{split}\mathbf{x}_{i}(t+\delta
t)&=\mathbf{x}(t)+\mathbf{v}_{i}(t)\delta t+\frac{1}{2m}\mathbf{f}_{a}\delta
t^{2},\\\ \mathbf{v}_{i}(t+\delta
t)&=\mathbf{v}_{i}(t)+\frac{1}{m}\mathbf{f}_{a}\delta t.\end{split}$ (27)

where $\mathbf{f}_{a}=f_{a}\hat{\mathbf{z}}$. 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 $\eta^{-1}$ [65, 25],

$\eta=\frac{nf_{a}S_{X}^{2}}{8a_{0}^{3}v_{max}}$ (28)

In the absence of polymers we measure $v_{max}^{0}=0.070a_{0}/t_{0}$ which
corresponds to viscosity $\eta_{0}=16.04\sqrt{mk_{B}T/a_{0}^{4}}$, as obtained
in previous work [38, 51]. In the presence of polymers $v_{max}<v_{max}^{0}$
in the presence of polymers which enables us to determine all the viscosities
$\eta>\eta_{0}$ presented in Fig. 2(g) in the main text.

##  9 Supplementary Figures

###  9.1 Rotational diffusion

We measure the orientational correlation $C_{e}(t)$ of the squirmer,
$C_{e}(t)=\langle\mathbf{e}(0)\cdot\mathbf{e}(t)\rangle$, where time $t=0$ is
set after $50\%$ of the simulation time. From $C_{e}(t)$ we fit the
theoretical curve $C_{e}^{th}(t)=e^{-2D_{r}t}$ to obtain the rotational
diffusion constants $D_{r}$. We show the results in Fig. S1 where we compare
$D_{r}$ to the theoretical value in the absence of polymers,
$D_{r}^{0}=k_{B}T/(8\pi\eta_{0}R^{3})$.

###  9.2 Steric friction

We plot the effective steric friction coefficient $\gamma_{st}$ devided by the
bulk friction coefficient $\gamma_{b}$, similar as in Fig. 5(b) in the main
text but now depending on polymer length, stiffness and density in Fig. S2.

###  9.3 Local polymer properties

In Fig. S3(a,b) we show the local polymer alignment angle $\alpha$ compared to
the average angle $\alpha_{0}$ defined as the angle between the polymer end-
to-end vector and the squirmer orientation for two different fluids. Here
$\alpha=0$ corresponds to alignement in the direction of the squirmer, and
$\alpha=\pi/2$ perpendicular to it.

In Fig. S3(c,d) we show the local polymer end-to-end distance $l$ compared to
the average value $l_{0}$ for the same two fluids, which indicates regions of
small stretching (typically $l/l_{0}<1.3$) for flexible polymers and small
compression for stiff filaments.

##  10 Supplementary Movies

###  10.1 Movie M1

M1 shows the dynamics of a pusher ($\beta=-3$) and its local polymeric
environment, characterized by polymer length $N=12$, stiffness $k_{b}=0$ and
density $\rho=0.1$, in the reference frame of the moving squirmer.

###  10.2 Movie M2

M2 shows the dynamics of a puller ($\beta=+3$) and its local polymeric
environment, characterized by polymer length $N=12$, stiffness $k_{b}=0$ and
density $\rho=0.1$, in the reference frame of the moving squirmer.

###  10.3 Movie M3

M3 shows the dynamics of a pusher ($\beta=-3$) and its local polymeric
environment, characterized by polymer length $N=30$, stiffness
$k_{b}=3000k_{B}T$ and density $\rho=0.2$, in the reference frame of the
moving squirmer.

###  10.4 Movie M4

M4 shows the dynamics of a puller ($\beta=+3$) and its local polymeric
environment, characterized by polymer length $N=30$, stiffness
$k_{b}=3000k_{B}T$ and density $\rho=0.2$, in the reference frame of the
moving squirmer.

Figure S1: Measured rotational diffusion constant $D_{r}$ compared to the
theoretical value in the absence of polymers $D_{r}^{0}$, for different fluids
at different densities $\rho$. blue: pushers, green: pullers, red: neutral
squirmers. Figure S2: Effective steric friction coefficient $\gamma_{st}$ for
different polymer types depending on density. The data points are the same as
used in Fig. 5(b) in the main text. Figure S3: (a,b) Local polymer orientation
angle $\alpha$ compared to the average bulk value $\alpha_{0}$. (c,d) Local
polymer end-to-end distance compared to the average value $l_{0}$. Shown are
results for differents quirmer types and different fluids.

## References

  * [1] Lauga E. Powers T. R. Rep. Prog. Phys.722009096601.
  * [2] Elgeti J., Winkler R. G. Gompper G. Rep. Prog. Phys.782015056601.
  * [3] Zöttl A. Stark H. J. Phys. Condens. Matter282016253001.
  * [4] Bechinger C., Di Leonardo R., Löwen H., Reichhardt C., Volpe G. Volpe G. Rev. Mod. Phys.882016045006.
  * [5] Lighthill J. M. Commun. Pure Appl. Math.51952109.
  * [6] Blake J. R. J. Fluid Mech.461971199.
  * [7] Pedley T. J. IMA J. Appl. Math. (Institute Math. Its Appl.812016488.
  * [8] Li G., Lauga E. Ardekani A. M. J. Nonnewton. Fluid Mech.2972021104655.
  * [9] Reigh S. Y. Lauga E. Phys. Rev. Fluids22017.
  * [10] Datt C. Elfring G. J. Phys. Rev. Lett.1232019158006.
  * [11] Eastham P. S. Shoele K. Phys. Rev. Fluids520201.
  * [12] Zhu L., Do-Quang M., Lauga E. Brandt L. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys.8320111.
  * [13] Zhu L., Lauga E. Brandt L. Phys. Fluids24201251902.
  * [14] De Corato M., Greco F. Maffettone P. L. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys.9220151.
  * [15] Montenegro-Johnson T. D., Smith D. J. Loghin D. Phys. Fluids252013.
  * [16] Datt C., Zhu L., Elfring G. J. Pak O. S. J. Fluid Mech.7842015R1.
  * [17] Ouyang Z., Lin J. Ku X. Phys. Fluids302018.
  * [18] Leshansky A. M. Phys. Rev. E8020091.
  * [19] Nganguia H. Pak O. S. J. Fluid Mech.8552018554.
  * [20] Nganguia H., Zheng K., Chen Y., Pak O. S. Zhu L. Phys. Fluids322020.
  * [21] Kirch J., Schneider A., Abou B., Hopf A., Schaefer U. F., Schneider M., Schall C., Wagner C. Lehr C. M. Proc. Natl. Acad. Sci. U. S. A.109201218355.
  * [22] Mickel W., Münster S., Jawerth L. M., Vader D. A., Weitz D. A., Sheppard A. P., Mecke K., Fabry B. Schröder-Turk G. E. Biophys. J.9520086072.
  * [23] Gutsche C., Kremer F., Krüger M., Rauscher M., Weeber R. Harting J. J. Chem. Phys.1292008.
  * [24] Zöttl A. Yeomans J. M. J. Phys. Condens. Matter312019234001.
  * [25] Zöttl A. Yeomans J. M. Nat. Phys.152019554.
  * [26] Nganguia H., Zhu L., Palaniappan D. Pak O. S. Phys. Rev. E10120201.
  * [27] Malevanets A. Kapral R. J. Chem. Phys.11019998605.
  * [28] Malevanets A. Yeomans J. M. Eur. Lett.522000231.
  * [29] Gompper G., Ihle T., Kroll D. M. Winkler R. G. Adv. Poly. Sci.22120091.
  * [30] Zöttl A. Chin. Phys. B292020.
  * [31] Hu J., Yang M., Gompper G. Winkler R. G. Soft Matter1120157867.
  * [32] Elgeti J., Kaupp U. B. Gompper G. Biophys. J.9920101018.
  * [33] Babu S. B. Stark H. New J. Phys.14201285012.
  * [34] Huang M. J., Schofield J. Kapral R. Soft Matter1220165581.
  * [35] Zöttl A. Stark H. Eur. Phys. J. E41201861.
  * [36] Downton M. T. Stark H. J. Phys. Condens. Matter212009204101.
  * [37] Götze I. O. Gompper G. Phys. Rev. E822010041921.
  * [38] Zöttl A. Stark H. Phys. Rev. Lett.1082012218104.
  * [39] Schaar K., Zöttl A. Stark H. Phys. Rev. Lett.1152015038101.
  * [40] Rühle F., Blaschke J., Kuhr J. T. Stark H. New J. Phys.202018025003.
  * [41] Mandal S. Mazza M. G. Eur. Phys. J. E4420211.
  * [42] Zöttl A. Stark H. Phys. Rev. Lett.1122014118101.
  * [43] Blaschke J., Maurer M., Menon K., Zöttl A. Stark H. Soft Matter1220169821.
  * [44] Theers M., Westphal E., Qi K., Winkler R. G. Gompper G. Soft Matter1420188590.
  * [45] Zantop A. W. Stark H. J. Chem. Phys.1552021.
  * [46] Qi K., Westphal E., Gompper G. Winkler R. G. Commun. Phys.5202249.
  * [47] Qi K., Westphal E., Gompper G. Winkler R. G. Phys. Rev. Lett.124202068001.
  * [48] Gomez-Solano J. R., Blokhuis A. Bechinger C. Phys. Rev. Lett.11620162.
  * [49] Kapral R. Adv. Chem. Phys.140200889.
  * [50] Padding J. T. Louis A. A. Phys. Rev. E742006031402.
  * [51] Zöttl A. Hydrodynamics of Microswimmers in Confinement and in Poiseuille Flow Ph.D. thesis TU Berlin (2014).
  * [52] Noguchi H. Gompper G. Phys. Rev. E782008.
  * [53] Weeks J. D., Chandler D. Anderson H. C. J. Chem. Phys.5419715237.
  * [54] Zöttl A. Stark H. Annu. Rev. Condens. Matter Phys.142023109.
  * [55] Allen M. P., Frenkel D. Talbot J. Comput. Phys. Reports91989301.
  * [56] Berke A., Turner L., Berg H. Lauga E. Phys. Rev. Lett.101200838102.
  * [57] Wang S. Ardekani a. Phys. Fluids242012.
  * [58] Khair A. S. Chisholm N. G. Phys. Fluids262014.
  * [59] Fan T. H., Dhont J. K. Tuinier R. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys.7520071.
  * [60] Du Y., Jiang H. Hou Z. Soft Matter1520192020.
  * [61] Kim Y., Joo S., Kim W. K. Jeon J. H. Macromolecules5520227136.
  * [62] Figueroa-Morales N., Dominguez-Rubio L., Ott T. L. Aranson I. S. Sci. Rep.920199713.
  * [63] Lamura A., Gompper G., Ihle T. Kroll D. M. Eur. Lett.562001319.
  * [64] Noguchi H., Kikuchi N. Gompper G. Eur. Lett.78200710005.
  * [65] Backer J. A., Lowe C. P., Hoefsloot H. C. J. Iedema P. D. J. Chem. Phys.1222005154503.
  * [66] Allahyarov E. Gompper G. Phys. Rev. E662002036702.