Three-body correlations and conditional forces in suspensions of active hard disks

Self-propelled Brownian particles show rich out-of-equilibrium physics, for instance, the motility-induced phase separation (MIPS). While decades of studying the structure of liquids have established a deep understanding of passive systems, not much is known about correlations in active suspensions. In this work we derive an approximate analytic theory for three-body correlations and forces in systems of active Brownian disks starting from the many-body Smoluchowski equation. We use our theory to predict the conditional forces that act on a tagged particle and their dependence on the propulsion speed of self-propelled disks. We identify preferred directions of these forces in relation to the direction of propulsion and the positions of the surrounding particles. We further relate our theory to the effective swimming speed of the active disks, which is relevant for the physics of MIPS. To test and validate our theory, we additionally run particle-resolved computer simulations, for which we explicitly calculate the three-body forces. In this context, we discuss the modeling of active Brownian swimmers with nearly hard interaction potentials. We find very good agreement between our simulations and numerical solutions of our theory, especially for the nonequilibrium pair-distribution function. For our analytical results, we carefully discuss their range of validity in the context of the different levels of approximation we applied. This discussion allows us to study the individual contribution of particles to three-body forces and to the emerging structure. Thus, our work sheds light on the collective behavior, provides the basis for further studies of correlations in active suspensions, and makes a step towards an emerging liquid state theory.

Figure 1

Simulation snapshot of 4096 self-propelled disks at number density $\overline{\rho }=0.3$ and constant propulsion speed $v_0/d_{\mathrm{eff}}=5$. The system size is $L\times L$, with $L\approx 116.85$, and the directions of propulsion ${\widehat{e}}_i$ for each particle $i$ are shown by arrows. (b) Situation from within the snapshot in (a) with two tagged particles and our corresponding relative coordinates. The origin is fixed at the position of the first particle with the $x$ direction along its direction of propulsion. The second particle is located at the position $\vec{r}=(rcos\theta ,rsin\theta )$. The normalized basis vectors ${\widehat{e}}_r$ and ${\widehat{e}}_{\theta }$ are shown for the position of particle 2. (c) Situation from within the snapshot in (a) with two tagged particles interacting via two intermediate particles 3 and 4.

Figure 2

Sketch for two fixed hard disks labeled 1 and 2 and additional hard disks (dashed) in contact with the first one. The shaded area around particle 2 (red) is not accessible to a third particle due to the presence of the second particle. (a) Angles ${\theta }^{*}$ and ${\theta }_{\mathrm{\Delta }}$ and (b) example of the decomposition of the conditional force ${\vec{F}}_1$ acting on the first particle. The projected components are ${\left({\vec{F}}_1\right)}_r=({\widehat{e}}_r·{\vec{F}}_1){\widehat{e}}_r$ and ${\left({\vec{F}}_1\right)}_1=({\widehat{e}}_1·{\vec{F}}_1){\widehat{e}}_1$.

Figure 3

Pair distribution $g(r,\theta ,{\phi }_2)$ (first column) and conditional force ${\vec{F}}_1(r,\theta ,{\phi }_2)$ (second and third columns) from BD simulations at number density $\overline{\rho }=0.3$ and propulsion speed $v_0/d_{\mathrm{eff}}=5$. The first column shows the distribution of a second particle around the first one, the second column shows the projection of ${\vec{F}}_1$ onto the radial direction ${\widehat{e}}_{\mathrm{r}}$, and the third column shows the projection of ${\vec{F}}_1$ onto the orientation ${\widehat{e}}_1$. The relative position of the second particle with respect to the first one is given by the separations (a)–(c) $r\approx 1$, (d)–(f) $r\approx 1.5$, and (g)–(i) $r\approx 4$ and by the angle $\theta$, where $\theta =0$ corresponds to the position in front of the tagged first particle as sketched in Fig. 1. The relative orientation of the second particle with respect to the first one is given by ${\phi }_2$. To help interpret these plots both the relative position $\theta$ and the orientation ${\phi }_2$ are sketched in (g) for certain settings of particles 1 (black) and 2 (red) at the corresponding position in the plot.

Figure 4

Projected conditional force ${\widehat{e}}_1·{\vec{F}}_1$ on a tagged particle dependent on the relative position $(r,\theta )$ of a second particle as obtained from our BD simulations with a number density $\overline{\rho }=0.3$. Data are shown for (a) propulsion speed $v_0/d_{\mathrm{eff}}=5$ and (b)–(d) speeds 5, 15, 25, and 35. As indicated in (a), the data of (b)–(d) are shown along the cutting lines (b) along the positive $x$ axis (in the direction of propulsion), (c) along the positive $y$ axis (equal to the negative $y$ axis), and (d) along the negative $x$ axis.

Figure 5

(a) Pair distribution ${\langle g(r,\theta ,{\phi }_2)\rangle }_{{\phi }_2}$ and (b) and (c) projected conditional force ${\vec{F}}_1(r,\theta ,{\phi }_2)$ as shown in Fig. 3 ($v_0/d_{\mathrm{eff}}=5$ and $\overline{\rho }=0.3$), but averaged over the angle ${\phi }_2$ of the relative orientation of the second particle. The plots show the averaged simulation data from Fig. 3 (symbols), least-squares fits to the data in (b) and (c) as noted in the respective legend (dotted lines), and theoretical predictions (solid lines) from (a) Eq. (34), (b) Eq. (35), and (c) Eq. (37). For the theoretical predictions we use the parameters ${\zeta }_0=2\pi g_0\left(1\right)$ and ${\zeta }_1=\pi g_1\left(1\right)$, which we have calculated from our simulations via Eqs. (41) and (42).

Figure 6

Coefficients (a) and (b) $f_{\mathrm{a}}\pi /{\zeta }_1$, (c) and (d) $f_{\mathrm{b}}2\pi /{\zeta }_0$, and (e) and (f) $f_{\mathrm{c}}\pi /{\zeta }_1$ from theory and simulations as defined in Eq. (37). Theoretical curves are given by Eqs. (38, 39, 40) and simulation results follow from least-squares fits as shown in Fig. 5. Coefficients are shown at (a), (c), and (e) constant propulsion speed $v_0$ dependent on the number density $\overline{\rho }$ and (b), (d), and (f) constant number density dependent on the propulsion speed.

Figure 7

Coefficient $f_{\mathrm{b}}$ from theory and simulations as defined in Eq. (37) for a system of passive disks. Similar to Fig. 6, we show the coefficient $f_{\mathrm{b}}$ dependent on the number density $\overline{\rho }$. Note that data are shifted to enhance readability. The dashed lines mark zero for each number density $\overline{\rho }$.

Figure 8

Pair-distribution functions $g(r,\theta )$ around a self-propelled particle that is located at $(0,0)$ and swims in the positive $x$ direction. We show data at number density $\overline{\rho }=0.3$ and at propulsion speeds (a)–(c) $v_0=0$, (d)–(f) $v_0=5$, and (g)–(i) $v_0=20$. The full function is shown in (a), (d), and (g) and has the symmetry $g(r,\theta )=g(r,-\theta )$. Accordingly, we show data obtained from our BD simulations (BD) on the left half of the plot and numerical results from our theory (NUM) on the right half of the same plot. (b), (e), and (f) Values at particle contact, where we use the effective diameter as the particle-particle separation in our simulations. In addition, we show the theoretical predictions according to the parameters ${\zeta }_0$ and ${\zeta }_1$, which we have calculated from the respective data. Values are presented in Table 1. (c), (f), and (i) Function $g(x,0)$ along the positive $x$ axis in front of the tagged particle.