Kinetic Theory of Self-Propelled Particles with Nematic Alignment

Abstract

We present the results from kinetic theory for a system of self-propelled particles with alignment interactions of higher-order symmetry, particularly nematic ones. To this end, we use the Landau equation approach, a systematic approximation to the BBGKY hierarchy for small effective couplings. Our calculations are presented in a pedagogical way with the explicit goal of serving as a tutorial from a physicists’ perspective into applying kinetic theory ideas beyond mean-field to active matter systems with essentially no prerequisites and yield predictions without free parameters that are in quantitative agreement with direct agent-based simulations

1. Introduction

The seminal work on the Vicsek model [1,2] showing a phase transition to a state with long-ranged order in a two-dimensional model of inter-aligning self-propelled particles, a transition which would be in opposition to the Mermin–Wagner theorem [3] for equilibrium systems, can reasonably be used to mark the beginning of an ongoing foray from a theoretical perspective into understanding the “new” physics in models of active matter. In general, active matter is a subfield of nonequilibrium systems that are driven out of equilibrium on the smallest scale considered. The prototypical realization of this are self-propelled particles, as in the case of the Vicsek model. In a real-world version of such a system, there will be some microscopic mechanism (the usage of fuel, the metabolism of a living creature, external driving forces, etc.) that falls into the broader class of nonequilibrium physics, but if taken as fundamental units, these self-propelled particles are active.

In the decades since the first publication of the Vicsek model, considerable progress has been made in active matter physics, both on an analytical–computational level [4,5,6,7,8] and in the controlled experimental realization [9,10] of active matter systems. A somewhat comprehensive list of related review articles can be found in ref. [11]. The status of Vicsek-like models (an alternative nomenclature is that of dry active aligning matter [8]) as a workhorse in the research, particularly for more fundamental insights, has prevailed. This is not only limited to the original aim of studying flocking [12], the onset of collective motion or a global polar order in the directions of motion, but also extends to exploring systems with anti-aligning behavior [13,14,15,16,17,18,19], as well as systems with higher-order symmetry.

The most direct extension from polar order is the case of nematic order. Active nematics [20,21,22,23,24,25,26,27,28,29] are of high practical relevance, as many abundant instances of active-matter systems, such as bacteria and polymers, are based on particles with one distinct axis (rod-shaped). Typically, the alignment interaction in these systems is the result of actual collisions (in the sense that forces are acting on the spatial degrees of freedom) and the resulting symmetry is less clear, as inelastic collisions can generally induce a polar order [17,30,31]. Here, we consider the case of microscopically nematic pure alignment interactions. In line with the active matter concept, this can be seen as an abstraction to more microscopic interactions, such as steric interactions or by other means, say hydrodynamically [32].

On the technical side, interest in kinetic approaches to statistical physics has been renewed in the context of active matter [15,18,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50]. The extension of maximum entropy concepts (a discussion with some historical notes regarding the different conceptual schools can, for example, be found in ref. [51]) to these systems is not straightforward, whereas kinetic theory is systematically built up starting from the equations of motion and can therefore be adapted to a great variety of problems with few underlying assumptions. Specific Langevin dynamics being given, a systematic marginalization of the corresponding Fokker–Planck equation establishes the Born–Bogoliubov–Green–Kirkwood–Yvon (BBGKY) hierarchy of n-particle functions. Most physical quantities of interest can be written in terms of the one-particle function, making this a natural starting point. Closure of the hierarchy has been discussed in the literature [42,45] on polar models in detail. Often employed is the assumption of molecular chaos, which leads to mean-field dynamics (which can also be derived as the formal $N\to \infty$ limit of the equations of motion [35,52,53,54,55,56]) which have been referred to as the Vlasov (or McKean–Vlasov [52]) equation [57,58] (see also the earlier work by Jeans [59] and ref. [60] for an argument to call it the collisionless Boltzmann equation).The emergence of mean-field dynamics in certain continuum limits is the subject of active mathematical research, often being referenced as the propagation of chaos [43,52,61,62,63,64]. A full treatment of dynamics beyond mean-field is intricate and usually limited to perturbative approaches in either the density or the coupling strength. The density route, neglecting contributions from non-pair collisions, leads to the proper Boltzmann equation, which is valid for arbitrary interaction strengths [65]. As there are several technical intricacies related to this method (see the discussion of polar models in ref. [15]), we follow the route of weak interactions here. This leads [66] to the Landau equation, which arguably could also be named after Bogoliubov and/or Balescu, who provided more careful systematic derivations [67,68,69]. In particular, it becomes apparent under diligent derivation that the Landau equation is valid at arbitrary densities for sufficiently small coupling (the limit of higher densities, but weak interaction does only exist in models with soft interactions; relevant steric interactions, for example, are not suitable for a perturbative treatment in interaction strength). The Landau equation has also seen usage in the case of polar models [46,70]. There exist rigorous statements that the spatially homogeneous Landau equation that we will derive is indeed the proper limit of the Boltzmann equation with weak interactions [63,71,72].

The bottom-up nature of this kind of kinetic theory is complementary to more top-down mesoscopic hydrodynamic approaches that try to formulate the effective equations of motion for slow, collective variables from general principles [12,36,43,73,74,75].

In this work, we present some (in parts rather old) ideas on kinetic theory, particularly the Landau equation approach, in a concise and self-contained way that is tailored to the (rather new in the grand scheme of things) application in Vicsek-like models and, along the way, we derive novel results for Vicsek-like systems with nematic (and generally higher-order symmetry) coupling.

2. Model

We consider a class of generalized Vicsek-like models [34,35] wherein particles are described by two-dimensional positions ${\mathit{x}}_i$ ($i=1,\dots ,N$) and orientational angles ${\theta }_i$ whose dynamics are given by

$$\begin{array}{cc}\hfill \dot{{\mathit{x}}_i}& =v_0{\widehat{\mathit{n}}}_i\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill \dot{{\theta }_i}& =\mathsf{\Gamma }\sum _{j\in {\mathsf{\Omega }}_i}sin\left(n({\theta }_j\left(t\right)-{\theta }_i\left(t\right))\right)+{\eta }_i\hfill \end{array}$$

(1b)

with $\widehat{{\mathit{n}}_i}={\left(cos{\theta }_i,sin{\theta }_i\right)}^{\mathrm{T}}$, $\langle {\eta }_i\rangle =0$ and $\langle {\eta }_i\left(t\right){\eta }_j\left(t^{\prime }\right)\rangle =2D_r{\delta }_{i,j}\delta (t-t^{\prime })$. It is possible to formulate analogous dynamics in higher spatial dimensions [35,76]. Many aspects of our discussion are model-independent, but we limit the discussion to the two-dimensional case as the simplest non-trivial realization of orientational interactions whenever necessary. A schematic of our model is depicted in Figure 1. Crucially, this specific model does not include a prefactor that scales the interaction according to the number of interaction agents. This is essential in the limit of many interacting partners, as will be discussed below.

Thus, we consider ensembles of self-propelled particles that move with a fixed propulsion speed $v_0$ in the plane along the direction given by their orientational angle (relative to some fixed axis). These angles change via purely orientational or alignment interactions, whose strength is given by $\mathsf{\Gamma }$, and via explicit noise, whose strength is given by the resulting rotational diffusion constant $D_r$. Crucially, the interaction topology changes over time as the set of particles ${\mathsf{\Omega }}_i$ considered adjacent to particle i change. Specifically, throughout this work, we use

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_i=\{k\ne i||{\mathit{x}}_i-{\mathit{x}}_k|\le R\},\end{array}$$

(1c)

i.e., every particle has an identical fixed interaction radius R and interacts with every particle within a circle of radius R around it. For specific applications, other models of adjacency, such as metric-free models [39,77,78,79], can be more appropriate.

The interaction radius introduces a relevant length scale into the system, which gives rise to a suitable control parameter of the system, the non-dimensional density $M=\pi R^2\rho =\pi R^{2}N/L^2$. The intuitive meaning of this is that there are on average about M particles within the interaction radius of one particle. Similarly, we define the non-dimensional interaction strength $\mathrm{Sc}=\mathsf{\Gamma }R/v_0$ and the non-dimensional inverse diffusion strength or Peclet number $\mathrm{Pe}=v_0/\left(R{D}_r\right)$. This work focuses on the case of weakly interacting particles (to be specified later on) without any requirement regarding diluteness (i.e., arbitrary M).

Finally, the kind of interaction between any pair of angles is crucially affected by the chelation number $n\in \mathbb{N}$ in the argument of the interaction kernel. This allows us to control the symmetry of the alignment interaction, cf. Figure 2, with $n=1$ corresponding to the case of polar interactions and $n=2$ to nematic interactions. Also shown is the case of triatic interactions for $n=3$ as a representative for higher n-atic interaction (as observed, for example, in dihedral liquid crystals [80]) types with n-fold symmetry. A suitable order parameter is given by the bond order parameter $m_l={\sum }_i{e}^{il{\theta }_i}$ [81]. Generally, there are n stable and n unstable equilibria. While we are focused on $n=2$ and the case of $n=1$ has been treated [46,70] in the literature, we provide a general treatment for the sake of self-containedness. The symmetry of the interaction kernel directly affects the positional statistics observed in agent-based simulations (to be introduced in Section 4), see Figure 3.

3. Landau Kinetic Theory

By means of standard stochastic calculus [82,83], we find the Fokker–Planck equation for the full N-particle distribution function $P_N$ to the Langevin Equation (1) after rescaling to units such that $R=v_0=1$,

$$\begin{array}{cc}\hfill {\partial }_t{P}_N\left(\{{\mathit{x}}_i,{\theta }_i\}\right)& =-\sum _{i=1}^N\left(\left({\widehat{\mathit{n}}}_i·{\nabla }_i\right)P_N+\mathrm{Sc}{\partial }_{{\theta }_i}\left(\sum _{j=1}^N{a}_{i,j}sin\left(n({\theta }_j-{\theta }_i)\right)P_N\right)-{\mathrm{Pe}}^{-1}{\partial }_{{\theta }_i}^2{P}_N\right)\hfill \end{array}$$

(2)

wherein the omitted arguments of $P_N$ are the same on the right as on the left, which we will use throughout to increase legibility, and we introduce an indicator function

$$\begin{array}{cc}\hfill a_{i,j}& =\mathsf{\Theta }(1-|{\mathit{x}}_i-{\mathit{x}}_j|).\hfill \end{array}$$

(3)

The distribution $P_N$ contains the full information for all $3N$ microscopic degrees of freedom, which makes it unwieldy to use. As we consider all the particles to be identical, most relevant quantities, in particular the order parameter, should be reducible to integrals over a one-particle distribution. Therefore, we follow the textbook [69,84,85,86,87,88] Born–Bogoliubov–Green–Kirkwood–Yvon (BBGKY) hierarchy [67,89,90,91] approach and introduce reduced k-particle distribution functions by integrating out all the degrees of freedom except for those of the first k particles,

$$\begin{array}{cc}\hfill p_k& ={\displaystyle \frac{N!}{(N-k)!}}{\left({\displaystyle \frac{1}{M}}\right)}^k\int \dots \int P_{N}d{\mathit{x}}_{k+1}d{\theta }_{x+1}\dots d{\mathit{x}}_{N}d{\theta }_N.\hfill \end{array}$$

(4)

The combinatorial pre-factor accounts for possible other selections of the k pivot particles. The non-standard pre-factor ${(L^2/\left(\pi R^{2}N\right))}^k={(1/M)}^k$ is chosen such that all $p_k$ are of the order $\mathcal{O}\left(1\right)$ for a uniform distribution $P_N$. Using the shorthand notation $\mathit{i}={({\mathit{x}}_i,{\theta }_i)}^{\mathrm{T}}$, the evolution of the one-particle distribution, i.e., the first member of the BBGKY hierarchy, is given by

$$\begin{array}{c}\hfill {\partial }_t{p}_1(\mathbf{1},t)=\left({\mathcal{L}}_1+{\mathrm{Pe}}^{-1}{\partial }_{{\theta }_1}^2\right)p_1(\mathbf{1},t)-\gamma {\partial }_{{\theta }_1}\int \int a_{1,2}sin\left(n({\theta }_2-{\theta }_1)\right)p_2(\mathbf{1},\mathbf{2},t)d\mathbf{2}.\end{array}$$

(5a)

Here, we introduce the notation $\gamma =\mathrm{Sc}M$ for the effective interaction strength, the parameter that controls applicability of the theory, and make use of the one-particle Liouville operator ${\mathcal{L}}_i=-{\widehat{\mathit{n}}}_i{\nabla }_i$ determining the evolution in the free deterministic ($\mathrm{Sc}={\mathrm{Pe}}^{-1}=0$) problem. As for any system with additive pair interactions, the evolution of the k-particle function relies on the $k+1$-particle function establishing the BBGKY hierarchy. In particular, the two-particle function is found to evolve as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\partial }_t{p}_2(\mathbf{1},\mathbf{2},t)=\left({\mathcal{L}}_1+{\mathcal{L}}_2+{\mathrm{Pe}}^{-1}\left({\partial }_{{\theta }_1}^2+{\partial }_{{\theta }_2}^2\right)\right)p_2\hfill \\ \hfill & -\gamma [M^{-1}\left({\partial }_{{\theta }_1}\left(a_{1,2}sin\left(n\left({\theta }_2-{\theta }_1\right)\right)p_2\right)+{\partial }_{{\theta }_2}\left(a_{2,1}sin\left(n\left({\theta }_1-{\theta }_2\right)\right)\right)p_2\right)\hfill \\ \hfill & +{\partial }_{{\theta }_1}\int \int a_{1,3}sin\left(n\left({\theta }_3-{\theta }_1\right)\right)p_3(\mathbf{1},\mathbf{2},\mathbf{3},t)d\mathbf{3}+{\partial }_{{\theta }_2}\int \int a_{2,3}sin\left(n\left({\theta }_3-{\theta }_2\right)\right)p_3(\mathbf{1},\mathbf{2},\mathbf{3},t)d\mathbf{3}].\hfill \end{array}\end{array}$$

(5b)

The (abridged) hierarchy of Equation (5) is the starting point for any meaningful kinetic theory. Before moving on conceptually, it is very advantageous to separate off some contributions to the higher-order functions. In the free problem, all the particles are uncorrelated and one therefore has factorizing k-particle functions $p_{k,\mathrm{free}}={\prod }_{i=1^k}{p}_1(\mathit{i},t)$. In light of this, we perform an Ursell (or Mayer or cluster) expansion [92,93] introducing reduced correlation functions $g_i$. The correlation function $g_k$ contains the novel contributions on the k-particle level such that it contains the information about $f_k$ that is not contained in the lower-order $f_j$ and $g_j$ with $j<k$. In particular, the first three correlation functions are (omitting time dependencies for clarity here and in the following) given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_1\left(\mathbf{1}\right)& =g_1\left(\mathbf{1}\right)\hfill \\ \hfill p_2(\mathbf{1},\mathbf{2})& =p_1\left(\mathbf{1}\right)p_1\left(\mathbf{2}\right)+g_2(\mathbf{1},\mathbf{2})\hfill \\ \hfill p_3(\mathbf{1},\mathbf{2},\mathbf{3})& =p_1\left(\mathbf{1}\right)p_1\left(\mathbf{2}\right)p_1\left(\mathbf{3}\right)+p_1\left(\mathbf{1}\right)g_2(\mathbf{2},\mathbf{3})+p_1\left(\mathbf{2}\right)g_2(\mathbf{1},\mathbf{3})+p_1\left(\mathbf{3}\right)g_2(\mathbf{1},\mathbf{2})+g_3(\mathbf{1},\mathbf{2},\mathbf{3}).\hfill \end{array}\end{array}$$

(6)

We follow the diagrammatic approach of Balescu [69] to reduce the resulting equations for ${\partial }_t{g}_2$ to terms that actually contribute, utilizing a kinetic equivalent of the linked cluster theorem. We represent reduced one-particle distribution functions by unconnected horizontal lines that are tagged by their respective particle index; a k-tuple of horizontal lines connected by a vertical line represents a reduced k-particle correlation function. Each diagram contains one vertex that corresponds to the interaction kernel, $\sim \gamma {\partial }_{\theta }sin\left(n({\theta }^{\prime }-\theta )\right)$ (with a convention that the derivative is to be taken with respect to the (upper, where ambiguous) incoming line) and, finally, the lines corresponding to the given degrees of freedom start on the left and all other degrees of freedom (red lines in the figure) are to be integrated over. The important part to understand is that only fully connected (in the sense that every point on any line can be reached on lines from any other point) diagrams can contribute to the correlation function. In Figure 4, we show all the diagrams corresponding to the result of naively inserting Equation (6) into Equation (5b). The diagrams contained within the blue box are not fully connected and, therefore, do not contribute to $g_2$. Unaffected by these considerations is the free ($\mathrm{Sc}=0$) part of the time evolution. Reconstituting these considerations into equations, we find that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & ({\partial }_t-{\mathcal{F}}_{1,2})g_2(\mathbf{1},\mathbf{2})=-\gamma [M^{-1}\{{\partial }_{{\theta }_1}(a_{1,2}sin\left(n\left({\theta }_2-{\theta }_1\right)\right))(p_1\left(\mathbf{1}\right)p_1\left(\mathbf{2}\right)+g_2)\hfill \\ \hfill & +{\partial }_{{\theta }_2}(a_{2,1}sin\left(n\left({\theta }_1-{\theta }_2\right)\right))(p_1\left(\mathbf{1}\right)p_1\left(\mathbf{2}\right)+g_2)\}\hfill \\ \hfill & +{\partial }_{{\theta }_1}\int \int a_{1,3}sin\left(n\left({\theta }_3-{\theta }_1\right)\right)(p_1\left(\mathbf{1}\right)g_2(\mathbf{2},\mathbf{3})+p_1\left(\mathbf{3}\right)g_2+g_3(\mathbf{1},\mathbf{2},\mathbf{3}))d\mathbf{3}\hfill \\ \hfill & +{\partial }_{{\theta }_2}\int \int a_{2,3}sin\left(n\left({\theta }_3-{\theta }_2\right)\right)(p_1\left(\mathbf{2}\right)g_2(\mathbf{1},\mathbf{3})+p_1\left(\mathbf{3}\right)g_2+g_3(\mathbf{1},\mathbf{2},\mathbf{3}))d\mathbf{3}]\hfill \end{array}\end{array}$$

(7)

wherein ${\mathcal{F}}_{1,2}=\left({\mathcal{L}}_1+{\mathcal{L}}_2+{\mathrm{Pe}}^{-1}\left({\partial }_{{\theta }_1}^2+{\partial }_{{\theta }_2}^2\right)\right)$ is the free two-particle Fokker–Planck operator. As stated earlier, we are interested in the weakly interacting limit, which in light of the resulting equations, we can specify to be $\gamma =\mathrm{Sc}M\ll 1$. To this end, there is a fruitful insight to be drawn from the fact that only fully connected diagrams contribute to the $g_k$: any fully connected k-particle diagram will have to have one vertex and at least $k-2$ lines linked vertically by correlation functions, meaning that we can conclude that $g_k=\mathcal{O}\left({\mathrm{Sc}}^{k-1}\right)$. Thus, all of the hierarchical terms in any $g_k$ evolution, i.e., those originally belonging to $p_{k+1}$-terms and containing integrals, are of a higher order in $\gamma$. To order $\mathcal{O}\left({\mathrm{Sc}}^2\right)$, we therefore obtain the following Landau kinetic theory

$$\begin{array}{cc}\hfill {\partial }_t{p}_1\left(\mathbf{1}\right)& =\left({\mathcal{L}}_1+{\mathrm{Pe}}^{-1}{\partial }_{{\theta }_1}^2\right)p_1-\gamma {\partial }_{{\theta }_1}\int \int a_{1,2}sin\left(n\left({\theta }_2-{\theta }_1\right)\right)\left(p_1{p}_1\left(\mathbf{2}\right)+g_2(\mathbf{1},\mathbf{2})\right)d\mathbf{2}\hfill \\ \hfill & \equiv \left({\mathcal{L}}_1+{\mathrm{Pe}}^{-1}{\partial }_{{\theta }_1}^2\right)p_1-J_{\mathrm{coll}}\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill {\partial }_t{g}_2(\mathbf{1},\mathbf{2})& ={\mathcal{F}}_{1,2}{g}_2(\mathbf{1},\mathbf{2})+Q(\mathbf{1},\mathbf{2})\hfill \end{array}$$

(8b)

with the source term (using $a_{1,2}=a_{2,1}$)

$$\begin{array}{cc}\hfill Q(\mathbf{1},\mathbf{2})& =-\mathrm{Sc}{a}_{1,2}({\partial }_{{\theta }_1}-{\partial }_{{\theta }_2})sin\left(n\left({\theta }_2-{\theta }_1\right)\right)p_1\left(\mathbf{1}\right)p_2\left(\mathbf{2}\right).\hfill \end{array}$$

(8c)

It is a great strength of this approach that it is directly apparent that the collision integral $J_{\mathrm{coll}}=J_{\mathrm{coll}}\left(\mathbf{1}\right)$ consists of two contributions, a mean-field one, $J_{\mathrm{coll}}^{\mathrm{MF}}$, which is obtained by neglecting correlations and is of the order $\mathcal{O}\left(\gamma \right)$, and a genuinely Landau contribution, $J_{\mathrm{coll}}^{\mathrm{L}}$, of the order $\mathcal{O}\left(M{\mathrm{Sc}}^2\right)$ for which we have to know $g_2$ during interactions by solving Equation (8b). That is, we split

$$\begin{array}{cc}\hfill J_{\mathrm{coll}}& =J_{\mathrm{coll}}^{\mathrm{MF}}+J_{\mathrm{coll}}^{\mathrm{L}}\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill J_{\mathrm{coll}}^{\mathrm{MF}}& =\gamma {\partial }_{{\theta }_1}\int \int a_{1,2}sin\left(n\left({\theta }_2-{\theta }_1\right)\right)p_1{p}_1\left(\mathbf{2}\right)d\mathbf{2}\hfill \end{array}$$

(9b)

$$\begin{array}{cc}\hfill J_{\mathrm{coll}}^{\mathrm{L}}& =\gamma {\partial }_{{\theta }_1}\int \int a_{1,2}sin\left(n\left({\theta }_2-{\theta }_1\right)\right)g_2(\mathbf{1},\mathbf{2})d\mathbf{2}.\hfill \end{array}$$

(9c)

Enslaving $g_2$ to $p_1$, i.e., pretending we know $p_1\left(\mathbf{1}\right)$ and $p_1\left(\mathbf{2}\right)$ as functions of time, we can directly give this solution formally as a variation of constants

$$\begin{array}{c}\hfill g_2(\mathbf{1},\mathbf{2},t)={\mathcal{U}}_{1,2}(t_0,t,{\mathrm{Pe}}^{-1})g_2(\mathbf{1},\mathbf{2},t_0)+{\mathcal{U}}_{\mathbf{1},2}(t_0,t,{\mathrm{Pe}}^{-1}){\int }_{t_0}^t{\mathcal{U}}_{1,2}^{-1}(t_0,t^{\prime },{\mathrm{Pe}}^{-1})Q\left(t^{\prime }\right)d{t}^{\prime }\end{array}$$

(10)

with the free two-particle propagator

$$\begin{array}{c}\hfill {\mathcal{U}}_{1,2}(t_0,t,{\mathrm{Pe}}^{-1})=e^{{\mathcal{F}}_{1,2}(t-t_0)}.\end{array}$$

(11)

Operating in the weak noise limit, we set ${\mathcal{U}}_{1,2}(t_0,t,{\mathrm{Pe}}^{-1})\approx {\mathcal{U}}_{1,2}(t_0,t,0)$, neglecting noise contributions during the (typically short) timescale of interactions but keeping the diffusive terms arising from the (typically large) timescale of free flight. Additionally, the one-particle operators trivially commute and the system is symmetric under translation in time, so that we find that

$$\begin{array}{cc}\hfill {\mathcal{U}}_{1,2}(t_0,t,{\mathrm{Pe}}^{-1})& \approx {\mathcal{U}}_{1,2}(t_0,t,0)=e^{({\mathcal{L}}_1+{\mathcal{L}}_2)(t-t_0)}=U_1(t_0,t)U_2(t_0,t).\hfill \end{array}$$

(12)

Recalling the definition of the one-particle Liouville operator, ${\mathcal{L}}_i=-{\widehat{\mathit{n}}}_i{\nabla }_i$, which is simply the generator of infinitesimal translations, we see that the time propagation mediated by ${\mathcal{U}}_{1,2}$ is indeed a displacement by $\Delta \mathit{x}={\mathit{n}}_i\Delta t$. This is the physical expectation: in Landau kinetic theory, we compute the contributions from interactions assuming that particles move freely during the interaction, which is free ballistic transport in this problem. With this in mind, the back propagation mediated by ${\mathcal{U}}^{-1}(t_0,t)={\mathcal{U}}_{1,2}(t,t_0)$ is just a translation to the location where the agents had been at an earlier time.

To progress, two physically informed approximations are necessary. Firstly, we have to specify the meaning of the time $t_0$ and the value of the correlation function at this time. As stated previously, we only have good control over the propagator on short timescales and, by construction, the correlation function $g_2(\mathbf{1},\mathbf{2},t)$ only contains information if particles 1 and 2 have interacted. More importantly, as these particles subsequently interact with other (screened out) particles, any built-up correlation will decrease over time. This leads to the natural choice of assuming one-sided molecular chaos: we treat particles as uncorrelated at the beginning of a pair interaction [84]. The distinction of one-sidedness is important here, as the molecular chaos assumption (originally coined as the Boltzmann property by Kac [61]) usually refers to a mean-field description, i.e., $g_2\equiv 0$ throughout. Mathematically, these considerations provide $t_0$ with the meaning of the onset of an interaction (i.e., the point in time when the two particles enter each other’s interaction radii) and we can impose $g_2(\mathbf{1},\mathbf{2},t_0)=0$ as the particles being assumed to be uncorrelated at this point in time. We will address $t_0$ later on.

As an aside, we use physics terminology here, wherein mean-field is used to describe the effective dynamics derived from a neglection of correlations. It is an ongoing effort within the mathematics literature to show that the correlations are vanishing and that the long-term N-particle probability distribution functions factorize when considering the direct continuum ($N\to \infty$) limit of Equation (1). This, again following Kac [61], is referred to as the propagation of chaos [55,62,64,94]. For Vicsek-type models, this is discussed in ref. [55]. While we do omit contributions $\mathcal{O}(1/N)$ (from the combinatorial pre-factors), our perspective is that N is finite throughout (most ensembles of active matter particles such as flocks of birds or schools of fish are tiny, perhaps ${10}^{2–4}$, compared to the typical particle numbers considered in, say, gases, which are of the order of moles, so ${10}^{23}$), and if at all, the formal limit $N\to \infty$ should be performed at the end for actual observables such that the control parameters $M,\mathrm{Pe},\mathrm{Sc}$ remain finite.

Secondly, we know that in the relevant limits, ${\mathrm{Pe}}^{-1},\gamma \ll 1$, the one-particle distributions are spatially homogeneous on the scales of interactions reducing the action of the propagator solely to that which it has on the indicator function $a_{1,2}$, which evaluates to one during the interaction. In summation, we can write Equation (10) as

$$\begin{array}{cc}\hfill g_2(\mathbf{1},\mathbf{2},t)& =-\mathrm{Sc}{\int }_{t_0}^t{U}_{\mathbf{1},\mathbf{2}}(t^{\prime },t)a_{1,2}({\partial }_{{\theta }_1}-{\partial }_{{\theta }_2})sin\left(n\left({\theta }_2-{\theta }_1\right)\right)p_1(\mathbf{1},t)p_1(\mathbf{2},t)d{t}^{\prime }\hfill \\ & =-\mathrm{Sc}({\partial }_{{\theta }_1}-{\partial }_{{\theta }_2})sin\left(n\left({\theta }_2-{\theta }_1\right)\right)p_1(\mathbf{1},t)p_1(\mathbf{2},t){\int }_{t_0}^t{U}_{\mathbf{1},\mathbf{2}}(t^{\prime },t)a_{1,2}d{t}^{\prime }\hfill \\ & =-\mathrm{Sc}({\partial }_{{\theta }_1}-{\partial }_{{\theta }_2})sin\left(n\left({\theta }_2-{\theta }_1\right)\right)p_1(\mathbf{1},t)p_1(\mathbf{2},t)t_{\mathrm{passed}}\hfill \end{array}$$

(13)

Herein, $t_{\mathrm{passed}}=t_{\mathrm{passed}}(\mathbf{1}\left(t_0\right),\mathbf{2}\left(t_0\right))=t-t_0$ is the time passed during the encounter (up to time t), which depends on the initial coordinates. While the presentation here is focused on the specific case of pure alignment in two dimensions, the machinery is easily adapted to the case of higher dimensions or more complex interactions. However, the specific expressions as well as the evaluation of the passed time later on are highly dimensionally dependent.

As we consider purely orientational interactions, the dependency on the relative orientation $\Delta ={\theta }_2-{\theta }_1$ will be more critical than that on the relative positions $\Delta x={\mathit{x}}_2-{\mathit{x}}_1$. Delaying the specifics for a moment, it is helpful to keep track of the proceedings. The result of Equation (13), which was found by means of one-sided molecular chaos considerations to be the solution to Equation (8b), is the missing link to phrase Equation (8a) as a closed evolution equation for the one-particle function $p_1$. Specifically, it will be used to evaluate the Landau contribution to the collision integral.

3.1. Mode Equations

Thus, the relevant equation is the now closed first BBGKY-equation, which we write as

$$\begin{array}{cc}\hfill {\partial }_t{p}_1\left(\mathbf{1}\right)& =\left({\mathcal{L}}_1+{\mathrm{Pe}}^{-1}{\partial }_{{\theta }_1}^2\right)p_1\left(\mathbf{1}\right)-J_{\mathrm{coll}}^{\mathrm{MF}}\left[p_1\right]\left(\mathbf{1}\right)-J_{\mathrm{coll}}^{\mathrm{L}}\left[p_1\right]\left(\mathbf{1}\right).\hfill \end{array}$$

(14)

While we continue to postpone it, evaluation of the collision integrals is possible for a known $p_1$. Still, Equation (14) remains an integro-differential equation that is difficult to handle. The strategy that we pursue here is to make use of a Fourier transform. As we are particularly interested in the emergence or relaxation of angular order and generally working in a parameter regime close to the free problem, we only perform this transform with respect to the angular degrees of freedom and assume spatial homogeneity of $p_1$ for the rest of this work. This is to keep the presentation concise, and taking the full spatial dependencies into account poses no fundamental problem [15,48,95]. The convention of Fourier transformation we use is that the mode ${\widehat{p}}_m$ of $p_1$ to the wavenumber m is given by

$$\begin{array}{cc}\hfill {\widehat{p}}_m\left(t\right)& ={\displaystyle \frac{1}{2\pi }}{\int }_0^{2\pi }{e}^{-im\theta }{p}_1(\theta ,t)d\theta .\hfill \end{array}$$

(15a)

Consequently, the inverse transformation is given by

$$\begin{array}{cc}\hfill p_1(\theta ,t)& =\sum _{m=-\infty }^{\infty }{\widehat{p}}_m\left(t\right)e^{im\theta }.\hfill \end{array}$$

(15b)

As the one-particle function is real-valued, we have to have ${\widehat{p}}_m={\widehat{p}}_{-m}^{\ast }.$

We can generally establish an expectation for the mode equations, i.e., the time evolution of the Fourier modes, which is helpful for guiding one through the proceedings. The kinetic theory we derived is up to the order $\mathcal{O}\left({\gamma }^2\right)\sim \mathcal{O}\left(M^2\right)$, or more physically up to pair terms. Thus, the time evolution of a mode will be a quadratic form in all modes

$$\begin{array}{cc}\hfill {\partial }_t{\widehat{p}}_m& =-m^2{\mathrm{Pe}}^{-1}{\widehat{p}}_m+{\tilde{K}}_m^{\phantom{m}vw}{\widehat{p}}_v{\widehat{p}}_w,\hfill \end{array}$$

(16)

where we directly write the free part of the evolution due to rotational diffusion and use Einstein’s summation convention. We also know that our models are rotationally symmetric: any transformation ${\theta }_i\to {\theta }_i^{\prime }+\alpha$ amounts only to a different choice of reference axis. Therefore, both sides of the mode equations have to display identical behavior under this symmetry operation, making the gauge choice $\alpha$ immaterial. This leads to a sum rule, $w+v=m$, for the indices in the general approach, which also means that we can reduce the general tensor $\tilde{K}$ to a mode coupling matrix K. In summation, we find

$$\begin{array}{cc}\hfill {\partial }_t{\widehat{p}}_m& =-m^2{\mathrm{Pe}}^{-1}{\widehat{p}}_m+K_m^w{\widehat{p}}_w{\widehat{p}}_{m-w}.\hfill \end{array}$$

(17)

The matrix elements $K_m^w=V_m^w+L_m^w$ each will have contributions $V_m^w$ from the mean-field or Vlasov description and Landau corrections $L_m^w$ to it.

3.1.1. The Mean-Field Collision Integral

We recall that the mean-field collision integral was found to be, cp. Equation (9b),

$$\begin{array}{cc}\hfill J_{\mathrm{coll}}^{\mathrm{MF}}\left(\mathbf{1}\right)& =\gamma {\partial }_{{\theta }_1}\int a_{1,2}sin\left(n\left({\theta }_2-{\theta }_1\right)\right)p_1\left({\theta }_1\right)p_1\left({\theta }_2\right)d{\theta }_2\hfill \end{array}$$

(18)

where we used spatial homogeneity. Inserting Fourier transforms yields

$$\begin{array}{cc}\hfill V_m^w{\widehat{p}}_w{\widehat{p}}_{m-w}& =-\gamma {\int }_0^{2\pi }{\displaystyle \frac{d{\theta }_1}{2\pi }}{e}^{-im{\theta }_1}{\partial }_{{\theta }_1}{\int }_0^{2\pi }d{\theta }_{2}sin\left(n\left({\theta }_2-{\theta }_1\right)\right)\sum _{k,l=-\infty }^{\infty }{\widehat{p}}_k{\widehat{p}}_l{e}^{ik{\theta }_1}{e}^{il{\theta }_2}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & =-2\pi mi\gamma \sum _{l=-\infty }^{\infty }{\int }_0^{2\pi }d\Delta sin(n\Delta ){\widehat{p}}_{m-l}{\widehat{p}}_l{e}^{il\Delta }\hfill \end{array}$$

(20)

by integration by parts and changing integration variables. Inspecting the interaction kernel $sin(n\Delta )={\displaystyle \frac{1}{2i}}(e^{in\Delta }-e^{-in\Delta })$, performing the integral is trivial, and we can directly infer that

$$\begin{array}{cc}\hfill V_m^w& =m\pi \gamma ({\delta }_n^w-{\delta }_{-n}^w)\hfill \end{array}$$

(21)

Two conclusions are immediate from this result. One, the modes selected by the global chelation number n couple to the conserved mode ${\widehat{p}}_0$. Thus, they will be dominant hydrodynamically, which was to be expected. Secondly, there is an implied timescale ${\tau }_n^{-1}=n\pi \gamma$ for the relaxation dynamics of this dominant mode. This is also to be expected, as it reflects the perturbative behavior $sin\left(n({\Delta }_0+\epsilon )\right)\approx n\epsilon$ around the equilibrium positions.

The mean-field description should intuitively be the correct description for $M\gg 1$, as there will be no relevant correlations in the fully-connected model. This propagation of chaos property is indeed correct [35,55]. However, this limit should be approached carefully in this model. Inspection of the interaction in Equation (1) shows that $M\to \infty$ (i.e., $R\to \infty$ and effectively crossing over into a moving, fully-connected $XY$-model) is only sensible if we rescale $\mathsf{\Gamma }\to \mathsf{\Gamma }/M$, i.e., $\mathrm{Sc}\to \mathrm{Sc}/M$. We find corrections to the mean-field of order $M{\mathrm{Sc}}^2$, which are therefore negligible in this limit, validating earlier propagation of chaos results within our Landau approach.

3.1.2. The Landau Collision Integral

Collecting results, we can write the Landau collision integral as

$$\begin{array}{cc}\hfill J_{\mathrm{coll}}^{\mathrm{L}}(\mathbf{1},t)& =M{\mathrm{Sc}}^2{\partial }_{{\theta }_1}\int \int sin(n\Delta )t_{\mathrm{passed}}({\partial }_{{\theta }_1}-{\partial }_{{\theta }_2})sin(n\Delta )p_1{p}_1(\mathbf{2},t)d{\mathit{x}}_{2}d{\theta }_2.\hfill \end{array}$$

(22)

In this case, we cannot directly proceed by evaluating the spatial integration, as the passed time $t_{\mathrm{passed}}$ depends on it. We recall that this is meant to be the time an encounter lasted up to time t that started at $t_0$ with particles coordinates at time $t_0$ given by $\mathbf{1}$ and $\mathbf{2}$, respectively. Encounters have to end at a relative distance given by the interaction radius, $|{\mathit{x}}_2-{\mathit{x}}_1|=R$. Choosing a frame of reference such that ${\mathit{x}}_2\left(t_0\right)={\mathit{x}}_1\left(t_0\right)+{(-cos\varphi ,sin\varphi )}^T$, with $\varphi$ being the angle relative to the relative velocity ${\mathit{v}}_{\mathrm{rel}}={\mathit{v}}_2-{\mathit{v}}_1=({\widehat{\mathit{n}}}_2-{\widehat{\mathit{n}}}_1)$ whose modulus computes to $v_{\mathrm{rel}}=2|sin{\displaystyle \frac{\Delta }{2}}|$, we find from geometry that the distance traversed at time t is

$$\begin{array}{cc}\hfill s(t,t_0)& =x\left(t\right)+cos\varphi \left(t_0\right).\hfill \end{array}$$

(23)

wherein $cos\varphi \ge x\left(t\right)\ge -cos\varphi \ge 0$. We have to consider any possibility for $t_0$, which simply means that ${\mathit{x}}_2$ is uniformly distributed within the interaction radius around ${\mathit{x}}_1$ at time t. Thus, we can compute the spatially averaged passed time

$$\begin{array}{cc}\hfill t_{\odot }& ={\displaystyle \frac{2}{\pi }}{\int }_0^{1}dy{\int }_0^{\sqrt{1-y^2}}dx{\displaystyle \frac{x+\sqrt{1-y^2}}{v_{\mathrm{rel}}}}={\displaystyle \frac{4}{\pi v_{\mathrm{rel}}}}{\int }_0^{1}dy(1-y^2)={\displaystyle \frac{8}{3\pi v_{\mathrm{rel}}}}\hfill \end{array}$$

(24)

Plugging the averaged passed time into the Landau collision integral yields

$$\begin{array}{cc}\hfill J_{\mathrm{coll}}^{\mathrm{L}}({\theta }_1,t)& =M{\mathrm{Sc}}^2{\partial }_{{\theta }_1}\int sin(n\Delta )\pi t_{\odot }({\partial }_{{\theta }_1}-{\partial }_{{\theta }_2})sin(n\Delta )p_1\left({\theta }_1\right)p_1\left({\theta }_2\right)d{\theta }_2\hfill \\ \hfill & ={\displaystyle \frac{4}{3}}M{\mathrm{Sc}}^2{\partial }_{{\theta }_1}\int {\displaystyle \frac{sin(n\Delta )}{|sin\frac{\Delta }{2}|}}({\partial }_{{\theta }_1}-{\partial }_{{\theta }_2})sin(n\Delta )p_1\left({\theta }_1\right)p_1\left({\theta }_2\right)d{\theta }_2.\hfill \end{array}$$

(25)

Thus, we receive after Fourier transformation

$$\begin{array}{cc}\hfill & L_m^w{\widehat{p}}_w{\widehat{p}}_{m-w}={\displaystyle \frac{4miM{\mathrm{Sc}}^2}{3}}\sum _{v=-\infty }^{\infty }{\widehat{p}}_{m-v}{\widehat{p}}_v{\int }_0^{2\pi }d\Delta {\displaystyle \frac{sin(n\Delta )}{|sin\frac{\Delta }{2}|}}(-2ncos(n\Delta )+i(m-2v)sin(n\Delta ))e^{il\Delta }\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{4mM{\mathrm{Sc}}^2}{3}}\sum _{v=-\infty }^{\infty }{\widehat{p}}_{m-v}{\widehat{p}}_v{\int }_0^{2\pi }d\Delta {\displaystyle \frac{sin(n\Delta )}{|sin\frac{\Delta }{2}|}}(2ncos(n\Delta )sin(v\Delta )-(m-2v)sin(n\Delta )cos(v\Delta ))\hfill \end{array}$$

(27)

We outline the calculation of this ultimately trivial integral for the case of general n and, specifically, perform the calculation for $n=1,2$ in the Appendix A. After collecting terms, we end up at

$$\begin{array}{cc}\hfill L_m^w{|}_{n=1}=& mM{\mathrm{Sc}}^2{\displaystyle \frac{64}{3}}{\displaystyle \frac{\left(\left(4w^2-3\right)m-4w\right)}{(2w-1)(2w+1)(2w-3)(2w+3)}}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill L_m^w{|}_{n=2}=& mM{\mathrm{Sc}}^2{\displaystyle \frac{256}{3}}{\displaystyle \frac{\left(\left(64w^6-848w^4+2524w^2-1155\right)m-320w^5+2848w^3-4148w\right)}{(2w-1)(2w+1)(2w-3)(2w+3)(2w-5)(2w+5)(2w-7)(2w+7)}}.\hfill \end{array}$$

(29)

We show a graphical representation of these Landau contributions to the mode coupling matrix in Figure 5. The manifestation of the symmetry imposed onto the system via the chelation number n is less intuitive than for the mean-field contributions that connect the chelation mode to the density. In the second order, the structure of the corrections, $K_m^w{\widehat{p}}_{m-w}{\widehat{p}}_w$, means that there are two elements that contribute to the coupling of the mth mode to the density, $K_m^0$ and $K_m^m$. The second order part of the former is strongly negative, while that of the latter is weakly positive. Their sum, however, is negative, meaning that taking into account second-order contributions has (for $n=2$ and $n=1$) an effect that counteracts spontaneous alignment.

The coupling matrix elements $K_m^w$, or rather the full mode equations of Equation (17), contain the full dynamical information within the made assumptions (spatial homogeneity, cutoff at order ${\gamma }^2$). As such, they can be used to infer predictions on many relevant quantities. For example, we can consider the hydrodynamic expansion around an unordered state that we alluded to earlier, i.e., we assume that all $p_m$ for $m>0$ are small quantities. The zeroth mode $p_0={\displaystyle \frac{1}{2}}$ corresponds to the density and is a conserved quantity. This yields

$$\begin{array}{c}\hfill {\partial }_t{\widehat{p}}_m\approx -m^2{\mathrm{Pe}}^{-1}{\widehat{p}}_m+{\displaystyle \frac{1}{2}}(K_m^m+K_0^m){\widehat{p}}_m.\end{array}$$

(30)

For the mode that is selected by the symmetry encoded into the chelation number n, this means

$$\begin{array}{c}\hfill {\partial }_t{\widehat{p}}_n\approx -n^2{\mathrm{Pe}}^{-1}{\widehat{p}}_n+\gamma {\alpha }_n-{\gamma }^2{\beta }_n\end{array}$$

(31)

Here, we separate off the $\gamma$-dependencies from the Vlasov ($\alpha >0$) and Landau ($\beta >0$) contributions. Thus, it is directly evident that the onset of instability of this unordered state shifts towards lower values of ${\mathrm{Pe}}^{-1}$, i.e., lower noise.

4. Comparison with Agent-Based Simulation

The big advantage of kinetic theory based on the microscopic equations of motion versus more global hydrodynamic approaches is that we can validate our findings directly by comparison to agent-based simulation without any free parameters. Using an Euler scheme to numerically integrate Equation (1) [96,97,98], we measure the Fourier modes. This is done for the deterministic system (${\mathrm{Pe}}^{-1}=0$), starting from a polarized, spatially uniform state—see Figure 3 for some exemplary snapshots. Relaxation from an ordered state for $\gamma <0$ is a direct way of assessing the correctness of the angular mode couplings, as there will be no spatial patterns emerging (as would be the case for $\gamma >0$ or with strong noise). In particular, we initialize the directions uniformly $\theta \in [-\alpha /2,\alpha /2]$ with $\alpha =5\pi /6$. A secondary benefit of this choice is that the system is polarized along the x-axis, which in our convention means that ensemble-averaged Fourier modes will always be real-valued functions.

In Figure 6, we show that there is quantitative agreement between the measured mode dynamics in agent-based simulation and the results obtained from numerical integration of the mode equations, cp. Equation (17). The mode equations are numerically integrated by setting $p_m=0$ for all $\left|m\right|>10$. Here, we use rather small system sizes of $N=493$, but we do not expect qualitative differences with respect to the mode dynamics. For larger systems, spatial instabilities can become relevant, making the assumption of spatial homogeneity false [18,99,100]. However, omitting the spatial derivatives was done to keep the presentation concise, and fully accounting for them is not a principle difficulty. Thus, there are relevant contributions from terms beyond the mean-field at any N [42,50,101]. As discussed before, this is the result of fixing M and $\mathrm{Sc}$ to finite values as control parameters for any N. If the limit $N,M\to \infty$, $\mathrm{Sc}\to 0$ with $M\mathrm{Sc}$ fixed were to be considered instead, then we would expect the corrections to the mean-field to vanish for large N, establishing the propagation of chaos [35,55].

An important transport quantity that is suited to determination within this approach is the self-diffusion constant $D_{\mathrm{self}}$ at negative coupling when the system is self-mixing. Labeling a pivot particle with index 1, we define it as

$$\begin{array}{cc}\hfill D_{\mathrm{self}}& =\underset{t,\tau \to \infty }{lim}{\displaystyle \frac{1}{\tau }}\langle {({\mathit{x}}_1(t+\tau )-{\mathit{x}}_1\left(t\right))}^2\rangle \hfill \end{array}$$

(32)

By means of a Green–Kubo relation, this is directly linked to the velocity correlation time ${\tau }_c$, which determines the decay of the velocity autocorrelation function $\langle \mathit{v}(t+\tau )·\mathit{v}\left(t\right)\rangle \sim e^{-\tau /{\tau }_c}$ as ${\tau }_c=2D_{\mathrm{self}}$. As has been pointed out before in ref. [15], this quantity does not exist in mean-field calculations. The calculation via a Boltzmann–Lorentz-type theory is easily adapted, and we do not reproduce it in full here, but a rough argument would be as follows: if we label a very small fraction as tracer particles with modes ${\widehat{q}}_m$, then their interaction partners will essentially always be with regular particles whose modes we denote by ${\widehat{p}}_m$. Going through the derivation of Equation (17), it is apparent that the analogue for modes of the tracer particles, which explicitly takes into account spatial derivatives as is necessary for diffusive phenomena, would be

$$\begin{array}{cc}\hfill {\partial }_t{\widehat{q}}_m+{\displaystyle \frac{1}{2}}({\nabla }^{-}{\widehat{q}}_{m-1}+{\nabla }^{+}{\widehat{q}}_{m+1})& =-m^2{\mathrm{Pe}}^{-1}{\widehat{q}}_m+K_m^w{\widehat{q}}_{m-w}{\widehat{p}}_w\hfill \end{array}$$

(33)

with the shorthand ${\nabla }^{\pm }={\partial }_x\pm i{\partial }_y$. Again taking the limit of an unordered background, ${\widehat{p}}_m\propto {\delta }_{m,0}$, we see that the mean velocity of the tracer particles, ${\widehat{q}}_1$, only couples with the background density via $K_1^0$, an element that does not exist in the mean-field irrespective of n. In particular, we find a decay of the velocity autocorrelation function consistent with

$$\begin{array}{cc}\hfill D_{\mathrm{self}}& =-{\displaystyle \frac{1}{K_1^0-{\mathrm{Pe}}^{-1}}}.\hfill \end{array}$$

(34)

We present numerical results for $D_{\mathrm{self}}$ as inferred from the velocity auto-correlation function in Figure 7. The data are found in quantitative agreement for sufficiently small values of $M{\mathrm{Sc}}^2$, which is the order first correction beyond mean-field theory. Similar to the discussion of the polar case [15,48], the scaling itself (and thus the validity of only considering Landau equation contributions) and, at least qualitatively, corrections to scaling can be understood by a random telegraph model (labeled semi-analytical in Figure 7), which we reproduce in the Appendix B. Interestingly and not explicitly captured in the telegraph process description, the diffusion coefficients’ dependency on the chelation number n is non-trivial. At ${\mathrm{Pe}}^{-1}=0$, we find $D_{\mathrm{self}}^{n=2}/D_{\mathrm{self}}^{n=1}={\displaystyle \frac{35}{44}}\approx 0.8$. This is indicative of self-diffusion being a result of the full-interaction kernel, whereas other properties could be understood from the perspective of linearization around equilibria leading to power-law scaling in n.

5. Discussion

In summation, we gave a concise introduction into Landau equation-based methods of kinetic theory, as they can be applied to active matter with few prerequisites. This was performed for a specific Vicsek-like model with interactions of n-fold symmetry. By means of careful inspection of the relevant scalings, it becomes apparent that the Landau equation is not the mere weak coupling limit of the Boltzmann equation, but rather a separate approach to finding a solution to the BBGKY hierarchy. The Boltzmann equation is fundamentally an expansion in the density, whereas the Landau equation is an expansion in the coupling strength, valid at any density (as long as the effective interactions $M\mathrm{Sc}$ are sufficiently small).

We derived the full mode equations, establishing the connection between the mean-field (or Vlasov or McKean) and the Landau contributions to the assumptions of molecular chaos and the more physical assumption of one-sided molecular chaos, respectively. The former neglects all inter-particle correlations, whereas the latter takes correlations during interactions into account. The main technical results are the elements beyond the mean-field of the mode coupling matrix (see Figure 5).

As a result of this, we are able to analytically determine an asymptotic expression for the self-diffusion coefficient of the self-propelled particles, a measure for the self-mixing they display at anti-aligning couplings even in a deterministic system. This result is then independently corroborated by a random-telegraph process that gives identical asymptotics and offers some insight into corrections to scaling at larger couplings. As such, we give an example demonstrating that first-principle kinetic theories for active matter, when properly derived from the exact Liouville– or Fokker–Planck equation using the weaker, generic assumption of one-sided molecular chaos, lead to quantitative agreement, even in transient/dynamic situations far from stationary states. This is easily underappreciated in favor of more qualitative descriptions based on strong, problem-specific assumptions [8].

We hope that this work serves as one element refamiliarizing a broader community with the ideas and merits of approaching statistical physics kinetically starting from microscopic equations of motion. Going forward, these equations will provide a valuable set of tools in the understanding of active matter.