Disorder-induced persistent random motion and trapping of microswimmers

Abstract

Microorganisms often move in confined, disordered environments, where hydrodynamic couplings can modify their transport behavior. Using extensive finite-element simulations, we investigate the dynamics of microswimmers – modeled as squirmers – in two-dimensional disordered porous media by resolving the full hydrodynamic interactions. We reveal that the deterministic coupling between activity, hydrodynamics, and disorder is sufficient to generate effective diffusive transport. Strong pushers and pullers become localized in the porous medium either by trapping at corners or dynamic trapping, depending on swimmer type and the obstacle packing fraction. Squirmers can escape from dynamic traps, leading to prominent “hop-and-trap'' dynamics. Strikingly, we find a pusher-puller asymmetry in the trapping probability that can be reversed by short-range swimmer–obstacle interactions, highlighting the sensitivity of transport to near-field effects.

---

Converting chemical energy into directed motion allows microorganisms to optimize their survival strategies and to perform biological functions,^1–5 while synthetic microrobots – equipped with a self-propulsion mechanism – are expected to play a major role for future nano-technological applications, ranging from targeted-drug-delivery^6,7 to bioremediation.^8,9 At microscopic scales, the dynamics of these active microswimmers are dominated by viscous forces,¹⁰ imposing constraints on the swimming kinematics. Geometrical confinement of their surroundings further shape their transport behavior through hydrodynamic couplings and steric interactions.^11–13 Complex disordered media, omnipresent in natural and technological contexts, have been demonstrated to modify the motion of bacteria,¹⁴ worms,¹⁵ and active colloids¹⁶ to intriguing hop-and-trap dynamics, while truncating their directional persistence and controlling large-scale bacterial transport through repeated surface-scattering events.¹⁷ Revealing the physical ingredients underlying these non-equilibrium processes not only provides fundamental insights into microbiological systems but also paves the way towards novel synthetic smart microswimmers capable to adapt to their surrounding environment.^18–22

From a theory perspective, in the absence of hydrodynamic interactions transport of “dry” active particles through porous media has revealed localization in geometric pockets and corners driven by steric interactions.^23–29 Recent studies have demonstrated that self-propelled microswimmers optimize their diffusion through porous media via tuning their reorientation mechanisms, such as run-and-tumble or run-reverse patterns, allowing them to truncate their trapping phases in geometric pockets and to hop through the pore space.^{24–28,30,31} While this body of work highlights the importance of geometry, it typically neglects hydrodynamic couplings between a swimming microswimmer and nearby surfaces, which could play an important role in confined, complex environments.

Hydrodynamic interactions have been studied, both experimentally and theoretically, primarily for planar walls^32–39 and isolated spherical obstacles,^40,41 which have been demonstrated to profoundly influence swimmer dynamics, leading to scattering, stable swimming states, and trapping near boundaries. These effects depend strongly on the swimmer's shape, its hydrodynamic (such as pusher/puller) signature, and tumbling behavior, and can be strongly amplified by the properties of the confinement, such as repulsive potentials,^37,38 surface structure,⁴² or boundary deformability.⁴³ A paradigmatic model for describing microswimmers in this context is the squirmer model, which was first introduced by Lighthill⁴⁴ and later refined by Blake.^45,46 Within this framework, swimmer–surface interactions are known to depend sensitively on the prescribed surface-slip modes. In particular, it has been shown that hydrodynamic interactions alone can result in squirmer trapping at boundaries, alignment, or scattering from planar walls, as well as forward and backward orbiting around single obstacles.^37,47–53 Extensions to periodic lattices have further demonstrated how geometric confinement guides swimmer trajectories and induces localization.^54–57 While these studies offer valuable insights into swimmer-surface interactions in controlled settings, microorganisms are typically found in much more complex environments, where disorder is an intrinsic property of the system. In such settings, swimmer–boundary interactions, geometric disorder, and activity are intricately coupled, making it difficult to anticipate the transport behavior of swimmers from single-obstacle studies.

In this work, we isolate the effects of hydrodynamic interactions and spatial disorder on swimmer dynamics, while neglecting rotational noise or tumbling mechanisms of active microswimmers. Combining a phase-field approach with extensive finite-element simulations, we study the motion of squirmers through a disordered porous medium. Our results reveal that deterministic hydrodynamic interactions with disordered obstacles give rise to multiple distinct trapping mechanisms, including static confinement, quasi-periodic orbital states, and geometric pocket trapping. Moreover, we find a pronounced asymmetry between pushers and pullers, leading to qualitatively different survival probabilities and transport regimes. Even in the absence of noise, swimmers exhibit diffusive, sub-diffusive, or super-diffusive motion, demonstrating that disorder-induced hydrodynamic reorientation alone is sufficient to generat anomalous transport and localization.

1 The model

We consider a microswimmer navigating a disordered environment at low Reynolds number in two dimensions (2D). The swimmer is modeled as a disk-shaped squirmer of radius a, self-propelled by a prescribed tangential slip-velocity at its surface ∂S^44–46:where n denotes the unit surface normal and p is the swimming direction. Further, B₁ represents the first squirming mode, which sets the swim speed in free space U = B₁/2.⁵⁸ The swimmer type is characterized by the dimensionless squirming parameter β, measuring the stresslet strength, where β < 0 corresponds to pushers, β = 0 to neutral swimmers, and β > 0 to pullers. A sketch illustrating the setup is provided in Fig. 1.

Fig. 1 Squirmer sketch. Here, x_C is the squirmer center of mass, p denotes the squirmer orientation, θ corresponds to the angle between the squirmer orientation and the x-axis, and n is the unit normal to the squirmer surface pointing outward. Further, the slip velocity u_S is indicated in blue at the squirmer surface for β = 0. The line thickness indicates the slip velocity magnitude rescaled as 5||u_S||/U.

The disordered medium is represented by N randomly distributed non-overlapping disks of radius R⁵⁹ and characterized by the packing fraction ϕ = NπR²/L², where L is the length of the domain. To prevent overlap between the squirmer and the obstacles we apply a short-ranged repulsive potential ensuring that the swimmer-obstacle distance exceeds the minimal distance h_cut.

The spatially-varying fluid velocity and pressure fields, u(x, t) and Π(x, t), are governed by the quasi-steady incompressible Stokes equations:

μ∇²u = ∇Π and ∇·u = 0, (2)

with viscosity μ and no-slip boundary conditions on the obstacle surfaces. The velocity on the squirmer surface (∂S) entails its translational and angular velocities, U and ω, via: u = u_S + U + ωe_z × r (with e_z the unit normal to the 2D plane and r pointing from the swimmer center to its surface). The velocity and pressure fields and the velocities of the swimmer are obtained by solving the Stokes equations, subject to the force- and torque-free conditions on the swimmer:where is the stress tensor. We employ a phase-field approach for the squirmer and perform finite-element simulations with adaptive grid refinement,^60–65 taking into account periodic boundary conditions to model an effectively infinitely-large porous medium. The swimmer's instantaneous center of mass x_C(t) and orientation θ(t) are then updated according to the equations of motion: ẋ_C(t) = U and [small theta, Greek, dot above](t) = ω. Details of the numerical implementation and its validation with ref. 47 and 53 are provided in the Appendix.

2 Results

We perform extensive numerical simulations to investigate the dynamics of microswimmers, characterized by different squirming parameters β = −8, …8. Specifically, we set R = 4a and L = 128a, and select two cut-off values δ = h_cut/a = 1/20, 1/4, which can lead to different squirmer-obstacle interactions.³⁸ To resemble dilute and dense porous media, we consider two packing fractions, ϕ = 0.15 and 0.45. We remark that from the perspective of a squirmer of radius a interacting via a short-range repulsion of range δ, the relevant excluded radius becomes R_eff = R + a + δ. The corresponding effective packing fraction, accounting for overlaps of the enlarged disks, can be extracted numerically and yields ϕ_eff ≈ 0.69 for the dense systems considered here. Thus, although ϕ = 0.45 may appear moderately dilute, the swimmer effectively experiences a medium close to the percolation threshold.⁵⁹ We collect representative statistics by sampling 128 trajectories, starting at different locations with different initial orientations, for four different porous media realizations and each parameter set. The initial position is restricted to lie within the percolating cluster of the medium.

2.1 Microswimmer dynamics: from swimming to trapping

To elucidate the impact of the squirming parameter on active motion through a dense disordered porous medium (ϕ = 0.45), we consider trajectories starting from the same initial position and orientation. Distinct swimming patterns emerge, as shown in Fig. 2 (left) and in the supplementary movies. Neutral squirmers (β = 0) undergo frequent scattering from obstacles, leading to long, curved trajectories through the pore space. Increasing |β| > 0 results in trapping of squirmers in localized regions, whose nature depends on both the squirmer strength and the pusher/puller flow signature. For instance, a puller (β = 2) is scattered by several obstacles before it ends up statically trapped at a fixed position due to geometric confinement. Differently, a pusher with β = −4 becomes trapped and orbits around a single obstacle. Decreasing the squirmer parameter (β = −8) leads to its localization between two obstacles, where it follows a quasi-periodic orbit, as the increase in squirming parameter enables hydrodynamic interactions with a second nearby obstacle. Here, it is important to clarify that a quasi-periodic orbit is referred to a trajectory that remains spatially confined and repeatedly follows a similar path, but does not close exactly on itself. Instead, the swimmer returns slightly displaced from its initial position after each cycle. For instance, in the single-obstacle case the orbital motion is accompanied by small oscillations, so successive revolutions do not exactly repeat the same trajectory. Similarly, between two obstacles the swimmer undergoes multiple scattering between them while following a nearly identical path, yet each cycle terminates at a slightly shifted location, see Fig. 2 (right).

Fig. 2 Representative trajectories of squirmers in a disordered porous medium. (Left) Four swimmers with different squirming parameters β start at the same position and orientation (dark gray disk). The end points of the trajectory are marked by the red stars. Here, we set the packing fraction to ϕ = 0.45 and use the cut-off δ = 1/20. (Inset) A trajectory transitioning between dynamic trapping phases and free swimming motion. (Right) Examples of dynamically trapped trajectories exhibiting quasi-periodic motion for β = 0 (top) and β = −4 (bottom). Trajectory opacity makes repeated overlaps appear darker. In the former case, trapping is of topological origin, whereas in the latter it arises from hydrodynamic interactions. The shaded regions around the obstacles indicate the effective disk radius, accounting for both the squirmer size and the short-range repulsive interaction.

Inspection of our full dataset thus reveals two distinct states: exploration of the pore space and localization, either in the form of static trapping at obstacle boundaries or dynamic trapping in confined pockets. Transitions between the two states are also observed, see Fig. 2 (left-inset). Here, the squirmer first becomes dynamically trapped and orbits between two obstacles. As the orbits are not closed, the swimmer eventually escapes and explores the medium. Subsequently, it enters a second dynamic trap formed by more closely spaced obstacles, where it remains for the remainder of the simulation, spending significantly longer times there than in the first trap.

These different behaviors are also observed for the dilute porous medium, yet with different statistical signatures depending on the squirmer parameter, the flow signature, and the short-range repulsive potential.

2.2 A note on steric vs. hydrodynamic interactions

It is worth noting that in the absence of hydrodynamic interactions, the microswimmer cannot reorient in our model as the steric interactions are only repulsive. To isolate the role of hydrodynamic interactions from purely steric effects, we performed simulations in which hydrodynamics was switched off while keeping the porous geometry and short-range repulsive obstacle interaction unchanged. In this case, swimmers move with fixed speed B₁/2 and retain their initial swimming direction. The resulting dynamics is therefore only affected by the geometry. In dense environments (ϕ = 0.45), swimmers become statically trapped with probability one. In dilute environments (ϕ = 0.15), by contrast, the trapping probability is P ≈ 0.65 and the survival probability approaches a plateau S(t) ≈ 0.35, for t → ∞. The free trajectories display ballistic mean-squared displacement scaling, indicating that swimmers continue along nearly straight paths. The dependence on δ is negligible in both dense and dilute systems, demonstrating that in the no-hydrodynamics limit the trapping probability is controlled mainly by the arrangement of obstacles and the available pore space, rather than by the precise range of the steric repulsion. These control simulations are therefore consistent with a picture in which steric interactions alone produce geometric trapping, while hydrodynamic interactions are required for nontrivial reorientation and transport dynamics.

2.3 Exploring microswimmers: disorder-induced persistent random motion

To quantify the statistical properties of the different states emerging from hydrodynamic interactions, we first identify parameter combinations (ϕ, δ, β) that allow squirmers to explore the environment without being trapped. We analyze the motion of these microswimmers in terms of their mean-square displacement (MSD), 〈|Δr(t)|²〉, where Δr(t) = x_C(t) − x_C(0) is the displacement at time t. Irrespective of the parameter set, our results show that at short times the MSDs obey 〈|Δr(t)|²〉 ≈ U²t² due to their purely persistent swimming motion (Fig. 3). Once the microswimmers have traversed the center-to-center contact distance, t ≈ τ: = (R + a)/U, the MSDs start to deviate from a purely persistent behavior as hydrodynamic interactions with the randomly distributed obstacles impact their dynamics through modifying the swimming direction. Details of these interactions depend on the specific parameters.

Fig. 3 Exploring microswimmers. Mean-squared displacement 〈|Δr(t)|²〉/a² for microswimmers exploring the disordered medium. The left and right panel show results for δ = 1/20 and δ = 1/4, respectively. Further, ϕ denotes the packing fraction and β is the squirming parameter. The triangle marks the crossover time τ = (R + a)/U, corresponding to the time it takes the microswimmer to move the center-to-center contact distance.

In particular, for a small cutoff, δ = 1/20, and a dense porous medium, ϕ = 0.45, free trajectories are observed only for neutral squirmers (β = 0) (Fig. 3, left). These display a diffusive behavior at long times 〈|Δr(t)|²〉 ∼ t, indicating that successive encounters with obstacles effectively randomize their swimming direction. For dilute systems (ϕ = 0.15), the MSD displays a super-diffusive regime at long times for β ∈ {−1, 1, 2} with a clear trend: pullers exhibit more pronounced super-diffusive behavior than pushers. This occurs because puller interactions with obstacles lead to fewer tangential sliding events and less frequent reorientations, in contrast to their pusher-counterpart, allowing them to preserve an effectively persistent motion for extended times, see Appendix A.6.

By increasing the potential range (δ = 1/4), we find that moderate pushers (β = −2) in dilute systems (ϕ = 0.15) also remain untrapped. This occurs because the repulsive layer increases the minimum squirmer–obstacle gap: since hydrodynamic interactions are strongest at small gaps, preventing close approach weakens the near-field coupling between the swimmer and the obstacles (Fig. 3, right). The same qualitative trend is observed as for δ = 1/20: pullers exhibit more pronounced super-diffusive behavior than pushers, with 〈|Δr(t)|²〉 ∼ t² for β = 2, indicating a longer persistence of the swimming direction for pullers. In dense systems (ϕ = 0.45), neutral squirmers also display diffusive behavior while weak pushers exhibit slight sub-diffusion. This arises because hydrodynamic interactions with obstacles tend to align pushers tangentially to nearby surfaces. As a result, pushers spend longer times navigating close to obstacles, which reduces their effective displacement and slows their spatial exploration compared to a neutral squirmer.

We stress that the MSDs reported here characterize finite-time dynamics, not necessarily an asymptotic diffusive regime. Trajectories classified as free are those that do not become trapped within the simulation window. The coexistence, for the same values of β, δ, and ϕ, of trapped trajectories and trajectories that remain mobile over the observation window suggests that the latter get trapped at later times. Since the domain is finite and contains absorbing trapping configurations, continued exploration of the system can eventually lead to trapping. In the absence of noise or other escape mechanisms, the ultimate long-time behavior may therefore be localization rather than diffusion. Thus, the MSDs should be interpreted as finite-time measures of microswimmers mobility.

2.4 Localized microswimmers: dynamic and static trapping

Many microswimmers are unable to explore their environment for long times and instead become localized due to hydrodynamic coupling with obstacles and geometric confinement. Upon localization, we further characterize swimmer states through a confinement radius ρ (i.e. the radius of gyration), computed from the spatial extent of each trajectory after its trapping time t_trap(see Appendix). Swimmers with ρ < a are classified as statically trapped, corresponding to localization at a fixed position, whereas swimmers with ρ ≥ a are classified as dynamically trapped, exhibiting persistent motion within a confined region. Using this criterion, we measure the trapping probability P(β) as a function of the squirming parameter, shown in Fig. 4(a). To quantify the temporal aspects of this trapping process, we analyze the survival probability S(t), defined as the probability that a swimmer is not localized in the porous medium up to time t.

Fig. 4 Localized microswimmers. (a) and (d) Trapping probability P(β). Here, 's', 'd', and 't' denote static, dynamic, and total trapping events. (b) and (e) Survival probability S(t) as a function of time t for different squirming parameters β. (Inset) shows the rescaled median of the trapping time T̃. (c) and (f) Histogram of the normalized confinement radius ρ/(2(a + R)). (a)–(f) Different panels correspond to different packing fractions ϕ and repulsive potential strengths δ.

We identify three generic features that are robust across packing fractions ϕ and cutoff distances δ: (i) Static trapping occurs, in almost all cases, more frequently than dynamic trapping. This difference may follow from the different constraints underlying the two trapping mechanisms. Static trapping is mainly a geometric arrest: it occurs whenever neighboring obstacles form a constriction whose surface-to-surface gap is smaller than the effective squirmer diameter 2(a + δ). Such configurations are common at the packing fractions studied here. Dynamic trapping is more selective. It requires the squirmer to remain mobile while its center of mass follows a bounded, quasi-periodic trajectory for long times. This can happen through persistent orbiting around an isolated obstacle, repeated interactions between two suitably spaced obstacles, or motion inside a pocket-like obstacle cluster with only narrow escape routes. These situations require a more restricted range of obstacle separations and orientations than simple steric blockage. (ii) Higher squirmer strength leads to faster trapping. (iii) Neutral squirmers are the least likely to become trapped.

In fact, neutral squirmers can become statically trapped only by approaching a surface with an orthogonal orientation. Such a configuration is, however, unstable, since even a small perturbation in the swimming direction would cause the squirmer to scatter away from the obstacle, thus making trapping of neutral squirmers highly unlikely.

2.4.1 Trapping in dense environments. In dense media (ϕ = 0.45), confinement is significantly enhanced and pullers are trapped with high probability for β ≠ 0. Weak pushers, in contrast, can remain mobile provided the cutoff is sufficiently large (δ = 1/4), highlighting the screening effect of short-ranged repulsion (Fig. 4a). Consistently, the survival probability S(t) together with its median T show that pushers persist longer than pullers at the same squirmer parameter. A higher cutoff markedly increases the survival probability, with a particularly strong effect for weak pushers, whose S(t) become comparable to those of neutral squirmers (Fig. 4b).

Analysis of dynamically trapped trajectories reveals two distinct localization mechanisms (Fig. 4c). For large squirmer parameters, trapping typically involves only two nearby obstacles and is characterized by small confinement radii, indicating hydrodynamically stabilized states in narrow gaps (Fig. 2). Neutral squirmers instead display a qualitatively different behavior: their confinement radii are broadly distributed and trapping generally involves three or more obstacles, consistent with repeated scattering within the pore space; an example is displayed in Fig. 2 (right-top). Thus, while strong pushers and pullers are localized by hydrodynamic coupling to closely spaced obstacles, neutral swimmers and weak pushers become trapped only through purely geometric confinement.

2.4.2 Trapping in dilute environments. In dilute porous media (ϕ = 0.15), trapping depends sensitively on both the squirmer parameter β and the short-range cutoff δ (Fig. 4d). Here, δ controls the effective range of the short-range repulsive potential; the two values considered therefore represent distinct model regimes of near-field interaction. For small cut-offs (δ = 1/20), strong pushers and pullers (|β| ∈ {4,8}) are almost certainly localized, but through distinct mechanisms: pushers predominantly undergo dynamic trapping in quasi-periodic orbital states, whereas pullers are mainly statically trapped with vanishing translational and angular velocities. Pullers that are dynamically trapped are either stuck between two obstacles or orbit around one. For moderate squirmer strengths, trapping is primarily static and occurs more frequently for pushers than for pullers at the same squirmer magnitude.

Increasing the cutoff to δ = 1/4 screens near-field hydrodynamic interactions by increasing the minimum squirmer–obstacle separation and makes dynamic trapping between two obstacles rare. Strikingly, it also reverses the pusher–puller asymmetry: at fixed |β|, pullers become more likely localized than pushers. The trend is directly reflected in the median trapping times extracted from the survival probability (Fig. 4e, inset), which are larger for pullers than pushers at δ = 1/20 and larger for pushers than pullers at δ = 1/4. This behavior reflects the tendency of pullers to align orthogonally towards obstacle surfaces so that self-propulsion is counterbalanced by repulsion, whereas pushers preferentially adopt tangential sliding states that facilitate escape for larger cut-off where hydrodynamic interactions are effectively weaker. We emphasize that this does not imply that the sign of β is irrelevant at the level of individual swimmer–obstacle interactions. Indeed, our single-contact analysis for moderate squirmer strengths, β ∈ {−2, −1, 0, 1, 2}, shows that pushers and pullers approach and leave obstacle surfaces differently, with distinct contact orientations, contact times, and scattering angles; see Appendix A.6. The trapping statistics discussed here, however, are global observables resulting from repeated encounters with different obstacle configurations. Therefore, the pusher–puller asymmetry known from simpler geometries, such as channels³⁷ and obstacles,³⁸ can be modified by the disordered porous environment.

In dilute environments geometric confinement is very weak; therefore, hydrodynamic trapping occurs only for large values of |β|. The associated confinement radii distribution for β ∈ {−8, 4, 8} is shown in Fig. 4f, indicating frequent single-obstacle traps and trapping between two obstacles for small short-range cutoff δ. Increasing the latter leads to weaker near-field hydrodynamic interactions, which in turn reduces the probability of two-obstacle confinement.

2.4.3 Fluid dynamics of a dynamic trapping event. To elucidate the microscopic origin of dynamic trapping, we examine a representative trajectory of a strong pusher confined between two obstacles (Fig. 5). This example illustrates a two-obstacle localization mechanism that occurs repeatedly in our simulations and contributes to the multi-obstacle trapping events reported in Fig. 4c and f. For an isolated obstacle, a pusher swimming directly towards the surface cannot approach closer than a finite gap, where self-propulsion is balanced by hydrodynamic repulsion generated by the fluid pushed against the boundary. A slight tilt of the swimming direction destabilizes this head-on configuration, reorients the swimmer, and induces a lateral hydrodynamic attraction towards the obstacle, leading to a stable orbital state in which frontal repulsion and side-wise attraction balance. As a consequence of these hydrodynamic interactions, the squirmer orientation (p, blue arrows) and its swimming velocity (U, red arrows) can differ, see Fig. 5.

Fig. 5 Velocity fields during a dynamic trapping event. Snapshots of a strong pusher (β = −8) and small cutoff (δ_cut = a/20) illustrating dynamic trapping between two obstacles. The blue arrow denotes the squirmer' s orientation p, while the red arrow indicates the direction of its instantaneous velocity U. The small arrows indicate the fluid velocity direction.

The resultant orbit is not strictly circular but exhibits persistent oscillations. When the swimmer approaches the surface too closely and reorients toward it, near-field hydrodynamic repulsion causes it to bounce away. Since the swimmer is not perfectly orthogonal to the surface, lateral attraction subsequently draws it back, resulting in a quasi-periodic oscillatory motion around the obstacle. During such an orbit, the swimmer may enter the hydrodynamic attraction range of a second nearby obstacle. In this case, while still oriented towards the first obstacle, the rear flow of the pusher generates a repulsive interaction with the second surface, causing the swimmer to slide along it in the opposite direction. Once sufficient separation is regained, the swimmer detaches and returns towards the first obstacle, thereby resulting in persistent confinement between the two obstacles, with alternating interactions with both boundaries.

Summary and conclusions

Combining a phase-field approach with extensive finite-element simulations, we have studied the dynamics of squirmers in dilute and dense disordered porous media, accounting for the full hydrodynamic interactions between the swimmer and the obstacles. We have revealed the emergence of persistent random motion of neutral squirmers due to hydrodynamically-induced reorientation by the nearby obstacles. We have further identified both static and dynamic trapping events, where microswimmers quasi-periodically trace certain paths within confined pockets for significant amounts of time. Examples range from single-obstacle orbiting to orbiting between two obstacles to moving in pockets confined by multiple obstacles. We find that these trapping characteristics strongly depend on the squirmer strength, the repulsive potential between the squirmer and the obstacles, and the packing fraction of the medium. Interestingly, our results show transient trapping phases, leading to hydrodynamically-induced “hop-and-trap” motion, reminiscent to 'dry' active dynamics in porous media.²⁷

Our approach, focused on disk-shaped microswimmers, can be readily extended to other swimmer shapes. This may allow us to understand the interplay of cell morphology and hydrodynamics for transport through porous media, where experiments have revealed that cell length and deformability can selectively enhance exploration or promote dead-end trapping depending on pore disorder.^15,66 It further provides a promising framework for unraveling why the motion of E. coli becomes rectified in ordered obstacle arrays, while spherical Janus colloids orbit around obstacles and stochastically switch to neighboring ones.⁵⁴

Our study neglects rotational noise, which could promote escape from localized states⁴¹ and induce hop-and-trap dynamics. Incorporating this effect in crowded media, where nearby obstacles may modify the rotational diffusivity relative to its bulk value, represents an interesting direction for future work. Including tumbling mechanisms into our modeling approach, may allow us to elucidate the interplay of stochasticity in cell behavior and hydrodynamic couplings with the surrounding confinement.³⁹

An additional open research endeavor entails the change in behaviors from 2D systems, to quasi-2D channels, often used in microfluidic experiments, to three-dimensional dynamics. We speculate that our findings should remain qualitatively relevant in quasi-2D geometries, although quantitative differences are expected. In microfluidic channels, nearby boundaries can modify, and possibly screen, the hydrodynamic interactions responsible for scattering and trapping. By contrast, fully three-dimensional systems may differ qualitatively, since swimmers can evade planar trapping mechanisms, such as two-obstacle trapping, by escaping through the additional spatial dimension.

While we have considered a quiescent fluidic environment, where flows are solely generated by the propulsion mechanism of active microswimmers, many natural environments are subject to externally-imposed hydrodynamic flows.^3,67 Controlled microfluidic set-ups of swimming microorganisms and theoretical studies of active microswimmers in channel-geometries have revealed intriguing phenomena, such as rheotaxis (i.e. motion along/against a flow),^68–72 butterfly-like trajectories at channel constrictions,⁷³ and accumulation of microswimmers behind spherical obstacles.^74,75 In complex environments, however, the resultant flow and vorticity fields are highly heterogeneous,⁷⁶ leading to new effects. For example, these can induce curly swimmer trajectories,^77,78 trap point-like active microswimmers via vorticity-induced reorientation in open channels,⁷⁸ or generate shear-induced Taylor dispersion.⁷⁹ Yet, the ramifications of the swimmer's own size and hydrodynamic signature on transport across porous channels in the presence of an externally-imposed flow remain largely unexplored.

Beyond single-swimmer dynamics, our framework can be naturally generalized to account for hydrodynamic interactions between multiple squirmers. Collective effects are known to qualitatively modify transport properties in confinement, giving rise to clustering, enhanced diffusion, hydrodynamic bound states, and large-scale coherent motion.^80–82 In porous environments, microswimmers could clog channels and thereby redirect flows, which could impact other active microswimmers. Thus, the interplay of disordered porous media, hydrodynamic flows, and swimmer properties may produce exciting novel physical phenomena – bridging clogging and self-organization – which could have substantial implications for the formation of microbial communities and biofilms^83,84 and could inspire new tools for microswimmer sorting and filtration.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available upon reasonable request from the authors.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d6sm00308g.

Appendix

The appendix contains a detailed description of the phase-field model, its asymptotic analysis, and our numerical approach. We further present the validation of the model with known analytical results, provide details of the classification algorithm for the different trapping phases, and discuss squirmer-obstacle contact statistics.

A.1 Phase-field model

Let denote the region occupied by the swimmer and Ω\S the surrounding fluid domain. The governing equations read

∇·σ = 0, in Ω\S, (4)

∇·u = 0, in Ω\S, (5)

u = u_S + U + ωe_z × r, on ∂S, (6)

u = 0, on ∂Ω, (7)

where u and p denote the velocity and pressure fields, respectively, U and ω are the translational and angular velocities of the swimmer, and u_S is a prescribed tangential slip velocity on the swimmer surface. The stress tensor is with viscosity μ. We adopt a phase-field penalization formulation to avoid explicit tracking of the swimmer boundary, as required in immersed boundary approaches.⁸⁵ Our formulation follows the general framework of diffuse-domain and phase-field approaches, in which sharp boundaries are replaced by smooth indicator functions and interfacial conditions are imposed in a distributed manner; see ref. 86–90. In our model, the sharp interface ∂S is replaced by a smooth indicator function ψ, and the equations are extended to the entire domain Ω in the following way:

∇·σ = η(U + ωe_z × r + u_S − u), (8)

∇·u = 0, (9)

with the boundary condition u = 0 on ∂Ω. The penalization function is defined as

η = Kψ, (12)

where K ≫ 1 is a large constant. The indicator function is constructed from a smooth phase fieldwhere a is the swimmer radius and ξ the interface thickness. In practice one chooses K ∼ ξ⁻².

A.2 Asymptotic analysis

In this section, we use matched asymptotic analysis to show that eqn (8)–(11) converge to the sharp-interface eqn (3) and (4)–(7) as ξ → 0. In this approach, we expand the variables in powers of the interface thickness ξ in regions close to (inner expansion) and far (outer expansion) from the interface. The two expansions are matched in an intermediate region where both expansions are presumed to be valid (e.g., see^91,92 for a general description of the procedure).

Outer expansion. Let S: = {x:ϕ(x) ≥ 0.5}. We assume that the stress σ has a regular expansion in ξ in the fluid domain Ω\S, given by σ(ξ) = σ⁽⁰⁾ + ξσ⁽¹⁾ + ξ²σ⁽²⁾ + …. Integrating eqn (8) over Ω and equating it with eqn (10) gives

Using the divergence theorem and inserting the expansion of σ gives to leading order

Moreover, away from the squirmer, we have ϕ = 0 to all orders in ξ, and so ∇ϕ = 0 to all orders. Therefore, the contribution of η in eqn (8) becomes negligible outside the squirmer and the leading order of this equation yields ∇·σ⁽⁰⁾= 0 in Ω\S. We therefore readily recover the sharp momentum balance eqn (4).

Integrating this expression over Ω\S and using the divergence theorem gives

Subtracting eqn (16), we recover the sharp interface force balance (3). Similarly, for the torque balance, we may take the cross product of eqn (8) with r and integrate over Ω. Equating the result with eqn (11) gives

Using the divergence theorem and the symmetry of σ, and inserting its expansion gives to leading order

Moreover, away from the squirmer, a cross product of r with eqn (8) provides to highest order r × ∇·σ⁽⁰⁾ = 0 in Ω\S. Integrating over Ω\S and using the divergence theorem gives

Subtracting eqn (19), we recover the sharp interface torque balance (eqn (3)).

Inner expansion. It remains to be shown that the phase-field approach recovers the velocity boundary condition on the squirmer surface [eqn (6)]. Near this boundary, we introduce a local coordinate system

x(s,z;ξ) = X(s) + ξzn(s) (21)

where X(s) is a parametrization of the interface, z = r(x)/ξ is the stretched variable, and r is the signed distance from the point x to ∂S, which is taken to be positive inside the squirmer. We then assume that all variables can be written as functions of z and s and that in these coordinates the variables have regular expansions in ξ. That is, for the velocity field

û(z,s;ξ) ≡ u(x;ξ) = u(X(s) + ξzn(s);ξ) (22)

and the inner expansion is

û(z,s;ξ) = û⁽⁰⁾(z,s) + ξû⁽¹⁾(z,s) + ξ²û⁽²⁾(z,s) + … (23)

The nabla operator in the local coordinate system becomes

where ∇_Γ is the covariant surface derivative.

Inserting these local derivatives, we can consider the inner expansions of eqn (8). Assuming that the penalization constant K scales as ξ⁻², we find at highest order 1/ξ³:

∂_z[small phi, Greek, circumflex]⁽⁰⁾(Û⁽⁰⁾ + [small omega, Greek, circumflex]⁽⁰⁾e_z × r̂ + û_S − û⁽⁰⁾) = 0. (24)

Since ∂_z[small phi, Greek, circumflex]⁽⁰⁾ > 0 for all z, we recover eqn (6):

Û⁽⁰⁾ + [small omega, Greek, circumflex]⁽⁰⁾e_z × r̂ + û_S = û⁽⁰⁾. (25)

We may further take the derivative of this expression with respect to z to obtain to highest order

∂_zû_S = ∂_zû⁽⁰⁾, (26)

where we have used that ∂_zû⁽⁰⁾ = 0, ∂_z[small omega, Greek, circumflex]⁽⁰⁾ = 0, ∂_zr̂ = n̂ξ. Matching inner and outer solution in an overlap region where both expansions are valid, leads to the typical matching condition: u is continuous across the surface if ∂_zû⁽⁰⁾ = 0 (see⁹²). Hence, if the imposed slip velocity u_S is continuous across the squirmer surface, then we have global continuity of the resulting velocity field u.

A.3 Numerical approach

The phase-field equations are discretized using the finite-element method. The diffuse interface has thickness ξ and introduces steep gradients in the indicator function ψ and the penalization field η. Accurate resolution of this region requires locally refined meshes. Finite elements allow straightforward adaptive mesh refinement, enabling the interface region to be resolved while maintaining a moderate number of degrees of freedom in the bulk fluid, see Fig. 6.

Fig. 6 (Left) Adaptive mesh refinement used to resolve steep gradients and hydrodynamic interactions in the phase-field formulation. (Right) Flow field generated by a pusher (β = −8) in a confined geometry, rescaled by the free-space swimming speed.

We next derive the weak formulation and the corresponding discrete system. Let H₀¹(Ω) denote the usual Sobolev space. In the weak formulation we aim at finding (u,Π,U,ω) ∈ H₀¹(Ω) × L²(Ω) × [Doublestruck R]² × [Doublestruck R] such that

for all test functions (v, q). Here, we denote r^⊥ = e_z × r = (−r₂, r₁, 0).

Let denote the velocity basis functions and the pressure basis functions. The discrete velocity and pressure fields are

The resulting algebraic system takes the block form

where s = (U₁,U₂,ω)^T. The matrix entries are defined as

The rigid-body coupling matrix reads

where ê₁ and e₂ are the Cartesian basis vectors. The matrix M is given byand the right-hand side vectors are

The resulting linear system is implemented using the AMDiS framework⁶⁵ and solved with the PETSc linear algebra backend,^93,94 using the MUMPS solver.^95,96

A.3.1 Parallel implementation. The formulation is well suited to distributed-memory parallelization by MPI. Under domain decomposition, the standard Stokes blocks A and B are assembled entirely from element-wise contributions and therefore require only the usual local finite-element assembly. The rigid-body coupling block C can likewise be assembled from local contributions on each subdomain, since its entries are of the formwhere w_i denotes one of the rigid-body modes associated with translation or rotation. Hence C is obtained by standard parallel assembly of locally computed integrals.

The only genuinely global object is the small matrix M, whose entries involve integrals over the whole domain, e.g.:

Its evaluation therefore requires a global reduction across MPI ranks. However, M is of size 3 × 3 for a single squirmer in 2D (and remains block-structured and very small for multiple squirmers), so both its assembly and inversion are negligible compared to the solution of the Stokes system. In practice, the parallel cost is therefore dominated by the fluid solve, while the additional overhead associated with the rigid-body coupling is minimal.

A.3.2 Short-range hard-wall repulsive potential. Before advancing the squirmer center of mass x_c, a short-range contact correction is applied to avoid overlap and to produce physically meaningful near-wall motion. For each nearby obstacle, the closest point p to the squirmer center x_c is computed, withwhere a is the squirmer radius and h is the signed gap.

A short-range repulsive potential with prescribed cutoff h_cut is applied only when the squirmer is sufficiently close to the obstacle and approaching it, i.e. when h < h_cut and U·ĝ < 0.

The center velocity U is decomposed into normal and tangential parts:

U_⊥ = (U·ĝ)ĝ, U_‖ = U − U_⊥. (41)

When the above cutoff condition is satisfied, the update suppresses/limits the inward normal contribution and keeps the tangential component, so that motion transitions to sliding along the obstacle instead of penetrating it. In practice, this yields: (i) no geometric overlap, (ii) robust behavior in narrow gaps, and (iii) near-wall kinematics dominated by tangential motion.

A.3.3 Adaptive time-step update. After each solve, the time step is updated to keep the squirmer motion stable near obstacles, while avoiding unnecessary slowdown far from walls. The timestep Δt is computed from two ingredients: (i) an interface-based bound proportional to ξ/|U·ĝ|, and (ii) a gap-based bound proportional to g/|U·ĝ|. Let N be the number of neighboring obstacles to the squirmer and let Δt_ξ,k and Δt_g,k be the timestep restrictions from (i) and (ii), respectively. The final step is the most restrictive admissible value over nearby obstacles, then clipped by global lower/upper safety bounds:with additional damping that limits abrupt changes relative to Δtⁿ (e.g. growth/shrink factors). This provides small steps during close approach and larger steps in free-swimming regimes.

A.4 Validation

A.4.1 Squirmer in free space. As a first validation of the numerical method, we consider a single squirmer in free space. An analytical solution of the Stokes equations for a 2D squirmer can be derived using a stream-function formulation in polar coordinates. The solution provides explicit expressions for the velocity field both inside and outside the squirmer disk. Typically, for squirmer motions only the first two modes are retained (i.e. B_n = 0 for n > 2). In that case, the analytical solution inside the squirmer disk in the squirmer's reference frame is

whereas outside the squirmer disk the analytical solution is given by

Direct comparison with the analytical solution in the unbounded domain is not appropriate for numerical simulations, which are necessarily performed in bounded computational domains. Boundary conditions imposed at the outer boundary introduce deviations from the infinite-domain solution. For this reason, we validate the method by comparing the numerical velocity field with the analytical solution inside the squirmer disk, where both descriptions are consistent. Therefore, we define the relative error as

where u denotes the numerical velocity field and u_exact the analytical solution. Further, we compare the free space speed of the squirmer U = B₁/2 with the one computed numerically. The results are shown in Fig. 7a, showing nice convergence and relative errors of order for small interface thickness ξ.

Fig. 7 (a) Relative error with respect to the analytical solution for selected values of β. The results show slightly better than first-order convergence, consistent with the phase-field approximation. (Inset) Relative error of the squirmer speed in free space compared with the exact value U = B₁/2. Linear convergence with respect to the interface thickness ξ is observed, in agreement with the phase-field model. (b) Translational and angular velocities of a squirmer near a wall for different interface thicknesses ξ. The top row shows B₁ = 1, B₂ = 0, while the bottom row shows B₁ = 0, B₂ = 1. Black lines denote the analytical solution.⁵³ The numerical results converge toward the analytical prediction as ξ decreases.

A.4.2 Squirmer near a wall. We further validate our model by computing the instantaneous translational and rotational velocities of a squirmer near a no-slip wall and comparing them to the analytical results presented in ref. 53:where ρ = d/a−(d²/a²−1)^1/2, and d is the distance of the center of mass x_C to the wall.

We fix the distance of the squirmer' s center of mass to d = 1.2a from the wall and vary its orientations as θ ∈ {0, π/6, π/4, π/3, π/2}. For each orientation, we compute the translational and angular velocities for two sets of squirmer modes (B₁ = 1, B₂ = 0) and (B₁ = 0, B₂ = 1). These computations are repeated for decreasing values of the interface thickness ξ = a/10, a/20 and a/40. The results are shown in Fig. 7b, where we observe that decreasing ξ leads to numerical solutions that closely match the analytical predictions.

A.4.3 Squirmer interacting with one obstacle. We investigate the interaction between a squirmer of radius a and a circular obstacle of radius R = 4a. Simulations are performed in a square domain of size L = 64a, with no-slip boundary conditions imposed both on the obstacle surface and on the domain boundaries. To prevent overlap between the squirmer and the obstacle, we employ a short-range hard-wall repulsive potential characterized by a cutoff gap h_cut. Two values are considered: h_cut = a/20 and h_cut = a/4.

Consistent with the results reported in ref. 38, we observe two distinct dynamical regimes. For small values of |β|, the squirmer scatters from the obstacle, whereas for sufficiently large |β| it becomes hydrodynamically trapped and orbits around it (results not shown). In the orbiting regime, the dynamics depend on the swimmer type. Pullers (β > 0) exhibit forward orbits, with their orientation aligned with the direction of motion, and display periodic oscillations whose amplitude gradually decays. Pushers (β < 0), in contrast, show a behavior that depends on the cutoff gap h_cut. For the smaller gap h_cut = a/20, the swimmer tends to orient opposite to its direction of motion, resulting in a backward orbit. For the larger gap h_cut = a/4, the motion resembles that of pullers: after a short transient backward motion, the swimmer settles into a forward orbit. An example of these behaviors is shown in Fig. 8.

Fig. 8 Interaction of a pusher (β = −10) with a circular obstacle of radius R = 4a. (Left) Forward orbit obtained for a repulsive cutoff gap h_cut = a/4. (Right) Backward orbit obtained for h_cut = a/20. The squirmer is color coded according to the rescaled time tU/a, with a the radius of the squirmer and U = B₁/2 the free-space velocity. The black arrows indicate the squirmer orientation. Qualitatively similar behaviors were reported in ref. 38.

A.5 Classification of squirmer trajectories

Squirmer trajectories are classified according to the late-time behavior of their displacement over the simulation window, see Fig. 9. For each trajectory, we compute the displacement from the initial position,

d(t) = |x_C(t) − x_C(0)|,

and its upper envelope

Fig. 9 Top row. One-obstacle dynamic trapping (left), multi-obstacle dynamic trapping (middle), free trajectory (right). The dark gray small disk corresponds to the initial position of the squirmer. The shaded regions around the obstacles indicate the effective disk radius, accounting for both the squirmer size and the short-range repulsive interaction. Bottom row. Corresponding displacement relative to the initial position as a function of time for the trajectories shown above. Note the quasi-periodic pattern displayed by the displacement in the dynamical trapping states, which is absent in the free trajectory.

We identify candidate confined trajectories from the emergence of a plateau in D(t), detected from the local scaling exponent of the envelope. However, since D(t) is measured relative to the initial position, a plateau alone does not necessarily imply localization: a trajectory may stop increasing its distance from x_C(0) because it turns back or changes direction, while still exploring new regions of the domain. To avoid such false positives, once a candidate plateau onset time t_p is identified, we discard the earlier part of the trajectory and repeat the displacement-envelope analysis on the late-time segment t ∈ [t_p,t_f], using x_C(t_p) as the new reference position. The procedure is then repeated on progressively shorter final segments. A trajectory is classified as confined only if a plateau persists under these changes of the reference point, indicating that the motion remains bounded in space rather than merely bounded relative to the initial position. Trajectories for which no persistent plateau is detected are classified as free.

For confined trajectories, we further distinguish two regimes according to the spatial extent of the residual motion after the onset of confinement. Let t_p denote the detected onset time of the plateau and define

If d_max < a/10, where a is the squirmer radius, the squirmer effectively remains localized and the trajectory is classified as static trapping. Larger bounded amplitudes correspond to dynamic trapping, where the squirmer remains confined but continues to move within a finite region.

A.6 Contact event statistics

In this section we provide a direct comparison between pushers and pullers during contact events with obstacles. A squirmer is considered to be in contact with obstacles when the gap between them is sufficiently small:where x_Oi is the center of the i-th obstacle.

We focus on free trajectories in dilute systems. Along each trajectory, we identify contact events as the contiguous time intervals over which eqn (47) holds. For each contact event, we measure:

(1) The contact time T_contact, defined as the duration of the contact interval,

(2) The reorientation associated with the encounter, quantified by 1 − cosΔθ, where Δθ is the scattering angle between the incoming and outgoing orientations,

(3) The event-averaged contact angle α_contact between the squirmer orientation and the local tangent to the obstacle surface.

The results are summarized in Fig. 10, where we report the median of the distributions obtained by pooling all contact events over the free trajectories for each value of β and δ. Note that no free trajectories are available for δ = 1/20 and β = −2, and therefore this data point is omitted. We emphasize here that the results are for squirmer-obstacle interactions that do not lead to an immediate trapping event.

Fig. 10 Squirmer-obstacle interactions. (Top) Median of contact times T_contact, (middle) reorientation Δθ, and (bottom) average squirmer alignment with respect to the obstacle surface α_contact as a function of the squirming parameter β. Error bars indicate in interquartile range, from the 25th to the 75th percentile.

Our results show that the median contact times T_contact are of comparable magnitude for pushers and pullers. They are, however, systematically larger for pushers, indicating that pushers typically spend longer times in the vicinity of obstacles. We further find that the reorientation associated with obstacle encounters is stronger for pushers. In particular, for β = 2, pullers are only weakly reoriented by contacts with obstacles. Finally, the event-averaged contact angle provides direct evidence of tangential alignment. For pushers, the median contact angle lies between approximately 20° and 30°, where 0° corresponds to perfect tangential alignment. Pullers, in contrast, display larger median contact angles, typically between 40° and 50°, indicating a less tangential orientation with respect to obstacle surfaces.

Acknowledgements

C.K. and M.R. acknowledge discussions with Sagnik Garai and Axel Voigt for theoretical insights and Simon Praetorious for the code development. S.A. acknowledges support from DFG (grant 568855070). Open Access funding provided by the Max Planck Society.

References

1. R. Tecon and D. Or, FEMS Microbiol. Rev., 2017, 41, 599–623.
2. S. Suarez and A. A. Pacey, Hum. Reprod. Update, 2006, 12, 23–37.
3. J. S. Guasto, R. Rusconi and R. Stocker, Annu. Rev. Fluid Mech., 2012, 44, 373–400.
4. Y. F. Dufrene and A. Persat, Nat. Rev. Microbiol., 2020, 18, 227–240.
5. J.-B. Raina, V. Fernandez, B. Lambert, R. Stocker and J. R. Seymour, Nat. Rev. Microbiol., 2019, 17, 284–294.
6. P. Erkoc, I. C. Yasa, H. Ceylan, O. Yasa, Y. Alapan and M. Sitti, Adv. Ther., 2019, 2, 1800064.
7. M. Sedighi, A. Zahedi Bialvaei, M. R. Hamblin, E. Ohadi, A. Asadi, M. Halajzadeh, V. Lohrasbi, N. Mohammadzadeh, T. Amiriani and M. Krutova, et al., Cancer Med., 2019, 8, 3167.
8. W. Gao and J. Wang, ACS Nano, 2014, 8, 3170–3180.
9. J. S. T. Adadevoh, S. Triolo, C. A. Ramsburg and R. M. Ford, Environ. Sci. Technol., 2016, 50, 181.
10. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.
11. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
12. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
13. S. E. Spagnolie and P. T. Underhill, Annu. Rev. Condens. Matter Phys., 2023, 14, 381–415.
14. T. Bhattacharjee and S. S. Datta, Nat. Commun., 2019, 10, 2075.
15. R. Sinaasappel, M. Fazelzadeh, T. Hooijschuur, Q. Di, S. Jabbari-Farouji and A. Deblais, Phys. Rev. Lett., 2025, 134, 128303.
16. H. Wu, B. Greydanus and D. K. Schwartz, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2101807118.
17. A. Dehkharghani, N. Waisbord and J. S. Guasto, Commun. Phys., 2023, 6, 18.
18. F. Cichos, K. Gustavsson, B. Mehlig and G. Volpe, Nat. Mach. Intell., 2020, 2, 94–103.
19. A. Doostmohammadi, M. G. Mazza, T. N. Shendruk and H. Stark, Front. Phys., 2023, 11, 1231716.
20. L. Yang, J. Jiang, F. Ji, Y. Li, K.-L. Yung, A. Ferreira and L. Zhang, Nat. Mach. Intell., 2024, 6, 605–618.
21. V. A. Baulin, A. Giacometti, D. A. Fedosov, S. Ebbens, N. R. Varela-Rosales, N. Feliu, M. Chowdhury, M. Hu, R. Fuchslin, M. Dijkstra, M. Mussel, R. V. Roij, D. Xie, V. Tzanov, M. Zu, S. Hidalgo-Caballero, Y. Yuan, L. Cocconi, C.-M. Ghim, C. Cottin-Bizonne, M. Carmen Miguel, M. Jose Esplandiu, J. Simmchen, W. J. Parak, M. Werner, G. Gompper and M. M. Hanczyc, Soft Matter, 2025, 21, 4129–4145.
22. H. Lowen and B. Liebchen, in Towards Intelligent Active Particles, ed. M. te Vrugt, Springer Nature Switzerland, Cham, 2026, pp. 257–271.
23. O. Chepizhko and F. Peruani, Phys. Rev. Lett., 2013, 111, 160604.
24. T. Bertrand, Y. Zhao, O. Bénichou, J. Tailleur and R. Voituriez, Phys. Rev. Lett., 2018, 120, 198103.
25. N. A. Licata, B. Mohari, C. Fuqua and S. Setayeshgar, Biophys. J., 2016, 110, 247–257.
26. G. Volpe and G. Volpe, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 11350–11355.
27. C. Kurzthaler, S. Mandal, T. Bhattacharjee, G. Frangipane, N. Klongvessa, H. Lowen and G. Volpe, Nat. Commun., 2021, 12, 7088.
28. C. Reichhardt and C. J. Olson Reichhardt, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2014, 90, 012701.
29. M. Zeitz, K. Wolff and H. Stark, Eur. Phys. J. E, 2017, 40, 23.
30. H. H. Mattingly, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2407313122.
31. T. Pietrangeli, R. Foffi, R. Stocker, C. Ybert, C. Cottin-Bizonne and F. Detcheverry, Phys. Rev. Lett., 2025, 134, 188303.
32. E. Lauga, W. R. DiLuzio, G. M. Whitesides and H. A. Stone, Biophys. J., 2006, 90, 400–412.
33. A. P. Berke, L. Turner, H. C. Berg and E. Lauga, Phys. Rev. Lett., 2008, 101, 038102.
34. S. E. Spagnolie and E. Lauga, J. Fluid Mech., 2012, 700, 105–147.
35. M. Contino, E. Lushi, I. Tuval, V. Kantsler and M. Polin, Phys. Rev. Lett., 2015, 115, 258102.
36. V. Kantsler, J. Dunkel, M. Blayney and R. E. Goldstein, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 1187–1192.
37. J. S. Lintuvuori, A. T. Brown, K. Stratford and D. Marenduzzo, Soft Matter, 2016, 12, 7959–7968.
38. M. Kuron, P. Stark, C. Holm and J. de Graaf, Soft Matter, 2019, 15, 5908–5920.
39. G. Junot, Phys. Rev. Lett., 2022, 128, 248101.
40. D. Takagi, J. Palacci, A. B. Braunschweig, M. J. Shelley and J. Zhang, Soft Matter, 2014, 10, 1784.
41. S. E. Spagnolie, G. R. Moreno-Flores, D. Bartolo and E. Lauga, Soft Matter, 2015, 11, 3396–3411.
42. C. Kurzthaler and H. A. Stone, Soft Matter, 2021, 17, 3322–3332.
43. S. Garai, U. Makanga, A. Varma and C. Kurzthaler, Phys. Rev. Res., 2025, 7, L042005.
44. M. J. Lighthill, Commun. Pure Appl. Math., 1952, 5, 109–118.
45. J. R. Blake, J. Fluid Mech., 1971, 46, 199–208.
46. J. R. Blake, Bull. Austr. Math. Soc., 1971, 5, 255–264.
47. P. Ahana and S. P. Thampi, Fluid Dyn. Res., 2019, 51, 065504.
48. C. Kvs and S. P. Thampi, Phys. Rev. Fluids, 2021, 6, 083101.
49. G. Li and A. M. Ardekani, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2011, 83, 046309.
50. D. Nie, Y. Ying, G. Guan, J. Lin and Z. Ouyang, J. Fluid Mech., 2023, 960, A31.
51. T. Ishikawa, M. P. Simmonds and T. J. Pedley, Phys. Rev. Lett., 2004, 93, 124501.
52. K. Ishimoto and E. A. Gaffney, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062702.
53. K. Ishimoto and D. G. Crowdy, J. Fluid Mech., 2017, 821, 647–667.
54. A. T. Brown, I. D. Vladescu, A. Dawson, T. Vissers, J. Schwarz-Linek, J. S. Lintuvuori and W. C. K. Poon, Soft Matter, 2016, 12, 131–140.
55. A. Chamolly, T. Ishikawa and E. Lauga, New J. Phys., 2017, 19, 115001.
56. K. Ishimoto, E. A. Gaffney and D. J. Smith, Front. Cell Dev. Biol., 2023, 11, 1123446.
57. M. Ramprasad, S. Mandal and P. Sinha Mahapatra, Phys. Rev. E, 2025, 111, 065402.
58. J. J. Molina, Y. Nakayama and R. Yamamoto, Soft Matter, 2013, 9, 4923–4936.
59. S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, New York, 2002, vol. 16.
60. M. Alkämper, A. Dedner, R. Klöfkorn and M. Nolte, Arch. Num. Softw., 2016, 4, 1–28.
61. S. Vey and A. Voigt, Comput. Visualiz. Sci., 2007, 10, 57–67.
62. T. Witkowski, S. Ling, S. Praetorius and A. Voigt, Adv. Comput. Math., 2015, 41, 1145–1177.
63. C. Engwer, C. Graser, S. Muthing, S. Praetorius and O. Sander, Concepts for Composing Finite Element Function Space Bases, 2025, https://arxiv.org/abs/2508.10125.
64. P. Bastian, M. Blatt, A. Dedner, N.-A. Dreier, C. Engwer, R. Fritze, C. Gräser, C. Grüninger, D. Kempf, R. Klöfkorn, M. Ohlberger and O. Sander, Comput. Math. Appl., 2021, 81, 75–112.
65. Adaptive Multi-Dimensional Simulations, https://gitlab.com/amdis/amdis, v2.11.
66. D. Gao, Z. Wang, M. Jain, A. J. T. M. Mathijssen and R. Tao, Selective trapping of bacteria in porous media by cell length, 2025.
67. J. D. Wheeler, E. Secchi, R. Rusconi and R. Stocker, Annu. Rev. Cell Dev. Biol., 2019, 35, 213–237.
68. R. Tao, A. Thery, S. Que and A. J. T. M. Mathijssen, Newton, 2026, 2, 100337.
69. G. Jing, A. Zottl, É. Clément and A. Lindner, Sci. Adv., 2020, 6, eabb2012.
70. A. J. T. M. Mathijssen, N. Figueroa-Morales, G. Junot, É. Clément, A. Lindner and A. Zottl, Nat. Commun., 2019, 10, 3434.
71. Z. Peng and J. F. Brady, Phys. Rev. Fluids, 2020, 5, 073102.
72. H. L. Gertack, P. A. Hampshire, C. Wohlgemuth, R. Alert and S. Aland, Soft Matter, 2026, 22, 907–925.
73. N. Waisbord, A. Dehkharghani and J. S. Guasto, Nat. Commun., 2021, 12, 5949.
74. G. L. Mino, M. Baabour, R. Chertcoff, G. Gutkind, E. Clement, H. Auradou and I. Ippolito, Adv. Microbiol., 2018, 08, 451–464.
75. E. Secchi, A. Vitale, G. L. Mino, V. Kantsler, L. Eberl, R. Rusconi and R. Stocker, Nat. Commun., 2020, 11, 2851.
76. M. Residori, S. Mandal, A. Voigt and C. Kurzthaler, Phys. Rev. Res., 2025, 7, L012032.
77. A. Creppy, E. Clement, C. Douarche, M. V. D'Angelo and H. Auradou, Phys. Rev. Fluids, 2019, 4, 013102.
78. P. Das, M. Residori, A. Voigt, S. Mandal and C. Kurzthaler, Vorticity-induced surfing and trapping in porous media, 2025, https://arxiv.org/abs/2511.02471.
79. R. Alonso-Matilla, B. Chakrabarti and D. Saintillan, Phys. Rev. Fluids, 2019, 4, 043101.
80. M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189.
81. A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 118101.
82. D. Saintillan and M. J. Shelley, C. R. Phys., 2013, 14, 497–517.
83. K. Drescher, Y. Shen, B. L. Bassler and H. A. Stone, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 4345–4350.
84. D. L. Kurz, E. Secchi, F. J. Carrillo, I. C. Bourg, R. Stocker and J. Jimenez-Martinez, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2122202119.
85. N. Aguillon, A. Decoene, B. Fabrèges, B. Maury and B. Semin, ESAIM: Proc., 2012, 38, 36–53.
86. X. Li, J. Lowengrub, A. Rätz and A. Voigt, Commun. Math. Sci., 2009, 7, 81–107.
87. S. Aland, J. Lowengrub and A. Voigt, Comput. Model. Eng. Sci., 2010, 57, 77–106.
88. W. Marth, S. Aland and A. Voigt, J. Fluid Mech., 2016, 790, 389–406.
89. M. Reder, P. W. Hoffrogge, D. Schneider and B. Nestler, Int. J. Numer. Meth. Eng., 2022, 123, 3757–3780.
90. D. Mokbel, H. Abels and S. Aland, J. Comput. Phys., 2018, 372, 823–840.
91. G. Caginalp and P. C. Fife, SIAM J. Appl. Math., 1988, 48, 506–518.
92. R. L. Pego, Proc. R. Soc. London, Ser. A, 1989, 422, 261–278.
93. S. Balay, W. D. Gropp, L. C. McInnes and B. F. Smith, Modern Software Tools in Scientific Computing, 1997, pp. 163–202.
94. S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang and J. Zhang, PETSc Web page, 2025, https://petsc.org/.
95. P. Amestoy, I. S. Duff, J. Koster and J.-Y. L'Excellent, SIAM J. Matrix Anal. Appl., 2001, 23, 15–41.
96. P. Amestoy, A. Buttari, J.-Y. L'Excellent and T. Mary, ACM Trans. Math. Softw., 2019, 45, 1–26.

---