# Modeling complex particle suspensions: Perspectives on the rigid multiblob method

Blaise Delmotte, Florencio Balboa Usabiaga

► To cite this version:

Blaise Delmotte, Florencio Balboa Usabiaga. Modeling complex particle suspensions: Perspectives on the rigid multiblob method. Physical Review Fluids, 2025, 10, pp.100701. 10.1103/64b1-lfmc. hal-05347137

HAL Id: hal-05347137

https://hal.science/hal-05347137v1

Submitted on 17 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling complex particle suspensions: perspectives on the rigid multiblob method

Blaise Delmotte¹, * and Florencio Balboa Usabiaga², †

¹LadHyX, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

²BCAM - Basque Center for Applied Mathematics,
Mazarredo 14, Bilbao, E48009, Basque Country - Spain

Many suspensions contain particles with complex shapes that are affected not only by hydrodynamics, but also by thermal fluctuations, internal activity, kinematic constraints and other long-range non-hydrodynamic interactions. Modeling these systems represents a significant numerical challenge due to the interplay between different effects and the need to accurately capture multiscale phenomena. In this article we review recent developments to model large suspensions of particles of arbitrary shapes and multiple couplings with controllable accuracy within the rigid multiblob framework. We discuss the governing equations, highlight key numerical developments, and illustrate applications ranging from microswimmers to complex colloidal suspensions. This review illustrates the effectiveness and versatility of the rigid multiblob method in tackling a wide range of physical problems in fluid mechanics, soft matter physics, biophysics, materials and colloidal science.

## CONTENTS

I. Introduction 2

II. Governing equations and modeling challenges 3

A. Governing equations for the hydrodynamic problem 3

1. The mobility problem 4

B. Accounting for thermal fluctuations 4

C. Constrained kinematics 6

D. Additional couplings: example with reactive particles 7

III. The rigid multiblob framework 8

A. Rigid multiblob model for the Stokes problem 8

B. Including thermal fluctuations 10

1. Stochastic particle velocities 10

2. Drift term 11

C. How to scale big and compute fast? 11

1. Fast hydrodynamics with Green's functions 12

2. Fast hydrodynamics with grid-based methods 12

D. Handling kinematic constraints in large suspensions 13

E. The multi-physics multiblob method 14

IV. Applications 15

1. Active particles 15

2. Complex colloidal suspensions 18

V. Perspectives 18

Acknowledgments 20

References 20

<footnote>* blaise.delmotte@cnrs.fr</footnote>

<footnote>† fbalboa@bcamath.org</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Illustration of various effects at play in complex suspensions. Each box represents one of these effects and shows experimental illustrations. a) Patchy colloids used to form gels or suspensions with complex structures [1]. Scale bars: 1µm. b) Fumed silica particle used to tune the rheology of colloidal suspensions [2]. Scale bar: 150nm. c) Diatom chain *Bacillaria Paxillifer* made of sliding rod-like cells [3]. Scale bar: 100µm. d) Bacterium *Vibrio Cholerae* swimming with its flagellum [4]. e) Lock and key colloids [5]. Scale bars: 2µm. f) Chains of freely-jointed droplets [6]. g) Self assembly of catalytic rods [7]. h) Spontaneous rotation of chiral catalytic particles [8]. Scale bar: 20µm.</figcaption></figure>

# I. INTRODUCTION

Suspensions of small, micron-sized, particles are ubiquitous in biological, environmental and industrial systems. Their typical size and speed are such that, at their scale, the Reynolds number is very small (Re ≪ 1) and the surrounding flow follows the Stokes equations. In this limit, the hydrodynamic interactions between the particles are long-ranged, meaning that all particle motions are instantaneously coupled through the fluid.

In many applications, suspended particles have complex shapes and are subject to additional effects than hydrodynamic interactions, such as thermal fluctuations (i.e. Brownian motion), internal activity, kinematic constraints, or non-hydrodynamic long-ranged interactions (e.g. chemical, magnetic, electrostatic...), see Fig. 1. These effects, in combination with the particle shape, lead to complex dynamics at the scale of the particles, such as non-trivial diffusion, self-assembly, rotation, or propulsion, and provide intricate macroscopic properties at the scale of the suspensions, such as complex rheological and elastic responses, collective diffusion, collective motion, flow generation and mixing.

The interplay between the individual behavior of complex particles and their suspension dynamics has attracted significant attention in the past decades. Experimental methods for synthesizing complex suspensions and numerical tools for modeling them have developed rapidly. However, significant hurdles remain in the numerical domain. Indeed, the modeling of complex suspensions poses multiple challenges due to particle geometry, thermal fluctuations with multiplicative noise, nontrivial boundary conditions, kinematic constraints, active deformations or multiphysics couplings (see Fig. 1). In addition, these effects must be accurately captured to bridge the gap between the particle scale and the macroscopic suspension scale. Traditional numerical methods, such as Boundary Element Methods, are accurate and can now handle large collections of spherical [9], or even spheroidal particles [10], but they remain restricted to few particles of complex shape due to their high computational cost and cannot handle thermal fluctuations (see Ref. [11] for an exception in 2D).

An alternative approach is the rigid multiblob framework. In the rigid multiblob (RMB) method, also called composite bead model [12], or Stokesian Dynamics of rigid assemblies [13], rigid bodies are discretized using a collection of minimally resolved spherical markers, also called *blobs*, or beads, which are constrained to move as a rigid body. This idea first emerged nearly 50 years ago [14] to compute the friction coefficients of non-spherical aggregates. Unlike point-forces, i.e. Stokeslets, which are singular, these blobs have a finite size which allows to compute their hydrodynamic interactions, and thus the hydrodynamic coefficients of the aggregates, with regularized Green's functions, i.e. mobility matrices, which were first derived by Rotne & Prager [15] and Yamakawa [16] in the early 70s. This seminal multiblob idea was further formalized and developed in the late 1970s and early 1980s [17–19], see [20] for a review of these
pioneering works. Nowadays, the multiblob method is still widely used to compute hydrodynamic coefficients (mobility, diffusion and resistance matrices) of rigid bodies such as colloids, molecules or clusters [12, 21–37]. However, since the 1990s, and more particularly in the past decade, several groups across the world have extensively developed the method to extend its domains of application beyond the motion of passive colloidal clusters, but also to improve its numerical performances, using efficient summation techniques, modern linear algebra algorithms and grid-based numerical methods [12, 38–48].

While it shares some similarities with the well-known method of Regularized Stokeslets [35, 49], in which bodies are discretized with regularized point forces, the rigid multiblob method additionally handles thermal fluctuations and effectively couples with grid-based methods (such as Immersed Boundary methods [50, 51]) to solve hydrodynamic interactions between blobs, whereas the former can only be used for deterministic problems and relies exclusively on analytical Green’s functions.

In this Perspective, intended for a wide audience of fluid dynamicists and soft matter physicists, we present recent developments of the rigid multiblob method to model suspensions containing a large number of particles of arbitrary shapes and subject to multiple physical couplings, with controllable accuracy and high scalability. These recent efforts open new perspectives for the modeling of complex particle suspensions and their intricate phenomena. We first introduce the governing equations of complex particle suspensions in Section II, with an emphasis on the modeling challenges that they represent. Then we introduce the rigid multiblob framework in Section III and present recent developments that successfully address these modeling challenges. Finally we showcase some examples and applications where the method has been used in Section IV and outline future directions regarding the development of the method and its use to address new physical phenomena in Section V.

The framework presented in this article is implemented in a code which is free software, user friendly for non-experts and continuously enriched by a group of collaborators in a GitHub repository (https://github.com/stochasticHydroTools/RigidMultiblobsWall). The repository contains many hands-on examples that users can run on a basic laptop. This code has been actively used in the community, in combination with experiments and theory, to study a variety of systems involving (micro-)swimmers in various geometries [44, 47, 52, 53], passive, active and reactive colloids [46, 54–60], diffusing molecules [61, 62], suspensions of sedimenting cubes [63], viruses [64], functionalized nanoparticles [65], Möbius strips [66], branched particles suspensions [45, 67], fibers in structured environments [68, 69], coated drops [70], magnetic spheres [71], chains [72] or sheets [73].

## II. GOVERNING EQUATIONS AND MODELING CHALLENGES

This section introduces the equations that govern the motion of complex suspensions in viscous flows. Complexity is introduced at multiple levels of the system: in the shape of the particles, in their surface velocity and deformation, in the fluid through thermal fluctuations, in the relative motion between particles due to kinematic constraints, or with non-hydrodynamic couplings such as chemical interactions. We first introduce the governing equations of the hydrodynamic problem and then extend them to include thermal fluctuations. Then we explain how the equations of motion are modified in the presence of kinematic constraints. Finally, we tackle the situation where particles interact through additional physical couplings and provide the governing equations for the case of chemical and phoretic interactions between the suspended particles. At each step, we briefly discuss the main modeling challenges that will be addressed in Section III.

### A. Governing equations for the hydrodynamic problem

Let {$\mathcal{B}_p$}$_{p=1}^M$ be a set of $M$ bodies immersed in a Stokes flow. The configuration of each body $p$ is described by the location of a tracking point, $q_p$, and its orientation represented by the unit quaternion $\theta_p$, or in compact notation $x_p = \{q_p, \theta_p\}$. The linear and angular velocities of the tracking point are $u_p$ and $\omega_p$. The external force and torque applied to a body are $f_p$ and $\tau_p$. We will use the concise notation $U_p = \{u_p, \omega_p\}$ and $F_p = \{f_p, \tau_p\}$ as well, while unscripted vectors will refer to the composite vector formed by the variables of all the bodies, e.g. $U = \{U_p\}_{p=1}^M$. The velocity and pressure, $v$ and $p$, of the fluid with viscosity $\eta$ are governed by the Stokes equations [74]

$$-\nabla p + \eta \nabla^2 v = 0, \qquad (1)$$

$$\nabla \cdot v = 0, \qquad (2)$$

while at the bodies surface the fluid obeys the boundary condition

$$v(r) = (\mathcal{K}U_p)(r) + u_s(r) = u_p + \omega_p \times (r - q_p) + u_s(r) \text{ for } r \in \partial \mathcal{B}_p, \qquad (3)$$
where we have introduced the geometric operator $\mathcal{K}$ of rigid body motion. $u_s$ is the surface velocity of the particle relative to its rigid body motion. This velocity can correspond to the displacement of the surface due to active deformations, also called swimming gait [39], it can also be interpreted as the velocity of an ambient flow (e.g. shear, vortical or quadratic flows), or as an active slip velocity due to phoresis [75, 76] or ciliary motion [77] evaluated at the surface of a rigid particle. Boundary conditions also apply at the bounding surface of the domain considered.

Since inertia does not play a role in Stokes flows the conservation of linear and angular momentum reduces to the balance of force and torque. For every body $p$ the balance between hydrodynamic and external stresses is given by

$$ \int_{\partial \mathcal{B}_p} \lambda(\mathbf{r}) \, dS_r = \mathbf{f}_p, \qquad (4) $$

$$ \int_{\partial \mathcal{B}_p} (\mathbf{r} - \mathbf{q}_p) \times \lambda(\mathbf{r}) \, dS_r = \tau_p, \qquad (5) $$

where $\lambda(\mathbf{r})$ is a constraint force per unit area that enforces the boundary condition (3). In absence of surface velocity ($u_s = 0$), $-\lambda$ is the hydrodynamic traction exerted on the bodies by the fluid, i.e. $-\lambda(\mathbf{r}) = \boldsymbol{\sigma} \cdot \mathbf{n}$ where $\boldsymbol{\sigma} = -p\mathbf{I} + 2\eta(\nabla\mathbf{v} + \nabla\mathbf{v}^T)$ is the fluid stress tensor [78, 79]. The adjoint of the geometric operator $\mathcal{K}$ can be used to write the balance of force and torque for all bodies as [80]

$$ \mathcal{K}^T \lambda = \mathbf{F}. \qquad (6) $$

In order to compute the trajectories of the immersed bodies, the equations of motion are integrated in time

$$ \frac{d\mathbf{x}}{dt} = \mathbf{U}. \qquad (7) $$

The integration of unit quaternions requires some care to correctly represent rotations, details are discussed in [40, 45].

If the particle surface moves relative to the rigid body motion, then the surface shape can be integrated in time using (3).

## 1. The mobility problem

Given the nonhydrodynamic forces and torques acting on the rigid bodies, $\mathbf{F}$, and their surface velocity $u_s$, the equations (1)-(5) can be solved for the bodies velocities $\mathbf{U}$ and the constraint forces $\lambda$. Since the Stokes equations are linear we can write the velocity of $\mathbf{M}$ rigid bodies as

$$ \mathbf{U} = \mathbf{N}\mathbf{F} - \widetilde{\mathbf{N}}u_s, \qquad (8) $$

where $\mathbf{N} = \mathbf{N}(\mathbf{x})$ is the $6M \times 6M$ body mobility matrix that couple the forces and torques acting on the rigid bodies with their velocities, and $\widetilde{\mathbf{N}}$ is a linear operator ($\widetilde{\mathbf{N}} : L^2(\partial \mathcal{B}, \mathbb{R}^3) \to \mathbb{R}^{6M}$)[81], where $\partial \mathcal{B}$ is the particle surface, which we will call "surface mobility operator", mapping the surface motion to the particle velocities. As shown by Swan et al. [39], which generalizes the work of Stone and Samuel [82], $\widetilde{\mathbf{N}}$ can be obtained for a single particle using the Lorentz reciprocal theorem.

### Modeling challenges:

In order to model complex suspensions, it is crucial to develop methods that can capture short and long-ranged hydrodynamic interactions (HI) between tens or hundreds of thousands of particles.

However, complex particle shapes, non-trivial surface boundary conditions, the slow decay of HI in Stokes flow and their singular nature in the near-field make this task challenging.

Traditional methods like the Boundary Element Methods [78] have been widely used to solve the many-body mobility problem in particle suspensions. However, they suffer from convergence issues at small interparticle distances, require special, and complex, techniques to deal with the singularities distributed on the particle surfaces and are therefore limited to very few particles in dynamic simulations.

## B. Accounting for thermal fluctuations

If the suspended particles are small enough, they will also move randomly due to thermal fluctuations in the surrounding fluid. In order to satisfy the correct diffusion in the suspension, the statistics of the resulting particle
fluctuating velocities, $U_{\text{th}}$, must satisfy the fluctuation-dissipation theorem [83, 84]

$$ \langle U_{\text{th}} \rangle = 0, \qquad (9) $$

$$ \langle U_{\text{th}}(t)U_{\text{th}}(t')^T \rangle = D\delta(t-t') = 2k_B T N \delta(t-t'), \qquad (10) $$

where the operator $\langle \cdot \rangle$ is the ensemble average, $D$ the short-time diffusion tensor of the particles in the suspension, and $k_B T$ the thermal energy. Eq. (10) is a generalization of the famous Stokes-Einstein relation for a single particle to suspensions with multiple particles, where the $6M \times 6M$ short-time diffusion tensor of the particles $D$ is proportional to the mobility matrix. The fluctuations have the same origin as the viscous dissipation one must do work against when perturbing the system [85], and the fluctuation-dissipation theorem predicts that the response of a system in thermodynamic equilibrium to a small applied force is the same as its response to a spontaneous fluctuation [86].

In order to incorporate the effects of Brownian motion in the fluid, a stochastic stress tensor $Z$ [87, 88] is added into the right hand side of the Stokes equations

$$ -\nabla p + \eta \nabla^2 v = \nabla \cdot Z, \qquad (11) $$

$$ \nabla \cdot v = 0, \qquad (12) $$

The tensor $Z$ encodes the stress generated by the random motion and collisions of the fluid molecules, which, at scales much larger than the molecular scale, can be modeled as delta correlated in space and time

$$ \langle Z_{ij}(x, t) \rangle = 0, \qquad (13) $$

$$ \langle Z_{ij}(x, t) Z_{kl}(x', t') \rangle = 2\eta k_B T (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \delta(x - x')\delta(t - t'), \qquad (14) $$

The ensemble average of $Z$ is zero because $Z$ represents the random fluctuations around the mean deterministic flow [85, 87]. The magnitude of the correlations, proportional to the thermal energy $k_B T$, and their symmetries ensure that the fluctuation-dissipation theorem (10) is satisfied by Eqs. (11)-(14) [87, 88].

The mobility problem formed by (11)-(12) together with the boundary conditions (3) and force-torque balance (4)-(5) leads to the (Itô) overdamped Langevin equations for the particle positions and orientations [89–91]

$$ U = NF - \tilde{N}u_s + \sqrt{2k_B T N^{1/2} Z} + k_B T \partial_x \cdot N, \qquad (15) $$

where two additional terms appear on the right hand side.

The first is the random particle velocity increments due thermal agitation given by $U_{\text{th}} = \sqrt{2k_B T N^{1/2} Z}$, with $N^{1/2}$ being the square root of the mobility matrix (i.e. $N^{1/2} (N^{1/2})^T = N$), and $Z$ a vector of independent white noise processes. The noise dependence on the mobility matrix ensures that the fluctuation-dissipation theorem (10) is satisfied, a necessary condition for the Boltzmann distribution to be recovered at equilibrium [92]. Note that while both $Z$ and $N$ present in Eqs. (11) and (15) represent white noise terms they are not equivalent because the first is a continuous field defined on the fluid while the second is a discrete vector on the particles.

Along with the random velocities, the overdamped equations of motion require the inclusion of the thermal drift term, $k_B T \partial_x \cdot N$, which is the divergence of the mobility matrix with respect to the particle positions and orientations [93]. This ensures that the stochastic differential equation with an Itô interpretation of the stochastic integral yields dynamics consistent with those of Smoluchowski's equation for the corresponding probability distribution.

Using uncontrolled approximations in the random particle velocities and drift terms results in incorrect equilibrium particle distributions and distorted dynamics which cannot always be predicted a priori. For example, simple approximations in the random particle velocities can lead to errors of about 15% in the diffusion coefficient of polymers and up to 50% for the diffusion of particles in dense suspensions [94] (these errors were measured a posteriori because it was not possible to predict their magnitude beforehand). Neglecting the drift term leads to ~40% errors in the probability distribution function a single colloid diffusing inside a channel [92, 95, 96] and result in incorrect dynamics of active particles near walls [97].
### Modeling challenges:

Computing the square root of the mobility matrix $N^{1/2}$ with classical linear algebra, such as Cholesky decomposition, is costly as it scales badly with the particle number $O(M^3)$. Though several methods have been introduced to accelerate the matrix square root computation [94, 98–100], the inclusion of thermal fluctuations has often limited simulations to having very small particle numbers or ignoring completely hydrodynamic interactions between the particles.

The additional drift term $k_B T \partial_x \cdot N$ requires computing the divergence of the mobility matrix with respect to all the particle positions and orientations, which incurs a high computational cost and cannot be done analytically for complex particle shapes.

In the past decade, significant advances have been made to compute these terms at a nearly linear computational cost ($O(M)$) in suspensions of spherical particles [92, 94, 95, 99, 101], but their generalization to nonspherical shapes has been lacking, and traditional methods for more complex shapes, such as the Boundary Element Method [78] or the Method of Regularized Stokeslet [49], cannot handle thermal fluctuations.

## C. Constrained kinematics

In a variety of scenarios, the movement of certain bodies is either wholly or partially constrained, while others are free to move. This phenomenon can be observed in systems comprising freely-flowing particles and immobile obstacles, such as suspensions flowing through porous media or in microfluidic chips structured by pillar arrays. In other cases, which encompass a diverse array of microscopic systems, small objects are locked by geometric constraints or connected by rigid bonds with exceedingly short relaxation timescales. Stiff bonds can be found in polymer chains, freely-jointed molecules [5, 6] or in molecules with bond angles. They can also be found in the hinges of multilink artificial microswimmers [102–104], in the hook connecting the flagellum to the cell body of bacteria [105, 106], or in the form of adhesive forces between sliding cells in diatom colonies [107–109]. The use of kinematic constraints allows for the elimination of these fast degrees of freedom, thereby enabling larger time steps in numerical simulations.

In the majority of cases, the kinematic constraints mentioned above are holonomic, meaning that they depend only on the particle positions, orientations (and time $t$) and can be written in the form

$$g(x, t) = 0 \qquad (16)$$

where we recall that $x = \{q, \theta\}$ is a vector that collects all the particle positions and orientations. Typical holonomic, scalar or vectorial, constraints are for instance position constraints

$$g(x, t) = q_p - q_0(t) = 0, \qquad (17)$$

where $q_0(t)$ is the prescribed position of a body $p$; prescribed distance $d(t)$ between two bodies

$$g(x, t) = \frac{1}{2} \|q_p - q_q\|_2^2 - d(t)^2 = 0, \qquad (18)$$

or ball-and-socket joints between two bodies

$$g(x) = q_p + R(\theta_p) \Delta l_p(t) - q_q - R(\theta_q) \Delta l_q(t) = 0, \qquad (19)$$

where $\Delta l_p(t)$ is the vector from the body $p$ to the hinge of the joint in the body frame of reference. This vector is then rotated to the laboratory frame of reference by the rotation matrix $R(\theta_p)$. A constant vector $\Delta l_p = c$ represents a passive link, while $\Delta l_p = \Delta l_p(t)$ represents a moving link.

Following the classical framework of Lagrangian mechanics, a Lagrange multiplier is associated to each kinematic constraint in the equations of motion. In the context of immersed bodies at low Reynolds number, the inclusion of Lagrange multipliers results in the introduction of constraint forces $F^C$ in the mobility problem (8)

$$U = N(F + F^C) - \tilde{N}u_s. \qquad (20)$$

In order to model suspensions with kinematic constraints, one has to solve the augmented mobility problem (20) together with the equations of motion (7) subject to a diversity of kinematic constraints (e.g. (17)-(19)) acting on a portion, or the totality, of the degrees of freedom in the system. Due to the nonlinearity of the kinematic constraints, this yields a nonlinear system of equations for the positions, orientations, as well as the constraint forces associated with each constraint.
**Modeling challenges:**

Nonlinearities pose a challenge to integrate the equations of motion (7) for large suspensions [110, 111]. They require using nonlinear solvers, usually relying on Newton's method, where one needs to compute a large and complex inverse Jacobian matrix at each iteration. Until now, this challenging task has limited numerical studies to small numbers of particles [112–114] or led to methods restricted to a specific type of constraint [111].

## D. Additional couplings: example with reactive particles

In addition to hydrodynamic interactions, suspended particles can interact through other interactions such as chemical, electrostatic or magnetic couplings. Just like the Stokes equations, the equations governing these fields are elliptic, which means that they can be solved with similar methods. Below we outline the example of reactive phoretic particle.

The individual and collective dynamics of self-propelled reactive particles have attracted significant attention in recent decades [115–119]. Chemically-active phoretic colloids catalyze surface reactions to modify the concentration of chemical solutes surrounding them in order to self-propel, a phenomenon known as self-diffusiophoresis. In doing so, they generate long-ranged hydrodynamic flows and chemical gradients that modify the trajectories of other particles. As a result, the dynamics of reactive suspensions is fundamentally governed by hydro-chemical interactions.

Assuming that the magnitude of the phoretic flows is small and the solute diffusivity $D$ is large, the Péclet number characterizing the motion of the solute particles is small enough to neglect solute advection. Thus, in the fluid domain the concentration of solute $c = c(r)$ diffuses according to the Laplace equation

$$D \nabla^2 c = 0. \qquad (21)$$

The immersed colloids activate chemical reactions on their surfaces to produce or consume solute. They are introduced as boundary conditions of the Laplace equation. The boundary conditions are imposed on the concentration fluxes, i.e. Robin boundary conditions, and to simplify the mathematical problem we only consider linear boundary conditions. In particular, we consider sinks that consume solute at a rate $k$ proportional to the local concentration, sources that produce solute at a constant flux or a linear combination of both [76, 120, 121]

$$D \frac{\partial c}{\partial n} = D n \cdot \nabla c = kc - \alpha \quad \text{on } \partial B_p, \qquad (22)$$

where $n = n(r)$ is the surface normal on the colloids, $k = k(r) > 0$ the surface reaction rate and $\alpha = \alpha(r) > 0$ the surface production flux, both of which can vary over the colloidal surface. Since the reaction flux is proportional to the local concentration while the concentration diffuses with a finite diffusion coefficient $D$, this boundary condition allows to model either diffusive-limited or reaction-limited systems. The balance between diffusion and reaction is characterized by the Damkohler number, $Da = kR/D$, which represents the ratio between the diffusive time over one colloidal radius, $R^2/D$, and the reaction time scale $R/k$ [121, 122]. Additional boundary conditions at the domain boundary, such as vanishing concentration at infinity, or zero flux at the domain walls, must be included to close the system.

Due to diffusiophoresis the concentration gradients near the colloid surface induce a slip velocity [75, 123]

$$u_s = \mu \nabla_{\parallel} c = \mu (I - nn^T) \nabla c \quad \text{on } \partial B_p, \qquad (23)$$

where the coefficient of proportionality $\mu = \mu(r)$ is the surface mobility and $\nabla_{\parallel}$ the surface gradient, i.e. the gradient projected to the surface tangent plane.

Modeling reactive suspensions requires solving first the Laplace problem (21), with Robin boundary conditions on the particle surfaces (22), in order to determine the concentration and then its gradients on the particle surfaces $\nabla c$, which are used to compute the active slip (23) that appears in the boundary conditions of the (fluctuating) Stokes mobility problem (15).

**Modeling challenges:**

The Laplace problem bearsome similarities with the Stokes problem: interactions are long-ranged and boundary conditions are non-trivial on the particles' surface. Solving for the surface concentration gradients in large collections of (colloidal) reactive particles with arbitrary shapes and non-uniform catalytic surface properties is therefore challenging. Similarly to the Stokes problem, the preferred method for such systems is the Boundary Element, which does not handle thermal fluctuations and is restricted to very few particles in dynamic simulations.
<figure><img src="image_3.png" /><figcaption>FIG. 2. Diagram of the rigid multiblob (RMB) framework. The colored boxes refer to a methodological challenge that is addressed in a subsection of the manuscript.</figcaption></figure>

### III. THE RIGID MULTIBLOB FRAMEWORK

As explained in Section II the particle shape, hydrodynamic interactions, thermal fluctuations, constrained kinematics and additional non-hydrodynamic couplings pose multiple challenges to the modeling of complex particle suspensions.

In this section we present recent efforts to develop a unified, yet versatile, and scalable framework, based on the multiblob model, for modeling hydrodynamic and multiphysics interactions in large suspensions of arbitrarily shaped passive and active particles, taking into account thermal fluctuations and kinematic constraints.

Figure 2 shows a diagram of the rigid multiblob framework for modeling complex suspensions. Each colored box represents a bottleneck in the method which has recently undergone significant development, and the text next to it refers to the corresponding section in this article.

#### A. Rigid multiblob model for the Stokes problem

In the rigid multiblob (RMB) method rigid bodies are discretized using a collection of minimally resolved spherical nodes of hydrodynamic radius $a$, also called blobs, which are constrained to move as a rigid body. Figure 3 shows examples of various geometries discretized with the RMB.

Each blob is subject to a force, $\lambda_i$, that enforces the rigid motion of the body. Thus, the integrals in the balance of force and torque (4)-(5) become sums over the blobs

$$ \sum_{i \in \mathcal{B}_p} \lambda_i = f_p, \qquad (24) $$

$$ \sum_{i \in \mathcal{B}_p} (\boldsymbol{r}_i - \boldsymbol{q}_p) \times \boldsymbol{\lambda}_i = \boldsymbol{\tau}_p, \qquad (25) $$
<figure><img src="image_4.png" /><figcaption>FIG. 3. Examples of particles of various shapes discretized with the RMB model. (a) Nonuniform patchy sphere [124]. (b) Chemically active rod [46]. (c) Branched particle [67]. (d) Globular protein called Lysozyme [26]. (e) Möbius strip [66]. (f) HIV [64]. (e) Bacterium in a pipe [125].</figcaption></figure>

where $\mathbf{r}_i$ is the position of blob $i$ belonging to body $p$ ($\mathcal{B}_p$). Note that here $\lambda_i$ now represents finite forces and not density forces as in the continuum formulation since it includes the quadrature weights of the blobs. The boundary condition on the particle surfaces (3), as in collocation methods [78], is evaluated at each blob

$$v(\mathbf{r}_i) = \sum_j M_{ij} \lambda_j = \mathbf{u}_p + \omega_p \times (\mathbf{r}_i - \mathbf{q}_p) + \mathbf{u}_s(\mathbf{r}_i) \quad \text{for all } i \in \mathcal{B}_p, \qquad (26)$$

where $v(\mathbf{r}_i)$ is the fluid velocity at the position of blob $i$, $M$ is the $3N_b \times 3N_b$ mobility matrix that mediates the hydrodynamic interactions between blobs, and $N_b$ is the total number of blobs in the system. The term $M_{ij}$ couples the force acting on the blob $j$ to the velocity of the blob $i$. Depending on the geometry of the domain under consideration, this matrix can be written explicitly if the Green's function of the problem is known, or its action can be computed with matrix-free methods relying on grid-based Stokes solvers, such as Immersed Boundary Methods [50] or the Force Coupling Method [51] (see Section III C 2).

A well-known example of regularized Green's function is the Rotne-Prager-Yamakawa (RPY) approximation [15, 126], where the regularization is done by a double convolution of the Stokes Green's function, $\mathcal{G}(\mathbf{r}, \mathbf{r}')$, with Dirac delta functions defined on the surface of a sphere of radius $a$:

$$M_{ij} = M(\mathbf{r}_i, \mathbf{r}_j) = \frac{1}{(4\pi a^2)^2} \int \delta(|\mathbf{r}' - \mathbf{r}_i| - a) \mathcal{G}(\mathbf{r}', \mathbf{r}'') \delta(|\mathbf{r}'' - \mathbf{r}_j| - a) d^3\mathbf{r}' d^3\mathbf{r}'' \qquad (27)$$

Wajnryb et al. showed that, for non-overlapping blobs, an equivalent definition of the RPY mobility is [126]

$$M_{ij} = M(\mathbf{r}_i, \mathbf{r}_j) = \left[ \mathbf{I} + \frac{a^2}{6} \nabla_{\mathbf{r}'}^2 \right] \left[ \mathbf{I} + \frac{a^2}{6} \nabla_{\mathbf{r}''}^2 \right] \mathcal{G}(\mathbf{r}', \mathbf{r}'')|_{\substack{\mathbf{r}'=\mathbf{r}_i \\ \mathbf{r}''=\mathbf{r}_j}} \qquad (28)$$

which makes evident the role of the blob radius as a regularization parameter.

The well-known RPY mobility matrix was developed to provide a symmetric positive definite, pairwise approximation. The use of a symmetric positive definite matrix $M$ is one of the main advantages of the RMB as it allows to incorporate Brownian fluctuations efficiently as shown in Sec. III B. In contrast, the Regularized Stokeslets method uses a mobility $M$ that is not positive definite by design [35] and, current Boundary Element Methods use singular and near singular quadrature rules to compute the action of the Stokes Green's function, $\mathcal{G}(\mathbf{r}, \mathbf{r}')$ that generate non-symmetric mobility matrices $M$ and thus, cannot incorporate Brownian fluctuations efficiently (see Ref. [11] for an exception in 2D). As explained in Sec. III C 2, the mobility matrix of grid-based- methods is symmetric positive definite by design and has the same structure (28) as the RPY mobility to order $O(a^4)$.

Since $M$ is positive definite it is not necessary to employ special quadrature rules to evaluate near or self interactions. However, the formula (27) introduces a regularization error proportional to the blob radius [46]. As a typical
discretization of a rigid body surface employs $N_b \sim 1/a^2$ blobs, the convergence rate scales with the number of blobs as $\sim N_b^{-1/2}$. Despite this slow convergence it is possible to use optimized grids to obtain results accurate up to a few percent as far as blobs of different bodies do not overlap [46, 127]. We illustrate the accuracy of the RMB method in a couple of examples in Secs. IIID-III E and discuss optimized grids in Sec. V.

The original RPY mobility computes hydrodynamic interactions between spherical particles of equal radii in an unbounded domain [126]. Extensions of the RPY matrix for particles of different radii [128], in a background shear flow [126], above a fluid-fluid interface [129] or a no-slip boundary [130] with arbitrary periodicities [131, 132], or in a spherical cavity [133] are available in the literature. Section IIIC 1 enumerates fast methods to compute the action of the RPY tensor.

The whole linear system to solve the mobility problem is found by combining (24)-(26),

$$ \begin{bmatrix} M & -K \\ -K^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \\ U \end{bmatrix} = \begin{bmatrix} u_s \\ -F \end{bmatrix}. \qquad (29) $$

In the linear system the matrix $K$ is a discretization of the operator $\mathcal{K}$ introduced in Section II A: it transforms the rigid body velocities to blob velocities. The solution of (29) is

$$ U = NF - NK^T M^{-1} u_s \qquad (30) $$

where $N = (K^T M^{-1} K)^{-1}$ and $\tilde{N} = N K^T M^{-1}$ are the $6M \times 6M$ body mobility matrix and $6M \times 3N_b$ surface mobility operator given by the rigid multiblob model.

**Key contributions:**

The RMB allows to simulate suspensions of bodies of arbitrary shape, with surface velocities, without dealing with singularities, in a robust, simple and versatile framework. The RMB naturally interfaces with any fast method, either grid-based or relying on Green's functions, for computing hydrodynamic interactions between blobs (see Section III C). The linear system (29) can be solved efficiently with preconditioned iterative methods (see Section III C).

## B. Including thermal fluctuations

As explained in Section II B, two new terms appear in the equations of motion when thermal fluctuations are included: the stochastic particle velocities $U_{th} = \sqrt{2k_B T N^{1/2}} Z$ and the drift term $k_B T \partial_x N$. Both of these terms are challenging to account for with non-spherical particles and are computationally expensive. Below we present the formalisms recently developed to handle these terms correctly and efficiently.

### 1. Stochastic particle velocities

The computation of the stochastic particle velocities $U_{th} = \sqrt{2k_B T N^{1/2}} Z$ can be readily incorporated into the system by including a random velocity at the blob level, $u_{th}$, in the boundary conditions (26) [43, 45]. Specifically, we consider the fluctuating mobility problem

$$ \begin{bmatrix} M & -K \\ -K^T & 0 \end{bmatrix} \begin{bmatrix} \lambda \\ U \end{bmatrix} = \begin{bmatrix} u_s - u_{th} \\ -F \end{bmatrix}, \qquad (31) $$

with the random blob velocity,

$$ u_{th} = \sqrt{2k_B T M^{1/2}} Z, \qquad (32) $$

where $M^{1/2}$ is the square root of the blob mobility matrix, i.e. $M^{1/2} (M^{1/2})^T = M$, and $Z$ is a $3N_b \times 1$ white noise vector acting on the blobs, with covariance $\langle Z(t)Z^T(t') \rangle = I\delta(t - t')$. This form of the noise term is equivalent to including a stochastic stress directly into the Stokes equation as in (11). Since with the RPY approximation and grid-based approaches the mobility $M$ is positive definite [92, 95], the stochastic velocity (32) is well defined. The covariance of the resulting rigid body velocities generated by the stochastic blob velocities, setting $F = 0$ and $u_s = 0$,
is

$$ 
\begin{aligned}
\langle UU^T \rangle &= \langle NK^T M^{-1} u_{\text{th}} u_{\text{th}}^T M^{-1} KN \rangle \\
&= (NK^T M^{-1}) \sqrt{2k_B T} M^{1/2} \langle ZZ^T \rangle M^{1/2} \sqrt{2k_B T} (NK^T M^{-1})^T \\
&= 2k_B T N K^T M^{-1} M M^{-1} K N \delta(t-t') \\
&= 2k_B T N (KM^{-1}K^T) N \delta(t-t') \\
&= 2k_B T N \delta(t-t') = \langle U_{\text{th}} U_{\text{th}}^T \rangle.
\end{aligned} 
\quad (33) $$

Eq. (33) shows that the thermal noise is balanced by the viscous dissipation as required by the fluctuation-dissipation theorem (10). Thus the stochastic velocity imposed at the blob level (32) generates Brownian body velocities consistent with the Stokes equation. Equivalently, the square root of the body mobility matrix $N^{1/2} = NK^T M^{-1/2}$ satisfies the fluctuation-dissipation theorem with $N^{1/2}(N^{1/2})^T = N$.

## 2. Drift term

The second aspect of the computation is to account for Brownian drift when advancing the particle positions. It has been shown [43, 92, 96] that the drift can be incorporated using random finite differences (RFD), which a randomized version of the standard finite differences: given a Gaussian random vector, $W$, such that $\langle WW^T \rangle = I$, the divergence of the mobility matrix can be approximated as follows

$$ k_B T \partial_x \cdot N = k_B T \left\langle \frac{N(x + \delta W/2) - N(x - \delta W/2)}{\delta} \cdot W \right\rangle + O(\delta^2), \quad (34) $$

where $\delta$ a small number. This involves applying the mobility matrix $N$, and thus solving two mobility problems (31), evaluated at randomly displaced positions/orientations $\pm \delta W/2$ to a random force/torque vector $k_B T W$. Since the parameter $\delta$ is small the solution of one mobility problem can be used as initial guess to the second mobility problem reducing the computational cost [41].

In addition to RFD, a new midpoint scheme, called the “generalized Drifter-Corrector” (gDC), was recently developed [45]. This scheme, which is ideally suited for grid-based solvers, does not require additional mobility solves compared to the deterministic case and exhibit the same accuracy as RFD.

## Key contributions:

Introducing the noise at the blob level allows to generate the correct Brownian velocity increments at the level of the bodies with arbitrary shapes. The random blob velocities can be computed with any fast method, either grid-based or Green's function-based (see Section III C). Thanks to recent mathematical advances, including the drift term only requires one to two solves of the mobility problem (31) per time step. The resulting computational cost is negligible compared to traditional schemes, such as Fixman's method [98, 134].

## C. How to scale big and compute fast?

The linear system (31) needs to be solved efficiently to carry out dynamic simulations with large numbers of particles. To do so, two ingredients are crucial:

1. **An efficient iterative solver.** The size and condition number of the linear system (31) grows with the number of blobs $N_b$, it is therefore essential to use an iterative solver with a convergence rate that depends weakly on the number of particles. For that purpose block diagonal preconditioners applied to the iterative GMRES method have been developed [41, 42, 44]. The resulting solvers exhibit a convergence within a few iterations that is independent of the particle number $M$.

2. **Fast hydrodynamic interactions between blobs.** Hydrodynamic interactions between bodies are obtained directly from the matrix $M$. The action of $M$ on a vector is required to compute the deterministic velocities $U = NF - \tilde{N}u_s$, while the product $M^{1/2}Z$ is necessary to obtain the Brownian velocities $U_{\text{th}}$. It is therefore essential to use and develop efficient methods to compute the action of $M$ and its square root. In the following, we describe the approaches that have been developed for the fast evaluation of $M$ using Green's function and grid-based methods.
1. Fast hydrodynamics with Green's functions

Using direct linear algebra, the computational cost of a mobility-vector products $M\lambda$ scales quadratically with the number of blobs ($O(N_b^2)$) with Green's function based methods. In order to reduce this cost, matrix-vector products can be parallelized, either on CPU's or on the shared memory of GPU's, which yields a linear scaling up to $N_b \sim 10^4 - 10^5$ [43, 135-140]. More sophisticated methods are also used, such as fast multipole methods (FMM) [141], to achieve a linear scaling. For periodic domains, the positively split Ewald method [142, 143], the fast FCM [144] and periodic FMM [131, 132] are the most efficient and scalable methods for computing the action of $M$ on a vector.

To compute the action of $M^{1/2}$ on the random vector $Z$, iterative methods such as the Lanczos algorithm are best suited [94, 100]. This method iteratively approximates $M^{1/2}Z$ with a set of basis vectors of a vectorial space $K$ that are linear combinations of the powers of the mobility matrix times the random vector: $K = \text{span}\{Z, MZ, M^2Z, ...\}$. The cost of the method therefore depends on the number of basis vectors, and thus mobility-vector products, required to reach a given tolerance $\epsilon$. Using a block-diagonal preconditioner, the convergence rate is independent of the number of blobs and a small number of iterations (typically 3-5) is required to reach a reasonable tolerance.

2. Fast hydrodynamics with grid-based methods

When there is no explicit expression for the Green's function of the system, matrix-free methods can be used. They rely on grid-based Stokes solvers with fluctuating hydrodynamics combined with a large family of Immersed-Boundary (IBM) methods [92, 96, 101, 139, 145-148].

With fluctuating IBM methods, the coupling between the blobs and the fluid is achieved through a forcing term $f$ added to the fluctuating Stokes equations (11),

$$ \nabla p - \eta \nabla^2 v = \nabla \cdot \mathcal{Z} + f, \quad (35) $$

$$ \nabla \cdot v = 0. \quad (36) $$

The second term on the RHS of (35) is the blob forcing transferred to the fluid with a spreading operator $S$,

$$ f(r) = S[\lambda](r) = \sum_{i=1}^{N_b} \lambda_i \Delta_i(r), \quad (37) $$

where the finite size of the blobs is accounted for with a spreading envelope $\Delta_i(r)$. In the context of IBM, $\Delta_i(r)$ is a piece-wise regular compact kernel [50]. Other IBM-like methods, such as the Force Coupling Method (FCM), use Gaussian kernels [51].

The fluid velocity, obtained after solving (35) and (36), is the sum of a deterministic part $v^D$, due to the forcing $\lambda$, and a fluctuating term $v_{\text{th}}$, stemming from the fluctuating stress. We may express the total fluid velocity as

$$ v = \mathcal{L}^{-1}(S\lambda + \mathcal{D}\mathcal{Z}) \quad (38) $$

$$ = v^D + v_{\text{th}}, \quad (39) $$

where $S$ is the spreading operator in (37), $\mathcal{L}^{-1}$ is the inverse Stokes operator (i.e. the fluid solver), and $\mathcal{D}$ the divergence operator applied to the fluctuating stress [92].

The velocities of the blobs are then obtained from the fluid velocity using an averaging operator, $\mathcal{J}$, such that

$$ \mathbf{u}_i = (\mathcal{J}[v])_i = \int v \Delta_i(\mathbf{r}) d^3\mathbf{r} = \mathbf{u}_i^D + \mathbf{u}_{\text{th},i}, \quad i = 1, \dots, N_b, \quad (40) $$

where $\mathcal{J} = S^T$ is adjoint to the spreading operator.

We see then that the IBM mobility matrix can be written as the composition of three linear operators

$$ M_{\text{IBM}} = \mathcal{J} \mathcal{L}^{-1} S. \quad (41) $$

The IBM mobility differs from the RPY mobility defined in Eq. (27). However, one can perform a second order Taylor expansion of the flow field inside the kernels to show that $M_{\text{IBM}}$ obeys (28) to order $O(a^4)$. Thus, different kernels define different regularizations that are all essentially equivalent to the RPY mobility [92, 149].
Additionally, as demonstrated in [92, 95, 101] the velocity $\mathbf{u}_{th}$ satisfies the fluctuation-dissipation theorem with the covariance given by the IBM approximation of the mobility matrix,

$$ \langle \mathbf{u}_{th}(t) \mathbf{u}_{th}(t') \rangle = 2k_B T M_{IBM} \delta(t-t'), \quad (42) $$

where $\mathbf{u}_{th} = \sqrt{2k_B T} M_{IBM}^{1/2} Z(t)$, with $M_{IBM}^{1/2} = \mathcal{J} \mathcal{L}^{-1} \mathcal{D}$. To fulfill (42), it is crucial that the averaging and spreading operators are adjoint ($\mathcal{J} = S^T$) [150], which is not the case for the method of Regularized Stokeslets [35].

Thus, we see that fluctuating IBM methods simultaneously compute the actions of $M$ and $M^{1/2}$ on $\lambda$ and $Z$, respectively, by considering solutions to the forced fluctuating Stokes equations.

The fluctuating IBM can be implemented with any kind of Stokes solver. This involves first evaluating the forcing, $f(r)$, and the fluctuating stress, $Z$, on the grid. Next, the Stokes equations are solved with a Stokes solver (e.g. using Finite Volumes, Finite Elements, Finite Differences or a Spectral method). Finally, a quadrature rule is used to numerically integrate (40) to obtain the translational velocity for each of the blob making up the rigid bodies. Since many Stokes solvers can be well-parallelized on processor cores or GPUs, the computational cost scales linearly with the number of blobs $O(N_b)$ [139, 148, 151].

### Key contributions:

For both Green's function- and grid-based approaches, there are fast, parallel, methods for computing hydrodynamical interactions between blobs and their random noise, with linear scaling $O(N_b)$. Thanks to efficient preconditioned iterative solvers, the rigid multiblob linear system (31) can be solved with just a few iterations, i.e. mobility-vector products $M\lambda$, independently of the number of blobs in the system. Therefore the linear scaling of the computational cost is preserved at the level of the rigid bodies $O(M)$, which allows for simulations of complex suspensions at scales that had not been reached before.

## D. Handling kinematic constraints in large suspensions

As explained in Section II C, handling nonlinear holonomic constraints in large suspensions of bodies with arbitrary shape is challenging. However, when the constraints are holonomic their time derivatives are linear in the velocities of the rigid bodies. Therefore the time derivative of these constraints can be written compactly as a constraint equation on the particle velocities

$$ \dot{\mathbf{g}} = \mathbf{C}\mathbf{U} - \mathbf{B}(t) = \mathbf{0}, \quad (43) $$

where the matrix $\mathbf{C} = (\partial\dot{\mathbf{g}}/\partial\mathbf{U})$ is the constraints' Jacobian. For instance, the matrix $\mathbf{C}$ and vector $\mathbf{B}$ of the time derivative of the ball-and-socket constraint (19) are given by

$$ \begin{aligned} \mathbf{C} &= [\mathbf{I}_{3 \times 3} - \Delta \mathbf{l}_p^{\times} - \mathbf{I}_{3 \times 3} \ \Delta \mathbf{l}_q^{\times}], \\ \mathbf{B}(t) &= \mathbf{R}(\theta_q) \dot{\Delta} \mathbf{l}_q(t) - \mathbf{R}(\theta_p) \dot{\Delta} \mathbf{l}_p(t), \end{aligned} $$

where $\Delta \mathbf{l}_q^{\times} \mathbf{y} = \Delta \mathbf{l}_q \times \mathbf{y}$ for any vector $\mathbf{y}$ and $\dot{\Delta} \mathbf{l}_p(t) = \frac{\Delta \mathbf{l}_p}{\dot{\Delta} t}$.

Each constraint exerts a force and torque on the rigid bodies which in the case of passive links, $\mathbf{C}\mathbf{U} = \mathbf{B}(t) = \mathbf{0}$, generate no work [153]. This condition is enough to determine the structure of the constraint generalized forces,

$$ \mathbf{F}^{C} = \mathbf{C}^{T} \phi \quad (44) $$

where $\phi$ is a Lagrange multiplier homogeneous to a force [154].

Adding the constraint forces, $\mathbf{F}^C = \mathbf{C}^T \phi$, in the mobility relation (30), and inserting the velocity constraint (43) yields the linear system for the Lagrange multipliers $\phi$ and particle velocities $\mathbf{U}$ [44]

$$ \begin{bmatrix} -N \mathbf{C}^T & \mathbf{I} \\ \mathbf{0} & \mathbf{C} \end{bmatrix} \begin{bmatrix} \phi \\ \mathbf{U} \end{bmatrix} = \begin{bmatrix} N \mathbf{F} - \tilde{N} \mathbf{u}_s \\ \mathbf{B} \end{bmatrix}, \quad (45) $$

where $\mathbf{I}$ is a $6M \times 6M$ identity matrix. Any method to solve the Stokes problem, and therefore to apply the mobility matrix $\mathbf{N}$, can be used with (45). A similar framework to handle holonomic constraints in particle suspensions has been developed by [155] but with a different formulation of the linear system.

Solving (45) directly would require, in general, using nested linear solvers; an inner solver to compute the action of $\mathbf{N}$ and an outer solver to solve (45). Such approach would be computationally expensive. However, with the
<figure><img src="image_5.png" /><figcaption>FIG. 4. Illustration and validation of the RMB method with kinematic constraints. (a) Streamlines and fluid velocity magnitude around a swimming bacterium in free space discretized with 200 blobs. (b) Swimming speed of the bacterium normalized by the helical wave velocity $\omega/k$, where $k$ is the wavenumber of the helical wave and $\omega$ the rotation frequency, as a function of the number of wavelengths $N_{\lambda}$ along the helical tail for different ratios of flagellar length $L$ to head radius $R$. Dashed line: RMB method with kinematic constraints. Symbols: simulations from Higdon [152]. Figure adapted from [44].</figcaption></figure>

rigid multiblob model presented above, the whole linear system to solve the constrained mobility problem is found by combining (24)-(26), (43) and (44) [44]

$$ \begin{bmatrix} M & -K & 0 \\ -K^T & 0 & C^T \\ 0 & C & 0 \end{bmatrix} \begin{bmatrix} \lambda \\ U \\ \phi \end{bmatrix} = \begin{bmatrix} u_s \\ -F \\ B \end{bmatrix}. \qquad (46) $$

This linear system has the same structure as the unconstrained problem (29), which allows to use the same preconditioned iterative solvers (see Section III C). While resistance problems are known to scale poorly with the particle number, the iterative solver for the mixed mobility-resistance problem (46) is not sensitive to the system size, therefore allowing to study large suspensions of constrained bodies with quasilinear computational cost.

Additionally, constraint violations, e.g. due to discrete time-integration errors, are prevented with a body reconstruction that uses the constraints and a correction of the particles’ positions and orientations at the end of each time-step. The correction procedure, based on a nonlinear minimization algorithm, has negligible computational cost and preserves the accuracy of the time-integration scheme [44].

As a validation example, Figure 4a shows the flow field around a bacterium assembled with a spherical head and a helicoidal tail, both discretized with the multiblob method, that are attached and constrained to rotate in opposite directions along the same rotation axis. The resulting swimming speed shown in Figure 4b compares very well with other simulation results from the literature [152].

### Key contributions:

A velocity formulation of the kinematic constraints between bodies enables to write the constrained problem as a single linear system. This linear system is solved with a preconditioned iterative solver that couples effectively with any numerical method to compute hydrodynamic interactions between rigid bodies. The solver convergence is independent of the system size and constraint type, therefore allowing to simulate large suspensions in a scalable fashion. The resulting tool can handle different populations of constrained bodies simultaneously where both passive and time-dependent constraints are directly read from a simple input file. Thanks to its simplicity and flexibility, the user can readily use it to study physical and biological systems involving large collections of constrained bodies.

## E. The multi-physics multiblob method

Similarly to the Stokes equations, the Laplace equation (21) is an elliptic partial differential equations that can be written with an integral formulation at the particles’ surfaces [157] and discretized with the multiblob grid. The method has recently been extended to account for hydrochemical and phoretic interactions in suspensions of reactive particles with arbitrary shapes [46].
<figure><img src="image_6.png" /><figcaption>FIG. 5. Illustration and validation of the RMB method for reactive particles. a) Swimming speed for a Janus sphere as of function of the Damkohler number. Results for resolutions with different number of blobs, $N_b$, and comparison with the results of Michelin and Lauga [121]. The results with 42 blobs agree well for high Da numbers and shows an error of about 5% for Da ≪ 1; for finer resolutions the agreement is good for all Da numbers. Inset shows a the concentration field around a Janus particle (adapted from [121]). b) Velocities along the x-axis for two spheres with uniform surface activity as a function of the gap $d$ between the spheres. The motion is symmetric, thus only the velocity magnitude is shown. The results for discretizations with $N_b = 12, 42$ and $162$ are compared with multipolar expansions [156] and the exact solution computed in bispherical coordinates. Inset: sketch of the configuration. Figure adapted from [46].</figcaption></figure>

Just like the hydrodynamic problem, the multiblob method for the Laplace problem naturally interfaces with fast methods, such as GPU parallelization or the Fast Multipole Method [158], for computing convolutions with the Laplace Green's function and its derivatives (see Section III C). Canonical simulations of reactive particles show that the resulting method compares very well with analytical results from the literature with *only a few blobs* on the particle surface (see Fig. 5).

**Key contributions:**

The multiblob method can simulate large suspensions of particles of arbitrary shape with additional coupling, such as chemical reactions, while accounting for thermal fluctuations. It employs regularized kernels and a grid optimization strategy to solve the coupled Laplace-Stokes equations with reasonable accuracy at a fraction of the computational cost associated with traditional numerical methods.

The framework is not limited to diffusiophoresis, i.e. solute concentration fields. It broadly applies to any scalar field that satisfies elliptic equations, such as temperature or electric potential. With minor modifications, the code can be used to study the dynamics of thermophoretic or electrophoretic particle suspensions.

## IV. APPLICATIONS

In this section, we give a brief overview of applications that have been studied with the RMB method, sometimes in close collaboration with experimentalists. First, we present applications involving suspensions of active particles. These particles are driven either by external magnetic fields, by internal activity, or by hydrochemical interactions. Then we present systems involving large suspensions of colloids, where the particle shape and potential interactions are varied to tune the microstructure and rheological response of the suspension. All of these examples demonstrate the versatility of the RMB to study a wide range of systems in the fields of biophysics, soft and active matter, rheology, materials science and colloids.

### 1. *Active particles*

**Microrollers.** A microroller is a spherical colloid with a permanent magnetic moment that can be rotated by an external magnetic field. When rotated above a substrate they can *roll* in a direction determined by the external magnetic field. Systems formed by microrollers show an interesting interplay between the external drive, hydrodynamic interactions and thermal fluctuations. For example, a dense monolayer of colloids rolling above a flat substrate develops
<figure><img src="image_7.png" /><figcaption>FIG. 6. Illustrations of active particle systems simulated with the RMB framework. (a) Superposition of experimental (top) and numerical (bottom) images of a microroller interacting with an obstacle [59]. In both cases the microroller remains close to the obstacle for a while before the Brownian motion pushes the microroller away. (b) Lateral and top view of 100 bacteria swimming above a no-slip wall [44]. (c) Streamlines around a moving *Bacillaria* colony formed by 16 cells. (d) Lateral view of a self-propelled nanorod swimming along a tilted substrate [58]. The blobs of different colloids represent two different metals, the experiments used gold-platinum and gold-rhodium nanorods. The red (blue) shaded areas show the regions of high (low) pressure while the continuous lines represent the streamlines and the white arrows the active slip imposed in the nanorod. (e) Flow field around phoretic chiral particles [46]. Color-scale from white to dark red: magnitude of the velocity field. Black lines with red arrows: streamlines. Inset: experimental reactive stirrer, adapted from [161].</figcaption></figure>

a fingering instability that eventually seeds motile clusters of colloids which are stable for long periods of time [54]. The origin of the instability is the strong flows generated by the rotating microrollers, however, its growth rate and the size of the seeded clusters is controlled by the height of the microrollers above the substrate [159]. The height, in turn, is affected by the Brownian motion that lifts the colloids above the substrate [55, 160]. Therefore, to accurately capture the fingering instability and the emergence of clusters it is necessary to solve the hydrodynamic problem with thermal fluctuations.

The interplay between Brownian noise and hydrodynamic interactions affects the dynamics of microrollers in other ways. For example, a microroller can become trapped by a cylindrical obstacle, with the trapping time increasing with the obstacle size, see Fig. 6a [59]. The RMB method was used to understand this phenomenon by discretizing both the microroller and the obstacle with a small number of blobs. and by solving a mixed mobility-resistance problem in which the obstacle was fixed while the microroller was free to move under the action of an external magnetic field. Thanks to simulations, it was found that the rotation of the microroller creates flows that can be one or two orders of magnitude greater than the self-induced velocity, and that the hydrodynamic response of the obstacle to these flows strongly attracts the microroller towards the obstacle. The balance between this hydrodynamic attraction and the roller self-induced velocity delimits a basin of attraction that microrollers can enter or exit only through the diffusion generated by their Brownian motion. The trapping mechanism is therefore quite unique: noise, here thermal fluctuations, is the only way to reach an hydrodynamic attractor [59]. The basin of attraction size increases with the obstacle size which explains the escape times, as Brownian displacements are less likely to move a microroller out of a large basin of attraction. By solving the hydrodynamic and Brownian problems simultaneously, the RMB method successfully reproduced and explained the experimental results [59].

**Microorganisms.** The versatility of the RMB method, combined with fast methods to evaluate hydrodynamic interactions (cf. Section III C), allows to study a plethora of large constrained systems within a unified framework. For example, some microorganisms have developed mechanisms that break the kinematic reversibility of Stokes flows
to self-propel. The classical examples are flagellated bacteria (see Fig. 1d) and cilia covered organisms, but there are others, such as diatom colonies, that move by their relative displacements [162]. These organisms can be modeled as articulated rigid bodies with the formalism discussed in Secs. IIC and IIID, and the RMB method has been used to study some of them.

For instance, a bacterium can be modeled as a rigid body with a rigid flagellum attached by inextensible links so that the flagellum is free to rotate around its main axis but otherwise moves with the bacterium body. Within this formalism it is possible to apply equal but opposite torques to the bacterium body and flagellum to mimic the effect of the molecular motor acting on the flagellum, so the whole bacterium swims like a free-force swimmer (see Fig. 4). This approach was used with the RMB method to model the swimming dynamics of 100 bacteria above a no-slip wall, see Fig. 6b [44]. The bacteria were hydrodynamically attracted by the no-slip wall [163] and formed a large cluster, as shown in the snapshot in Fig. 6b, which eventually broke up. Such a large simulation, with 61,800 unknowns, was possible because the mobility problem was solved in a moderate number of iterations, because the hydrodynamic interactions were computed with a fast method [131, 158], and because it was possible to use large time step sizes as the inextensible links did not introduce any stiffness into the equations of motion.

The articulated body formalism was also used to study the dynamics of *Bacillaria* colonies. *Bacillaria* is a rod-shaped diatom about 70µm long and 10µm thick that habits marine and fresh water ecosystems [162]. *Bacillaria* cells can attach to each other along their long axes forming a stack of rods, see Figs. 1c and 6c. Once attached, the cells have the ability to slide with respect their neighbors in a coordinated fashion. In Ref. [44] *Bacillaria* cells were modeled as rods formed by only 14 blobs and the dynamics of a colony formed by 16 cells was studied. The cells were forced to slide against their first neighbors with a sinusoidal pattern and thus, the whole colony deformed following a wave form. Figure 6c shows a typical snapshot of the colony and the induced flow field. Unlike typical beating flagella, it was found the swimming direction and speed of the colony varied non-monotonically with the wavelength of the deformation wave [44]. These nontrivial results are currently investigated and experiments are carried out in parallel.

**Phoretic nanorods.** Another problem, studied with the RMB method, where the interplay between hydrodynamic interactions and Brownian motion is important is the dynamics of phoretic nanorods. Bimetallic nanorods can self-propel through electrophoresis in aqueous H₂O₂ solutions [164]. These nanorods, similar to the ones shown in Fig. 1g, are formed by two segments of different metals, typically gold and platinum, and decompose the H₂O₂ with asymmetric reactions on each metallic segment [165]. The chemical reactions produce a flux of ions near the nanomotor surface that drags the fluid and by momentum conservation propels the nanorod in the opposite direction. These nanorods have a tendency to swim against flow currents and also against gravity along tilted walls, i.e. they show rheotaxis and gravitaxis [57, 58]. RMB simulations were used to understand these behaviors. In those simulations the nanorods were discretized with a moderated number of blobs (~100). Instead of solving the complex electrochemical reactions an active slip velocity was set *by hand* near the bimetallic interface, see Fig. 6d.

The solution of the hydrodynamic problem showed that the nanorods swim with a tilt towards the substrate, see Fig. 6d. This tilt makes the nanorods behave like a weather vane and it is the origin of their rheotaxis. In the presence of a shear flow above the substrate a tilted rod reorients itself against the flow and once reoriented it tends to swim upstream. Pure deterministic simulations predict that a nanorod can swim upstream whenever its intrinsic swimming speed is larger than the flow speed at the height of the nanorod above the substrate. However, Brownian motion randomizes the nanorod orientation which makes it a less efficient rheotactor. RMB simulations were able to quantitatively reproduce the experimental upstream swimming velocities by solving the hydrodynamic and Brownian problems simultaneously [57]. A similar effect was observed for heavy-bottom nanorods swimming on inclined substrates [58]. Heavy-bottom nanorods reorient against gravity and swim up showing gravitaxis while Brownian motion randomizes their orientations and thus reduce their gravitactic efficiency. Interestingly, RMB simulations predicted that density homogeneous rods also experience a reorienting torque. A nanorod with a tilt towards the substrate rotates around a point displaced towards its head, therefore, even if the center of mass is at its geometric center the nanorod experiences a reorienting torque. These numerical findings were confirmed experimentally [58].

**Phoretic chiral particles.** Thanks to the scalability and versatility of the method, we have been able to explore new physical problems that, to the best of our knowledge, have not been addressed with numerical simulations. Fig. 6e shows an example, inspired by recent experimental work [8, 166, 167], consisting of a large collection of chiral particles, called *reactive stirrers*. We used the approach described in Secs. IID and IIIE to model these particles as reacting particles with uniform surface properties that emit solute at a fixed rate [46]. We confined the particles with a harmonic potential above a no-slip wall. Thanks to its chirality, a single particle rotates spontaneously due to self-phoresis and stirs the surrounding fluid. We found that a large collection of such particles exhibit nontrivial collective behaviors, such as the coexistence of a rotating outer rim and an inner suspension of particles rotating in opposite directions, see Fig. 6e [46].
<figure><img src="image_8.png" /><figcaption>FIG. 7. Examples of complex colloidal suspensions simulated with the RMB framework. (a) Snapshot of star shaped colloids under shear flow. A selected number of colloids are shown with full colors while the rest are shown with a transparency to easy the visualization. The colloids are colored according to their initial position along the x-axis to reveal the shear flow [67]. (b) Snapshot of the sheared branched particles in a periodic domain. Colors are chosen randomly [45]. (c) Suspension of patchy colloids with patch depended attraction-repulsion interactions [124].</figcaption></figure>

## 2. Complex colloidal suspensions

**Rheology of complex shaped colloids.** Some colloids are formed by the aggregation of smaller particles. For example, spherical silica particles can fuse to form complex colloids [2], see Fig. 1b. Suspensions formed by such colloids have an interesting rheology, exhibiting high viscosities at low volume fractions and continuous or discontinuous shear thickening [2]. Modeling such systems with traditional methods is too costly due to the complex shapes of the colloids and the role of Brownian motion in their rheological response. However, within the RMB's framework each silica particle can be modeled as a single blob which reduces the numerical complexity and cost, see Fig. 7a-b. RMB studies of star colloids suspensions have shown that their viscosity can be quite large at low and moderate volume fractions in agreement with experimental results. A hundred fold increase in the viscosity at 20% volume fraction was observed when the colloids have complex shapes, like hooks, that allow them to link to each other and form a network that resists stress [67]. However, in contrast with experimental results that found a shear thickening response [2], RMB studies found shear thinning [45, 67]. Such difference may come from the lack of friction and lubrication forces in the simulations, although both of them can be incorporated into the RMB formulation [41].

Friction has already been used in the context of Stokesian Dynamics and it can be included as local pairwise force between blobs of different colloids [168]. Since friction forces do not modify the linear systems discussed in Sec. ?? their implementation is simple, although, it may introduce stiffness into the equations of motion. Including lubrication corrections modifies the mobility problem (29), however, the additional complexity is small. The idea of including lubrication corrections was first introduced in the Stokesian dynamics method [13, 169] and it has been discussed in the context of rigid assemblies of blobs [41, 55].

**Percolation and gelation of patchy colloids.** The same numerical approach discussed in this manuscript was used by Wang & Swan to study the percolation and gelation of patchy colloids [124] (see Fig. 1a for an experimental example). They studied spherical colloids whose surfaces were covered with two types of patches exhibiting different attractive and repulsive interactions, see Fig. 7c. Each colloid was discretized with 42 blobs, which were assigned to one of the two patch types. Using this moderate resolution, they were able to study the dynamics of suspensions consisting of 500 colloids under thermal fluctuations and shear forcing. They found that, at moderate concentrations, patchy colloids can form microstructures quite different from those of homogeneous particles, which affect the percolation transition points and the suspensions rheology.

# V. PERSPECTIVES

Thanks to its simplicity and flexibility, the rigid multiblob (RMB) method has been used to simulate a wide variety of systems involving complex suspended particles. However, several challenges remain in generalizing the formulation to cover additional systems of interest and improving its efficiency and accuracy. We discuss some of these challenges here.
<figure><img src="image_9.png" /><figcaption>FIG. 8. a) Sketch of the grid optimization procedure for a body with arbitrary shape (in gray) and position $\mathbf{q}_p$. Left: surface discretization of a body with arbitrary shape with $N_b$ blobs of radius $a$ (yellow disks) and size scale $S = 1$. Right: optimized grid with $N_b$ larger blobs (red disks) positioned inside the actual body surface ($S < 1$). b) Relative error of the propulsion speed of a squirmer (sketched in the inset) as a function of the grid resolution $N_b \in [12; 2562]$. Disks: empirical grid. Triangles: optimized grid.</figcaption></figure>

The formulation discussed in Sec. II A can be generalized in several ways. For example, the no-slip boundary conditions can be extended to partial-slip or Navier boundary conditions [170]. Navier boundary conditions assume that the flow velocity at solid surfaces does not match the surface velocity but instead *slips* with a velocity proportional to the local velocity gradient. The proportionality factor between the velocity gradient and the slip velocity is the slip length, which represents the distance across the boundary where the linearly extrapolated flow velocity would vanish. At nanoscales, even small slip lengths can significantly affect the dynamics of suspended particles [171]. Boundary integral formulations for solving the Stokes equations with Navier boundary conditions exist [79, 171, 172], and the RMB can incorporate them into its formalism [173]. These formulations involve the double-layer operator of the Stokes equations, which is not positive definite. Therefore, the approach used in Sec. IIIB to generate Brownian velocities is not directly applicable, and how to efficiently generate Brownian velocities in this setting remains an open problem.

In this manuscript, we have only discussed problems where the suspended particles are rigid objects. However, blobs can also be used to model deformable objects, such as fibers. Assemblies of rigid and elastic objects can thus be built with blobs without any additional complexity. Elastic interactions can introduce stiffness into the equations of motion, which may require specialized integrators [174]. This approach has been used, for example, to study the coordinated motion of cilia on a spherical microorganism [175].

The methods used to compute hydrodynamic interactions in Sec. III C 1 rely on the Green's function of the Stokes equations, which requires discretizing only the particle surfaces rather than the full three-dimensional fluid domain. However, these methods are difficult to generalize to flows not governed by elliptic partial differential equations, i.e., flows beyond the Stokes regime, such as viscoelastic flows or flows at finite Reynolds numbers. In such cases, IBM-like methods, as discussed in Sec. III C 2, can be used to solve the flow dynamics. Examples modeling complex particles at finite Reynolds numbers can be found in Refs. [96, 149].

Hydrodynamic interactions between spherical blobs are introduced through the mobility matrix $M$, and the efficiency in computing these interactions impacts the overall performance of the RMB. Numerous codes are available to efficiently compute the action of $M$ (and of $M^{1/2}$ in some cases) [41, 92, 132, 137–140, 144, 176]. These codes differ in the domain geometry they support (periodic, bounded, infinite, etc.) and in the programming languages used (C, C++, CUDA, Fortran90, Python, JAX). To be invoked by the multiblob method, these codes must be linked as external libraries. A recent initiative called LibMobility aims to centralize various methods for computing hydrodynamic interactions between blobs in a user-friendly framework. The resulting library is designed for streamlined installation and easy integration into existing codes [177]. Such efforts are arduous but invaluable to the community and should be encouraged, even if they remain *poorly rewarded* in terms of scientific recognition.

To simulate large numbers of complex particles, it is necessary to use as few blobs per particle as possible while maintaining reasonable accuracy in the hydrodynamic interactions. Despite the regularization error introduced by the regularized Green's function (27), which scales with the blob radius $a$ and thus as $N_b^{-1/2}$ [46], it is possible to enhance numerical accuracy for a given resolution through a grid optimization procedure. Once the number of blobs per particle and the surface mesh are selected, two additional grid parameters remain to be tuned: the blob radius $a$ and the position of the blobs relative to the actual particle surface. Due to their finite size, blobs do not need to lie exactly on the particle surface. The distance of blob centers to the particle surface is described by a geometric parameter $S$, representing the scale of the discretized surface relative to the actual particle surface: when $S = 1$, the surface of the
grid coincides with the surface of the actual particle; if $S < 1$ ($S > 1$), the blobs are inside (outside) the surface (see Fig. 8a). In the early works on the multiblob method, this tuning was done empirically and only for the mobility matrix $N$ [12, 42]. More recently, Broms et al. [127] formalized the optimization of $N$ through a minimization problem. Building on their approach, we proposed a new optimization strategy based on the singular value decomposition of $\bar{N}$ to match the surface mobility operator as well [46]. Figure 8b shows a comparison of the propulsion speed of a spherical squirmer using the empirically chosen grid [42] and the optimized grid with the exact solution $U_{\text{th}}$. The optimized grid (triangles) significantly reduces the error compared to the original discretization (disks), enabling the simulation of slip-driven particles with coarse grids while achieving 2–3 digits of accuracy—comparable to that of other methods using much finer grids [79]. This grid optimization strategy is applicable to particles of any shape. One of the key challenges remaining is the development of quadrature rules that can balance, or at least further reduce, the regularization error already present in the continuous setting, with the goal of achieving a convergence rate faster than $\sim N_b^{-1/2}$ for the rigid multiblob method. An interesting alternative lies in recent developments of the method of fundamental solutions, which, although still much more computationally expensive than the RMB, show promising results in terms of accuracy and scalability [178].

In conclusion, the RMB is a versatile and robust framework for simulating complex particles in Stokes flows. Although there is still room for improvement in its accuracy, efficiency and range of application, it has now reached a stage of maturity that allows it to tackle a variety of physical problems that cannot be addressed using traditional methods.

ACKNOWLEDGMENTS

The authors would like to pay tribute to Aleks Donev, who was a key contributor to the rigid multiblob method. We also thank the anonymous reviewers for their insightful comments. B.D. acknowledges support from the French National Research Agency (ANR), under award ANR-20-CE30-0006. F.B.U. acknowledges funding by the Basque Government through the BERC 2022-2025 program and by the Ministry of Science, Innovation and Universities: BCAM Severo Ochoa accreditation CEX2021-001142-S/MICIN/AEI/10.13039/501100011033.

[1] W. Li, H. Palis, R. Mérindol, J. Majimel, S. Ravaine, and E. Duguet, Colloidal molecules and patchy particles: Complementary concepts, synthesis and self-assembly, Chemical Society Reviews **49**, 1955 (2020).

[2] P. Bourrianne, V. Niggel, G. Polly, T. Divoux, and G. H. McKinley, Tuning the shear thickening of suspensions through surface roughness and physico-chemical interactions, Physical Review Research **4**, 033062 (2022).

[3] iNaturalist, Photo 11624974, https://www.inaturalist.org/photos/11624974 (2017).

[4] K. A. Syed and K. E. Klose, Vibrio cholerae flagellar synthesis and virulence, in *Epidemiological and Molecular Aspects on Cholera* (Springer, 2010) pp. 203–212.

[5] S. Sacanna, W. T. Irvine, P. M. Chaikin, and D. J. Pine, Lock and key colloids, Nature **464**, 575 (2010).

[6] A. McMullen, M. Holmes-Cerfon, F. Sciortino, A. Y. Grosberg, and J. Brujic, Freely Jointed Polymers Made of Droplets, Physical Review Letters **121**, 138002 (2018).

[7] M. S. D. Wykes, J. Palacci, T. Adachi, L. Ristroph, X. Zhong, M. D. Ward, J. Zhang, and M. J. Shelley, Dynamic self-assembly of microscale rotors and swimmers, Soft matter **12**, 4584 (2016).

[8] A. M. Brooks, M. Tasinkevych, S. Sabrina, D. Velegol, A. Sen, and K. J. Bishop, Shape-directed rotation of homogeneous micromotors via catalytic self-electrophoresis, Nature communications **10**, 495 (2019).

[9] W. Yan, E. Corona, D. Malhotra, S. Veerapaneni, and M. Shelley, A scalable computational platform for particulate Stokes suspensions, Journal of Computational Physics **416**, 109524 (2020).

[10] L. Crowder, T. Li, E. Corona, and S. Veerapaneni, Boundary integral equation analysis for spheroidal suspensions, arXiv preprint arXiv:2506.20809 (2025).

[11] Y. Bao, M. Rachh, E. E. Keaveny, L. Greengard, and A. Donev, A fluctuating boundary integral method for Brownian suspensions, Journal of Computational Physics **374**, 1094 (2018).

[12] J. W. Swan and G. Wang, Rapid calculation of hydrodynamic and transport properties in concentrated solutions of colloidal particles and macromolecules, Physics of Fluids **28**, 011902 (2016).

[13] J. F. Brady and G. Bossis, Stokesian dynamics, Annu. Rev. Fluid Mech. **20**, 111 (1988).

[14] J. McCammon and J. Deutch, Frictional properties of nonspherical multisubunit structures. application to tubules and cylinders, Biopolymers: Original Research on Biomolecules **15**, 1397 (1976).

[15] J. Rotne and S. Prager, Variational treatment of hydrodynamic interaction in polymers, Journal of Chemical Physics **50**, 4831 (1969).

[16] H. Yamakawa, Transport properties of polymer chains in dilute solution: hydrodynamic interaction, The Journal of Chemical Physics **53**, 436 (1970).
[17] H. Nakajima and Y. Wada, A general method for evaluation of diffusion constants, dilute-solution viscoelasticity, and the dielectric property of a rigid macromolecule with an arbitrary configuration. i, Biopolymers: Original Research on Biomolecules **16**, 875 (1977).

[18] J. G. De La Torre and V. A. Bloomfield, Hydrodynamic properties of macromolecular complexes. i. translation, Biopolymers: Original Research on Biomolecules **16**, 1747 (1977).

[19] J. G. De La Torre and V. A. Bloomfield, Hydrodynamics of macromolecular complexes. ii. rotation, Biopolymers: Original Research on Biomolecules **16**, 1765 (1977).

[20] J. G. de la Torre and V. A. Bloomfield, Hydrodynamic properties of complex, rigid, biological macromolecules: theory and applications, Quarterly reviews of biophysics **14**, 81 (1981).

[21] M. Zurita-Gotor, J. Bławzdziewicz, and E. Wajnryb, Motion of a rod-like particle between parallel walls with application to suspension rheology, Journal of Rheology **51**, 71 (2007).

[22] T. T. Bringley and C. S. Peskin, Validation of a simple method for representing spheres and slender bodies in an immersed boundary method for stokes flow on an unbounded domain, Journal of Computational Physics **227**, 5397 (2008).

[23] A. Pandey, P. S. Kumar, and R. Adhikari, Flow-induced nonequilibrium self-assembly in suspensions of stiff, apolar, active filaments, Soft Matter **12**, 9068 (2016).

[24] B. Cichocki and K. Hinsen, Stokes drag on conglomerates of spheres, Physics of Fluids **7**, 285 (1995).

[25] M. X. Fernandes and J. G. de la Torre, Brownian dynamics simulation of rigid particles of arbitrary shape in external fields, Biophysical journal **83**, 3039 (2002).

[26] J. G. de la Torre, M. L. Huertas, and B. Carrasco, Calculation of hydrodynamic properties of globular proteins from their atomic-level structure, Biophysical journal **78**, 719 (2000).

[27] V. Lobaskin and B. Dünweg, A new model for simulating colloidal dynamics, New Journal of Physics **6**, 54 (2004).

[28] R. Kutteh, Rigid body dynamics approach to Stokesian dynamics simulations of nonspherical particles, The Journal of chemical physics **132** (2010).

[29] M. Długosz and J. M. Antosiewicz, Hydrodynamic effects on the relative rotational velocity of associating proteins, The Journal of Physical Chemistry B **117**, 6165 (2013).

[30] S. Poblete, A. Wysocki, G. Gompper, and R. G. Winkler, Hydrodynamics of discrete-particle models of spherical colloids: A multiparticle collision dynamics simulation study, Physical Review E **90**, 033314 (2014).

[31] A. Ortega, D. Amorós, and J. G. De La Torre, Prediction of hydrodynamic and other solution properties of rigid proteins from atomic-and residue-level models, Biophysical journal **101**, 892 (2011).

[32] A. Gastaldi and M. Vanni, The distribution of stresses in rigid fractal-like aggregates in a uniform flow field, Journal of colloid and interface science **357**, 18 (2011).

[33] J. J. Molina and R. Yamamoto, Direct numerical simulations of rigid body dispersions. i. mobility/friction tensors of assemblies of spheres, The Journal of chemical physics **139** (2013).

[34] M. Reichert and H. Stark, Synchronization of rotating helices by hydrodynamic interactions, The European Physical Journal E **17**, 493 (2005).

[35] R. Cortez, L. Fauci, and A. Medovikov, The method of regularized stokeslets in three dimensions: Analysis, validation, and application to helical swimming, Physics of Fluids **17**, 031504 (2005).

[36] A. Codutti, F. Bachmann, D. Faivre, and S. Klumpp, Bead-based hydrodynamic simulations of rigid magnetic micropro-ellers, Frontiers in Robotics and AI **5**, 109 (2018).

[37] S. Bianchi, V. Carmona Sosa, G. Vizsnyiczai, and R. Di Leonardo, Brownian fluctuations and hydrodynamics of a microhelix near a solid wall, Scientific Reports **10**, 4609 (2020).

[38] K. Hinsen, Hydrolib: a library for the evaluation of hydrodynamic interactions in colloidal suspensions, Computer physics communications **88**, 327 (1995).

[39] J. W. Swan, J. F. Brady, R. S. Moore, et al., Modeling hydrodynamic self-propulsion with Stokesian dynamics. or teaching Stokesian dynamics to swim, Physics of Fluids **23** (2011).

[40] S. Delong, F. B. Usabiaga, and A. Donev, Brownian dynamics of confined rigid bodies, The Journal of chemical physics **143**, 144107 (2015).

[41] A. M. Fiore and J. W. Swan, Fast Stokesian dynamics, Journal of Fluid Mechanics **878**, 544 (2019).

[42] F. Balboa Usabiaga, B. Kallemov, B. Delmotte, A. P. S. Bhalla, B. E. Griffith, and A. Donev, Hydrodynamics of suspensions of passive and active rigid particles: a rigid multiblob approach, Communications in Applied Mathematics and Computational Science **11**, 217 (2016).

[43] B. Sprinkle, F. Balboa Usabiaga, N. A. Patankar, and A. Donev, Large scale Brownian dynamics of confined suspensions of rigid particles, The Journal of Chemical Physics **147**, 244103 (2017), https://doi.org/10.1063/1.5003833.

[44] F. Balboa Usabiaga and B. Delmotte, A numerical method for suspensions of articulated bodies in viscous flows, Journal of Computational Physics **464**, 111365 (2022).

[45] T. A. Westwood, B. Delmotte, and E. E. Keaveny, A generalised drift-correcting time integration scheme for Brownian suspensions of rigid particles with arbitrary shape, Journal of Computational Physics **467**, 111437 (2022).

[46] B. Delmotte and F. B. Usabiaga, A scalable method to model large suspensions of colloidal phoretic particles with arbitrary shapes, Journal of Computational Physics **518**, 113321 (2024).

[47] H. Gidituri, G. Kabacaoğlu, M. Ellero, and F. B. Usabiaga, Mapping flagellated swimmers to surface-slip driven swimmers, Journal of Computational Physics **510**, 113081 (2024).

[48] E. Krucker-Velasquez, J. W. Swan, and Z. Sherman, Immersed boundary method for dynamic simulation of polarizable colloids of arbitrary shape in explicit ion electrolytes, arXiv preprint arXiv:2401.04865 (2024).

[49] R. Cortez, The method of regularized stokeslets, SIAM Journal on Scientific Computing **23**, 1204 (2001).
[50] C. S. Peskin, The immersed boundary method, Acta numerica 11, 479 (2002).

[51] M. Maxey and B. Patel, Localized force representations for particles sedimenting in Stokes flow, International journal of multiphase flow 27, 1603 (2001).

[52] T. Dombrowski, S. K. Jones, G. Katsikis, A. P. S. Bhalla, B. E. Griffith, and D. Klotsa, Transition in swimming direction in a model self-propelled inertial swimmer, Physical Review Fluids 4, 021101 (2019).

[53] H. Gidituri, M. Ellero, and F. Balboa Usabiaga, Swimming efficiently by wrapping, Journal of Fluid Mechanics 993, A7 (2024).

[54] M. Driscoll, B. Delmotte, M. Youssef, S. Sacanna, A. Donev, and P. Chaikin, Unstable fronts and motile structures formed by microrollers, Nature Physics 13, 375 (2017).

[55] B. Sprinkle, E. B. Van Der Wee, Y. Luo, M. M. Driscoll, and A. Donev, Driven dynamics in dense suspensions of microrollers, Soft Matter 16, 7982 (2020).

[56] B. Sprinkle, S. Wilken, S. Karapetyan, M. Tanaka, Z. Chen, J. R. Cruise, B. Delmotte, M. M. Driscoll, P. Chaikin, and A. Donev, Sedimentation of a colloidal monolayer down an inclined plane, Physical Review Fluids 6, 034202 (2021).

[57] Q. Brosseau, F. Balboa Usabiaga, E. Lushi, Y. Wu, L. Ristroph, J. Zhang, M. Ward, and M. J. Shelley, Relating rheotaxis and hydrodynamic actuation using asymmetric gold-platinum phoretic rods, Phys. Rev. Lett. 123, 178004 (2019).

[58] Q. Brosseau, F. B. Usabiaga, E. Lushi, Y. Wu, L. Ristroph, M. D. Ward, M. J. Shelley, and J. Zhang, Metallic microswimmers driven up the wall by gravity, Soft Matter 17, 6597 (2021).

[59] E. B. van der Wee, B. C. Blackwell, F. Balboa Usabiaga, A. Sokolov, I. T. Katz, B. Delmotte, and M. M. Driscoll, A simple catch: Fluctuations enable hydrodynamic trapping of microrollers by obstacles, Science Advances 9, eade0320 (2023).

[60] S.-Y. Chen, H. M. L. Rios, M. O. de la Cruz, and M. Driscoll, Restructuring a passive colloidal suspension using a rotationally driven particle, Soft Matter 20, 2151 (2024).

[61] M. Dlugosz, B. Cichocki, and P. Szymczak, Estimating near-wall diffusion coefficients of arbitrarily shaped rigid macromolecules, Physical Review E 106, 014407 (2022).

[62] M. Dlugosz, B. Cichocki, and P. Szymczak, First coarse grain then scale: How to estimate diffusion coefficients of confined molecules, The Journal of Chemical Physics 159 (2023).

[63] D. Kundu, F. Balboa Usabiaga, and M. Ellero, Settling dynamics of a non-Brownian suspension of spherical and cubic particles in Stokes flow, Journal of Fluid Mechanics 1010, A63 (2025).

[64] N. Moreno, D. Moreno-Chaparro, F. B. Usabiaga, and M. Ellero, Hydrodynamics of spike proteins dictate a transport-affinity competition for sars-cov-2 and other enveloped viruses, Scientific Reports 12, 11080 (2022).

[65] D. Moreno-Chaparro, N. Moreno, F. B. Usabiaga, and M. Ellero, Computational modeling of passive transport of functionalized nanoparticles, The Journal of Chemical Physics 158 (2023).

[66] N. Moreno, D. Vazquez-Cortes, and E. Fried, Sedimentation dynamics of triply-twisted möbius bands: Geometry versus topology, arXiv preprint arXiv:2402.13017 (2024).

[67] F. Balboa Usabiaga and M. Ellero, Rheology of moderated dilute suspensions of star colloids: The shape factor, Physics of Fluids 36 (2024).

[68] U. Makanga, M. Sepahi, C. Duprat, and B. Delmotte, Obstacle-induced lateral dispersion and nontrivial trapping of flexible fibers settling in a viscous fluid, Physical Review Fluids 8, 044303 (2023).

[69] U. Makanga, C. Duprat, and B. Delmotte, Sedimentation of flexible fibers through ordered arrays of pillars, To be submitted (2024).

[70] Y. Gao, B. Sprinkle, E. Springer, D. W. Marr, and N. Wu, Rolling of soft microbots with tunable traction, Science Advances 9, eadg0919 (2023).

[71] E. S. Bililign, F. Balboa Usabiaga, Y. A. Ganan, A. Poncet, V. Soni, S. Magkiriadou, M. J. Shelley, D. Bartolo, and W. T. Irvine, Motile dislocations knead odd crystals into whorls, Nature Physics 18, 212 (2022).

[72] T. Yang, B. Sprinkle, Y. Guo, J. Qian, D. Hua, A. Donev, D. W. Marr, and N. Wu, Reconfigurable microbots folded from simple colloidal chains, Proceedings of the National Academy of Sciences 117, 18186 (2020).

[73] X. Wang, B. Sprinkle, H. K. Bisoyi, T. Yang, L. Chen, S. Huang, and Q. Li, Colloidal tubular microrobots for cargo transport and compression, Proceedings of the National Academy of Sciences 120, e2304685120 (2023).

[74] S. Kim and S. J. Karrila, Microhydrodynamics: principles and selected applications, Butterworth-Heinemann series in chemical engineering (Butterworth-Heinemann, Boston, 1991) oCLC: 246965391.

[75] J. L. Anderson, Colloid transport by interfacial forces, Annu. Rev. Fluid Mech. 21, 61 (1989).

[76] R. Golestanian, T. Liverpool, and A. Ajdari, Designing phoretic micro-and nano-swimmers, New Journal of Physics 9, 126 (2007).

[77] J. Blake, A spherical envelope approach to ciliary propulsion, Journal of Fluid Mechanics 46, 199 (1971).

[78] C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow, Cambridge Texts in Applied Mathematics (Cambridge University Press, 1992).

[79] D. J. Smith, M. T. Gallagher, R. Schuech, and T. D. Montenegro-Johnson, The role of the double-layer potential in regularised stokeslet models of self-propulsion, Fluids 6, 411 (2021).

[80] The adjoint property ensures that the power dissipated by the motion of the particle surface integrated over the whole surface matches the power dissipated by the motion of a rigid particle. For a particle without surface velocity ($u_s = 0$)

$$P = \int_{\partial \mathcal{B}_p} \lambda^T(\mathbf{r}) v(\mathbf{r}) \, dS_r = \int_{\partial \mathcal{B}_p} \lambda^T(\mathbf{r}) (\mathcal{K} \mathcal{U})(\mathbf{r}) \, dS_r = \mathcal{U}^T \mathcal{K}^T \lambda = \mathcal{U}^T \mathcal{F} = \mathbf{u} \cdot \mathbf{f} + \tau \cdot \omega. \quad (47)$$
[81] $L^2(\partial\mathcal{B}, \mathbb{R}^3)$ is the space of vector-valued functions in $\mathbb{R}^3$, i.e. $\mathbf{u}_s$, defined on the manifold $\partial\mathcal{B}$ (the particle surface) whose euclidean norm is square integrable on $\partial\mathcal{B}$, i.e. $\int_{\partial\mathcal{B}} \|\mathbf{u}_s\|_2^2 dS < \infty$. We think this function space is a safe choice for $\mathbf{u}_s$, but it might not be the most appropriate.

[82] H. A. Stone and A. D. T. Samuel, Propulsion of microorganisms by surface distortions, Phys. Rev. Lett. **77**, 4102 (1996).

[83] R. Kubo, The fluctuation-dissipation theorem, Reports on progress in physics **29**, 255 (1966).

[84] W. B. Russel, Brownian motion of small particles suspended in liquids, Annual Review of Fluid Mechanics **13**, 425 (1981).

[85] S. R. de Groot and P. Mazur, *Non-Equilibrium thermodynamics* (Dover Publications Inc., 1984).

[86] E. Guazzelli and J. F. Morris, *A physical introduction to suspension dynamics*, Vol. 45 (Cambridge University Press, 2011).

[87] L. Landau and E. Lifshitz, *Fluid Mechanics* (Pergamon Press, 1959).

[88] G. D. Fabritiis, M. Serrano, R. Delgado-Buscalioni, and P. V. Coveney, Fluctuating hydrodynamic modeling of fluids at the nanoscale, Physical Review E **75**, 026307 (2007).

[89] E. J. Hinch, Application of the Langevin equation to fluid suspensions, Journal of Fluid Mechanics **72**, 499 (1975).

[90] B. Noetinger, Fluctuating hydrodynamics and Brownian motion, Physica A: Statistical Mechanics and its Applications **163**, 545 (1990).

[91] J.-N. Roux, Brownian particles at different times scales: a new derivation of the Smoluchowski equation, Physica A: Statistical Mechanics and its Applications **188**, 526 (1992).

[92] S. Delong, F. B. Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, and A. Donev, Brownian dynamics without Green's functions, The Journal of chemical physics **140**, 134110 (2014).

[93] D. L. Ermak and J. A. McCammon, Brownian dynamics with hydrodynamic interactions, The Journal of Chemical Physics **69**, 1352 (1978).

[94] T. Ando, E. Chow, Y. Saad, and J. Skolnick, Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations, The Journal of Chemical Physics **137**, 064106 (2012).

[95] B. Delmotte and E. E. Keaveny, Simulating Brownian suspensions with fluctuating hydrodynamics, The Journal of chemical physics **143**, 244109 (2015).

[96] B. Sprinkle, A. Donev, A. P. S. Bhalla, and N. Patankar, Brownian dynamics of fully confined suspensions of rigid particles without Green's functions, The Journal of Chemical Physics **150**, 164116 (2019), https://doi.org/10.1063/1.5090114.

[97] P. Leishangthem and X. Xu, Thermodynamic effects are essential for surface entrapment of bacteria, Physical Review Letters **132**, 238302 (2024).

[98] M. Fixman, Construction of Langevin forces in the simulation of hydrodynamic interaction, Macromolecules **19**, 1204 (1986).

[99] T. Ando, E. Chow, and J. Skolnick, Dynamic simulation of concentrated macromolecular solutions with screened long-range hydrodynamic interactions: Algorithm and limitations, The Journal of chemical physics **139**, 121922 (2013).

[100] E. Chow and Y. Saad, Preconditioned Krylov Subspace Methods for Sampling Multivariate Gaussian Distributions, SIAM Journal on Scientific Computing **36**, A588 (2014).

[101] E. E. Keaveny, Fluctuating force-coupling method for simulations of colloidal suspensions, Journal of Computational Physics **269**, 61 (2014).

[102] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, and J. Bibette, Microscopic artificial swimmers, Nature **437**, 862 (2005).

[103] B. Jang, E. Gutman, N. Stucki, B. F. Seitz, P. D. Wendel-García, T. Newton, J. Pokki, O. Ergeneman, S. Pané, Y. Or, et al., Undulatory locomotion of magnetic multilink nanoswimmers, Nano letters **15**, 4829 (2015).

[104] P. Liao, L. Xing, S. Zhang, and D. Sun, Magnetically driven undulatory microswimmers integrating multiple rigid segments, Small **15**, 1901197 (2019).

[105] H. C. Berg and R. A. Anderson, Bacteria Swim by Rotating their Flagellar Filaments, Nature **245**, 380 (1973), bandiera_abtest: a Cg_type: Nature Research Journals Number: 5425 Primary_atype: Research Publisher: Nature Publishing Group.

[106] S. Trachtenberg and I. Hammel, The rigidity of bacterial flagellar filaments and its relation to filament polymorphism, Journal of structural biology **109**, 18 (1992).

[107] O. F. Müller, *Kleine Schriften aus der Naturhistorie*, Vol. 1 (Buchhandlung der Gelehrten, 1782).

[108] M. R. M. Kapinga and R. Gordon, Cell attachment in the motile colonial diatom bacillaria paxillifer, Diatom Research **7**, 215 (1992), publisher: Taylor & Francis.

[109] R. Gordon, Partial synchronization of the colonial diatom Bacillaria "paradoxa", Research Ideas and Outcomes **2**, e7869 (2016), publisher: Pensoft Publishers.

[110] R. Featherstone, *Robot Dynamics Algorithms, Robotics: Vision, Manipulation and Sensors*, Vol. 22 (Springer US, 1987).

[111] S. F. Schoeller, A. K. Townsend, T. A. Westwood, and E. E. Keaveny, Methods for suspensions of passive and active filaments, Journal of Computational Physics **424**, 109846 (2021).

[112] H. Shum, E. A. Gaffney, and D. J. Smith, Modelling bacterial behaviour close to a no-slip plane boundary: the influence of bacterial geometry, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences **466**, 1725 (2010).

[113] H. Shum and J. M. Yeomans, Entrainment and scattering in microswimmer-colloid interactions, Physical Review Fluids **2**, 113101 (2017).

[114] B. J. Walker, R. J. Wheeler, K. Ishimoto, and E. A. Gaffney, Boundary behaviours of Leishmania mexicana: A hydrodynamic simulation study, Journal of Theoretical Biology **462**, 311 (2019).

[115] J. L. Moran and J. D. Posner, Phoretic self-propulsion, Annu. Rev. Fluid Mech. **49**, 511 (2017).
[116] H. Stark, Artificial chemotaxis of self-phoretic active colloids: collective behavior, Accounts of chemical research **51**, 2681 (2018).

[117] P. Illien, R. Golestanian, and A. Sen, 'Fuelled' motion: phoretic motility and collective behaviour of active colloids, Chemical Society Reviews **46**, 5508 (2017).

[118] A. Domínguez and M. N. Popescu, A fresh view on phoresis and self-phoresis, Current Opinion in Colloid & Interface Science, 101610 (2022).

[119] A. Zöttl and H. Stark, Modeling active colloids: From active brownian particles to hydrodynamic and chemical fields, Annual Review of Condensed Matter Physics **14**, 109 (2023).

[120] S.-Y. Lu, Diffusion and reaction in regular arrays of spheres, The Journal of Chemical Physics **109**, 4985 (1998), https://doi.org/10.1063/1.477110.

[121] S. Michelin and E. Lauga, Phoretic self-propulsion at finite Péclet numbers, J. Fluid Mech. **747**, 572–604 (2014).

[122] A. P. S. Bhalla, B. E. Griffith, N. A. Patankar, and A. Donev, A minimally-resolved immersed boundary model for reaction-diffusion problems, The Journal of Chemical Physics **139**, 214112 (2013).

[123] S. Marbach and L. Bocquet, Osmosis, from molecular insights to large-scale applications, Chemical Society Reviews **48**, 3102 (2019).

[124] G. Wang and J. W. Swan, Surface heterogeneity affects percolation and gelation of colloids: dynamic simulations with random patchy spheres, Soft Matter **15**, 5094 (2019).

[125] G. Vizsnyiczai, G. Frangipane, S. Bianchi, F. Saglimbeni, D. Dell'Arciprete, and R. Di Leonardo, A transition to stable one-dimensional swimming enhances e. coli motility through narrow channels, Nature communications **11**, 2340 (2020).

[126] E. Wajnryb, K. A. Mizerski, P. J. Zuk, and P. Szymczak, Generalization of the Rotne–Prager–Yamakawa mobility and shear disturbance tensors, Journal of Fluid Mechanics **731**, R3 (2013).

[127] A. Broms, M. Sandberg, and A.-K. Tornberg, A locally corrected multiblob method with hydrodynamically matched grids for the Stokes mobility problem, Journal of Computational Physics **487**, 112172 (2023).

[128] P. J. Zuk, E. Wajnryb, K. A. Mizerski, and P. Szymczak, Rotne–Prager–Yamakawa approximation for different-sized particles in application to macromolecular bead models, Journal of Fluid Mechanics **741**, 10.1017/jfm.2013.668 (2014).

[129] B. Delmotte, Viscosity ratio across interfaces controls the stability and self-assembly of microrollers, Physical Review Fluids **8**, L062302 (2023).

[130] J. W. Swan and J. F. Brady, Simulation of hydrodynamically interacting particles near a no-slip boundary, Physics of Fluids (1994-present) **19**, 113306 (2007).

[131] W. Yan and M. Shelley, Universal image system for non-periodic and periodic stokes flows above a no-slip wall, Journal of Computational Physics **375**, 263 (2018).

[132] W. Yan and R. Blackwell, Kernel aggregated fast multipole method: Efficient summation of laplace and stokes kernel functions, Advances in Computational Mathematics **47**, 69 (2021).

[133] C. Aponte-Rivera and R. N. Zia, Simulation of hydrodynamically interacting particles confined by a spherical cavity, Physical Review Fluids **1**, 023301 (2016).

[134] M. Fixman, Simulation of polymer dynamics. I. General theory, The Journal of Chemical Physics **69**, 1527 (1978).

[135] F. B. Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, and A. Donev, Inertial coupling method for particles in an incompressible fluctuating fluid, Computer Methods in Applied Mechanics and Engineering **269**, 139 (2014).

[136] X. Liu and E. Chow, Large-scale hydrodynamic Brownian simulations on multicore and manycore architectures, in *2014 IEEE 28th International Parallel and Distributed Processing Symposium* (IEEE, 2014) pp. 563–572.

[137] R. Singh and R. Adhikari, Pystokes: phoresis and Stokesian hydrodynamics in python, Journal of Open Source Software **5**, 2318 (2020).

[138] A. K. Townsend, Stokesian dynamics in python, Journal of Open Source Software **9**, 6011 (2024).

[139] R. P. Peláez, P. Ibáñez-Freire, P. Palacios-Alonso, A. Donev, and R. Delgado-Buscalioni, Universally adaptable multiscale molecular dynamics (UAMMD), a native-gpu software ecosystem for complex fluids, soft matter, and beyond, Computer Physics Communications **306**, 109363 (2025).

[140] K. W. Torre, R. D. Schram, and J. de Graaf, Python-jax-based fast Stokesian dynamics, arXiv preprint arXiv:2503.07847 (2025).

[141] Z. Liang, Z. Gimbutas, L. Greengard, J. Huang, and S. Jiang, A fast multipole method for the Rotne–Prager–Yamakawa tensor and its applications, Journal of Computational Physics **234**, 133 (2013).

[142] A. M. Fiore, F. Balboa Usabiaga, A. Donev, and J. W. Swan, Rapid sampling of stochastic displacements in Brownian dynamics simulations, The Journal of Chemical Physics **146**, 124116 (2017), http://aip.scitation.org/doi/pdf/10.1063/1.4978242.

[143] A. M. Fiore and J. W. Swan, Rapid sampling of stochastic displacements in Brownian dynamics simulations with stresslet constraints, The Journal of chemical physics **148**, 044114 (2018).

[144] H. Su and E. E. Keaveny, Accelerating the force-coupling method for hydrodynamic interactions in periodic domains, Journal of Computational Physics **510**, 113060 (2024).

[145] P. J. Atzberger, P. R. Kramer, and C. S. Peskin, A stochastic immersed boundary method for fluid-structure dynamics at microscopic length scales, Journal of Computational Physics **224**, 1255 (2007).

[146] P. J. Atzberger, Stochastic eulerian lagrangian methods for fluid-structure interactions with thermal fluctuations, Journal of Computational Physics **230**, 2821 (2011).

[147] P. Plunkett, J. Hu, C. Siefert, and P. J. Atzberger, Spatially adaptive stochastic methods for fluid-structure interactions subject to thermal fluctuations in domains with complex geometries, Journal of Computational Physics **277**, 121 (2014).
[148] A. Hashemi, R. P. Peláez, S. Natesh, B. Sprinkle, O. Maxian, Z. Gan, and A. Donev, Computing hydrodynamic interactions in confined doubly periodic geometries in linear time, The Journal of Chemical Physics 158 (2023).

[149] B. Kallemov, A. Bhalla, B. Griffith, and A. Donev, An immersed boundary method for rigid bodies, Communications in Applied Mathematics and Computational Science 11, 79 (2016).

[150] It is also necessary that the discretization of the Stokes equations preserves some of the original symmetries, specifically $L^{-1}DD^T L^{-1} = L^{-1}$, which is true in the continuum setting but it is not respected by all discretizations [92].

[151] B. Delmotte, E. E. Keaveny, F. Plouraboue, and E. Climent, Large-scale simulation of steady and time-dependent active suspensions with the force-coupling method, Journal of Computational Physics 302, 524 (2015).

[152] J. Higdon, The hydrodynamics of flagellar propulsion: helical waves, Journal of Fluid Mechanics 94, 331 (1979).

[153] H. Goldstein, C. Poole, and J. Safko, *Classical mechanics (3rd ed.)* (Addison-Wesley, 2001).

[154] B. Delmotte, E. Climent, and F. Plouraboué, A general formulation of bead models applied to flexible fibers and active filaments at low reynolds number, Journal of Computational Physics 286, 14 (2015).

[155] W. T. Funkenbusch, K. S. Silmore, and J. W. Swan, Approaches for fast Brownian dynamics simulation with constraints, Journal of Computational Physics 509, 113043 (2024).

[156] F. Rojas-Pérez, B. Delmotte, and S. Michelin, Hydrochemical interactions of phoretic particles: a regularized multipole framework, Journal of Fluid Mechanics 919, A22 (2021).

[157] T. D. Montenegro-Johnson, S. Michelin, and E. Lauga, A regularised singularity approach to phoretic problems, The European Physical Journal E 38, 139 (2015).

[158] W. Yan and R. Blackwell, Kernel aggregated fast multipole method: Efficient summation of laplace and Stokes kernel functions, arXiv preprint arXiv:2010.15155 (2020).

[159] B. Delmotte, A. Donev, M. Driscoll, and P. Chaikin, Minimal model for a hydrodynamic fingering instability in microroller suspensions, Physical Review Fluids 2, 114301 (2017).

[160] F. Balboa Usabiaga, B. Delmotte, and A. Donev, Brownian dynamics of confined suspensions of active microrollers, The Journal of Chemical Physics 146, 134104 (2017), software available at https://github.com/stochasticHydroTools/RigidMultiblobsWall, http://dx.doi.org/10.1063/1.4979494.

[161] Y. Zhang, D. A. Gregory, Y. Zhang, P. J. Smith, S. J. Ebbens, and X. Zhao, Reactive inkjet printing of functional silk stirrers for enhanced mixing and sensing, Small 15, 1804213 (2019).

[162] R. Jahn and A.-M. M. Schmid, Revision of the brackish-freshwater diatom genus *Bacillaria* Gmelin (Bacillariophyta) with the description of a new variety and two new species, European Journal of Phycology 42, 295 (2007).

[163] A. P. Berke, L. Turner, H. C. Berg, and E. Lauga, Hydrodynamic attraction of swimming microorganisms by surfaces, Phys. Rev. Lett. 101, 038102 (2008).

[164] W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St. Angelo, Y. Cao, T. E. Mallouk, P. E. Lammert, and V. H. Crespi, Catalytic nanomotors: Autonomous movement of striped nanorods, Journal of the American Chemical Society 126, 13424 (2004).

[165] J. L. Moran and J. D. Posner, Electrokinetic locomotion due to reaction-induced charge auto-electrophoresis, J. Fluid Mech. 680, 31 (2011).

[166] Y. Zhang, D. A. Gregory, Y. Zhang, P. J. Smith, S. J. Ebbens, and X. Zhao, Reactive inkjet printing of functional silk stirrers for enhanced mixing and sensing, Small 15, 1804213 (2019).

[167] P. Kumar, Y. Zhang, S. J. Ebbens, and X. Zhao, 3d inkjet printed self-propelled motors for micro-stirring, Journal of Colloid and Interface Science 623, 96 (2022).

[168] R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions, Journal of Rheology 58, 1693 (2014), https://doi.org/10.1122/1.4890747.

[169] L. Durlofsky, J. F. Brady, and G. Bossis, Dynamic simulation of hydrodynamically interacting particles, Journal of Fluid Mechanics 180, 21 (1987).

[170] C. Neto, D. R. Evans, E. Bonacurso, H.-J. Butt, and V. S. Craig, Boundary slip in newtonian liquids: a review of experimental studies, Reports on progress in physics 68, 2859 (2005).

[171] C. Kamal, S. Gravelle, and L. Botto, Alignment of a flexible platelike particle in shear flow: Effect of surface slip and edges, Physical Review Fluids 6, 084102 (2021).

[172] J. W. Swan and A. S. Khair, On the hydrodynamics of ‘slip-stick’ spheres, Journal of Fluid Mechanics 606, 115 (2008).

[173] D. Moreno-Chaparro, F. B. Usabiaga, N. Moreno, and M. Ellero, Simulating non-brownian suspensions with non-homogeneous navier slip boundary conditions (2025), arXiv:2505.13505 [cond-mat.soft].

[174] O. Maxian, B. Sprinkle, and A. Donev, Bending fluctuations in semiflexible, inextensible, slender filaments in Stokes flow: Toward a spectral discretization, The Journal of Chemical Physics 158, 154114 (2023), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0144242/16843155/154114_1.5.0144242.pdf.

[175] T. A. Westwood and E. E. Keaveny, Coordinated motion of active filaments on spherical surfaces, Physical Review Fluids 6, L121101 (2021).

[176] E. E. Keaveny, Fluctuating force-coupling method for simulations of colloidal suspensions, Journal of Computational Physics 269, 61 (2014).

[177] R. Fish, B. Sprinkle, R. Pérez Peláez, and A. Donev, libmobility: a fast Stokesian dynamics library, in *APS Division of Fluid Dynamics Meeting Abstracts* (2024) pp. A19-006.

[178] A. Broms, A. H. Barnett, and A.-K. Tornberg, Accurate close interactions of stokes spheres using lubrication-adapted image systems, Journal of Computational Physics 523, 113636 (2025).