Numerical Treatment of the Boltzmann Equation for Self-Propelled Particle Systems

Kinetic theories constitute one of the most promising tools to decipher the characteristic spatiotemporal dynamics in systems of actively propelled particles. In this context, the Boltzmann equation plays a pivotal role, since it provides a natural translation between a particle-level description of the system’s dynamics and the corresponding hydrodynamic fields. Yet, the intricate mathematical structure of the Boltzmann equation substantially limits the progress toward a full understanding of this equation by solely analytical means. Here, we propose a general framework to numerically solve the Boltzmann equation for self-propelled particle systems in two spatial dimensions and with arbitrary boundary conditions. We discuss potential applications of this numerical framework to active matter systems and use the algorithm to give a detailed analysis to a model system of self-propelled particles with polar interactions. In accordance with previous studies, we find that spatially homogeneous isotropic and broken-symmetry states populate two distinct regions in parameter space, which are separated by a narrow region of spatially inhomogeneous, density-segregated moving patterns. We find clear evidence that these three regions in parameter space are connected by first-order phase transitions and that the transition between the spatially homogeneous isotropic and polar ordered phases bears striking similarities to liquid-gas phase transitions in equilibrium systems. Within the density-segregated parameter regime, we find a novel stable limit-cycle solution of the Boltzmann equation, which consists of parallel lanes of polar clusters moving in opposite directions, so as to render the overall symmetry of the system’s ordered state nematic, despite purely polar interactions on the level of single particles.

Popular Summary

In the second half of the 19th century, Ludwig Boltzmann took a major leap to reconcile two apparently disparate notions—the atomistic and the hydrodynamic paradigms—and capture the physics of gaseous systems. His remarkable equation of motion for the one-particle distribution, which incorporates a particle-level description of binary interactions and lends itself to the derivation of a hydrodynamic description of gaseous systems, thereby connects the macroscopic kinetic coefficients to the microphysics of particle interactions.

A new class of gaseous systems, so-called active gases, has piqued the interest of many scientists. Unlike gas molecules in thermal equilibrium, the constituent particles of active gases are self-propelled, and mutual interactions are typically highly dissipative. Such systems are far from thermal equilibrium and exhibit a wealth of patterns on macroscopic scales, which ultimately arise as a consequence of activity and dissipation at the microscopic level. The Boltzmann equation is a conceptually appealing tool for understanding the unusual collective dynamics of such active systems. Previous Boltzmann-equation-based approaches to assessing the macroscopic dynamics of specific active systems largely relied on expansion techniques that substantially limited the range of applicability of the ensuing hydrodynamic equations. We propose an alternative approach to solve the Boltzmann equation numerically in two dimensions; the computation scheme we develop allows us to study parameter spaces even well beyond the onset of collective phenomena and also provides a highly flexible means to investigate spatially confined systems with arbitrary boundary geometries. In particular, we investigate an archetypical model system of self-propelled particles subject to binary, polar particle interactions. We find striking similarities between the ordering transition in this active model system and the liquid-gas phase transitions in equilibrium fluids. Moreover, we observe a novel stable limit-cycle solution of the Boltzmann equation, which consists of parallel lanes of polar clusters, moving in opposite directions.

Our findings, which reproduce hysteresis effects and density-segregated patterns, shed light on macroscopic collective motion and the transition between ordered and disordered states. We anticipate that future studies will build on our results to investigate systems with different interaction symmetries.

Figure 1

Spatial and angular tessellations: The system’s spatial domain $\mathrm{\Omega }$ is discretized using an arbitrary tessellation $\{{\mathfrak{s}}^{\alpha }{\}}_{\alpha }$ obeying Eqs. (7). A set of angular channels according to the tessellation (9) is attached to each spatial cell ${\mathfrak{s}}^{\gamma }$. Spatial transport: Convection in the bulk of the system along the direction ${\theta }_n$ (indicated by the dashed arrow) amounts to translating each spatial cell ${\mathfrak{s}}^{\gamma }$ by the convection length $\tau$. The convective transport operator ${\mathcal{T}}_{n\gamma }^{\alpha k}$ [Eq. (16b)] then redistributes the current value of $f_n^{\gamma }$ among all neighboring cells according to their intersection with the convected cell.

Figure 2

(a) Stationary solution for the one-particle distribution function $f(\theta )$ according to Eq. (21) for a spatially homogeneous system with $L=1$, $\sigma ={\sigma }_0=0.5$, and $K=32$, 64, and 128 angular slices for various overall densities ${\rho }_0$ (indicated in the figure). (b) Representative snapshot of a stationary traveling-wave pattern for ${\rho }_0=0.25$ and $\sigma ={\sigma }_0=0.5$. The color code indicates local density; the length and direction of arrows indicate the magnitude and orientation of the local particle flow. (c) Phase diagram for different values of the overall density ${\rho }_0$ and noise amplitude $\sigma ={\sigma }_0$. The analytic result for the phase boundary [refer to Eq. (34)] is indicated by the dashed line. A DSP is observed all along the phase boundary between the IHP and the SHP. The system size $L=100$, ${\epsilon }_x=5$, and $K=32$ for (b) and (c).

Figure 3

(a) Average momentum against density $\rho$ for two noise levels $\sigma ={\sigma }_0$ and three system sizes $L$; specific values are indicated in the figure (${\epsilon }_x=5$). The results indicate that the transition is first order. (b) Hysteresis study. The average momentum against overall density ${\rho }_0$ for $\sigma ={\sigma }_0=0.35$ and two system sizes $L\in \{50,100\}$. Solid symbols (with the paths indicated by solid arrows) indicate the system’s stationary state when starting from random initial conditions. Open symbols (with the paths indicated by open arrows) indicate the system’s stationary states upon lowering the density quasistatically starting from a stationary wavelike pattern at ${\rho }_0>{\rho }_t$. The system exhibits hysteresis, corroborating the assertion that the underlying ordering transition is first order.

Figure 4

Density variation inside the DSP parameter regime ($\sigma ={\sigma }_0=0.5$). (a) Typical profiles of wavelike patterns inside the density-segregated phase for various overall densities ${\rho }_0$. The asymmetry of the wave profiles increases with density. (${\rho }_t\approx 0.21$ for the noise values considered.) (b) Density segregation as a function of the overall density ${\rho }_0$. Data points have been obtained by quasistatically changing ${\rho }_0$ starting at a stationary traveling-wave pattern at ${\rho }_0=0.25$. As a consequence of hysteresis effects [see Figs. 3 and 5], phase-space coexistence can be observed between the DSP and the IHP at low densities, as well as between the DSP and the SHP at higher densities. (c) Surface fraction of the HD phase as a function of the overall density ${\rho }_0$. Within the DSP, the surface fraction varies virtually linearly with ${\rho }_0$ in agreement with the mean-field picture [Eq. (38)]. At the transitions $\mathrm{IHP}\to \mathrm{DSP}$ and $\mathrm{DSP}\to \mathrm{SHP}$, the surface fraction exhibits a discontinuous jump.

Figure 5

(a) Hysteresis study ($\sigma ={\sigma }_0=0.25$). Surface fraction of the HD phase against overall density ${\rho }_0$. Solid symbols (with their paths indicated by solid arrows) indicate the system’s stationary state upon increasing ${\rho }_0$ quasistatically, starting at a wavelike pattern at ${\rho }_0=0.25$. Open symbols (with their paths indicated by open arrows) indicate the system’s stationary states upon quenching the density to successively lower values, starting from a stationary, spatially homogeneous state computed at ${\rho }_0=0.415$. The system exhibits hysteresis, corroborating the assertion that the underlying ordering transition is first order. Inset: Dynamics of the density-segregation parameter $\eta$ (open data points) follows the subcritical bifurcation model [Eq. (39)] (solid line). (b) Sketch of the bifurcation diagram for the bifurcation parameter ${\rho }_0$ (at $\sigma ={\sigma }_0=0.5$). Solid lines indicate stable fixed points of the system’s dynamics and have been determined numerically. Dashed lines are drawn “by hand” and indicate the estimated position of unstable fixed points.

Figure 6

Cluster-lane patterns. (a) Particle currents in the cluster-lane solution of the active Boltzmann equation. Black curves represent instantaneous realizations of the longitudinally averaged, transverse particle currents $j_y(y,t)$ as functions of the (transverse) $y$ coordinate taken at different instances of time; see the snapshots in (b) and (c). The red curve represents the time-averaged transverse particle currents $\overline{j_y}(y)$. The gray shaded areas indicate the depletion zones where cluster grazing occurs, “feeding” the cluster zones (white areas). (b) Snapshot of a cluster-lane patten taken at the instant $t_1$ where the clusters of both lanes collide. (c) Snapshot of the same cluster-lane patten taken at some later instant $t_2>t_1$ where the clusters of both lanes spread due to transverse diffusion. (b),(c) The color code indicates local density levels, white arrows indicate the magnitude and direction of the local momentum, and red streamlines indicate macroscopic flow patterns.

Figure 7

Particle exchange across the cluster-lane boundaries in the long-time limit. In the plots, $t=0$ corresponds to a simulation time of $t_{\text{sim}}=10000{\lambda }^{-1}$, during which the cluster-lane pattern has developed from random initial conditions at $t_{\text{sim}}=0$. Particle exchanges across the cluster-lane boundaries are driven by cluster grazing (spikes) and lateral diffusion of clusters (bumps); see the inset. In the long-time limit, the magnitude of the instantaneous value of the net particle exchange $\dot{n}$ across the cluster-lane boundaries converges to 0.

Figure 8

Scaling of the coarsening dynamics inside the DSP regime at intermediate and long times. (a) Temporal evolution of ${\ell }_{\perp }(t)$ for a set of four different seeds (data points) of wave patterns. The length scale ${\ell }_{\perp }(t)$ roughly grows linearly with time and saturates at ${\xi }_{\phi }\equiv {\ell }_{\perp }(t\to \infty )=L$. (b) Temporal evolution of ${\ell }_{\perp }(t)$ for a set of six different seeds (data points) of cluster-lane patterns. The length scale ${\ell }_{\perp }(t)$ roughly grows linearly with time and saturates at ${\xi }_{\phi }\equiv {\ell }_{\perp }(t\to \infty )\approx 0.4L$. The solid red lines in (a) and (b) represent the averaged time evolution of ${\ell }_{\perp }$.

Figure 9

Illustration of the bending instability in the vicinity of the transition $\mathrm{DSP}\to \mathrm{IHP}$. Snapshots from left to right are taken at different times $T$. The left panel indicates the initial extended wave pattern, which becomes unstable toward a bending instability (middle panel); the right panel shows the reinforced stable wave pattern, which is contracted along the lateral dimension, as compared to the initial pattern.

Figure 10

Momentum and velocity for wavelike patterns within the DSP regime. (a) Momentum profiles corresponding to the density profiles displayed in Fig. 4. Velocity profiles morphologically resemble the corresponding density profiles and exhibit the same increase in asymmetry with increasing density. (b) Characteristic velocities inside the LD phase ($v_{\min }$) and the HD phase ($v_{\max }$), respectively. The middle curve gives the spatially averaged velocity $⟨v⟩$. All velocities are measured in units of $v_0$. The velocity characterizing the LD phase $v_{\min }$ attains finite values for ${\rho }_0\gtrsim {\rho }_v\approx 0.3$.

Figure 11

Cluster-lane geometry in periodic systems. (a) Schematic of the position of a cluster lane’s midline for winding number $w=2$ and length parameter $l=3$. (b) Illustration of a cluster-lane pattern with orientation $\alpha ={\alpha }_{1,1}=\pi /4$, observed in a system of linear extent $L=125$. Cluster lanes are indicated by dashed white lines, cluster lanes’ midlines indicated by thin solid lines. The direction of motion is indicated by arrows. (c) Asymmetric cluster-lane pattern observed in a system of linear extent $L=150$. The parameters for (b) and (c) are $\sigma ={\sigma }_0=0.5$ and $\rho =0.25$.

Figure 12

Illustration of the winding number $w$ and the length parameter $l$; see Eq. (b1). Each plot indicates the projection of the cluster lane’s midline onto the torus’s equatorial plane. (The system’s extent in the equatorial plane is indicated by a pair of dashed, concentric circles.) The winding number $w$ (rows) gives the number of full revolutions about the torus’s midpoint in the equatorial plane. The special case $w=0$ is indicated by the straight solid lines in the top row. The length parameter $l$ (columns) gives the total length of the closed contour representing the cluster lane’s midline.