Hydrodynamic interactions dominate the structure of active swimmers' pair distribution functions

TL;DR

This study shows that hydrodynamic interactions significantly influence the pair distribution functions of microswimmer suspensions, reducing overall structure and creating orientation-dependent patterns, advancing understanding of active matter systems.

Contribution

It provides the first detailed analysis of how hydrodynamics affect the statistical structure of microswimmer liquids through particle-based simulations.

Findings

Hydrodynamic interactions alter the orientation-dependent pair distribution function.

Hydrodynamics decrease the overall liquid structure compared to passive models.

At higher densities, a stronger first neighbor shell forms in specific directions.

Abstract

Microswimmers often exhibit surprising patterns due to the nonequilibrium nature of their dynamics. Collectively, suspensions of microswimmers appear as a liquid whose properties set it apart from its passive counterpart. To understand the impact of hydrodynamic interactions on the basic statistical features of a microswimmer's liquid, we investigate its structure by means of the pair distribution function. We perform particle-based simulations of microswimmers that include steric effects, shape anisotropy, and hydrodynamic interactions. We find that hydrodynamic interactions considerably alter the orientation-dependent pair distribution function compared to purely excluded-volume models like active Brownian particles, and generally decrease the structure of the liquid. Depletion regions are dominant at lower filling fractions, while at larger filling fraction the microswimmer liquid develops a stronger first shell of neighbors in specific directions, while losing structure at larger distances. Our work is a first step towards a statistico-mechanical treatment of the structure of microswimmer suspensions.

Figures

Equations

$$\frac{\text{d}\mathbf{r}}{\text{d}t}=v_0\mathbf{e}+\mathbf{F}/\gamma +\mathbf{\eta },$$

$$\frac{\text{d}\mathbf{e}}{\text{d}t}=\mathbf{\zeta }\times \mathbf{e},$$

$$g(r_b,{\theta }_b)=\frac{h(r_b,{\theta }_b)}{N^{2}n(r_b,{\theta }_b)},$$

$$n(r_b,{\theta }_b)=$$

$$\times [\cos ({\theta }_b)-\cos ({\theta }_b+\delta \theta )],$$

$$g(r)={\int }_0^{\pi }g(r,\theta )\mathrm{d}\theta .$$

$$C={\int }_0^{R_{\min }}g(r)\mathrm{d}r,$$

Full text

I Introduction

Collections of microorganisms such as bacteria or microalgae show a myriad of self-organized structures and patterns that are also governed by physical processes stemming from their nonequilibrium nature and mutual interactions. Some examples are self-concentration Schwarzendahl and Mazza (2018); Dombrowski et al. (2004), complex interaction with solid surfaces Ostapenko et al. (2018); Lauga et al. (2006); Berke et al. (2008), bacterial turbulence Wensink et al. (2012), swarming Copeland and Weibel (2009), spontaneous formation of spiral vortices Wioland et al. (2013), and directed motion Wioland, Lushi, and Goldstein (2016). Many of these systems collectively resemble the behavior and structure of a fluid and therefore it has been found useful to model these collections of organisms as an active or living fluid Koch and Subramanian (2011); Marchetti et al. (2013). Furthermore, the methods of statistical mechanics of classical fluids have been found to be powerful tools also in nonequilibrium situations, such as the analysis of velocity correlations and energy spectra for mesoscale turbulence of bacteria Wensink et al. (2012), and even a local equilibrium Maxwell–Boltzmann approximation for sedimenting active colloids Enculescu and Stark (2011).

The structure of an equilibrium fluid is typically analyzed in terms of the pair distribution function Allen and Tildesley (1989); Gray and Gubbins (1984), and is a cornerstone of liquid state theory Hansen and McDonald (1990). For active fluids, some work on the structural properties in terms of pair distribution functions has been put forward Bialké, Löwen, and Speck (2013); Härtel, Richard, and Speck (2018); Hoell, Löwen, and Menzel (2018); Wittkowski, Stenhammar, and Cates (2017); Riedel, Kruse, and Howard (2005); Alarcón and Pagonabarraga (2013); Kogler and Klapp (2015); Yang et al. (2017); Pessot, Löwen, and Menzel (2018); Locatelli et al. (2015); Theers et al. (2018); de Macedo Biniossek et al. (2018); Pessot, Löwen, and Menzel (2018). An important finding is the anisotropic form of the pair correlation function of active Brownian particles (ABP) presented in Bialké, Löwen, and Speck (2013); Härtel, Richard, and Speck (2018); Hoell, Löwen, and Menzel (2018); Wittkowski, Stenhammar, and Cates (2017), which arises from their self-propulsion.

Typically, pair correlation functions are strongly influenced by the mutual interactions between particles, which for microswimmers such as bacteria or algae moving in an aqueous milieu include hydrodynamic interactions. In this work we analyze the influence of the hydrodynamic interactions on the orientation-dependent pair distribution function of active swimmers, by considering a recent model for biological swimmers which also includes shape anisotropy and steric effects Schwarzendahl and Mazza (2018). In order to disentangle the effects of activity, hydrodynamic and steric interactions, we compare the pair distribution functions for our model to the case of active Brownian particles. We find that hydrodynamic interactions introduce considerably more structure to the pair distribution function at low filling fractions and suppress structures at high filling fractions.

The rest of this work is organized as follows. Section II introduces our model and the setup of the numerical simulations. In Sec. III we discuss the orientation-dependent pair distribution function, and its angular average: the radial distribution function. Section IV presents the distribution functions for both our swimmer’s model and for ABP. Finally, in Sec. V we collect our conclusions.

II Models and parameters

II.1 Stroke-averaged biological microswimmer

We employ the stroke-averaged biological microswimmer model presented in Schwarzendahl and Mazza (2018), which consists of an asymmetric dumbbell. Figure 1(a) shows a schematic representation of the puller-type model swimmer mimicking a Chlamydomonas reinhardtii cell, while Fig. 1(b) shows a schematic of the pusher-type model swimmer representing an Escherichia coli cell. In both cases, the small green sphere, marked with (1), represents the swimmer body and the large sphere, marked with (2), represents the stroke-average region spanned by the flagellar motion. The rigid motion of the dumbbells is governed by Newton’s equations and implemented via quaternion dynamics (for details see Schwarzendahl and Mazza (2018)).

As hydrodynamic interactions between microswimmers can play an important role in their collective behavior, we also include an explicit model for the surrounding fluid by means of the multiparticle collision dynamics (MPCD) technique. MPCD is a particle-based method that reproduces hydrodynamic modes up to the Navier–Stokes level of description Gompper et al. (2009).

Experimental measurements Drescher et al. (2011) have shown that the flow field of the pusher-type swimmer can be modeled by a force dipole, whereas puller-type swimmers’ flow fields are well described by a combination of three-Stokeslets Drescher et al. (2010). We use regularized force regions, as shown by the red spheres with embedded arrows (3) in Fig. 1, to account for the hydrodynamic flow fields. Additionally, the dumbbells are coupled to the fluid by using a no-slip boundary condition on their surface (for more details see Schwarzendahl and Mazza (2018)).

The temperature $T$ of the fluid, the size of an MPCD grid cell $a$, and the mass of an MPCD particle $m$, are the reference values for temperature, length, and mass, respectively. We perform simulations in a cubic domain of size $100a$ with periodic boundary conditions in all three directions; each MPCD cell contains an average of $\langle N_{\mathsf{C}}\rangle=20$ MPCD particles, giving a total of $2\times 10^{7}$ MPCD particles in the system. The effective radius of an individual swimmer is $\sigma\approx 5a$. For further details on the numerical implementation the reader is referred to Schwarzendahl and Mazza (2018).

II.2 Active Brownian particles

In addition to the our stroke-averaged biological microswimmer, we also carry out simulations of ABP, to disentangle the effects of activity, hydrodynamic and steric interactions. ABP are spherical particles with radius $\sigma$ that self-propel along their orientation $\bm{e}$ with a typical speed $v_{0}$ (see also Wysocki, Winkler, and Gompper (2014); Fily and Marchetti (2012); Redner, Hagan, and Baskaran (2013); Bialké, Speck, and Löwen (2012)). The translational motion of the particle’s position $\bm{r}$ is governed by the equation

[TABLE]

where the force $\bm{F}$ between the particles is purely steric and is calculated from a Weeks–Chandler–Anderson potential Weeks, Chandler, and Andersen (1971). Here, $\bm{\eta}$ is a Gaussian white noise with zero mean and $\langle\bm{\eta}(t)\bm{\eta}(t^{\prime})\rangle=2D\bm{1}\delta(t-t^{\prime})$, where $D$ is the translational diffusion coefficient. The motion of the orientational degrees of freedom is determined by the equation

[TABLE]

where rotational diffusion is included by the Gaussian white noise $\bm{\zeta}$ with $\langle\bm{\zeta}(t)\bm{\zeta}(t^{\prime})\rangle=2D_{r}\bm{1}\delta(t-t^{\prime})$ and the rotational diffusion coefficient is given by $D_{r}=3D/(2\sigma)^{2}$. Also for the ABP model, we use a three dimensional cubic domain with periodic boundary conditions.

II.3 Dimensionless parameters

For both our stroke-averaged biological microswimmer and ABP we simulate $N=300-1350$ active swimmers, resulting in filling fractions from $\phi=0.05$ to $\phi=0.25$, at a fixed volume $(100a)^{3}$. A dimensionless quantification of the relative importance of the self-propulsion speed to diffusive processes is the Péclet number, which we define as $\mathcal{P}={v_{\mathrm{eff}}\sigma}/{D}$, where $v_{\mathrm{eff}}$ and $\sigma$ are the effective speed and the effective size of the swimmer, respectively, and $D$ the diffusion constant of the swimmer (for ABP $v_{\mathrm{eff}}=v_{0}$.) We present results at $\mathcal{P}\approx 1.9\times 10^{3}$ and $\mathcal{P}\approx 6\times 10^{2}$. The Reynolds number of the flow around our swimmers (measuring the ratio of inertial to viscous forces) is given by $\mathcal{R}=\sigma v_{\mathrm{eff}}\rho/\eta$, where $\rho$ and $\eta$ are the MPCD fluid’s density and viscosity, respectively. The simulated Reynolds numbers are $\mathcal{R}=0.06$ for high Péclet number ($\mathcal{P}\approx 1.9\times 10^{3}$) and $\mathcal{R}=0.02$ at lower Péclet number ($\mathcal{P}\approx 6\times 10^{2}$).

III Pair distribution function

In general, the complete information of a liquid’s structure is encapsulated in the angular pair distribution function $g(\bm{r}_{\mu\nu},\bm{e}_{\mu},\bm{e}_{\nu})$ of two particles $\mu$ and $\nu$, which depends on their relative position $\bm{r}_{\mu\nu}=\bm{r}_{\mu}-\bm{r}_{\nu}$ as well as on their orientations $\bm{e}_{\mu}$, $\bm{e}_{\nu}$ Gray and Gubbins (1984); Allen and Tildesley (1989). However, in the general case the arguments of the pair distribution function are high dimensional and it has thus been proven useful to restrict oneself to a partial descriptions of the orientational ordering Allen and Tildesley (1989); Gray and Gubbins (1984). Since our active swimmers self-propel into the direction of their orientation $\bm{e}_{\mu}$, we are expecting non-trivial correlations mainly with respect to this axis. Therefore, we compute correlation of a particles orientation with respect to the relative position of all other particles, as illustrated in Fig. 2.

In practice, we use the following definition of our orientation-dependent pair distribution function: given two particles $\mu$ and $\nu$ with relative position $\bm{r}_{\mu\nu}=\bm{r}_{\mu}-\bm{r}_{\nu}$ and orientations $\bm{e}_{\mu}$, $\bm{e}_{\nu}$, we compute the distance $r=|\bm{r}_{\mu\nu}|$ and the polar angle $\theta=\bm{e}_{\mu}\cdot\bm{r}_{\mu\nu}/|\bm{r}_{\mu\nu}|$ with respect to the orientation of particle $\mu$ (Fig. 2). The average over the ensemble, polar angle, and orientation of the neighboring swimmer $\bm{e}_{\nu}$ is obtained from a histogram $h(r_{b},\theta_{b})$ for all $N^{2}$ pairs of swimmers, with radial bins $r_{b}$ and polar bins $\theta_{b}$. The orientation-dependent pair distribution function is then given by

[TABLE]

where $n(r_{b},\theta_{b})$ is the normalization

[TABLE]

where $\delta r$ is the bin size in the radial component and $\delta\theta$ is the bin size in the polar component. It should be noted that $g(r,\theta)$ for a homogeneous configuration is equal to one. Physically, $g(r,\theta)$ is related to the probability of finding a neighboring particle at distance $r$ and with a polar angle $\theta$ with respect to the reference particle and its orientation.

In addition to the orientation-dependent pair distribution function we calculate the radial distribution function Allen and Tildesley (1989); Hansen and McDonald (1990) given by

[TABLE]

Furthermore, as an effective measure of the local structure around a microswimmer, we consider the coordination number Hansen and McDonald (1990), which is the average number of nearest neighbours of a particle

[TABLE]

where $R_{\mathrm{min}}$ is the position of the first minimum of $g(r)$.

IV Results

Figure 3 shows a matrix chart of our calculations of [Eq. (3)]. Different rows represent different Péclet numbers for both puller- and pusher-type swimmers, whereas the filling fraction changes with the column. First of all, low filling fractions are characterized by a complex texture of peaks and minima of $g(r,\theta)$. These are stationary state structures, and not transient, as explicitly verified in our simulations.

Pushers have increased probability to find neighbors in the frontal and rear polar regions of their bodies, while around their equator neighboring swimmers tend to be depleted. As the Péclet number increases, the asymmetry between front and back poles increases, and the rear tail region grows more populated than the front. Pullers present the most complex texture at low filling fractions. A preference for an accumulation of probability in the front pole and in the rear region at intermediate latitudes is, however, discernible. Furthermore, at low filling fraction, $g(r,\theta)$ for both pullers and pushers shows a rather long-ranged structure, extending more than 5 times the radius of the individual swimmer.

As the filling fraction increases, $g(r,\theta)$ loses structure for both pushers and pullers at all Péclet numbers. For pushers, perpendicular to the swimming direction the pair distribution function first increases to $g(r,\theta)>1$ and then decreases to values $g(r,\theta)<1$. Furthermore, on the sides of the pushers one observes a band like structure for $\phi\gtrsim 0.10$ and a quadrupolar structure for $\phi=0.06$. The increased values of $g(r,\theta)$ in front and back poles of pushers become more concentrated in smaller areas as the filling fraction is increased.

For pullers, we observe that as $\phi$ increases, $g(r,\theta)>1$ exhibits a ring-like structure of increased probability close to the swimmers. Furthermore, we find two distinct regions at intermediate latitudes in the rear region of the swimmer, within the ring-like structure, that show very large values of $g(r,\theta)$.

To disentangle the effects of activity, steric and hydrodynamic interactions on the pair distribution function, we performed additional simulations of ABP Wysocki, Winkler, and Gompper (2014); Fily and Marchetti (2012); Redner, Hagan, and Baskaran (2013); Bialké, Speck, and Löwen (2012), which only include steric interactions among swimmers. Figure 4 shows the orientation-dependent pair distribution function of ABP at different Péclet numbers and filling fractions. Even from a cursory inspection it is visible that the complex texture present at lower $\phi$ has nearly vanished. In its stead, a partial ring (interrupted in the lower pole) of pronounced peak of $g(r,\theta)$ emerges. A clear increase to $g(r,\theta)>1$ in front and a clear decrease to to $g(r,\theta)<1$ behind the active Brownian particle can be seen. This asymmetry is an effect of the self-propulsion of the active Brownian particle and has also been studied in Bialké, Löwen, and Speck (2013); Härtel, Richard, and Speck (2018); Hoell, Löwen, and Menzel (2018); Wittkowski, Stenhammar, and Cates (2017).

The comparison of active Brownian particles to our active swimmer model shows that the hydrodynamic interactions have a strong impact on the pair distribution function. The long-ranged structures that are observed in Fig. 3 at low filling fractions are not present for active Brownian particles in Fig. 4, and thus we conclude that these are a consequence of the hydrodynamic interactions between swimmers.

At high filling fractions, however, the opposite effect is observed: active Brownian particles show more structure than our biological swimmer model, with a second and third order of rings of $g(r,\theta)>1$ in front and extending to the sides of the ABP. Again, this a consequence of hydrodynamic interactions.

The regions of increased probability [$g(r,\theta)>1$] around the equator of pushers, and in the rear region at intermediate latitudes for pullers strongly correlate with the structure of the respective flow fields. For pushers, the streamlines will bring a neighboring swimmer to the side of the one under consideration [see Fig. 5(a)]. For pullers, the stagnation point in the rear region at intermediate latitudes correlates well with the region of enhanced $g(r,\theta)$ [see Fig. 5(b)]. The strong influence from the specific flow field has also been recognized in Pessot, Löwen, and Menzel (2018) in the context of a minimal microswimmer model.

Although the orientation-dependent pair distribution function already contains a coarse-grained level of description, complex structures are visible. It proves useful to observe an even more coarse-grained quantity $g(r)$, obtained from averaging over the polar angle $\theta$ and that will capture the average environment experienced by a swimmer over a long period of time. Figure 6 shows the radial distribution function $g(r)$ for pushers, pullers, and ABP. ABP exhibit a remarkably pronounced first peak and a much weaker second peak; a similar structure for ABP was also found by de Macedo Biniossek et al. (2018); Farage, Krinninger, and Brader (2015).

The first shell of neighbors for both pushers and pullers is considerably weaker than the ABP’s. This is a combination of the hydrodynamic interactions that quickly reorient particles, and the typically long reorientation times of ABP combined with steric effects Matas-Navarro et al. (2014). Even at this considerable level of coarse-graining, it is visible how the different flow fields of pushers and pullers generate different structures in $g(r)$. Pushers show only a single peak, whereas pullers have a weak first and a stronger second peak. Comparing pushers and pullers with ABP reveals that hydrodynamic interactions drastically reduce the height of the first peak but enhance its width.

To study the nearest-neighbor environment in more detail, we report the dependence of the coordination number $C$ on the filling fraction $\phi$ in Fig. 7. As $\phi$ increases, $C$ increases monotonically for pushers, pullers, and ABP at both Péclet numbers investigated. However, for ABPs the coordination number is always larger than for pushers and pullers, confirming the result that the hydrodynamic interactions decrease the number of nearest neighbors.

V Conclusions

We have presented the orientation-dependent pair distribution function $g(r,\theta)$ for hydrodynamically interacting active swimmers at different Péclet numbers and a range of filling fractions. To disentangle the effects for activity, steric and hydrodynamic interactions we compared our results to simulations of active Brownian particles. We found that the hydrodynamic interactions introduce long-ranged structures into the orientation-dependent pair correlation function at low filling fractions. As the filling fraction increases, hydrodynamic interactions suppress the complex texture of $g(r,\theta)$. In fact, the influence of hydrodynamic interactions on the orientation-dependent pair correlation function is strong and the specific flow field introduces a non-trivial asymmetry in the $g(r,\theta)$ for both pushers and pullers.

Investigation of the radial distribution function for pushers, pullers, and ABP shows that hydrodynamic interactions strongly both suppress the height of the first peak and broaden it. As an effective measure of the average immediate neighborhood of a given microswimmer, the coordination number reveals that hydrodynamic interactions among microswimmers suppress the strong first shell of neighbors.

Our present work lays the foundations for a theoretical investigation of the complex interplay of hydrodynamics, steric forces, and active motion. An open —and ambitious— question is whether knowledge of the active liquid’s structure is sufficient for the development of a nonequilibrium statistico-mechanical theory.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

We gratefully acknowledge support from the Deutsche Forschungsgemeinschaft (SFB 937, project A20).