Pair-distribution function of active Brownian spheres in three spatial dimensions: simulation results and analytical representation

The gist The authors analyze the full pair-distribution function of spherical active Brownian particles in three spatial dimensions using Brownian dynamics simulations and derive an analytical representation parametrized by activity and concentration, which advances quantitative models of active Brownian particles.

Simulation Methodology

The study employs modified LAMMPS software for Brownian dynamics simulations to investigate the dynamics of N spherical ABPs described by overdamped Langevin equations. The translational motion is governed by equation (1), which includes a self-propulsion force and noise term ft,i(t). Rotational motion is given by equation (3), characterized by the rotational diffusion coefficient Dr, which relates to the translational diffusion coefficient Dt. The Péclet number Pe = v0σ/Dt is varied between 50 and 500 by changing temperature, and the average packing density Φ0 is varied from 0.05 to 0.7 by changing the particle number N.

Pair-Distribution Function Analysis

The pair-distribution function g(r1, r2,ub1,ub2, t) is defined as the two-body density derived from the probability distribution P(AB). In a stationary and homogeneous system, this function simplifies to depend only on the Péclet number Pe, packing density Φ0, relative position r = r ubd, and orientations ub1 and ub2. The authors define coordinates such that the pair-distribution function depends on four variables: g = g(r, θ1, θ2, ϕ2). The angles are defined by θ1 = arccos(ub1 · ubd) and θ2 = arccos(ub1 · ub2).

State Diagram and Homogeneity

The analysis involves determining the characteristic length Lc, which quantifies density inhomogeneities and is used to distinguish between a homogeneous system and one undergoing phase separation. The state diagram in Fig. 2 shows that the system remains homogeneous over time if either Φ0 < 0.35 or Pe < 100, indicating that phase separation occurs outside this region. In three dimensions, the packing density of particles in a cluster is lower than the perfect packing density of an fcc grid, suggesting no phase separation based on the characteristic length definition in this case.

Analytical Representation and Fitting

To enable field-theoretical models, an analytical expression for the product function p(Pe, Φ0; r, θ1, θ2, ϕ2) = −U′(r) g(Pe, Φ0; r, θ1, θ2, ϕ2) is developed. This product function p is nonzero only for 0.9σ < r < 2/6σ because g vanishes for r less than approximately σ due to strong repulsion and U′ vanishes for r > 2/6σ due to the interaction force cutoff. The function p is approximated using a Fourier expansion (Eq. 19) and truncated at second order, resulting in an expression involving coefficients αh,j,k and βh,j,1. These coefficients are fitted using an Exponentially Modified Gaussian distribution (EMG) function combined with a linear cutoff term to describe the dependence on r.

Results and Discussion

The pair-distribution function g is visualized for various fixed distances r and angles θ1, θ2, ϕ2 in Figs. 4, 6, and 7. Maxima of g are found for configurations that are very stable or easy to reach, such as configuration B where particles are oriented towards each other. The product function p is shown in Fig. 3, where the mean absolute error MAE/MAV varies between 2 and 55 percent with high errors occurring only for very small packing densities and very high Péclet numbers. The results confirm the state diagram by Stenhammar et al. [27] and allow determination of regions where symmetry assumptions are valid. The work provides an analytical expression for the product of the pair-distribution function and the derivative of the interaction potential that provides an excellent fit to simulation data.

Conclusion

The study confirms and improves state diagrams obtained in previous works [27–29] and extracts a fully orientation-resolved pair-distribution function for homogeneous particle distributions across a wide range of parameters. The results provide insights into the collective dynamics of ABPs in three dimensions and can be exploited in the derivation of active field theories and for obtaining microscopic predictions for state boundaries in active systems. Possible extensions include investigations of mixtures of active and passive particles and particles with more complex shapes. The consistency with two-dimensional cross sections is demonstrated by agreement with Ref. [26].

Supplementary Material

The Supplementary Material contains a spreadsheet with the values of the fit parameters (as shown in Appendix A) needed to recreate the analytical representation of the product function, a Python script abp.spherical3d.pairdistribution that recreates the approximation of the product function p using these values, and Python scripts and raw data needed to recreate Figs. 1–8.

Conflicts of Interest

There are no conflicts of interest to declare.

Acknowledgements

R.W. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – ProjectID 433682494 – SFB 1459 The simulations for this work were performed on the computer cluster PALMA II of the University of Münster.

References

[1] C. Bechinger, R. Di Leonardo, H. L¨owen, C. Reichhardt, G. Volpe, and G. Volpe, Active particles in complex and crowded environments, Reviews of Modern Physics 88, 045006 (2016)

[2] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A Simha, Hydrodynamics of soft active matter, Reviews of Modern Physics 85, 1143 (2013)

[3] K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly, and R. E. Goldstein, Fluid dynamics and noise in bacterial cell–cell and cell–surface scattering, Proceedings of the National Academy of Sciences U.S.A. 108, 10940 (2011)

[4] H. C. Berg, E. coli in Motion (Springer-Verlag, New York, 2008)

[5] A. P. Petroff, X.-L. Wu, and A. Libchaber, Fast-moving bacteria self-organize into active two-dimensional crystals of rotating cells, Physical Review Letters 114, 158102 (2015)

[6] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Statistical mechanics for natural flocks of birds, Proceedings of the National Academy of Sciences U.S.A. 109, 4786 (2012)

[7] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A Procaccini, M Viale and V Zdravkovic Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study Proceedings of the National Academy of Sciences U S A 105, 1232 (2008)

[8] A. Walther and A. H. E. M¨uller, Janus particles, Soft Matter 4, 663 (2008)

[9] J. Tailleur and M. E. Cates, Motility-induced phase separation, Physical Review Letters 100, 218103 (2008)

[10] M. E.

1. textbfIncrease predictive accuracy for MIPS phase boundaries in 3D systems through analytical modeling of pair interactions. The improved system can derive an analytical expression for the product of the interaction force and the pair-distribution function [18], which has been already used in Ref. [39] to find the spinodal and critical point of a system of ABPs. This allows for microscopic predictions beyond purely numerical simulation results, specifically by providing a predictive field theory [39].

2. textbfEnhance structural characterization by extracting orientation-resolved pair distributions for homogeneous states. The improved system can provide the fully orientation-resolved pair-distribution function [9], which is parametrized as g(r, θ1, θ2, ϕ2) and allows researchers to determine the regions where the symmetry assumptions introduced in Sec. II B are valid (namely the regions where the particle distribution is homogeneous) by analyzing a highresolution state diagram [5].

3. textbf Improve simulation efficiency and convergence for complex parameter spaces using analytical approximations of interaction kernels. By developing an analytical representation for the product function p(Pe, Φ0; r, θ1, θ2, ϕ2) [18], the system can utilize a Fourier expansion truncated at second order to approximate g with high accuracy, reducing statistical errors and allowing for faster fitting of simulation data [71].

4. textbf Enable high-fidelity modeling of structural features through optimized fitting procedures. The improved system employs an exponentially modified Gaussian distribution (EMG function) [25] combined with a linear cutoff term to fit Fourier coefficients, ensuring that all coefficients can be fitted by a product of the EMG function, the linear cutoff term, and a polynomial of a degree less than four [26]. This rigorous fitting procedure leads to reduced relative errors in approximating the simulation data.