2.1. Soft microswimmer and problem settings

Consider microswimmers with a radius of $a$ swimming freely in an incompressible Newtonian liquid with a viscosity of $\mu$ and a density of $\rho$ . Due to the small size of the microswimmers, the particle Reynolds number is much smaller than unity, as such the inertial effect of the fluid flow can be neglected. The flow field is then governed by the Stokes equation. The microswimmer consists of a two-dimensional hyperelastic thin membrane enclosing a Newtonian liquid with the same viscosity and density as the external liquid. A torque swimmer model (Ishikawa et al. Reference Ishikawa, Tanaka, Imai, Omori and Matsunaga2016) was originally developed to represent a ciliate. The model expresses the thrust force at the tip of a cilium and the reaction force at the root of the cilium as a torque on the surface $A_t$ at a distance of $\varepsilon$ from the membrane (cf. figure 1 a). The magnitude of the propulsive torque density per unit area in the reference state $L_r$ is given by

(2.1) \begin{align} L_r = L_0 + \kappa \left ( \theta - \frac {\pi }{2} \right ), \end{align}

where $\kappa$ determines the torque strength of the asymmetry and $\theta$ is the polar angle with respect to the orientation vector $\boldsymbol {e}$ .

The torque density $L_0$ is defined as the value at which the swimming speed of the rigid body is equal to $U_0$ . The swimmer in the reference state satisfies the torque-free condition on the overall torque surface due to the axisymmetric torque distribution. The instantaneous torque density per unit area in the deformed state $L$ is given by

(2.2) \begin{align} L = L_r\frac { \mathrm{d} \it A_{t,r}}{ \mathrm{d} \it A_t} , \end{align}

where $ \mathrm{d}\it A_{t,r}$ and $ \mathrm{d} \it A_t$ are the local torque surface area in the reference and deformed states, respectively. This correction means that the number of cilia does not change with deformation. The propulsive torque in vector form is given by $\boldsymbol {L} = L\boldsymbol {n} \times \boldsymbol {t}$ , where $\boldsymbol {n}$ is the unit outward normal vector and $\boldsymbol {t}$ is the unit tangential vector inthe $\theta$ direction on the membrane surface. The soft swimmer is deformed from its reference shape by the propulsive torque distribution, resulting in a swimming-induced stresslet $\unicode{x1D64E}$ (Matsui et al. Reference Matsui, Omori and Ishikawa2020b ):

(2.3) \begin{align} \unicode{x1D64E} = \int { \left [ \tfrac {1}{2}(\boldsymbol {q}\boldsymbol {x} + \boldsymbol {x}\boldsymbol {q}) - \tfrac {1}{3}\boldsymbol {q} \boldsymbol{\cdot} \boldsymbol {x}\unicode{x1D644} \right ]} \mathrm{d}\kern1pt \textit{A}, \end{align}

where $\boldsymbol {q}$ is the traction at the point $\boldsymbol {x}$ on the membrane $A$ and $\unicode{x1D644}$ is the identity matrix. By changing $\kappa$ in (2.1), various swimming modes can be expressed as puller, neutral and pusher types (figure 1 b–d), respectively. In the rigid squirmer (Blake Reference Blake1971), the swimming mode is often expressed by positive and negative values of the second mode $\beta$ :

(2.4) \begin{align} v_{\theta } = U_0 \left ( \tfrac {3}{2}\sin \theta + \tfrac {3}{2}\beta \sin \theta \cos \theta \right ). \end{align}

Here $v_{\theta }$ is the squirming velocity on the squirmer surface. In this study the swimming mode is also defined using $\beta$ obtained from $\unicode{x1D64E}$ with reference to the rigid squirmer model, i.e.

(2.5) \begin{align} {S}_{ee}=4\pi \mu a^2\beta , \end{align}

where ${S}_{ee}$ is the stresslet component in the swimming direction.Here $\kappa$ is adjusted to ensure that, under the same $\beta$ , the stresslet produced by the rigid squirmer and the soft swimmer are the same. The correspondences between $\kappa$ and $\beta$ are almost linear as shown in figure 2. In this study, $\beta$ is set to −0.9 for pushers, 0 for neutral swimmers and 0.9 for pullers.

Figure 2.

The relationship between the asymmetry of torque strength $\kappa$ , the stresslet component in the swimming direction ${S}_{ee}$ and the swimming mode $\beta$ of the soft swimmer with each $Ca$ . Here $Ca$ is the ratio between the viscous force due to swimming and the elastic force of the membrane defined by (2.22).

Assuming an infinite suspension, triply periodic boundary conditions were applied, as shown in figure 1(e). Twenty-seven soft swimmers with random orientations were placed in random initial positions in the computational unit domain. The volume fraction $\phi$ was varied from 0.03 to 0.2 by changing the domain size. According to the scaling by Batchelor & Green (Reference Batchelor and Green1972), the frequency of $n$ -body interference is proportional to $\phi ^n$ . Thus, for $\phi = 0.2$ , the influence of three-body interference is expected to be 20 % of two-body interference.

2.2. Membrane mechanics

The coupling of thin-shell theory of capsules with fluid mech anics has been developed by the work of Barthès-Biesel, Pozrikidis and their co-workers. Therefore, we briefly explain the outline of thin-shell theory in the main text. For further details, please refer to Barthès-Biesel et al. (Reference Barthès-Biesel, Diaz and Dhenin2002), Lac et al. (Reference Lac, Barthes-Biesel, Pelekasis and Tsamopoulos2004), Walter et al. (Reference Walter, Salsac, Barthes-Biesel and Tallec2010) and Pozrikidis (Reference Pozrikidis2010). The swimmer’s membrane is modelled as a hyperelastic thin material. A material point of the membrane in the reference state is represented by $\boldsymbol {X}$ ; $\boldsymbol {x}(\boldsymbol {X},t)$ expresses the deformed state. Using a surface deformation gradient tensor $\unicode{x1D641}_s$ , reference and deformed states have the following relationship (Pozrikidis Reference Pozrikidis2010):

(2.6) \begin{align} \boldsymbol {dx} = \unicode{x1D641}_s \boldsymbol{\cdot} \boldsymbol {d \it X}, \end{align}

where $\unicode{x1D641}_s = \partial \boldsymbol{x}/\partial \boldsymbol{X}$ . The Green--Lagrange strain tensor $\unicode{x1D640}$ describing local deformation of the membrane is as follows:

(2.7) \begin{align} \unicode{x1D640} = \tfrac {1}{2} \big( \unicode{x1D641}_s^T \unicode{x1D641}_s - \unicode{x1D644} \big). \end{align}

The first and second invariants of the strain tensor are given by the principal extension ratios $\lambda _1$ and $\lambda _2$ as

(2.8) \begin{align} I_1= \lambda _1^2 + \lambda _2^2 -2 \end{align}

and

(2.9) \begin{align} I_2 = \lambda _1^2 \lambda _2^2 -1 = J^2_s - 1 , \end{align}

where $J_s$ = det( $\unicode{x1D641}_s)=\lambda _1 \lambda _2$ is the Jacobian expressing the area dilatation. Because the membrane is modelled as a two-dimensional hyperelastic material, the Cauchy tension $\unicode{x1D64F}_{\boldsymbol{s}}$ and the strain energy function $w_s(I_1, I_2)$ have the following relationship:

(2.10) \begin{align} \unicode{x1D64F}_{\boldsymbol {s}} &= \frac {1}{J_s} \unicode{x1D641}_s \boldsymbol{\cdot} \frac {\partial w_s}{\partial \unicode{x1D640}} \boldsymbol{\cdot} \unicode{x1D641}^T_s . \end{align}

In order to express strain hardening behaviour and the incompressible areal dilatation of a biological membrane, we utilise the constitutive law proposed by Skalak et al. (Reference Skalak, Tozeren, Zarada and Chien1973). The strain energy function of the Skalak model is given by

(2.11) \begin{align} w_s = \frac {G_s}{4} \big( I^2_1 + 2I_1 - 2I_2 + CI^2_2 \big) , \end{align}

where $G_s$ is the surface shear elastic modulus and $C$ is the non-dimensional areal dilatation constant. The Poisson ratio of the Skalak law is $\nu _s = C/(1+C)$ . The value of $C$ is set to 10 with reference to the erythrocyte model (Omori et al. Reference Omori, Ishikawa, Imai and Yamaguchi2013).

The bending energy of the biological membrane, derived by Helfrich (Reference Helfrich1973), is

(2.12) \begin{align} w_b = \frac {E_b}{2}\int \left ( 2H-c_0 \right ) \mathrm{d}\kern1pt \textit{A}, \end{align}

where $H$ is the mean curvature of the membrane, $c_0$ is the spontaneous curvature and $E_b$ is the bending modulus. By assuming a symmetric membrane property (i.e. $c_0 = 0$ ), the energy can be converted to the force density $\boldsymbol {q}_b$ as (Matsui et al. Reference Matsui, Omori and Ishikawa2020a )

(2.13) \begin{align} \boldsymbol {q}_b = -2E_b\big[\Delta _sH + 2H(H^2 - K)\big]\boldsymbol {n}, \end{align}

where $\Delta _s$ is the Laplace--Beltrami operator and $K$ is the Gaussian curvature on the membrane surface.

2.3. Fluid mechanics

Due to the small size of the swimmer, the inertial effect of fluid flow is negligible. Therefore, the flow field is governed by the Stokes equation, and the fluid velocity $\boldsymbol {v}$ at an arbitrary point $\boldsymbol {y}$ is given by the boundary integral equation (Pozrikidis Reference Pozrikidis1992) as

(2.14) \begin{align} \boldsymbol {v}(\boldsymbol {y}) = - \frac {1}{8\pi \mu }\sum _i^N \left [ \int _{A} \unicode{x1D645}(\boldsymbol {x},\boldsymbol {y}) \boldsymbol{\cdot} \boldsymbol {q}(\boldsymbol {x}) \mathrm{d}\kern1pt \textit{A}^i + \int _{A_{t}} \unicode{x1D64D}(\boldsymbol {x},\boldsymbol {y}) \boldsymbol{\cdot} \boldsymbol {L}(\boldsymbol {x}) \mathrm{d} \kern1pt A^i_t \right ] , \end{align}

where $N$ is the number of microswimmers in the fluid domain. Here, $\unicode{x1D645}$ and $\unicode{x1D64D}$ are the Green’s function of the Stokeslet and torque in the triply periodic boundaries, respectively (Beenakker Reference Beenakker1986):

(2.15) \begin{align} \unicode{x1D645} = \unicode{x1D645}^{1} + \unicode{x1D645}^{\kern1pt 2}, \end{align}

(2.16) \begin{align} \unicode{x1D64D} = \unicode{x1D64D}^{1} + \unicode{x1D64D}^{\kern1pt 2}, \end{align}

(2.17) \begin{align} \begin{bmatrix} {J}^{1}_{i j} \\[4pt] {J}^{\kern1pt 2}_{i j} \end{bmatrix} = (\delta _{i j} \nabla ^{2}-\nabla _{i} \nabla _{j}) \begin{bmatrix} r\:\mathrm {erfc}(\xi r) \\[4pt] r\:\mathrm {erf}(\xi r) \end{bmatrix}, \end{align}

(2.18) \begin{align} \begin{bmatrix} {R}^{1}_{i j} \\[4pt] {R}^{\kern1pt 2}_{i j} \end{bmatrix} = \tfrac {1}{4} \varepsilon _{l k j}(\nabla _{k} - \nabla _{l}) \begin{bmatrix} {J}^{1}_{i l} \\[4pt] {J}^{\kern1pt 2}_{i k} \end{bmatrix}, \end{align}

where $\boldsymbol {\delta }$ is the Kronecker delta, $\boldsymbol {\varepsilon }$ is the Levi-Civita epsilon and $\xi$ is an arbitrary positive constant with dimensions of inverse length for controlling the convergence rate of the error function and complementary error function. For a simple cubic lattice, the optimal value of the convergence rate that minimises the computational cost for the periodic expansion is given by $\sqrt {\pi }/V^{1/3}$ (Pozrikidis Reference Pozrikidis1992), where $V$ is the volume of the computational domain. The traction $\boldsymbol {q}$ is given by the stress jump across the thin membrane as

(2.19) \begin{align} \boldsymbol {q}(\boldsymbol {x}) = \left [ \boldsymbol {\sigma }_{out}(\boldsymbol {x}) - \boldsymbol {\sigma }_{in}(\boldsymbol {x}) \right ] \boldsymbol{\cdot} \boldsymbol {n}(\boldsymbol {x}) , \end{align}

where $\boldsymbol {\sigma }_{in}$ and $\boldsymbol {\sigma }_{out}$ represent viscous stress on the surface from inside and outside of the membrane, respectively.

2.4. Numerical methods

In this section we briefly explain the numerical method. The swimmer’s membrane and torque surface are respectively discretised by 1280 triangular elements, with 642 vertices, as shown in figure 1(a). Each material point is tracked in terms of its associated Lagrangian, and the strain at any time $t$ is determined. The tension is calculated based on Skalak’s constitutive law, and the load $\boldsymbol {q}$ acting on the membrane is obtained from the weak-form equilibrium equation, given as

(2.20) \begin{align} \int _A \hat {\boldsymbol {u}} \boldsymbol{\cdot} \boldsymbol {q}_s \mathrm{d}\kern1pt \textit{A} = \int _A \boldsymbol{\hat{\varepsilon }} : \unicode{x1D64F}_s \, \mathrm{d}\kern1pt \textit{A}, \end{align}

where $\hat {\boldsymbol {u}}$ is the virtual displacement and ${\boldsymbol{\hat{\varepsilon}}} = \ ({1}/{2})(\nabla _s \hat {\boldsymbol {u}} + \nabla _s \hat {\boldsymbol {u}}^T)$ is the virtual strain. The above equation is solved using a finite-element method (Walter et al. Reference Walter, Salsac, Barthes-Biesel and Tallec2010) with respect to unknown $\boldsymbol {q}_s$ . Assuming a linear sum of the loads due to in-plane deformation and the loads due to bending deformation, the load on the membrane surface is determined by $\boldsymbol {q} = \boldsymbol {q}_s + \boldsymbol {q}_b$ . When $\boldsymbol {q}$ is obtained, it is substituted into the boundary integral equation; the velocity $\boldsymbol {v}$ is computed using a Gaussian quadrature method with six numerical integration points in each element. A material point on the torque surface $\boldsymbol {x}_t$ is given by $\boldsymbol {x}_t = \boldsymbol {x}+\varepsilon \boldsymbol {n}(\boldsymbol {x})$ , and the flow due to the torque $\boldsymbol {L}$ is calculated in the same manner. In order to coincide with the half-length of the cilia, $\varepsilon$ is set to $0.05a$ , where $a$ is the radius of the swimmer. Assuming a no-slip boundary at the membrane, $\mathrm{d}\kern1pt \boldsymbol{x}/ \mathrm{d}\it t = \boldsymbol {v}$ , the material point $\boldsymbol {x}$ is updated by a second-order Runge–Kutta method. The above cycles are iteratively calculated to simulate hydrodynamic interactions in the suspension.

A multipole expansion technique is introduced to accelerate the computation (Pozrikidis Reference Pozrikidis1992). Depending on the distance between swimmers’ centres, three regions corresponding to (i) short-range, (ii) middle-range and (iii) long-range interactions are considered in the calculation. For the short-range interaction ( $r \leqslant r_{sm}$ ), the boundary integral equations are numerically integrated using all 1280 elements. In the middle-range interaction ( $r_{sm} \lt r \leqslant r_{ml}$ ), a coarse-grained 20 triangle mesh is used for the calculation; the point stresslet and torque defined at each mass centre are utilised for long-range interactions ( $r \gt r_{ml}$ ). Threshold values $r_{sm} = 4a$ and $r_{ml} =V^{1/3}$ were set to preserve computational accuracy.

The overlap between meshes becomes problematic when the swimmers are very close together (i.e. of the order of the mesh size). Assuming repulsion between cilia, the repulsive force at very short distances is defined as

(2.21) \begin{align} \boldsymbol {F}_{rep}(\boldsymbol {r}^{\prime}) = \begin{cases} k (2.4\varepsilon - |\boldsymbol {r}^{\prime}|) \frac {\boldsymbol {r}^{\prime}}{|\boldsymbol {r}^{\prime}|} & (|\boldsymbol {r}^{\prime}| \leqslant 2.4\varepsilon ), \\[4pt] \boldsymbol {0} &(|\boldsymbol {r}^{\prime}|\gt 2.4\varepsilon ) , \end{cases} \end{align}

where $\boldsymbol {r}^{\prime}$ is the distance between node points on the torque surface of each swimmer and $k$ represents the magnitude of the repulsion. In this study, $k/\mu U_0 = 500$ was set to be consistent with the ciliary elasticity estimated from the bending stiffness of the cilia (Katoh et al. Reference Katoh2023).

Membrane deformability is normalised by the capillary number $Ca$ :

(2.22) \begin{align} \textit {Ca} = \frac {\mu U_0}{G_s} . \end{align}

This dimensionless number represents the ratio between the viscous force due to swimming and the elastic force of the membrane. All computations were normalised by using the radius $a$ , swimming speed $U_0$ and surface shear elastic modulus $G_s$ . Here $\phi$ was changed by varying the length, whereas the number of swimmers remained fixed at 27.

3.1. Orientational correlations

For a more quantitative analysis, we define self-correlation of the swimming direction from the initial direction as

(3.1) \begin{align} I_{e_0} (t) = \left \langle \boldsymbol {e}_i(t) \boldsymbol{\cdot} \boldsymbol {e}_i(0) \right \rangle _{N} , \end{align}

where the bracket $\left \langle \right \rangle _N$ represents an ensemble average over the number of swimmers $N$ , $\boldsymbol {e}_i$ is the orientation vector of the $i$ th swimmer defined from the posterior and anterior node points. The result is shown in figure 4(c), in which the correlation gradually decreases over time. The relaxation time, corresponding to the time when $I_{e_0} = 1/e$ ( $e$ is Napier’s constant), was approximately estimated in the range of $100 \leqslant tU_0/a \leqslant 200$ for $Ca = 0.03$ , $0.1$ and $0.2$ , whereas that of the rigid puller (Ishikawa & Pedley Reference Ishikawa and Pedley2007) was approximately $tU_0/a \simeq 50$ . With respect to soft swimmers, the scattering angles after face-to-face interaction are reduced according to deformability (Matsui et al. Reference Matsui, Omori and Ishikawa2020a ). Thus, the self-correlation to the initial swimming direction persists for a longer duration as $Ca$ increases, and the swimmers retain their initial orientations longer compared with rigid swimmers. We also investigated the self-correlation of pusher and neutral swimmers (cf. figures 4 a and 4 b). Pusher and neutral soft swimmers tend to have longer relaxation times relative to rigid squirmers, as do pullers.

Figure 4.

Self-correlations of swimmers’ orientations $I_{e_0}(t) (=\langle \boldsymbol {e}_i(t) \boldsymbol{\cdot} \boldsymbol {e}_i(0) \rangle )$ in the suspension with $\phi = 0.1$ . Here $\boldsymbol {e}_i$ is the swimming direction of the $i$ th swimmer. The coloured lines and area shows the averaged values and standard deviations for the three independent cases throughout the paper. The black lines show the result of rigid squirmmers reported by Ishikawa & Pedley (Reference Ishikawa and Pedley2007).

In order to investigate the time scale for microstructure development in the suspension, the global order correlation is defined as

(3.2) \begin{align} I_e(t)=\left \langle \boldsymbol {e}_i(t) \boldsymbol{\cdot} \boldsymbol {e}_j(t) \right \rangle _N. \end{align}

To ensure that they would not affect the final structure, two initial conditions were set: (i) an initial polar order and (ii) an initial random configuration. The global correlation of the swimming directions reached a developed state when $tU_0/a \gt 100$ for pushers (cf. figure 5 a) and $tU_0/a \gt 200$ for pullers (cf. figure 5 b). These times were close to the relaxation times of the self-correlation; the converged values of $I_e$ were almost 0 for pushers and 0.2 for pullers. These results suggest that pullers achieve collective swimming, in which the swimmers are more aligned with each other’s swimming direction relative to pushers, similar to the behaviour of rigid squirmers (Ishikawa & Pedley Reference Ishikawa and Pedley2007).

Figure 5.

Time change of the global correlation of swimmers’ orientations $Ie(t)(=\langle \boldsymbol {e}_i(t) \boldsymbol{\cdot} \boldsymbol {e}_j(t) \rangle )$ with $\phi = 0.1$ and $Ca = 0.1$ .

Figure 6.

Correlations of swimmers’ orientations as a function of relative distance $I_e(t)(=\langle \boldsymbol {e}_i(r_i) \boldsymbol{\cdot} \boldsymbol {e}_j(r_j \rangle )$ with the $\phi = 0.1$ in the developed state. The black lines show the result of rigid squirmers.

We then investigated the spatial correlations of the swimmers in the developed state. Based on the results in figure 5, the criterion for determining a sufficiently developed microstructure was set to $tU_0/a = 200$ , and the spatial correlation considered the time average between $tU_0/a = 200$ and $400$ . For neutral swimmers, structural development was very slow and convergence values could not be obtained beyond $tU_0/a=1000$ . Thus, in the following sections we focus on the comparison between pushers and pullers; the results for neutral swimmers can be found in the Appendix. The spatial correlation function of the swimmers’ orientations within the radial interval $r_g -\Delta r$ to $r_g$ is defined as

(3.3) \begin{align} I_e(r_g) = \langle \boldsymbol {e}_i(\boldsymbol {r}_i) \boldsymbol{\cdot} \boldsymbol {e}_j(\boldsymbol {r}_j) \rangle _{N,t} , \end{align}

where $r_g$ is the relative centre--centre distance of the swimmers, $\boldsymbol {r}_i$ is the mass centre of the $i$ th swimmer and the mass centre of the $j$ th swimmer within the interval is given by $\boldsymbol {r}_j \in (r_g - \Delta r) \leqslant |\boldsymbol {r}_j-\boldsymbol {r}_i| \leqslant r_g$ . In this study, $\Delta r$ is defined as a constant that divides the domain length into 60 sampling points. The results are shown in figure 6. For the comparison, the results of rigid squirmers (Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008) are also presented. Ishikawa et al. (Reference Ishikawa, Locsei and Pedley2008) simulated three-dimensional suspensions with 64 rigid squirmers applied with periodic boundary conditions and then reported that hydrodynamic interactions form long-ranged orientational correlations of the domain size level ( $r_g/a \approx 7$ ). Negatively correlated interactions were not observed among rigid swimmers, whereas high $Ca$ swimmers showed negative correlations. Negative areas extended to the long range, especially among pushers. The correlations of soft swimmers maintain up to the domain size level, similar to the rigid swimmers. Since it is difficult to increase the domain size further due to computational constraints, innovations in computers and computational methods are needed to discuss larger scales. These results suggest that deformability enables swimmers to easily pass each other (e.g. during face-to-face interactions). The effects of membrane deformation are therefore expected to weaken the swimming-induced microstructure formation.

3.2. Swimmer aggregation

In order to quantify swimmer aggregation, the radial distribution function $g(r_g)$ is introduced as

(3.4) \begin{align} g(r_g) = \frac {\langle n(r_g)\rangle _{N,t}}{n_{mean}} , \end{align}

where $n_{mean} = N/V$ is the global mean of the number density. The local number density $n(r_g)$ is defined by the number of swimmers within the spherical shell domain of radius $r_g - \Delta r$ to $r_g$ .

The aggregation ratio $g$ of pushers and pullers is shown in figure 7. Both pushers and pullers tended to aggregate within $r_g \approx 3a$ in all cases, a tendency similar to the rigid squirmers (Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008). However, due to the membrane deformability, soft swimmers aggregated closer than $r_g=2a$ , which is equivalent to the contact distance of rigid swimmers. Higher capillary numbers allowed swimmers to aggregate more closely, but the peak value of $g$ was considerably reduced. Concerning rigid swimmers, the values were high very close to swimmers with $r_g\simeq 2a$ and converged to an average value of 1 at a distance of $r_g \geqslant 3a$ . The number of swimmers forming a cluster is estimated by integrating $\left \langle n(r_g)\right \rangle$ in the region where the swimmer density is grater than the mean. The results are summarised in table 1; in all cases, the clusters were formed by approximately two to three swimmers.

Figure 7.

Swimmers’ aggregations in the suspension with $\phi = 0.1$ . The black broken line indicates $g = 1.0$ , where the local number density is equivalent to the global number density. The black sold lines show the result of rigid squirmers (Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008).

Table 1.

Average number of swimmers in a unit cluster. The marker ‘*’ shows reference to the previous results of rigid squirmer models provided by Ishikawa et al. (Reference Ishikawa, Locsei and Pedley2008).

3.3. Swimming speed

The ensemble average of the swimming speed $\left \langle |\boldsymbol {U}_i|\right \rangle _N(t)$ was computed from $200 \leqslant tU_0/a \leqslant 400$ . The time-averaged value is represented by $\bar {\langle U\rangle }$ . Figure 8 shows $\bar {\langle U\rangle }$ as a function of $\phi$ , with $Ca = 0.03, 0.1$ and $0.2$ . For both pushers and pullers, $\bar {\langle U\rangle }$ monotonically decreased according to $\phi$ . For the same number density and $Ca$ , the average swimming speed of pushers exceeded that of pullers. This is because pushers can more easily pass each other. Among pusher swimmers, a higher $Ca$ value leads to faster swimming under the same number density conditions.

Figure 8.

Time-averaged swimming speeds of the swimmers with various $Ca$ . The error bar indicates the standard deviations computed from independent three initial conditions.

Translational dispersion is a measure of microswimmer spreading. The mean square displacement $R^2$ was used to determine the swimmer’s self-dispersion, as given by

(3.5) \begin{align} R^2(\Delta t) = \big\langle | \boldsymbol {r}_i(t + \Delta t) - \boldsymbol {r}_i(t)|^2 \big\rangle _N, \end{align}

where $\boldsymbol {r}$ is the swimmer’s displacement, $\Delta t$ is the time interval and $R^2$ is a function of $\Delta t$ with $\phi = 0.20$ , as shown in figure 9. We observed the different effects of $Ca$ between pushers and pullers. Among pushers, the value was higher for larger $Ca$ in the high $\Delta t$ regime, whereas it was almost invariant for $Ca$ among pullers. In the double logarithmic axes (cf. inset graph in figure9), both pushers and pullers had a slope of 2. Thus, the swimmers’ movement was ballistic, rather than diffusive, in the range of $\Delta t U_0/a \leqslant 100$ . A longer $\Delta t$ analysis may ultimately lead to a transition to diffusion; however, the computational time interval $\Delta tU_0/a \leqslant 100$ is not sufficient to confirm this transition. Since the time scale of the relaxation time for orientational self-correlation is estimated as $O(10^2)$ , a time scale of $O(10^3)$ may be required to observe a diffusion phenomenon, which is difficult to achieve due to the long computational time. The converged values could not be confirmed in this study.

Figure 9.

Mean square displacements $R^2$ in the suspension with $\phi = 0.2$ during the time interval $\Delta t$ . The inserted graph is drawn to a double logarithmic scale.

3.4. Rotational diffusivity

Rotational diffusivity $\unicode{x1D63F}_R$ is defined as follows to describe scattering in the swimming direction:

(3.6) \begin{align} \unicode{x1D63F}_R(\Delta t) = \lim _{\Delta t\rightarrow \infty } \frac {\left \langle ( \boldsymbol {\omega }(t + \Delta t) - \boldsymbol {\omega }(t))\otimes ( \boldsymbol {\omega }(t + \Delta t) - \boldsymbol {\omega }(t))\right \rangle _N}{2\Delta t},\\[-24pt]\nonumber \end{align}

(3.7) \begin{align} \boldsymbol {\omega }(t + \Delta t) - \boldsymbol {\omega }(t) = \int ^{\Delta t}_{t} \boldsymbol {\Omega } \mathrm{d}\it t. \end{align}

Here $\boldsymbol {\Omega }$ is the swimmer’s rotational velocity and $\boldsymbol {\omega }$ is the rotational displacement. The average of the diagonal components of $\unicode{x1D63F}_R$ represents the isotropic contribution:

(3.8) \begin{align} \unicode{x1D63F}_R = ({D}_{R,xx} + {D}_{R,yy} + {D}_{R,zz})/3. \end{align}

The time change of $\unicode{x1D63F}_R$ with $\phi = 0.10$ and $0.20$ is shown in figure 10. Here $\unicode{x1D63F}_R$ linearly increases in the small $\Delta t$ regime ( $\Delta t U_0/a \lt 1$ ) and then gradually converges to a specific value at $\Delta t U_0/a \geqslant 10$ . This tendency is similar to that of rigid swimmers (Ishikawa & Pedley Reference Ishikawa and Pedley2007). Regarding soft pushers, the converged value was slightly smaller with increasing $Ca$ , which suggests that softer pushers can maintain their swimming direction during the interaction, (i.e. the softer pusher has highly direct swimmability). Among pullers, $\unicode{x1D63F}_R$ was slightly higher with $Ca$ . In addition, we observed long and stable two-body clusters during face-to-face collisions, and $\unicode{x1D63F}_R$ convergence became worse (cf. figure 10 b).

Figure 10.

(a–d) Time change of swimmers’ rotational diffusivities during the time interval $\Delta t$ ; (a,b) $\phi = 0.1$ and (c,d) $\phi = 0.2$ . Converged values of $D_R$ in the suspension with (e) $\phi = 0.1$ and (f) $\phi = 0.2$ . (Note that the result of the puller type with $Ca = 0.2$ is excluded because the convergent value has not been obtained during the time scale that has been simulated.)

To further analyse the short-range interactions, the effects of swimming mode on the nearby swimmer distribution were examined. Because most clusters were formed by two or three swimmers (cf. table 1), we focused on the closest two-body interaction. The distributions of the closest swimmers during the short-range interaction were quantified, and the probability density $P$ was derived as a function of $\alpha$ . Here $\alpha$ is the relative angle between the closest swimmers $i$ and $j$ , (i.e. the angle in which the closest swimmer $j$ is from the direction of $\boldsymbol {e}_i$ ):

(3.9) \begin{align} \alpha = \cos ^{-1} \left ( \frac {\boldsymbol {e}_i \boldsymbol{\cdot} (\boldsymbol {r}_j - \boldsymbol {r}_i)}{|\boldsymbol {r}_j - \boldsymbol {r}_i|} \right ). \end{align}

As shown in figure 11, pushers with low $Ca$ exhibited asymmetric distributions of $P(\alpha )$ , and the peak appeared in the range of $\pi /3 \lt \alpha \lt \pi /2$ . As $Ca$ increases, the peak shifted to $\pi /2$ and the distribution became almost symmetric. With respect to pullers with high $Ca$ , a large peak was observed at $\alpha \lt \pi /12$ , indicating that a soft puller is likely to collide with the head. Among pullers, structures with heads coalesced due to the retracted flow have also been reported for the ellipsoidal squirmers (Kyoya et al. Reference Kyoya, Matunaga, Omori and Ishikawa2015); the results of the present study are consistent with this finding.

Figure 11.

Probability density of the nearest neighbours as a function of relative angle $\alpha$ from the swimming direction ( $\phi = 0.1$ ).

4.1. Strain

In order to quantify membrane deformation, we utilised the second invariant of the strain tensor $I_2$ within the short-range interaction ( $r \lt 3a$ ). The second invariant represents the area dilatation of the membrane, which controls the mechanosensitive channel of prokaryotes during severe osmotic challenges (Sukharev et al. Reference Sukharev, Blount, Martinac and Kung1997; Perozo et al. Reference Perozo, Cortes, Sompornpisut, Kloda and Martinac2002). The time- and ensemble-averaged strain increment $\Delta I_{2}$ relative to the solitary swimming is presented in figure 12. As indicated in figure 12(a), $\Delta I_{2}$ was nearly proportional to $Ca$ for both pushers and pullers. This phenomenon represents a direct consequence of an increase in $Ca$ , resulting in reduced membrane rigidity. Here $\Delta I_{2}$ was greater for pullers than for pushers, due to the pulling effect in the swimming direction. We also investigated the effect of the volume fraction $\phi$ on deformation (cf. figure 12 b). A linear trend was observed for the volume fraction. A greater volume fraction led to more frequent microswimmer collisions, resulting in larger deformation.

Figure 12.

Difference of the second strain invariant from the solitary swimming: (a) effect of $Ca$ ( $\phi = 0.1$ ) and (b) effect of $\phi$ ( $Ca = 0.1$ ).

4.2. Tension

We evaluated the principal tension during the interaction. The principal tension $\tau _1$ of the Skalak law is given by

(4.1) \begin{align} \tau _{1} = \frac {G_s \lambda _1}{\lambda _2}(\lambda ^2_1 - 1 + C\lambda ^2_2 I_2) \end{align}

and $\tau _{2}$ is obtained by switching $\lambda _1$ and $\lambda _2$ . Because the in-plane isotropic tension $T$ , defined as $T = (\tau _1 + \tau _2)/2$ , could be relevant for activation of the mechanosensing ion channel (Matsui et al. Reference Matsui, Omori and Ishikawa2020a ), we focused on the in-plane isotropic tension $T$ .

In this study the reference shape of the elastic membrane was assumed to be a sphere. Membrane tension was then generated, even in dilute suspensions without interparticle interactions. The mean value of $T$ in the solitary swimming, $T_{mean}^{sol} = {\int (T { \mathrm{d}}A)}/{A}$ , is shown in figure 13(a), where $A$ is the membrane surface area of the swimmer. Here $T_{mean}^{sol}/\mu U_0$ increased according to $Ca$ . The tension of pushers was always higher than that of pullers because pusher swimmers displayed a dented shape in the back (cf. figure 1), and tended to have a higher $T_{mean}^{sol}$ compared with other swimmer types.

Figure 13.

In-plane isotropic tension of swimmers $T (= (\tau _1 + \tau _2)/2)$ . Tensions are derived in each triangular element and then averaged over swimmers and times. (a) Mean tension of the solitary swimming. (b,c) Mean tension normalised by the solitary swimming. (d,e) Time-averaged maximum tension in the suspension. The bar indicates instantaneous maximum and minimum values of the maximum tension in the suspension.

Considering variation in the tension due to the hydrodynamic interactions within the suspension, the averaged incremental tension from the solitary swimming $\langle T_{mean} \rangle _{N,t}/T_{mean}^{sol}$ was lower for pushers (cf. figure 13 b) because they are more likely to pass each other. As shown in figure 13(c), the slip-through effect of pushers was effective even at higher number densities; among the pullers, tension increased according to number density. Regarding pullers, tensions were strongly influenced by $Ca$ . When $Ca = 0.03$ , $\langle T_{mean} \rangle _{N,t}/T_{mean}^{sol}$ was approximately 3.2, and the value monotonically decreased according to $Ca$ .

The instantaneous maximum tension in the suspension was determined; the time averages are shown in figure 13. Both swimmers showed similar trends in $Ca$ and $\phi$ . As $Ca$ increased with fixed $\phi$ , the maximum tension slightly decreased. Softer swimmers were able to avoid hard collisions, potentially explaining why the maximum tension was reduced by $Ca$ , whereas the maximum tension linearly increased according to $\phi$ . The frequency of many-body interferences increased according to number density. Accordingly, the capsular membrane experienced greater tension.

The tension distribution changed with pushing and pulling hydrodynamic interactions. The incremental tension during the short-range interaction was evaluated as a function of the azimuth angle $\theta$ relative to the orientation vector, as defined by $\Delta T(\theta )$ in the following:

(4.2) \begin{align} \Delta T(\theta ) = \big\langle T(\theta ) - T^{sol}(\theta )\big\rangle _{N,t}. \end{align}

Here $T(\theta )$ is the circumferentially averaged isotropic tension at $\theta$ in the suspension and $T^{sol}(\theta )$ is that of solitary swimming. The result is shown in figure 14. In all cases, incremental tension $\Delta T(\theta )$ decreases according to $Ca$ (cf. figures 14 a and 14 b). Pushers showed a relatively flat distribution of $\Delta T(\theta )$ , except at the $\theta = \pi$ pole, whereas pullers tended to have a decreasing $\Delta T(\theta )$ in the $\theta \leqslant \pi /2$ region.

Figure 14.

Incremental tension $\Delta T$ as a function of the polar angle on a swimmer surface $\theta$ ( $\phi = 0.1$ ). (a,b) Distributions of local incremental tension $\Delta T$ from the solitary swimming. (c,d) Distributions of the ratio $\Delta T/T_{sol}$ .

Upon further examination of directional distributions that caused the largest tension fluctuations, the incremental tension was normalised according to the value used for solitary swimming. The results are shown in figures 14(c) and 14(d). Among pushers, the ratio of increase was high for $3\pi /4\lt \theta \lt \pi$ ; among pullers, the ratio of increase was high around $\theta = \pi /2$ . These tendencies were more pronounced at low capillary numbers. These results indicate that for pushers (pullers), the greatest change in tension occurs in the posterior (side) during the interaction.

Figure 15.

Orientational correlation of the neutral swimmers with $\phi = 0.1$ . Each line indicates three different initial conditions with (a,c,e) showing initially random configurations and (b,d,f) initially polar order configurations.