2.1. Squirmer model

We consider spherical swimmers that propel themselves by generating tangential velocities along their surface, known as a squirmer. For each swimmer, the surface velocity is given by (Lighthill Reference Lighthill1952; Blake Reference Blake1971; Ishikawa et al. Reference Ishikawa, Simmonds and Pedley2006)

(2.1) \begin{equation} \boldsymbol{u}_{\!s} = B_1\left (\sin \theta { + \beta }\sin \theta \cos \theta \right )\boldsymbol{e}_\theta , \end{equation}

where $\theta =\cos ^{-1}(\boldsymbol{p} \boldsymbol{\cdot }\boldsymbol{r} / r)$ represents the polar angle between the position vector $\boldsymbol{r}$ , measured from the particle’s centre, and the swimming direction $\boldsymbol{p}$ , while ${\boldsymbol{e}}_\theta$ is the unit tangent vector along the $\theta$ direction. Throughout this work, we fix $B_1$ which determines the swimming speed of a solitary squirmer in the Stokes regime, where $V_s=2B_1/3$ , which is independent of $\beta$ . The squirmer parameter $\beta$ regulates the strength of the force dipole generated by the swimmer and determines whether the squirmer behaves as a puller ( $\beta \gt 0$ ) or a pusher ( $\beta \lt 0$ ). Furthermore, in the present study, we neglect bottom-heaviness of the squirmer and assume that the centre of mass coincides with the centre of the sphere.

2.2. Basic equations

The transport of swimmers in a fluid is simulated by direct numerical solution of the Navier–Stokes equations using the fictitious domain method (Yu & Shao Reference Yu and Shao2007). Here, we briefly describe the basic equations; details of the method can be found in Lin & Gao (Reference Lin and Gao2019). Let $\varOmega$ and $P$ represent the entire domain and the domain occupied by the swimmers, respectively ( $P \subset \varOmega$ ). Here $P$ is bounded by the swimmer surface $S$ .

The flow is governed by the Navier–Stokes equations,

(2.2) \begin{equation} \rho _f \frac {\mathrm{d} \boldsymbol{u}}{\mathrm{d} t}= \boldsymbol{\nabla } \boldsymbol{\cdot }\left (-p \boldsymbol{I}+\boldsymbol{\tau }\right ) + \boldsymbol{f}\quad {\textrm {in}} \quad {\varOmega }. \end{equation}

(2.3) \begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}=\boldsymbol{0} \quad {\textrm {in}} \quad {\varOmega }. \end{equation}

The velocity boundary condition on swimmer surface is given by

(2.4) \begin{equation} \boldsymbol{u}=(\boldsymbol{U}+\boldsymbol{\omega } \times \boldsymbol{r})+\boldsymbol{u}_{\!s}\quad {\textrm {in}} \quad {S}. \end{equation}

Here, $\boldsymbol u$ , $p$ and $\rho _f$ represent the fluid velocity, pressure and density, respectively. The stress tensor $\boldsymbol{\tau }=2 \mu \boldsymbol{\dot \gamma }$ , where $\mu$ is the viscosity and $\boldsymbol{\dot \gamma }$ is the strain rate tensor. The body force $\boldsymbol{f}$ accounts for the fluid–particle interaction defined in the particle domain.

Next, we describe the motion of spherical swimmers with density $\rho _s$ and radius $a$ . The swimmer velocity $\boldsymbol U$ and angular velocity $\boldsymbol \omega$ obey Newton’s equations of motion,

(2.5) \begin{equation} \rho _s V_{\!p} \frac {\mathrm{d} \boldsymbol{U}}{\mathrm{d} t}=\int _S \boldsymbol{n} \boldsymbol{\cdot }\left (-p \boldsymbol{I}+\boldsymbol{\tau }\right ) \,\mathrm{d} s+(\rho _s-\rho _f)V_{\!p}\boldsymbol g, \end{equation}

(2.6) \begin{equation} \frac {\mathrm{d}\left (\boldsymbol{J} \boldsymbol{\cdot }{\boldsymbol{\omega }}\right )}{\mathrm{d} t}=\int _S \boldsymbol{r} \times \left [\boldsymbol{n} \boldsymbol{\cdot }\left (-p \boldsymbol{I}+\boldsymbol{\tau }\right )\right ] \,\mathrm{d} s. \end{equation}

Here, $V_{\!p}$ and $\boldsymbol J$ denote the swimmer volume and the moment of inertia tensor, respectively. The vector $\boldsymbol{r}$ is measured relative to the particle centre of mass, and $\boldsymbol{g}$ denotes the gravitational acceleration.

The position of each swimmer, $\boldsymbol{x}=(x,y,z)^T$ , is determined from the following kinematic equation:

(2.7) \begin{equation} \frac {{{\textrm {d}}{\boldsymbol{x}}}}{{{\textrm {d}}t}} = {\boldsymbol{U}}. \end{equation}

To accurately capture the orientation of swimmers, we employ a body-fixed coordinate system ( $x'$ , $y'$ , $z'$ ), which is aligned with the swimmer’s centre, with the $z'$ -axis pointing in the direction of swimming, denoted as $\boldsymbol p$ . The swimmer orientation is determined by

(2.8) \begin{equation} \frac {\textrm {d} \boldsymbol q}{\textrm {d}t} = \frac {1}{2}\left (\begin{array}{cccc} q_4 & -q_3 & q_2 & q_1 \\[5pt] q_3 & q_4 & -q_1 & q_2 \\[5pt] -q_2 & q_1 & q_4 & q_3 \\[5pt] -q_1 & -q_2 & -q_3 & q_4 \end{array}\right )\left (\begin{array}{c} \omega '_x \\[5pt] \omega '_y \\[5pt] \omega '_z \\[5pt] 0 \end{array}\right ) \!. \end{equation}

In (2.8), $q_1$ , $q_2$ , $q_3$ and $q_4$ are the components of the unit quaternion $\boldsymbol q$ , which can be defined in terms of the Euler angles (Yu & Shao Reference Yu and Shao2007). The components of the angular velocity, $\omega '_x$ , $\omega '_y$ and $\omega '_z$ , are measured in the body-fixed frame. With the quaternion $\boldsymbol q$ , the rotation matrix from the space-fixed frame to the body-fixed frame is

(2.9) \begin{equation} \boldsymbol{R}=2\left (\begin{array}{ccc} q_1^2+q_4^2-\frac {1}{2} & q_1 q_2+q_3 q_4 & q_1 q_3-q_2 q_4 \\[5pt] q_1 q_2-q_3 q_4 & q_2^2+q_4^2-\frac {1}{2} & q_2 q_3+q_1 q_4 \\[5pt] q_1 q_3+q_2 q_4 & q_2 q_3-q_1 q_4 & q_3^2+q_4^2-\frac {1}{2} \end{array}\right ) \!. \end{equation}

Since the swimmer is rotationally symmetric, its orientation $\boldsymbol{p}$ is fully determined by the direction of its symmetry axis, the $z'$ axis, in the space-fixed frame. orientation vector is given by $\boldsymbol p = \boldsymbol{R} \boldsymbol{\cdot }(0,0,1)^T$ .

Figure 1.

Schematic representation of the model system. (a) Simulation set-up of the suspension. (b) Sketch of a squirmer, where $\gamma$ is the angle between the orientation vector $\boldsymbol p$ and the positive $x$ -axis.

The computational code has been extensively validated for the dynamics of squirmers at low Reynolds numbers (Xia et al. Reference Xia, Yu, Hu, Zhu and Lin2025b , Reference Xia, Yu, Zhang, Lin and Ouyangc ), as well as for the sedimentation of a single sphere (Yu & Shao Reference Yu and Shao2007; Xia et al. Reference Xia, Yu, Pan, Lin and Guo2022). Validations covering squirmer dynamics in the inertial regime and the hindered settling of passive spheres are presented in § 3.

2.3. Simulation set-up

As delineated in figure 1, we consider a suspension of $N$ spherical swimmers with radius $a$ in a periodic cubic domain of side length $L$ , yielding a volume fraction of swimmers $\phi =4\pi Na^3/(3L^3)$ . For the simulations presented in this work, the computational domain of size $L = 32a$ is discretised on a fixed Cartesian grid with $N_m = 512$ mesh points in each direction, corresponding to 32 grid points across the swimmer diameter (i.e. $2a/\varDelta = 32$ ), which is sufficient to accurately resolve the hydrodynamic interactions between squirmers. The volume fraction of swimmers ( $\phi$ ) ranges from 0.025 to 0.1, corresponding to a total number of swimmers ranging from 196 to 782. Initially, the swimmers are randomly distributed in space with random orientations. The effect of the domain size $L$ has been examined numerically in Appendix B. In the simulations, zero-volume-flux condition ( $\phi \langle \boldsymbol U \rangle +(1-\phi ) \langle \boldsymbol u \rangle =\boldsymbol 0$ , with $\left \langle \boldsymbol{\cdot }\right \rangle$ denoting an instantaneous average over all particles or fluid nodes) is imposed via a dynamically adjusted pressure gradient. All computations were normalised by using the swimmer radius $a$ for length, $V_s$ for velocity, $a/V_s$ for time, $\mu V_s/a$ for stress and $\mu V_s/a^2$ for body force. We treat the steric interaction between swimmers by adding an inter-particle repulsive force (Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008) in order to avoid the overlap of swimmers. The time step in all simulations is set to $\Delta t = 10^{-3} a/V_s$ . Each simulation is run for a total duration of $2500a/V_s$ to ensure that the measured quantities were fully converged.

Fluid inertia is characterised by the swimming Reynolds number

(2.10) \begin{equation} \textit{Re}_S=\frac {\rho _fV_sa}{\mu }, \end{equation}

which represents the ratio between the inertial force of the fluid and the viscous force due to swimming. Owing to their small size, most microorganisms swim in the Stokes regime, where fluid inertia is negligible. However, in nature, there exists a wide variety of swimming organisms. For instance, the millimetre-sized copepods and Pleurobrachia commonly achieve $\textit{Re}_S$ up to $O(10)$ (Kiørboe et al. Reference Kiørboe, Jiang and Colin2010). Motivated by such examples, we investigate the hydrodynamic diffusion of swimmers with $0.01 \leqslant \textit{Re}_S \leqslant 10$ , covering regimes of both negligible and finite fluid inertia.

The numerical cases are further classified into neutrally buoyant swimmers and settling swimmers, depending on whether a density mismatch exists between the swimmer and the surrounding fluid. To isolate the effects of fluid inertia and sedimentation, Brownian motion and bottom-heaviness are neglected in the present study. The density of microswimmers is usually close to that of the background fluid. For neutrally buoyant swimmers, the density ratio $\rho _s/\rho _f$ is equal to unity; therefore, they do not settle under gravity. For settling swimmers, $\rho _s/\rho _f$ is chosen as 1.02 to reflect the small density mismatch commonly found in natural aquatic microorganisms. The gravitational direction is taken in the negative $x$ direction, i.e. $\boldsymbol{g}=-\boldsymbol{e}_{\!x}g$ , with $g=\|\boldsymbol g\|$ . The ratio of the settling velocity of a sphere to the swimming velocity of a solitary squirmer in the Stokes limit will be an important parameter to quantify the settling effect. For a spherical particle, the Stokes settling velocity is given by

(2.11) \begin{equation} V_g=\frac {2\left (\rho _s-\rho _f\right )ga^2}{9\mu }. \end{equation}

It should be noted that the actual fall velocity of dilute spherical particles in a periodic domain is smaller than $V_g$ due to the fluid inertial and finite-domain effect. We refer to $\alpha =V_g/V_s$ as the settling coefficient, and $\alpha =0$ corresponds to the neutrally buoyant cases. Finally, the settling Reynolds number is defined as

(2.12) \begin{equation} \textit{Re}_G=\frac {\rho _f V_g a}{\mu }, \end{equation}

which can also be expressed as $\textit{Re}_G=\alpha \textit{Re}_S$ . We restrict $\alpha \leqslant 1$ (i.e. $V_s \geqslant V_g$ ) to ensure that self-propulsion dominates over gravitational settling. When $\alpha \gg 1$ , gravitational settling overwhelms the swimming; such conditions are beyond the scope of the present study. The squirmer parameter was set to $\beta = \pm 5$ , motivated by previous studies reporting qualitative differences between strong pullers and pushers. The definitions and values of non-dimensional parameters are summarised in table 1.

Table 1.

List of the non-dimensional parameters.

Table 1. Long description

A table with six rows and four columns. The columns are labeled Parameter name, Symbol, Definition, and Values. The rows list the following parameters: Swimming Reynolds number with symbol Re_S, Volume fraction with symbol phi, Squirmer parameter with symbol beta, Density ratio with symbol rho_S/rho_f, and Settling coefficient with symbol alpha. Each row provides the definition and corresponding values for each parameter.

Figure 2.

Swimming speed of a solitary squirmer, $U_{\textit{sol}}$ , as a function of $\textit{Re}_S$ . Here, $U_{\textit{sol}}$ is normalised by $V_s$ . The symbols represent our simulation results and lines show the result of Chisholm et al. (Reference Chisholm, Legendre, Lauga and Khair2016).

3.1. Single-swimmer speed

As reported by Wang & Ardekani (Reference Wang and Ardekani2012) and Chisholm et al. (Reference Chisholm, Legendre, Lauga and Khair2016), below a critical Reynolds number, a single squirmer exhibits a steady swimming speed $U_{\textit{sol}}$ , which depends on both $\beta$ and $\textit{Re}_S$ . The calculated steady swimming speed $U_{\textit{sol}}$ versus $\textit{Re}_S$ with $\beta =\pm 5$ for $\textit{Re}_S \leqslant 10$ is shown in figure 2. Pushers exhibit a monotonic increase in swimming speed with increasing $\textit{Re}_S$ , exceeding twice the Stokes-regime value at $\textit{Re}_S = 10$ . In contrast, pullers show a decrease in speed for $\textit{Re}_S \lesssim O(1)$ . As $\textit{Re}_S$ increases beyond $O(1)$ , the speed increases. Our results are in good agreement with the numerical results of Chisholm et al. (Reference Chisholm, Legendre, Lauga and Khair2016) using the spectral element method.

3.2. Trajectories of a swimmer pair

Two-body hydrodynamic interactions play a dominant role in transport of swimmers in semi-dilute suspensions. Figure 3 compares the trajectories of two pullers ( $\beta =5$ ) obtained from our simulations with the results reported by Li et al. (Reference Li, Ostace and Ardekani2016). The two swimmers are initially located in the same plane ( $z = 0$ ). They initially swim towards each other, and their swimming orientations are initially parallel. They are initially positioned with a centre-to-centre distance of $10a$ along the $x$ -axis and $a$ along the $y$ -axis. As shown in figure 3(a,b), at $\textit{Re}_S = 0.1$ and 1, the pullers rotate away from each other during the interaction and eventually escape with positive scattering angles. The scattering angle at $\textit{Re}_S = 1$ is larger than that at $\textit{Re}_S = 0.1$ . In figure 3(c), the pullers become hydrodynamically entrapped near each other, rotating along circular trajectories at $\textit{Re}_S = 10$ . Our simulation results are in agreement with the numerical values reported by Li et al. (Reference Li, Ostace and Ardekani2016). The validation is limited to pullers because, at Reynolds numbers of order unity, pushers develop intrinsically three-dimensional trajectories after close contact due to the amplification of out-of-plane perturbations. Consequently, the reference results do not provide a suitable benchmark for quantitative comparison.

Figure 3.

Trajectories of two pullers ( $\beta =5$ ) swimming in opposite directions at different Reynolds numbers: (a) $\textit{Re}_S=0.1$ , (b) $\textit{Re}_S=1$ and (c) $\textit{Re}_S=10$ .

3.3. Hindered settling velocity of passive spheres

The hindered settling of spherical particles is a classical and extensively studied problem. We compare our simulations with previous models and simulations for the average settling speed of randomly distributed non-Brownian spheres in the low-Reynolds-number regime. The Richardson–Zaki correlation $U_{sed}/{V_g}=(1-\phi )^{n}$ (Richardson & Zaki Reference Richardson and Zaki1997) summarises experimental data, and its validity has been confirmed in numerical simulations (Cunha et al. Reference Cunha, Abade, Sousa and Hinch2002) and experiments (Nicolai et al. Reference Nicolai, Herzhaft, Hinch, Oger and Guazzelli1995); the Hayakawa–Ichiki model $U_{sed}/{V_g}={(1-\phi )^3} /{ (1+2 \phi +1.429 \phi (1-\phi )^3 )}$ (Hayakawa & Ichiki Reference Hayakawa and Ichiki1995) has been favourably compared against simulations for a wide range of volume fractions (Padding & Louis Reference Padding and Louis2004).

In figure 4, numerical results from current simulations for passive spheres at low $\textit{Re}_G$ in a triply periodic cubic box of size $L=32a$ are plotted as a function of the volume fraction. The average settling speed decreases as the particle volume fraction increases due to the hindered effect. Our simulation results lie between the predictions of the Richardson–Zaki correlation and the Hayakawa–Ichiki model. The simulated speed for $\phi =0.025$ is slightly smaller than the values predicted by the models. This is expected because of the use of triply periodic boundary conditions in a finite domain (Phillips et al. Reference Phillips, Brady and Bossis1988). Our results also agree well with the simulation results of Li & Botto (Reference Li and Botto2024) using the Stokesian dynamics method with domain size of $L=56a$ .

Figure 4.

Average settling speed of randomly distributed spheres versus volume fraction.

4.1. Effects of swimming Reynolds number on diffusion

The diffusion coefficient is a measure of the increasing displacements of particles, which is calculated based on the mean-square displacement (MSD) of particles. In the case of neutrally buoyant swimmers, there is no preferred direction of motion, so particle motion is expected to be isotropic. Thus, we calculate the total MSD by

(4.1) \begin{equation} \big \langle \Delta r^2(t)\big \rangle \equiv \big \langle [\boldsymbol{x}(\tau +t)-\boldsymbol{x}(\tau )]^2\big \rangle , \end{equation}

where $\boldsymbol{x}(t)=\boldsymbol{x}(t_0)+\int ^t_{t_0} \boldsymbol{U}(t^\prime )\, \mathrm{d} t^\prime$ is the position of a particle at time $t$ and $\left \langle \boldsymbol{\cdot }\right \rangle$ is an average over all particles and reference times. If the MSD grows linearly with time, then the spreading is diffusive. The time-dependent diffusivity is defined as

(4.2) \begin{equation} D(t)= \frac {\left \langle \Delta r^2(t)\right \rangle }{6t}. \end{equation}

Then it allows us to evaluate the long-time diffusion coefficient as $\displaystyle D_t=\lim _{t \rightarrow \infty }D(t)$ .

Before discussing the diffusion coefficient, the MSDs of pullers and pushers at different swimming Reynolds numbers and $\phi =0.05$ are shown as a log–log plot in figure 5. We see that MSDs grow as ${\sim} t^2$ for short times and as ${\sim} t$ for long times in all cases, which confirms a ballistic-to-diffusive transition. Interestingly, the MSDs of pullers decrease with increasing $\textit{Re}_S$ , whereas for pushers the MSDs increase as $\textit{Re}_S$ is increased. For the puller suspension at $\textit{Re}_S=10$ , the MSD oscillates strongly in the intermediate-time regime, resulting in much smaller spreading. This oscillation in the MSD arises from pairwise interactions between pullers (see figure 3 c). The pullers undergo continuous rotation during near interactions, resulting in a small net displacement over long times. The modified diffusion process is also evident from trajectories of sample swimmers. As illustrated by the insets in figure 5, increasing $\textit{Re}_S$ renders puller trajectories more circular, while pusher trajectories become more extended and persistent. We note that the damped MSD oscillations share similarities with those observed for swimmers subject to an external torque (Sandoval Reference Sandoval2013) or shear flow (Ishikawa Reference Ishikawa2012) and bubbles rising in quiescent water (Huang et al. Reference Huang, Hessenkemper, Tan, Ni, Sommer, Bragg and Ma2025). However, the underlying physical mechanisms differ substantially. In those studies, the oscillations originate from forced rotations or wake instabilities, respectively, whereas in the present case they arise from pairwise interactions.

Figure 5.

Mean-square displacement of (a) pullers and (b) pushers with $\phi =0.05$ at various Reynolds numbers. Here $\tau _s=a/V_s$ is the swimming time scale. The insets show a few sample trajectories at $\textit{Re}_S=0.01$ and 10.0, respectively.

Figure 6.

Velocity autocorrelation $C_U(t)$ of (a) pullers ( $\beta =5$ ) and (b) pushers ( $\beta =-5$ ) at $\phi =0.05$ over a range of swimming Reynolds numbers.

To gain further insight into the origin of the contrasting behaviours observed for pullers and pushers, we examine the velocity autocorrelation function, which characterises the temporal persistence of swimmer motion. The normalised velocity autocorrelation function is defined as

(4.3) \begin{equation} C_U(t)=\frac {{{{\left \langle {{\boldsymbol{U}}\left ( \tau \right ) \boldsymbol{\cdot }{\boldsymbol{U}}\left ( {\tau + t} \right )} \right \rangle } }}}{{{{\left \langle {{{\boldsymbol{U}}^2}\left ( \tau \right )} \right \rangle } }}}. \end{equation}

Figure 6 presents $C_U(t)$ of swimmers under different $\textit{Re}_S$ . When fluid inertia becomes non-negligible, the fore–aft antisymmetric vorticity becomes unstable, leading to a qualitative modification of hydrodynamic interactions. Consequently, near-field encounters induce stronger angular deflections, thereby accelerating the decay of velocity autocorrelation. This interpretation is consistent with previous studies on pairwise interactions (Li et al. Reference Li, Ostace and Ardekani2016), which showed that pullers exhibit increasingly larger scattering angles with increasing Reynolds number. As a result, $C_U(t)$ exhibits a markedly faster decay with increasing $\textit{Re}_S$ as shown in figure 6(a). Notably, at $\textit{Re}_S=10$ , $C_U(t)$ displays a pronounced damped oscillation around zero, which gives rise to an oscillating MSD in the intermediate-time regime in figure 5(a). Unlike pullers, figure 6(b) shows that $C_U(t)$ for pushers decays more slowly as $\textit{Re}_S$ increases. This trend can be understood from the inertia-modified wake dynamics of pushers. As the Reynolds number increases, the pusher generates a more coherent and stable wake structure, characterised by increased wake velocity and a reduced wake width (Li et al. Reference Li, Ostace and Ardekani2016). Unlike pullers, whose wake tends to destabilise under inertia, pushers maintain a long-lived vorticity wake and do not exhibit immediate wake separation or a standing eddy even at moderate $\textit{Re}_S$ (Chisholm et al. Reference Chisholm, Legendre, Lauga and Khair2016). This persistent wake stabilises pusher’s orientation and thereby slowing the decay of velocity correlations.

Figure 7.

(a) Diffusion coefficient, (b) mean-square swimmer speed and (c) decorrelation time as functions of the swimming Reynolds number for different values of $\beta$ at $\phi = 0.05$ . All quantities are normalised by their corresponding values when $\textit{Re}_S = 0.01$ , i.e. $D_{t,0}$ .

Figure 7. Long description

Panel A: A scatter plot shows the diffusion coefficient normalized by its value at beta equals zero (D_t/D_t,0) as a function of the swimming Reynolds number (Re_S). The plot includes data points for beta values of 5, 3, −3, and −5, represented by different symbols. Panel B: A scatter plot displays the mean-square swimmer speed normalized by its value at beta equals zero (U_rms^2/U_rms,0) against the swimming Reynolds number (Re_S). Data points for beta values of 5, 3, −3, and −5 are shown with distinct symbols. Panel C: A scatter plot illustrates the decorrelation time normalized by its value at beta equals zero (t_p/t_p,0) as a function of the swimming Reynolds number (Re_S). The plot includes data points for beta values of 5, 3, −3, and −5, represented by different symbols.

We extract the diffusion coefficient $D_t$ at $\phi =0.05$ from the long-time MSDs shown in figure 7(a). Clearly, the results confirm that $D_t$ decreases with $\textit{Re}_S$ for pullers, whereas it increases with $\textit{Re}_S$ for pushers. The diffusion coefficient for pullers at $\textit{Re}_S = 10$ is more than an order of magnitude smaller than its value in the low- $Re$ limit, whereas for pushers it increases several times over the same range. The figure further illustrates the effect of $\beta$ on the diffusion coefficient, revealing that the sensitivity of $\textit{Re}_S$ on $D_t$ becomes stronger as $\left |\beta \right |$ increases. The diffusion coefficient $D_t$ can be expressed in terms of the velocity autocorrelation function through the Green–Kubo relation

(4.4) \begin{equation} {D}_t=\frac {1}{3}\int _0^{\infty }\langle \boldsymbol{U}(\tau ) \boldsymbol{\cdot }\boldsymbol{U}(\tau +t)\rangle \,\mathrm{d} t. \end{equation}

Then we have ${D}_t=U^2_{rms}t_p$ , where $U^2_{rms}={\left \langle {{{\boldsymbol{U}}^2} ( \tau )} \right \rangle }/3$ is the mean-square speed of swimmers and $t_p=\int _0^{\infty } C_U(t)\, \mathrm{d} t$ is a time scale of the velocity decorrelation. Figure 7(b,c) shows $U^2_{rms}$ and $t_p$ of swimmers at $\phi =0.05$ , normalised by their corresponding values when $\textit{Re}_S=0.01$ . For pullers, both $U^2_{rms}$ and $t_p$ decrease with increasing $\textit{Re}_S$ , whereas for pushers, both $U^2_{rms}$ and $t_p$ increase with $\textit{Re}_S$ . Consequently, fluid inertia reduces the swimming speed and accelerates the velocity decorrelation of pullers, whereas it has the opposite effect on pushers, enhancing their swimming speed and prolonging the velocity correlation.

Figure 8.

Correlation between $D_t$ and $\phi$ at different $\textit{Re}_S$ and $\beta$ . The dashed line $D_t/(aV_s)= 0.24\phi ^{-1}$ (Ishikawa & Pedley Reference Ishikawa and Pedley2007a ) for $\beta =5$ is drawn to guide the eye.

4.2. Scaling of diffusion coefficient

Dimensionally, $D_t$ has units of ${\textrm {m}^2}/{\rm s}^{-1}$ , and so can be expressed as $U^2_cT_c$ , with $U_c$ and $T_c$ being the velocity and time scales. In semi-dilute suspensions of squirmers without inertia, it has been reported that $D_t$ scales inversely with the volume fraction, i.e. $D_t \sim \phi ^{-1}$ (Ishikawa & Pedley Reference Ishikawa and Pedley2007a ). We first revisit the scaling argument. In a semi-dilute suspension, the diffusive behaviour of squirmers is primarily governed by short-range two-body hydrodynamic interactions that induce significant particle reorientation. This motivates a random-walk model, yielding, for example, $D_t\sim U^2_{\textit{sol}}t_{\textit{mfp}}/ (1-\langle \cos \Delta \theta \rangle )$ (Berg Reference Berg1993), where $U_{\textit{sol}}$ is the swimming speed of an isolated swimmer in the absence of interactions, $t_{\textit{mfp}}$ is the effective time duration between successive collisions and $\Delta \theta$ represents the angular deflection between successive collisions. The time scale $t_{\textit{mfp}}$ can be scaled using a simple kinetic-theory argument. A swimmer sweeps out a volume $V_{sw\textit{ept}} \sim \pi a^2 U_{\textit{sol}} t_{\textit{mfp}}$ in time $t_{\textit{mfp}}$ , and the number of particles within this volume is $N \sim 3\phi U_{\textit{sol}} t_{\textit{mfp}}/(4a)$ . Requiring that $N\sim 1$ gives $t_{\textit{mfp}} \sim (a/U_{\textit{sol}})\phi ^{-1}$ . As $U_{\textit{sol}}$ and $\langle \cos \Delta \theta \rangle$ are independent on $\phi$ , the random-walk model results in

(4.5) \begin{equation} D_t \sim a U_{\textit{sol}} \phi ^{-1}. \end{equation}

It is important to note that the motion of non-Brownian swimmers is intrinsically deterministic. However, repeated reorientations progressively erase their orientational memory, such that the random-walk framework provides an effective description of their diffusive behaviour. In the present study, we examine whether this scaling law remains valid when finite fluid inertia is included. To this end, we plot the diffusion coefficient $D_t$ as a function of $\phi$ for $\textit{Re}_S=0.01$ and 1 on a log–log scale in figure 8. The dashed line, representing the scaling $D_t/(aV_s) = 0.24\phi ^{-1}$ for pullers ( $\beta = 5$ ) in the zero- $Re$ limit proposed by Ishikawa & Pedley (Reference Ishikawa and Pedley2007a ), is included for reference. The numerical data for $\textit{Re}_S = 0.01$ and $\beta = 5$ closely follow the reference line, while the results for $\textit{Re}_S = 1$ continue to exhibit the inverse scaling of $D_t$ with $\phi$ . Our results demonstrate that, although fluid inertia significantly modifies the diffusion coefficient, it does not alter the underlying $\phi ^{-1}$ scaling in the presence of weak fluid inertia.

Figure 9.

Effect of volume fraction $\phi$ on the MSDs of (a) pullers and (b) pushers in the presence of non-negligible inertial effects ( $\beta = \pm 5$ , $\textit{Re}_S=10$ ). The insets show the corresponding long-time diffusion coefficient as a function of $\phi$ .

However, such scaling is expected to break down at higher $\textit{Re}_S$ . Figure 9 shows the MSD of swimmers ( $\beta \pm 5$ , $\textit{Re}_S=10$ ) with $\phi$ ranges from 0.025 to 0.1. The MSD of pushers decreases as $\phi$ is increased as expected, whereas it increases as $\phi$ is increased for pullers. In the insets, we plot the corresponding long-time diffusion coefficient as a function of $\phi$ , confirming these opposite trends. This difference in behaviour stems from the qualitative differences in the two-body interactions between pullers and pushers. In the case of pushers, the two swimmers immediately separate and move in different directions upon collision. In contrast, pullers tend to move together as a stable pair without separating, as shown in figure 3(c). This stable paired motion is disrupted by collisions with a third swimmer, causing the swimmers to scatter in different directions. Consequently, as the volume fraction increases, pairwise entrapment is weakened, leading to enhanced swimmer transport. These results indicate that the introduction of inertia qualitatively alters the diffusion characteristics of pullers and pushers.

In addition, we note that the transition between these two regimes for pullers is non-trivial. For intermediate Reynolds numbers ( $1\lt \textit{Re}_S\lt 10$ ), pairwise scattering and pairwise entrapment are both active, and their relative importance depends on both $\textit{Re}_S$ and $\phi$ . Consequently, no simple universal dependence of $D_t$ on $\phi$ is expected in this regime. A systematic investigation of this Reynolds-number-dependent transition requires extensive simulations over a broad range of Reynolds numbers and volume fractions, and is therefore left for future work.

5.1. Settling dynamics in the weak-to-moderate inertia regime

When swimmers are denser than the surrounding fluid, they experience vertical migration. We first focus on the weak-to-moderate inertia regime ( $\textit{Re}_S \leqslant 1$ ), which is representative of most micro- and meso-scale systems. As is shown in the following, even such modest fluid inertia, when coupled with sedimentation, is sufficient to fundamentally alter the collective migration and diffusion. The influence of stronger fluid inertia ( $\textit{Re}_S = 10$ ) is examined separately in § 5.4.

The vertical drift velocity is defined as

(5.1) \begin{equation} {U}_d=\lim _{t \rightarrow \infty } \frac {\langle {x}(\tau +t)-{x}(\tau )\rangle }{t}. \end{equation}

Figure 10 presents $U_d$ as a function of settling coefficient $\alpha$ for $\beta =\pm 5$ . The inclined dashed lines represent $U_d=-V_g$ , corresponding to a simple superposition of isotropic swimming and directional sedimentation. At low $\textit{Re}_S$ , both pullers and pushers settle downwards under gravity, with settling speeds slightly smaller than the Stokes settling velocity $V_s$ due to the hindered effect. Contrary to the common results at low Reynolds numbers, as shown in figure 10(a), at $\textit{Re}_S = 1$ the pullers exhibit almost no net vertical migration ( $U_d \simeq 0$ ) at low values of $\alpha \leqslant 0.5$ , indicating that settling pullers are, on average, able to counteract gravitational settling in the finite-inertia regime. However, as $\alpha$ further increases, the downwards drift velocity increases markedly, and at $\alpha = 1$ the settling speed far exceeds that observed in the $\textit{Re}_S = 0.01$ cases. This observation is clearly anomalous, as fluid inertia typically enhances hydrodynamic drag and therefore decreases the settling velocity. Pushers, on the other hand, display a more conventional behaviour. As shown in figure 10(b), at $\textit{Re}_S = 1$ the increased hydrodynamic drag reduces their settling speed relative to the low- $\textit{Re}_S$ cases.

Figure 10.

Vertical drift velocity of settling swimmers as a function of the settling coefficient: (a) $\beta =5$ and (b) $\beta =-5$ . The inclined dashed lines represent $U_d=-V_g$ .

Beyond their vertical drift velocity, settling swimmers exhibit long-time diffusive behaviours. We discuss separately the diffusion coefficients parallel ( $D_t^\parallel$ ) and perpendicular ( $D_t^ \bot$ ) to the gravity direction, which are computed as

(5.2) \begin{equation} D_t^\parallel =\lim _{t \rightarrow \infty } \frac {\big \langle \Delta r^2_\parallel (t)\big \rangle }{2t}, \quad D_t^ \bot =\lim _{t \rightarrow \infty } \frac {\big \langle \Delta r^2_ \bot (t)\big \rangle }{4t}. \end{equation}

Here, $\Delta r_\parallel ^2(t) = [\hat {x}(\tau + t) - \hat {x}(\tau )]^2$ , where $\hat {x}(t) = x - \int U_d\,\mathrm{d}t$ denotes the MSD in the vertical direction relative to the mean vertical displacement. Meanwhile, $\Delta r_\bot ^2(t) = \left \langle [y(\tau + t) - y(\tau )]^2 + [z(\tau + t) - z(\tau )]^2 \right \rangle$ represents the horizontal MSD.

The results for the diffusion coefficients are shown in figure 11. At low $\textit{Re}_S$ , both the vertical and horizontal diffusion coefficients vary little with the settling coefficient, indicating that sedimentation has only a weak influence on the diffusive dynamics in the low- $\textit{Re}_S$ regime. At $\textit{Re}_S=1$ , however, both $D_t^\parallel$ and $D_t^\bot$ deviate from the neutrally buoyant cases. For pushers, both the vertical and horizontal diffusion coefficients decrease as $\alpha$ increases. Pullers, however, exhibit a markedly different behaviour: the vertical diffusion coefficient increases significantly for $\alpha \gt 0.5$ , while the horizontal diffusion coefficient decreases rapidly with increasing $\alpha$ . In particular, $D_t^\parallel$ becomes more than one order of magnitude larger than that of the neutrally buoyant case for $\alpha \geqslant 0.75$ . This result reveals an extreme anisotropy in the swimmer diffusion for settling pullers with finite fluid inertia.

Our observations indicate that sedimentation plays a key role in swimmer transport in the finite-inertia regime, particularly for pullers. The origin of these behaviours is examined in the following subsections.

Figure 11.

Relative diffusion coefficients as a function of the settling coefficient for (a,b) pullers ( $\beta =5$ ) and (c,d) pushers ( $\beta =-5$ ): (a,c) $D^{\parallel }_t/D^{\parallel }_{t,0}$ and (b,d) $D^{\bot }_t/D^{\bot }_{t,0}$ . Here $D^{\parallel }_{t,0}$ and $D^{\bot }_{t,0}$ represent corresponding values at $\alpha =0$ .

Figure 12.

Dependence of the preferential orientation of swimmers on the settling coefficient. Mean orientation angle for (a) $\beta = 5$ and (b) $\beta = -5$ . (c) Normalised angular probability density function of the orientation for inertial swimmers at $\textit{Re}_S = 1$ and $\alpha = 0.5$ .

5.2. Preferential orientation

To gain insight into the mechanisms underlying the observed transport behaviour, we first examine the preferential orientation of swimmers in the presence of gravity. We define $\gamma$ as the angle between the orientation vector $\boldsymbol p$ and the positive $x$ -axis (see figure 1 b). Figure 12(a) shows the dependence of the mean orientation of pullers, quantified by $\langle \cos \gamma \rangle$ , on the settling coefficient $\alpha$ . At low $\textit{Re}_S$ , no clear preferential orientation is observed. At $\textit{Re}_S=1$ , however, pullers exhibit a positive $\langle \cos \gamma \rangle$ , indicating a tendency to swim upwards against gravity. The corresponding results for pushers are shown in figure 12(b), where only a much weaker upwards preferential orientation is observed at $\textit{Re}_S=1$ .

We introduce a normalised angular probability density function, defined as (Ishikawa & Pedley Reference Ishikawa and Pedley2007b )

(5.3) \begin{equation} p^{\prime }(\gamma )=\frac {2}{n V \sin (\gamma )} \iint _{\gamma =\textit{const}} P(\boldsymbol{p}) \,\mathrm{d} A_e \,\mathrm{d} V, \end{equation}

where $P(\boldsymbol{p})$ is the probability that there is a swimmer with an orientation vector $\boldsymbol p$ within the solid angle $\textrm {d} A_e$ . If one assumes isotropic orientation of swimmers, then $p^{\prime }(\gamma )=1$ for all $\gamma$ . The results of $p^{\prime }(\gamma )$ for swimmers ( $\textit{Re}_S=1$ , $\alpha =0.5$ ) are shown in figure 12(c). Our results reveal that the $p^{\prime }(\gamma )$ distribution has a peak at $\gamma =0$ for $\textit{Re}_S=1$ . Meanwhile, $p^{\prime }(\gamma )$ has a larger peak for smaller $\phi$ due to weaker interactions among swimmers.

The gravitactic behaviour observed in figure 12 for non-bottom-heavy swimmers with sedimentation is reminiscent of that reported in bottom-heavy swimmers, which align vertically due to a gravitational torque arising from mass asymmetry (Ishikawa & Pedley Reference Ishikawa and Pedley2007a ). However, the underlying mechanisms are fundamentally different. In the finite- $\textit{Re}_S$ regime, the alignment of swimmers is governed by a hydrodynamic torque exerted on them during sedimentation, rather than from an intrinsic gravitational torque. For isolated spherical squirmers, Candelier et al. (Reference Candelier, Qiu, Zhao, Voth and Mehlig2022) has demonstrated that sedimentation of squirmer generates hydrodynamic torques arising from convective fluid inertia. These torques tend to align swimmers in the opposite direction of gravity, offering a unique mechanism of gyrotaxis. In Appendix A, we perform additional simulations of individual settling swimmers, which confirm that at $\textit{Re}_S=1$ , the puller experiences a substantially larger torque than the pusher.

The preferential orientation explains why puller suspensions with $\textit{Re}_S=1$ exhibit little net settlement in relatively small $\alpha$ , as seen in figure 10(a). The preferential orientation also explains the hindered diffusion of settling pushers at $\textit{Re}_S=1$ shown in figure 11(c,d). As the pushers become preferentially aligned against gravity, their orientations become much less dispersed. Consequently, neighbouring pushers tend to swim in nearly the same direction, which reduces their relative velocity and the frequency of close hydrodynamic encounters. The weakened hydrodynamic interactions lead to less-frequent reorientation and therefore suppress the long-time translational diffusion. However, preferential orientation fails to explain the anomalous vertical migration and diffusion of pullers at $\textit{Re}_S = 1$ and $\alpha \gt 0.5$ . In § 5.3, we provide a detailed analysis of the microstructure of the swimmer suspension, which offers further insights into these unusual transport phenomena.

5.3. Inhomogeneous structures

Sedimentation of swimmers will induce various microstructures and collective organisation, which can substantially alter the transport phenomena of swimmers. In order to quantify the microstructure of suspensions, we examine the pair correlation function $g(\boldsymbol{r}_{\!g})$ (Kulkarni & Morris Reference Kulkarni and Morris2008), which represents the normalised probability of finding a particle within a finite-volume element centred at the relative separation vector $\boldsymbol r_g$ ,

(5.4) \begin{equation} g(\boldsymbol{r}_{\!g})=\frac {P_{1 \mid 1}(\boldsymbol{r}_{\!g} \mid \mathbf{0})}{n}, \end{equation}

where $n=N/L^3$ is the average number density and $P_{1 \mid 1}(\boldsymbol{r}_{\!g} \mid \mathbf{0})$ is the conditional probability density of finding a particle at $\boldsymbol r_g$ , given that a reference particle is located at the origin $\boldsymbol 0$ . In the numerical implementation, $P_{1 \mid 1}(\boldsymbol{r}_{\!g} \mid \mathbf{0})$ is estimated by counting swimmer pairs within finite volume elements.

Figure 13.

Pair correlation functions $g({r}_{\perp },x_g)$ of swimmers at volume fraction of $\phi =0.025$ : (a–d) $\alpha =0$ ; (e–h) $\alpha =1.0$ . For reference, the pair correlation functions of passive suspensions are shown: (i) $\textit{Re}_G=0.01$ and ( j) $\textit{Re}_G=1$ .

For suspensions in the absence of gravity, particle pairs are expected to be isotropically distributed, resulting in a radially symmetric pair correlation function. In contrast, in sedimenting suspensions, previous numerical results indicated that the settling of passive particles has anisotropic microstructures (Yin & Koch Reference Yin and Koch2007; Xia et al. 2022, 2024, Reference Xia, Lin and Yu2025a ). Since the horizontal distribution remains approximately isotropic, the relative position vector $\boldsymbol r_g = (x_g,y_g,z_g)^T$ can be decomposed into vertical and horizontal components,

(5.5) \begin{equation} \boldsymbol{r}_{\!g}=\boldsymbol{r}_{\perp }+x_g \boldsymbol{e}_{\!x}, \quad \text{with } \quad \boldsymbol{r}_{\perp }=y_g \boldsymbol{e}_y+z_g \boldsymbol{e}_z, \quad r_{\perp }=\left \|\boldsymbol{r}_{\perp }\right \|=\sqrt {y_g^2+z_g^2} . \end{equation}

As a result, the three-dimensional pair correlation function $g(\boldsymbol{r}_{\!g})$ is averaged over the azimuthal direction and represented as $g({r}_{\perp },x_g)$ . Figure 13 shows the pair correlation functions $g({r}_{\perp },x_g)$ for suspensions of both neutrally buoyant swimmers ( $\alpha =0$ ) and settling swimmers ( $\alpha =1$ ) for $\phi =0.025$ . Figure 13(a–d) confirm that the microstructure of neutrally buoyant swimmers is nearly isotropic, with a high probability of finding swimmer pairs in near contact in all directions. In contrast, a higher pair correlation is observed for settling swimmer pairs with vertical separation (figure 13 e–h). In particular, as shown in figure 13( f), $g({r}_{\perp },x_g)$ exhibits pronounced values along the vertical direction for pullers with $\textit{Re}_S=1$ , indicating a column-like structure. Such columnar clustering has been previously reported in dilute passive suspensions at high settling Reynolds numbers ( $\textit{Re}_G \gg O(10)$ ), where it originates from wake-induced hydrodynamic attraction (Uhlmann & Doychev Reference Uhlmann and Doychev2014; Moriche et al. Reference Moriche, Hettmann, García-Villalba and Uhlmann2023; Jiang et al. Reference Jiang, Xu and Zhao2024). However, this mechanism is not applicable here, since the present regime is restricted to $\textit{Re}_G \leqslant 1$ , where wake effects are weak. To further isolate the underlying mechanism, we examine passive particles at $\textit{Re}_G = 0.01$ and 1, respectively (figure 13 i, j). In both cases, no columnar structures are observed. At $\textit{Re}_G=1$ , we observe a decreased pair distribution function for pairs with vertical separation, which is consistent with previous studies (Yin & Koch Reference Yin and Koch2007; Xia et al. Reference Xia, Yu, Pan, Lin and Guo2022). Taken together, these results indicate that the observed column-like structures arise uniquely in suspensions of settling pullers, and cannot be attributed to either classical wake-induced attraction.

Figure 14.

Velocity correlation functions $I_{U'}({r}_{\perp },x_g)$ of settling pullers ( $\beta =5$ , $\textit{Re}_S=1.0$ , $\phi =0.025$ ): (a) $\alpha =0.5$ , (b) $\alpha =0.75$ and (c) $\alpha =1.0$ .

Pullers exhibits a much stronger tendency to reorient against gravity than pushers owing to a stronger fluid-inertial torque. Although the restoring torque originates from fluid inertia rather than bottom-heaviness, the resulting orientational dynamics are qualitatively analogous to those of gyrotactic swimmers. An upwards-oriented puller simultaneously settles under gravity, generating a disturbance flow associated with its sedimentation. Equivalently, the sedimentation-induced disturbance flow establishes a local hydrodynamic shear that tends to rotate neighbouring swimmers. Due to a balance between fluid-inertial torque and shear-induced torques, this self-generated shear favours the reorientation and gradual migration of pullers towards the descending flow region, thereby promoting vertically aligned clustering.

While the above results demonstrate the emergence of anisotropic, column-like microstructures for settling pullers in the finite-inertia regime, they do not fully capture the coherence of the swimmer motion. To quantify the coherent motion, we compute the following pair correlation functions for the velocity fluctuations:

(5.6) \begin{equation} I_{U^{\prime }}\left (\boldsymbol{r}_{\!g}\right )=\left \langle U_x^{\prime }\left (\boldsymbol{x}_i\right ) U_x^{\prime }\left (\boldsymbol{x}_j\right ) \mid \boldsymbol{r}_{i j}=\boldsymbol{r}_{\!g}\right \rangle , \end{equation}

where $U'_x(\boldsymbol x_i)=U_x(\boldsymbol x_i)-\langle U_x\rangle$ is the $x$ component of the velocity fluctuation of swimmer $i$ . Figure 14 shows $I_{U'}({r}_{\perp },x_g)$ for settling pullers at $\textit{Re}_S=1$ with increasing settling. Figure 14(a) indicates that the velocity fluctuations are nearly decorrelated at all relative positions between swimmer pairs at $\alpha =0.5$ . As $\alpha$ increases, a pronounced positive correlation emerges at small horizontal separations, while a negative correlation is observed at separations larger than $10a$ (figure 14 c). These results indicate that stronger sedimentation and column-like structure will induce large-scale, vertically coherent flows in puller suspensions with finite inertia. The velocity correlations further demonstrate that the columnar microstructures are accompanied by coherent downwelling motions extending over many swimmer diameters, indicating that the swimmer accumulation is dynamically coupled with the surrounding disturbance flow. As shown in figure 15(a), such a downwelling flow generates vorticity that reorients the surrounding swimmers towards the direction of the downwelling flow. The increased local concentration within the columns further strengthens the descending disturbance flow. This mutual reinforcement establishes a positive feedback between swimmer accumulation and coherent downwelling, thereby sustaining the columnar organisation. Consequently, pullers preferentially reside within coherent downwelling regions than in a spatially homogeneous suspension, resulting in a significantly enhanced mean settling velocity. Taken together, the resulting collective dynamics can be interpreted within the framework of gyrotactic focusing (Kessler Reference Kessler1985, Reference Kessler1986; Pedley & Kessler Reference Pedley and Kessler1992; Croze et al. Reference Croze, Bearon and Bees2017; Rühle et al. Reference Rühle, Zantop and Stark2022; Ishikawa et al. Reference Ishikawa2025a ), in which the competition between fluid shear and a restoring torque drives the preferential accumulation of negatively buoyant swimmers. Although the restoring torque in the present system originates from fluid inertia rather than bottom-heaviness, the resulting collective organisation is qualitatively analogous to classical gyrotactic focusing.

Figure 15.

(a) Illustration for the columnar structure and coherent flow in settling puller suspensions. The background colour represents the streamwise fluid velocity fluctuation, where blue corresponds to downwelling and red corresponds to upwelling. Probability density functions of instantaneous swimmer ( $\beta =5$ , $\textit{Re}_S=1.0$ , $\phi =0.025$ ) velocities in the (b) $x$ and (c) $y$ directions for different values of $\alpha$ .

The same coherent organisation also accounts for the anisotropic diffusion. As shown in figure 15, the coexistence of coherent downwelling within the columns and compensating upwards flow in the surrounding dilute regions broadens the distribution of streamwise swimmer velocities (figure 15 b), thereby increasing the vertical velocity variance. Conversely, because swimmers remain laterally confined within columns due to gyrotactic focusing, transverse relative motion is suppressed, leading to a reduced horizontal velocity variance (figure 15 c). Consequently, the diffusion becomes increasingly anisotropic with increasing settling strength, characterised by enhanced transport in the vertical direction and reduced transport in the horizontal directions. The coherent collective dynamics associated with this gyrotactic-focusing-like state therefore provide a unified explanation for the enhanced settling velocity and strongly anisotropic diffusion observed in finite-inertia puller suspensions.

We note that the column-like structures and large-scale disturbance flows are constrained by the finite computational domain with periodic boundary conditions. Consequently, when inhomogeneous structures emerge (see Appendix B), the quantitative values of the drift velocity and diffusion coefficient exhibit some dependence on the domain size. Nevertheless, the qualitative conclusions regarding the formation of microstructures, the enhanced drift, as well as the anisotropic diffusion remain unaffected.

Figure 16.

Dynamics of settling pullers ( $\beta =5$ and $\phi =0.025$ ). (a) Vertical drift velocity as a function of the settling coefficient at $\textit{Re}_S=10$ . (b) Velocity autocorrelation. (c) Pair correlation function with $\alpha =1$ . For reference, (d) shows the pair correlation function of passive suspension with $\textit{Re}_G=10$ .

5.4. Settling dynamics in the high-inertia regime

The preceding subsections demonstrated that moderate fluid inertia ( $\textit{Re}_S = 1$ ) promotes settling-induced preferential orientation and collective organisation for pullers. A natural question is whether these collective states are further enhanced at higher Reynolds numbers, where isolated pullers experience an even stronger inertial torque. To address this issue, we additionally consider settling pullers ( $\beta =5$ and $\phi =0.025$ ) at $\textit{Re}_S = 10$ . In figure 16(a) we plot the vertical drift velocity and the mean orientation angle as a function of settling coefficient. Surprisingly, in contrast to the moderate-inertia regime, the suspension no longer develops any collective ordering, and the settling speed increases with increasing settling coefficient. The disappearance of preferential orientation at $\textit{Re}_S = 10$ cannot be attributed to a weaker fluid-inertial aligning torque. On the contrary, isolated pullers exhibit even stronger upwards alignment at this Reynolds number. We attribute the suppression of collective migration to the persistent pairwise entrapment characteristic of pullers at $\textit{Re}_S \sim O(10)$ . Figure 16(b) compares the velocity autocorrelation function $C_U(t)$ of settling pullers at $\textit{Re}_S=1$ and 10. At $\textit{Re}_S=10$ , the oscillatory velocity autocorrelation remains prominent when settling is introduced, indicating that the pairwise dynamics identified in § 4 remain the dominant hydrodynamic interaction. Although the fluid-inertial torque continuously drives swimmers toward an upwards-swimming state, repeated pairwise encounters prevent this alignment from being maintained over long times. Consequently, sustained gyrotactic focusing cannot be established, leading to the disappearance of collective organisation.

Figure 16(c) further displays the pair correlation function of pullers with $\textit{Re}_S=10$ and $\alpha =1$ , confirming that the pronounced vertically aligned structure observed in figure 13( f) disappear at higher Reynolds numbers. For comparison, we also include the pair correlation function of passive particles at $\textit{Re}_G=10$ in figure 16(d). In both cases, the probability of finding particle pairs in near contact is reduced, which suggests that the pairwise interactions destroy the vertical structure, so that, after averaging, the resulting distribution resembles that of passive particles settling under gravity.

Consequently, the collective migration identified at moderate Reynolds numbers no longer develops in the high-inertia regime. The non-trivial phenomenon highlights the importance of investigating the Reynolds-number dependence of collective migration of swimmers under gravity. We plan to explore the state transition in future work.