Long-time Brownian dynamics of hard and sticky colloids by geometric projection

For colloidal suspensions, Brownian dynamics (BD) provides a particle-scale description of how microscopic interactions drive collective phenomena such as aggregation, gelation, sedimentation and deposition. However, the accessible physical duration of these simulations is often set by the fastest time scale retained in the numerical model, leaving many slower processes of interest beyond the reach of a direct particle-resolved simulation. Brownian dynamics already removes one such time scale by eliminating inertial relaxation and retaining only the overdamped configurational dynamics [1]. Nevertheless, this does not eliminate all microscopic stiffness: explicit near-contact potentials can introduce short configurational relaxation times. A steep repulsive wall and a narrow attractive well differ in the shape of the interaction potential, but both generate rapid positional relaxation that must be resolved by the integrator. For collective observables on time scales much longer than near-contact relaxations, such as aggregation kinetics or structural correlations, the detailed relaxation pathway is often irrelevant: particles cannot overlap, and irreversibly adhesive contacts remain bonded. The corresponding limiting description is therefore naturally reformulated in terms of geometric constraints rather than time-resolved near-contact forces [2, 3]. For hard spheres, this idea underlies several potential-free Brownian schemes, from event-driven [4] and rejection-based [5] methods to overlap-correction [6, 7] and reflected-transition-probability [8] approaches. At high density, however, their efficiency can be challenged by frequent encounters and many-body overlap events. This difficulty is further amplified by irreversible aggregation: sticky contacts persist over many Brownian steps and, in dense clusters, gels, or deposits, form strongly coupled constraint networks. Here we take a complementary route centered on a simple, well-defined geometric prescription designed to remain flexible and robust in dense, strongly coupled constraint networks. We formulate hard exclusion, sticky contacts, and wall interactions within a common constraint-projection framework for overdamped BD. The goal is to remove irrelevant fast modes associated with near-contact interactions while retaining longer-range smooth forces, Brownian diffusion, external driving, and the dense hard or sticky contact networks that control long-time aggregation and deposition pathways.

Once stiff near-contact interactions are replaced by hard or sticky contact rules, their role in the dynamics shifts from providing forces to restricting the admissible domain of configuration space. We exploit this by enforcing these geometric constraints on the configuration obtained after each free BD step, borrowing from the computer-graphics framework of position-based dynamics [9, 10], but integrated within an overdamped stochastic integrator. The instantaneous configuration of a system is described by the set of coordinates $\mathbf{X}$. Although the method can be applied to various types of degrees of freedom, such as rigid-body transformations, in this work we restrict ourselves to the Cartesian coordinates of the centers of $N$ spherical particles, $\mathbf{X} = \{\mathbf{r}_i\}_{i=1}^{N}$. The admissible configurations are expressed through a set of $N_{\mathrm{cons}}$ coupled constraint residuals, $\mathcal{C} = \{C_\alpha(\mathbf{X})\}_{\alpha=1}^{N_{\mathrm{cons}}}$, with $C_\alpha(\mathbf{X}) = 0$ for all satisfied constraints. The interaction of a particle $i$ of radius $a_i$ with a lower hard wall is modeled by

where $z_i^0 = z^{\mathrm{wall}} + a_i + h^0$ is the lowest $z$-coordinate allowed for the particle, with a minimum gap $h^0$ imposed on all the contacts. We model the bonded interaction of particles $i$ and $j$ by

where $d_{ij}^0 = a_i+a_j+h^0$, and $r_{ij} = \lVert\mathbf{r}_j-\mathbf{r}_i\rVert$. Similarly, for non-bonded particles $i$ and $j$, we model their hard-sphere interaction by

With the discontinuous interactions encoded in $\mathcal{C}$, time integration alternates between a free BD update under smooth forces and a geometric projection. The free update yields an unconstrained configuration $\mathbf{X}^{(0)}$, from which the projection solver constructs a sequence $\mathbf{X}^{(k)}$ until the maximum absolute residual falls below the tolerance $\varepsilon$. An iteration of the solver updates the configuration according to $\mathbf{X}^{(k+1)} = \mathbf{X}^{(k)} + \omega\Delta\mathbf{X}^{(k)}$, where $\omega\le1$ is a relaxation factor and the unrelaxed correction is given by

where $\mathbf{M}$ is the mobility matrix. In the present Cartesian representation, hydrodynamic interactions are neglected and $\mathbf{M}$ is therefore a $3N \times3N$ diagonal matrix, with each coordinate belonging to particle $i$ assigned the Stokes mobility $\mu_i = 1/\gamma_i$. All pair-constraint corrections are evaluated from $\mathbf{X}^{(k)}$, accumulated independently, and applied simultaneously. This simultaneous update is naturally parallel, but can lead to persistent correction oscillations in dense bonded systems; under-relaxation with $\omega< 1$ was therefore used to improve convergence robustness.

Adhesive aggregation is introduced as an update of the constraint set: a non-bonded pair is assigned an irreversible bond constraint when $r_{ij} < d_{ij}^0+\varepsilon$. The same $\varepsilon$ is used as the projection tolerance, ensuring consistency between bond creation and bonded-contact projection. Together, the contact distance $h^0$ and the tolerance $\varepsilon$ set the geometric resolution of near-contact dynamics and are chosen to balance stable convergence with physically consistent long-time behavior. The resulting scheme explicitly separates smooth interactions, which generate the free Brownian step, from singular hard and sticky interactions, which define the admissible geometry onto which the configuration is projected. The method is implemented as a LAMMPS fix [11, 12] and is available, with input examples, in a zenodo archive (DOI: 10.5281/zenodo.23802331). Further details of the Brownian update, constraint projection, and bond handling are provided in the ESI†.

We validate the constraint formulation in the two elementary limits on which the method relies: non-aggregating hard-sphere Brownian motion and irreversible sticky aggregation. In the non-aggregating case, the projection should enforce excluded volume without altering either the equilibrium structure or the long-time diffusive dynamics. We therefore compare hard-sphere constrained simulations with explicit-potential BD that uses the repulsive branch of the integrated Lennard-Jones potential [13], truncated at its minimum. The radial distribution function is measured at $\phi= 0.55$, where steric structure is strongly constrained, while the time-dependent diffusion coefficient $D(t)$ is measured at $\phi= 0.25$, where the long-time diffusive regime remains accessible to explicit-potential BD. As shown in Figure 1, the projected dynamics preserves the tested equilibrium structure and diffusive dynamics within the tested time-step range. We then test the adhesive contact rule at $\phi= 0.25$ by comparing the normalized cluster number $N_{\mathrm{c}}(t)/N$ to explicit BD simulations with attractive integrated Lennard-Jones potentials of different ranges and well shapes, potential 1, 2 and 3, shown in Figure 2. These explicit potentials retain the near-contact attraction but still require the small time step associated with fast bonded-pair relaxation. As the attraction becomes short-ranged but still effectively irreversible, the explicit-potential kinetics approaches the sticky-contact result, showing that the constraint rule captures the contact-limited aggregation dynamics targeted by the coarse-grained model. These comparisons were also used to select the largest time step for which the relevant observables remained unchanged, thereby defining the validated time steps used in the performance measurements below. Time-step sensitivity tests and complete simulation parameters are reported in the ESI†.

The validated time steps already reveal the origin of the computational gain: explicit-potential BD required $\Delta t = 50\ \mathrm{ps}$ for non-aggregating suspensions and $\Delta t = 25\ \mathrm{ps}$ for aggregating ones, whereas the constraint-based simulations reached $\Delta t = 50\ \mathrm{ns}$ for both non-aggregating hard spheres and sticky hard spheres. Be-

cause this three-order-of-magnitude increase in physical time per step comes with the additional cost of completing the projection, we assess the net speed increase through physical time throughput $t_{\mathrm{phys}}/t_{\mathrm{wall}}$ where $t_{\mathrm{phys}}$ is the simulated physical time and $t_{\mathrm{wall}}$ is the corresponding wall-clock time. fig:3a shows a throughput one to two orders of magnitude larger than explicit-potential BD over the tested volume fractions, despite the increasing projection cost at high density and in bonded systems. The solver was deliberately chosen for its simplicity and parallel implementation: constraint corrections are evaluated independently before being applied simultaneously. This design is naturally parallel but can require many iterations to converge, especially in dense bonded networks. Thus, even with this simple projection scheme, the gain obtained by removing the stiff near-contact time scale is not lost to the solver. This computational margin is what makes it possible to address low-field deposition regimes that remain impractical with explicit potentials. Benchmark definitions, simulation parameters, and parallel-scaling results are provided in the ESI$^\dagger$.

To demonstrate the new dynamical regimes made accessible by the method, we apply it to electrophoretic deposition (EPD), a coating technique in which charged particles are driven by an electric field toward an electrode, where densification and irreversible coagulation jointly lead to deposit formation[14, 15]. Despite extensive study, direct experimental access to particle-scale dynamics during EPD remains limited, particularly in concentrated suspensions and growing deposits, which complicates mechanistic interpretation and makes particle-resolved simulations valuable for linking prescribed effective interactions to deposition pathways[16]. Lowering the electric field directly lowers the Péclet number,

and increases the drift time over a distance $L$,

Therefore, at fixed drift-scaled deposition time, lowering Pe increases the physical duration to be simulated, while explicit DLVO Brownian dynamics remains constrained by the short time step imposed by the stiff near-contact potential. A finite wall-clock budget therefore sets a lower bound on the Péclet numbers that can be reached. In a favorable setup, with $N = 10{,}000$ particles, initial volume fraction $\phi= 0.1$ and box length $L_z = 167a$ along the field direction, explicit DLVO simulations would already require several weeks to reach $\mathrm{Pe} = 10$ and would remain impractical at $\mathrm{Pe} = 1$. In contrast, the present constraint-based method reaches $\mathrm{Pe} = 10$ in about 1.5 days and $\mathrm{Pe} = 1$ in about 6 days on the same computational setup, enabling the comparison of the deposition pathways shown in Figure 4. The effective interaction, simulation setup, and time-step validation for the EPD calculations are detailed in the ESI$^\dagger$.

To characterize the local formation of the bonded network, we measure the surface density $n(z)$ of bonds crossing a plane at height $z$, and their orientational order parameter $S(z)$ along the deposition direction,

where $u_{i,z}$ is the $z$-component of the bond unit vector, and $\langle i,j\rangle_z$ denotes bonds intersecting the plane at height $z$. Because this plane-crossing selection weights bonds by their probability of intersecting the plane, an isotropic bond network gives $S^{\mathrm{iso}} = 1/4$. Details of the density, bond-orientation, and pore-size analyses are provided in the ESI$^\dagger$. The accessible low-Pe regime already displays signatures that are captured by a simple extrapolation from high-field deposition. At high Pe, particles are rapidly driven into the growing front: at a given height, the local volume fraction and the surface density of crossing bonds rise sharply, and the first bonds are preferentially oriented along the field direction. Accordingly, $S(t)$, measured at $z = 15a$, starts above the isotropic reference before relaxing toward it as the deposit densifies and rearranges. At $\mathrm{Pe} = 1$, by contrast, the local density builds up more progressively and bond formation starts close to the isotropic reference, consistent with barrier-limited coagulation occurring upon densification rather than immediate front-induced bonding. This change in deposition pathway also leaves a structural signature in the final deposit: the core pore-size distribution shifts from small interstitial voids at high Pe toward a larger population of particle-scale pores at Pe = 1. These dynamical and structural signatures show that the method opens the way to systematic particle-resolved studies of low-field EPD regimes that were previously impractical with explicit potentials.

Overall, these results show that stiff near-contact interactions can be removed from time integration when their relevant effect can be represented by geometric constraints and contact-capture rules. The constraint-based Brownian dynamics scheme preserves the targeted hard-sphere and contact-limited aggregation observables while extending the accessible physical time, as illustrated here by low-Pe EPD. Beyond the spherical hard and sticky contacts considered in this work, the formulation is modular: additional degrees of freedom and constraint functions, for example for non-spherical particles, can in principle be introduced within the same projection framework. Extensions to reversible contacts or configuration-dependent hydrodynamic mobilities would, however, require additional ingredients beyond the present formulation. More broadly, the approach provides a practical way to align the numerical time scales of Brownian dynamics simulations with the collective time scales of interest.

Author Contributions

Rémi Martin: Methodology, Software, Validation, Investigation, Formal analysis, Writing – Original Draft, Writing – Review & Editing. Sandrine Duluard: Methodology, Writing – Review & Editing, Supervision, Resources, Funding acquisition. Céline Merlet: Conceptualization, Formal analysis, Methodology, Validation, Writing – Review & Editing, Supervision, Resources, Funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The Zenodo archive (DOI: 10.5281/zenodo.23802331) provides the LAMMPS implementation of the geometric projection method together with example input files for the simulations.

Acknowledgment

The authors thank the French «Ministère de l’Enseignement Supérieur et de la Recherche et de l’Espace (MESRE)» for the PhD grant of Rémi Martin. This work was granted access to the HPC resources of the CALMIP supercomputing centre under the allocation P21014. The authors acknowledge Florence Ansart, Pierre-Louis Taberna, Vincent Dahirel and Marie Jardat for useful discussions.

Notes and references