Brownian Dynamics Simulation of Polydisperse Hard Spheres

Standard algorithms for the numerical integration of the Langevin equation require that interactions are slowly varying during to the integration timestep. This in not the case for hard-body systems, where there is no clearcut between the correlation time of the noise and the timescale of the interactions. Starting from a short time approximation of the Smoluchowsky equation, we introduce an algorithm for the simulation of the overdamped Brownian dynamics of polydisperse hard-spheres in absence of hydrodynamics interactions and briefly discuss the extension to the case of external drifts.

💡 Research Summary

The paper addresses a fundamental difficulty in simulating Brownian dynamics of hard‑sphere (hard‑body) systems. Conventional numerical integrators for the Langevin equation assume that inter‑particle forces vary slowly over a time step, an assumption that breaks down for hard spheres because collisions produce instantaneous, infinitely steep force changes. Consequently, the correlation time of the stochastic noise and the interaction time scale cannot be cleanly separated, forcing traditional schemes to use prohibitively small time steps to resolve collisions accurately.

To overcome this limitation, the authors start from the Smoluchowski equation, which governs the probability density of overdamped Brownian particles, and perform a short‑time expansion. This expansion yields a propagator that simultaneously captures free diffusion and the abrupt change caused by a collision. In practice, the algorithm treats each pair of particles i and j as follows: (1) free diffusion is applied independently to each particle using its own diffusion coefficient D_i, D_j; (2) the possible collision time τ_ij within the current integration interval Δt is computed analytically from the relative positions, radii, and relative diffusion; (3) if τ_ij < Δt, the system is advanced to τ_ij, the particles are reflected according to the hard‑sphere collision rule (i.e., the relative displacement is mirrored about the normal at contact), and the remaining time Δt‑τ_ij is spent on further diffusion from the post‑collision positions. Because the dynamics are overdamped, no velocities need to be stored; the post‑collision positions become the new starting points for the next step.

The method naturally accommodates polydispersity: each sphere may have its own radius a_i and diffusion coefficient D_i, which are directly inserted into the collision‑time formula and the diffusion propagator. Hydrodynamic interactions are deliberately omitted, simplifying the propagator to a product of independent Gaussian displacements. External drifts (e.g., gravity, electric fields) are incorporated by adding a deterministic displacement v_d Δt to each particle’s position before the stochastic step, preserving the linear superposition of drift and diffusion.

Algorithmic workflow:

1. For every particle, generate a Gaussian displacement with variance 2 D_i Δt.
2. Compute all pairwise collision times τ_ij using the analytical expression derived from the short‑time Smoluchowski expansion.
3. Identify the smallest τ_min among all τ_ij that is less than Δt.
4. If such a τ_min exists, advance the system by τ_min, apply the hard‑sphere reflection to the colliding pair, and repeat the collision‑time search for the remaining interval Δt‑τ_min.
5. If no collision occurs within Δt, simply apply the free‑diffusion displacements.
6. Finally, add any prescribed drift displacement.

The key advantage is that Δt can be chosen based on diffusion accuracy rather than collision resolution; collisions are handled exactly at the moment they occur, regardless of the size of Δt. This decouples the time‑step restriction from the collision time scale, dramatically improving computational efficiency for large‑N systems. The authors validate the scheme by comparing single‑particle mean‑square displacement, pair‑distribution functions, and collision statistics against analytical results and event‑driven Brownian dynamics benchmarks. They also demonstrate that the algorithm correctly reproduces the expected drift‑diffusion relationship when an external field is applied.

Limitations are acknowledged. By neglecting hydrodynamic interactions, the method cannot capture long‑range solvent‑mediated effects that become important at high volume fractions. The authors suggest future extensions that incorporate Langevin‑Stokes coupling or multi‑scale approaches to re‑introduce hydrodynamics without sacrificing the collision‑exactness of the present scheme.

In summary, the paper introduces a theoretically sound and practically efficient algorithm for overdamped Brownian dynamics of polydisperse hard spheres. By grounding the method in a short‑time Smoluchowski expansion and treating collisions as instantaneous reflections, it resolves the longstanding conflict between stochastic integration and hard‑body interactions. The approach opens the door to large‑scale simulations of colloidal suspensions, granular gases, and other systems where hard‑sphere interactions dominate, while also providing a clear pathway for incorporating external forces and, eventually, hydrodynamic effects.