A new universality class describes Vicsek’s flocking phase in physical dimensions (2024)

Abstract

The Vicsek simulation model of flocking together with its theoretical treatment by Toner and Tu in 1995 were two foundational cornerstones of active matter physics. However, despite the field’s tremendous progress, the actual universality class (UC) governing the scaling behavior of Viscek’s “flocking” phase remains elusive. Here, we use nonperturbative, functional renormalization group methods to analyze, numerically and analytically, a simplified version of the Toner-Tu model, and uncover a novel UC with scaling exponents that agree remarkably well with the values obtained in a recent simulation study by Mahault et al.[Phys.Rev.Lett.123, 218001 (2019)], in both two and three spatial dimensions. We therefore believe that there is strong evidence that the UC uncovered here describes Vicsek’s flocking phase.

Using the rescaling convention,

where, without loss of generality, the flocking direction is chosen to be along the $x$-axis,these novel exponents are:

for $d<11/3$ where $d$ is the spatial dimension. Remarkably, the values of these exponents agree very well with the simulation results in both two and three dimensions (falling within the given simulation errors, see Table 1). Therefore, we believe that the new UC uncovered here describes the ordered phase of the Vicsek model.

Simplified Toner-Tu model.—We start with the celebrated TT EOM that describe generic compressible polar active fluids, derived simply from considering the underlying conservation law and symmetries of the system [2, 19, 4]:

where $\rho$ is the mass density field and $\mathbf{g}$ is the momentum density field. Note that instead of using the velocity density field as one of the two hydrodynamic variables in the original formulation [2, 19], we have opted for the momentum field. The physics of course remains the same but this choice has the virtue of simplifying the continuity equation (3) by rendering it linear.In the EOM of $\mathbf{g}$ (A new universality class describes Vicsek’s flocking phase in physical dimensions), all coefficients are generic functions of $\rho$ and $|\mathbf{g}|$, the “pressure” terms $P$’s are functions of $\rho$:

where $\rho_{0}$ is the mean density, and the coefficients $\kappa_{n}$’s and $\nu_{n}$’s are themselves functions of $|\mathbf{g}|$.Furthermore, “h.o.t.” in Eq.(A new universality class describes Vicsek’s flocking phase in physical dimensions) denotes higher order terms in spatial derivatives (e.g., $\nabla^{4}\mathbf{g}$, etc) that are irrelevant to our discussion, and the noise term $\mathbf{f}$ is Gaussian with vanishing mean and statistics:

Linear Theory.—In the flocking phase, the mean magnitude of the momentum field, $g=|\mathbf{g}|$, is nonzero and we are interested in the fluctuating fields around this flocking state:

where hats denote normalized vectors.We now further partition $\delta\mathbf{g}$ into three components that are more natural in our analysis: $\delta\mathbf{g}=\delta\mathbf{g}_{x}+\delta\mathbf{g}_{L}+\delta\mathbf{g}_{T}$, where $\delta\mathbf{g}_{x}=\hat{\mathbf{x}}(\hat{\mathbf{x}}\cdot\delta\mathbf{g})$, $\delta\mathbf{g}_{L}=\hat{\mathbf{q}}_{\bot}(\hat{\mathbf{q}}_{\bot}\cdot$, where$\mathbf{q}_{\bot}$ denotes the wavevector (in spatially transformed Fourier space) perpendicular to the $\mathbf{x}$-direction, i.e.,$\mathbf{q}_{\bot}=\mathbf{q}-q_{x}\hat{\mathbf{x}}$ and $q_{x}=\hat{\mathbf{x}}\cdot\mathbf{q}$. Namely, the three components of $\delta\mathbf{g}$ correspond to its component along the flocking direction, along the direction of the wavevector (with the $\mathbf{x}$-component subtracted), and along the direction perpendicular to both wavevector and flocking direction. Note that Simplification 1 enforces that $\delta\mathbf{g}_{x}=0$ here.

We now analyze the scaling behavior of the ordered phase at the linear level, i.e., by first truncating the TT EOM to linear order in $\delta\rho,\delta\mathbf{g}_{L}$, and $\delta\mathbf{g}_{T}$. The propagators can thus be obtained by inverting the “dynamical matrix” constructed from the linear TT EOM [20]:

where we have defined $\tilde{\mathbf{q}}=(\omega_{q},\mathbf{q})$, ${\cal P}_{L,{ij}}(\mathbf{q})=\hat{q}_{\bot,i}\hat{q}_{\bot,j}$, and ${\cal P}_{T,{ij}}(\mathbf{q})=\delta_{ij}-\hat{x}_{i}\hat{x}_{j}-\hat{q}_{\bot$ is the projector transverse to the $x$ and the longitudinal direction. Further we have defined $\lambda_{g}=\lambda_{1}g_{0}$, $\mu_{\bot}=\mu_{1}$, $\mu_{x}=\mu_{1}+\mu_{3}g_{0}^{2}$ and $\mu_{\bot}^{L}=\mu_{1}+\mu_{2}$, which are all evaluated at $\delta\rho=\delta\mathbf{g}=0$. Since $\delta\mathbf{g}_{x}=0$, ${\bf G}$ is perpendicular to $\hat{\bf x}$.

The equal-time correlation functions can then be obtained in the usual way, giving,

where $\int_{\mathbf{q}}=\int d^{d}\mathbf{q}/(2\pi)^{d}$.In particular, our linear analysis identifies the following scaling exponents $\chi^{\rm lin}_{\rho}=\chi^{\rm lin}=(2-d)/2,z^{\rm lin}=2$, and $\zeta^{\rm lin}=1$,which, as expected, are identical to those in previous works [2, 19, 4].

Nonlinear analysis using FRG.— Applying Simplification 2 to eliminate all nonlinearities involving $\delta\rho$, the only nonlinearities left are terms involving the $\lambda$’s and $U$, which become independent of $\delta\rho$.The standard power counting method (e.g., see [21]) shows that below $d=4$, the leading order contributions of these nonlinearities (i.e., the $\lambda$’s, which are no longer functions of $\delta\mathbf{g}$, and $U=\beta|\delta\mathbf{g}|^{2}/2$) can modify the scaling behavior and thus have to be incorporated into the analysis.RG methods provide a systematic way to accomplish this task and we will use here the functional version of the renormalization group based on the exact Wetterich equation [11, 12, 13]:

where $\Gamma_{k}$ is the wavelength ($k$) dependent effective average actionand $R_{k}$ is a regulator that serves to control the length scale ($\sim 1/k$) beyond which fluctuations are averaged over. The exact flow equation (14) serves to interpolate between the microscopic action $\Gamma_{\Lambda}$ (where all model details are encoded) and the macroscopic effective average action $\Gamma_{0}$, from which the EOM for the averages of the fields can be obtained.The trace is a sum over all degrees of freedom, i.e., over all field indices, wavevectors and frequencies, and $\Gamma_{k}^{(2)}$ is the matrix containing the second order functional derivatives of $\Gamma_{k}$ with respect to the fields. The boundary conditions for $\Gamma_{k}$ described above are enforced by requiring that $R_{\Lambda}\sim\infty$ and $R_{0}=0$. Otherwise, it can be chosen freely, and in principle $\Gamma_{0}$ is independent of that choice. In practice, we typically constrain the form of the microscopic action $\Gamma_{\Lambda}$in order to close the flow equations with a finite number of coefficients, which are now also dependent on $k$.

To proceed with our NPRG analysis, we use the Martin-Siggia-Rose-de Dominicis-Janssen formalism [22, 23, 24, 18], introducing the response fields $\bar{\mathbf{g}}_{\bot}$ and $\bar{\rho}$, to obtain a scalar action that describes our theory at the microscopic scale $\Lambda$. Making all microscopic couplings dependent on $k$ (not written explicitly), we obtain an Ansatz for the scale-dependent effective average action,

where we have defined $\mathbf{g}_{\bot}=\delta\mathbf{g}_{L}+\delta\mathbf{g}_{T}$.We have also introduced the coefficient $\gamma$ to allow for the potential renormalization of the time-derivative term. Due to its linear structure, the “density sector” (proportional to $\bar{\rho}$) does not get renormalized [25, 26], and therefore its coefficients remain unity.

The last ingredient in the FRG formulation is the regulator, which we choose to be [27],

where $\Theta_{\epsilon}$ is a smooth, nonzero function that approaches the Heaviside function in the limit of $\epsilon\rightarrow 0$, which is to be taken at the end of the calculation. Note that the regulator effectively modifies the propagators with a factor independent of frequency and wavenumber in $x$-direction $q_{x}$. With this property, we are able to evaluate all integrals in the trace of Eq.(14) analytically, except the $q_{x}$-integral. We further note that since the regulator has the same structure as the propagator, even though it is frequency dependent, causality is preserved [18].

Integrating the flow equations numerically, we always find a nontrivial stable fixed point. The associated scaling exponents are shown in Fig.1 (blue squares). For dimensions $11/3\lesssim d<4$, the scaling exponents agree with those obtained by Toner and Tu in Refs.[2, 19]:

Intriguingly, below $d\approx 11/3$, the values of the exponents are found to agree with the new formula shown inEq.(2), until for $d\lesssim 2.4$ when the RG flow seems to become divergent. We hypothesize that the divergence is due to our truncation/simplification of the scale-dependent average action, as we reason that the density-dependent couplings could become more important in lower dimensions and potentially stabilize the RG flow.

We interpret our findings as follows. The flocking phase of our simplified TT model is generically described by the TT UC for $11/3<d<4$. Below $d<11/3$, a new stable RG fixed point, with scaling behavior given by Eq.(2), emerges. And comparing the values of the exponents obtained using Eq.(2) to a recent simulation study [5] (red stars in Fig.1), we believe that the UC uncovered here describes the ordered phase of the Vicsek model. Schematics of the RG flows illustrating the stability exchange between the TT UC and the UC described here are shown in Fig.2 in terms of the anomalous dimensions $\eta_{\bot}$ and $\eta_{x}$, defined as the graphical corrections of $\mu_{\bot}$ and $\mu_{x}$,

Unexpectedly, different behaviour is observed for the graphical corrections of $\mu_{x}$,

where $\int_{\tilde{h}}d^{d}hd\omega_{h}/(2\pi)^{(d+1)}$ and $V_{ijn}$ is the three-point momentum density vertex (shown in Ref.[20]). First, contrarily to the other corrections, the limit of large $\bar{\kappa}_{1}$ and $\bar{\lambda}_{g}$ does not commute with the integral over $h_{x}$ in Eq.(19). If the limit is performed before the integral, one would incorrectly conclude that this correction vanishes. Evaluated in the correct order however, we find via numerical evaluation of Eq.(19) that, as both $\bar{\kappa}_{1}$ and $\bar{\lambda}_{g}$ diverge, and $\bar{\lambda}_{1}$, $\bar{\lambda}_{2}$ and $\bar{\lambda}_{3}$ approach a fixed point value (marked by an asterisk), the function $F_{x}$ behaves asymptotically as

where asterisks denote fixed point values and $\bar{F}_{x}$ has the asymptotic properties, $\bar{F}_{x}(\bar{\lambda}_{1},\bar{\lambda}_{2},\bar{\lambda}_{3},0)=0$ and $\lim_{x\rightarrow\infty}\bar{F}_{x}(\bar{\lambda}_{1},\bar{\lambda}_{2},\bar{$ [20]. The flow equations therefore only allow a fixed point for two possible scenarios. Either $\bar{\kappa}_{1}$ tends to infinity fast enough such that $\bar{\lambda}_{g}^{13}\ll\bar{\kappa}_{1}^{7}$, in which case $F_{x}=\eta_{x}=0$ and the TT UC is recovered, or $\bar{\kappa}_{1}$ and $\bar{\lambda}_{g}$ diverge in such a way that the ratio $\bar{\lambda}_{g}^{13}/\bar{\kappa}_{1}^{7}$ approaches a constant value, leading to the new UC uncovered here. Specifically, since neither $\bar{\lambda}_{g}$ nor $\bar{\kappa}_{1}$ receive any graphical corrections in this limit, using the standard rescaling, we can write down an effective flow equation for $\bar{\lambda}_{\kappa}=\bar{\lambda}_{g}^{13}/\bar{\kappa}_{1}^{7}$:

with precisely two fixed points: either $\bar{\lambda}_{\kappa}=0$ or $\bar{\lambda}_{\kappa}={\rm const}$. In the latter case, the exponents have to fulfill the hyperscaling relation:

Together with the scaling relations from the vanishing graphical corrections for $\lambda_{1}$ and $D$ (also obtained in Refs [2, 19]):

we can determine the analytical expressions of the exponents in Eq.(2). Further, our calculation here applies below $d\sim 2.4$, thus extending our numerical FRG results to d=2.

Finally, we can support the scenario regarding the exchange of stabilities between the TT UC and our novel UC (Fig.2) by expanding Eq.(21) around the TT fixed point at $\bar{\lambda}_{\kappa}=0$,

which clearly becomes unstable below $d<d_{c}^{\prime}=11/3$. Further supporting evidence is that the exponents of both TT UC (16) and our UC (2) coincide at $d=11/3$, as expected from such an exchange of fixed point stabilities.