Simulating squirmers with multiparticle collision dynamics

Multiparticle collision dynamics is a modern coarse-grained simulation technique to treat the hydrodynamics of Newtonian fluids by solving the Navier-Stokes equations. Naturally, it also includes thermal noise. Initially it has been applied extensively to spherical colloids or bead-spring polymers immersed in a fluid. Here, we review and discuss the use of multiparticle collision dynamics for studying the motion of spherical model microswimmers called squirmers moving in viscous fluids.


Introduction
Recently, there has been a great interest in understanding the physical principles of the locomotion mechanisms of microswimmers.Important examples are biological microorganisms such as bacteria, sperm cells or algae, and artificial self-propelled swimmers such as active colloids [1,2] or active emulsion droplets [3].Their locomotion is mainly governed by low Reynolds number hydrodynamics and thermal or biological noise.While swimming speed and flow fields of a single noiseless microswimmer in a bulk fluid may be calculated analytically in some simple cases by solving Stokes equations [4][5][6][7][8], analytical methods are less suitable for swimmers in complex environments, in the presence of other swimmers, or when Brownian fluctuations become important.Hence several numerical methods are used nowadays to overcome these limitations.Very succesful tools are coarse-grained hydrodynamic simulation methods such as lattice Boltzmann (LB), dissipative particle dynamics (DPD), and multiparticle collision dynamics (MPCD) -a special version of which is stochastic rotation dynamics (SRD).These methods allow to simulate the time evolution of one or many swimmers, in bulk as well as in confinement.While MPCD and DPD are particle-based methods and naturally include thermal fluctuations, LB is based on the time evolution of smooth velocity distribution functions and noise may be included on top.
We first introduce the basics of the MPCD method in section 2 and discuss the squirmer model briefly in section 3. We then explain in detail how the motion of the squirmer can be coupled to the MPCD fluid using a hybrid molecular dynamics -multiparticle collision dynamics (MD-MPCD) scheme in section 4. We then present results for swimming in a narrow slit geometry in section 5, and finish with a conclusion in section 6.

Multiparticle collision dynamics
Multiparticle collision dynamics (MPCD) is a coarse-grained solver of the Navier-Stokes equations [54,55], which naturally includes thermal fluctuations because it simulates the fluid explicitly via pointlike, effective fluid particles kept at temperature T 0 .Transport coefficients such as the shear viscosity η can be calculated analytically for a wide range of simulation parameters but depend on the particular collision rule between fluid particles [55].
The MPCD method was originally proposed by Malevanets and Kapral and termed stochastic rotation dynamics (SRD) [56,57].The fluid is modeled by N p effective fluid particles, which are pointlike and have mass m 0 .They perform alternating streaming and collision steps and thereby solve the Navier-Stokes equations on a coarsegrained level since the collision step preserves local momentum.Initially, the particles are randomly placed in the simulation domain at initial positions r i (t = 0) and at velocities v i (t = 0) that are drawn from a normal distribution with zero mean and standard deviation σ 0 = k B T 0 /m 0 .
In the streaming step [see Fig. 1(a)] the fluid particles are moved ballistically for a time ∆t such that their positions are updated to Noteworthy, here ∆t is a simulation parameter, which determines fluid properties such as the viscosity.This is in contrast to molecular dynamics (MD) simulations, where ∆t only determines the accuracy of a simulation run.Typically, ∆t = 10 −2 −10 0 in MPCD units.After the streaming step all the fluid particles are sorted into cubic cells with edge length a 0 .To perform the collision step, all particles in a cell interact with each other and their velocities are updated to [see also Fig. 1(b)] Here ξ = 1, . . ., N c ist the cell index, N c is the total number of cells, and u ξ (t) is the average particle velocity in cell ξ.The collision operator Ξ(r j (t), v j (t)) in general depends on positions and velocities of the fluid particles in cell ξ.It has to be chosen such that momentum is conserved locally in each cell, in order to reproduce a flow field, which solves the Navier-Stokes equations.Depending on its concrete realization it either conserves energy or keeps temperature and thus the mean kinetic energy constant.In the original SRD method the fluid particles are part of a microcanonical ensemble and both momentum and energy is conserved locally.The collision operator reads where R α is a rotation operator, which rotates the relative velocities v i − u ξ about a randomly chosen axis by a fixed angle α.In order to fulfill Galilean invariance but also to reduce memory effects and correlations in the collision step, the cell grid is shifted randomly at each time step [58] by a random vector s = {s 1 , s 2 , s 3 }, where the components s i are drawn from a uniform distribution In their original paper Malevanets and Kapral showed [56] that by averaging out the particle dynamics on short length and time scales using a Chapman-Enskog expansion of the corresponding Boltzmann equation, the Navier-Stokes equations can be derived.Furthermore, the H-theorem guarantees fluid relaxation towards equilibrium.The authors also demonstrated that the number of fluid particles per cell is Poisson distributed and the fluid particle velocity distribution converges fast towards a Maxwell-Boltzmann distribution.In addition, transport coefficients like the fluid viscosity can be calculated from the relevant model parameters m 0 , a 0 , T 0 , γ, α and ∆t [56][57][58][59][60][61][62], where γ = Nc ξ=1 n ξ /N c is the average number of fluid particles n ξ per cell.
Noguchi et al. developed alternative collision operators for MPCD fluids.They use a Langevin or Anderson thermostat to keep temperature constant but also to conserve momentum locally [63,64].Thus, a canonical ensemble of fluid particles is modeled, where energy is conserved on average.An advantage of these alternative collision rules is that local angular momentum conservation, which is missing in the original SRD formulation, can be included in a rather simple way.Although angular momentum conservation does not directly alter the Navier Stokes equations, it ensures that the stress tensor is symmetric, as it should be for the Newtonian fluid with its isotropy.Indeed, it has been shown that neglecting angular momentum conservation in the fluid may lead to incorrect physical effects [65].One popular approach for the collision operator is to use the Anderson thermostat and, in addition, conserve angular momentum.This type of collision rule is abbreviated by MPC-AT+a [63].It has widely been used to study the dynamics of squirmers in MPCD fluids [46-50, 52, 53].

Squirmer model
The squirmer is a simple spherical model swimmer.In its simplest form, a single parameter allows to tune the flow field initiated by the squirmer in the surrounding fluid.Thereby it mimics typical types of microswimmers ranging from pushers via neutral swimmer to pullers.The squirmer has been introduced originally as a model for ciliated microorganisms such as the spherical Volvox or the elongated Opalina [5,6].Their cell bodies are covered entirely by synchronously beating cilia, which propel them forward by propagating so-called metachronal waves along the cell surface [66].After averaging over a full beating cycle, characterized by a power stroke and recovery stroke, the resulting net flow close to the surface is nonzero.To account for these flows, the squirmer model was developed by Lighthill [5] and extended by Blake [6].Cilia-induced flows close to the surface of the swimmer is here modeled by an effective slip velocity v s directly at the surface of the sphere.Considering only axisymmetric flows, the timedependent surface velocity fields in spherical coordinates read v s (θ, t) = v θ (θ, t)e θ + v r (θ, t)e r .
Without any thermal noise, a squirmer in bulk simply swims along its symmetry axis e, and θ is the angle between this axis and the point r s on the surface [see Fig. 2(a)].
The surface velocities can be expanded in Legendre polynomials P n , where The expansion coefficients A n and B n are sometimes called surface velocity modes.Solving the governing equations of low-Reynolds-number hydrodynamics, the Stokes equations [67,68], together with the boundary conditions [Eq. ( 5)] gives analytic expressions for the flow fields created by a squirmer in the surrounding fluid.In many cases the surface velocity is assumed to be only tangential to the squirmer surface and time-independent.This describes in a good approximation the hydrodynamic flows close to spherical artificial microswimmers.Examples are active colloids [1,2] and droplets [3,69,70] which move themselves forward realizing static tangential stresses close to their surface, for example induced by selfphoretic mechanisms or self-induced Marangoni flows, respectively.Other examples are biological microswimmers such as Paramecia [71] or Volvox [72] which propel themselves by synchronously beating thousands of tiny cilia attached to their surface.Averaged over a beating cycle, they induce strong tangential flows to move forward.Setting v r (θ, t) = 0, the velocity field v(r) around a squirmer of radius R with orientation e and located at position r 0 can be calculated to [73] v(r) = − 1 3 Here r = |r − r 0 | is the distance from the squirmer center, r = (r − r 0 )/r is the radial unit vector, and The swimming speed v 0 of a noiseless squirmer in bulk solely depends on the first mode B 1 [5], and its velocity is constant, V 0 = v 0 e.It can be shown that the velocity of a sphere with tangential surface velocity v s is equivalent to the average velocity on the surface [74,75], Indeed, it can be shown directly that only the first mode gives a nonzero contribution to the integral resulting in Eq. ( 9).It is interesting to have a closer look at the flow field presented in Eq. (7), in particular at the behavior far away from the squirmer.The flow field created by the first mode B 1 is that of a source dipole.It decays as ∼ r −3 and is often the dominant contribution for spherical microswimmers [1].The flow field created by the second mode B 2 is a force dipole field.It decays as ∼ r −2 and is the dominant field far away from the squirmer, because all other modes B i with i 3 decay faster, see Eq. (7).For B 2 < 0 the far-field is the same as for a so-called pusher (or contractile swimmer) and for B 2 > 0 a so-called puller (or extensile swimmer) is realized.Typical pushers have the propelling apparatus at the back of their cell bodies, as observable in many bacteria or sperm cells.In contrast, pullers like Chlamydomonas have their flagella in front of the cell body.By tuning B 2 , the squirmer can therefore be used as a simple model to capture different types of swimmers: pushers (B 2 < 0), pullers (B 2 > 0), and source-dipole swimmers (B 2 = 0).
In order to capture the flow fields of microswimmers in leading order correctly, it is sufficient to take the first two modes B 1 and B 2 into account.Experiments with active emulsion droplets have indeed demonstrated that the flow field can be approxiamted by these two modes only [76].The flow field of a swimming Paramecium has been measured as well.It can also be approximated well by considering only the leading modes with B 1 being the dominant mode [71].

Details of MPCD simulations of squirmers
In this section we discuss in detail how the dynamics of squirmers moving in the MPCD fluid is implemented.We also explain how the dynamics of multiple squirmers is treated, how to include solid boundaries, and how pressuredriven channel flow is realized in the simulation.
In the streaming step both squirmers and fluid particles are moved forward.In addition, squirmers collide with each other and squirmers and fluid particles hit bounding walls.In the collision step virtual particles contribute when fluid particles exchange momentum.They are inside squirmers and walls that overlap with collision cells.In both streaming and collision step momentum and angular momentum are transferred between fluid, squirmers, and walls.In the following we describe in detail the procedure how to perform the actual simulations.
Typical simulation geometries in 3D are bulk systems implemented by periodic boundary conditions in x, y, and z direction, or a fluid confined between two infinitely extended parallel walls separated by a distance H.For the latter case the number of collision cells in the direction perpendicular to the walls is H/a 0 + 1.This ensures that despite the random grid shifts all fluid particles and a sufficiently thick layer of virtual particles are covered by the grid all the time.Initially, the effective fluid particles are randomly distributed in the simulation domain but are not allowed to overlap with the squirmers.The velocities are normal distributed with zero mean and standard deviation σ 0 = k B T 0 /m 0 .The mass density of the fluid is given by ρ f = γm 0 /a 3 0 and the number density by n f = γ/a 3 0 with γ the average number of fluid particles per cell.One possibility is to choose cell length a 0 , mass m 0 , and energy k B T 0 as the basic units [128].Then the unit of time becomes τ 0 = a 0 m 0 /k B T 0 and the unit of viscosity is η 0 = √ m 0 k B T 0 /a 2 0 .The remaining parameters are ∆t, γ, and, in addition, the rotation angle α, when the SRD collision operator is used.These parameters deter-mine transport coefficients of the fluid such as the viscosity [55].
In our simulations we consider the dynamics of one or N S identical squirmers of radius R, with surface velocity modes B 1 and B 2 , and with the same mass density as the fluid, ρ S = ρ f .Note, when B 1 = B 2 = 0 passive colloids can be simulated as well.Squirmers are initially placed such that they do not overlap with other squirmers, walls, or fluid particles.The concrete initial conditions for their positions and orientations depend on the specific problem of interest.Usually, initial values for velocities and angular velocities are set to V i = 0 and Ω i = 0.The mass M S of a squirmer is given by with volume V S = 4πR 3 /3.The moment of inertia tensor of a squirmer is isotropic, I S = I S 1, with

Streaming step
In the streaming step squirmers and effective fluid particles move ballistically and they interact with each other and with confining walls.There are different ways how to integrate the streaming of the system.One way is to use interaction potentials, such as a cut-off Lennard-Jones or WCA potential [129], between squirmers to ensure they do not overlap.Then, the streaming step is divided into N MD time steps, which defines the MD time interval δt = ∆t/N MD .It has to be chosen sufficiently small in order to correctly integrate motion in possibly very steep potentials.Another way is to perform event-driven MD simulations during the time intervall δt = ∆t.
Motion of the squirmers Squirmers are moved ballistically during time δt in order to update their positions R i and orientations e i .In some situations external forces F i (t) and/or torques T i (t) act on the squirmers.They can also include the interaction forces between squirmers.Then, one possibility to move squirmers in the streaming step is to use a simple Velocity-Verlet integration [130,131].It updates positions and orientations according to and translational and angular velocities according to When hard potentials between squirmers are used, an alternative is to perform event-driven simulations.Squirmers may overlap when they move for a time δt.Then, all collision times t ij between squirmer i and another squirmer j or a wall j have to be calculated.Instead of moving during time δt, the squirmers can only move for a time t m < δt with t m = min(t ij ).If squirmer i hits a wall and the wall is assumed to be smooth so that it does not act with a frictional torque on the squirmer, the squirmer velocity is simply updated to with n the wall normal.If two squirmers i and j are the collision partners, which have smooth surfaces so that they do not exchange angular momentum, their velocities are updated to [132] with and After the collision is performed, all squirmers are moved during time interval δt − t m and again the collision partners and times are determined and their velocities updated.This procedure is repeated until the total time δt consumed.
Motion of the fluid particles In the absence of external forces acting on the fluid and fluid-squirmer interactions, the streaming of the fluid particles is simply given by Eq. ( 1).However, to implement a Poiseuille flow in a channel geometry or between two infinitely extended parallel plates, one needs to apply a constant force f i to them representing the constant pressure gradient [133].Then, again the Velocity-Verlet algorithm can be used to update positions and velocities, Fluid particles also collide with squirmers or bounding walls during the streaming step.In contrast to the eventdriven implementation of squirmer dynamics, one may not want to calculate exact collision times between fluid particles and squirmers or walls since this slows down the simulation.One way to overcome this is to first stream all the squirmers during time δt neglecting that some fluid particles may enter them and only then move the fluid particles.Padding et al. [134] have proposed an efficient method to do so for simulating colloidal suspensions, which has been used for squirmer-fluid interactions as well [45,[47][48][49].We start with a fluid particle hitting a wall.Once this happens, it is moved half a time step back and its velocity is updated as follows.In order to fulfill the no-slip boundary condition at the wall, the bounce back rule is applied, where the velocity is simply reversed, meaning that surfaces are rough for the fluid.Then, the fluid particle is moved forward for half a time step with the new velocity.A similar procedure is applied if fluid particle i hits squirmer j, but one has to account for the fact that the surface of the squirmer carries a slip velocity field as well [45], ) Here, r * ij is the collision point of the particle on the surface of the squirmer and v s (e j , r * ij ) is the slip velocity [Eq.( 11)] of squirmer j located at r * ij .When passive colloids are simulated, v s is set to zero.Note that during time interval δt a fluid particle is allowed to interact several times with different squirmers or bounding walls.These multiple reflections help to reduce possible depletion interactions between squirmers.
When a fluid particle interacts with a squirmer or a wall, its momentum p i = m 0 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 squirmer j is simply which is then transferred to the squirmer.In general, if N * j fluid particles hit squirmer j during time interval δt, its velocity V j and angular velocity Ω j are updated to where are the total linear and angular momentum transferred to squirmer j in the streaming step.

Collision step
In the collision step fluid particles interact with each other but also interact with squirmers and bounding walls via virtual particles [135].The collision step starts with shifting the cell grid by a random vector s as discussed in Sec. 2.Then, collision cells partly overlapping with squirmers and walls are filled with virtual particles, which are placed in the walls and in the squirmers, see Fig. 3.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 γ and hence a smaller viscosity [55].The virtual particles are placed randomly distributed at density ρ f in a layer of thickness d W = a 0 inside a bounding wall and in layers of thickness d S = √ 3a 0 inside the squirmers, see Fig. 3.It is always guaranteed that grid cells overlapping with walls and squirmers are always completely filled.The positions and velocities of the virtual particles are denoted by ri 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 0 /m 0 , and zero mean.In addition, virtual particle i located in squirmer j also assumes the local slip velocity of the point r * ij on the squirmer surface, which is closest to ri , plus the velocity of this point due to the translation and rotation of the squirmer.So, one sets [46] vi = vr where vr i are the random velocities.Then, as usual, fluid and virtual particles are sorted into the cells.We denote the number of fluid and virtual particles in cell ξ by N f ξ and N v ξ , respectively.The total number of particles in the cell is In order to perform now the collision step, the mean velocity in a cell is computed, and the new velocities of the particles in a cell are set by the mean velocity u ξ plus a term given by the collision operator, see Eq. ( 2).When the Anderson thermostat is used, random velocities have to be computed, which are denoted by v r i and vr i , and which are again drawn from a Gaussian distribution with width σ 0 = k B T 0 /m 0 .Then, in order to conserve linear momentum the change of total velocity due to the added random velocities, V ξ , has to be subtracted.If conservation of angular momentum is required, by using the collision operator MPC-AT+a, a further term needs to be added.To implement the collision operator, first the center of mass has to be calculated, and the relative positions r s i = r i − r s ξ and rs i = ri − r s ξ .They are needed to calculate the inverse of the moment of inertia tensor I −1 ξ in each cell, where (33) The random velocities added in the collision step change angular momentum by (34) To compensate for ∆L ξ and conserve angular momentum in each cell, the following angular velocity is computed [63], ω and used to rotate the particle velocities in a cell.Thus, by adding the extra terms ω AMC ξ × r s i or ω AMC ξ × rs i to the new particle velocities, angular momentum is conserved without changing linear momentum.To summarize, in the collision step the particle velocities are updated according to [63] Finally, momentum and angular momentum is transferred to the squirmers, no matter which collision operator has been used.In each collision cell momentum (and angular momentum) is conserved locally but exchanged between fluid particles and virtual particles.The change of momentum for a virtual particle is then The momenta and angular momenta of all virtual particles i located in squirmer j are then added up and assigned to squirmer j.Thus, after a collision step each squirmer assumes an additional momentum ∆P c j and angular momentum ∆L c j , where the sum goes over all virtual particles located in squirmer j with total number N v j .So, after each collision step the squirmer velocities and angular velocities are updated according to where I S is the moment of inertia of the squirmer.

Dynamics of a squirmer in a narrow slit geometry
Finally, we apply the discussed simulation procedure to the case of a single squirmer moving in a narrow-slit or Hele-Shaw geometry.We use a squirmer of radius R = 3a 0 , swimming speed determined by B 1 = 0.1a 0 /τ 0 , and choose different values of β = B 2 /B 1 .We place planar walls at positions x = ±4a 0 such that there is only a small gap between squirmer surface and walls.Then, we perform several simulations at various β and initial conditions.All the recorded trajectories are shown in Fig. 4(a) where the position x is plotted versus the angle θ, which is the orientation of the squirmer relative to the upper wall normal.For strongly negative β values (strong pushers) the squirmer tends to orient perpendicular to the walls (shown for β = −3 in Fig. 4(a)).This is also sketched in Fig. 4(b).Increasing β above a certain value induces oscillations of the swimmer between the walls: The squirmer hits a surface, spends some time there, reorients away from the surface, hits the second surface and so on (shown for β = −1.5 in Fig. 4(a) and sketched in Fig. 4(c)).
neutral squirmer still performs oscillations but mainly avoids touching the walls by swimming with an orientation performing small oscillations about θ = π/2 along the bounding walls, see also Fig. 4(d).In contrast, pullers (β > 0) swim stable between the two walls, as sketched in Fig. 4(e).Interestingly, the β-dependent dynamics can be captured by considering the lubrication approximation of a squirmer in front of a wall [49,73], which predicts the correct dynamics for all considered β values [48].
We also measured the swimming velocity and the flow field around a squirmer with radius R = 3a 0 and surface velocity modes B 1 = 0.1a 0 /τ 0 and β = 3, which swims stable between the plates, by averaging over many realisations and over time.The swimming velocity results in v MPCD 0 = 0.0657a 0 /τ 0 , which is very close to the analytic expression for the swimming speed in bulk v 0 = 0.0667 [Eq.( 9)].To obtain the flow fields, the simulation box around the swimmer is divided into small cubic boxes, for which we determine the flow velocity vector.The boxes do not have to coincide with collision cell boxes and may even have higher resolution.Figure 5 shows the flow field (blue arrows) and stream lines (red) for a slice of the field parallel to the walls (top) and perpendicular to the walls (bottom) in the lab fame.Due to the presence of the walls the force dipole fields develop closed loops [136,137].We use a grid size with edge length a 0 /2 to monitor the velocity field, meaning that MPCD is able to resolve flow fields below the the cell size a 0 of the collision grid.
In Fig. 6 we show the decay of the flow field perpendicular (v ⊥ ) and parallel (v || ) to the swimming direction, which we compare to the decay in bulk using the analytic expressions [Eq.(13)].Due to the leading force dipole field, they decay with ∼ r −2 .For intermediate distances away from the squirmer the flow fields obtained from our simulations in strong confinement show a similar decay.Indeed it has been shown analytically that for very strong con-finement flow fields of microswimmers decay in general with ∼ r −2 [50,138].However, in our simulations the decay of the flow fields far away from the squirmer is faster, which probably stems from the periodic boundary conditions [38,40].Also, the fact that the strength of the flow field in front of the swimmer is weaker than at the back may be a signature of the periodic boundary conditions or due to the fact that MPCD does not model a perfect incompressible Stokes flow, but finite Reynolds and Schmidt numbers which account for finite fluid inertia and compressibility, respectively [55,128].

Conclusion
Here we have reviewed and discussed how simple spherical model swimmers -squirmers -can be modeled and simulated using the method of multiparticle collision dynamics.We have shown how the dynamics of the fluid and the squirmers is coupled both in the streaming step and in the collision step, where a correct momentum and angular momentum transfer between particles ensures correct hydrodynamics [45] and thermal fluctuations [46].In order to enhance the accuracy of the flow fields around the squirmer, the use of virtual particles inside squirmers, which contribute to the collision step, is necessary.Also bounding surfaces and pressure-driven flows can be simulated straightforward.Since the method models the fluid explicitely via effective fluid particles, thermal noise is included naturally in the simulation.In order to obtain smooth flow fields averaging over several realisations and over time is necessary.

Figure 1 :
Figure 1: Sketch of the MPCD method.(a) Streaming step: Pointlike fluid particles move ballistically for a time ∆t with velocities v i .(b) Collision step: Particles are sorted into cubic cells of edge length a 0 and exchange momentum with other particles in the cell such that the total momentum is conserved.

Figure 3 :
Figure 3: 2D sketch of collision step.Squirmers are shown in red and a wall in grey.Small dots in the white region indicate the fluid particles and dots in the squirmers and in the wall are virtual particles, which are located in a thin layer of thickness d W = a 0 in the wall and in a thin shell of thickness d S = √ 3a 0 in the squirmers.The dashed lines show the collision cell grid.

Figure 4 :
Figure 4: (a) Typical phase portraits from MPCD simulations for a squirmer moving in a narrow-slit geometry.Walls are located at x = ±4a 0 and the radius of the squirmer is R = 3a 0 .Thus, minimum and maximum positions of its center are at x = ±a 0 .The angle θ indicates the orientation of the squirmer relative to the wall normal.So, θ = 0 (θ = π/2) is the orientation perpendicular (parallel) to the walls.The long-time dynamics is shown in green.(b-e): Sketch of typical squirmer trajectories for different β.The initial positions are shown in blue, the final positions in green, and intermediate positions in orange.The orientation vectors are sketched by black arrows.

Figure 5 :
Figure 5: Flow field for a puller swimming in a narrowslit geometry.The velocity field is indicated in blue and streamlines in red.Top: Top view of the flow field.Bottom: side view.

Figure 6 :
Figure 6: (a) Decay of the flow field perpendicular (v ⊥ ) to the swimmer direction.The flow field in the −z direction is shown in blue, and in the +z direction in orange.(b) Decay of the flow field parallel (v || ) to the swimmer direction.The flow field in the +y direction (front of the swimmer) is shown in blue, and in the −y direction (back of the swimmer) in orange.Solid lines are curves obtained from the analytic model of a squirmer in bulk [Eq.(13)].