Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

Alignment-Induced Self-Organization of Autonomously Steering Microswimmers: Turbulence, Vortices, and Jets

Segun Goh Theoretical Physics of Living Matter, Institute for Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany    Elmar Westphal Peter Grünberg Institute and Jülich Centre for Neutron Science, Forschungszentrum Jülich, 52425 Jülich, Germany    Roland G. Winkler Theoretical Physics of Living Matter, Institute for Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany    Gerhard Gompper g.gompper@fz-juelich.de Theoretical Physics of Living Matter, Institute for Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany
(June 25, 2024)
Abstract

Systems of motile microorganisms exhibit a multitude of collective phenomena, including motility-induced phase separation and turbulence. Sensing of the environment and adaptation of movement plays an essential role in the emergent behavior. We study the collective motion of wet self-steering polar microswimmers, which align their propulsion direction hydrodynamically with that of their neighbors, by mesoscale hydrodynamics simulations. The simulations of the employed squirmer model reveal a distinct dependence on the swimmer flow field, i.e., pullers versus pushers. The collective motion of pushers is characterized by active turbulence, with nearly homogeneous density and a Gaussian velocity distribution. Pullers exhibit a strong tendency for clustering and display velocity and vorticity distributions with fat exponential tails; their dynamics is chaotic, with a temporal appearance of vortex rings and fluid jets. Our results show that the collective behavior of intelligent microswimmers is very diverse and still offers many surprises to be discovered.

I Introduction

Emergence of dynamic structures and patterns is an essential feature of biological active motile systems. Examples include microbial swarms [1, 2] on the cellular level, to schools of fish [3, 4], flocks of birds [5, 6], and collective motion in human crowds [7] on a macroscopic level. Also, in artificial active systems consisting of synthetic self-propelled particles and microrobots, the collective dynamics of constituent objects is of fundamental importance for their application in engineering and medicine to achieve a large spectrum of functionalities [8, 9, 10, 11].

A fundamental aspect in such systems is the active and autonomous motion of the constituting particles. While activity and self-propulsion can give rise to several novel types of collective behaviors, such as motility-induced phase separation (MIPS) [12, 13] and active turbulence [14, 15, 16], the fact that biological microswimmers are not only motile, but also gather information about their environment and adapt their motion by self-steering, remains largely unexplored and yet to be elucidated [17, 18, 19].

Many living organisms are immersed in a fluid medium, and their collective behavior is strongly affected or even dominated by hydrodynamics [20, 21]. The hydrodynamic environment is not just the background medium in which aquatic microorganisms are based, but it is rather essential for locomotion on the individual level as well as inter-organism interactions [1]. On the mesoscale, life has adapted to low-Reynolds hydrodynamics, e.g., in the emergence of bacterial turbulence [22, 23, 16] and in coordinated cell migration during embryogenesis [24, 25]; similar physical laws govern the swarming microrobots [8, 26]. However, studying large-scale hydrodynamic systems is challenging as the fluid adds a large number of degrees of freedom, while the consideration of large-scale systems is unavoidable to accurately capture long-range hydrodynamics interactions and to minimize potential finite-size effects. So far, only limited studies have been conducted in this direction, particularly in three dimensions.

The goal of our current endeavor is to unravel the emergent collective behavior of systems, which combine two essential components of living and artificial active systems, self-steering and hydrodynamics. In the context of active matter, the Vicsek model [27, 28] is among the earliest and simplest models for the collective directional motion of self-propulsion and self-steering particles due to mutual alignment. While the role of alignment interactions has been extensively investigated for dry active systems [28], hydrodynamic interactions themselves have rather been viewed as a physical alignment mechanism in wet systems [29, 30]. However, hydrodynamic propulsion can in fact also destabilize polar order [31, 32, 33], which implies the necessity for an additional stabilization mechanism for alignment. To achieve stable alignment, we consider a hydrodynamic extension of the Vicsek model. Our active agents, modeled as squirmers [34, 35, 36, 37, 38], sense the propulsion direction of neighboring agents and adapt their propulsion direction accordingly by hydrodynamic self-steering [39, 40], with a “slow” temporal response due to limited maneuverability. We perform large-scale simulations in three-dimensional systems, capturing the fluid environment by the multiparticle collision dynamics (MPC) technique [41, 42, 43], a particle-resolved mesoscale hydrodynamic simulation approach.

Refer to caption
Figure 1: Illustration of the alignment interaction. (a) The microswimmer with the orientation 𝒆𝒆{\bm{e}}bold_italic_e (petrol) senses the orientations of neighboring microswimmers (purple) within the sensing range Rasubscript𝑅𝑎R_{a}italic_R start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT (green dashed circle) and reorients toward 𝒆aimsubscript𝒆aim{\bm{e}}_{\rm aim}bold_italic_e start_POSTSUBSCRIPT roman_aim end_POSTSUBSCRIPT (black arrow), which is the average orientation of the neighbors determined via Eq. (3). (b) and (c) Non-axisymmetric surface flow fields for hydrodynamic self-steering of (b) puller- and (c) pusher-type microswimmers with angular velocity 𝝎𝝎\bm{\omega}bold_italic_ω. Adapted from Ref. [40]. (d) Global polar order parameters ΨΨ\Psiroman_Ψ for pullers (purple squares), pushers (green circles), and ABPs (yellow triangles) as a function of the maneuverability ΩΩ\Omegaroman_Ω. Here, Pe=128Pe128{\rm Pe}=128roman_Pe = 128.

We observe and characterize the emergence of swarming dynamics in polar active fluids, consisting either of pusher or puller microswimmers. Our results show that hydrodynamic interactions destabilize polar order, giving rise to rich collective spatio-temporal behavior beyond the simple symmetry breaking of the dry Vicsek model. Pusher systems feature active turbulence with non-universal scaling exponents in the kinetic energy spectrum, revealing a route toward active turbulence via self-steering. Systems of pullers are accompanied by formation of dense, swarming clusters driven by hydrodynamic interactions. In particular, formation of toroidal structures are observed in the vorticity field, which are characterized by enhanced spatial vorticity-velocity cross correlations at short distances. This demonstrates that the formation of vortex rings is a direct consequence of strong active jets caused by propulsion and alignment.

II Result

II.1 Hydrodynamic Vicsek model and Polar Order

We consider a system of N𝑁Nitalic_N spherical active microswimmers with radius Rsqsubscript𝑅sqR_{\rm sq}italic_R start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT and instantaneous orientation 𝒆isubscript𝒆𝑖{\bm{e}}_{i}bold_italic_e start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT, i{1,,N}𝑖1𝑁i\in\{1,\dots,N\}italic_i ∈ { 1 , … , italic_N }. For self-propulsion and self-steering, the squirmer model is employed, where the self-propulsion is achieved via an axisymmetric surface-slip boundary condition, characterized by the speed v0subscript𝑣0v_{0}italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and the active stress β𝛽\betaitalic_β [see Materials and Methods, Eq. (4)]. Self-steering of microswimmers is modeled via an adaptive non-axisymmetric surface-slip boundary condition [see Materials and Methods, Eqs. (4) and (5)), and Fig. 1(b),(c)], which enables the squirmer to rotate and consequently reorient toward a desired direction 𝒆aim,isubscript𝒆aim𝑖{\bm{e}}_{{\rm aim},i}bold_italic_e start_POSTSUBSCRIPT roman_aim , italic_i end_POSTSUBSCRIPT with the limited angular velocity [44]

𝝎i=C0𝒆i×𝒆aim,i,subscript𝝎𝑖subscript𝐶0subscript𝒆𝑖subscript𝒆aim𝑖\displaystyle{\bm{\omega}}_{i}=C_{0}{\bm{e}}_{i}\times{\bm{e}}_{{\rm aim},i},bold_italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = italic_C start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT bold_italic_e start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT × bold_italic_e start_POSTSUBSCRIPT roman_aim , italic_i end_POSTSUBSCRIPT , (1)

where C0subscript𝐶0C_{0}italic_C start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT characterizes the strength of adaptation. Accordingly, we introduce two dimensionless parameters, the Péclet number PePe{\rm Pe}roman_Pe and the maneuverability ΩΩ\Omegaroman_Ω in the form

Pe=v0σDR,Ω=C0DR,formulae-sequencePesubscript𝑣0𝜎subscript𝐷𝑅Ωsubscript𝐶0subscript𝐷𝑅\displaystyle{\rm Pe}=\frac{v_{0}}{\sigma D_{R}},\quad\Omega=\frac{C_{0}}{D_{R% }}\,,roman_Pe = divide start_ARG italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_σ italic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT end_ARG , roman_Ω = divide start_ARG italic_C start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT end_ARG , (2)

where σ=2Rsq𝜎2subscript𝑅sq\sigma=2R_{\rm sq}italic_σ = 2 italic_R start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT and DRsubscript𝐷𝑅D_{R}italic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT are the diameter and the (thermal) rotational diffusion coefficient of a squirmer, respectively.

Refer to caption
Figure 2: Active turbulence in pusher systems. (a)-I Snapshot of squirmers for Pe=128Pe128{\rm Pe}=128roman_Pe = 128 and Ω=2048Ω2048\Omega=2048roman_Ω = 2048. The orientation of each particle is indicated by a petrol hemisphere. (a)-II Two-dimensional projection of the streamlines of the velocity field (white lines) and the magnitude of the vorticity field (heat map) for Pe=128Pe128{\rm Pe}=128roman_Pe = 128 and Ω=2048Ω2048\Omega=2048roman_Ω = 2048. (b) Equal-time spatial velocity correlation function as a function of the squirmer distance for various ΩΩ\Omegaroman_Ω as indicated in the legend. An enlarged view for small Csqsubscript𝐶sqC_{\rm sq}italic_C start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT is shown in the bottom panel. (c) Squirmer density distribution as a function of the local packing fraction for various ΩΩ\Omegaroman_Ω. The global packing fraction is represented by a black dashed line at ρloc0.093subscript𝜌loc0.093\rho_{\rm loc}\approx 0.093italic_ρ start_POSTSUBSCRIPT roman_loc end_POSTSUBSCRIPT ≈ 0.093. (d) Energy spectrum as a function of the wave number k𝑘kitalic_k for Ω/Pe=1ΩPe1\Omega/{\rm Pe}=1roman_Ω / roman_Pe = 1, 4444, 16161616, and 64646464, and Pe=32Pe32{\rm Pe}=32roman_Pe = 32, 128128128128, and 512512512512. kσ=2π/σsubscript𝑘𝜎2𝜋𝜎k_{\sigma}=2\pi/\sigmaitalic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT = 2 italic_π / italic_σ is the wave number for the squirmer diameter. (e) Energy spectrum as a function of the wave number for Pe=128Pe128{\rm Pe}=128roman_Pe = 128, Ω=512Ω512\Omega=512roman_Ω = 512, and various densities (legend); here, L/a=256𝐿𝑎256L/a=256italic_L / italic_a = 256. In (d) and (e), solid and dashed lines represent energy spectra of fluids and squirmers, respectively.

In Eq. (1), sensed information of each particle about the orientation of other neighboring particles is represented by a vector 𝒆aim,i({𝐞j})subscript𝒆aim𝑖subscript𝐞𝑗{\bm{e}}_{{\rm aim},i}(\{{\bf e}_{j}\})bold_italic_e start_POSTSUBSCRIPT roman_aim , italic_i end_POSTSUBSCRIPT ( { bold_e start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT } ). In the interest of the investigation and understanding of the generic collective behavior of an active polar fluid, we focus here on a “minimal” model, where only the swimming and steering mechanism via the adaptive surface flow field is explicitly considered, whereas sensing and information processing is taken into account implicitly by the sensed-information vector 𝒆aimsubscript𝒆aim{\bm{e}}_{\rm aim}bold_italic_e start_POSTSUBSCRIPT roman_aim end_POSTSUBSCRIPT in Eq. (1), see Ref. [40] for a more detailed discussion. As a representative example of information exchanges between intelligent microswimmers, we consider a Vicsek-type alignment interaction, where each microswimmer aims at adapting its orientation and propulsion direction to the average orientation of neighboring particles, see Fig. 1(a) for illustration. Specifically, we employ a non-additive rule of orientation adaptation, which results in non-reciprocal interactions between microswimmers, where [45]

𝒆aim,i=1NijΓi𝒆j.subscript𝒆aim𝑖1subscript𝑁𝑖subscript𝑗subscriptΓ𝑖subscript𝒆𝑗\displaystyle{\bm{e}}_{{\rm aim},i}=\frac{1}{N_{i}}\sum_{j\in\Gamma_{i}}{\bm{e% }}_{j}.bold_italic_e start_POSTSUBSCRIPT roman_aim , italic_i end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG italic_N start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_ARG ∑ start_POSTSUBSCRIPT italic_j ∈ roman_Γ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT bold_italic_e start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT . (3)

Here, ΓisubscriptΓ𝑖\Gamma_{i}roman_Γ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is the set of neighbors of the i𝑖iitalic_i-th particle in its alignment range Rasubscript𝑅𝑎R_{a}italic_R start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT, and Nisubscript𝑁𝑖N_{i}italic_N start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is the number of neighbors. As apparent from the definition in Eq. (3), 𝒆aimsubscript𝒆aim{\bm{e}}_{\rm aim}bold_italic_e start_POSTSUBSCRIPT roman_aim end_POSTSUBSCRIPT, which serves as an input signal triggering adaptive surface flows according to Eqs. (4) and (5), is typically not a unit vector. Therefore, the magnitude of the adaptation force depends on the strength of orientational order of the neighboring microswimmers. We emphasize that in our approach the steering is achieved solely via the modification of the surface flow fields, mimicking the autonomous behaviors of microorganisms, in contrast to external driving forces, see, e.g., Ref. [46]. This difference leads to very different behaviors in wet and dry systems.

For the MPC fluid dynamics simulations, we consider a MPC variant with angular momentum conservation [47], see Supplementary Materials (SM), Sec. S-I for more details. Our highly parallelized, GPU-accelerated implementation employed for the simulations is based on the framework proposed in Ref. [48]. For an accurate characterization of emergent behaviors, we consider large system sizes up to L/a=1024𝐿𝑎1024L/a=1024italic_L / italic_a = 1024, where L𝐿Litalic_L is the length of the cubic simulation box and a𝑎aitalic_a is the side length of a MPC collision cell, and up to N=884,736𝑁884736N=884,736italic_N = 884 , 736 squirmers. For the squirmers, we choose a sensing range Ra=4Rsqsubscript𝑅𝑎4subscript𝑅sqR_{a}=4R_{\rm sq}italic_R start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT = 4 italic_R start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT, and a strength of the active stress β=3𝛽3\beta=-3italic_β = - 3 and β=3𝛽3\beta=3italic_β = 3 for pushers and pullers, respectively. For most simulations, we consider the packing fraction ρ(4πRsq3/3)N/L3=0.093𝜌4𝜋superscriptsubscript𝑅sq33𝑁superscript𝐿30.093\rho\equiv(4\pi R_{\rm sq}^{3}/3)N/L^{3}=0.093italic_ρ ≡ ( 4 italic_π italic_R start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT / 3 ) italic_N / italic_L start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT = 0.093 (based on squirmer radius), or ρa(4πRa3/3)N/L36.0subscript𝜌𝑎4𝜋superscriptsubscript𝑅a33𝑁superscript𝐿36.0\rho_{a}\equiv(4\pi R_{\rm a}^{3}/3)N/L^{3}\approx 6.0italic_ρ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ≡ ( 4 italic_π italic_R start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT / 3 ) italic_N / italic_L start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ≈ 6.0 (based on sensing range), if not explicitly stated otherwise. For comparison, we also perform simulations of a dry system of aligning self-steering “intelligent” active Brownian particles (iABPs) of the same packing fraction.

Results for the global polar order parameter Ψ|i𝒆i|/NΨsubscript𝑖subscript𝒆𝑖𝑁\Psi\equiv|\sum_{i}{\bm{e}}_{i}|/Nroman_Ψ ≡ | ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_italic_e start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | / italic_N are shown in Fig. 1(d). While the systems are disordered for small maneuverability ΩΩ\Omegaroman_Ω in both dry and wet cases, a global polar order emerges only in dry systems for large ΩΩ\Omegaroman_Ω, indicating that hydrodynamic interactions destabilize the polar order. Instead of simple polar ordering, our squirmer systems show swarming dynamics with and without density modulation depending on the swimming mechanism (pusher or puller).

Refer to caption
Figure 3: Mean-square displacement of pushers. The results are shown as a function of time for various ΩΩ\Omegaroman_Ω (legend). The black dashed line indicates the ballistic dynamics of non-interacting squirmers, (Δr)2=v02(Δt)2superscriptΔ𝑟2superscriptsubscript𝑣02superscriptΔ𝑡2(\Delta r)^{2}={v_{0}}^{2}(\Delta t)^{2}( roman_Δ italic_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( roman_Δ italic_t ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, and DRsubscript𝐷𝑅D_{R}italic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT is the rotational diffusion coefficient.

II.2 Pushers: Active Turbulence via Self-Steering

Systems of rearly propelled squirmers (pushers) exhibit active turbulence, see Fig. 2(a)-I for a snapshot (also Movie S1), and their collective motion features vortical flows. Figure 2(a)-II displays the fluid velocity field reflecting vortical structures and fluctuations in the magnitude of the vorticity. Accordingly, equal-time spatial velocity correlations become negative at large distances (less than L/2𝐿2L/2italic_L / 2) for large ΩΩ\Omegaroman_Ω, see Fig. 2(b). However, no pronounced density fluctuations are visible, see Fig. 2(a)-I. This is confirmed by the analysis of the local density distribution via from Voronoi tessellation [49]. As shown in Fig. 2(c), the distribution exhibits a peak near the global squirmer density independent of ΩΩ\Omegaroman_Ω.

To characterize the dynamics, we examine the kinetic energy spectrum as a representative indicator of active turbulence [15]. We determine the energy spectra for both the squirmer and the fluid. For the squirmers, we first calculate spatial velocity correlation functions and then perform a Fourier transform to obtain the energy spectrum. For the fluid, we calculate the energy spectrum directly from the (Eulerian) velocity field.

Before proceeding to a more detailed analysis of systems of self-steering squirmers, we briefly discuss as reference case the dynamics of pushers without self-steering, i.e., squirmers with Ω=0Ω0\Omega=0roman_Ω = 0. Interestingly, we do not observe any significant collective behavior without self-steering for Pe=128Pe128{\rm Pe}=128roman_Pe = 128 as well as 256256256256, with β=3𝛽3\beta=-3italic_β = - 3 and β=5𝛽5\beta=-5italic_β = - 5, in our simulations (see SM, Fig. S1). It seems that thermal fluctuations and consequently rotational diffusion of the microswimmers suppress any structure formation at the considered small packing fraction, compare also Ref. [50] for a related lattice Boltzmann simulation study of pushers. Moreover, the energy spectrum of squirmers strongly deviates from the fluid energy spectrum for large and intermediate k𝑘kitalic_k (Fig. S1), indicating that the fluid flows only mildly affect the squirmer dynamics. As also pointed out in Ref. [23], presumably even higher densities of force-dipoles are necessary to promote the emergence of collective motion engaging both the fluid and microswimmers of pushers without self-steering, solely based on hydrodynamic alignment. Indeed, the number density explored in Refs. [51, 52] is typically one or two orders of magnitude larger than in our case due to different geometry of active particles and their use of force dipoles.

Refer to caption
Figure 4: Swarming dynamics in puller systems. (a)-I Snapshot of squirmers with their orientation indicated by a petrol hemisphere. (a)-II Two-dimensional projection of the streamlines of the velocity field (white lines) and the magnitude of the vorticity field (heat map). The parameters are Pe=128Pe128{\rm Pe}=128roman_Pe = 128 and Ω=8192Ω8192\Omega=8192roman_Ω = 8192. (b) Squirmer density distribution as a function of the local packing fraction for various ΩΩ\Omegaroman_Ω. (c) Mean-square displacement as a function of time for various ΩΩ\Omegaroman_Ω (legend). The black dashed line indicates the ballistic dynamics of non-interacting squirmers, (Δr)2=v02(Δt)2superscriptΔ𝑟2superscriptsubscript𝑣02superscriptΔ𝑡2(\Delta r)^{2}=v_{0}^{2}(\Delta t)^{2}( roman_Δ italic_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( roman_Δ italic_t ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (d) Alignment distribution as a function of scalar product of a squirmer’s orientation and velocity vector. Data for pullers (solid lines) as well as pusher (dashed line) are displayed for comparison. (e) Energy spectrum as a function of the wave number k𝑘kitalic_k for various ΩΩ\Omegaroman_Ω (legend). Energy spectra of fluids and squirmers are represented by solid and dashed lines, respectively. Various power-laws are indicated by short black lines. In all figures Pe=128Pe128{\rm Pe}=128roman_Pe = 128.

In sharp contrast, systems of self-steering squirmers display pronounced self-organization, as demonstrated in Fig. 2. With increasing maneuverability ΩΩ\Omegaroman_Ω, a scaling regime emerges in the fluid power spectrum for wave numbers k/kσ0.3less-than-or-similar-to𝑘subscript𝑘𝜎0.3k/k_{\sigma}\lesssim 0.3italic_k / italic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ≲ 0.3, see Fig. 2(d). Corresponding energy spectra of the squirmer motion exhibit the same scaling behavior with the same exponents – with scaling extending even to lower wave numbers k/kσ0.02𝑘subscript𝑘𝜎0.02k/k_{\sigma}\approx 0.02italic_k / italic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ≈ 0.02 for Ω=8192Ω8192\Omega=8192roman_Ω = 8192 – confirming the emergence of active turbulence, where the dynamics of squirmers and fluid are strongly correlated on larger length scales. The exponents ν𝜈\nuitalic_ν of the energy spectrum, E(k)(k/kσ)νsimilar-to𝐸𝑘superscript𝑘subscript𝑘𝜎𝜈E(k)\sim(k/k_{\sigma})^{-\nu}italic_E ( italic_k ) ∼ ( italic_k / italic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT - italic_ν end_POSTSUPERSCRIPT are found to be non-universal, depending on and increasing with ΩΩ\Omegaroman_Ω, roughly in the range 2.8ν4.0less-than-or-similar-to2.8𝜈less-than-or-similar-to4.02.8\lesssim\nu\lesssim 4.02.8 ≲ italic_ν ≲ 4.0 for 4Pe/Ω644PeΩ644\leq{\rm Pe}/\Omega\leq 644 ≤ roman_Pe / roman_Ω ≤ 64, see Fig. S2. Moreover, with increase ΩΩ\Omegaroman_Ω, not only ν𝜈\nuitalic_ν increases, but also the scaling regimes extends to large length scales (smaller k𝑘kitalic_k, as more squirmers participate in the self-organized vortex structures, and the peak height in the energy spectrum increases (up to |E|200v02𝐸200superscriptsubscript𝑣02|E|\approx 200v_{0}^{2}| italic_E | ≈ 200 italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT for Ω=8192Ω8192\Omega=8192roman_Ω = 8192). This indicates that squirmers attain much higher velocities, which is quantitatively confirmed by the mean-square displacement (MSD), see Fig. 3. In the ballistic regime, the MSD (Δr)2v2ttsimilar-todelimited-⟨⟩superscriptΔ𝑟2superscript𝑣2superscript𝑡𝑡\langle(\Delta r)^{2}\rangle\sim v^{2}t^{t}⟨ ( roman_Δ italic_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ⟩ ∼ italic_v start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_t start_POSTSUPERSCRIPT italic_t end_POSTSUPERSCRIPT, hence, the increasing amplitude in Fig. 3 is a direct measure of the increasing speed as ΩΩ\Omegaroman_Ω increases for Ω512Ω512\Omega\geq 512roman_Ω ≥ 512, i.e., in the regime where a broad scaling region can be identified in the energy spectrum (Fig. 2(d)). The Péclet number PePe{\rm Pe}roman_Pe merely affects turbulent dynamics on large length scales, but the ratio Pe/ΩPeΩ{\rm Pe}/\Omegaroman_Pe / roman_Ω determines the scaling exponent ν𝜈\nuitalic_ν as shown in Fig. S2. We also probe density effects varying the number of squirmers. We first notice that the fluid energy spectra seem to converge for high densities, as shown in Fig. 2(e), in accordance with the previous report [52]. As the packing fraction ρ𝜌\rhoitalic_ρ decreases from 0.1860.1860.1860.186 to 0.0120.0120.0120.012 corresponding to ρa=11.9subscript𝜌𝑎11.9\rho_{a}=11.9italic_ρ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT = 11.9 and 0.7680.7680.7680.768, downward shifts in the fluid energy spectra are observed, indicating that fluid stirring by squirmers are not strong at low densities. Consequently, the scaling regime in the energy spectra of squirmers shrinks as the density decreases, e.g., to a narrow range of 0.2k/kσ0.4less-than-or-similar-to0.2𝑘subscript𝑘𝜎less-than-or-similar-to0.40.2\lesssim k/k_{\sigma}\lesssim 0.40.2 ≲ italic_k / italic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ≲ 0.4 at ρ=0.012𝜌0.012\rho=0.012italic_ρ = 0.012.

For k/kσ0.5greater-than-or-equivalent-to𝑘subscript𝑘𝜎0.5k/k_{\sigma}\gtrsim 0.5italic_k / italic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ≳ 0.5, or length scales smaller than about 2σ2𝜎2\sigma2 italic_σ, differences between the fluid and squirmer energy spectra appear. In this regime, near-field hydrodynamic flows play a significant role. Moreover, the properties of the energy spectra depend on PePe{\rm Pe}roman_Pe, indicating that noise effects are significant at these small length scales, see Fig. S3 for fluid thermal energy spectra [53]. Therefore, in wet systems and small densities, a strong self-steering of active particles is crucial for the emergence of a large-scale coherent collective motion. Otherwise disorder and hydrodynamic instabilities on small length-scales may prevail [31]. Such an observation should also apply to systems where the alignment is mediated via steric repulsion of elongated body shapes [15] or strong hydrodynamic force dipoles [52]. In any cases, in active turbulence of wet systems, collective fluid flows induced by microswimmers build up the collective behavior of microswimmers on large length-scales with fast dynamics. In sharp contrast, active turbulence in dry systems [54, 55] requires densely packed active particles, because a speedup mechanism as in wet systems is lacking as steric repulsions may only result in a slow-down of active particles. Instead, in dry active turbulence, the scaling regime develops in large k𝑘kitalic_k, or equivalently, small |E|𝐸|E|| italic_E | regimes due to chaotic interparticle collisions on small length scales.

Refer to caption
Figure 5: Formation of vortex rings. (a) System cutout (100z/a300100𝑧𝑎300100\leq z/a\leq 300100 ≤ italic_z / italic_a ≤ 300) of the surface plot for the vorticity field with Pe=128Pe128{\rm Pe}=128roman_Pe = 128, Ω=2048Ω2048\Omega=2048roman_Ω = 2048, L/a=768𝐿𝑎768L/a=768italic_L / italic_a = 768, and ω0.33ωmax𝜔0.33subscript𝜔max\omega\approx 0.33\omega_{\rm max}italic_ω ≈ 0.33 italic_ω start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, where ωmaxsubscript𝜔max\omega_{\rm max}italic_ω start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT indicates the maximum value of ω𝜔\omegaitalic_ω in the system at this particular time step. (b) Quasi two-dimensional slice of the vorticity field, projected onto the corresponding two dimensions (white lines), together with the magnitude of the vorticity field (heat maps) for Pe=128,Ω=8192formulae-sequencePe128Ω8192{\rm Pe}=128,\Omega=8192roman_Pe = 128 , roman_Ω = 8192, and L/a=1024𝐿𝑎1024L/a=1024italic_L / italic_a = 1024. (c) Time evolution of squirmer configurations during the emergence and destruction of a vortex ring for I t=t0𝑡subscript𝑡0t=t_{0}italic_t = italic_t start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT (formation of cluster), II t1=t0+0.041/DRsubscript𝑡1subscript𝑡00.041subscript𝐷𝑅t_{1}=t_{0}+0.041/D_{R}italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_t start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + 0.041 / italic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT (formation of jet flow), III t2=t0+0.10/DRsubscript𝑡2subscript𝑡00.10subscript𝐷𝑅t_{2}=t_{0}+0.10/D_{R}italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_t start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + 0.10 / italic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT (emergence of vortex ring), and IV t3=t0+0.13/DRsubscript𝑡3subscript𝑡00.13subscript𝐷𝑅t_{3}=t_{0}+0.13/D_{R}italic_t start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = italic_t start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + 0.13 / italic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT (jelly-fish like spreading of squirmers). (d) Vortex ring (blue torus, ω0.4ωmax𝜔0.4subscript𝜔max\omega\approx 0.4\omega_{\rm max}italic_ω ≈ 0.4 italic_ω start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT) and fluid velocity field (tubes) together with the squirmer positions (red bullets) and orientations (indicated by red bars) at t=t2𝑡subscript𝑡2t=t_{2}italic_t = italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT. I Top view and II side-view of the vortex ring, which depict the cluster of squirmers on the bottom-left pulling fluid on the top-right toward the cluster. Squirmers in the frontal region of the cluster are moving along the outer surface of the vortex ring, see also (c)-III to identify the locations of the vortex ring and squirmers. (e)-I Surface plot of the jet flow (|𝒗fl|1.6v0subscript𝒗fl1.6subscript𝑣0|{\bm{v}}_{\rm fl}|\approx 1.6v_{0}| bold_italic_v start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT | ≈ 1.6 italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT) generated at t=t1𝑡subscript𝑡1t=t_{1}italic_t = italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT, together with the vortex ring subsequently formed at t=t2𝑡subscript𝑡2t=t_{2}italic_t = italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT. Note that the jet flow is formed at the frontal region of the aligned cluster ((c)-II), and directed toward the cluster. (e)-II Positions and velocities of the squirmers (bullets) at t=t3𝑡subscript𝑡3t=t_{3}italic_t = italic_t start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT ((c)-IV), together with the vortex ring previously formed at t=t2𝑡subscript𝑡2t=t_{2}italic_t = italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT. The colors of the squirmers indicate the inner product of their orientation and propulsion direction. In (c), (d), and (e), Pe=128Pe128{\rm Pe}=128roman_Pe = 128, Ω=2048Ω2048\Omega=2048roman_Ω = 2048.

II.3 Pullers: Swarming Dynamics via Self-Steering

In a systems of self-steering pullers, a rich swarming dynamics develops, as shown in Fig. 4. The self-organization is characterized by the formation of morphologically complex clusters of microswimmers, which, on larger length scales, exhibit visually chaotic movements and exchange constitutive squirmers with each others, see Fig. 4(a)-I (also Movie S2). Still, the puller system exhibits a velocity field with vortical structures (see Fig. 4(a)-II), surprisingly similar to pushers.

Refer to caption
Figure 6: Velocity-vorticity coupling. (a) Vorticity field of pullers (white lines) combined with the magnitude of the velocity field (heat maps) (see Fig. 5(b)). (b) Distribution of the Cartesian components of the vorticity field ω¯αsubscript¯𝜔𝛼\bar{\omega}_{\alpha}over¯ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT as a function of the vorticity ω𝜔\omegaitalic_ω. The black solid line shows the Gaussian distribution. (c) Cross correlation function between the velocity and the vorticity field as a function of the radial squirmer distance.

Density Modulation.

The local density distribution is calculated by employing a Voronoi tessellation [49]. In sharp contrast to pusher systems (Fig. 2(c)), puller systems exhibit significant density modulations, see Fig. 4(b). For large ΩΩ\Omegaroman_Ω, Ω4096greater-than-or-equivalent-toΩ4096\Omega\gtrsim 4096roman_Ω ≳ 4096, the density distribution is broad, while for Ω=2048Ω2048\Omega=2048roman_Ω = 2048, a clear tendency of segregation is observed (Movie S3) with a low-density peak at ρloc0subscript𝜌loc0\rho_{\rm loc}\approx 0italic_ρ start_POSTSUBSCRIPT roman_loc end_POSTSUBSCRIPT ≈ 0, and a high-density peak at ρloc0.5subscript𝜌loc0.5\rho_{\rm loc}\approx 0.5italic_ρ start_POSTSUBSCRIPT roman_loc end_POSTSUBSCRIPT ≈ 0.5. For small ΩΩ\Omegaroman_Ω in the range 128Ω512128Ω512128\leq\Omega\leq 512128 ≤ roman_Ω ≤ 512, non-mobile clusters appear (see Fig. S4). We focus here on the dynamic clusters, as formation of a dense static cluster may be affected by a depletion of MPC fluid particles inside the cluster – which is related to the (weak) compressibility of the MPC fluid, and artificially enhance cluster stability [56]. For Ω32Ω32\Omega\leq 32roman_Ω ≤ 32, unimodal distributions are recovered, but with a peak at a density smaller than the global density of ρ=0.093𝜌0.093\rho=0.093italic_ρ = 0.093, and with fatter tails than those in pusher systems, which reflects a clustering tendency reported for pullers [56, 57].

Decoupling of Puller Orientation and Velocity.

The emergence of a high-density regime for Ω2048Ω2048\Omega\geq 2048roman_Ω ≥ 2048 needs to be distinguished from motility-induced phase separation (MIPS) in dry systems. First, squirmers in a (small) cluster exhibit local polar order, giving rise to a coherent directional motion of the cluster. Second, the pullers actually swim faster on average than their bare self-propulsion speed v0subscript𝑣0v_{0}italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, as shown in Fig. 4(c) – as in the system of pushers (see Fig. 2(e)). In particular, decoupling of self-propulsion and velocity of squirmers is far more significant than in MIPS. Even situations, where pullers are driven backward occur frequently (Movie S4). Indeed, as shown in Fig. 4(d), the probability that the velocity of a squirmer is anti-parallel to its orientation, i.e., 𝒗/|𝒗|=𝒆𝒗𝒗𝒆{\bm{v}}/|{\bm{v}}|=-{\bm{e}}bold_italic_v / | bold_italic_v | = - bold_italic_e, is even higher than for the parallel case, i.e., 𝒗/|𝒗|=𝒆𝒗𝒗𝒆{\bm{v}}/|{\bm{v}}|={\bm{e}}bold_italic_v / | bold_italic_v | = bold_italic_e. Therefore, we conclude that hydrodynamic interactions between self-steering pullers dominate over the self-propulsion forces. The attractive hydrodynamic interactions between aligned pullers in a head-to-tail configuration promotes the formation of dense clusters. Yet, as we will demonstrate below, dense clusters are not static, but exhibit a highly dynamical morphology.

Decoupling of Puller and Fluid Dynamics.

The kinetic energy spectra E(k)𝐸𝑘E(k)italic_E ( italic_k ) for squirmer motion and the fluid in puller systems are presented in Fig. 3. While for the swarming of pullers, E(k)𝐸𝑘E(k)italic_E ( italic_k ) exhibits a power-law decay in the intermediate wave-number regime (0.04<k/kσ<0.20.04𝑘subscript𝑘𝜎0.20.04<k/k_{\sigma}<0.20.04 < italic_k / italic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT < 0.2) for the fluid, as in pusher systems, several features significantly deviate from those of pushers, which again demonstrates the uniqueness of puller swarming. Most of all, a pronounced mismatch between squirmer and fluid energy spectra is observed in the intermediate regime of k/kσ0.1𝑘subscript𝑘𝜎0.1k/k_{\sigma}\approx 0.1italic_k / italic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ≈ 0.1 (see Fig. S5(a) and Fig. S7(a) for the corresponding correlation functions), which indicates that pullers are not simply driven by the fluid flow. The exponents of the three scaling regimes in Fig. 4(e) show universal behavior for the investigated range of ΩΩ\Omegaroman_Ω. Notably, in the case of the squirmer energy spectra for 0.2<k/kσ<0.60.2𝑘subscript𝑘𝜎0.60.2<k/k_{\sigma}<0.60.2 < italic_k / italic_k start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT < 0.6, which correspond the length scales of the high-density intra-cluster regime where steric repulsion and near-field hydrodynamic interactions dominate, the obtained value of the exponent ν5/3similar-to-or-equals𝜈53\nu\simeq 5/3italic_ν ≃ 5 / 3 is in accordance with the Kolmogorov exponent for classical hydrodynamic turbulence and with active turbulence of dense quasi-two-dimensional pusher systems [15]. For larger length scales, an exponent ν11/3similar-to-or-equals𝜈113\nu\simeq 11/3italic_ν ≃ 11 / 3 is obtained, which persists under variation of global densities, see Fig. S6. A similar value has been found previously in Lattice Boltzmann simulations of extended force dipoles [52], although at higher densities and for pushers.

Formation of vortex ring.

An even more striking feature is observed in configurations of the vorticity field. Visual observation of the time evolution of the system indicates a typical dynamical behavior in the morphology of clusters, which involves the pulsatile transformation of aggregates from a spherical shape into jellyfish-like arrangements. In terms of fluid mechanics, the jellyfish-like morphology suggests formation of a vortex ring, which is indeed confirmed by emergence of toroidal structures in the vorticity fields extracted from the simulations, see Fig. 5(a). As shown in Fig. 5(b), the vorticity field exhibits whirling patterns within regions where the magnitude of the vorticity field is large.

For a detailed illustration, we consider the small system size L/a=128𝐿𝑎128L/a=128italic_L / italic_a = 128, where only a single cluster emerges (Movie S5). As shown in Fig. 5(c)-I, the dynamics initiates with formation of a cluster. Then, due to alignment, squirmers rotate, and a polar order emerges within the cluster, see Fig. 5(c)-II, and also Fig. S5(b) for the corresponding order parameter. Such an ordered structure gives rise to a strong collective fluid flow, which generates a pronounced jet in front of the cluster, see the yellow surface in Fig. 5(e)-I. Notably, the jet flow is self-generated via active stirring of microswimmers in this case, instead of external perturbation as in passive hydrodynamic fluids. Subsequently, a spread-out motion of squirmers is initiated (Fig. 5(c)-III). Simultaneously, a vortex ring is formed around the cluster (blue ring in Fig. 5(e)-I), while squirmers are moving forward. As shown in Fig. 5, the velocity field is indeed wrapping around the vortex ring, in accordance with Fig. 5(b). Then, the pullers continue to spread out, rolling about the region where the vortex ring forms, see Fig. 5(c)-IV. While swimmers at the cluster center swim forward, they are dragged backward at the periphery, as shown in Fig. 5(e)-II, which contributes to the anomalous behavior in the distribution of 𝒗𝒆/|𝒗|𝒗𝒆𝒗{\bm{v}}\cdot{\bm{e}}/|{\bm{v}}|bold_italic_v ⋅ bold_italic_e / | bold_italic_v | (see Fig. 4(d)). Eventually, the cluster dissolves, ending its life cycle.

II.4 Velocity-Vorticity Coupling

The sequential time evolution described so far indicates a strong coupling between the fluid velocity and the vorticity field in a puller cluster. As shown in Fig. 6(a), the rotation of the vorticity field for pullers is indeed centered at regions with a strong velocity field.

For a more quantitative characterization, we examine the distribution of the Cartesian components of the velocity and vorticity fields. In pusher systems, both the velocity and vorticity fields (Fig. S7(b) and Fig. 6(b), respectively) exhibit a Gaussian distribution, which is an indicator of active turbulence [15]. For pullers, both deviate from a Gaussian, demonstrating that the swarming dynamics of pullers is not active turbulence. Specifically, the velocity distribution of pullers exhibits “fat” exponential tails, as shown in Fig. S7(b), in line with the emergence of stronger fluid flows than “expected”, i.e., the occurrence of jet plumes induced by aligned pullers. Also the vorticity field for pullers shows a broader distribution than that of pushers as displayed in Fig. 6(b). Furthermore, the vorticity distribution for pullers in Fig. 6(b) shows a rather sharp peak at ω¯αsubscript¯𝜔𝛼\bar{\omega}_{\alpha}over¯ start_ARG italic_ω end_ARG start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT, the average of three Cartesian components ωαsubscript𝜔𝛼\omega_{\alpha}italic_ω start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT (α=x,y,z𝛼𝑥𝑦𝑧\alpha=x,y,zitalic_α = italic_x , italic_y , italic_z) of 𝝎𝝎\bm{\omega}bold_italic_ω, indicating a weak separation between regions with strong and weak vorticity.

A more fundamental difference in velocity-vorticity coupling is revealed by a cross correlations between the magnitude of vorticity and that of the velocity field (|ω(𝒓)|𝜔𝒓|\omega(\bm{r})|| italic_ω ( bold_italic_r ) | and |𝒗(𝒓)|𝒗𝒓|{\bm{v}}(\bm{r})|| bold_italic_v ( bold_italic_r ) |), as defined in Eq. (15). In Fig. 6(c), the various cross correlations for pushers exhibit a pronounced peak at r/a=20𝑟𝑎20r/a=20italic_r / italic_a = 20, 60606060, and 120120120120 for Ω=512Ω512\Omega=512roman_Ω = 512, 2048204820482048, and 8192819281928192, respectively, while the cross correlation at r=0𝑟0r=0italic_r = 0 is not strong. Hence, for pushers, the velocity field of a vortex is weak at the center but strong at intermediate regions between the center and periphery. In sharp contrast, for pullers, the cross correlation between vorticity and velocity fields is already strong at small distances, which demonstrates that a strong velocity field generates a strong vorticity in the immediate vicinity of the jet flow. Moreover, the cross correlation decays faster than for pushers, assuming negative values before approaching zero.

III Discussion and Conclusions

We have studied self-organization and dynamics in three-dimensional wet systems of self-steering squirmers, which aim for alignment of their orientation with that of their neighbors. We demonstrate that alignment via hydrodynamic self-steering gives rise to a rich collective behavior in such polar active fluids, depending on the type of active stresses, i.e., whether the microswimmers are pushers or pullers. In both cases, an essential role of hydrodynamics is the breaking of long-range polar order, which causes the emergence of chaotic behavior.

For pushers, the particle distribution is quite homogeneous, the distribution of the Cartesian velocity components is Gaussian, and the kinetic energy spectrum displays a peak and subsequent power-law decay with increasing wave vector, which indicates active-turbulent behavior. An intriguing feature is that strong self-steering enhances the coherent movement of microswimmers, which leads to collective particle motion with speeds much faster then the individual swim speed – as can be seen in the increasing magnitude of the peak in the energy spectrum, combined with an extension of the scaling regime toward large length scales. This implies that large-scale flows are induced by the collective motion, which drag the microswimmers along and supersede their individual motion. Thus, the polarity field and the fluid flow field are strongly coupled.

For pullers, another type of self-organization emerges, which is strictly distinguished from motility-induced phase separation (MIPS) of dry ABP systems as well as active turbulence of pushers. The particle density is now found to be very inhomogeneous, as the pullers tend to form clusters. However, these clusters are not static, but appear as quite unstable. The particle alignment inside the cluster generates a strong fluid jet and a vortex ring, which pulls apart the cluster and leads to its disintegration. These strong flows imply that the probability of fast fluid flows is enhanced, which is reflected in the emergence of fat tails in the velocity distribution.

In contrast, wet systems of self-propelled squirmers without self-steering display no interesting collective behavior in three dimensions, not even at high squirmer volume fractions.

Our numerical observations for ensembles of self-steering pullers challenges current theoretical views on collective behaviors in wet active systems. So far, it has been typically assumed that the polarity and velocity fields of active fluids are essentially identical, based on the assumption of a nearly homogeneous distribution of active particles. Heterogeneous densities have been observed recently in models of compressible polar active fluids for bacterial suspensions [58, 59]; however, the mechanism is entirely different in this case, as hydrodynamic interactions are not considered, and clustering is driven by a strong dependence of self-propulsion speed on the local density.

Our results demonstrate that in wet systems of self-steering microswimmers in three dimensions, the interplay of the particle density and polarity, and of the fluid velocity field can give rise to a surprisingly rich variety of emergent behaviors – already for a highly simplified model systems with only a single particle type.

IV Materials and Methods

IV.1 Mesoscale Fluid Model: Multi-Particle Collision Dynamics

We adopt the multiparticle collision dynamics (MPC) method [42, 43], a particle-based mesoscale simulation approach, as model for the fluid. Specifically, we employ the stochastic rotation variant of MPC with angular momentum conservation (MPC-SRD+a) [47] and the cell-level Maxwell-Boltzmann scaling thermostat [60]. The algorithm consists of alternating streaming and collision steps. In the streaming step, the MPC point particles of mass m𝑚mitalic_m propagate ballistically over the collision time interval hhitalic_h, denoted as collision time. In the collision step, fluid particles are sorted into the cells of a cubic lattice of lattice constant a𝑎aitalic_a defining the collision environment. Then, their relative velocities, with respect to the center-of-mass velocity of the collision cell, are rotated around a randomly oriented axes by a fixed angle α𝛼\alphaitalic_α. The algorithm conserves mass, linear, and angular momentum on the collision-cell level, while thermal fluctuations are automatically incorporated. More details are described in Refs. [38, 15] and the SM, which refers to Refs. [61, 62, 63, 64] additionally. Our GPU-based, highly parallelized implementation employed for the simulations is proposed in Ref. [48].

We use the average MPC fluid density (particles per collision cell) Nc=20delimited-⟨⟩subscript𝑁𝑐20\langle N_{c}\rangle=20⟨ italic_N start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ⟩ = 20, the collision time h=0.02am/(kBT)0.02𝑎𝑚subscript𝑘B𝑇h=0.02a\sqrt{m/(k_{\rm B}T)}italic_h = 0.02 italic_a square-root start_ARG italic_m / ( italic_k start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT italic_T ) end_ARG and the rotation angle α=130𝛼superscript130\alpha=130^{\circ}italic_α = 130 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT, which yield the fluid viscosity of η=42.6mkBT/a2𝜂42.6𝑚subscript𝑘B𝑇superscript𝑎2\eta=42.6\sqrt{mk_{\rm B}T}/a^{2}italic_η = 42.6 square-root start_ARG italic_m italic_k start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT italic_T end_ARG / italic_a start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [47, 65]. With the squirmer radius Rc=3asubscript𝑅c3𝑎R_{\rm c}=3aitalic_R start_POSTSUBSCRIPT roman_c end_POSTSUBSCRIPT = 3 italic_a, these MPC parameters yield the rotational diffusion coefficient DR=4.1×105kBT/m/asubscript𝐷𝑅4.1superscript105subscript𝑘B𝑇𝑚𝑎D_{R}=4.1\times 10^{-5}\sqrt{k_{\rm B}T/m}/aitalic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT = 4.1 × 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT square-root start_ARG italic_k start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT italic_T / italic_m end_ARG / italic_a.

IV.2 Hydrodynamic Self-Steering

Hydrodynamic self-steering of squirmers is achieved via adaptive surface flow fields, given as

uϕ=subscript𝑢italic-ϕabsent\displaystyle u_{\phi}=italic_u start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT = 32v0sinθ(1+βcosθ)1Rsq2(C~11cosϕC11sinϕ)32subscript𝑣0𝜃1𝛽𝜃1superscriptsubscript𝑅sq2subscript~𝐶11italic-ϕsubscript𝐶11italic-ϕ\displaystyle\frac{3}{2}v_{0}\sin{\theta}(1+\beta\cos{\theta})-\frac{1}{R_{\rm sq% }^{2}}(\tilde{C}_{11}\cos{\phi}-C_{11}\sin{\phi})divide start_ARG 3 end_ARG start_ARG 2 end_ARG italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT roman_sin italic_θ ( 1 + italic_β roman_cos italic_θ ) - divide start_ARG 1 end_ARG start_ARG italic_R start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ( over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT roman_cos italic_ϕ - italic_C start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT roman_sin italic_ϕ ) (4)
uθ=subscript𝑢𝜃absent\displaystyle u_{\theta}=italic_u start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT = cosθRsq2(C11cosϕ+C~11sinϕ),𝜃superscriptsubscript𝑅sq2subscript𝐶11italic-ϕsubscript~𝐶11italic-ϕ\displaystyle\frac{\cos{\theta}}{R_{\rm sq}^{2}}(C_{11}\cos{\phi}+\tilde{C}_{1% 1}\sin{\phi}),divide start_ARG roman_cos italic_θ end_ARG start_ARG italic_R start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ( italic_C start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT roman_cos italic_ϕ + over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT roman_sin italic_ϕ ) , (5)

where θ𝜃\thetaitalic_θ and ϕitalic-ϕ\phiitalic_ϕ are the polar and azimuthal angles in a body-fixed reference frame. The parameter v0subscript𝑣0v_{0}italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT characterizes the swim speed and β𝛽\betaitalic_β the strength of the force dipole, where β<0𝛽0\beta<0italic_β < 0 for pullers and β>0𝛽0\beta>0italic_β > 0 for pushers [21, 66]. C11subscript𝐶11C_{11}italic_C start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT and C~11subscript~𝐶11\tilde{C}_{11}over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT control the magnitudes of non-axisymmetric surface-flow components, leading to rotational motion of the body [39, 40]. Specifically, the non-axisymmetric flow fields enable the adaptive motion of Eq. (1) via

C11subscript𝐶11\displaystyle C_{11}italic_C start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT =C0Rsq3(𝒆×𝒆aim)𝒆x,absentsubscript𝐶0superscriptsubscript𝑅sq3𝒆subscript𝒆aimsubscript𝒆𝑥\displaystyle=C_{0}R_{\rm sq}^{3}({\bm{e}}\times{\bm{e}}_{\rm aim})\cdot{\bm{e% }}_{x},= italic_C start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_R start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ( bold_italic_e × bold_italic_e start_POSTSUBSCRIPT roman_aim end_POSTSUBSCRIPT ) ⋅ bold_italic_e start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT , (6)
C~11subscript~𝐶11\displaystyle\tilde{C}_{11}over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT =C0Rsq3(𝒆×𝒆aim)𝒆y,absentsubscript𝐶0superscriptsubscript𝑅sq3𝒆subscript𝒆aimsubscript𝒆𝑦\displaystyle=C_{0}R_{\rm sq}^{3}({\bm{e}}\times{\bm{e}}_{\rm aim})\cdot{\bm{e% }}_{y},= italic_C start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_R start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ( bold_italic_e × bold_italic_e start_POSTSUBSCRIPT roman_aim end_POSTSUBSCRIPT ) ⋅ bold_italic_e start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT , (7)

where 𝒆xsubscript𝒆𝑥{\bm{e}}_{x}bold_italic_e start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT and 𝒆ysubscript𝒆𝑦{\bm{e}}_{y}bold_italic_e start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT are the unit vectors along the axis of the body-fixed reference frame. For the self-propulsion, we use v0/kBT/m=0.007872subscript𝑣0subscript𝑘B𝑇𝑚0.007872v_{0}/\sqrt{k_{\rm B}T/m}=0.007872italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT / square-root start_ARG italic_k start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT italic_T / italic_m end_ARG = 0.007872 and 0.0314880.0314880.0314880.031488, which correspond to Pe=v0/(σDR)=32Pesubscript𝑣0𝜎subscript𝐷𝑅32{\rm Pe}=v_{0}/(\sigma D_{R})=32roman_Pe = italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT / ( italic_σ italic_D start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT ) = 32 and 128128128128, and Re=0.022Re0.022{\rm Re}=0.022roman_Re = 0.022 and 0.0890.0890.0890.089, respectively. The values of the self-steering strength C0subscript𝐶0C_{0}italic_C start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT are varied from 00 to 0.335872kBT/m/a0.335872subscript𝑘B𝑇𝑚𝑎0.335872\sqrt{k_{\rm B}T/m}/a0.335872 square-root start_ARG italic_k start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT italic_T / italic_m end_ARG / italic_a yielding 0Ω81920Ω81920\leq\Omega\leq 81920 ≤ roman_Ω ≤ 8192.

IV.3 Steric squirmer interaction

Steric repulsion between two squirmers is described by the separation-shifted Lennard-Jones potential

ULJ(ds)=4ϵ0[(σ0ds+σ0)12(σ0ds+σ0)6+14],subscript𝑈LJsubscript𝑑𝑠4subscriptitalic-ϵ0delimited-[]superscriptsubscript𝜎0subscript𝑑𝑠subscript𝜎012superscriptsubscript𝜎0subscript𝑑𝑠subscript𝜎0614\displaystyle U_{\rm LJ}(d_{s})=4\epsilon_{0}\left[\left(\frac{\sigma_{0}}{d_{% s}+\sigma_{0}}\right)^{12}-\left(\frac{\sigma_{0}}{d_{s}+\sigma_{0}}\right)^{6% }+\frac{1}{4}\right],italic_U start_POSTSUBSCRIPT roman_LJ end_POSTSUBSCRIPT ( italic_d start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ) = 4 italic_ϵ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT [ ( divide start_ARG italic_σ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_d start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT + italic_σ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT - ( divide start_ARG italic_σ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_d start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT + italic_σ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT + divide start_ARG 1 end_ARG start_ARG 4 end_ARG ] , (8)

for ds<(21/61)σ0subscript𝑑𝑠superscript2161subscript𝜎0d_{s}<(2^{1/6}-1)\sigma_{0}italic_d start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT < ( 2 start_POSTSUPERSCRIPT 1 / 6 end_POSTSUPERSCRIPT - 1 ) italic_σ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and zero otherwise, where dssubscript𝑑𝑠d_{s}italic_d start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT indicates the surface-to-surface distance between the two squirmer. To avoid loss of hydrodynamic interactions when two squirmers contact with each other, we also include a virtual safety distance dvsubscript𝑑𝑣d_{v}italic_d start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT [67, 15], which leads to the effective distance ds=rcσ2dvsubscript𝑑𝑠subscript𝑟𝑐𝜎2subscript𝑑𝑣d_{s}=r_{c}-\sigma-2d_{v}italic_d start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = italic_r start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT - italic_σ - 2 italic_d start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT, where rcsubscript𝑟𝑐r_{c}italic_r start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT denotes the center-to-center distance and σ𝜎\sigmaitalic_σ is the squirmer diameter. We choose σ0=2dvsubscript𝜎02subscript𝑑𝑣\sigma_{0}=2d_{v}italic_σ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 2 italic_d start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT. Numerically, the equations of motion for the rigid-body dynamics of the squirmers are solved by the velocity-Verlet algorithm, see SM for more details.

IV.4 Spatial correlation and Energy Spectrum

Squirmers: Particle-based approach.

The spatial velocity correlation function of the squirmers is defined as

Csq(r)=1v¯2ij𝒗i𝒗jδ(r|𝒓i𝒓j|)ijδ(r|𝒓i𝒓j|).subscript𝐶sq𝑟1superscript¯𝑣2delimited-⟨⟩subscript𝑖𝑗subscript𝒗𝑖subscript𝒗𝑗𝛿𝑟subscript𝒓𝑖subscript𝒓𝑗delimited-⟨⟩subscript𝑖𝑗𝛿𝑟subscript𝒓𝑖subscript𝒓𝑗\displaystyle C_{\rm sq}(r)=\frac{1}{\bar{v}^{2}}\frac{\left\langle\sum_{i\neq j% }{\bm{v}}_{i}\cdot{\bm{v}}_{j}\delta(r-|{\bm{r}}_{i}-{\bm{r}}_{j}|)\right% \rangle}{\left\langle\sum_{i\neq j}\delta(r-|{\bm{r}}_{i}-{\bm{r}}_{j}|)\right% \rangle}.italic_C start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT ( italic_r ) = divide start_ARG 1 end_ARG start_ARG over¯ start_ARG italic_v end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG divide start_ARG ⟨ ∑ start_POSTSUBSCRIPT italic_i ≠ italic_j end_POSTSUBSCRIPT bold_italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ⋅ bold_italic_v start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT italic_δ ( italic_r - | bold_italic_r start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - bold_italic_r start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT | ) ⟩ end_ARG start_ARG ⟨ ∑ start_POSTSUBSCRIPT italic_i ≠ italic_j end_POSTSUBSCRIPT italic_δ ( italic_r - | bold_italic_r start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - bold_italic_r start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT | ) ⟩ end_ARG . (9)

where v¯2(i|𝒗i|2/N)superscript¯𝑣2subscript𝑖superscriptsubscript𝒗𝑖2𝑁\bar{v}^{2}\equiv\left(\sum_{i}|{\bm{v}}_{i}|^{2}/N\right)over¯ start_ARG italic_v end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ≡ ( ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | bold_italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_N ). Here, we use the velocity averaged over a short time interval ΔtΔ𝑡\Delta troman_Δ italic_t, for which we consider Δt=0.26σ/v0Δ𝑡0.26𝜎subscript𝑣0\Delta t=0.26\sigma/v_{0}roman_Δ italic_t = 0.26 italic_σ / italic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. The energy spectrum can be calculated via Fourier transformation. Here we consider the Fourier sine transform [68]

Esq(k)=kπdrrsinkrv¯2Csq(r).subscript𝐸sq𝑘𝑘𝜋differential-d𝑟𝑟𝑘𝑟superscript¯𝑣2subscript𝐶sq𝑟\displaystyle E_{\rm sq}(k)=\frac{k}{\pi}\int{\rm d}r\,r\sin{kr}\,\bar{v}^{2}C% _{\rm sq}(r).italic_E start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT ( italic_k ) = divide start_ARG italic_k end_ARG start_ARG italic_π end_ARG ∫ roman_d italic_r italic_r roman_sin italic_k italic_r over¯ start_ARG italic_v end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_C start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT ( italic_r ) . (10)

Fluid: Field-based approach.

We first extract the fluid velocity field 𝒗flsubscript𝒗fl{\bm{v}}_{\rm fl}bold_italic_v start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT from the simulation data by introducing a grid dividing the whole system into Ng3superscriptsubscript𝑁𝑔3N_{g}^{3}italic_N start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT cells. The velocities of all MPC particles, averaged over a short time interval ΔtΔ𝑡\Delta troman_Δ italic_t as for squirmers, are additionally averaged over each cell to obtain 𝒗𝒏subscript𝒗𝒏{\bm{v}}_{\bm{n}}bold_italic_v start_POSTSUBSCRIPT bold_italic_n end_POSTSUBSCRIPT with 𝒏=(nx,ny,nz)T𝒏superscriptsubscript𝑛𝑥subscript𝑛𝑦subscript𝑛𝑧𝑇{\bm{n}}=(n_{x},n_{y},n_{z})^{T}bold_italic_n = ( italic_n start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT , italic_n start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT , italic_n start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT for ni=0,,Ng1subscript𝑛𝑖0subscript𝑁𝑔1n_{i}=0,\ldots,N_{g}-1italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 , … , italic_N start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT - 1 and i{x,y,z}𝑖𝑥𝑦𝑧i\in\{x,y,z\}italic_i ∈ { italic_x , italic_y , italic_z }. Then the discrete Fourier transform

𝒗fl(𝒌)=1Ng3𝒏𝒗sq(𝒏)ei2πa𝒌𝒏/Ngsubscript𝒗fl𝒌1superscriptsubscript𝑁𝑔3subscript𝒏subscript𝒗sq𝒏superscript𝑒𝑖2𝜋𝑎𝒌𝒏subscript𝑁𝑔\displaystyle{\bm{v}}_{\rm fl}({\bm{k}})=\frac{1}{N_{g}^{3}}\sum_{\bm{n}}{\bm{% v}}_{\rm sq}({\bm{n}})e^{i2\pi a{\bm{k}}\cdot{\bm{n}}/N_{g}}bold_italic_v start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT ( bold_italic_k ) = divide start_ARG 1 end_ARG start_ARG italic_N start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG ∑ start_POSTSUBSCRIPT bold_italic_n end_POSTSUBSCRIPT bold_italic_v start_POSTSUBSCRIPT roman_sq end_POSTSUBSCRIPT ( bold_italic_n ) italic_e start_POSTSUPERSCRIPT italic_i 2 italic_π italic_a bold_italic_k ⋅ bold_italic_n / italic_N start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT end_POSTSUPERSCRIPT (11)

is performed. The energy spectrum is calculated straightforwardly via [69]

Efl(𝒌)=12|𝒗fl(𝒌)|2,subscript𝐸fl𝒌12superscriptsubscript𝒗fl𝒌2\displaystyle E_{\rm fl}({\bm{k}})=\frac{1}{2}|{\bm{v}}_{\rm fl}({\bm{k}})|^{2},italic_E start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT ( bold_italic_k ) = divide start_ARG 1 end_ARG start_ARG 2 end_ARG | bold_italic_v start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT ( bold_italic_k ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (12)

which is then averaged over all directions of 𝒌𝒌\bm{k}bold_italic_k to obtain Efl(k)subscript𝐸fl𝑘E_{\rm fl}(k)italic_E start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT ( italic_k ). Then, the spatial velocity correlation function is obtained via the Fourier transformation

Cv(𝒏)=𝒏Efl(𝒌)ei2πa𝒌𝒏/Ng,subscript𝐶𝑣𝒏subscript𝒏delimited-⟨⟩subscript𝐸fl𝒌superscript𝑒𝑖2𝜋𝑎𝒌𝒏subscript𝑁𝑔\displaystyle C_{v}({\bm{n}})=\sum_{\bm{n}}\left\langle E_{\rm fl}({\bm{k}})% \right\rangle e^{-i2\pi a{\bm{k}}\cdot{\bm{n}}/N_{g}},italic_C start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT ( bold_italic_n ) = ∑ start_POSTSUBSCRIPT bold_italic_n end_POSTSUBSCRIPT ⟨ italic_E start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT ( bold_italic_k ) ⟩ italic_e start_POSTSUPERSCRIPT - italic_i 2 italic_π italic_a bold_italic_k ⋅ bold_italic_n / italic_N start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT end_POSTSUPERSCRIPT , (13)

from which we calculate Cfl(r)subscript𝐶fl𝑟C_{\rm fl}(r)italic_C start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT ( italic_r ) by averaging over all directions of 𝒏𝒏\bm{n}bold_italic_n. To reduce noise effects at small length scales, velocity fields are averaged over boxes with side lengths σ𝜎\sigmaitalic_σ, 2σ2𝜎2\sigma2 italic_σ, or 3σ3𝜎3\sigma3 italic_σ.

The vorticity field is then defined as

𝝎(𝒏)𝒏×𝒗fl(𝒏).𝝎𝒏subscriptbold-∇𝒏subscript𝒗fl𝒏\displaystyle{\bm{\omega}}({\bm{n}})\equiv{\bm{\nabla}}_{\bm{n}}\times{\bm{v}}% _{\rm fl}({\bm{n}}).bold_italic_ω ( bold_italic_n ) ≡ bold_∇ start_POSTSUBSCRIPT bold_italic_n end_POSTSUBSCRIPT × bold_italic_v start_POSTSUBSCRIPT roman_fl end_POSTSUBSCRIPT ( bold_italic_n ) . (14)

Numerically, the vorticity field is calculated by the five-point stencil method from the velocity field. The vorticity spatial correlation Cωsubscript𝐶𝜔C_{\omega}italic_C start_POSTSUBSCRIPT italic_ω end_POSTSUBSCRIPT are also obtained from Eqs. (11)-(13). Moreover, a gliding time average is performed for the velocity with a time window ΔtΔ𝑡\Delta troman_Δ italic_t, corresponding to v0Δt0.26σsubscript𝑣0Δ𝑡0.26𝜎v_{0}\Delta t\approx 0.26\sigmaitalic_v start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT roman_Δ italic_t ≈ 0.26 italic_σ.

Cross correlation.

We again utilize Fourier transformation to calculate cross correlations between velocity and vorticity fields. Specifically, we first calculate the magnitudes of velocity and vorticity fields, which are then shifted by their average values, i.e., v~(𝒏)v(𝒏)𝒏v(𝒏)/Ng3~𝑣𝒏𝑣𝒏subscript𝒏𝑣𝒏superscriptsubscript𝑁𝑔3\tilde{v}(\bm{n})\equiv v(\bm{n})-\sum_{\bm{n}}v(\bm{n})/N_{g}^{3}over~ start_ARG italic_v end_ARG ( bold_italic_n ) ≡ italic_v ( bold_italic_n ) - ∑ start_POSTSUBSCRIPT bold_italic_n end_POSTSUBSCRIPT italic_v ( bold_italic_n ) / italic_N start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT and ω~(𝒏)ω(𝒏)𝒏ω(𝒏)/Ng3~𝜔𝒏𝜔𝒏subscript𝒏𝜔𝒏superscriptsubscript𝑁𝑔3\tilde{\omega}(\bm{n})\equiv\omega(\bm{n})-\sum_{\bm{n}}\omega(\bm{n})/N_{g}^{3}over~ start_ARG italic_ω end_ARG ( bold_italic_n ) ≡ italic_ω ( bold_italic_n ) - ∑ start_POSTSUBSCRIPT bold_italic_n end_POSTSUBSCRIPT italic_ω ( bold_italic_n ) / italic_N start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT. Then, from the Fourier transforms of the fields v~(𝒌)~𝑣𝒌\tilde{v}({\bm{k}})over~ start_ARG italic_v end_ARG ( bold_italic_k ) and ω~(𝒌)~𝜔𝒌\tilde{\omega}({\bm{k}})over~ start_ARG italic_ω end_ARG ( bold_italic_k ), the cross correlation is obtained via

Cvω(𝒏)=𝒏v~(𝒌)ω~(𝒌)ei2πa𝒌𝒏/Ng,subscript𝐶𝑣𝜔𝒏subscript𝒏delimited-⟨⟩~𝑣𝒌superscript~𝜔𝒌superscript𝑒𝑖2𝜋𝑎𝒌𝒏subscript𝑁𝑔\displaystyle C_{v\omega}({\bm{n}})=\sum_{\bm{n}}\left\langle\tilde{v}({\bm{k}% })\tilde{\omega}^{*}({\bm{k}})\right\rangle e^{-i2\pi a{\bm{k}}\cdot{\bm{n}}/N% _{g}},italic_C start_POSTSUBSCRIPT italic_v italic_ω end_POSTSUBSCRIPT ( bold_italic_n ) = ∑ start_POSTSUBSCRIPT bold_italic_n end_POSTSUBSCRIPT ⟨ over~ start_ARG italic_v end_ARG ( bold_italic_k ) over~ start_ARG italic_ω end_ARG start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ( bold_italic_k ) ⟩ italic_e start_POSTSUPERSCRIPT - italic_i 2 italic_π italic_a bold_italic_k ⋅ bold_italic_n / italic_N start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT end_POSTSUPERSCRIPT , (15)

where superscript * indicates the complex conjugate. Cvω(r)subscript𝐶𝑣𝜔𝑟C_{v\omega}(r)italic_C start_POSTSUBSCRIPT italic_v italic_ω end_POSTSUBSCRIPT ( italic_r ) is obtained by averaging over all 𝒏𝒏\bm{n}bold_italic_n directions.

References

  • Pedley and Kessler [1992] T. J. Pedley and J. O. Kessler, Hydrodynamic phenomena in suspensions of swimming microorganisms, Annu. Rev. Fluid Mech. 24, 313 (1992).
  • Be’er and Ariel [2019] A. Be’er and G. Ariel, A statistical physics view of swarming bacteria, Mov. Ecol. 7, 9 (2019).
  • Katz et al. [2011] Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin, Inferring the structure and dynamics of interactions in schooling fish, Proc. Natl. Acad. Sci. USA 108, 18720 (2011).
  • Gautrais et al. [2012] J. Gautrais, F. Ginelli, R. Fournier, S. Blanco, M. Soria, H. Chaté, and G. Theraulaz, Deciphering interactions in moving animal groups, PLOS Comput. Biol. 8, 1 (2012).
  • Ballerini et al. [2008] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Empirical investigation of starling flocks: a benchmark study in collective animal behaviour, Anim. Behav. 76, 201 (2008).
  • Bialek et al. [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Statistical mechanics for natural flocks of birds, Proc. Natl. Acad. Sci. USA 109, 4786 (2012).
  • Helbing et al. [2000] D. Helbing, I. Farkas, and T. Vicsek, Simulating dynamical features of escape panic, Nature , 487 (2000).
  • Xie et al. [2019] H. Xie, M. Sun, X. Fan, Z. Lin, W. Chen, L. Wang, L. Dong, and Q. He, Reconfigurable magnetic microrobot swarm: Multimode transformation, locomotion, and manipulation, Sci. Robot. 4, eaav8006 (2019).
  • Wang et al. [2021] Q. Wang, K. F. Chan, K. Schweizer, X. Du, D. Jin, S. C. H. Yu, B. J. Nelson, and L. Zhang, Ultrasound doppler-guided real-time navigation of a magnetic microswarm for active endovascular delivery, Sci. Adv. 7, eabe5914 (2021).
  • Chen and Bechinger [2022] C.-J. Chen and C. Bechinger, Collective response of microrobotic swarms to external threats, N. J. Phys. 24, 033001 (2022).
  • Wang et al. [2023] X. Wang, P.-C. Chen, K. Kroy, V. Holubec, and F. Cichos, Spontaneous vortex formation by microswimmers with retarded attractions, Nat. Commun. 14, 56 (2023).
  • Fily and Marchetti [2012] Y. Fily and M. C. Marchetti, Athermal phase separation of self-propelled particles with no alignment, Phys. Rev. Lett. 108, 235702 (2012).
  • Buttinoni et al. [2013] I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles, Phys. Rev. Lett. 110, 238301 (2013).
  • Wensink et al. [2012] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Meso-scale turbulence in living fluids, Proc. Natl. Acad. Sci. USA 109, 14308 (2012).
  • Qi et al. [2022] K. Qi, E. Westphal, G. Gompper, and R. G. Winkler, Emergence of active turbulence in microswimmer suspensions due to active hydrodynamic stress and volume exclusion, Commun. Phys. 5, 49 (2022).
  • Aranson [2022] I. S. Aranson, Bacterial active matter, Rep. Prog. Phys. 85, 076601 (2022).
  • Ziepke et al. [2022] A. Ziepke, I. Maryshev, I. S. Aranson, and E. Frey, Multi-scale organization in communicating active matter, Nat. Commun. 13, 6727 (2022).
  • Sawicki et al. [2023] J. Sawicki, R. Berner, S. A. M. Loos, M. Anvari, R. Bader, W. Barfuss, N. Botta, N. Brede, I. Franović, D. J. Gauthier, S. Goldt, A. Hajizadeh, P. Hövel, O. Karin, P. Lorenz-Spreen, C. Miehl, J. Mölter, S. Olmi, E. Schöll, A. Seif, P. A. Tass, G. Volpe, S. Yanchuk, and J. Kurths, Perspectives on adaptive dynamical systems, Chaos 33, 071501 (2023).
  • Negi et al. [2024] R. S. Negi, R. G. Winkler, and G. Gompper, Collective behavior of self-steering active particles with velocity alignment and visual perception, Phys. Rev. Res. 6, 013118 (2024).
  • Koch and Subramanian [2011] D. L. Koch and G. Subramanian, Collective hydrodynamics of swimming microorganisms: Living fluids, Annu. Rev. Fluid Mech. 43, 637 (2011).
  • Elgeti et al. [2015] J. Elgeti, R. G. Winkler, and G. Gompper, Physics of microswimmers—single particle motion and collective behavior: a review, Rep. Prog. Phys. 78, 056601 (2015).
  • Dombrowski et al. [2004] C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein, and J. O. Kessler, Self-concentration and large-scale coherence in bacterial dynamics, Phys. Rev. Lett. 93, 098103 (2004).
  • Drescher et al. [2011] K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly, and R. E. Goldstein, Fluid dynamics and noise in bacterial cell-cell and cell-surface scattering, Proc. Natl. Acad. Sci. USA 108, 10940 (2011).
  • Deneke et al. [2019] V. E. Deneke, A. Puliafito, D. Krueger, A. V. Narla, A. De Simone, L. Primo, M. Vergassola, S. De Renzis, and S. Di Talia, Self-organized nuclear positioning synchronizes the cell cycle in drosophila embryos, Cell 177, 925 (2019).
  • Hernández-López et al. [2023] C. Hernández-López, A. Puliafito, Y. Xu, Z. Lu, S. D. Talia, and M. Vergassola, Two-fluid dynamics and micron-thin boundary layers shape cytoplasmic flows in early Drosophila embryos, Proc. Natl. Acad. Sci. USA 120, e2302879120 (2023).
  • Ceron et al. [2023] S. Ceron, G. Gardi, K. Petersen, and M. Sitti, Programmable self-organization of heterogeneous microrobot collectives, Proc. Natl. Acad. Sci. USA 120, e2221913120 (2023).
  • Vicsek et al. [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75, 1226 (1995).
  • Chaté [2020] H. Chaté, Dry aligning dilute active matter, Annu. Rev. Condens. Matter Phys. 11, 189 (2020).
  • Reinken et al. [2018] H. Reinken, S. H. L. Klapp, M. Bär, and S. Heidenreich, Derivation of a hydrodynamic theory for mesoscale dynamics in microswimmer suspensions, Phys. Rev. E 97, 022613 (2018).
  • Liu et al. [2021a] Z. T. Liu, Y. Shi, Y. Zhao, H. Chaté, X. qing Shi, and T. H. Zhang, Activity waves and freestanding vortices in populations of subcritical quincke rollers, Proc. Natl. Acad. Sci. USA 118, e2104724118 (2021a).
  • Aditi Simha and Ramaswamy [2002] R. Aditi Simha and S. Ramaswamy, Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles, Phys. Rev. Lett. 89, 058101 (2002).
  • Llopis and Pagonabarraga [2010] I. Llopis and I. Pagonabarraga, Hydrodynamic interactions in squirmer motion: Swimming with a neighbour and close to a wall, J. Non-Newtonian Fluid Mech. 165, 946 (2010).
  • Clopés et al. [2020] J. Clopés, G. Gompper, and R. G. Winkler, Hydrodynamic interactions in squirmer dumbbells: active stress-induced alignment and locomotion, Soft Matter 16, 10676 (2020).
  • Lighthill [1952] M. J. Lighthill, On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers, Comm. Pure Appl. Math. 5, 109 (1952).
  • Blake [1971] J. R. Blake, A spherical envelope approach to ciliary propulsion, J. Fluid Mech. 46, 199 (1971).
  • Ishikawa et al. [2006] T. Ishikawa, M. P. Simmonds, and T. J. Pedley, Hydrodynamic interaction of two swimming model micro-organisms, J. Fluid Mech. 568, 119 (2006).
  • Götze and Gompper [2010] I. O. Götze and G. Gompper, Mesoscale simulations of hydrodynamic squirmer interactions, Phys. Rev. E 82, 041921 (2010).
  • Theers et al. [2016a] M. Theers, E. Westphal, G. Gompper, and R. G. Winkler, Modeling a spheroidal microswimmer and cooperative swimming in a narrow slit, Soft Matter 12, 7372 (2016a).
  • Pak and Lauga [2014] O. S. Pak and E. Lauga, Generalized squirming motion of a sphere, J. Eng. Math. 88, 1 (2014).
  • Goh et al. [2023] S. Goh, R. G. Winkler, and G. Gompper, Hydrodynamic pursuit by cognitive self-steering microswimmers, Commun. Phys. 6, 310 (2023).
  • Malevanets and Kapral [1999] A. Malevanets and R. Kapral, Mesoscopic model for solvent dynamics, J. Chem. Phys. 110, 8605 (1999).
  • Kapral [2008] R. Kapral, Multiparticle collision dynamics: Simulations of complex systems on mesoscale, Adv. Chem. Phys. 140, 89 (2008).
  • Gompper et al. [2009] G. Gompper, T. Ihle, D. M. Kroll, and R. G. Winkler, Multi-particle collision dynamics: A particle-based mesoscale simulation approach to the hydrodynamics of complex fluids, Adv. Polym. Sci. 221, 1 (2009).
  • Goh et al. [2022] S. Goh, R. G. Winkler, and G. Gompper, Noisy pursuit and pattern formation of self-steering active particles, New J. Phys. 24, 093039 (2022).
  • Chepizhko et al. [2021] O. Chepizhko, D. Saintillan, and F. Peruani, Revisiting the emergence of order in active matter, Soft Matter 17, 3113 (2021).
  • Bera et al. [2022] A. Bera, S. Sahoo, S. Thakur, and S. K. Das, Active particles in explicit solvent: Dynamics of clustering for alignment interaction, Phys. Rev. E 105, 014606 (2022).
  • Noguchi and Gompper [2008] H. Noguchi and G. Gompper, Transport coefficients of off-lattice mesoscale-hydrodynamics simulation techniques, Phys. Rev. E 78, 016706 (2008).
  • Westphal et al. [2024] E. Westphal, S. Goh, R. G. Winkler, and G. Gompper, Htmpc: A heavily templated C++ library for large scale particle-based mesoscale hydrodynamics simulations using multiparticle collision dynamics, arXiv preprint arXiv:2406.15236 10.48550/arXiv.2406.15236 (2024).
  • Rycroft [2009] C. H. Rycroft, VORO++: A three-dimensional Voronoi cell library in C++, Chaos 19, 041111 (2009).
  • Gascó et al. [2023] A. Gascó, I. Pagonabarraga, and A. Scagliarini, Three-dimensional active turbulence in microswimmer suspensions: simulations and modelling, arXiv preprint arXiv:2304.03662 10.48550/arXiv.2304.03662 (2023).
  • Stenhammar et al. [2017] J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo, and A. Morozov, Role of correlations in the collective behavior of microswimmer suspensions, Phys. Rev. Lett. 119, 028005 (2017).
  • Bárdfalvy et al. [2019] D. Bárdfalvy, H. Nordanger, C. Nardini, A. Morozov, and J. Stenhammar, Particle-resolved lattice Boltzmann simulations of 3-dimensional active turbulence, Soft Matter 15, 7747 (2019).
  • Huang et al. [2012] C.-C. Huang, G. Gompper, and R. G. Winkler, Hydrodynamic correlations in multiparticle collision dynamics fluids, Phys. Rev. E 86, 056711 (2012).
  • Wensink and Löwen [2012] H. H. Wensink and H. Löwen, Emergent states in dense systems of active rods: from swarming to turbulence, J. Phys.: Condens. Matter 24, 460130 (2012).
  • Keta et al. [2024] Y.-E. Keta, J. U. Klamser, R. L. Jack, and L. Berthier, Emerging mesoscale flows and chaotic advection in dense active matter, Phys. Rev. Lett. 132, 218301 (2024).
  • Theers et al. [2018] M. Theers, E. Westphal, K. Qi, R. G. Winkler, and G. Gompper, Clustering of microswimmers: interplay of shape and hydrodynamics, Soft Matter 14, 8590 (2018).
  • Zantop and Stark [2022] A. W. Zantop and H. Stark, Emergent collective dynamics of pusher and puller squirmer rods: swarming, clustering, and turbulence, Soft Matter 18, 6179 (2022).
  • Worlitzer et al. [2021a] V. M. Worlitzer, G. Ariel, A. Be’er, H. Stark, M. Bär, and S. Heidenreich, Motility-induced clustering and meso-scale turbulence in active polar fluids, New J. Phys. 23, 033012 (2021a).
  • Worlitzer et al. [2021b] V. M. Worlitzer, G. Ariel, A. Be’er, H. Stark, M. Bär, and S. Heidenreich, Turbulence-induced clustering in compressible active fluids, Soft Matter 17, 10447 (2021b).
  • Huang et al. [2010] C.-C. Huang, A. Chatterji, G. Sutmann, G. Gompper, and R. G. Winkler, Cell-level canonical sampling by velocity scaling for multiparticle collision dynamics simulations, J. Comput. Phys. 229, 168 (2010).
  • Omelyan [1998] I. P. Omelyan, Algorithm for numerical integration of the rigid-body equations of motion, Phys. Rev. E 58, 1169 (1998).
  • Lamura et al. [2001] A. Lamura, G. Gompper, T. Ihle, and D. M. Kroll, Multiparticle collision dynamics: Flow around a circular and a square cylinder,  Europhys. Lett. 56, 319 (2001).
  • Padding et al. [2005] J. T. Padding, A. Wysocki, H. Löwen, and A. A. Louis, Stick boundary conditions and rotational velocity auto-correlation functions for colloidal particles in a coarse-grained representation of the solvent, J. Phys.: Condens. Matter 17, S3393 (2005).
  • Theers and Winkler [2014] M. Theers and R. G. Winkler, Effects of thermal fluctuations and fluid compressibility on hydrodynamic synchronization of microrotors at finite oscillatory Reynolds number: A multiparticle collision dynamics simulation study, Soft Matter 10, 5894 (2014).
  • Theers and Winkler [2015] M. Theers and R. G. Winkler, Bulk viscosity of multiparticle collision dynamics fluids, Phys. Rev. E 91, 033309 (2015).
  • Shaebani et al. [2020] M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper, and H. Rieger, Computational models for active matter, Nat. Rev. Phys. 2, 181 (2020).
  • Theers et al. [2016b] M. Theers, E. Westphal, G. Gompper, and R. G. Winkler, From local to hydrodynamic friction in Brownian motion: A multiparticle collision dynamics simulation study, Phys. Rev. E 93, 032604 (2016b).
  • Batchelor [1959] G. K. Batchelor, The theory of homogeneous turbulence (University Press, Cambridge, 1959).
  • Liu et al. [2021b] Z. Liu, W. Zeng, X. Ma, and X. Cheng, Density fluctuations and energy spectra of 3D bacterial suspensions, Soft Matter 17, 10806 (2021b).

Acknowledgements

Funding: The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). Author contributions: G.G. and R.G.W. designed the research. E.W. wrote the simulation code. S.G. performed the simulations and analyzed the data. S.G., R.G.W., and G.G. discussed the results and wrote the manuscript together. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are available from the corresponding author upon reasonable request.