Transport coefficients of self-propelled particles: Reverse perturbations and transverse current correlations

TL;DR

This paper extends the reverse perturbation method to active particles modeled by the Vicsek model, deriving transport coefficients and comparing theoretical predictions with simulations, revealing deviations from mean-field assumptions.

Contribution

It introduces a novel approach to measure transport coefficients in active particle systems and derives a new collisional contribution to viscosity using kinetic theory.

Findings

Transport coefficients agree well near the transition to collective motion.

Measured viscosity and amplification coefficients are slightly higher than mean-field predictions.

The mean-field assumption of molecular chaos is less valid in velocity-alignment systems.

Abstract

The reverse perturbation method [Phys. Rev. E 59, 4894 (1999)] for shearing simple liquids and measuring their viscosity is extended to the Vicsek-model (VM) of active particles [Phys. Rev. Lett. 75, 1226 (1995)] and its metric-free version. The sheared systems exhibit a phenomenon that is similar to the skin effect of an alternating electric current: momentum that is fed into the boundaries of a layer decays mostly exponentially towards the center of the layer. It is shown how two transport coefficients, i.e. the shear viscosity $\nu$ and the momentum amplification coefficient $\lambda$, can be obtained by fitting this decay with an analytical solution of the hydrodynamic equations for the VM. The viscosity of the VM consists of two parts, a kinetic and a collisional contribution. While analytical predictions already exist for the former, a novel expression for the collisional part is derived by an Enskog-like kinetic theory. To verify the predictions for the transport coefficients, Green-Kubo relations were evaluated and transverse current correlations were measured in independent simulations. Not too far to the transition to collective motion, we find excellent agreement between the different measurements of the transport coefficients. However, the measured values of $\nu$ and $1-\lambda$ are always slightly higher than the mean-field predictions, even at large mean free paths and at state points quite far from the threshold to collective motion, that is, far in the disordered phase. These findings seem to indicate that the mean-field assumption of molecular chaos is much less reliable in systems with velocity-alignment rules such as the VM, compared to models obeying detailed balance such as Multi-Particle Collision Dynamics.

Equations

$${\mathbf{x}}_i(t+\tau )={\mathbf{x}}_i(t)+\tau {\mathbf{v}}_i(t).$$

$${\Phi }_i=\mathrm{Arg}\left(\sum _{\{j\}}{\mathrm{e}}^{i{\theta }_j}\right)$$

$${\theta }_i(t+\tau )={\Phi }_i+{\xi }_i,$$

$$v_{\mathrm{a}}=\frac{1}{N{v}_0}\left|\sum _{i=1}^N{\mathbf{v}}_i\right|.$$

$${\eta }_{\mathrm{c}}=\sqrt{48M\left(\frac{2}{\pi }-\frac{1}{2}\right)}.$$

$$⟨\sigma ⟩=\frac{⟨\Delta p_x⟩}{2\Delta t{L}_x},$$

$$\nu ={\nu }_{\mathrm{kin}}+{\nu }_{\mathrm{coll}}.$$

$${\nu }_{\mathrm{kin}}=\frac{v_0^2\tau }{8}\frac{1+p}{1-p}.$$

$$p=\frac{4}{\eta }\sin \eta \sum _{n=1}^N\frac{{\mathrm{e}}^{-M}}{n!}{n}^2{M}^{n-1}{K}_{2c}^{11}(n),$$

$$p\approx \frac{\sin \eta }{\eta }\frac{1+0.327M^2+0.072M^3}{1+M+M^2/2+M^3/6}.$$

$$p\approx \frac{\sin \eta }{2\eta }.$$

$${\nu }_{\mathrm{coll}}=\frac{R^2}{\tau }\frac{\sin (\eta /2)}{2\eta }\sum _{n=1}^{\infty }\frac{{\mathrm{e}}^{-M}}{(n-1)!}{M}^n{K}_C^1(n+1).$$

$${\nu }_{\mathrm{coll}}\approx \frac{M{R}^2}{\tau }\frac{\sin (\eta /2)}{8\eta }\sqrt{\frac{\pi }{M+2}},\left[1+\frac{3M}{8(M+2{)}^2}\right]$$

$${\nu }_{\mathrm{coll}}\approx \frac{M{R}^2}{\tau }\frac{\sin (\eta /2)}{2\eta }\left[\frac{1/\pi +0.2624M+0.11245M^2+0.03347M^3}{1+M+M^2/2+M^3/6}\right].$$

$${\mu }_1{u}_y{\partial }_y{u}_x+{\mu }_2{u}_x{\partial }_y{u}_y$$

$$({\mu }_1+{\mu }_2+2{\mu }_3)u_y{\partial }_y{u}_x+2{\mu }_3{u}_x{\partial }_y{u}_x$$

$$\kappa$$

$$u_x=d_0\sinh (d_{1}y),$$

$$d_1=\sqrt{\frac{1-\lambda }{\nu \tau }}$$

$$\sigma ={\nu {\rho }_m\frac{\partial u_x}{\partial y}|}_{y=L_y/4}$$

$$\nu =\frac{\sigma }{{\rho }_m{d}_0{d}_1\cosh (d_1{L}_y/4)}$$

$$1-\lambda =\tau d_1^2\nu$$

$$\nu u_x^{{}^{\prime \prime }}-(\kappa +q{u}^2)u_x$$

$$\nu u_y^{{}^{\prime \prime }}-(\kappa +q{u}^2)u_y$$

$$0=\frac{\partial }{\partial y}\left[-\frac{\kappa }{2}{u}^2-\frac{q}{4}{u}^4+\frac{\nu }{2}{\left(\frac{\partial \mathbf{u}}{\partial y}\right)}^2\right].$$

$$\mathbf{u}$$

$$d_1$$

$$⟨u_x⟩=d_2+d_0\sinh (d_{1}y).$$

$$\lambda -1=\frac{\tau d_1^2\nu }{2}.$$

$${\partial }_t{w}_{\alpha }=-{\partial }_{\alpha }P-\kappa w_{\alpha }+{\partial }_{\beta }{\sigma }_{\alpha \beta }+{\Psi }_{\alpha }.$$

$${\sigma }_{\alpha \beta }=\nu \left({\partial }_{\alpha }{w}_{\beta }+{\partial }_{\beta }{w}_{\alpha }-\frac{2}{d}{\delta }_{\alpha \beta }{\partial }_{\gamma }{w}_{\gamma }\right),$$

$${\hat{w}}_{\alpha }(\mathbf{k},t)=\frac{1}{V}\int w_{\alpha }(\mathbf{x},t){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}}\mathrm{d}\mathbf{x},$$

Full text

PACS numbers:87.10.-e, 05.20.Dd, 64.60.Cn, 02.70.Ns

I Introduction

During the past two decades, there has been a large interest in active matter systems, such as bird flocks animal_flocks , swarming bacteria kearns_11 ; copeland_09 , active colloids zoettl_14 ; ginot_15 , microtubule mixtures nedelec_02 and actin networks actin_net driven by molecular motors. These systems display interesting behaviors such as pattern formation, collective motion and non-equilibrium phase transitions vicsek_12 ; marchetti_13 . Some of these features already occur in one of the simplest models for active matter, the Vicsek-model (VM) of self-propelled particles vicsek_95 ; czirok_97 ; nagy_07 and its variants peruani_08 ; aldana_09 ; peng_09 ; barbaro_09 ; ginelli_10a ; chou_12 ; mishra_12 ; romensky_14 . Because of the simplicity of its interaction rules and the existence of a non-standard transition to a collective state of polar order, the VM became an archetype of active matter.

Due to the many degrees of freedom, theoretical studies of active matter systems are often based on coarse-grained macroscopic transport equations for the slow variables such as density or momentum. Originally, the general forms of these equations were postulated by symmetry and renormalization group arguments, such as in the seminal Toner-Tu theory toner_95 ; toner_98 ; toner_12 for polar active matter. However, this approach leaves the coefficients of the terms in the transport equation largely undetermined. Furthermore, memory and other nonlocal terms are usually not considered, although for particular models there is evidence on their relevance kuersten_17 . These shortcomings motivated many researchers to derive macroscopic transport equations directly from the microscopic interactions and to obtain explicit expressions for the occurring coefficients bertin_06 ; baskaran_08a ; baskaran_08b ; bertin_09 ; ihle_11 ; peshkov_12 ; farrell_12 ; grossmann_13 ; thueroff_13 ; hanke_13 ; ihle_14_a ; peshkov_14_b ; chepizhko_14 . Even though most of these coarse-graining approaches are based on some type of mean-field assumption, they can still be rather involved and, in addition, rely on further approximations such as time scale separation, the thermodynamic limit or the irrelevance of higher order spatial gradients.

In principle, different kinetic theory approaches for the same microscopic model can lead to different macroscopic expressions, see for example Ref. 38. Thus, the validity of the derived transport coefficients and of the macroscopic description in general is often questionable and has led to debates aldana_09 ; ihle_13 ; peshkov_14_b_C ; ihle_14_b ; ihle_14_c . To the best of our knowledge, so far there has been no comprehensive work on the verification of transport coefficients in polar active matter. In this article, we start filling this void, at least for two transport coefficients of the standard and the metric-free (or topological) VM. In particular, we perform extensive agent-based simulations of the VM, measure the kinematic viscosity, $\nu$, and the momentum amplification coefficient, $\lambda$, by means of several complementary methods, and compare them to predictions from kinetic theory. To test our numerical tools, we perform additional measurements on a momentum-conserving, particle-based model for computational fluid dynamics, called Multi-Particle Collision Dynamics (MPCD) fluid malevanets_99 ; malevanets_00 ; gompper_08 ; kapral_08 ; howard_18 ; howard_19 .

First, we employ a non-equilibrium approach and measure the response to shear, generated through the reverse perturbation (RP) method mueller_plathe_99 . When applying RP to the VM, one observes a phenomenon that is similar to the skin effect of an alternating electric current: momentum that is fed into the boundaries of a channel decays mostly exponentially towards the center of the channel. We show how $\nu$ and $\lambda$, can be obtained by fitting this decay with an analytical solution of the hydrodynamic equations for the VM. In analogy to the MPCD case, the viscosity of the VM consists of two parts, the kinetic ($\nu_{\rm kin}$) and the collisional viscosity ($\nu_{\rm coll}$). The latter contribution was missing in previous theories of the VM, and we derive here an analytical expression for $\nu_{\rm coll}$ by extending the Enskog-like kinetic theory from Refs. 29; 50.

Furthermore, we introduce and apply two other methods to evaluate $\nu$ and $\lambda$ – the transverse current correlation method (TC) hoheisel_88 ; palmer_94 and the Green-Kubo (GK) method green_54 ; kubo_57 ; zwanzig_65 ; forster_75 , which, in contrast to RP, operate without introducing velocity gradients. We discuss colored noise and estimate memory effects in the theory for TC, which explains limitations of this method at very small mean free path. We find excellent agreement between the measurements of the RP and the TC method. Reasonable quantitative agreement between agent-based simulations (using the RP, GK and TC method) and predictions by kinetic theory is observed. This supports previous concerns on the validity of the mean-field assumption of molecular chaos in systems without detailed balance and underlines the need for a theory that includes correlation effects.

The rest of the manuscript is organized as follows. In Sec. II we present the standard and metric-free version of the VM, and briefly discuss the RP method for applying shear flow. In Sec. III we develop analytic expressions for the transport coefficients $\lambda$ and $\nu$, and discuss the theory of the TC method. In Sec. IV we shows our numerical results and compare them to the theoretical predictions. Finally, in Sec. V we provide a final discussion of our results and briefly outline open questions. The Appendix contains additional theoretical details related to the kinetic theory of the VM, as well as results for an MPCD fluid under shear.

II Model and Methods

II.1 The standard Vicsek model

The VM vicsek_95 ; czirok_97 ; nagy_07 consists of $N$ point particles at global number density $\rho_{0}$, which move at constant speed $v_{0}$ in two dimensions. The positions and velocities of the particles at time $t$ are given by ${\bf x}_{i}(t)$ and ${\bf v}_{i}(t)$, respectively. In the VM, the particles are propagated via sequential streaming and collision steps with time step $\tau$. (The term “collision” should not to be taken literally, but instead it just denotes any action that changes the momentum of a particle.) During the streaming step, the particles move ballistically

[TABLE]

Because the speeds of the particles stay the same at all times, the velocities are parameterized by the “flying” angles, $\theta_{i}$, i.e. ${\bf v}_{i}=v_{0}(\cos{\theta_{i}},\sin{\theta_{i}})$.

In the collision step, the directions $\theta_{i}$ are changed so that the particles align with their neighbors within a fixed distance $R$ plus some external noise. In practice, a circle of radius $R$ is drawn around the focal particle $i$, and the average direction $\Phi_{i}$ of motion of the particles within the circle is determined according to

[TABLE]

where the sum goes over all particles within the interaction range $R$ (including particle $i$). Once all average directions $\Phi_{i}$ are known, the new directions follow as

[TABLE]

where $\xi_{i}$ is the so-called angular noise. The random numbers $\xi_{i}$ are uniformly distributed in the interval $[-\eta/2,\eta/2]$, with noise strength $\eta$. The model uses parallel updating, and in this paper we will also assume the so-called standard VM which uses a forward-updating rule. Thus, the already updated positions ${\bf x}_{i}(t+\tau)$ are used for determining the average directions $\Phi_{i}$ at time $t$.

Another relevant model parameter is the average particle number $M$ that can be found inside a circle of radius $R$, i.e. $M=\rho_{0}\pi R^{2}$ with the global number density $\rho_{0}$. The dimensionless parameter $M$ measures the ratio of the interaction range $R$ to the average particle distance $1/\sqrt{\rho_{0}}$. By increasing $M$ and/or decreasing the noise $\eta$, the VM can be driven from a disordered phase to a phase of collective motion. The degree of alignment of the particle velocities can be quantified through the polar order parameter

[TABLE]

Assuming a spatially homogeneous system, the threshold condition for this non-equilibrium phase transition can be calculated within mean-field kinetic theory (see Appendix A for details). For sufficiently small $M\ll 1$, the threshold noise $\eta_{\rm c}$ is predicted as

[TABLE]

For parameters where the molecular chaos assumption is strongly violated, $\eta_{\rm c}$ can be much lower than this theoretical prediction, sometimes by a factor between two and three. For more details on the calculations and for a discussion of this transition, see Refs. 29; 39; 57; 58; 50.

II.2 The Vicsek model with topological interactions

Recent experiments by Ballerini et al. ballerini_08 ; cavagna_10 on flocks of starlings indicated that a Vicsek-like interaction rule with a fixed interaction range might not be appropriate for animal flocks. Instead, a statistical analysis revealed that, on average, each bird interacts with a fixed number of neighbors, typically six to seven. This constitutes a topological or metric-free interaction because the metric distance is not relevant; rather, it is a question of who the closest neighbors are. Ballerini et al. argued further that, due to evolutionary pressure, the main goal of interaction among individuals is to maintain cohesion. By comparing simulations with the regular VM and a modified VM with metric-free interactions, they found that flocks, when facing predators, kept cohesion much better in the metric-free model. These observations inspired several other groups to study versions of the VM with topological interactions.

In this paper we will focus on a simple modification of the VM, which was suggested by one of us chou_12 , because it allows an analytical description by a similar Enskog-like kinetic theory as the one outlined in Appendix. A. In this model, the alignment rule of the regular VM, given by Eqs. (2) and (3), is slightly modified such that the number of particles in every collision circle is kept constant and equal to $M$ at all times by locally adjusting the interaction radius. Thus, only the $M-1$ closest neighbors together with particle $i$ itself are included in the calculation of the average angle $\Phi_{i}$ of a particular particle $i$. This procedure leads to large interaction ranges in areas with sparse populations, whereas the interaction radius becomes small at locations with a high particle number density. We introduce an effective interaction radius for this metric-free model, $R_{\rm eff}=\sqrt{M/(\pi\rho_{0})}$, which is always set to one by appropriately choosing the particle number density $\rho_{0}$.

In the regular VM, a larger local particle density leads to more robust alignment and stronger local order. This behavior can be seen in the phase diagram of the VM, for example Fig. 1 in Ref. 29. This coupling between density and order is the main reason behind the occurrence of soliton-like density waves near the order/disorder threshold in the regular VM ihle_13 . In the topological VM, however, density and order are decoupled because it is always the same number of particles that participate in the alignment interaction. Therefore, the long-wave length instability of the regular VM as well as the density waves are absent in the topological VM chou_12 ; peshkov_12_b .

II.3 The reverse perturbation method

We performed non-equilibrium simulations to compute the shear viscosity from the simulations. These approaches often provide significantly better signal to noise ratios compared to equilibrium methods, such as the GK relation green_54 ; kubo_57 ; zwanzig_65 ; forster_75 . To generate shear flow in our system, we employed the RP method mueller_plathe_99 , where the shear stress on the system is imposed externally, by generating a momentum flux through a slab perpendicular to the flow direction. This flux is achieved by swapping the particle velocities in the following way: first, the periodic simulation box is subdivided into equally sized slabs with thickness $a$ along the gradient direction of the flow ($y$). Then, particle $i$ in the $y=0$ slab with the largest positive $e_{x}$ value and particle $j$ in the $y=L_{y}/2$ slab with the largest negative $e_{x}$ value are identified, and their velocities are swapped. This swapping procedure artificially generates a momentum flux, which gives rise to a physical flow.

If both particles have the same mass, as is the case in all our models, swapping conserves both the linear momentum and the global kinetic energy. In our implementation, momentum swaps were applied to the system with equal probability either before or after the collision step. Note, that when using the RP method for the VM, it is crucial to also swap $e_{y}$ of the particle pair so that the particle speed $v_{0}$ is conserved. We verified that this additional swapping does not introduce an unwanted momentum flux in the $x$ direction.

The imposed shear stress can be controlled by the amount of mome*•*ntum swaps in one step and by the time between swaps, $\Delta t$. For the chosen geometry of our two-dimensional systems, the average shear stress can be computed as:

[TABLE]

where $\langle\Delta p_{x}\rangle$ is the $x$ component of the average total momentum exchanged during one time step. Figure 1 shows a schematic view of the shear procedure and the emerging flow profile.

III Theory

III.1 The collisional viscosity $\nu_{\rm coll}$

It has been shown by several groups tuzel_03 ; kikuchi_03 ; pooley_05 ; ihle_05 ; noguchi_08 , that the kinematic shear viscosity of particle-based models, which consist of subsequent streaming and collision steps, is a sum of two terms, namely the kinetic part, $\nu_{\rm kin}$, and the collisional part, $\nu_{\rm coll}$,

[TABLE]

Thus, it is plausible that such a decomposition is also valid for the VM. The kinetic part is due to the momentum that is carried by a particle moving ballistically and can, for example, be calculated by a Boltzmann-like kinetic equation. For the standard VM, this calculation has been done in Refs. 29 and 50, resulting in

[TABLE]

The auxiliary quantity $p$ involves an infinite sum,

[TABLE]

where the coefficients $K_{2c}^{11}$ are given in Table I of Ref. 50. Expression (9) can be evaluated approximately at small and large partner number $M$. For small normalized densities $M\ll 1$, a good approximation is

[TABLE]

In the opposite limit, $M\gg 1$, one finds to leading order:

[TABLE]

In Ref. 50 it was demonstrated that, like in regular fluids, the same expression for $\nu_{\rm kin}$ can be obt‘ined by evaluating a simple Green-Kubo relation by means of the molecular chaos approximation.

The collisional contribution, $\nu_{\rm coll}$ in Eq. (7), stems from collisional transfer of momentum across the finite interaction range $R$, and is therefore outside the scope of Boltzmann-like equations. In contrast, an Enskog-like theory should be able to capture this contribution (see Appendix A). However, in previous calculations ihle_11 ; ihle_16 , a large mean-free path $\Lambda=v_{0}\tau\gg R$ was assumed, where $\nu_{\rm coll}$ becomes negligible. In this manuscript, we show how to calculate $\nu_{\rm coll}$ for the standard VM within mean-field kinetic theory. (Note, that Boltzmann approaches such as those of Refs. 25 and 67 are unable to obtain this important contribution to the viscosity.) The details of this calculation have been moved to the Appendix A for conciseness, and we summarize here only the final result for $\nu_{\rm coll}$

[TABLE]

This term has been neglected in previous publications, and it is clear by dimensional analysis, that the collisional part dominates the viscosity in the typical regime of the VM, such as originally used by Vicsek et al. vicsek_95 , because $\Lambda\ll R$. This is because $\nu_{\rm kin}$ scales with time step $\tau$ and the effective temperature, $k_{B}T/m\sim v_{0}^{2}/2$, whereas $\nu_{\rm coll}$ is proportional to $R^{2}/\tau$, thus $\nu_{\rm coll}/\nu_{\rm kin}\propto(R/\Lambda)^{2}$.

For $M\gg 1$, Eq. (12) can be approximated as

[TABLE]

by means of a saddle point expansion inside the infinite sum of Eq. (12). In the opposite limit $M\ll 1$, we keep only the first terms in the sum and find

[TABLE]

Figure 2 shows the predicted collisional viscosity, Eq. (12), as a function of $M$ in comparison to the two approximations, Eqs. (13) and (14). Interestingly, it turns out that the asymptotic expansion, Eq. (13), is not only excellent for $M\geq 1$ but remains a very good approximation for $M<1$ with an error of around one to two percent. In contrast, the approximative expression Eq. (14) which was obtained by truncating an infinite series becomes very accurate at small $M$ but should not be used for $M>1$.

Figures 3 and 4 show both (kinetic and collisional) contributions to the viscosity as a function of noise, $\eta$, and normalized density, $M$. The kinetic contribution is largest at both small $\eta$ and small $M$, whereas the collisional contribution increases with $M$ and decreases with $\eta$. Figure 5 shows the total viscosity $\nu=\nu_{\rm kin}+\nu_{\rm coll}$ for two particular sets of parameters in comparison with $\nu_{\rm kin}$ and $\nu_{\rm coll}$. Clearly, for these parameters, neglecting the collisional part leads to a large error.

Finally, we consider the system used in Vicsek’s original paper, Ref. 10. Translating the parameters from their Fig. 2(a) into our notation leads to $M=12.57$, $R=1$, and $\Lambda=v_{0}\tau=0.03$. Choosing $\eta=3.5$, which is slightly above $\eta_{\rm c}$, and applying expressions (8) and (12), we predict $\nu_{\rm coll}=1.7$ and $\nu_{\rm kin}=5.4\times 10^{-5}$. This finding confirms the expectation that the kinetic part of the viscosity is negligible here. Of course, these are predictions within the mean-field approximation, which are not expected to be valid at this small ratio $\Lambda/R=0.03$. For improved results, pre-collisional correlations as discussed in Ref. 68 need to be taken into account.

III.2 Hydrodynamic theory for Vicsek-like models

There is shared belief, that on a macroscopic level, polar active systems are described by a minimal set of equations for mass and momentum density – the well-known Toner-Tu equations toner_95 ; toner_98 ; toner_12 . These equations were first postulated on the basis of symmetry and renormalization group arguments. Within the mean-field assumption of molecular chaos, they have also been derived from first principles, for a VM-like model with binary collisions by Bertin et al. bertin_06 ; bertin_09 , and by Ihle ihle_11 ; ihle_14_a for the standard VM with discrete time evolution as considered here. In the latter approach, additional nonlinear gradient terms which were not part of the original Toner-Tu theory, were found ihle_16 .

We would like to apply the Toner-Tu theory to the shear setup given in Fig. 1, in order to extract the values of several transport coefficients from simulation data. Measurements are taken in the stationary state, and thus, all time derivatives in the hydrodynamic equations are set to zero. The stationary state is established by feeding a small amount of $x$-momentum into the top layer of the channel, and by extracting $x$-momentum from the layer in the middle (see Sec. II.3 for details). This procedure leads to a shear flow of small size at not too low noise $\eta$. Therefore, we can neglect most nonlinear terms in the flow velocity. Furthermore, there is no pressure gradient in the $x$-direction. Hence, we assume translational invariance for that direction and neglect all spatial derivatives with respect to $x$. Because of the particular way particle velocities are swapped, no net $y$-momentum is transferred between the feeding layers. Analyzing the continuity equation, $\partial_{t}\rho+\partial_{\alpha}(\rho u_{\alpha})=0$, under the previous assumptions shows that the transversal derivative, $\partial_{y}w_{y}$ of the $y$-component of the momentum density, ${\bf w}=\rho{\bf u}=(w_{x},w_{y})$, should be zero. In shear flow of a regular Newtonian fluid, the density is constant and the transversal velocity vanishes, $u_{y}=0$. This is not quite the case for the active fluid considered here. Instead, due to a lack of Galilean invariance, there are additional convective terms which prevent such a simple shear solution. Nevertheless, agent-based simulations of the VM (see Sec. IV.1) showed that the density variations across the channel are less than one percent, at least outside the parameter range where the well-known density instability of the regular VM occurs bertin_06 ; ihle_11 ; ihle_13 . Therefore, density gradients will be ignored in our theory.

Under these circumstances, the Toner-Tu equations for the components of the macroscopic velocity ${\bf u}=(u_{x},u_{y})$ take a simplified form:

[TABLE]

The kinematic viscosity $\nu$, the coefficients of the convective terms, $\mu_{1}$, $\mu_{2}$, $\mu_{3}$, and the strength of the cubic nonlinearity $q$ depend on the time step $\tau$, density $\rho$, noise $\eta$ and the interaction radius $R$. The main difference to a regular fluid is the linear term in ${\bf u}$ which results from the violation of momentum conservation. The coefficient $\lambda$ describes, whether on average, momentum is lost or gained in a collision. The cubic term, $\propto u^{2}{\bf u}$, becomes relevant below the threshold noise $\eta_{\rm c}$, where $\lambda>1$. The threshold noise is defined by the condition $\lambda(\eta_{\rm c})=1$.

III.2.1 Analysis for the disordered state, $\lambda<1$

In the disordered state, $\eta>\eta_{\rm c}$, the momentum amplification factor $\lambda$ is smaller than one, which means that, on average, momentum is lost in collisions. If momentum is “fed” into the boundary layer, it can only penetrate into the bulk of the channel within a certain distance $l_{S}$ due to the interplay of momentum-diffusion and “-evaporation”. Since this behavior appears to be similar to the skin effect in electrodynamics, $l_{S}$ will be called skin depth. In this scenario, we can neglect the cubic term in Eq. (15). Since at $\eta>\eta_{\rm c}$ there is no spontaneous symmetry breaking, we also neglect the transversal component $u_{y}$. Both assumptions have been justified numerically, and they allow us to obtain an analytical solution for the velocity profile across the channel:

[TABLE]

where $u_{x}$ is the $x$-component of the macroscopic velocity. This profile is to be applied to the upper (or lower) half of the channel with the $y$-coordinate set to zero in the middle of the considered half-channel. The coefficient $d_{1}$ is given by

[TABLE]

As shown in Sec. IV.1 of this paper, velocity profiles from agent-based simulations which were averaged in time and over the length of the channel, show excellent agreement with this ${\rm sinh}$ profile. Fitting data to this profile enables the determination of the constants $d_{0}$ and $d_{1}$.

To recover both transport coefficients $\lambda$ and $\nu$, an additional quantity – the momentum flux – is needed. The momentum flux $\sigma$ is determined by measuring the amount of momentum which is fed into the top layer per time and length in the simulations. This flux is linked to the velocity gradient by

[TABLE]

where $\rho_{m}$ is the mass density, and the gradient is to be evaluated at the top of the channel, at $y=L_{y}/4$ ($y=0$ is defined in the middle of the upper half-channel, see Fig. 1). Inserting the solution, Eq. (17) into Eq. (19) gives an equation for the viscosity $\nu$:

[TABLE]

Note, that using the coefficients $d_{0}$ and $d_{1}$ from a fit of the velocity profile, instead of applying Eq. (19) directly, circumvents the problem of numerically evaluating a velocity gradient in the fluctuating top layer of the channel. Once $\nu$ has been determined, it can be inserted in the relation for $d_{1}$, Eq. (18), yielding an expression for the coefficient $1-\lambda$:

[TABLE]

It is possible to formally integrate Eq. (15) with $\mu_{1}=\mu_{2}=0$, but with the cubic nonlinearity on the right hand side included. However, fitting this solution to numerically obtained velocity profiles failed, in the sense, that it did not give reliable estimates for the coefficient $q$. The reason is that the averaged velocities in our simulation data were too small for the nonlinearity to be relevant. This was verified independently by using the mean-field prediction for this coefficient from Refs. ihle_11 ; ihle_16 , evaluating the cubic term by hand and observing that it is negligible compared to the linear terms.

III.2.2 Analysis for the ordered state, $\lambda>1$

The situation in the ordered state is more complicated than the one at $\lambda<1$ because, (i) at noise values slightly below the threshold noise, soliton-like density waves occur in the regular VM chate_08 ; bertin_06 ; ihle_13 . That means, density gradients are large and derivatives with respect to $x$ cannot be neglected. (ii) Spontaneous symmetry breaking occurs, leading to large macroscopic velocities that are not necessarily parallel to the walls of the channel.

The former issue will be ignored, because there are models such as the metric-free VM cavagna_10 ; ginelli_10a ; chou_12 or the incompressible active liquid chen_15 , where such density waves do not occur. Thus, for simplicity, in our analysis we still omit density gradients. The latter issue means that nonlinear terms are relevant and that, depending on the situation, the transversal component of the velocity $u_{y}$ could be larger than $u_{x}$. Here, as a first step, we assume to be close to the threshold, $\lambda-1\ll 1$. Furthermore, assuming small velocity gradients and small momentum transfer rates, we still ignore the convective nonlinearities but keep the cubic nonlinear term with coefficient $q$ to stabilize the solution. With these considerations, Eqs. (15) and (16) now become,

[TABLE]

where the definition $u_{\alpha}^{\prime}\equiv{\partial u_{\alpha}\over\partial y}$ was used. Multiplying the first equation (22) by $u_{x}^{\prime}$, the second one by $u_{y}^{\prime}$, and adding both equations yields:

[TABLE]

We define the flow velocity of a homogeneous ordered state, $u_{0}=\sqrt{|\kappa|/q}$, and the normal vector $\mathbf{\hat{n}}$ for its flow direction. Near the threshold to collective motion, $\lambda-1\ll 1$, and for a particular constant direction $\mathbf{\hat{n}}$, using Eq. (24), we expand the solution around the homogeneous ordered state, and obtain the approximate result

[TABLE]

where the unit vector $\mathbf{\hat{t}}$ and the constant $A$ are arbitrary. In finite, not too large systems, both directions $\mathbf{\hat{n}}$ and $\mathbf{\hat{t}}$ fluctuate over time. Because our simulation data are time-averaged, such an average is also performed over Eq. (25). For the $x$-component of the flow velocity, one obtains

[TABLE]

Apart from the constant $d_{2}$, the solution has the same form as the one in the disordered state. Note, however, that the coefficients $d_{0}$ and $d_{2}$ originate from the time-average of the fluctuating unit vectors, $d_{0}\equiv\left\langle A\,\hat{t}_{x}\right\rangle$, $d_{2}\equiv\left\langle u_{0}\,\hat{n}_{x}\right\rangle$, and therefore strongly depend on the system parameters and the details of the time-average. Thus, for example, it is possible to observe an averaged flow profile with $d_{2}\approx 0$ which deceivingly looks like the one found in the disordered phase, even though particles have strong orientational order at any given time. The procedures and formulas to obtain the viscosity, Eq. (20), for both the ordered and the disordered phase are identical. However, there is a difference between Eqs. (18) and (26) for the fit parameter $d_{1}$. Because of that, for $\eta<\eta_{\rm c}$ one finds,

[TABLE]

III.3 Transverse current fluctuations

The most common methods of calculating the shear viscosity from simulations are the GK approach and nonequilibrium Molecular Dynamics. A third, less popular, approach is the use of transverse-current auto-correlation functions. This method relies on the fact that in molecular liquids in thermal equilibrium, long-wavelength fluctuations in transverse momentum fields decay exponentially with a decay constant $\nu k^{2}$, where $\mathbf{k}$ is the wave vector of the fluctuation, and $\nu$ is the kinematic shear viscosity. This approach was used, for example, to calculate the shear viscosity for mono-atomic liquids in Ref. 51, and for liquid carbon dioxide, a molecular fluid, in Ref. 52. In an “artificial” fluid without momentum conservation such as the VM, the decay constant should contain an additional term which does not depend on $k$, because even at zero wave number, momentum fluctuations still decay. This additional term should contain information about the momentum amplification factor $\lambda$.

To describe small fluctuations in a stationary state, we start with the linearized Toner-Tu equations for the momentum density ${\bf w}=(w_{x},w_{y})$ which, in the spirit of Landau-Lifschitz’s theory on “fluctuating hydrodynamics” landau_59 , are augmented with a noise source ${\bf\Psi}=(\Psi_{x},\Psi_{y})$,

[TABLE]

Here, $P$ is the pressure, $\kappa\equiv(1-\lambda)/\tau$, and $\sigma_{\alpha\beta}$ is the viscous stress tensor,

[TABLE]

where $d=2$ is the spatial dimension. The bulk viscosity is not included in Eq. (29) because it is irrelevant for the TC fluctuations. Higher order gradient terms and nonlinear terms were neglected in this equation. We also assume that we are in the disordered phase, that is $\kappa\geq 0$. As usual, the average of the noise can be chosen to vanish $\left\langle{\bf\Psi}\right\rangle=0$. However, not much is known about its correlations; in general, we can neither assume that the noise is white nor that its components are uncorrelated.

Defining the Fourier transform of the momentum density,

[TABLE]

Eq. (29) reads in Fourier space as

[TABLE]

where $\hat{P}$ and $\hat{\Psi}_{\alpha}$ are the Fourier transforms of pressure and noise, respectively. The simplest way to model a colored noise is to assume an exponential form

[TABLE]

where $C$ models the unknown (but irrelevant) strength of the noise and $\gamma$ is the decay rate of the noise correlations. This leads to the definition of the memory time of the noise, $\tau_{\rm N}\equiv 1/\gamma$. The tensor $g$, in particular its off-diagonal elements, describe possible correlations between the different spatial components of the noise. For simplicity and for symmetry reasons, one has $g_{xx}=g_{yy}=1$, and $g_{xy}=g^{*}_{yx}$. Without those correlations, i.e. with $g_{xy}=0$, and in the limit $\gamma\to\infty$, the white noise behavior of a regular fluid is recovered, where the correlations become equal to $\delta_{\alpha\beta}C\delta(t-\tilde{t})$. All three quantities $\gamma$, $C$, and $g$ are likely to depend on the wave vector $\mathbf{k}$.

A general result of the Mori-Zwanzig projector operator formalism mori-zwanzig is, that the correlations of the (internal) noise are equal to the memory kernel in the corresponding generalized Langevin-Equation. However, our Langevin-equation, Eq. (32), contains only local terms. Therefore, even though we did not apply this formalism explicitly, for consistency FOOT_MORI we assume a white noise, i.e. $\gamma\rightarrow\infty$ in the following calculations followed by an estimate of when this assumption is likely to fail.

One way to extract the TC fluctuations is to focus on the vorticity $\bm{\omega}\equiv\nabla\times{\bf w}$ of the flow, whose $z$-component is given by

[TABLE]

In Fourier space, Eq. (34) becomes

[TABLE]

Multiplying the $x$-component of Eq. (29) by $k_{y}$ as well as multiplying the $y$-component by $k_{x}$, and subtracting both equation from each other leads to a closed equation for $\Omega$,

[TABLE]

where

[TABLE]

Thus, by using the vorticity, one has managed to effectively “project out” the longitudinal modes which contain pressure and bulk viscosity. The correlations of the noise $\hat{\phi}$ follow from Eq. (33) in the limit $\gamma\rightarrow\infty$ as

[TABLE]

where

[TABLE]

The stochastic differential equation, Eq. (36), is solved, and in the stationary limit, $t\to\infty$ and $\tilde{t}\to\infty$, we obtain the vorticity correlations

[TABLE]

Let us define the viscosity-related decay time $\tau_{\mu}\equiv 1/\mu$ and the ratio of the two characteristic times scales,

[TABLE]

Once the decay time $\tau_{\mu}$ and the memory time of the noise $\tau_{\rm N}$ are approximately equal, memory effects should matter for the decay of the vorticity correlations, and we expect deviations from the white noise prediction, Eq. (40).

In the collision-dominated regime of the VM, where $v_{0}\tau\ll R$, $\nu_{coll}/\nu_{kin}\gg 1$ and $\nu\propto 1/\tau$, the total viscosity, $\nu$, can become very large at small time steps $\tau$. Thus, $\mu$ would also become large. This trend to small $\tau_{\mu}=1/\mu$ is intensified if one is deep in the disordered phase, where $\kappa$ is large too, and also in systems with a small linear dimension $L$ since the smallest useful wavenumber is equal to $2\pi/L$ and thus can be rather large. In the same limit of small time steps, a particle needs more iterations to move out of the collision circle of current collision partners. Hypothesizing, that, at least far in the disordered phase, the memory time of the noise $\tau_{\rm N}$ is approximately given by the time two particles diffuse away from each other by a distance $R$, we obtain a rough estimate for $\tau_{\rm N}$,

[TABLE]

For small time steps, we also have $\mu\propto 1/\tau$, and therefore

[TABLE]

Hence, we predict that at sufficiently small mean free path $v_{0}\tau$ (compared to the radius of the interaction circle), the ratio of the time scales, $\delta$, will become larger than 1, meaning that the memory time of the noise cannot be neglected. Therefore, we apply the TC method only at sufficiently large mean free paths where we can assume decent accuracy of Eq. (40). The corresponding numerical results are presented in Sec. IV.2

IV Numerical results

IV.1 Velocity profiles, polar order, and transverse current fluctuations

We performed two-dimensional agent-based simulations of the regular and the metric-free (topological) VM with system sizes ranging from $16\times 16$ to $128\times 32$, and followed the RP protocol with swap-times between $\Delta t=1$ an $\Delta t=4$ (see Sec. II.3). The simulations usually ran for $5\times 10^{6}$ to $3\times 10^{7}$ iterations after the stationary state has been reached in order to ensure sufficient accuracy in the time averages. Since momentum is not conserved in the VM, the stationary velocity profile across each of the half channels (see Fig. 1) is usually not linear. Hence, unlike as for momentum-conserving fluids such as MPCD (see Appendix B), the shear viscosity cannot simply be obtained as the proportionality factor between the measured velocity gradient and the applied shear stress. Instead, the theory outlined in Sec. III.2 is used to evaluate the simulation data. In particular, the velocity profile was fitted with a ${\rm sinh}$ profile, according to Eq. (27). The extracted fitting coefficients $d_{0}$ and $d_{1}$ were inserted in Eq. (20) to obtain the viscosity $\nu$. If the polar order parameter was above about $0.15$ and the coefficient $d_{2}$ significantly deviated from zero, expression (28) for the ordered state was used to obtain the momentum gain coefficient $\lambda$. Otherwise, Eq. (21) for the disordered state was applied to extract $\lambda$.

Figure 6 shows the measured velocity profiles as a function of height $y$ for the metric-free VM for three different noise values $\eta$ in the disordered phase. Only the lower half of the channel is shown and the velocities inside the bottom and top layer were discarded to obtain a better fit. At the border of these “feeding” layers the profiles change abruptly, as shown by the dotted red line for $\eta=4$, which is included in the plot for illustration but was not used in the fitting procedure. One sees, that the ${\rm sinh}$ function provides a perfect fit, at least for this set of parameters.

We have also investigated whether shearing the system has an effect on the (average) ordering of the particle velocities, quantified by the polar order parameter $v_{\rm a}$ [see Eq. (4)]. Figure 7 shows typical results of $v_{\rm a}$ as a function of noise, $\eta$, at rest and under shear for the regular VM in a $16\times 16$ box with $M=5$, $\tau=2$, $R=1$, and $v_{0}=1$. The data for these three cases are virtually identical, indicating that shear has a negligible impact on the overall ordering of the particles for the applied shear stresses.

To study the transverse current (or vorticity) correlations, we performed agent-based simulations of the regular and the metric-free VM with periodic boundary conditions with system sizes $L_{x}\times L_{y}$ ranging from $16\times 16$ to $64\times 64$. Neither external forces nor the RP-swapping procedure were applied. Once a system reached its stationary state, the momentum density was measured in our simulations through

[TABLE]

which is the Fourier transformation of the microscopic expression for the momentum density of a system of $N$ point particles at a particular position $\mathbf{x}$, given by (see for example Refs. 56)

[TABLE]

Here, $\mathbf{v}_{j}=(v_{j,x},v_{j,y})$ and $\mathbf{r}_{j}=(r_{j,x},r_{j,y})$ are the velocity and position of particle $j$, respectively. Inserting the components of the transformed momentum density, $\mathbf{\hat{w}}=(\hat{w}_{x},\hat{w}_{y})$, into Eq. (35), the $z$-component of the vorticity, denoted by $\Omega$, can be recorded. Typically, the simulations ran for $10^{5}$ to $10^{6}$ iterations after the stationary state has been reached, and $\Omega(\mathbf{k},t)$ was recorded for a set of small wave vectors. After completion of the simulations, the stored time series was used to calculate the time-averaged vorticity correlations, $\left\langle\Omega({\bf k},t)\Omega^{*}({\bf k},\tilde{t})\right\rangle$. We found that in the disordered phase, for small $k$ and for not too small time steps, these fluctuations decayed exponentially, $\propto{\rm exp}(-\mu|t-\tilde{t}|)$, as shown in Fig. 8.

The decay rate $\mu$ was extracted from exponential fits with different values for the wave vector $\mathbf{k}$ (see Fig. 8). We then extracted the transport coefficients $\nu$ and $\lambda$ from these data by fitting the obtained $\mu$ by the theoretical expectation, $\mu=(1-\lambda)/\tau+\nu k^{2}$ [see Eq. (37)]. Figure 9 shows exemplary data for simulations of the metric-free VM at $M=5$, $\tau=2$, $R_{\rm eff}=1$, and $v_{0}=1$ in a quadratic simulation box with $L_{x}=L_{y}=64$.

IV.2 Transport coefficients

Figure 10 shows the total shear viscosity, $\nu$, obtained by the RP and TC methods as a function of noise, $\eta$. The trend of a decreasing viscosity with increasing noise is the same as in the theoretical prediction. However, the measured values lie consistently by about $15\%$ to $18\%$ above the mean-field prediction, given by Eqs. (7), (8), and (12), even at parameter ranges where one naively would expect mean-field theory to hold. This discrepancy is confirmed by the viscosity measurements through the TC method, which agree rather well with the results from the RP simulations at noise values in the disordered phase.

Figure 11 shows the difference $1-\lambda$ ($\lambda=1$ denotes the threshold condition for the order/disorder transition) as a function of noise $\eta$ for different system sizes and time steps. Apart from the values obtained by the RP method, the figure also shows the mean-field prediction and values obtained by the TC method. Although both the scaled density $M=\rho\pi R^{2}=5$ and the ratio between mean free path and interaction radius, $v_{0}\tau/R=1$ and $2$, are rather large, there is a significant deviation between the numerical results and the theory. This discrepancy indicates that mean-field theory gives quite an inaccurate prediction for the threshold condition at these parameters. One notices that doubling the time step, $\tau$, reduces the deviation to mean-field theory only slightly.

Furthermore, the parameter $\lambda$ seems to jump from a positive value to a negative one around some critical value if the noise is decreased. This behavior could be due to the appearance of density waves right at the onset of collective motion, that render the order/disorder transition discontinuous. A quantitative, mean-field theory of this mechanism in the regular VM is presented in Ref. 39. However, density waves only occur in sufficiently large systems; Figure 12 shows the density distribution of particles along the $x$-direction, taken relative to the center of mass of the system for the large $128\times 32$ systems as well as the small $16\times 16$ systems at $\tau=2$. Measurements have been taken for noise values in the ordered ($\eta=3.3<\eta_{\rm c}$) and disordered regime ($\eta=3.8>\eta_{\rm c}$). The large systems exhibit distinct density waves when $\eta<\eta_{\rm c}$, whereas the small systems did not develop any such density waves. Note that applying shear did not have any appreciable effect on the formation of density waves, but the wave fronts appear to be less sharp due to the emerging ${\rm sinh}$-shaped velocity profile (see Fig. 6).

Nevertheless, even in small systems there might be strong density fluctuations and/or transient clusters at the threshold to collective motion that are precursors of the discontinuous phase transition which is observed in larger systems. In order to test the hypothesis that the jump in the measured values of $1-\lambda$ is caused by those density fluctuations and are not artifacts of the RP method, we also performed measurements for the VM with metric-free interactions. As shown in Fig. 13 and comparing with Fig. 11, these jumps are about a factor of three smaller in the metric-free model and are hardly noticeable in the plot. Further note, that mean-field theory underestimates the coefficient $1-\lambda$ also for the metric-free model. However, increasing the time step from $\tau=2$ to $\tau=5.66$ leads to better agreement with mean-field theory, as expected. In general, it appears that the agreement with mean-field theory is better for the metric-free model than for the regular VM.

Still, we cannot completely rule out that the observed jump in $1-\lambda$ might be an artifact of the assumed hydrodynamic theory which we used to evaluate the RP measurements. In that case, the jump could be interpreted as an error bar in the determination of the momentum gain coefficient $\lambda$ or, alternatively, could simply mean that the RP method is not very reliable in the ordered phase. These interpretations are consistent with an alternative measurement of $\lambda$ by the TC method. In particular, comparing the green triangles to the orange triangles in Fig. 11, which correspond to the same set of parameters but different methods, we see good agreement in the disordered phase with a difference between the two curves that is smaller than the jump in the green curve.

Similar agreement between the RP and TC method is seen in Fig. 14, which shows measurements of the viscosity for the metric-free model (see Sec. II.1 for a definition of the model). Although near the threshold noise, $\eta_{\rm c}$, excellent agreement between the TC and the RP method occurs, deviations are observed at larger noise. At these larger noises, the TC method appears to be more accurate than the RP method or at least seems to require less fine-tuning of numerical parameters such as the appropriate thickness of the feeding layers, swap times $\Delta t$ and so on.

To quantify the errors of the RP method we fitted three different sections of the velocity profiles for the parameters of Figs. 13 and 14 at $\eta=4.0$ and for both $\tau=2$ and $\tau=5.66$. This leads to different fitting coefficients $d_{0}$ and $d_{1}$, which consequently lead to different predictions when plugged into Eqs. (20), (21) or (28). For the runs with $\tau=2$ we found rather small errors, about $5\%$ for $\nu$ and $2\%$ for $\lambda$. However, for the larger time step $\tau=5.66$ the errors become huge if sections of the profile are picked for fitting that either only include profile parts from near the center of the sample or only parts from the vicinity of the feeding layers. In this worst case scenario, one obtains a $100\%$ error in both $\nu$ and $\lambda$. However, by comparing the green with the red curve in Fig. 13, the error for $\lambda$ actually appears to be only around $15\%$ to $18\%$ at the largest noises and even smaller for $\nu$. Nevertheless, this serves as a warning that the thickness of the feeding layer and the channel height $L_{y}$ must be chosen carefully and large enough compared to the mean-free path.

Let us now focus on the effect of time step for otherwise identical conditions (purple circles and blue squares in Fig. 11). Here, we observe that the jump, which indicates the order/disorder transition, occurs at a smaller noise value in the system with the smaller time step $\tau=1$. This shift is a well-known effect in the standard VM, which has been reported for instance in Refs. 69; 29. It has been shown that at large mean free paths, corresponding to large time steps, the threshold noise converges to the mean field prediction ihle_11 , which follows from Eq. (67) by setting $\lambda$ equal to one. This mean-field prediction for $\eta_{\rm c}$ does not depend on particle velocity, time step or interaction radius. However, as observed in Ref. 69, the actual value of $\eta_{\rm c}$ becomes smaller by up to a factor between two and three if the mean free path (or $\tau$ in our case) is reduced. This deviation from mean-field theory is attributed to correlation effects which grow at decreasing mean free path. Although a ring-kinetic theory for correlation effects in the standard VM was attempted in Ref. 68, it fell short of explaining the dependence of $\eta_{\rm c}$ on the mean free path.

A similar dependence of the threshold noise on the mean free path can be seen in Fig. 13 which shows $1-\lambda$ for the metric-free model at two different mean free paths, $v_{0}\tau=2$ and $5.66$. Here, the threshold noises are $\eta_{\rm c}\approx 3.24$ and $\eta_{\rm c}\approx 3.63$, respectively. Around these noise values one observes a tiny jump of $1-\lambda$ in the figure. Additionally, one sees that the RP method becomes less accurate deep in the ordered phase but also further in the disordered phase, away from the threshold, see also the orange curve in Fig. 14 for $\tau=2$ and $\eta\geq 3.9$. This is partly because if $\lambda$ differs significantly from one, the velocity profile decays rapidly towards the middle of each half-channel, making a fit by a $\sinh$-function less reliable. Furthermore, at the larger mean free path, $v_{0}\tau=5.66$ and a “feeding layer” smaller than this length (as used in our simulations), the discreteness of the dynamics impacts the velocity profile which deviates from a $\sinh$-function. Fitting it anyway by such a function creates a rather large error. This is especially visible in Fig. 13 at noises $\eta\geq 4.0$ and large $\tau=5.66$.

We noticed a significant $k$-dependence of the viscosity near the order/disorder threshold at certain parameter values. In particular, fitting the decay constant $\mu$ by a function $\kappa+\tilde{\nu}k^{\beta}$ with three free parameters sometimes led to an exponent $\beta$ smaller than two, at least in the range of $k$-values we investigated. We checked that $\beta$ approaches the value two at larger mean free paths and further away from the threshold, i.e. at larger noise. Since this effect seems to be more pronounced than in regular fluids near criticality, we think that it is a result of the velocity alignment interaction. A detailed numerical and analytical investigation of this behavior is beyond the scope of the paper but is subject of current research.

The quantity, $\xi\equiv\sqrt{\nu/\kappa}$ has units of length, and within kinetic theory, it can be established ihle_19 as a mean-field approximation to the correlation length. At a continuous phase transition, the correlation length should diverge with the critical exponent $\bar{\nu}$ as

[TABLE]

with the usual mean-field exponent $\bar{\nu}=1/2$, see for example Ref. 75. Inserting our measured values of $\kappa$ and $\nu$ into the expression for $\xi$, and plotting this as a function of the relative distance to the threshold noise, $(\eta-\eta_{\rm c})/\eta_{\rm c}$, the divergence with the mean-field exponent of $1/2$ is reproduced rather well, as shown in Fig. 15. Of course, this result does not rule out that the actual correlation length diverges with an exponent different from $1/2$. Measuring of actual critical exponents requires careful finite size scaling and is beyond the scope of this paper.

IV.3 Evaluation of Green-Kubo relations

GK relations green_54 ; kubo_57 ; zwanzig_65 ; forster_75 provide a convenient way to measure transport coefficients in equilibrium Molecular Dynamics or other particle-based simulation methods of regular fluids. Typically, these relations are used to obtain the self-diffusion coefficient and the viscosity. The kinetic part of the viscosity describes the convection of the transverse momentum by a particle. More specifically, every particle that moves in the $x$-direction with some velocity $v_{x}$ carries its transverse momentum $p_{y}=mv_{y}$ with it; when it eventually collides with another particle, it has transferred $y$-momentum in the $x$-direction. This mechanism leads to the appearance of the off-diagonal element $\sigma_{xy}^{\rm kin}(t)=m\sum_{j=1}^{N}v_{j,x}(t)v_{j,y}(t)$ of the kinetic stress tensor in the derivation of the corresponding GK relation by the projector-operator method for regular fluids. A similar derivation for the MPCD-fluid can be found in Refs. 76; 77.

It seems plausible that the same mechanism of momentum transport acts also in generalized fluids, such as the VM, that neither respect momentum conservation nor detailed balance. Indeed, it was shown ihle_16 that the analytical evaluation of the usual GK relation for the kinetic part of the viscosity

[TABLE]

within the mean-field assumption of molecular chaos and setting the temperature $k_{B}T$ equal to $mv_{0}^{2}/2$, leads to an expression for $\nu_{\rm kin}$ which is identical to the one obtained from the Chapman-Enskog theory of the VM, Eq. (8). Thus, at least at the mean-field level, the validity of the GK relation, Eq. (47), has been proven. Although a microscopic derivation of the GK relation for the VM has not been performed yet, its correctness beyond mean-field seems likely, and we use it here anyway.

Figure 16 shows the kinetic viscosities measured in direct simulations of the VM without shear gradient by means of Eq. (47) as a function of noise, $\eta$, in comparison with the theoretical (mean-field) expression. For higher noise, where the theory should become more accurate, we find excellent agreement for a large range of particle velocities, $v_{0}$, as shown in Fig. 17. However, even at the threshold to ordered motion, the deviation is only around $15\%$. These results support the validity of the mean-field derivation of the analytical expression for $\nu_{\rm kin}$, Eq. (8), which was first reported in Ref. 29. Understanding and reducing the deviations between kinetic theory and agent-based simulations will be left for future studies.

V Conclusions

By now, one can find many derivations of hydrodynamic equations from the microscopic interactions of active particle systems in the literature. These derivations are often complicated and involve several more or less severe approximations. Therefore, the validity of the obtained expressions is not a priori clear, and it would be useful to verify them. In this manuscript, we developed a hydrodynamic theory for both the standard and the metric-free version of the Vicsek model (VM) of self-propelled particles, when shear is applied through Müller-Plathe’s reverse perturbation (RP) method. Feeding momentum into the boundaries of a channel filled with self-propelled particles led to an almost exponential decay of the flow speed towards the center of the channel due to the lack of momentum conservation. We demonstrated how fitting this decay with an analytical solution of the hydrodynamic equations for the VM allows extracting the two transport coefficients, namely the shear viscosity, $\nu$, and the momentum amplification coefficient, $\lambda$. In order to compare with existing kinetic theories, an improvement of a previous derivation of the viscosity from an Enskog-like kinetic theory was required. This calculation resulted in a new explicit formula for the missing contribution – the collisional part of the viscosity. For a typical choice of parameters from Vicsek’s original paper, we showed that this collisional contribution is larger by a factor of $\approx 10^{4}$ than the previous prediction for the viscosity.

To verify our theory, we performed agent-based simulations of both the standard and the metric-free version of the VM. We measured the transport coefficients $\nu$ and $\lambda$ using two different methods, namely the RP method and the transverse current fluctuation method (TC). In the disordered phase and not too far from the threshold to collective motion, excellent agreement between the measurements of $\nu$ was found. These findings verify our extension of Müller-Plathe’s RP method to active particle systems.

Further, we found reasonable agreement when comparing our measurements of $\nu$ with the predictions from mean-field kinetic theory. However, the measured viscosities were consistently higher by $15\%$ to $18\%$ than the predicted ones. To elucidate the origin of this systematic discrepancy, we also measured the kinetic part of the shear viscosity, $\nu_{\rm kin}$, using the Green-Kubo (GK) approach. In the GK calculations we observed very good agreement between theory and measurements at large noise, close to the maximum noise of $\eta_{\rm max}=2\pi$. However, close to the threshold to collective motion we found a similar difference of about $15\%$ as in the RP and TC measurements. Because most of our measurements were done at rather large time step $\tau\geq 1$, where the viscosity is dominated by its kinetic part, we hypothesize that most of the discrepancy in $\nu$ between theory and simulation is due to the invalidity of the mean-field assumption in the analytical calculation of $\nu_{\rm kin}$. Therefore, future theoretical efforts to improve the expression for the viscosity should focus on this contribution.

The agreement between the results for $\lambda$ obtained by the RP and TC method was also very good but still not as good than for the viscosities. Moreover, we observed that the mean-field prediction for $\lambda$ only became accurate at very large mean-free paths, where there is large mixing of particles and where the assumption of molecular chaos should become valid.

While our results support the correctness of the novel mean-field calculation of the collisional part of the viscosity as well as the validity of earlier results from an Enskog-like kinetic theory, they however underline previous concerns about mean-field assumptions and the relevance of correlation effects in active matter systems. It appears that even at a large average number of interaction partners, $M=5$, and a mean free path that is twice as large as the interaction range, pre-collisional correlations still significantly influence the transport coefficients.

VI Acknowledgments

A.N. acknowledges funding from the German Research Foundation (DFG) under the project number NI 1487/2-1. Computing time was granted on the supercomputer Mogon at Johannes Gutenberg University Mainz (www.hpc.uni-mainz.de).

Appendix A Kinetic theory for the Vicsek model

A.1 Introduction to Enskog-like kinetic theory

In the VM, a given particle $i$ is described by its location ${\bf x}_{i}$ and the angle $\theta_{i}$ of its velocity vector. Hence, the microstate of a system of $N$ particles corresponds to a point in $3N$-dimensional phase space. The time-evolution of the VM in this phase space is Markovian, since information about microstates from earlier times is irrelevant for further evolution. Hence, we can write down an exact evolution equation for the $N$-particle probability density $P$ of the corresponding Markov chain,

[TABLE]

which describes the transition from microscopic state ${\bf A}$ to state ${\bf B}$ during one time step with transition probability $W_{AB}$. The state of the system at time $t+\tau$ is given by the vector, ${\bf B}\equiv(\theta^{(N)},{\bf X}^{(N)})$, where $\theta^{(N)}\equiv(\theta_{1},\theta_{2},\ldots,\theta_{N})$ contains the flying directions of all $N$ particles, and ${\bf X}^{(N)}\equiv({\bf x}_{1},{\bf x}_{2},\ldots,{\bf x}_{N})$ describes all particle positions. The initial microscopic state at time $t$ is denoted as ${\bf A}\equiv(\tilde{\theta}^{(N)},{\bf\tilde{X}}^{(N)})$. The integral over the initial state translates to $\int\,{\rm d}{\bf A}\equiv\prod_{i=1}^{N}\int_{-\pi}^{\pi}\,{\rm d}\tilde{\theta}_{i}\int\,{\rm d}{\bf\tilde{x}}_{i}$, where pre-collisional angles and positions are given by $\tilde{\theta}_{i}$ and ${\bf\tilde{x}}_{i}$, respectively. The transition probability $W_{AB}$ encodes the microscopic collision rules,

[TABLE]

and consists of two parts: the first $\delta$-function describes the streaming step which changes particle positions. The second part contains the periodically continued delta function, $\hat{\delta}(x)=\sum_{m=-\infty}^{\infty}\delta(x+2\pi m)$, which accounts for the modification of angles in the collision step. The particle velocities ${\bf V}^{(N)}\equiv({\bf v}_{1},{\bf v}_{2},\ldots,{\bf v}_{N})$, are given in terms of angular variables $\theta_{i}$,

[TABLE]

For the standard VM, the noise distribution $w_{n}$ is given by

[TABLE]

with noise strength $\eta$. Solving Eq. (48) is intractable without major simplification. The common way to proceed is to use Boltzmann’s molecular chaos approximation by assuming that the particles are uncorrelated just prior to every microscopic interaction foot1 . This approximation amounts to a factorization of the $N$-particle probability into a product of one-particle probabilities, i.e. $P(\theta^{(N)},{\bf X}^{(N)})=\prod_{i=1}^{N}P_{1}(\theta_{i},{\bf x}_{i})$ on the right hand side of Eq. (48). Because molecular chaos neglects pre-collisional correlations, the resulting theory has a mean-field nature. By integrating out all particles except one – the so-called focal particle – in Eq. (48), an Enskog-like equation for the distribution function $f=NP_{1}$ is obtained,

[TABLE]

where $C$ is an Enskog collision operator for multi-particle collisions. In the thermodynamic limit, $N\to\infty$, $L\to\infty$, and $\rho_{0}=N\L^{2}=const.$, this operator is given by

[TABLE]

Here, $\langle\ldots\rangle_{\bf x}=\int_{\odot}\ldots\,{\rm d}{\bf x}_{2}\,{\rm d}{\bf x}_{3}\ldots{\rm d}{\bf x}_{n}$ denotes the integration over all positions of the particles $2,3,\ldots n$ inside the collision circle, and $\langle\ldots\rangle_{\tilde{\theta}}=\int_{0}^{2\pi}\ldots{\rm d}\tilde{\theta}_{1}\,{\rm d}\tilde{\theta}_{2}\ldots{\rm d}\tilde{\theta}_{n}$ refers to the integration over the pre-collisional angles of all $n$ particles inside the circle. The average angle of the focal particle $i=1$, $\Phi_{1}$ is defined in Eq. (2) and is a function of both the pre-collisional angles and the positions of all particles. For more details on the derivation of Eq. (52) and a discussion of the molecular chaos assumption, see Refs. 35 and 50.

A.2 Calculation of the collisional viscosity $\nu_{\rm coll}$

To calculate the collisional viscosity, we will heavily rely on the notations and equations presented in Ref. 50, which are too lengthy to be repeated here in full detail foot2 . There, a Chapman-Enskog expansion (CE) cercignani_88 ; enskog_21 ; chapman_52 , which is basically an elaborated gradient expansion, was constructed to obtain hydrodynamic equations of the VM. To systematically account for gradients in the hydrodynamic fields, a dimensionless ordering parameter $\epsilon$ had been introduced, which was set to unity at the end of the calculation. As a “byproduct” of the CE, expressions for the transport coefficients and the equation of state in terms of microscopic parameters were obtained.

The non-standard CE procedure of Ref. 50 starts with a Taylor expansion of the left hand side of Eq. (52), in which spatial gradients are scaled as $\partial_{\alpha}\to\epsilon\partial_{\alpha}$, and multiple time scales $t_{i}$, whose physical meaning is explained at the end of this Appendix, are introduced in the temporal gradients,

[TABLE]

In addition, the distribution function $f$ and the collision integral, e.g. the right hand side of Eq. (52), are expanded in powers of $\epsilon$,

[TABLE]

In Ref. 50 it was shown that the expansion of the distribution function $f$ in Eq. (55) can be identified as an angular Fourier series,

[TABLE]

with Fourier coefficients $a_{i}$ and $b_{i}$. Thus, the reference state $f_{0}$ of the CE, that is, the leading order contribution to $f$, coincides with the zero mode of the Fourier series.

To obtain Toner-Tu-like equations, the CE expansion has to be performed up to third order in $\epsilon$, given the chosen scaling of Eqs. (54) and (55). Collecting terms in orders of $\epsilon$ leads to a hierarchy of coupled equations for the temporal evolution of $f_{i}$, which are given by Eqs. (22-25) in Ref. 50. These equations contain the microscopic velocity vector, given in Eq. (50).

The goal is to obtain macroscopic equations for the first two moments of $f$, namely the particle density $\rho$ and the momentum density vector ${\bf w}=(w_{x},w_{y})$, which are the “slow” fields in this problem,

[TABLE]

where ${\bf u}=(u_{x},u_{y})={\bf w}/\rho$ denotes the macroscopic flow velocity. To proceed, velocity moments of the hierarchy equations are taken, that is, they are multiplied by products of $e_{x}$ and $e_{y}$ and integrated over the angle $\theta$. This calculation leads to evolution equations for density and momentum, however, split up for the different time scales. For example, there are separate equations for $\partial_{t_{0}}\rho$ and for $\partial_{t_{2}}\rho$. Successively inserting and partially solving the equations, and finally adding all pieces together, for example like $\partial_{t}\rho=\partial_{t_{0}}\rho+\partial_{t_{1}}\rho+\partial_{t_{2}}\rho+\ldots$ ($\epsilon$ has been set to one at this stage) leads to the desired hydrodynamic equations, Eqs. (94) and (130) in Ref. 50.

The microscopic collision rules enter this procedure through the velocity moments of the collision integral $C$, i.e. through quantities like $\langle e_{x}C_{1}\rangle$ or $\langle e_{x}e_{y}C_{2}\rangle$ with $\langle\ldots\rangle\equiv\int_{0}^{2\pi}\ldots{\rm d}\theta$. For example, the former quantity is the $O(\epsilon)$ contribution of the following moment,

[TABLE]

and is defined as

[TABLE]

In these moments of $C\circ f$, a crucial approximation was made in Refs. 29, 35 and 50, that led to the formal absence of collisional contributions to the transport coefficients. This approximation consists of neglecting spatial variations of the distribution $f$ across the interaction circle. This issue comes up because the Enskog-like collision term $C\circ f$ involves integrals with products of $f$ over the collision circle. Here, we abandon this approximation which is not justified if the interaction radius is of the same order or larger than the mean free path, i.e. $R\gtrsim\Lambda$.

Comparing Eqs. (56) and (57) with (58) leads to the identification of the Fourier coefficients $a_{1}$ and $b_{1}$ with the components of the momentum density, ${\bf w}=(a_{1},b_{1})$. Now, inserting $f_{0}$ and $f_{1}$ from Eqs. (56) and (57) into Eq. (59), and performing the integrations yields

[TABLE]

with

[TABLE]

The $n$-dimensional angular integral, $K_{C}^{1}$, has been evaluated before, see table I in Ref. 50.

Expanding the density and the $x$-component of the momentum density around ${\bf x}$ and decorating every spatial gradient with a power of $\epsilon$ gives

[TABLE]

Only the first term from Eq. (64) will contribute to $\langle e_{x}C_{1}\rangle$ since the gradient terms are higher order in $\epsilon$. Thus, we can replace $\int_{\odot}{\rm d}{\bf x}_{2}\,w_{x}({\bf x}_{2})$ by $Aw_{x}({\bf x})$ where $A=\pi R^{2}$ is the area of the collision circle. Inserting the expansion (63) into the defining equation for $M$,

[TABLE]

one finds $M({\bf x})=A\rho({\bf x})+O(\epsilon^{2})$. Thus, $\rho({\bf x})$ can be approximated by $M/A$ in Eq. (61) if one only cares about the first order contribution $\langle e_{x}C_{1}\rangle$. This result is identical to Eqs. (38) and (39) in Ref. 50,

[TABLE]

with the factor $\lambda$

[TABLE]

This factor was discussed in detail in Ref. 50 and it describes the ensemble-averaged amplification of the momentum density. The threshold condition for the transition to collective motion is given by $\lambda=1$ (assuming molecular chaos and a spatially homogeneous system). For $M\gg 1$, Eq. (67) can be approximated as

[TABLE]

whereas for $M\ll 1$ one finds

[TABLE]

Similar to the calculation above, we recalculated moments of the collision operator in second order in $\epsilon$ such as $\langle e_{x}^{2}C_{2}\rangle$ and $\langle e_{x}e_{y}C_{2}\rangle$ without the approximation of large mean free path and again did not see any difference to previous results. Thus, we conclude that, at least at a mean-field level, previous calculations of transport coefficients that depend solely on moments of $C\circ f$ in linear and quadratic order in $\epsilon$ remain correct at small mean free paths. However, in third order in $\epsilon$, additional terms arise that were neglected previously in the large mean free path approximation. Consider the third order contribution to the moment from Eq. (61),

[TABLE]

which, according to Eq. (61) contains an integration of the momentum density over the collision circle, $\int_{\odot}{\rm d}{\bf x}_{2}\,w_{x}({\bf x}_{2})$ where the expansion (64) is inserted, and the integration over the collision circle can be performed explicitly in every term of the series. This calculation gives,

[TABLE]

where $\hat{n}=(\hat{n}_{x},\hat{n}_{y})=(\cos{\alpha},\sin{\alpha})$ is the radial unit vector. Terms with odd powers of $\epsilon$ disappear because of the symmetric (circular) shape of the collision area.

Inserting Eq. (71) into Eq. (61) leads together with Eq. (70) to

[TABLE]

where the coefficients $\Gamma$ and $S$ are given in Eqs. (61) and (62) of Ref. 50, and $a_{2}$ and $b_{2}$ are Fourier coefficients defined in Eq. (57). The new result of the current paper is the third term whose coefficient $H$ is,

[TABLE]

We checked that relaxing the previous restriction on the mean free path only affects the moment $\langle e_{\beta}C_{3}\rangle$ and does not impact other relevant moments, at least in a third-order CE expansion. In order to obtain improved transport coefficients, it therefore suffices to formally replace all occurrences of $\Gamma w_{\beta}w^{2}$ by $\Gamma w_{\beta}w^{2}+H\nabla^{2}w_{\beta}$ in the calculations of Ref. 50 after Eq. (111) of that paper. As a result of this straightforward but technical exercise, we observed that, at least up to third order in $\epsilon$, all transport coefficients except the viscosity remain unchanged. In particular, we found the novel collisional contribution to the kinematic viscosity, $\nu_{\rm coll}$, as presented in Eq. (12) above.

The time scale $t_{0}$, introduced in Eq. (54), is the fast convective time that is associated with the Euler-equation (hence non-dissipative), and it measures the time, momentum is convected due to a pressure gradient. The multitime scale expansion, Eq. (54), allows that an approximation of a given order can vary rapidly with respect to one time scale but more slowly with respect to another. Here, the density $\rho$ does not vary on the time scale $t_{0}$, i.e. $\partial_{t_{0}}\rho=0$. It turns out that the time scale $t_{1}$ is spurious and physically irrelevant in the chosen vicinity to the transition threshold, where we assumed $1-\lambda=\mathcal{O}(\epsilon^{2})$. This is because both hydrodynamic variables, density $\rho$ and momentum density $\mathbf{w}$, do not change at all on this scale, so that $\partial_{t_{1}}\rho=0$ and $\partial_{t_{1}}\mathbf{w}=0$. The timescale $t_{1}$ is only present in the current formalism due to a systematic scaling Ansatz in powers of $\epsilon$, that is, for “historical” reasons. Finally, the time scale $t_{2}$ is a slower relaxation time scale, which is associated with the viscous processes that bring the system into its stationary state. For our chosen scaling, $1-\lambda=\mathcal{O}(\epsilon^{2})$, the local relaxation of momentum due to transferring it to and from the environment (encoded in the alignment interaction and the angular noise), as well as the nonlinear processes that lead to the cubic term $\propto w^{2}\mathbf{w}$ in the hydrodynamic equations, also occur on the same time scale $t_{2}$ as the momentum diffusion. This can be seen in Eq. (114) of Ref. 50.

Appendix B Verification of the RP method for an MPCD fluid

To validate our approach and the implementation of the shear algorithm, we used the RP method to compute the shear viscosity of an MPCD solvent in two dimensions malevanets_99 ; malevanets_00 ; gompper_08 ; kapral_08 ; howard_18 ; howard_19 . MPCD is a particle-based mesoscale technique, often used for simulating the dynamics of complex fluids such as polymeric suspensions and blood flow. The general idea behind MPCD is to adopt a computationally inexpensive, coarse-grained solvent model that faithfully reproduces the solvent-mediated hydrodynamic interactions. Similar as in the VM, the motion of the MPCD particles is governed by alternating streaming and collision steps. During the streaming step, the velocities of all fluid particles are updated according to Eq. (1). In the collision step, the fluid particles undergo stochastic collisions with particles in the same quadratic collision cell, where the edge length of these cells, $a$, dictates the spatial resolution of the hydrodynamic interactions huang:pre:2012 . Here, we employed both the stochastic rotation dynamics (SRD) malevanets_99 and the Andersen thermostat (AT) collision rule allahyarov:pre:2002 . For both variants of the MPCD algorithm, analytic expressions and previous numeric calculations for the transport coefficients are readily available in the literature tuzel_03 ; kikuchi_03 ; ihle_05 ; pooley_05 .

We achieved isothermal conditions in the MPCD-SRD simulation by employing a Monte Carlo style thermostat hecht_05 ; gompper_08 , which correctly conserves the local momentum in each collision cell and reproduces the desired Maxwell velocity distribution. In the MPCD-AT simulations, thermalization was achieved directly through the collision step. Galilean invariance was restored by applying a random shift of the collision cells before every collision step ihle_01 . In all MPCD simulations, the particle mass $m$ was set to unity, and a temperature of $k_{\rm B}T=1$ was used. Simulations were conducted in a quadratic simulation box with $L_{x}=L_{y}=16a$ and periodic boundary conditions in all directions. A particle number density of $\rho=10a^{-2}$ has been used throughout. For the SRD rule, we set the collision angle to $\alpha=110^{\circ}$. We determined the shear viscosity $\eta$ of the MPCD fluids for time steps $\tau=0.1$, $0.2$, $0.4$ and $1.0$, by conducting multiple simulations at different average shear stress $\left\langle\sigma\right\rangle$. Figure 18 shows $\left\langle\sigma\right\rangle$ vs. the measured shear rate, $\dot{\gamma}$ for the MPCD-AT simulations (the MPCD-SRD results are qualitatively similar), demonstrating that the MPCD fluid behaves like a Newtonian liquid, as expected. From these data, the shear viscosity can then be computed as $\nu=\left\langle\sigma\right\rangle/\dot{\gamma}$.

Figure 19 shows $\nu$ as a function of $\tau$ compared to the theoretical prediction for the MPCD-AT and MPCD-SRD algorithm. In both cases, the shear viscosity of the fluid, $\nu=\nu_{\rm coll}+\nu_{\rm kin}$, is dominated at small $\tau$ by the collisional contribution, $\nu_{\rm coll}$. However, as $\tau$ is increased (and thus the mean free path of the particles, $\tau\sqrt{k_{B}T/m}$, becomes larger), particle collisions become less important, and the shear viscosity of the fluid is dominated by the kinetic contribution, $\nu_{\rm kin}$, instead. The viscosity computed from the shear simulations follows these trends perfectly, and we achieved quantitative agreement with the theoretical expressions within $3\,\%$. This small difference of a few percent between theory and simulation is in the same range of errors which were observed previously by using other methods such as GK relations pooley_05 ; gompper_08 ; ihle_05 .