Flocking from a quantum analogy: spin–orbit coupling in an active fluid

In the quantum system, if $\kappa \to \infty$ then $\hat{{ \mathcal H }}$ becomes the two-dimensional Dirac Hamiltonian in the Weyl representation (up to an overall energy scale and shift):

$$\begin{eqnarray}&&{H}_{{\rm{D}}}={\boldsymbol{\sigma }}\cdot {\rm{\nabla }}+m(I-{\sigma }_{z}).\end{eqnarray} \tag{ 10 }$$

When trying to use this equation to build an analogy with classical systems, spatial dimensionality plays an important role; the one-dimensional Dirac equation is obtained by considering Ψ independent of (for example) y. In order to establish a probabilistic interpretation of the spinor, recall from equation (7) that in this case the (one-dimensional) velocity operator is ${\sigma }_{x}$. We take full advantage of this fact by rewriting H_D in the eigenbasis of ${\sigma }_{x}$. Under this change of basis, the Pauli matrices translate as ${\sigma }_{x}\to {\sigma }_{z}$ and ${\sigma }_{z}\to {\sigma }_{x}$. In the new basis, the upper (lower) component of the spinor corresponds to the eigenstate with velocity +1 (−1). Thus, in one dimension a general spinor has the structure

$$\begin{eqnarray}&&{\rm{\Psi }}=\left(\begin{array}{c}{p}_{\to }(x,t)\\ {p}_{\leftarrow }(x,t)\end{array}\right).\end{eqnarray} \tag{ 11 }$$

The interpretation is that ${p}_{\to }(x,t)$ $[{p}_{\leftarrow }(x,t)]$ is the probability density at time t for finding a particle at position x traveling to the right (left). Furthemore, in the new basis the Hamiltonian now reads:

$$\begin{eqnarray}&&{H}_{{\rm{D}}}={\sigma }_{z}{\partial }_{x}+m(I-{\sigma }_{x}).\end{eqnarray} \tag{ 12 }$$

By substituting this Hamiltonian into the equation of motion (9), we obtain the master equations:

$$\begin{eqnarray}&&{\partial }_{t}{p}_{\to }=-{\partial }_{x}{p}_{\to }+m({p}_{\leftarrow }-{p}_{\to }),\end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{\leftarrow }={\partial }_{x}{p}_{\leftarrow }+m({p}_{\to }-{p}_{\leftarrow }).\end{eqnarray} \tag{ 13b }$$

We can therefore interpret the one-dimensional Dirac equation, rotated into imaginary time, as the master equation for the probability distribution of a persistent random walker moving along a line at constant speed and with a turning rate (i.e., a rate for changing direction) given by m. The stochastic process corresponding to such a walker is a Poisson process, and the one-dimensional Dirac equation in imaginary time can be restated as the telegrapher's equations describing this process. This analogy has its origins in the path-integral formulation of the Dirac equation [33], also restated in imaginary time in [34–36].

Given the success in using the Dirac Hamiltonian in one dimension, one could expect that the imaginary-time Dirac equation naturally generalizes to higher dimensions. However, unlike in the one-dimensional case, the two-dimensional Dirac Hamiltonian does not consistently describe the probability density of a single, self-propelled particle: to ensure a physical description, κ must be constrained to be less than $8m$. We now proceed to derive this result.

We now proceed to demonstrate the link between (imaginary-time Schrödinger) equation (9) and the probability densities and currents of self-propelled particles. To do so, we decompose the spinor Ψ into real and imaginary parts:

$$\begin{eqnarray}&&{\rm{\Psi }}=\left(\begin{array}{c}\rho +{\rm{i}}\chi \\ {j}_{x}+{{\rm{i}}{j}}_{y}\end{array}\right),\end{eqnarray} \tag{ 14 }$$

where ρ, χ, j_x, and j_y are real-valued functions of the position ${\boldsymbol{r}}$ and time t (and are independent of θ). With this parametrization, equation (9) becomes

$$\begin{eqnarray}-2{\partial }_{t}\rho & = & {\rm{\nabla }}\cdot {\boldsymbol{j}}-\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}\rho ,\\ -2{\partial }_{t}\chi & = & -{{\rm{\nabla }}}_{\perp }\cdot {\boldsymbol{j}}-\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}\chi ,\\ -2{\partial }_{t}{\boldsymbol{j}} & = & {\rm{\nabla }}\rho -{{\rm{\nabla }}}_{\perp }\chi +2\,m{\boldsymbol{j}}-\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}{\boldsymbol{j}},\end{eqnarray} \tag{ 15 }$$

where we have introduced ${\boldsymbol{j}}\equiv ({j}_{x},{j}_{y})$ and, for any vector ${\boldsymbol{a}}$, ${{\boldsymbol{a}}}_{\perp }\equiv ({a}_{y},-{a}_{x})$. Note that the first of these equations can be interpreted as a continuity equation, with ρ taking the role of a density and with both ${\boldsymbol{j}}$ and ${\rm{\nabla }}\rho$ contributing to the current. Furthermore, we can interpret χ as a gauge degree of freedom for the orientation of the local coordinate frame. We make the simplest choice of gauge: $\chi =0$. Substituting this condition into equation (15), we find

$$\begin{eqnarray}&&-2{\partial }_{t}\rho ={\rm{\nabla }}\cdot {\boldsymbol{j}}-\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}\rho ,\,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\perp }\cdot {\boldsymbol{j}}=0,\,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&-2{\partial }_{t}{\boldsymbol{j}}={\rm{\nabla }}\rho +2\,m{\boldsymbol{j}}-\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}{\boldsymbol{j}}.\end{eqnarray} \tag{ 18 }$$

We check the consistency of our gauge choice by noting that if the initial conditions satisfy equation (17), the evolution given by equations (16) and (18) remains consistent with equation (17). Indeed, we find this consistency condition to hold by applying ${{\rm{\nabla }}}_{\perp }$ to equation (18):

$$\begin{eqnarray}&&-2{\partial }_{t}({{\rm{\nabla }}}_{\perp }\cdot {\boldsymbol{j}})=2\,m({{\rm{\nabla }}}_{\perp }\cdot {\boldsymbol{j}})-\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}({{\rm{\nabla }}}_{\perp }\cdot {\boldsymbol{j}}).\end{eqnarray} \tag{ 19 }$$

In what follows, we parametrize the velocity of self-propulsion via ${\boldsymbol{v}}(\theta )=(\cos \theta ,\sin \theta )$, and use equation (16) to show that the system of equations (16)–(18) is equivalent to a Fokker–Planck equation that describes the dynamics of the probability density $P({\boldsymbol{r}},\theta )\equiv \rho +{\boldsymbol{v}}(\theta )\cdot {\boldsymbol{j}}$. Physically, P describes the probability of having a particle near ${\boldsymbol{r}}$ and oriented at an angle near θ. We first decompose the current into components parallel to and perpendicular to the velocity ${\boldsymbol{v}}$ via ${\boldsymbol{j}}=({\boldsymbol{v}}\cdot {\boldsymbol{j}}){\boldsymbol{v}}+({{\boldsymbol{v}}}_{\perp }\cdot {\boldsymbol{j}}){{\boldsymbol{v}}}_{\perp }$ and substitute this identity into the continuity equation (16). To get the dynamics of the distribution of the orientation angle θ, we multiply equation (18) by ${\boldsymbol{v}}$. We add this equation for the time evolution of the current to the continuity equation to find an equation for the probability density P:

$$\begin{eqnarray}&&-2{\partial }_{t}P=({\boldsymbol{v}}\cdot {\rm{\nabla }})P+({{\boldsymbol{v}}}_{\perp }\cdot {\rm{\nabla }})({{\boldsymbol{v}}}_{\perp }\cdot {\boldsymbol{j}})+2m({\boldsymbol{v}}\cdot {\boldsymbol{j}})-\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}P.\end{eqnarray} \tag{ 20 }$$

The two terms that are not yet expressed in terms of P can be addressed in the following way. First, note that ${\boldsymbol{v}}\cdot {\boldsymbol{j}}$ encodes orientational diffusion: ${\partial }_{\theta }^{2}{\boldsymbol{v}}=-{\boldsymbol{v}}$, and thus

$$\begin{eqnarray}&&{\boldsymbol{v}}\cdot {\boldsymbol{j}}=-{\partial }_{\theta }^{2}({\boldsymbol{v}}\cdot {\boldsymbol{j}})=-{\partial }_{\theta }^{2}P.\end{eqnarray} \tag{ 21 }$$

We also find an extra component to diffusion that couples translational noise, rotational noise, and convection. This can be obtained using the identity ${{\boldsymbol{v}}}_{\perp }\cdot {\boldsymbol{j}}=-{\partial }_{\theta }P$:

$$\begin{eqnarray}&&({{\boldsymbol{v}}}_{\perp }\cdot {\rm{\nabla }})({{\boldsymbol{v}}}_{\perp }\cdot {\boldsymbol{j}})=({\boldsymbol{v}}\cdot {\rm{\nabla }})P-{\partial }_{\theta }{\partial }_{x}(P\sin \theta )+{\partial }_{\theta }{\partial }_{y}(P\cos \theta ).\end{eqnarray} \tag{ 22 }$$

By substituting all of the diffusive and convective terms in equations (21) and (22) into (20), we arrive at the Fokker–Planck equation for P:

$$\begin{eqnarray}&&{\partial }_{t}P=-({\boldsymbol{v}}\cdot {\rm{\nabla }})P+\displaystyle \frac{1}{2}\left[{\partial }_{\theta }\,{\partial }_{x}(P\sin \theta )+{\partial }_{\theta }\,{\partial }_{y}(-P\cos \theta )+2m\,{\partial }_{\theta }^{2}P+\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}P\right].\end{eqnarray} \tag{ 23 }$$

From this formulation, we may note that m plays the role of a diffusion constant for the orientation θ, and ${\kappa }^{-1}$ plays that role for the position ${\boldsymbol{r}}$. The derivation of equation (23) is one of our main results: we showed that the Fokker–Planck equation (23) is equivalent to the imaginary-time Schrödinger equation (9) with the Hamiltonian (8), which includes a spin–orbit coupling term.

We now derive the precise microscopic model that corresponds to the Fokker–Planck equation (23). Before doing so, first note that for a Fokker–Planck equation to represent a stochastic microscopic model, the associated diffusion matrix must be positive-definite. In this subsection, we derive and analyze the drift vector and diffusion matrix for the model defined by equations (8) and (9), and derive the conditions under which an underlying microscopic model exists. To begin this analysis, we first define $z\equiv a\theta$, where a is a particle lengthscale, and we re-express θ in terms of z to ensure that the different entries in the diffusion matrix have the same dimensionality. We compare equation (23) with the usual Fokker–Planck equation, viz.,

$$\begin{eqnarray}&&{\partial }_{t}P=-\displaystyle \sum _{i=1}^{3}{\partial }_{i}({\mu }_{i}P)+\displaystyle \frac{1}{2}\displaystyle \sum _{i,j=1}^{3}{\partial }_{i}{\partial }_{j}({D}_{{ij}}P),\end{eqnarray} \tag{ 24 }$$

in which ${\boldsymbol{\mu }}$ is the drift vector and D is the diffusion matrix. We then read off as follows:

$$\begin{eqnarray}&&{\boldsymbol{\mu }}=({\boldsymbol{v}},0)=(\cos \theta ,\sin \theta ,0),\end{eqnarray} \tag{ 25 }$$

using a vector notation in which the third component corresponds to θ and

$$\begin{eqnarray}D=\left(\begin{array}{ccc}1/\kappa & 0 & (a/2){\rm{\sin }}\,\theta \\ 0 & 1/\kappa & -(a/2){\rm{\cos }}\,\theta \\ (a/2){\rm{\sin }}\,\theta & -(a/2){\rm{\cos }}\,\theta & 2\,{{ma}}^{2}\end{array}\right).\end{eqnarray} \tag{ 26 }$$

As mentioned above, in order to construct a particle-based model for equation (23), it is necessary that D be positive-definite. To check this, we examine its eigenvalues $({\lambda }_{1},{\lambda }_{+},{\lambda }_{-})$:

$$\begin{eqnarray}&&{\lambda }_{1}=1/\kappa ,\,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\lambda }_{\pm }=\displaystyle \frac{1}{2}\left(2\,{{ma}}^{2}+\displaystyle \frac{1}{\kappa }\right)\pm {[{(2{{ma}}^{2}-1/\kappa )}^{2}+{a}^{2}]}^{-1/2}.\end{eqnarray} \tag{ 28 }$$

κ, a, and m real and positive guarantees that ${\lambda }_{1}$ and ${\lambda }_{+}$ are positive. From the expression for ${\lambda }_{-}$ we conclude that D is positive-definite if and only if

$$\begin{eqnarray}&&8m\gt \kappa .\end{eqnarray} \tag{ 29 }$$

Thus, in order for equation (23) with $m\gt 0$ to correspond to a microscopic model, translational noise must be present, because without translational noise (i.e., in the limit $\kappa \to \infty$) there isn't an orientational noise parameter m that satisfies equation (29). This implies that, in two dimensions, the Dirac equation alone cannot describe a self-propelled particle, in contrast to the one-dimensional case.

From equation (23), further insight into the relationship between m and κ and their physical interpretations can be gained by deriving the connection between this Fokker–Planck equation and the underlying microscopic process, i.e., the stochastic Langevin equation

$$\begin{eqnarray}&&{\rm{d}}{{\boldsymbol{R}}}_{t}={\boldsymbol{\mu }}({{\boldsymbol{R}}}_{t},t){\rm{d}}{t}+{\rm{\Sigma }}({{\boldsymbol{R}}}_{t},t){\rm{d}}{{\boldsymbol{W}}}_{t}.\end{eqnarray} \tag{ 30 }$$

In equation (30), ${{\boldsymbol{R}}}_{t}$ has N components (corresponding to the random variables), ${\boldsymbol{\mu }}$ is an N-component associated drift vector, ${\rm{\Sigma }}({{\boldsymbol{R}}}_{t},t)$ is an N × M matrix, and ${{\boldsymbol{W}}}_{t}$ is an M-dimensional Wiener process interpreted in either the Itô or Stratonovich sense [37, 38]. For example, in the present case of N = 3, ${{\boldsymbol{R}}}_{t}$ consists of the position vector ${\boldsymbol{R}}=(x,y)$ and orientation angle θ.

Formally, this interpretation can be established for an Itô process by considering a diffusion matrix D of the form ${\rm{\Sigma }}{{\rm{\Sigma }}}^{T}$. Notice that if one such ${{\rm{\Sigma }}}_{0}$ furnishes this decomposition then, for any orthogonal matrix R, ${{\rm{\Sigma }}}^{{\prime} }={{\rm{\Sigma }}}_{0}\,R$ also satisfies it, so there are many different Langevin equations that yield the same Fokker–Planck equation. By using a Cholesky decomposition [39], we find a particular solution for the case M = 3:

$$\begin{eqnarray*}{\rm{\Sigma }}=\left(\begin{array}{ccc}{\kappa }^{-1/2} & 0 & 0\\ 0 & {\kappa }^{-1/2} & 0\\ \tfrac{1}{2}\,a{\kappa }^{1/2}\sin \theta & -\tfrac{1}{2}\,a{\kappa }^{1/2}\cos \theta & \sqrt{2}\,a\,{\left[m-\tfrac{\kappa }{8}\right]}^{1/2}\end{array}\right),\end{eqnarray*}$$

which yields the following Langevin equation:

$$\begin{eqnarray}&&{\rm{d}}{\boldsymbol{R}}={\boldsymbol{v}}(\theta ){\rm{d}}{t}+{\kappa }^{-1/2}{\boldsymbol{\xi }}\,{\rm{d}}{t},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{\rm{d}}\theta =\displaystyle \frac{1}{2}{\kappa }^{1/2}({{\boldsymbol{v}}}_{\perp }(\theta )\cdot {\boldsymbol{\xi }}){\rm{d}}{t}+\sqrt{2}{\left[m-\displaystyle \frac{1}{8}\kappa \right]}^{1/2}{\xi }_{3}\,{\rm{d}}{t}.\end{eqnarray} \tag{ 32 }$$

Here, ${\boldsymbol{\xi }}=({\xi }_{1},{\xi }_{2})$ is the two-dimensional translational noise that acts on the position of the particle, whereas ${\xi }_{3}$ is a rotational noise influencing the polarization angle θ. Note that to interpret this microscopic model we have assumed that equations (31), (32) are Itô stochastic differential equations. Generally, this differs from a Stratonovich process by an extra, noise-induced, drift vector having components

$$\begin{eqnarray}&&{\mu }_{i}=\displaystyle \frac{1}{2}\displaystyle \sum _{k,j=1}^{3}({\partial }_{j}{{\rm{\Sigma }}}_{{ik}}){{\rm{\Sigma }}}_{{jk}}.\end{eqnarray} \tag{ 33 }$$

However, in the case we are considering, the corresponding term is identically zero, and thus equations (31) and (32) can also be seen as a Stratonovich stochastic differential equation.

Let us now discuss the physical picture of the single-particles dynamics described by equations (31) and (32). This microscopic model has similarities to the models of active particles used, e.g., in [1, 40]. At each instant in time, a particle is oriented at an angle θ and attempts to propagate in this direction at a constant speed. However, translational noise can change the direction of propagation away from the particle polarization. As a unique feature, the model we consider has feedback between translational and rotational noise: the larger the translational noise, the weaker the rotational noise. Quantitatively, if we define ${\alpha }_{\xi }$ to be the angle that the force from the translational noise, ${\boldsymbol{\xi }}=({\xi }_{1},{\xi }_{2})$, makes with respect to the x-axis, we can rewrite the rotational noise term ${\kappa }^{1/2}({{\boldsymbol{v}}}_{\perp }(\theta )\cdot {\boldsymbol{\xi }})/2$ in equation (32) as

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{\kappa }^{1/2}| {\boldsymbol{\xi }}| \sin (\theta -{\alpha }_{\xi }).\end{eqnarray} \tag{ 34 }$$

From equation (34), we observe that particles in effect try to oppose the translational noise, and prefer to align opposite to the direction of each kick. This coupling acts as a guidance system: in the absence of translational noise, the particle does not know which way to point. This is a consequence of how κ enters equation (32): the translational noise strength is inversely proportional κ, whereas the rotational noise strength is proportional to κ. Thus, the particle depends on feedback from translational noise to decide where to go. (A curious analogy emerges from the physics of hair cells in the inner ear, which depend on the presence of external noise to complete their function [41].)

To further examine this feedback feature, consider the extreme case in which m is only slightly bigger than the lower-bound of $\kappa /8$, i.e., $m=\kappa /8+\delta$, with $\delta /\kappa \ll 1$. If δ is sufficiently small then rotational noise becomes irrelevant, compared to the large translational noise, and equation (32) becomes

$$\begin{eqnarray}&&{\rm{d}}\theta =\displaystyle \frac{1}{2}{\kappa }^{1/2}| {\boldsymbol{\xi }}| \,\sin (\theta -{\alpha }_{\xi }){\rm{d}}{t}.\end{eqnarray} \tag{ 35 }$$

In this regime, the noise dominates over the self-propelled aspect of particle motion. The strength of this noise, quantified by ${\kappa }^{-1}$, is not subject to any restrictions, and both large- and small-noise regimes are physically accessible.

Curiously, the interplay between the strength of the translational noise and its effect on the polarization is an expression of an uncertainty principle in this model. To see this, disregard the drift, and consider the feedback on the angle as simple additive noise. One then finds $\langle {[{\boldsymbol{R}}(t)-{\boldsymbol{R}}(0)]}^{2}\rangle \sim 4t/\kappa$ and $\langle {[\theta (t)-\theta (0)]}^{2}\rangle \sim t\kappa /2$, which suggest the relation:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{t}^{2}}\langle {[{\boldsymbol{R}}(t)-{\boldsymbol{R}}(0)]}^{2}\rangle \langle {[\theta (t)-\theta (0)]}^{2}\rangle \sim 2.\end{eqnarray} \tag{ 36 }$$

This is a direct analog of the Heisenberg uncertainty principle, which relates the uncertainty of the position and velocity (here captured by orientation θ) of a quantum particle.

We now compare this microscopic model with others discussed in the literature. One related example involves a system composed of self-propelled hard rods that also experience translational noise, as examined in [42–44]. In these models, the translational and orientational noises are assumed to be uncorrelated. A situation closer to ours is explored in [45], in which the translational noise affects both orientational and spatial diffusion. In that case, the effects of these correlations have been examined in the inertial regime, in which Fokker–Planck dynamics are not equivalent to the imaginary-time Schrödinger equation that we examine here.

To conclude this section, let us generalize this interplay between translational and rotational noise and give it an arbitrary strength. In this case, equation (32) acquires an additional arbitrary (real) parameter λ via:

$$\begin{eqnarray*}&&{\rm{d}}\theta =\displaystyle \frac{1}{2}{\kappa }^{1/2}\lambda | {\boldsymbol{\xi }}| \sin (\theta -{\alpha }_{\xi }){\rm{d}}{t}+\sqrt{2}{\left[m-\displaystyle \frac{\kappa {\lambda }^{2}}{8}\right]}^{1/2}{\xi }_{3}\,{\rm{d}}{t}\end{eqnarray*}$$

(with $m\gt \kappa {\lambda }^{2}/8$). This parameter λ controls the response of a self-propelled particle to translational noise. The sign of λ determines the type of response: for $\lambda \lt 0$, the particle turns in the direction of any translational kick, whereas for $\lambda \gt 0$, as in the case above, the particle reacts in opposition to the kick. The Fokker–Planck equation associated with the Langevin dynamics of equation (31) is given by

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}P & = & -({\boldsymbol{v}}\cdot {\rm{\nabla }})P+\displaystyle \frac{1}{2}\lambda [{\partial }_{\theta }{\partial }_{x}(P\sin \theta )-{\partial }_{\theta }{\partial }_{y}(P\cos \theta )]+m{\partial }_{\theta }^{2}P+\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}P.\\ & & \end{array}\end{eqnarray} \tag{ 37 }$$

Although the additional parameter λ introduces more flexibility into the model, it destroys the uncertainty principle (36) and the bridge with the Schrödinger equation describing a two-component spinor.

Instead of trying to exactly solve equation (44), in the present section we introduce a self-consistent approximation. To do this, we rewrite the potential as

$$\begin{eqnarray}&&{V}_{{\rm{i}}}(\{{\boldsymbol{r}},\theta \})=g\,{\mathfrak{I}}(h({{\boldsymbol{r}}}_{{\rm{i}}}){{\rm{e}}}^{{\rm{i}}\alpha ({{\boldsymbol{r}}}_{{\rm{i}}})}\,{{\rm{e}}}^{-{\rm{i}}{\theta }_{{\rm{i}}}}),\end{eqnarray} \tag{ 45 }$$

where

$$\begin{eqnarray}&&h({{\boldsymbol{r}}}_{{\rm{i}}}){{\rm{e}}}^{{\rm{i}}\alpha ({{\boldsymbol{r}}}_{{\rm{i}}})}\equiv \displaystyle \sum _{j(\ne {\rm{i}})}R({{\boldsymbol{r}}}_{{\rm{i}}}-{{\boldsymbol{r}}}_{j}){{\rm{e}}}^{{\rm{i}}{\theta }_{j}}.\end{eqnarray} \tag{ 46 }$$

The defining assumption of the self-consistent approximation is that

$$\begin{eqnarray}&&h({{\boldsymbol{r}}}_{{\rm{i}}}){{\rm{e}}}^{{\rm{i}}\alpha ({{\boldsymbol{r}}}_{{\rm{i}}})}\approx {\langle h({{\boldsymbol{r}}}_{{\rm{i}}}){{\rm{e}}}^{{\rm{i}}\alpha ({{\boldsymbol{r}}}_{{\rm{i}}})}\rangle }_{C},\end{eqnarray} \tag{ 47 }$$

where the average $\langle \cdots {\rangle }_{C}$ is taken with respect to the conditional probability that one of the particles is at position ${{\boldsymbol{r}}}_{i}$ and oriented along ${\theta }_{i}$:

$$\begin{eqnarray}&&P({\{{{\boldsymbol{r}}}_{j},{\theta }_{j}\}}_{j(\ne i)}| {{\boldsymbol{r}}}_{i},{\theta }_{i})=\displaystyle \prod _{{\ell }(\ne i)}{P}_{{\ell }}({{\boldsymbol{r}}}_{{\ell }},{\theta }_{{\ell }}).\end{eqnarray} \tag{ 48 }$$

This approximation treats the inter-particle interaction as an external potential due to the average effect of all the other particles. An explicit computation of the conditional average in equation (47) leads to

$$\begin{eqnarray}&&{\langle h({{\boldsymbol{r}}}_{{\rm{i}}}){{\rm{e}}}^{{\rm{i}}\alpha ({{\boldsymbol{r}}}_{{\rm{i}}})}\rangle }_{C}=\displaystyle \sum _{j(\ne {\rm{i}})}{\int }_{0}^{2\pi }{\rm{d}}{\theta }_{j}\,{{\rm{e}}}^{{\rm{i}}{\theta }_{j}}{\langle R({{\boldsymbol{r}}}_{{\rm{i}}}-{{\boldsymbol{r}}}_{j})\rangle }_{{\theta }_{j}},\end{eqnarray} \tag{ 49 }$$

where

$$\begin{eqnarray}&&{\langle R({{\boldsymbol{r}}}_{i}-{{\boldsymbol{r}}}_{j})\rangle }_{{\theta }_{j}}\equiv {\int }_{A}{{\rm{d}}}^{2}{r}_{j}\,R({{\boldsymbol{r}}}_{i}-{{\boldsymbol{r}}}_{j}){P}_{j}({{\boldsymbol{r}}}_{j},{\theta }_{j})\end{eqnarray} \tag{ 50 }$$

and the integral is taken over the two-dimensional area A. For simplicity, we now consider a purely local interaction [that is to say, taking $R({{\boldsymbol{r}}}_{i}-{{\boldsymbol{r}}}_{j})\to \delta ({{\boldsymbol{r}}}_{i}-{{\boldsymbol{r}}}_{j})$], in which case the self-consistency condition reduces to the simple form:

$$\begin{eqnarray}&&{\langle h({{\boldsymbol{r}}}_{{\rm{i}}}){{\rm{e}}}^{{\rm{i}}\alpha ({{\boldsymbol{r}}}_{{\rm{i}}})}\rangle }_{C}=\displaystyle \sum _{j(\ne {\rm{i}})}{\int }_{0}^{2\pi }{\rm{d}}{\theta }_{j}\,{{\rm{e}}}^{{\rm{i}}{\theta }_{j}}{P}_{j}({{\boldsymbol{r}}}_{{\rm{i}}},{\theta }_{j}).\end{eqnarray} \tag{ 51 }$$

Within the self-consistent approximation, all particles are identical and experience the same forcing. This forcing is, in turn, determined by considering the effect of a particle on its neighbors. The assumption of identical particles leads to all particles having the same probability distributions for all observables. In terms of probabilities, we thus have ${P}_{j}({\boldsymbol{r}},\theta )=P({\boldsymbol{r}},\theta )$ for all j. By using equation (51) and taking the $N\gg 1$ limit, we obtain

$$\begin{eqnarray}&&{\langle h({{\boldsymbol{r}}}_{{\rm{i}}}){{\rm{e}}}^{{\rm{i}}\alpha ({{\boldsymbol{r}}}_{{\rm{i}}})}\rangle }_{C}=N{\int }_{0}^{2\pi }{\rm{d}}\theta \,{{\rm{e}}}^{{\rm{i}}\theta }{P}_{j}({{\boldsymbol{r}}}_{{\rm{i}}},\theta ).\end{eqnarray} \tag{ 52 }$$

Substituting equation (52) into the expression for the potential, we find that the self-consistent potential has the form:

$$\begin{eqnarray}&&{V}_{{\rm{SC}}}({{\boldsymbol{r}}}_{i},{\theta }_{i})=g\,h({{\boldsymbol{r}}}_{i})\sin [\alpha ({{\boldsymbol{r}}}_{i})-{\theta }_{i}].\end{eqnarray} \tag{ 53 }$$

For convenience, we rewrite this expression using an external alignment field ${\boldsymbol{h}}({\boldsymbol{r}})$, defined via

$$\begin{eqnarray}&&{\boldsymbol{h}}({\boldsymbol{r}})\equiv h({\boldsymbol{r}}){\boldsymbol{v}}(\alpha ).\end{eqnarray} \tag{ 54 }$$

In terms of ${\boldsymbol{h}}$, the potential has the form

$$\begin{eqnarray}&&{V}_{{SC}}({{\boldsymbol{r}}}_{i},{\theta }_{i})=g\,{{\boldsymbol{h}}}_{\perp }({{\boldsymbol{r}}}_{i})\cdot {\boldsymbol{v}}(\theta ).\end{eqnarray} \tag{ 55 }$$

The self-consistent alignment field satisfies $| {\boldsymbol{h}}({\boldsymbol{r}})| =h({\boldsymbol{r}})$ and ${\boldsymbol{h}}({\boldsymbol{r}})=\pi N{\boldsymbol{j}}$, i.e., it is a measure of the spontaneous alignment between the particle velocities. The advantage of using this self-consistent approximation is that it reduces the many-body Fokker–Planck equation to the one-particle nonlinear equation, i.e.,

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}P & = & -({\boldsymbol{v}}\cdot {\rm{\nabla }})P-g{\partial }_{\theta }[({{\boldsymbol{h}}}_{\perp }\cdot {\boldsymbol{v}})P]+\displaystyle \frac{1}{2}[{\partial }_{\theta }{\partial }_{x}(P\sin \theta )-{\partial }_{\theta }{\partial }_{y}(P\cos \theta )+2\,m{\partial }_{\theta }^{2}P+\displaystyle \frac{1}{\kappa }{{\rm{\nabla }}}^{2}P].\\ & & \\ & & \end{array}\end{eqnarray} \tag{ 56 }$$

In the present subsection we have restricted ourselves to considering a description of interacting, self-propelled particles in terms of the probability density P_N rather than in terms of the Schrödinger equation. This is done out of necessity: the external potential term $g\,{\partial }_{\theta }[({{\boldsymbol{h}}}_{\perp }\cdot {\boldsymbol{v}})P]$ in equation (56) cannot be captured within the Hamiltonian (8). To demonstrate this impossibility within a concrete example, let us consider a term in the Hamiltonian of the form ${{\boldsymbol{h}}}_{\perp }\cdot {\boldsymbol{\sigma }}$ as a possible candidate. Such a term presents two issues that cannot be overcome within the framework we are considering: (i) Such a term generates a nonzero value of χ (in the imaginary part of the spinor). This issue can be overcome if one considers a more general framework in which the space of quantum states includes four-spinors with the structure ${\rm{\Psi }}=(\phi ,\bar{\phi })$ as well as by including additional terms in the Hamiltonian (8). (ii) More significantly, the external potential term in equation (56) couples the lowest two Fourier modes of the orientation to higher Fourier modes. As a result, a description based on only the first two modes does not form a closed system of equations. We thus conclude that, in general, the Schrödinger equation in imaginary time with a spin–orbit coupling term describes single-particle dynamics only.

Although the spinorial description works for single-particle dynamics only, we can use the Fokker–Planck description to examine the stability of the interacting isotropic active gas. In this subsection, we explore the onset of alignment due to inter-particle interactions. We follow the standard approach based on the dynamics of Fourier modes of the distribution of orientations θ [21]. First, we expand the single-particle probability density in Fourier modes:

$$\begin{eqnarray}&&P({\boldsymbol{r}},\theta )=\rho +{\boldsymbol{j}}\cdot {\boldsymbol{v}}+\displaystyle \sum _{n\geqslant 2}{{\boldsymbol{j}}}_{n}\cdot {{\boldsymbol{v}}}_{n},\end{eqnarray} \tag{ 57 }$$

where ${{\boldsymbol{v}}}_{n}\equiv (\cos [n\theta ],\sin [n\theta ])$ and ${{\boldsymbol{j}}}_{n}$ are the vectors whose components are the distinct Fourier modes of P, i.e.,

$$\begin{eqnarray}&&{j}_{n,x}({\boldsymbol{r}})=\displaystyle \frac{1}{\pi }{\int }_{0}^{2\pi }{\rm{d}}\theta \,\cos (n\theta )P({\boldsymbol{r}},\theta ),\end{eqnarray} \tag{ 58a }$$

$$\begin{eqnarray}&&{j}_{n,y}({\boldsymbol{r}})=\displaystyle \frac{1}{\pi }{\int }_{0}^{2\pi }{\rm{d}}\theta \,\sin (n\theta )P({\boldsymbol{r}},\theta ).\end{eqnarray} \tag{ 58b }$$

Substituting equations (58a), (58b) into (56) and using the linear independence of the Fourier components leads to the following set of coupled equations describing the time-evolution of the 3 lowest Fourier modes:

$$\begin{eqnarray}&&{\partial }_{t}\rho =-\displaystyle \frac{1}{2}{\rm{\nabla }}\cdot \left({\boldsymbol{j}}-\displaystyle \frac{1}{\kappa }{\rm{\nabla }}\rho \right),\end{eqnarray} \tag{ 59a }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{x}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-m\right){j}_{x}-\displaystyle \frac{1}{2}{\partial }_{x}\rho -\displaystyle \frac{3}{4}{\rm{\nabla }}\cdot {{\boldsymbol{j}}}_{2}+g\left({h}_{x}\rho -\displaystyle \frac{1}{2}{\boldsymbol{h}}\cdot {{\boldsymbol{j}}}_{2}\right),\end{eqnarray} \tag{ 59b }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{y}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-m\right){j}_{y}-\displaystyle \frac{1}{2}{\partial }_{y}\rho -\displaystyle \frac{3}{4}{\rm{\nabla }}\cdot {{\boldsymbol{j}}}_{2\perp }+g\left({h}_{x}\rho -\displaystyle \frac{1}{2}{\boldsymbol{h}}\cdot {{\boldsymbol{j}}}_{2\perp }\right),\end{eqnarray} \tag{ 59c }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{2,x}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-4m\right){j}_{2,x}-{\rm{\nabla }}\cdot {{\boldsymbol{j}}}_{3}+g{\boldsymbol{h}}\ast {\boldsymbol{j}}-g{\boldsymbol{h}}\cdot {{\boldsymbol{j}}}_{3},\end{eqnarray} \tag{ 59d }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{2,y}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-4m\right){j}_{2,y}-{\rm{\nabla }}\cdot {{\boldsymbol{j}}}_{3\perp }+g{\boldsymbol{h}}\ast {{\boldsymbol{j}}}_{\perp }-g{\boldsymbol{h}}\cdot {{\boldsymbol{j}}}_{3\perp },\end{eqnarray} \tag{ 59e }$$

where for compactness we have introduced the notation for the ∗ product of two vectors, defined via ${\boldsymbol{a}}\ast {\boldsymbol{b}}\equiv {a}_{x}{b}_{x}-{a}_{y}{b}_{y}$. Note that equations (59a)–(59e) are similar to those in [21]. One noteworthy difference is the notation: whereas throughout the present work we represent currents as vectors having real components, in, e.g., [21, 46] these quantities are represented via complex numbers. We now translate between these two notations by providing an explicit dictionary. Let us first translate the Fourier components ${{\boldsymbol{j}}}_{n}$ (as defined in equations (58a), (58b)) to complex numbers f_n defined via

$$\begin{eqnarray}&&{f}_{n}({\boldsymbol{r}})=\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }d\theta \,{{\rm{e}}}^{{\rm{i}}{n}\theta }\,P({\boldsymbol{r}},\theta ),\end{eqnarray} \tag{ 60 }$$

where n are integers, positive, negative, or zero. From this definition, we see that ${f}_{0}=\rho$ and ${f}_{n}=(1/2)({j}_{n,x}+{{\rm{i}}{j}}_{n,y})$ for positive n. For negative n's, notice that, ${f}_{-n}=\bar{{f}_{n}}$. We also relate the gradient operator ∇to derivatives with respect to $z\equiv (x+{\rm{i}}{y})/2$ via $({\partial }_{x},{\partial }_{y})\to {\partial }_{z}\equiv {\partial }_{x}-{\rm{i}}{\partial }_{y}$. From this equivalence, ${\rm{\nabla }}\cdot {{\boldsymbol{j}}}_{n}=2{\mathfrak{R}}[{\partial }_{z}{f}_{n}]$, ${\rm{\nabla }}\cdot {{\boldsymbol{j}}}_{n\perp }=2{\mathfrak{I}}[{\partial }_{z}{f}_{n}]$, ${\rm{\nabla }}\ast {{\boldsymbol{j}}}_{n}=2{\mathfrak{R}}[{\partial }_{\bar{z}}{f}_{n}]$, and ${\rm{\nabla }}\ast {{\boldsymbol{j}}}_{n\perp }=2{\mathfrak{I}}[{\partial }_{\bar{z}}{f}_{n}]$. Finally, the product operators between vectors that we denote via · and ∗ are similarly translated: ${{\boldsymbol{j}}}_{m}\cdot {{\boldsymbol{j}}}_{n}=4{\mathfrak{R}}[{f}_{-m}{f}_{n}]$ and ${{\boldsymbol{j}}}_{m}\cdot {{\boldsymbol{j}}}_{n\perp }=4\,{\mathfrak{I}}[{f}_{-m}{f}_{n}];$ for ∗, ${{\boldsymbol{j}}}_{m}\ast {{\boldsymbol{j}}}_{n}=4{\mathfrak{R}}[{f}_{m}{f}_{n}]$ and ${{\boldsymbol{j}}}_{m}\ast {{\boldsymbol{j}}}_{n\perp }=4\,{\mathfrak{I}}[{f}_{m}{f}_{n}]$. These expressions can be substituted into equations (59b)–(59e) to transform these two pairs of real equations into two complex-number equations. In the present work, we continue with the real-vector notation, but the above dictionary can also be used to translate the rest of the equations to the complex convention.

From equations (59a)–(59e) we explicitly see that the interaction terms (which are proportional to g) couple the higher-order Fourier modes to the lowest ones. Nevertheless, notice that the dependence of ${{\boldsymbol{j}}}_{2}$ on ${\boldsymbol{j}}$ is of higher order in the nonlinearity: ${\boldsymbol{h}}\propto {\boldsymbol{j}}$, and the interaction term is quadratic in the currents. Thus, we may deduce whether the isotropic phase is stable by performing a linear stability analysis in which we assume that the current density $| {\boldsymbol{j}}|$ is small compared to the particle density ρ. Then, all higher Fourier modes, such as ${{\boldsymbol{j}}}_{2}$, may be neglected, and we obtain a linearized theory. As this approach neglects all stabilizing nonlinear terms, it does not yield a description of the polar active phase—we leave that task to the following subsections.

We thus proceed with examining the stability of the isotropic phase, while neglecting all nonlinear terms. The linearized equations are

$$\begin{eqnarray}&&{\partial }_{t}\rho =-\displaystyle \frac{1}{2}{\rm{\nabla }}\cdot \left({\boldsymbol{j}}-\displaystyle \frac{1}{\kappa }{\rm{\nabla }}\rho \right),\,\end{eqnarray} \tag{ 61a }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{j}}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-m\right){\boldsymbol{j}}-\displaystyle \frac{1}{2}{\rm{\nabla }}\rho +\displaystyle \frac{{gN}}{2A}{\boldsymbol{j}}.\end{eqnarray} \tag{ 61b }$$

In order to study the stability of the isotropic phase within equations (61a), (61b), we first look for solutions of the form $(\rho ,{\boldsymbol{j}})=({\rho }_{0},{{\boldsymbol{j}}}_{0}){{\rm{e}}}^{\lambda t}$. We take spatial Fourier transforms, which re-express the gradient terms through the wavevector ${\bf{k}}\equiv ({k}_{x},{k}_{y})$. These steps allow us to transform the above differential equations into an eigenvalue problem, wherein ${\rho }_{0}$ and ${{\boldsymbol{j}}}_{0}$ act as eigenvector components and λ as an eigenvalue. The stability of the solutions of this system can then be analyzed by looking at the sign of the eigenvalues for each value of ${\boldsymbol{k}}$. Specifically, there are three eigenvalues associated with the right-hand side of equations (61a), (61b):

$$\begin{eqnarray}&&{\lambda }_{1}=-\displaystyle \frac{{k}^{2}}{2\kappa }-\left(m-\displaystyle \frac{{gN}}{2A}\right),\,\end{eqnarray} \tag{ 62a }$$

$$\begin{eqnarray}&&{\lambda }_{\pm }=-\displaystyle \frac{{k}^{2}}{2\kappa }-\displaystyle \frac{1}{2}\left(m-\displaystyle \frac{{gN}}{2A}\right)\pm \displaystyle \frac{1}{2}{\left[{\left(m-\displaystyle \frac{{gN}}{2A}\right)}^{2}-{k}^{2}\right]}^{1/2},\end{eqnarray} \tag{ 62b }$$

where $k\equiv | {\bf{k}}|$. From equations (62a), (62b), we note that for small wavenumbers (or, equivalently, long wavelengths) the eigenvalues have a real part that is negative only if $m\gt {gN}/2A$. Defining the areal number density of the particle system via ${\rho }_{d}\equiv N/A$, we may rewrite this condition as

$$\begin{eqnarray}&&m\gt g\,\displaystyle \frac{{\rho }_{d}}{2}.\end{eqnarray} \tag{ 63 }$$

Physically, this condition corresponds to the regime in which the disordered (i.e., isotropic or non-polar) fluid phase is stable with respect to small fluctuations. In fact, the same condition appears in the analysis of the XY and Kuramoto models of, respectively, two-dimensional spins and synchronizing oscillators; see, e.g., [47, 48]. In the special case of maximally large orientational noise (i.e., $m\to \kappa /8$), this condition implies that

$$\begin{eqnarray}&&\displaystyle \frac{1}{4g{\rho }_{d}}\gt {\kappa }^{-1}.\end{eqnarray} \tag{ 64 }$$

This counterintuitive result states that the translational noise must be large in order for the disordered phase to be stable. This is a special feature of the model we are considering, in which the coupling between translational and orientational noise (and the resulting 'uncertainty relation') requires a large translational noise if the orientations of the particles are to be distributed narrowly.

In the following subsection, we examine the behavior of long-wavelength fluctuations in the vicinity of the transition between the polar and disordered phases of the active fluid. Let us first consider the case in which the disordered phase is highly unstable, and show that the system's preferred spatially homogeneous state has rotational-symmetry-broken polar order. As a polar fluid maintains translational symmetry, all spatial derivatives in the Fokker–Planck equation (56) vanish, and the equation reduces to

$$\begin{eqnarray}&&{\partial }_{t}P=g{\partial }_{\theta }\,[({{\boldsymbol{h}}}_{\perp }\cdot {\boldsymbol{v}})P]+m\,{\partial }_{\theta }^{2}P.\end{eqnarray} \tag{ 65 }$$

In this homogeneous case, pressure and translational noise are neglected, and the only obstacle to polar order is the orientational noise. The exact solution of this equation takes the form $P=\tfrac{1}{{AZ}}\exp \left(-\tfrac{g}{m}{\boldsymbol{h}}\cdot {\boldsymbol{v}}\right)$, in which Z is a normalization constant. Although this is an exact solution, we obtain a simpler expression by only considering the case in which $| {\boldsymbol{j}}|$ is small. It is more illustrative, though, to operate in terms of a Landau theory: in this regime, we may consider only up to the first three Fourier modes, which brings us to the following systems of equations:

$$\begin{eqnarray}{\partial }_{t}\rho & = & 0,\\ {\partial }_{t}{j}_{x} & = & -m\,{j}_{x}+g\left({h}_{x}\,\rho -\displaystyle \frac{1}{2}{\boldsymbol{h}}\cdot {{\boldsymbol{j}}}_{2}\right),\end{eqnarray} \tag{ 66a }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{y}=-m\,{j}_{y}+g\left({h}_{x}\,\rho -\displaystyle \frac{1}{2}{\boldsymbol{h}}\cdot {{\boldsymbol{j}}}_{2\perp }\right),\,\end{eqnarray} \tag{ 66b }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{2,x}=-4m\,{j}_{2,x}+g\,{\boldsymbol{h}}\ast {\boldsymbol{j}},\,\end{eqnarray} \tag{ 66c }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{2,y}=-4m\,{j}_{2,y}+g\,{\boldsymbol{h}}\ast {{\boldsymbol{j}}}_{\perp }.\,\end{eqnarray} \tag{ 66d }$$

Solving these equations, and recalling that ${\boldsymbol{h}}=\pi N{\boldsymbol{j}}$, leads us to the following self-consistent equation for the alignment vector ${\boldsymbol{h}}$:

$$\begin{eqnarray}&&{\boldsymbol{h}}=\displaystyle \frac{{\rho }_{d}}{2}\left[\displaystyle \frac{g}{m}-\displaystyle \frac{1}{8}{\left(\displaystyle \frac{g}{m}\right)}^{3}{h}^{2}\right]{\boldsymbol{h}},\end{eqnarray} \tag{ 67 }$$

where $h=| {\boldsymbol{h}}|$. Solutions of equation (67) are either the trivial ${\boldsymbol{h}}=0$ solution (which is unstable for $g{\rho }_{d}\lt 2m$), or the solution

$$\begin{eqnarray}&&h=4{\left(\displaystyle \frac{{m}^{3}}{{\rho }_{d}{g}^{3}}\right)}^{1/2}\sqrt{\displaystyle \frac{{\rho }_{d}g}{2m}-1}.\end{eqnarray} \tag{ 68 }$$

This solution is only physical for $g{\rho }_{d}\geqslant 2m$. Thus, this inequality presents the boundary between the disordered fluid phase and the ordered, polar, phase of the active fluid. We show this boundary in figure 2. Note that although in this subsection we have neglected all spatial derivatives, it is crucial that we have included the contributions of the higher Fourier modes in the orientational angle, i.e., ${{\boldsymbol{j}}}_{2}$. The coupling of this mode to the lower Fourier modes introduces nonlinearity in the system, which is necessary to stabilize the ordered polar phase.

Zoom In Zoom Out Reset image size

The self-consistent approximation is known not to hold for phase transitions in which fluctuations are strong. Fluctuations can shift the transition line in parameter space and, furthermore, can change the nature of the transition. We thus expect the prediction that the order-disorder transition in polar active fluids is continuous not to hold—in physical polar active fluids, the transition may be discontinuous.

Let us now go beyond the assumption of spatial homogeneity and examine the effects of temporally slow, long-wavelength fluctuations in the density and the current for the case of a polar active fluid. This yields the hydrodynamic theory which, as for any polar active fluid, reduces to a form of Toner–Tu theory [31, 49]. We proceed via the method introduced in [21]. Starting from equations (59a)–(59e), we go to the regime in which the order is weak and therefore $P({\boldsymbol{r}},\theta )$ depends only weakly on θ. Physically, this is the regime in which the active flow is much slower than the microscopic speed of each particle. Mathematically, this regime allows us to discard all Fourier modes higher than ${{\boldsymbol{j}}}_{2}$. (NB: for the sake of generality, let us also include the parameter λ from equation (37) which controls the feedback between the translational noise and the angle.) Then, we consider the equations

$$\begin{eqnarray}&&{\partial }_{t}\rho =-\displaystyle \frac{1}{2}{\rm{\nabla }}\cdot \left({\boldsymbol{j}}-\displaystyle \frac{1}{\kappa }{\rm{\nabla }}\rho \right),\end{eqnarray} \tag{ 69a }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{x}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-m\right){j}_{x}+\left(\displaystyle \frac{\lambda -2}{2}\right){\partial }_{x}\rho ,\end{eqnarray} \tag{ 69b }$$

$$\begin{eqnarray}&&-\left(\displaystyle \frac{\lambda +2}{4}\right){\rm{\nabla }}\cdot {{\boldsymbol{j}}}_{2}+g\left({h}_{x}\rho -\displaystyle \frac{1}{2}{\boldsymbol{h}}\cdot {{\boldsymbol{j}}}_{2}\right),\end{eqnarray} \tag{ 69c }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{y}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-m\right){j}_{y}+\left(\displaystyle \frac{\lambda -2}{2}\right){\partial }_{y}\rho ,\end{eqnarray} \tag{ 69d }$$

$$\begin{eqnarray}&&-\left(\displaystyle \frac{\lambda +2}{4}\right){\rm{\nabla }}\cdot {{\boldsymbol{j}}}_{2\perp }+g\left({h}_{x}\rho -\displaystyle \frac{1}{2}{\boldsymbol{h}}\cdot {{\boldsymbol{j}}}_{2\perp }\right),\end{eqnarray} \tag{ 69e }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{2,x}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-4m\right){j}_{2,x}+\left(\displaystyle \frac{\lambda -1}{2}\right){\rm{\nabla }}\ast {\boldsymbol{j}}+g{\boldsymbol{h}}\ast {\boldsymbol{j}},\end{eqnarray} \tag{ 69f }$$

$$\begin{eqnarray}&&{\partial }_{t}{j}_{2,y}=\left(\displaystyle \frac{1}{2\kappa }{{\rm{\nabla }}}^{2}-4m\right){j}_{2,y}+\left(\displaystyle \frac{\lambda -1}{2}\right){\rm{\nabla }}\ast {{\boldsymbol{j}}}_{\perp }+g{\boldsymbol{h}}\ast {{\boldsymbol{j}}}_{\perp }.\end{eqnarray} \tag{ 69g }$$

Note two important aspects of the associated approximation: (i) $| {{\boldsymbol{j}}}_{2}|$ is much smaller than $| {\boldsymbol{j}}| ;$ and (ii) as we are interested in hydrodynamics, we only consider time- and length-scales much larger than the microscopic ones. Therefore, we consider the regime characterized by

$$\begin{eqnarray}&&{\partial }_{t}{{\boldsymbol{j}}}_{2,x},{\partial }_{t}{{\boldsymbol{j}}}_{2,y}\ll m{{\boldsymbol{j}}}_{2,x},m{{\boldsymbol{j}}}_{2,y}\end{eqnarray} \tag{ 70 }$$

in which both the time-derivatives ${\partial }_{t}{{\boldsymbol{j}}}_{2}$ and the term $\tfrac{1}{\kappa }{{\rm{\nabla }}}^{2}{{\boldsymbol{j}}}_{2}$ are neglected. Re-expressing the current ${\boldsymbol{j}}$ via ${\boldsymbol{j}}={\boldsymbol{h}}/\pi N$ and ${\rho }_{d}=2\pi N\rho$, we can rewrite the above equations as

$$\begin{eqnarray}&&{\partial }_{t}{\rho }_{d}=-{\rm{\nabla }}\cdot \left({\boldsymbol{h}}-\displaystyle \frac{1}{2\kappa }{\rm{\nabla }}{\rho }_{d}\right)\end{eqnarray} \tag{ 71 }$$

and

$$\begin{eqnarray}{\partial }_{t}{\boldsymbol{h}} & = & \left[\left(\displaystyle \frac{{\rho }_{d}g}{2}-m\right)-\displaystyle \frac{{g}^{2}}{8m}{h}^{2}\right]{\boldsymbol{h}}+\left(\displaystyle \frac{\lambda -2}{4}\right){\rm{\nabla }}{\rho }_{d}\\ & & -(\lambda -1)\displaystyle \frac{g}{16m}[({\boldsymbol{h}}\cdot {\rm{\nabla }}){\boldsymbol{h}}+({{\boldsymbol{h}}}_{\perp }\cdot {\rm{\nabla }}){{\boldsymbol{h}}}_{\perp }]\\ & & -(\lambda +2)\displaystyle \frac{g}{8m}[{\boldsymbol{h}}({\rm{\nabla }}\cdot {\boldsymbol{h}})-{{\boldsymbol{h}}}_{\perp }({\rm{\nabla }}\cdot {{\boldsymbol{h}}}_{\perp })]\\ & & +\left[\displaystyle \frac{1}{2\kappa }+\displaystyle \frac{2-\lambda -{\lambda }^{2}}{32m}\right]{{\rm{\nabla }}}^{2}{\boldsymbol{h}}.\end{eqnarray} \tag{ 72 }$$

Finally, the unusual ${{\boldsymbol{h}}}_{\perp }$ terms may be rewritten in a more familiar way via the identities

$$\begin{eqnarray}{{\boldsymbol{h}}}_{\perp }({\rm{\nabla }}\cdot {{\boldsymbol{h}}}_{\perp }) & = & \displaystyle \frac{1}{2}{\rm{\nabla }}{h}^{2}-({\boldsymbol{h}}\cdot {\rm{\nabla }}){\boldsymbol{h}},\\ ({{\boldsymbol{h}}}_{\perp }\cdot {\rm{\nabla }}){{\boldsymbol{h}}}_{\perp } & = & \displaystyle \frac{1}{2}{\rm{\nabla }}{h}^{2}-{\boldsymbol{h}}({\rm{\nabla }}\cdot {\boldsymbol{h}}),\end{eqnarray} \tag{ 73 }$$

where, as above, $h\equiv | {\boldsymbol{h}}|$. The resulting hydrodynamic equations read:

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{h}}+(\lambda +1)\displaystyle \frac{3g}{16m}({\boldsymbol{h}}\cdot {\rm{\nabla }}){\boldsymbol{h}}=\left[\left(\displaystyle \frac{{\rho }_{d}g}{2}-m\right)-\displaystyle \frac{{g}^{2}}{8m}{h}^{2}\right]{\boldsymbol{h}}\\ &&\quad +\,\left(\displaystyle \frac{\lambda -2}{4}\right){\rm{\nabla }}{\rho }_{d}+(\lambda +5)\displaystyle \frac{g}{16m}\left[\displaystyle \frac{1}{2}{\rm{\nabla }}{h}^{2}-{\boldsymbol{h}}({\rm{\nabla }}\cdot {\boldsymbol{h}})\right]\\ &&\quad +\,\left[\displaystyle \frac{1}{2\kappa }+\displaystyle \frac{2-\lambda -{\lambda }^{2}}{32m}\right]{{\rm{\nabla }}}^{2}{\boldsymbol{h}}.\end{eqnarray} \tag{ 74 }$$

Equation (71) is the continuity equation that established the conservation of the number of self-propelled particles. On the other hand, equation (74) describes the hydrodynamics of the polarization order parameter for the active fluid. Note that the order parameter ${\boldsymbol{h}}$ is not identical to the current of the active fluid, which in addition includes the part of the motion due to translational diffusion.

Let us now discuss the result for the hydrodynamics, equation (74). The second term on the left-hand side is an advection term. In general, as consequence of the breaking of Galilean invariance in our system, its prefactor is not unity. This results from the presence of a special reference frame in the system: there is only one frame in which the particles can propagate with equal speed, regardless of their direction of movement. On the right-hand side, the first term is the one responsible for the spontaneous breaking of the rotational symmetry of the fluid in the polar state. The second and third terms are pressure-like terms. The third term includes the effects of a nonlinear compressibility.

The fourth term is a viscous damping term familiar from, e.g., the Navier–Stokes equations. In the present context, this term includes the effects of both the translational noise and the coupling between the translational and orientational noises. Curiously, the parameter λ has a strong effect on the effective viscosity. E.g., if $\lambda =1$ (the model we focus on for the single-particle dynamics), this viscosity reduces to $1/2\kappa$. In that case, the viscosity depends only on the translational noise. At first glance, the expression for the effective viscosity does not appear to be positive-definite. If the viscosity were to change sign, this would suggest that the homogeneous fluid state should become unstable. However, in the previous section we showed that the particle dynamics is physical only if the inequality $m\gt \kappa {\lambda }^{2}/8$ holds. We now re-express the mass in terms of a positive parameter , via $m=\kappa [\epsilon +({\lambda }^{2}/8)]$, and write the effective viscosity as

$$\begin{eqnarray}&&{\nu }_{{\rm{eff}}}=\displaystyle \frac{1}{2\kappa }\left[1+\displaystyle \frac{2-\lambda -{\lambda }^{2}}{16[\epsilon +({\lambda }^{2}/8)]}\right].\end{eqnarray} \tag{ 75 }$$

The effective viscosity as given by equation (75) is positive for λ real and nonnegative. To see this, note that when $2-\lambda -{\lambda }^{2}\lt 0$, equation (75) yields the minimal possible value for the effective viscosity for each value of λ in the limit $\epsilon \to 0$. This value is always positive and converges to $1/4\kappa$ as $\lambda \to \pm \infty$. On the other hand, when $2-\lambda -{\lambda }^{2}\gt 0$, the minimum value that the effective viscosity can achieve is $1/2\kappa$. We plot the minimum value of the effective viscosity for each value of λ in figure 3. Curiously, it is also worth noting that in the case $\lambda =-1$ (a noise coupling with the same strength but with the opposite sign of our original model), the advection term vanishes.

Zoom In Zoom Out Reset image size

The hydrodynamic equations (74), along with the continuity equation (71), are a version of the Toner–Tu equations, first derived in [31] based on symmetry considerations. By contrast, we obtain these equations based on the microscopic single-particle model that we introduced in the previous section along with inter-particle interactions. The coefficients that we obtained thus explicitly depend on the microscopic parameters of the model, which allows not only for the form but also for the precise numerical evaluation for the hydrodynamic coefficients in the Toner–Tu equations. In the special limit $(\kappa ,\lambda )\to (\infty ,0)$, with ${\lambda }^{2}\kappa \to 0$, the hydrodynamic equations that we derive coincide exactly with the ones derived in [49].