Active Fast Stokesian Dynamics

• Active Fast Stokesian Dynamics is an extension of classical Stokesian Dynamics that incorporates activity via prescribed slip, forces, and stresslets.
• It employs fast, matrix-free solvers with Ewald or FMM summation and pairwise lubrication to achieve efficient O(N) scaling for many-body simulations.
• The methodology preserves hydrodynamic fidelity while modeling complex active suspensions including squirmers, shakers, and odd-viscosity effects.

Active Fast Stokesian Dynamics is the extension of Stokesian Dynamics to suspensions of active particles together with the algorithmic machinery that makes many-body simulations tractable at large $N$. In this setting, the low-Reynolds-number hydrodynamics remain those of classical Stokes flow, but activity enters through prescribed boundary velocities, active forces, torques, stresslets, or slip moments, while fast matrix-free solvers, Ewald- or FMM-type summation, and pairwise lubrication retain the characteristic balance of hydrodynamic fidelity and computational efficiency associated with modern Fast Stokesian Dynamics (Elfring et al., 2021, Torre et al., 10 Mar 2025). The term is therefore used across a family of formulations: thermalized and athermal, periodic and open, particle- and aggregate-resolved, and with implementations ranging from JAX-based GPU codes to boundary-integral and HOOMD-integrated solvers (Ge et al., 2023, Yan et al., 2019).

1. Hydrodynamic and operator-theoretic foundations

Active Fast Stokesian Dynamics inherits the standard Stokesian assumptions: an incompressible Newtonian fluid, negligible inertia, and rigid particles whose surface velocity is the sum of rigid-body motion and an active contribution. In the notation used for active suspensions, the boundary condition on particle $i$ is

$$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$

with $\mathbf{u}^s_i$ the prescribed active slip or deformation velocity. The governing equations are the Stokes equations,

$$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$

or equivalently $\nabla p=\eta \nabla^2 \mathbf{u}$, $\nabla\cdot \mathbf{u}=0$, depending on notation (Elfring et al., 2021, Torre et al., 10 Mar 2025).

The central structural point is that activity changes the inhomogeneous terms but not the geometry-dependent grand mobility and grand resistance operators. In resistance form, one widely used SD expression writes

$$\begin{pmatrix}\mathcal{F}^H\ \mathcal{S}^H\end{pmatrix} = -\mathcal{R}\cdot \begin{pmatrix} (\mathcal{U}-\mathcal{U}^\infty)\ -\mathcal{E}^\infty \end{pmatrix},$$

with $\mathcal{U}$ the generalized rigid-body velocities and $\mathcal{S}^H$ the stresslets. The classical SD near-/far-field split is

$i$0

where $i$1 is the far-field grand mobility, $i$2 are pairwise lubrication resistances, and $i$3 removes double counting of pairwise far-field terms (Torre et al., 10 Mar 2025).

For active particles, the reciprocal-theorem formulation makes the same point in an exact way. The rigid-body dynamics can be written as

$i$4

or, with Brownian motion,

$i$5

so the passive resistance operator is reused while the active slip contributes through $i$6 and higher irreducible moments (Elfring et al., 2021).

This invariance of operator structure is the reason “active” and “fast” can be combined without changing the core numerical architecture. Fast solvers accelerate the same mobility or resistance action that appears in passive SD; activity modifies the right-hand side, the constraint terms, or the active velocity moments.

2. Fast formulations and computational scaling

The defining numerical move in modern Fast Stokesian Dynamics is to avoid explicit assembly and inversion of dense many-body mobility matrices. In the Fiore–Swan lineage adopted by JFSD, the SD update is recast as a linear saddle-point problem coupling generalized velocities, traction-like Lagrange multipliers, and near-field lubrication variables: $i$7 with overall size $i$8. The Schur complement reduces the task to repeated applications of the configuration-dependent mobility operator, so dense matrices are never formed or inverted. Because the near-field block uses neighbor lists and the far-field action uses fast summation, each block application can be performed in $i$9; with bounded Krylov iteration counts and suitable preconditioning, the overall step scales linearly with particle number (Torre et al., 10 Mar 2025).

The same fast philosophy appears in other SD accelerations. Spectral Ewald Accelerated Stokesian Dynamics evaluates the periodic far-field mobility with spectral accuracy and $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$0 complexity while retaining SD lubrication corrections, using GPU-accelerated spreading, FFTs, wave-space multiplication, and gathering. Its importance for active problems is explicit: when adding active forces, torques, or surface slips, the far-field mobility evaluation is unchanged and only the right-hand side of the linear system is modified (Wang et al., 2015).

A distinct but closely related route is a mobility-agnostic boundary-integral platform in which a second-kind boundary integral equation yields a symmetric positive-definite mobility, KIFMM provides $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$1 far-field evaluation, vector spherical harmonics resolve lubrication-dominated near-singular interactions, and non-overlap is enforced by a complementarity/QP contact solve. That framework is compatible with RPY, SD, and BI mobilities, and demonstrates weak and strong parallel scaling up to $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$2 particles with approximately $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$3 degrees of freedom on $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$4 cores (Yan et al., 2019).

These formulations differ in discretization and solver technology, but they share the same computational doctrine: matrix-free operator application, iterative Krylov solution, explicit or spectrally accurate treatment of near-field lubrication, and fast far-field summation. In that sense, Active Fast Stokesian Dynamics is best understood as an operator-centered methodology rather than as a single codebase.

3. Representations of activity

Within this framework, activity can be introduced in several algebraically equivalent ways. In JFSD, active forces and torques are added directly to the applied generalized force: $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$5 so the active extension is explicitly “parameteric”: one injects activity through the same blocks that already accept external forces, torques, stresslets, or background flows (Torre et al., 10 Mar 2025).

For squirmers and related swimmers, activity is more naturally represented by prescribed surface slip and its irreducible moments. In the irreducible-moment notation of Active Stokesian Dynamics,

$$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$6

so the first squirming mode sets self-propulsion, the second sets the active stresslet, and the geometry-dependent mobility blocks translate those moments into many-body motion (Elfring et al., 2021). In the AFSD liquid-film formulation for squirmers, the same information is written as a generalized swim velocity $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$7 and a symmetric slip moment $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$8 (Zhou et al., 30 Sep 2025).

Apolar active suspensions constitute a particularly important special case. In the shaker systems studied with AFSD, $$\mathbf{u}(\mathbf{x}\in\partial \mathcal{B}_i)=\mathbf{U}_i+\boldsymbol{\Omega}_i \times \mathbf{r}_i+\mathbf{u}^s_i(\mathbf{x}),$$9 and $\mathbf{u}^s_i$0, so there is no self-propulsive translation and activity is entirely stresslet-like; pullers have $\mathbf{u}^s_i$1 and pushers have $\mathbf{u}^s_i$2 (Ge et al., 2023). This representation is minimal but not trivial: it preserves the long-range $\mathbf{u}^s_i$3 stresslet field, near-field lubrication, and the many-body couplings between translation, rotation, and stresslet-induced flow.

A recurrent misconception is that “neutral” or “stealth” squirmers are hydrodynamically silent in suspension. AFSD studies of MIPS in liquid films show that this is false once collisions are present: even for neutral squirmers with no solitary disturbance flow, conservative steric collisions generate pusher force dipoles, and these collision-induced dipoles dominate the long-range hydrodynamic interactions in suspension (Zhou et al., 30 Sep 2025).

The activity concept can also be shifted from particles to the suspending medium. “Stokesian Dynamics with odd viscosity” generalizes the Green’s functions and mobility blocks of SD to quasi-2D active fluids with broken time-reversal symmetry, producing antisymmetric odd-mobility couplings and a chiral Oseen tensor. Fast Ewald and FMM acceleration are described as compatible in principle, but a rigorous odd lubrication theory is not yet available (Yuan et al., 2022).

4. Implementations and software realizations

Implementation	Hydrodynamic core	Distinctive scope
JFSD (Torre et al., 10 Mar 2025)	Python/JAX FSD saddle-point solver	JIT/XLA, GPU acceleration, open and periodic BCs
AFSD for squirmers (Ge et al., 28 Feb 2026)	Fast periodic summation with lubrication in HOOMD-blue	Sheared active suspensions up to $\mathbf{u}^s_i$4
Mobility-agnostic BI platform (Yan et al., 2019)	BI mobility, KIFMM, VSH, BBPGD contact solve	Dense rigid suspensions, torques or slip, $\mathbf{u}^s_i$5 particles
SHAPES (Gissinger et al., 11 Jun 2026)	Bead-level SD for aggregate tensors	Walls, flexible aggregates, active swimmers

JFSD is the clearest example of a general-purpose, portable AFSD implementation. It realizes the Fiore–Swan FSD method in a Python/JAX stack with JIT-compiled kernels and GPU acceleration, supports “open” and periodic hydrodynamics through scalar mobility functions, uses neighbor lists derived from JAX-MD, and provides a modular TOML-based configuration layer. It was tested on a mid-range NVIDIA GeForce RTX 4060 Ti (8 GB); the initial JIT compilation takes $\mathbf{u}^s_i$6–$\mathbf{u}^s_i$7 minutes, and the code can simulate on the order of $\mathbf{u}^s_i$8–$\mathbf{u}^s_i$9 particles on a desktop GPU within memory constraints (Torre et al., 10 Mar 2025).

The HOOMD-based AFSD used for sheared squirmers adopts the same far-/near-field decomposition and Krylov solution strategy but specializes to periodic 3D shear, athermal suspensions, and squirmer activity through $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$0 and $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$1. In practice it reaches $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$2 and resolves dissipation, stress, normal stress differences, and microstructure under simple shear (Ge et al., 28 Feb 2026).

The contact-aware boundary-integral platform occupies a complementary niche. It is not framed as FSD in the narrow Fiore–Swan sense, but it realizes large-scale active Stokes suspensions with matrix-free mobility solves, exact non-overlap constraints, and full contact-network reconstruction. Because the collision-resolution layer is mobility-agnostic, SD, RPY, and BI back ends can all be used (Yan et al., 2019).

SHAPES extends the SD viewpoint to rigid and flexible aggregates. It computes aggregate resistance and mobility tensors, including wall effects, from bead-level hydrodynamics and prescribed internal kinematics, and the paper explicitly links those tensor blocks to future Fast Stokesian Dynamics simulations of many interacting active aggregates treated as superparticles (Gissinger et al., 11 Jun 2026).

5. Validation, benchmark structure, and physical results

The passive and thermal benchmark base for AFSD is unusually extensive in JFSD. Physical validations include sedimenting three spheres, sedimentation velocity $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$3 for a simple cubic array of $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$4 spheres versus $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$5, pair dynamics in simple shear, shear viscosity $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$6 for $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$7, thermal MSD of a single particle in periodic boxes, equilibrium sampling under a linear pair potential $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$8, and short-time self-diffusion $$\mu \nabla^2 \mathbf{u} - \nabla p = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0,$$9 versus $\nabla p=\eta \nabla^2 \mathbf{u}$0 for $\nabla p=\eta \nabla^2 \mathbf{u}$1. These tests are used to verify both scaling and physical correctness, including Hasimoto-periodic corrections and the expected $\nabla p=\eta \nabla^2 \mathbf{u}$2 complexity (Torre et al., 10 Mar 2025).

For active suspensions, one major benchmark class concerns apolar shakers. AFSD simulations of spherical squirmers with $\nabla p=\eta \nabla^2 \mathbf{u}$3 and $\nabla p=\eta \nabla^2 \mathbf{u}$4 over volume fractions from $\nabla p=\eta \nabla^2 \mathbf{u}$5 to $\nabla p=\eta \nabla^2 \mathbf{u}$6 show little difference between pullers and pushers in their instantaneous and long-time dynamics. The translational diffusion coefficient varies non-monotonically with volume fraction, peaking around $\nabla p=\eta \nabla^2 \mathbf{u}$7–$\nabla p=\eta \nabla^2 \mathbf{u}$8, while rotational relaxation accelerates with $\nabla p=\eta \nabla^2 \mathbf{u}$9. The study attributes these trends to activity-induced hydrodynamic interactions together with weak local order seen in pair statistics, namely $\nabla\cdot \mathbf{u}=0$0 and weak local nematic ordering $\nabla\cdot \mathbf{u}=0$1 (Ge et al., 2023).

A second benchmark class concerns nonequilibrium phase behavior. In quasi-2D liquid films, AFSD simulations of up to $\nabla\cdot \mathbf{u}=0$2 squirmers show that MIPS disappears for stealth swimmers with full hydrodynamic interactions. Rotational HI suppress phase separation by increasing the effective rotational diffusivity and shortening the persistence length, while translational HI in thick films generate confinement-scale advective flows that sweep swimmers away from clusters. The same study distinguishes this from wall-induced aggregation: collisional dipoles between swimmers and boundaries can generate attractive flows and clustering, but that mechanism is not the same as MIPS (Zhou et al., 30 Sep 2025).

A third class concerns rheology and dissipation under flow. AFSD simulations of semi-dilute to concentrated suspensions of shakers in simple shear show that shear enhances the total dissipation but reduces the relative viscosity for both puller- and pusher-type shakers. At low shear rates where activity dominates, pushers dissipate more energy than pullers; at high shear rates, both behave effectively as passive spheres, producing negative normal stress differences due to shear-induced collision. The same simulations report enhanced nematic order and anisotropic pair correlation at intermediate rates, and much weaker shear-thinning for neutral self-propelled squirmers with $\nabla\cdot \mathbf{u}=0$3, indicating that internal activity rather than motility alone is the dominant driver of the strong rate dependence (Ge et al., 28 Feb 2026).

Taken together, these results show that AFSD is not merely a faster implementation of passive SD. It has become a quantitative tool for testing active-matter mechanisms that depend simultaneously on many-body HI, lubrication, steric contact, confinement, and, in some implementations, Brownian fluctuations.

6. Limitations, misconceptions, and current directions

The method’s strengths are coupled to clear limitations. JFSD retains the usual SD sensitivities: small $\nabla\cdot \mathbf{u}=0$4 is needed near contacts due to lubrication, Euler integration introduces $\nabla\cdot \mathbf{u}=0$5 bias in long tails unless $\nabla\cdot \mathbf{u}=0$6 is very small, GPU memory caps $\nabla\cdot \mathbf{u}=0$7, double precision may reduce throughput, and solver iteration counts increase with $\nabla\cdot \mathbf{u}=0$8 even though the normalized step time preserves linear complexity (Torre et al., 10 Mar 2025). In the apolar-squirmer and sheared-squirmer studies, Brownian motion is neglected, the particles are spherical, and activity is truncated to minimal slip moments; these are deliberate modeling choices, not universal features of AFSD (Ge et al., 2023, Ge et al., 28 Feb 2026).

Another important limitation is that different AFSD branches treat contacts differently. The Fiore–Swan/JFSD lineage typically prevents overlap with excluded-volume interactions while resolving lubrication hydrodynamics. The contact-aware BI platform instead enforces non-overlap by a complementarity/QP solve that imposes Newton’s Third Law, records the complete contact network, and permits large timesteps without overlap. These are distinct numerical philosophies rather than interchangeable implementation details (Yan et al., 2019).

Odd-viscosity extensions underscore an additional boundary of current knowledge. The odd-mobility kernels, generalized Faxén laws, and symmetry relations have been derived, and compatibility with Ewald, spectral Ewald, FMM, and Fast SD is described, but near-field odd lubrication theory is not yet available. For dense suspensions in odd media, this remains a substantive open problem rather than a minor implementation gap (Yuan et al., 2022).

Complex shapes pose a similar challenge. SHAPES shows that SD can be adapted to compute complete resistance and mobility tensors for aggregate particles, including confinement by walls and active internal actuation, but it omits lubrication by default for robustness and efficiency and notes that contact-rich or wall-gap-dominated phenomena still require additional modeling (Gissinger et al., 11 Jun 2026).

A further current direction is hybridization with learned surrogates. A 2025 perspective proposes a framework that replaces SD’s full many-body computations with a neural network trained on SD data, correcting the RPY approximation while avoiding costly matrix inversions. The same paper emphasizes that activity can still be represented through prescribed stresslets within the grand mobility/resistance formalism, suggesting a route toward matrix-free $\nabla\cdot \mathbf{u}=0$9 active simulations with SD-level correction structure (Torre et al., 20 Mar 2025).

The broad misconception that AFSD is a single monolithic method is therefore inaccurate. The literature instead supports a more precise view: AFSD is a class of operator-based active-suspension solvers that preserve SD’s many-body hydrodynamics and near-field structure while accelerating the far field through matrix-free fast summation. Its present frontier lies in extending that synthesis to richer activity models, more general shapes, stronger confinement, odd-viscosity media, and data-driven corrections without sacrificing the hydrodynamic structure that gives the method its predictive value.