Simulating Microswimmers Under Confinement With Dissipative Particle (Hydro) Dynamics

In this work we study microwimmers, whether colloids or polymers, embedded in bulk or in confinement. We explicitly consider hydrodynamic interactions and simulate the swimmers via an implementation inspired by the squirmer model. Concerning the surrounding fluid, we employ a Dissipative Particle Dynamics scheme. Differently from the Lattice-Boltzmann technique, on the one side this approach allows us to properly deal not only with hydrodynamics but also with thermal fluctuations. On the other side, this approach enables us to study microwimmers with complex shapes, ranging from spherical colloids to polymers. To start with, we study a simple spherical colloid. We analyze the features of the velocity fields of the surrounding solvent, when the colloid is a pusher, a puller or a neutral swimmer either in bulk or confined in a cylindrical channel. Next, we characterise its dynamical behaviour by computing the mean square displacement and the long time diffusion when the active colloid is in bulk or in a channel (varying its radius) and analyze the orientation autocorrelation function in the latter case. While the three studied squirmer types are characterised by the same bulk diffusion, the cylindrical confinement considerably modulates the diffusion and the orientation autocorrelation function. Finally, we focus our attention on a more complex shape: an active polymer. We first characterise the structural features computing its radius of gyration when in bulk or in cylindrical confinement, and compare to known results obtained without hydrodynamics. Next, we characterise the dynamical behaviour of the active polymer by computing its mean square displacement and the long time diffusion. On the one hand, both diffusion and radius of gyration decrease due to the hydrodynamic interaction when the system is in bulk. On the other hand, the effect of confinement is to decrease the radius of gyration, disturbing the motion of the polymer and thus reducing its diffusion.

1 Introduction

Active Matter is a branch of Physics that focuses on the study of intrinsically out-of-equilibrium systems due to energy being constantly supplied, converted into directed motion and dissipated by individual constituents. Active Matter is a field that has raised a lot of interest in the last decade, since it captures complex collective behaviours, often exclusively associated to living matter, and might enable a wide range of technological applications [1]. One of the paradigmatic systems of Active Matter consists of a suspension of active particles. Active particles can be living (such as bacteria) or synthetic (such as active colloids). Active colloids are micron-size particles which self-propel through a medium by converting energy extracted from their environment into directed motion [2, 3], with potential medical and technological applications [4–10]. The collective behaviour of systems constituted by a large number of these particles is rich and complex as shown by a series of recent numerical [11–17] and experimental [18] works, and in many cases cannot be ascribed solely to the particles motion since hydrodynamics due to the surrounding solvent might need to be taken into account [19]. This is the case for microswimmers [20], whose motion is an essential aspect of life.

Microswimmers are usually ciliated and/or flagellated microorganisms that achieve propulsion thanks to the movement of their cilia located on their outer surface: for this reason one can consider them as self-propelled microorganisms. In the last few years microswimmers have been intensively studied, being of interest in several interdisciplinary sciences. Examples of living microswimmers are Escherichia coli bacterium, Paramecium or sperm cells, or algae (such as Chlamidomonas or Volvox). Whereas examples of synthetic microswimmers are Janus colloidal particles. When considering the effect of hydrodynamic interactions, numerical studies of a two dimensional suspension of self-propelled repulsive swimmers have demonstrated that hydrodynamics affects, not only the phase behaviour of a dense suspension [21], as suggested by Ishikawa [22] in an early work, but also the dynamics of transient clusters at lower densities [23]. Moreover, other theoretical, numerical and experimental results have also revealed the importance of hydrodynamics in these systems [24–27].

To model microswimmers, Blake and Lighthill proposed the so called squirmer model [28, 29]. The squirmer model reproduces the induced hydrodynamic flow around a spherical swimmer while preserving the main features of the active stresses generated by it [30]. The spherical squirmer particle mimics the effect of the cilia on the fluid as a prescribed slip velocity tangential to the surface. The described mechanism is the one that leads to the swimmer’s propulsion. A squirmer is characterized by two modes accounting for its swimming velocity and its active stress. Depending on the active stress, it is possible to classify squirmers as pushers (e.g., E. coli, sperm), pullers (e.g., Chlamydomonas) and neutral (e.g., Paramecium) swimmers [31, 32]. The squirmer model has been expanded for complex swimmers, such as non-spherical swimmers [32] and explicitly ciliated microorganisms [31].

Besides mimicking the swimmer’s behaviour, it is important to choose a model to mimic the features of the surrounding fluid. The applicability of atomistic algorithms (Molecular Dynamics-like) to simulate the fluid is limited, since they only allow to study short time and length scales (few hundreds of nanoseconds and few tens of nanometers). To explore longer length/time scales, more relevant for living swimmers, atomistic methods become computationally inefficient. Thus, one might consider mesoscopic methods, that bridge the gap between the microscopic and the macroscopic continuum scale [33]. These methods span longer length and time scales: from several nanometers to micrometers and from nanoseconds to microseconds. The most renown mesocopic numerical models used to simulate fluids that fully consider hydrodynamic interactions are Lattice-Boltzmann [34], Multiparticle Collision Dynamics [35–37] and Dissipative Particle Dynamics [38]. The Lattice Boltzmann (LB) approach consists in describing the solvent in terms of the density of particles with a given velocity at a node of a given lattice. The discretized velocities join the nodes and prescribe the lattice connectivity [39]. The LB model reproduces the dynamics of a Newtonian liquid of a given shear viscosity η. Relevant hydrodynamic variables are recovered as moments of the one-particle velocity distribution functions. The total force and torque the fluid exerts on a particle embedded in it are obtained by imposing that the total momentum exchange between the particle and the fluid nodes vanishes. Since a Lattice Boltzmann code is computationally expensive, from a practical point of view it is possible to parallelize it using Message Passage Interface to exploit the excellent scalability of LB on supercomputing facilities [40]. In the Multiparticle Collision Dynamics (MPCD) approach [35–37] a fluid is represented by N point particles with continuous positions and velocities. The particle dynamics proceeds in two steps: streaming and collision. During the streaming step, particles move ballistically. Whereas in the collision step particles interact locally via an instantaneous stochastic process, that could be based on stochastic rotation dynamics with angular momentum conservation [32]. For this purpose, the simulation box is partitioned into cubic collision cells. Within MPCD Galilean invariance is ensured, together with thermal fluctuations. The algorithm conserves mass, linear, and angular momentum on the collision cell level, which gives rise to hydrodynamics on large length and long time scales. Dissipative particle dynamics (DPD) is one of the most efficient mesoscale coarse-grained approaches for modeling soft matter systems. DPD was originally proposed by Hoogerbrugge and Koelmann [38] as an off-lattice, momentum conserving, Galilean invariant mesoscopic method, the coarse-grained dynamics of which obeys the Navier-Stokes equations and preserve hydrodynamics. Later on, Espanol and Warren [41] reformulated the DPD model in terms of stochastic differential equations. DPD consists in modified Langevin equations that operate between pairs of particles interacting via three different forces: conservative, dissipative and random (thermal) forces. The DPD model has already been used to model complex colloidal suspensions, such as proteins [42] or red globules in blood [43].

Along with hydrodynamics, confinement also plays a major impact on the dynamics of microswimmers. Their interaction with bounding walls is different depending on the type of microswimmer we are dealing with and can lead to different transport and aggregation phenomena. As of today, this facts have been studied theoretically [44–46], numerically [47–49] and experimentally [31, 50]. While in these works the study focuses mainly in the interaction of microswimmers with plane boundaries, other types of confinement also display relevant features, for instance the effect of porous media in bacterial suspensions has also been reported [51]. Finally, of course it is worth noting that the effect of confinement is not only limited to active matter systems but also plays a role in a wide range of systems, such as passive hard spheres [52].

In the present work we propose to model suspensions of microwimmers with DPD hydrodynamics inspired by the squirmer model. When the agent is a sphere we choose a raspberry-like structure [53–55] and we will directly consider the squirmer model. Whereas when the agent is a polymer, we will build the polymer as a chain of monomers, and treat each monomer similarly, but not rigorously, as a squirmer. To properly deal with hydrodynamics, we will mimic the surrounding fluid via DPD interactions, using an in-house extension of the LAMMPS [56] open source package implementing appropriate reaction forces on the swimmer’s particles that balance the forces exerted on the fluid and enable its propulsion [57, 58]. Our choice is motivated by the fact that differently from LB [34], DPD easily allows to take into account thermal fluctuations and to simulate colloids with complex shapes (not only spherical). Moreover, it is also easier to control compressibility and Schmidt number in DPD than MCPD. Hence, in DPD it is easier to control the appropriate dynamic regime that couples the solvent and solute dynamics. Although this is not the focus of his paper, DPD also allows for a more thorough control of the phase diagram of the solvent and how to deal with fluid phase coexistence. Firstly we study the dynamical behaviour of either microwimmer in bulk. In the case of active colloids, we establish the flow fields surrounding the particle and compute their diffusion, comparing pushers, pullers and neutral swimmers. In the case of active polymers, besides the dynamics we also study its conformational features. Next, we confine either microwimmer in a cylindrical channel, and unravel the effect of hydrodynamics as compare to the equivalent systems where hydrodynamics is not present. For each system we explore different Reynolds and Péclet numbers. The Reynolds number [59] is the ratio of inertial to viscous forces within a fluid subjected to relative internal motion: this number measures the amount of turbulence of the solvent in the system. The Péclet number [11, 60] is defined as the ratio of the rate of advection of a physical quantity by the flow to the rate of diffusion of the same quantity. This number quantifies the degree of activity of microwimmers.

The manuscript is organised as follows. In Section 2 we describe the relevant physical quantities and the technical details of the implementation. We first describe the DPD method to simulate the solvent (Section 2.1), implemented within the LAMMPS open source numerical package [56]. Next, we present the two microwimmers under study: the active colloid (Section 3.1) and the active polymer (Section 3.2). In Section 3.1, we introduce the raspberry-like active colloid (Figure 1A) in bulk and when interacting with a cylindrical surface (Figure 1C). In Section 3.2, we report the active polymer (Figure 1B), as in ref. [11], in bulk and under cylindrical confinement (Figure 1D). The way we implemented hydrodynamics is reported in Section 2.4, being the same for both active objects embedded in a DPD solvent. In the same section we characterize the physical quantities of a fluid such as the kinematic viscosity (ν) and the solvent diffusion coefficient (D_sol), and parameters to quantify the activity of the colloid/polymer embedded in a fluid, such as the Reynolds number and the Péclet number. Finally, in Section 2.5 we report the analysis tools used to study the microwimmers in bulk or under confinement. In Section 3 we present the results obtained, first for the colloid (Section 3.1) and then for the polymer (Section 3.2). In Section 4 we discuss the results and comment on future avenues.

FIGURE 1

FIGURE 1. (A): A raspberry-like active colloid composed of 18 filler particles and one thruster particles at the center. Note that in these figures the filler particle radius is scaled down to 1 for better visibility, but in all simulations we use $r_c^{\text{sf}}=2$ as the solvent-filler DPD cutoff in order to obtain a more spherical colloid and to prevent the solvent particles from stepping into the colloid. (B): An active polymer composed of 20 beads. (C): An active colloids in a cylindrical confinement. (D): An active polymer in a cylindrical confinement. Red: thruster particles. Gray: filler particles. Blue: wall particles composing the confining channel. Yellow: solvent particles.

2 Materials and Methods

In this work we study an active colloid and an active polymer embedded in a fluid solvent either in bulk or confined inside a cylindrical channel. We simulate the active colloid as a spherically-shaped collection of particles merged together by rigid interactions. Whereas the active polymer is built as a chain of monomers glued together by harmonic interactions that enable their relative movement. The rest of the interactions are the DPD-like interactions between any two particles, the hydrodynamic force-field that enables the agents’ propulsion and, in the case of the confined polymer, a repulsive (WCA-like) potential between channel (particles) and polymer/solvent particles.

2.1 Modeling the Solvent With Dissipative Particle Dynamics

Our system consists of microwimmers embedded in a solvent, where hydrodynamics is explicitly taken into account. The fluid surrounding the microwimmer is simulated as a collection of individual particles interacting via Dissipative Particle Dynamics [38]. According to DPD, below a given cutoff r_c, the force acting on the i-th solvent particle consists of three contributions,

${\vec{F}}_i=\sum _j\left({\vec{F}}_{ij}^C+{\vec{F}}_{ij}^D+{\vec{F}}_{ij}^R\right){\widehat{r}}_{ij}\text{if}r<r_c,(1)$

being ${\widehat{r}}_{ij}=({\vec{r}}_i-{\vec{r}}_j)/r_{ij}$ the inter-particle unitary direction between the i-th and j-th particles. The conservative term is ${\vec{F}}_{ij}^C=Aw(r_{ij})$, where A is the amplitude and w (r_ij) = 1 − r_ij/r_c a weighting factor varying between 0 and 1 as in ref. [33]. The dissipative contribution reads ${\vec{F}}_{ij}^D=-\gamma w^2(r_{ij})\left({\widehat{r}}_{ij}\cdot {\vec{v}}_{ij}\right)$ with friction coefficient γ. Finally, the thermal contribution ${\vec{F}}_{ij}^R=\sigma \alpha w(r_{ij})/\sqrt{\mathrm{\Delta }t}$ is a random force, where α, is a Gaussian random number with zero mean and unit variance, Δt the chosen time-step for the time integration and $\sigma =\sqrt{2k_{B}T\gamma }$ is related to the mean of the random force via fluctuation-dissipation, being T the temperature of the system.

In the current work, we implement the DPD solvent via the LAMMPS open source package [56], setting the time step to Δt/τ = 10^–2 for the simulations of the active colloid and Δt/τ = 5 ⋅ 10^–3 for the simulations of the active polymer. In both systems, we choose an equilibration time of $\sim 10^4$ steps, while the production run is of the order of 10⁶ steps. The number of solvent particles for the system contaning the active colloid in bulk is N = 10,125, distributed in a cubic simulation box of L = 15. In cylindrical confinement, depending on the channel radius R_cyl = {3.5, 4.5, 5.5, 6.5} the number of solvent particles is N_sol = {3464, 5726, 8553, 11,946}, respectively, and the channel length is fixed to L_cyl = 30. The number of solvent particles for the polymer system in bulk is around N_sol = 24,000 distributed in a cubic simulation box of L = 20. In the polymer confined case, with channel radius R_cyl = 6 and length L = 50, the number of solvent particles is around N_sol = 17,000. For all simulations, the mass of the solvent particles is fixed to m = 1 and the numerical density to ρ_sol = 3. The characteristic length scale for all our simulations is the DPD cutoff distance between solvent particles $r_c\equiv r_c^{ss}=1$, the mass scale is fixed by the mass of one solvent particle m = 1, and for the time scale we fix τ = 1. Following ref. [33], the DPD interaction parameters between solvent-solvent particles are set to A^ss = 25.0, γ^ss = 4.5 and $r_c^{\text{ss}}=1.0$ (see Tables 1, 2 in the following subsections). The physical properties of a DPD fluid depend on its viscosity [33] that can be computed from the Green-Kubo relation [61] for the stress autocorrelation function (zero-shear viscosity). Later on we will discuss our choice for the fluid’s viscosity. Whereas the DPD parameters used for each microwimmer are reported in their corresponding sections.

TABLE 1

TABLE 1. DPD parameters used to study a spherical colloidal squirmer for (COL) embedded in a solvent (SOL), in bulk or in a cylindrical confinement (CYL). All parameters are in reduced DPD units.

TABLE 2

TABLE 2. DPD parameters used to study an active polymer (POL) embedded in a solvent (SOL), in bulk or in a cylindrical confinement (CYL). All parameters are in reduced DPD units.

It is worth noting that when confining the DPD fluid in a cylinder we observed concentric ring-shaped density fluctuations in the vicinity of the wall at T = 0.1. These fluctuation have been previously observed and can lead to undesirable effects [62]. For this reason we chose T = 1 for all our confined simulations in which these fluctuations were not observed.

2.2 Colloids in Bulk and in Confinement

To study a spherical squirmer, we build a raspberry-like [53–55] colloid made of 19 particles rigidly bonded. In Figure 1A, we represent the active colloids, consisting of one particle (the thruster particle) located at the center of the sphere and the remaining 18 (filler particles) evenly distributed on the surface of a sphere of radius R_col around the center particle. The reason for choosing this structure has been guided by simplicity, balancing the number of particles and sphericity, and is inspired by previously proposed models for complex colloids [53–55]. The reason to consider only one thruster particle is because such an approach is enough to generate activity with a minimal disturbance on the geometrical properties of the swimmer. The shell of passive particles is necessary to control the dimension and shape of the colloid, rendering it an extended body that has an orientation. The propulsion mechanism of the thruster particle will be explained in the next section, when detailing the implementation of the hydrodynamic interactions.

The orientation of the colloid is defined by the “active axis” identified by three chosen co-linear particles. This axis is also the symmetry axis of the force field we will apply to the solvent, and thus will define the colloid’s direction of propulsion. All particles belonging to each colloid interact via DPD: 1) with the solvent, 2) with particles belonging to other colloids and 3) with particles building the channel. However, particles belonging to each colloid do not interact between them (so their overlap does not cause any trouble), except for the rigid interactions that keep them glued together. In Table 1 we report the chosen DPD parameters for all interactions between particles: solvent-solvent, solvent-colloid, solvent-cylinder, colloid-colloid, colloid-cylinder.

In Section 2.4 we will describe different squirmer models, such as pushers, pullers and neutral swimmers, each one characterised by a different velocity field in the surrounding fluid. In order to check whether the raspberry-like colloid reproduces the features of the different squirmers, we compute the velocity fields and compared them to those reported for the different squirmers in ref. [49].

Having studied the active colloid in bulk, we study its physical behaviour when confined in cylindrical environments of different radii. The cylinder is composed of DPD overlapping particles, properly aligned along the x axis at given angles. Overlapped DPD particles are left out of the time integration and their DPD interactions are switched off thus they can be used to model a wall. Particles are first evenly distributed along a circumference in the yz-plane and then this circumference is repeated through the x-axis. The separation of the particles is chosen so that the roughness of the inner surface of the cylinder is the same along the angular and longitudinal directions. Periodic boundary condition (PBC) are applied along the longitudinal direction (x axis). Particles’ interaction parameters are reported in table 1. Choosing DPD interactions for modelling the collisions with the channel allows us to maintain a large time step Δt = 10^–2. Due to the softness of the DPD interactions, we have appropriately set the DPD parameters for the channel particles to avoid leaking of solvent particles through the channel wall. Moreover, DPD enables adding a friction between the solvent and the channel wall. In our case, we have tested that for high enough values of γ we are able to simulate Poiseuille flow. However, for our study we have decided to explore low values of γ, which correspond to the implementation of slip boundary conditions at the channel’s surface.

2.3 Polymers in Bulk and in Confinement

Following ref. [11], we model the active polymer as a chain of active monomers. As shown in Figure 1B, each of the N_b monomers is composed by a single thruster DPD particle, except the head and tail monomer. Since the first (last) particle of the polymer does not have previous (posterior) neighbors, no force is applied on them. Alternatively, one could consider that the activity direction is extrapolated from the neighbor monomer in the chain. However, this alternative approach will not affect the main results and features described in the manuscript. Previously proposed models for active polymers have also taken this approach [11]. Monomers are held together to their first neighbours via a harmonic potential $V_{\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{c}}(r)=K{(r-r_0)}^2$, acting between thruster particles of the connected beads separated by a distance r, with $K=30k_{B}T/r_c^2$, being r₀ ≡ 1.5 r_c. Since all interaction between particles are soft (DPD-like), we can choose dt = 10^–2 as the time step to integrate the equations of motion. As in ref. [11], we assume that all monomers are active apart from the first and the last (in grey in Figure 1). An active force F_a,i acts on each thruster monomer at r_i. The force is characterised by a constant magnitude F_a and a direction of r_i+1 − r_i−1 parallel to the polymer backbone tangent, being r_i+1 and r_i−1 the position vectors of the thruster particles first neighboring monomers.

To characterise the bulk properties of an active polymer, we study a dilute system of 4 active polymers in a box with edge L = 20 at a solvent density of ρ = 3. Care must be taken if the volume fraction of polymers is not low enough, since polymers might interact between each other via hydrodynamics. In our case we avoid this by working with a polymer volume fraction that is always lower than 5%.

To study the effects of confinement we embed the active polymer and the solvent in a cylindrical channel with periodic boundary conditions along the axial axis. The cylinder consists of N_c = 24,415 frozen WCA-like particles that interact with the DPD particles (solvent and polymers) via a WCA-like potential

$V_{\text{LJ}}\left(r\right)=\left\{\begin{array}{cc}4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]+ϵ;\hfill & \text{for}r<2^{1/6}\sigma ,\hfill \\ 0;\hfill & \text{for}r\ge 2^{1/6}\sigma ,\hfill \end{array}\right.(2)$

where ϵ is the unit of energy and σ represent the channel’s particle diameter set to σ = r_c = 1. In all simulations we set k_BT = 1.0 (Lennard-Jones units). Cylinder particles are located close enough to avoid DPD solvent particles to cross the cylinder’s wall.

The chosen values for the DPD parameters are reported in Table 2 for all interactions between particles: solvent-solvent, solvent-polymer, solvent-cylinder, polymer-polymer, polymer-cylinder.

We should stress the fact that when dealing with active colloids or active polymers we have chosen to simulate the cylindrical channel in a different way. In the former case, the channel has been simulated by means of particles interacting via DPD, as explained earlier. Whereas in the latter case, the channel has been built using particles interacting via a repulsive WCA potential, to compare with ref. [11].

2.4 Swimming Induced by Hydrodynamics

In order to numerically consider full hydrodynamic interactions between the microwimmers and the surrounding solvent, we prescribe a force field for the solvent particles surrounding the thruster particles (the red particles in Figure 1).

When dealing with a spherical squirmer, the usual approach consists in prescribing tangential velocities to the solvent particles at the swimmers surface [28]. Note that we have not followed the usual squirmer approach. In our case, tangential solvent forces, instead of velocities, are prescribed over a hydrodynamic active volume Γ_H around the colloid, instead of just at the colloid’s surface. This approach is more general since it enables the possibility of studying different agent shapes and inertial effects, which are present in many active systems [63].

In this study we only consider axisymmetric force fields (Eq. 3). We choose the hydrodynamic region Γ_H as a spherical shell around the thruster particles of inner and outer radii R_c and R_H, respectively. The expressions of the force fields considered for the colloid and the polymer are reported in what follows (see the Supplementary Appendix for more details). The general expression for an axisymmetric force field that vanishes everywhere except inside the aforementioned spherical shell is given by,

$\mathit{f}\left(r,\theta \right)=\left[f_r\left(r,\theta \right){\widehat{\mathit{e}}}_r+f_{\theta }\left(r,\theta \right){\widehat{\mathit{e}}}_{\theta }\right]P_{R_c,R_H}\left(r\right)(3)$

where r is the distance from the thruster to the solvent particle, θ the angle between the colloid’s orientation vector $\widehat{\mathit{e}}$ and the solvent position vector, ${\widehat{\mathit{e}}}_r$ and ${\widehat{\mathit{e}}}_{\theta }$ are the radial and tangential unitary vectors with respect to the colloid frame of reference and $P_{R_c,R_H}(r)=\mathrm{\Theta }(r-R_c)\mathrm{\Theta }(R_H-r)$ is a pulse function in the radial dimension which defines the spherical shell.

In order to have more control over the propelling force, we normalize the force field over the hydrodynamic region Γ_H, and multiply by a factor F_p. F_p is the input parameter for the magnitude of the self-propelling force. Thus the hydrodynamic force field that will be applied to the solvent reads,

${\mathit{f}}_H\left(r,\theta \right)=F_p\frac{\mathit{f}\left(r,\theta \right)}{N},\text{where}N=\left|{\int }_{{\mathrm{\Gamma }}_H}\mathit{f}\left(r,\theta \right)\text{d}V\right|.(4)$

Since in our case we are dealing with a discrete fluid (made of solvent particles), the i-th solvent particle will feel a force,

${\mathit{f}}_H^i=F_p\frac{\mathit{f}\left(r_i,{\theta }_i\right)}{N},\text{where}N=\left|\sum _{j\in {\mathrm{\Gamma }}_H}\mathit{f}\left(r_j,{\theta }_j\right)\mathrm{\Delta }V\right|(5)$

where the sum is taken over all the solvent particles that are inside Γ_H and $\mathrm{\Delta }V=r_c^3=1$. Similarly to the squirmer model, here we only consider the two first surface modes of the polar component, f_θ(r, θ), while the radial component is neglected f_r (r, θ) = 0 (see Supplementary Appendix Section 5.1 for the details). Thus, the force field becomes,

$\mathit{f}\left(r,\theta \right)=\left(B_1\sin \theta +B_2\sin \theta \cos \theta \right)P_{R_c,R_H}\left(r\right){\widehat{\mathit{e}}}_{\theta }\equiv {\mathit{f}}_{\text{col}},(6)$

for which $N=(R_H^3-R_c^3)B_1{\pi }^2/4$. In this way, the total propulsion force (which is precisely the integral appearing in Eq. 4) experimented by the colloid is just F_p. B₂ controls the force dipole contribution to the force field and thus wether we are dealing with pushers (B₂ < 0), neutrals (B₂ = 0) or pullers (B₂ > 0). Because of this formulation, B₁ plays no role and will be fixed to 1.0 from here on. As in the squirmer model, we define β = B₂/B₁ as the active stress parameter that controls the type of squirmer (see Figure 2). Under the assumption of Stokes flow (low Reynolds number), it is reasonable to think that the velocity field of the solvent particles will resemble that of the squirmer model¹ [28, 49].

FIGURE 2

FIGURE 2. Examples of 2D hydrodynamic redistribution force fields for all the studied cases: a pusher (A) with β = −5, a puller (B) with β = 5, a neutral squirmer (C) with β = 0 and a monomer of the active polymer (D). These correspond to a section of the 3D field passing through the equator of the redistribution sphere. The arrow inside the central colloid indicates the direction of the propulsion force F_p, with magnitude indicated in each panel. Here, the hydrodynamic region Γ_H would be the area contained between the solid and dashed circles.

Now we need to deal with the reaction force that is exerted on the colloid which will result in its thrust. Moreover, since an active colloid is an extended rigid object we would like to preserve the torque that may arise due to density fluctuations or interactions with other objects. The reaction thrust force, f_T, is applied on the nearest colloid particle (thruster or not) to each of the solvent particles and it is equal and opposite to the redistributed force on that solvent particle:

${\mathit{f}}_T^k\left(r,\theta \right)=-{\mathit{f}}_H^{i\left(k\right)}\left(r,\theta \right)(7)$

where i(k) represents the nearest solvent particle to the k-th colloid particle.

At each step, we implement the following algorithm:

1. Starting from a reference microswimmer, we identify the neighboring solvent particles around the swimmer’s thruster particles located between the swimmer’s radius (R_col for the colloid; r_c for the polymer) and the “hydrodynamic” radius R_H.

2. We compute the force field f_H in Eq. 4 at each of the neighbors positions, consistently with the swimmer’s orientation. The norm of the total distributed force is also computed.

3. For each neighbor:

3.1 We apply the corresponding normalized force.

3.2 We find the nearest agent particle and apply the same and opposite force.

In this way self-propulsion is achieved, while linear and angular momenta are locally conserved at each step. This procedure enables physically realistic modeling of the propulsion mechanism of a wide range of self-propelled systems, both living and artificial.

In case of the active polymer, we have considered a constant field modulated by $\cos \left(\theta /2\right)$ for each thruster monomer,

$\mathit{f}\left(r,\theta \right)=-\cos \left(\frac{\theta }{2}\right)P_{R_c,R_H}\left(r\right)\widehat{\mathit{e}}\equiv {\mathit{f}}_{\text{pol}},(8)$

where $\widehat{\mathit{e}}=\cos \theta {\widehat{\mathit{e}}}_r-\sin \theta {\widehat{\mathit{e}}}_{\theta }$ is the self-propulsion direction of the thruster particle. In this case, a reaction force that provides thrust to the agent is applied on each thruster particle. This force is equal and opposite to the total force distributed among the solvent particles in each step. Since in this case we are dealing with a flexible object that has many thruster particles, we need not to worry about finding the nearest agent particle, since this is already taken care of, as the force that each thruster particle redistributes is equal and opposite to the one that is exerted on it.

In Figure 2 we can see the hydrodynamic force fields, f, we have used in this work. The continuous and dashed circumferences represent the inner (R_c), and outer (R_H) radius respectively and define the region Γ_H where redistribution occurs. As an example, in this figure the propulsion force is computed as the surface integral of the vector field inside this region, ${\mathit{F}}_p=-{\int }_{{\mathrm{\Gamma }}_H}{\mathit{f}}_H(r,\theta )d\mathit{S}$, in this case computed in 2D. Since the force field is asymmetric there exists a net propulsion force that provides thrust to the agent. In the case of the polymer (d), each polymer bead (or monomer) acts as a small colloid with its own redistribution field, so in this case R_c = r_c would represent the beads radius, i.e. the thickness of the polymer.

To conclude, the total force experienced by an agent particle consists of the following contributions

$\mathit{F}={\mathit{F}}_C+{\mathit{F}}_D+{\mathit{F}}_R+{\mathit{F}}_T(9)$

where F_T is the total thrust force, computed as the sum of all the reaction forces on each colloid’s particle ${\mathit{F}}_T={\sum }_k{\mathit{f}}_T^k$.

It is worth noting that while in this study we have restricted ourselves to axisymmetric force fields, the code implementation is made for general force fields, allowing for example azimuthal flows, like those of the Volvox algae [64].

2.4.1 Quantifying Activity

To characterise a microwimmer in the solvent, we will define dimensionless numbers such as the Reynolds number and the Péclet number. For this, we will need to establish the viscosity η of the fluid. η can be numerically computed in a DPD fluid, as recently shown in ref. [65], or estimated via a mean field, as in Warren and Groot [33]. In our work, we follow the second approach, according to which the DPD solvent kinematic viscosity ν, defined as $\frac{\eta }{\rho }$, can be computed as

$\nu =\frac{\eta }{\rho }=\frac{D_{sol}}{2}+\frac{2\pi \gamma \rho r_c^5}{1575}(10)$

where the diffusion coefficient D_sol is

$D_{sol}=\frac{45k_{B}T}{2\pi \gamma \rho r_c^3}.(11)$

For more details, see Warren and Groot [33]. Note that in MPCD the viscosity can be computed as $\eta =16.05\sqrt{m_0{k}_{B}T}/a_0^2$, where m₀ and a₀ are the mass and the size of the cell used in MPCD algorithm. See [60, 66–68] for more details.

Once we know the viscosity, we compute the Reynolds number and the Péclet number. The Reynolds number quantifies the amount of inertial versus viscous forces acting on an object that moves in a fluid,

$\text{Re}=\frac{v_p{r}_c}{\nu }(12)$

where r_c = 1 is the solvent characteristic length and v_p is the microwimmer’s propulsion velocity. For both the active colloids and the active polymer, the velocity is the one of the center of mass.

The Péclet number is an adimensional number used to quantify the measure of the activity. It is directly proportional to the self-propulsion speed and to the reorientation time [69]. The Péclet number has been used in bulk suspensions of Active Brownian Particles [70, 71] to quantify particles’ activity. Differently from ABP, when dealing with swimmers (as in our case) the rotational dynamics of the colloid arises from interaction with the fluid and it is not prescribed (such as in ABP). Therefore we introduce a Péclet number based on the propulsion velocity, v_p. In this way one should expect that, for a given set of parameters, the reorientation time increases with the propulsion velocity and thus with the Péclet. In other words, the Péclet number is defined to describe the degree of activity in the system as the ratio between the self-propulsion of the microwimmer and a diffusion scale Pe = v_p τ_r/σ. However, different works have shown different definitions for this number. We define the Péclet number for the colloid following ref. [60] as,

${\text{Pe}}_{col}=\frac{v_p{R}_{\text{col}}}{D_{col}}(13)$

here D_col = k_BT/6πηR_col is the estimated diffusion coefficient of the colloid and R_col is the colloid radius that is fixed to 2 for all our simulations (with the exception of the flow fields shown in Figure 3 in which R_col = 3 was chosen for better visibility). We have studied the following ranges F_p ∈ [0, 50], v_p ∈ [0, 1.25], Pe ∈ [0, 156] and Re ∈ [0, 12]. As mentioned earlier, for some parameters we cannot assume that we are in Stokes flow conditions, so we should not use the relation v_p = F_p/6πηR_c for computing the colloids propulsion velocity, this is why we “measure” it as $v_p={⟨{\mathit{v}}_{\text{col}}\cdot \widehat{\mathit{e}}⟩}_t$. The values obtained are shown in Figure 1 of the Supplementary Material. Note that it was found that the propulsion velocity of the colloid, and thus the Pe and Re, are not always linear with the propulsion force F_p. Moreover, they change whether we are dealing with pusher, neutral or puller squirmers. In Figure 1 of the Supplementary Material we show the different propulsion velocities found (and their corresponding Pe and Re) for each type of squirmer. However, in all simulations presented, we remain in the range F_p ∈ [0, 50] where the separation between the v_p values for different squirmer types is not so dramatic and the behaviour does not depart too much from linearity.

FIGURE 3

FIGURE 3. Sections of the velocity fields in the lab frame (bottom row) and moving with the colloid (top row), for the pusher (A,D), neutral (B,E) and puller (C,F) squirmer. surrounded by $\sim 10^4$ fluid particles and swimming to the right. The colour of the background is the averaged density of fluid particles, the colour of the arrows shows the fields magnitude. Here the colloid radius is R_c = 3 and F_p = 100 in order to obtain a clearer flow field. These values produce Pe ≈ {818, 655, 573} and Re ≈ {25, 20, 17} for the pusher, neutral and puller squirmers respectively.

For the polymer we follow the ref. [11] and define the Péclet number as

${\text{Pe}}_{pol}=\frac{F_p{r}_c}{k_{B}T}(14)$

where r_c = 1 represents the characteristic length of the monomers. The polymer’s lengths, N_b, studied lay in a range between 40 and 100. For this particular cases we have explored values Pe = {0.01, 0.1, 1.0} that correspond with Reynolds numbers in the laminar regimen, around Re$\approx 0.3$.

2.5 Analysis Tools

In order to characterise our systems, we compute both structural and dynamical features. Concerning the active colloid, we first establish the velocity field of the solvent surrounding the swimmer to characterise the nature of each spherical squirmer (whether pusher, puller or neutral). Next, we study its dynamics by computing the mean square displacement, and estimate the effective diffusion coefficient from its long time behaviour. When dealing with the colloid in confinement we also analyze the orientation autocorrelation function (OACF) that supplies information about the reorientation time of the colloid. Concerning the active polymer, we first characterise how activity affects its structural features by computing the radius of gyration. Next, we study its dynamics by computing the mean square displacement of the center of mass and again estimate the effective diffusion coefficient from its long time behaviour.

2.5.1 Velocity Fields

For computing the solvent velocity fields around the colloid we run simulations of a fixed colloid in the center of the box pointing to the positive x-axis. Then, we perform a binning of the simulation box and average the velocities of the solvent particles inside each bin, finally we also take ensemble and time averages in the stationary state. The velocity fields shown in Figure 3 correspond to a slab that has the same height as the colloid (D_col). The arrows represent the xy-projection of the full 3D velocities.

2.5.2 MSD

Concerning dynamical features, we compute the mean square displacement

$\mathrm{M}\mathrm{S}\mathrm{D}\left(t\right)=⟨{\left[{\mathbf{r}}_{\text{cm}}\left(t\right)-{\mathbf{r}}_{\text{cm}}\left(0\right)\right]}^2⟩(15)$

Where r_cm indicates the position of the center of mass of the colloid/polymer. The average is taken over several colloids/polymers. The long time behaviour of the MSD, corresponds to the diffusion coefficient D, MSD(t) = 6Dt. It is worth noting that when confinement takes place in a cylinder with a small radius, it might be better to consider the system as one dimensional, thus MSD(t) = 2Dt. However, this is not our case since we consider that the agents have sufficient space to diffuse in the transverse directions. This leads to a more straightforward comparison between the different systems.

2.5.3 OACF

The orientation autocorrelation function is computed for the colloid in confinement to asses the impact of the confinement in the rotational diffusion (or equivalently, the reorientation time) of the colloid.

$\mathrm{O}\mathrm{A}\mathrm{C}\mathrm{F}\left(\mathrm{\Delta }t\right)=\frac{1}{N_{\mathrm{\Delta }t}}\sum _{i=0}^{N_{\mathrm{\Delta }t}}\widehat{\mathit{e}}\left(t_i\right)\cdot \widehat{\mathit{e}}\left(t_i+\mathrm{\Delta }t\right)(16)$

Here Δt = n dt where dt is our base time step. The scalar product of the orientation at a given time $\widehat{\mathit{e}}(t_i)$ with itself at a delayed time $\widehat{\mathit{e}}(t_i+\mathrm{\Delta }t)$ is averaged over the intervals of length Δt, starting at all the possible t_i’s, that fit into the total simulation time T_sim = N_totdt. So there would be N_Δt = N_tot − n + 1 intervals of the same length in the full simulation interval for a given n.

2.5.4 RoG

The radius of gyration R_g for the active polymer is computed according to the relation,

$R_g^2=\frac{1}{N}\sum _{k=1}^N{\left({\mathit{r}}_k-{\mathit{r}}_{\text{cm}}\right)}^2,(17)$

where r_cm is the position of the center of mass of the polymer, r_k is position of the k thruster particle and N is the number of bead of the polymer.

3 Results

In what follows we present the results obtained for both microwimmers, either in bulk or in cylindrical confinement. We start with the simplest object: the spherical squirmer (Section 3.1) characterising its hydrodynamic features (Section 3.1.1) and its dynamical properties (Section 3.1.2). When confined in a cylindrical channel, we also compute its orientation autocorrelation function (Section 3.1.3). Next, we study the more complex-shape active polymer (Section 3.2), characterizing its structural (Section 3.2.1) and dynamical (Section 3.2.2) properties, compare our results with the passive and Brownian counterpart.

3.1 Active Colloids

3.1.1 Flow Fields

To start with, we present our results for a spherical squirmer and study the velocity fields for the pusher, the puller, and the neutral swimmer.

Figure 3 displays the velocity flow fields for this three squirmers computed as explained in the previous section: pusher (a, d), neutral (b, e) and puller (c, f) squirmer. Comparing our results with the typical flow fields expected for squirmers (e.g. ref. [49]) the flow fields reported in Figure 3 are not so symmetrical, in the case of the puller and pusher lab frames (Figures 3D,F). The four characteristic vortices of the flow field when periodic boundary conditions are present [48] seem to be shifted to the negative x-direction, compressing the two at the front and stretching the two at the back. In the same way, in the relative frame, we can see smaller swirls than usual at the front of the pusher (Figure 3A) and somewhat elongated ones at the back of the puller (Figure 3C). In the case of the neutral swimmer, the characteristic source dipole of the lab frame (Figure 3E) is completely compressed against the swimmers surface, and some turbulent flow is appreciated at the edges of the y-dimension of the section. All these deviations from the usual flow fields are ascribed to inertial effects of the fluid stemming from the high Reynolds number present in our simulations [59]. In Supplementary Material we show the flow fields for a different set of parameters (lower Reynolds number, at Re ≈ 0.1) for which we find a more typical squirmer flow field [48, 49]. In ref. [48] the authors also comment on an analytical solution and state that it is indistinguishable from the one found in their simulations. An analytical solution for the flow field without PBCs in terms of a source dipole, a force dipole and a source quadrupole is also provided in ref. [46] following a different approach but compatible with the usual derivation by Blake followed by [48]. Ref. [59] studies in detail how the Reynolds number affects the flow fields around squirmers.

Confining the colloid inside a cylindrical channel has a drastic impact in the solvent flow fields (Figure 4), since the channel walls change the boundary conditions of the fluid.

FIGURE 4

FIGURE 4. Solvent flow fields of the colloid confined in an cylindrical channel (R_cyl = 3.5) in the lab frame (bottom row) and moving with the colloid (top row), for the pusher (A,D), neutral (B,E), and puller (C,F) squirmer. Here the colloid radius is R_c = 2 and F_p = 50. These values produce Pe ≈ {90, 76, 61} and Re ≈ {3.9, 3.3, 2.7} for the pusher, neutral and puller squirmer respectively.

For the pusher and the puller in the absolute frame (Figures 4A,D) we observe that the two vortices at the back and front respectively have disappeared, while the other two (at the front of the pusher and at the back of the puller) seem to have retracted to a closer position directly in front of the pusher and behind the puller. A similar damping of vorticity has already been reported in ref. [72] and may be attributed to the suppression of the fluid’s long-wavelength modes due to the confinement [73]. In the relative frame of reference, the swirls have also contracted further, and it is now difficult to distinguish them from just turbulent flow. In the case of the neutral swimmer (Figures 4B,E) the flow fields do not differ that much with respect to the ones encountered in bulk, with the exception that now there are no turbulent regions at the edges of the flow field.

3.1.2 Diffusion

When dealing with a colloidal squirmer in bulk, we study its dynamical features by estimating the long time diffusion coefficient normalised by the diffusion of a passive colloid in bulk via the center of mass mean square displacement, as explained in Section 2.5., for the three types of squirmers (Figures 5A–C).

FIGURE 5

FIGURE 5. MSDs for the three squirmer types studied. Left column: pusher [panels (A, D, and G)]. Center column: neutral [panels (B, E, and H)]. Right column: puller [panels (C, F, and I)]. Top row: bulk [panels (A, B, and C)]. Center row: in cylindrical confinement for different Pe’s for the smallest channel radius R_cyl = 3.5 [panels (D, E, and F)]. Bottom row: for different channel radii for the highest Péclet numbers available corresponding to the highest thrust force F_p = 50 [panels (G, H, and I)].

The top row of panels in Figure 5 represent the MSD for a bulk dilute suspension of pushers (a), neutrals (b) and pullers (c). Their long time behaviour corresponds to the diffusion coefficient reported in Figure 6A From the results presented, it is reasonable to conclude that the three types of squirmer diffuse almost the same for the ranges of Péclet numbers studied. As expected, the diffusion of the three of them increases when increasing their thrust force and thus their Péclet number.

FIGURE 6

FIGURE 6. Measured diffusion as a function of the Péclet number for the active colloid in bulk (A), in confinement for the narrowest channel R_cyl = 3.5 (B) and as a function of the cylinder radius for the highest Péclets (C) for the three squirmer types studied. We normalize by the diffusion of a passive colloid in bulk $D_b^0=0.032179$. In the two first plots a secondary horizontal axis shows the corresponding Reynols number.

In Figure 6A it is worth noting that as we increase the thrust force and thus the Péclet and Reynolds numbers, the diffusion behavior changes significantly. When we are in the range of Re ≪ 1 the diffusion increases significantly while we increase the Pe. When we approach Re ≈ 1 the increase in diffusion is dampened reaching what seems to be a saturation as Re ≫ 1.

The middle and bottom row of panels in Figure 5 represent the MSD for a confined dilute suspension of pushers (a), neutrals (b) and pullers (c). The middle panels study the dynamics of swimmers in a channel with the smallest radius, while varying the Peclet number for pushers (d), neutrals (e) and pullers (f). The bottom panels study the dynamics of swimmers at the highest Peclet in a channel with varying radius for pushers (g), neutrals (h) and pullers (i). When we confine the active colloid inside a cylindrical channel the symmetry between pushers and pullers is lost. Pushers will tend to reorient parallel to the wall, while pullers will do so perpendicularly. In this way, mobility of pushers should be increased while pullers should be more prone of getting “stuck” at the wall, which is reflected in the MSD curves in Figures 5D–F. This is a well known behavior for the interaction of squirmers with walls [45, 49], that has been used to explain the accumulation of certain microorganisms at surfaces [50]. The coupling between a microswimmer close to the wall and the solvent is affected because a portion of its hydrodynamics region (see Figure 2) lies outside the cylinder, where no solvent particles are present. Although this effect can contribute an additional torque, the volume of this excluded region is small compared with the rest of the hydrodynamic region, and it is not seen to affect the qualitative features of the hydrodynamic coupling between the microswimmer and the confining wall. Moreover, when studying the diffusion for different channel radii (Figure 6C), if this effect had a relevant contribution, we would expect to see a clear modulation of the diffusion varying radius for all types of squirmers, since the greater the radius, the smaller the excluded portion of the hydrodynamic region.

The same information can be recovered when plotting the OACF for each system (as will be shown in Figures 7E–G) curves).

FIGURE 7

FIGURE 7. Auto-correlation function for the colloids orientation vector. Top: for the smallest cylinder radius R_cyl = 3.5 and all the Péclet numbers studied for pusher (A), neutral (B) and puller (C). Bottom: for all the studied cylinder radii for a passive colloid (D) for which F_p = 0, Pe = 0.1 ≈ 0 and the highest Péclet numbers available (all corresponding to F_p = 50) for the three types of squirmers: pusher (E), puller (F), and neutral (G).

As shown in the MSD curves (Figures 5D–F), for the pusher we detect a slight increase at large times, while for the puller the curves collapse showing a significant decrease in its motility at large times for all studied Péclet numbers. This is due to the wall-facing effect described previously, which is consistent not only with the decrease of motility, but also with the apparent independence of the diffusion with the Péclet number. The shape of the MSDs curves for the neutral swimmer (Figure 5E) also follow from this argument. The neutral squirmer gains its thrust force symmetricaly between its front and back. Therefore, it propels on its front more than the pusher but less than the puller, and propels on its back more than the puller but less than the pusher. The fact that this system is between the two is confirmed by the MSDs curves. The diffusion curves (Figures 6A–C) show more clearly what we have just addressed.

In Figure 6C we report the normalized diffusion for the highest Péclet number of the three types of squirmers in confinement as a function of the channel radius. The major effect of varying the channel radius occurs for the pusher, while the puller and neutral squirmer’s diffusion seems to remain unaffected by it (in the studied range). This is coherent with the wall-facing argument previously described. The diffusion is a long time property, while for the studied radii the colloid reaches the channel wall at much shorter time scales. Therefore once the colloid has reached the wall, it might get stuck due to the wall-facing effect regardless of the channels radius.

3.1.3 Orientation Aturocorrelation Function

Finally, we compute the orientation auto-correlation function (OACF) when active colloids are confined in a cylindrical channel, as depicted in Figure 7. The OACF measures the rotational diffusion (or equivalently, the reorientation time) of a colloid, i.e., for how long the colloid retains its swimming direction before it is randomized by fluctuations.

The top row of Figure 7 represents the OACF for the system confined in the smallest cylinder, when varying the Péclet number. Whereas the bottom row represents the OACF for an active colloid propelling at the highest Péclet number and confined in cylinders with different radii. In the case of a pusher (Figure 7A) we detect a clear increase of the reorientation time with increasing Pe. This is expected for any non-chiral active particle which increases its Pe by increasing its propulsion force [69]. Moreover, due to the wall-rebound argument discussed previously, this effect could be amplified. When dealing with the neutral (Figure 7B) and the puller (Figure 7C) squirmers, the interpretation is less clear. It seems that in both cases starting from the lowest Pe the reorientation time increases until it reaches a point where the behaviours for both squirmers is different. For the neutral squirmer, as we keep increasing Pe the reorientation time decreases, reaching a minimum for the highest Pe. Whereas for the puller, at Pe = 14.5 there is a sharp decrease and then, as we keep increasing Pe, a slight recovery. Anyhow it is hard to draw solid conclusions in both cases. One reason could also be due to not enough statistics.

Figures 7D–G offers a much clearer interpretation. In these panels we show how the OACF changes as we vary the channel radius keeping in all cases the maximum Pe available, corresponding to the highest thrust force F_p = 50. As expected for a passive colloid (Figure 7D) the OACF is the same regardless of the channel radius. Moving now to the pusher (Figure 7E) we notice an increase of the reorientation time as we decrease the channel radius, consistent with the wall-rebound argument. For the puller (Figure 7G) we encounter the opposite behaviour, the reorientation time increases with increasing radius: this can be explained with the wall-facing argument plus the fact that when the puller is swimming against the wall it is in an unstable state, similar to when a pencil is left standing at its tip, so it will change its orientation. Some times this reorientation will lead him back to the center of the channel, but the narrower the channel, the sooner it will encounter again the wall and reorient again. For the neutral squirmer (Figure 7F) we are again in between pushers and pullers but since neutrals propel in their front side, as pullers, the behaviour observed is more similar to pullers than to pushers.

3.2 Active Polymer

In this section we present our results on structural and dynamical features of the active polymer in an explicit solvent. In particular, we focus on the radius of gyration R_g, and on the diffusion coefficient D, computed via the long-time behaviour of the mean square displacement of the polymer’s center of mass. We consider the active polymer first in bulk and then confined in a cylindrical channel, underlying the effect of the activity in comparison with the passive polymer behaviour in the same conditions. When in bulk, we unravel the effect of hydrodynamics comparing our results to the results obtained in ref. [11] for Active Brownian polymers (without hydrodynamics).

3.2.1 Radius of Gyration

Figures 8A,C shows the probability distribution function of the radius of gyration for the polymer in bulk (panel A, continuous lines) and confined in a channel (panel C, dashed lines), comparing the passive (thick lines) to the active (thin lines) case. We also sketch snapshots representing typical conformations observed in each case, both for the passive case (in blue) and for the active one (in magenta). Panel B represents the average radious of gyration as a function of the Péclet number for polymer length of N = 50 (blue), N = 100 (green) and N = 200 (red) in bulk (continous line) and in a channel (dashed line).

FIGURE 8

FIGURE 8. Probability distribution of the radius of gyration for passive (Pe = 0, thick line) and active (Pe = 10, thin line) polymers of N = 50 (blue), N = 100 (green), N = 200 (red) monomers in bulk (A) and in confinement (C). The snapshots illustrate typical conformations for the passive (blue) or active (magenta) polymers. (B) Average value of R_g as function of Pe for different polymer sizes. Continuous lines are for the polymer in bulk and dashed lines for the polymer under confinement.

When studying a passive polymer in bulk (thick lines in Figure 8A), our model recovers the expected increase of the radius of gyration with the polymer size. This behaviour is observed also in the presence of active forces as shown by the thin lines in Figure 8A. In order to underline the relevance of hydrodynamics it would be interesting to compare the results obtained for the active polymer with those for the Active Brownian Polymer reported in ref. [11]. However, a direct comparison is not possible, due to the different features of the chosen polymer’s model. In ref. [11] the authors used a bead-spring self-avoiding polymer, whereas in our study we have used an ideal polymer.

When studying the average of the radius of gyration as a function of the Péclet number (Figure 8B) we detect a non-monotonic behaviour.

For short polymers (N = 50), R_g remains constant at low activities (until Pe = 5): the same behaviour has been detected in ref. [11] for Active Brownian polymers, sign that hydrodynamics is not relevant for low activities.For relatively short polymers (N = 50 and 100) the radius of gyrations is almost constant when activity is low, and increases at high activity. This behaviour corresponds to what one would expect if the polymer behaved like a flexible polymer [17]. The collapse is a consequence of the time scale separation of the thermal and active contributions. The subsequent increase of the radius of gyration is due to the reduced influence of hydrodynamic interactions for larger values of Peclet number.However, when activity increases, the presence of hydrodynamics affects the polymer conformation since R_g increases. On the other side, without hydrodynamics, R_g decreases. For larger polymers (N > 50) R_g reaches a minimum value before increasing again. The same behaviour has been already reported in ref. [17] for active fully flexible Brownian self-avoiding polymer.

Larger polymers (N > 200) behave like a semi-flexible polymer [17], characterised by an initial decrease of R_g (more compact shape) for small values of the Péclet number, leading to an increase of R_g with the activity (more open shape). This non-monotonic behaviour resembles the behaviour observed for the end-to-end distance of active polymers in the presence of hydrodynamics [13, 14].

Even when an active polymer is confined in a cylindrical channel (Figure 8C), activity plays the same role on the probability distribution of the radius of gyration. The radius of gyration increases with the number of monomers N when hydrodynamics is taken into account. This is expected, as we increase the mass of a polymer. Moreover, comparing the active (thin) to the passive (thick) polymer, the increase is more stretched in the active than in the passive case. Interestingly, the confinement does not seem to affect R_g since same size active polymers (50, ≤ ,N ≤ 200) are characterised by the same radius of gyration when in bulk or in a channel. Probably, the reason for this is that we have chosen to study a channel whose diameter is relatively large, thus not differing too much from the bulk system.

3.2.2 Diffusion

In order to understand the dynamical features of an active polymer in bulk and in confinement, we compute the MSD of the active polymer’s center of mass.

As in ref. [11] for the system without hydrodynamics and like the squirmers studied in the previous section, the MSDs present at short time a ballistic regime (MSD ∝ t²), a diffusive dependence at long time (MSD ∝ t) and for intermediate times there is a crossover typical of a super-diffusive regime (MSD ∝ t^ν, with 1 < ν < 2). From the MSD long time dependence, we estimate the diffusion coefficient. Figure 9 represents the diffusion coefficient D_eff of the polymer’s center of mass as a function of the polymer size, when varying the Péclet number. While in panel A we have normalised the diffusion coefficient by the diffusion coefficient of the passive polymer (D₀), in panel B we have normalised the diffusion coefficient by the mean field diffusion of a DPD solvent particle (D_sol in Eq. 22).

FIGURE 9

FIGURE 9. Diffusion coefficient of the polymer center of mass D_eff as a function of the polymer length for three different Pe values 0 (red), 0.1 (orange) and 1 (yellow). (A) Normalised by the diffusion coefficient of the passive polymer (D₀); (B) normalised by the mean field diffusion of a DPD solvent particle Eq. 11. Square dots are results from ref. [11] triangles are results for the system in bulk and circles for system in confinement.

In Figure 9A. we show the results for the effective diffusion normalized by the diffusion coefficient of the passive case (D₀) as a function of the polymer size. We study values of activity ranging from the passive case (in red) to Pe = 1 (yellow case) and observe that activity increases the effective polymer diffusion. Meanwhile, if we compare the Brownian diffusion [11] (yellow square) for the same activity Pe = 1.0 with our results (yellow triangles) we detect the same dependence with N but approximately 10 times smaller. The effect of hydrodynamics is to slow down the polymers’ motion, as expected. Finally, if we compare the results obtained for the bulk system with the ones for the channel we observe how confinement does not seem to affect the small polymers (N < 50), but turns out to be relevant when the length of the polymer is increased. For longer polymers (N > 70) the confinement affects the polymer’s motion and the diffusion decreases when the polymer is too long.

On the other hand, in Figure 9B. the effective diffusion has been normalized by the mean-field diffusion given by Eq. (11). The idea to represent the data in this way was to be able to establish a power law dependence of the diffusion coefficient with the polymer size and compare it with the prediction expected for the diffusion by the Rouse and Zimm [74] theory of Gaussian chains. Within this theory, the chain center of mass diffusion D_Rouse ∝ N⁻¹ and D_Zimm ∝ N^−1/2. As shown in Figure 9B, when the polymer is passive (red line), the power law resembles that predicted by the Zimm model, which takes into account the hydrodynamic interactions between the beads of the polymer. While when the activity is relatively high (yellow line) the behavior is not similar to that expected in this model, since other effects appear due to the activity of the system.

In order to understand this behaviour, we compute the Schmidt number, defined as Sc = ν/D_sol, being ν = 1.25 the kinematic viscosity of the DPD fluid and D_sol the diffusion of the solvent particles. Estimating the kinematic viscosity via the mean field model of Groot and Warren’s [33], and measuring D_sol, we estimate Sc = 2.35 for the solvent in all simulations. Whenever the Schmidt number is larger than one, the momentum diffusion dominates and hydrodynamics is relevant. However, since in our case the Schmidt number is around one, we conclude that the hydrodynamic coupling is not too strong. Therefore, we observe both scaling regimes, Rouse and Zimm.

4 Discussion

In many paradigmatic examples of active matter such as biological microswimmers or synthetic active colloids, these are typically immersed in a solvent and the hydrodynamic interactions produced by the movement of the particles are relevant. Usually the introduction of these hydrodynamic interactions in active systems has been carried out through lattice models such as LB, that consider hydrodynamic effects but neglect thermal fluctuations, or MPCD that allow for the study of systems at low Reynolds number [75]. In this work we develop a new framework for the introduction of hydrodynamic interactions in mesoscopic molecular dynamics simulations. To do so, we have used the well known DPD model that has been showed to be a simple and well behaved coarse-grain model (for a specific set of parameters) for the implementation of hydrodynamics interactions in passive systems.

One of the main advantages of this new implementation is the possibility of easily taking into account thermal fluctuations for swimmers of complex shapes. Moreover, our implementation has been developed as an extension of the LAMMPSs open source package and will be sent to the LAMMPS developers (constantly maintaing the code), making our numerical approach readily available to everyone.

In active systems there are a plethora of different mechanisms that produce the propulsion of the microwimmers, such as beating of flagella or chemical reactions. In our approach, we focus on the fact that in all these cases, the agents exert a force on the solvent in which they are immersed in order to achieve thrust. Depending on the type of propulsion mechanism employed, the exerted force has its own distinct features but it always respects the conservation laws of the different physical quantities. The model and its implementation is described in details, mainly based on momentum conservation: this corresponds to the fact that the force experienced by the microwimmer in its propulsion must be compensated by the stresses induced in the solvent.

To verify the validity of our model, we study two particular cases whose phenomenology has been well characterised by other numerical methods. The first of these cases are spherical squirmers, which represent the simplest model of a microwimmer in which the hydrodynamics of the system is taken into account. The second example studied is an active polymer, which is nothing more than a first approximation to a slightly more complex structure: a chain of swimmers. In this case, our proposed method is applied in the same way for each of the monomers (swimmers) that form the polymer. As shown in the Results section, the proposed method leads to a phenomenology, such as flow fields, dynamical and structural features, consistent with the results obtained for the same systems studied with different numerical models.

Concerning the active colloid, we have been able to reproduce the solvent flow fields for the different types of swimmers (Figure 3), observing a characteristic deformation of the solvent flow fields due to the inertial effects present in the fluid at moderate Reynolds numbers Re ∼ 20. We have been able to asses the impact of these inertial effects on the dynamics and hydrodynamics of the swimmer and to conclude that pushers are the most efficient swimmers (in the sense that they develop a larger propulsion velocity for the same propulsion force, Figure 1 of the Supplementary Material), followed by neutrals and pullers when the Reynolds number is increased enough. For the ranges studied in the majority of our work, we showed that when swimming in bulk, diffusion is hardly affected by the choice of squirmer type (Figure 6A), and begins to saturate as we venture into higher Reynolds numbers. When the swimmer is confined inside a cylindrical channel, the flow fields changed dramatically in order for the fluid to adapt to the new boundary conditions (Figure 4). The confined geometry breaks this symmetry in the diffusion between the squirmer types (Figure 6B), as the behaviour of each swimmer near the channel wall is completely different: while pushers tend to rebound, aligning parallel to the wall and thus increasing their diffusion, pullers tend to get stuck, aligning perpendicularly with the wall and thus drastically decreasing their diffusion. Neutral squirmers lay in between both behaviours, but closer to pullers, as they slightly rely on the solvent ahead of them for achieving thrust. When varying the radius of the confining channel, diffusion of neutrals and pullers was hardly affected, while pushers enhanced their diffusion with decreasing channel radius (Figure 6C). Finally, we discussed the effect of the confinement in the reorientation time of the swimmers (Figure 7). We showed that pushers have slower reorientation dynamics the larger the Péclet number, but could not conclude anything too clear for neutrals and pullers. Although when the Péclet is the highest and we increase the channel radius, the reorientation behaviour of pushers and pullers is clearly opposite: pushers increase their reorientation time while pullers decrease it. Again, neutrals lay in-between both behaviours although a little closer to pullers than to pushers.

In the active polymer case, we have compared the radius of gyration R_g and diffusion to the system without hydrodynamics (Active Brownian Polymer). Concerning the radius of gyration R_g, the behavior of the polymer has been characterised as a function of both the polymer length and the Péclet number. Even though the radius of gyration monotonically increases with the polymer length (Figure 8), the dependence with the activity is not so straightforward. For short polymers R_g always increases with activity, whereas for long polymers it reaches a minimum value. This behaviour has been already detected in active fully flexible Brownian self-avoiding polymers. On the other hand, confinement always decreases R_g with respect to the system in bulk. When studying the dynamics of the polymer, we have compared our results with the analytical results for the Rouse and Zimm models, and concluded that our model (for the set of parameters used) is compatible with the prediction of the Rouse model at low Peclet (when hydrodynamics does not seem to play a relevant role), and is compatible to the Zimm model at higher Peclet number, when the monomers of the polymer chain can interact with each others due to hydrodynamics.

As in Ref. 76, having characterised the behaviour of individual swimmers in a solvent, we plan to use our numerical tool to study more dense suspensions of active colloids or active polymers. This will allow us to study their collective behaviour, their aggregation (if present) and the interplay played by hydrodynamics and activity, with the idea of comparing our numerical results to experiments on active synthetic colloids or active living swimmers (such as algae or bacteria), where hydrodynamics is relevant.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

CV acknowledges fundings from MINECO PID2019-105343GB-I00 and EUR2021-122001. I. Pagonabarraga acknowledges support from Ministerio de Ciencia, Innovación y Universidades MCIU/AEI/FEDER for financial support under grant agreement PGC2018-098373-B-100 AEI/FEDER-EU, from Generalitat de Catalunya under project 2017SGR-884, Swiss National Science Foundation Project No. 200021-175719 and the EU Horizon 2020 program through 766972-FET-OPEN NANOPHLOW.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2022.926609/full#supplementary-material

Footnotes

¹Under the same Stokes flow assumption, the colloid’s propulsion velocity v_p can be computed from the self-propulsion force F_p via the Stokes law F_p = 6πηR_cv_p. However, this assumption may not hold in some cases that are also worth studying. For these cases, we “measure” v_p as the time averaged projection of the colloids velocity v_col over its orientation axis $v_p={⟨{\mathit{v}}_{\text{col}}\cdot \widehat{\mathit{e}}⟩}_t$ as will be explained later on.

References

1. Gompper G, Winkler RG, Speck T, Solon A, Nardini C, Peruani F, et al. The 2020 Motile Active Matter Roadmap. J Phys Condens Matter (2020) 32(19):193001. doi:10.1088/1361-648x/ab6348

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Aranson IS. Active Colloids. Phys.-Usp. (2013) 56(1):79–92. doi:10.3367/ufne.0183.201301e.0087

CrossRef Full Text | Google Scholar

3. Zöttl A, Stark H. Emergent Behavior in Active Colloids. J Phys Condens Matter (2016) 28(25):253001. doi:10.1088/0953-8984/28/25/253001

CrossRef Full Text | Google Scholar

4. De Corato M, Arqué X, Patiño T, Arroyo M, Sánchez S, Pagonabarraga I. Self-propulsion of Active Colloids via Ion Release: Theory and Experiments. Phys Rev Lett (2020) 124(10):108001. doi:10.1103/physrevlett.124.108001

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Lucia A, Guzmán E. Emulsions Containing Essential Oils, Their Components or Volatile Semiochemicals as Promising Tools for Insect Pest and Pathogen Management. Adv Colloid Interf Sci (2021) 287:102330. doi:10.1016/j.cis.2020.102330

CrossRef Full Text | Google Scholar

6. Ebbens SJ. Active Colloids: Progress and Challenges towards Realising Autonomous Applications. Curr Opin Colloid Interf Sci (2016) 21:14–23. doi:10.1016/j.cocis.2015.10.003

CrossRef Full Text | Google Scholar

7. Hortelão AC, Patiño T, Perez-Jiménez A, Blanco À, Sánchez S. Enzyme-powered Nanobots Enhance Anticancer Drug Delivery. Adv Funct Mater (2018) 28(25):1705086. doi:10.1002/adfm.201705086

CrossRef Full Text | Google Scholar

8. Li J, Thamphiwatana S, Liu W, Esteban-Fernández de Ávila B, Angsantikul P, Sandraz E, et al. Enteric Micromotor Can Selectively Position and Spontaneously Propel in the Gastrointestinal Tract. ACS nano (2016) 10(10):9536–42. doi:10.1021/acsnano.6b04795

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Qian J, Kai G. Application of Micro/nanomaterials in Adsorption and Sensing of Active Ingredients in Traditional Chinese Medicine. J Pharm Biomed Anal (2020) 190:113548. doi:10.1016/j.jpba.2020.113548

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Palacci J, Sacanna S, Vatchinsky A, Chaikin PM, Pine DJ. Photoactivated Colloidal Dockers for Cargo Transportation. J Am Chem Soc (2013) 135(43):15978–81. doi:10.1021/ja406090s

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Bianco V, Locatelli E, Malgaretti P. Globulelike Conformation and Enhanced Diffusion of Active Polymers. Phys Rev Lett (2018) 121(21):217802. doi:10.1103/physrevlett.121.217802

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Isele-Holder RE, Elgeti J, Gompper G. Self-propelled Worm-like Filaments: Spontaneous Spiral Formation, Structure, and Dynamics. Soft matter (2015) 11(36):7181–90. doi:10.1039/c5sm01683e

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Eisenstecken T, Gompper G, Winkler R. Conformational Properties of Active Semiflexible Polymers. Polymers (2016) 8(8):304. doi:10.3390/polym8080304

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Winkler RG, Gompper G. The Physics of Active Polymers and Filaments. J Chem Phys (2020) 153(4):040901. doi:10.1063/5.0011466

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Michieletto D. Non-equilibrium Living Polymers. Entropy (2020) 22(10):1130. doi:10.3390/e22101130

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Locatelli E, Bianco V, Malgaretti P. Activity-induced Collapse and Arrest of Active Polymer Rings. Phys Rev Lett (2021) 126(9):097801. doi:10.1103/physrevlett.126.097801

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Das S, Kennedy N, Cacciuto A. The Coil-Globule Transition in Self-Avoiding Active Polymers. Soft Matter (2021) 17(1):160–4. doi:10.1039/d0sm01526a

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Deblais A, Woutersen S, Bonn D. Rheology of Entangled Active Polymer-like T. Tubifex Worms. Phys Rev Lett (2020) 124(18):188002. doi:10.1103/physrevlett.124.188002

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Martín-Gómez A, Eisenstecken T, Gompper G, Winkler RG. Active Brownian Filaments with Hydrodynamic Interactions: Conformations and Dynamics. Soft Matter (2019) 15(19):3957–69. doi:10.1039/c9sm00391f

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Elgeti J, Winkler RG, Gompper G. Physics of Microswimmers-Single Particle Motion and Collective Behavior: a Review. Rep Prog Phys (2015) 78(5):056601. doi:10.1088/0034-4885/78/5/056601

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Matas-Navarro R, Golestanian R, Liverpool TB, Fielding SM. Hydrodynamic Suppression of Phase Separation in Active Suspensions. Phys Rev E Stat Nonlin Soft Matter Phys (2014) 90(3):032304. doi:10.1103/PhysRevE.90.032304

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Ishikawa T, Pedley TJ. Coherent Structures in Monolayers of Swimming Particles. Phys Rev Lett (2008) 100(8):088103. doi:10.1103/PhysRevLett.100.088103

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Llopis I, Pagonabarraga I. Dynamic Regimes of Hydrodynamically Coupled Self-Propelling Particles. Europhys Lett (2006) 75(6):999–1005. doi:10.1209/epl/i2006-10201-y

CrossRef Full Text | Google Scholar

24. Schwarzendahl FJ, Mazza MG. Maximum in Density Heterogeneities of Active Swimmers. Soft Matter (2018) 14:4666–78. doi:10.1039/c7sm02301d

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Schwarzendahl FJ, Mazza MG. Hydrodynamic Interactions Dominate the Structure of Active Swimmers' Pair Distribution Functions. J Chem Phys (2019) 150(18):184902. doi:10.1063/1.5085755

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Stenhammar J, Nardini C, Nash RW, Marenduzzo D, Morozov A. Role of Correlations in the Collective Behavior of Microswimmer Suspensions. Phys Rev Lett (2017) 119:028005. doi:10.1103/PhysRevLett.119.028005

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Liu Z, Zeng W, Ma X, Cheng X. Density Fluctuations and Energy Spectra of 3d Bacterial Suspensions. Soft Matter (2021) 17:10806–17. doi:10.1039/d1sm01183a

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Lighthill MJ. On the Squirming Motion of Nearly Spherical Deformable Bodies through Liquids at Very Small reynolds Numbers. Comm Pure Appl Math (1952) 5(2):109–18. doi:10.1002/cpa.3160050201

CrossRef Full Text | Google Scholar

29. Blake JR. A Spherical Envelope Approach to Ciliary Propulsion. J Fluid Mech (1971) 46(1):199–208. doi:10.1017/s002211207100048x

CrossRef Full Text | Google Scholar

30. Ramachandran S, Sunil Kumar PB, Pagonabarraga I. A Lattice-Boltzmann Model for Suspensions of Self-Propelling Colloidal Particles. Eur Phys J E (2006) 20(2):151–8. doi:10.1140/epje/i2006-10009-1

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Lauga E, Powers TR. The Hydrodynamics of Swimming Microorganisms. Rep Prog Phys (2009) 72(9):096601. doi:10.1088/0034-4885/72/9/096601

CrossRef Full Text | Google Scholar

32. Theers M, Westphal E, Gompper G, Winkler RG. Modeling a Spheroidal Microswimmer and Cooperative Swimming in a Narrow Slit. Soft Matter (2016) 12(35):7372–85. doi:10.1039/c6sm01424k

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Groot RD, Warren PB. Dissipative Particle Dynamics: Bridging the gap between Atomistic and Mesoscopic Simulation. J Chem Phys (1997) 107(11):4423–35. doi:10.1063/1.474784

CrossRef Full Text | Google Scholar

34. Succi S. The Lattice Boltzmann Equation: For Fluid Dynamics and beyond. Oxford University Press (2001).

Google Scholar

35. Malevanets A, Kapral R. Mesoscopic Model for Solvent Dynamics. J Chem Phys (1999) 110(17):8605–13. doi:10.1063/1.478857

CrossRef Full Text | Google Scholar

36. Kapral R. Multiparticle Collision Dynamics: Simulation of Complex Systems on Mesoscales. Adv Chem Phys (2008) 140:89–146. doi:10.1002/9780470371572.ch2

CrossRef Full Text | Google Scholar

37. Kroll DM, Winkler RG, Gompper G, Ihle T. Multi-particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids. In: Advanced Computer Simulation Approaches for Soft Matter Sciences III (2009). p. 1–87.

Google Scholar

38. Hoogerbrugge PJ, Koelman JMVA. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhys Lett (1992) 19(3):155–60. doi:10.1209/0295-5075/19/3/001

CrossRef Full Text | Google Scholar

39. Cates ME, Stratford K, Adhikari R, Stansell P, Desplat J-C, Pagonabarraga I, et al. Simulating Colloid Hydrodynamics with Lattice Boltzmann Methods. J Phys Condens Matter (2004) 16(38):S3903–S3915. doi:10.1088/0953-8984/16/38/009

CrossRef Full Text | Google Scholar

40. Stratford K, Pagonabarraga I. Parallel Simulation of Particle Suspensions with the Lattice Boltzmann Method. Comput Mathematics Appl (2008) 55(7):1585–93. doi:10.1016/j.camwa.2007.08.018

CrossRef Full Text | Google Scholar

41. Español P, Warren P. Statistical Mechanics of Dissipative Particle Dynamics. Europhys Lett (1995) 30(4):191–6. doi:10.1209/0295-5075/30/4/001

CrossRef Full Text | Google Scholar

42. Wei J, Liu Y, Song F. Coarse-grained Simulation of the Translational and Rotational Diffusion of Globular Proteins by Dissipative Particle Dynamics. J Chem Phys (2020) 153(23):234902. doi:10.1063/5.0025620

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Fedosov DA, Pan W, Caswell B, Gompper G, Karniadakis GE. Predicting Human Blood Viscosity In Silico. Proc Natl Acad Sci U.S.A (2011) 108(29):11772–7. doi:10.1073/pnas.1101210108

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Blake JR. A Note on the Image System for a Stokeslet in a No-Slip Boundary. Math Proc Camb Phil Soc (1971) 70(2):303–10. doi:10.1017/s0305004100049902

CrossRef Full Text | Google Scholar

45. Spagnolie SE, Lauga E. Hydrodynamics of Self-Propulsion Near a Boundary: Predictions and Accuracy of Far-Field Approximations. J Fluid Mech (2012) 700:105–47. doi:10.1017/jfm.2012.101

CrossRef Full Text | Google Scholar

46. Pak OS, Lauga E. Chapter 4 Theoretical Models of low-reynolds-number Locomotion. In: Fluid-Structure Interactions in Low-Reynolds-Number Flows. London: Soft Matter Series (2016). p. 100–67.

Google Scholar

47. Zhu L, Lauga E, Brandt L. Low-reynolds-number Swimming in a Capillary Tube. J Fluid Mech (2013) 726:285–311. doi:10.1017/jfm.2013.225

CrossRef Full Text | Google Scholar

48. Kuron M, Stärk P, Burkard C, De Graaf J, Holm C. A Lattice Boltzmann Model for Squirmers. J Chem Phys (2019) 150(14):144110. doi:10.1063/1.5085765

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Zöttl A, Stark H. Simulating Squirmers with Multiparticle Collision Dynamics. Eur Phys J E Soft Matter (2018) 41(5):61–11. doi:10.1140/epje/i2018-11670-3

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Berke AP, Turner L, Berg HC, Lauga E. Hydrodynamic Attraction of Swimming Microorganisms by Surfaces. Phys Rev Lett (2008) 101:038102. doi:10.1103/PhysRevLett.101.038102

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Kurzthaler C, Mandal S, Bhattacharjee T, Löwen H, Datta SS, Stone HA. A Geometric Criterion for the Optimal Spreading of Active Polymers in Porous media. Nat Commun (2021) 12(1):7088. doi:10.1038/s41467-021-26942-0

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Mandal S, Lang S, Gross M, Oettel M, Raabe D, Franosch T, et al. Multiple Reentrant Glass Transitions in Confined Hard-Sphere Glasses. Nat Commun (2014) 5(1):4435. doi:10.1038/ncomms5435

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Lobaskin V, Dünweg B. A New Model for Simulating Colloidal Dynamics. New J Phys (2004) 6:54. doi:10.1088/1367-2630/6/1/054

CrossRef Full Text | Google Scholar

54. de Graaf J, Peter T, Fischer LP, Holm C. The Raspberry Model for Hydrodynamic Interactions Revisited. Ii. The Effect of Confinement. J Chem Phys (2015) 143(8):084108. doi:10.1063/1.4928503

PubMed Abstract | CrossRef Full Text | Google Scholar

55. de Graaf J, Menke H, Mathijssen AJTM, Fabritius M, Holm C, Shendruk TN. Lattice-boltzmann Hydrodynamics of Anisotropic Active Matter. J Chem Phys (2016) 144(13):134106. doi:10.1063/1.4944962

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Plimpton S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J Comput Phys (1995) 117(1):1–19. doi:10.1006/jcph.1995.1039

CrossRef Full Text | Google Scholar

57. Nash RW, Adhikari R, Tailleur J, Cates ME. Run-and-tumble Particles with Hydrodynamics: Sedimentation, Trapping, and Upstream Swimming. Phys Rev Lett (2010) 104:258101. doi:10.1103/PhysRevLett.104.258101

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Menzel AM, Saha A, Hoell C, Löwen H. Dynamical Density Functional Theory for Microswimmers. J Chem Phys (2016) 144(2):024115. doi:10.1063/1.4939630

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Chisholm NG, Legendre D, Lauga E, Khair AS. A Squirmer across reynolds Numbers. J Fluid Mech (2016) 796:233–56. doi:10.1017/jfm.2016.239

CrossRef Full Text | Google Scholar

60. Kuhr J-T, Blaschke J, Rühle F, Stark H. Collective Sedimentation of Squirmers under Gravity. Soft Matter (2017) 13(41):7548–55. doi:10.1039/c7sm01180f

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Kubo R. Statistical-mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems. J Phys Soc Jpn (1957) 12(6):570–86. doi:10.1143/jpsj.12.570

CrossRef Full Text | Google Scholar

62. Pivkin IV, Karniadakis GE. Controlling Density Fluctuations in wall-bounded Dissipative Particle Dynamics Systems. Phys Rev Lett (2006) 96:206001. doi:10.1103/PhysRevLett.96.206001

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Löwen H. Inertial Effects of Self-Propelled Particles: From Active Brownian to Active Langevin Motion. J Chem Phys (2020) 152(4):040901. doi:10.1063/1.5134455

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Drescher K, Leptos KC, Tuval I, Ishikawa T, Pedley TJ, Goldstein RE. DancingVolvox: Hydrodynamic Bound States of Swimming Algae. Phys Rev Lett (2009) 102(16):168101. doi:10.1103/physrevlett.102.168101

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Panoukidou M, Wand CR, Carbone P. Comparison of Equilibrium Techniques for the Viscosity Calculation from Dpd Simulations. Soft Matter (2021) 17(36):8343–53. doi:10.1039/d1sm00891a

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Noguchi H, Gompper G. Dynamics of Fluid Vesicles in Shear Flow: Effect of Membrane Viscosity and thermal Fluctuations. Phys Rev E Stat Nonlin Soft Matter Phys (2005) 72(1):011901. doi:10.1103/PhysRevE.72.011901

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Kikuchi N, Pooley CM, Ryder JF, Yeomans JM. Transport Coefficients of a Mesoscopic Fluid Dynamics Model. J Chem Phys (2003) 119(12):6388–95. doi:10.1063/1.1603721

CrossRef Full Text | Google Scholar

68. Tüzel E, Strauss M, Ihle T, Kroll DM. Transport Coefficients for Stochastic Rotation Dynamics in Three Dimensions. Phys Rev E Stat Nonlin Soft Matter Phys (2003) 68(3):036701. doi:10.1103/PhysRevE.68.036701

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Martin-Roca J, Martinez R, Alexander LC, Diez AL, Aarts DGAL, Alarcon F, et al. Characterization of Mips in a Suspension of Repulsive Active Brownian Particles through Dynamical Features. J Chem Phys (2021) 154(16):164901. doi:10.1063/5.0040141

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Stenhammar J, Marenduzzo D, Allen RJ, Cates ME. Phase Behaviour of Active Brownian Particles: the Role of Dimensionality. Soft Matter (2014) 10:1489–99. doi:10.1039/c3sm52813h

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Fang L, Li LL, Guo JS, Liu YW, Huang XR. Time Scale of Directional Change of Active Brownian Particles. Phys Lett A (2022) 427:127934. doi:10.1016/j.physleta.2022.127934

CrossRef Full Text | Google Scholar

72. Pagonabarraga I, Hagen MHJ, Lowe CP, Frenkel D. Short-time Dynamics of Colloidal Suspensions in Confined Geometries. Phys Rev E (1999) 59:4458–69. doi:10.1103/physreve.59.4458

CrossRef Full Text | Google Scholar

73. Bocquet L, Barrat J-L. Hydrodynamic Properties of Confined Fluids. J Phys Condens Matter (1996) 8:9297–300. doi:10.1088/0953-8984/8/47/019

CrossRef Full Text | Google Scholar

74. Zimm BH. Dynamics of Polymer Molecules in Dilute Solution: Viscoelasticity, Flow Birefringence and Dielectric Loss. J Chem Phys (1956) 24(2):269–78. doi:10.1063/1.1742462

CrossRef Full Text | Google Scholar

75. Michael PH, Arash N, Jeremy CP. Modeling Hydrodynamic Interactions in Soft Materials With Multiparticle Collision Dynamics. Curr Opin Chem Eng (2019) 23:34–43. doi:10.1016/j.coche.2019.02.007

CrossRef Full Text | Google Scholar

76. Qi K, Westphal E, Gompper G, Winkler RG. Emergence of Active Turbulence in Microswimmer Suspensions due to Active Hydrodynamic Stress and Volume Exclusion. Communications Physics (2022) 5:1–12.

CrossRef Full Text | Google Scholar