Critical phenomenon of the order–disorder transition in incompressible active fluids

Emergent properties of interacting non-equilibrium systems are of widespread and fundamental interest. One of the simplest, but most striking, of these is the self-organized phenomenon of 'flocking'—that is, collective motion (CM) in large groups of motile organisms [1–10]. This phemenon is fascinating in part because its occurrence in two spatial dimensions requires the spontaneous breaking of a continuous symmetry, which is forbidden in thermal equilibrium by the Mermin–Wagner theorem [11]. It was initially hoped [2] that the transition into this novel state could be a continuous one belonging to a new universality class. However, it was subsequently realized, from both simulations and theoretical analysis [12–18] of the hydrodynamic equations [3–6], that as this putative continuous transition is approached from the ordered side, but before it can be reached, the homogeneous CM state becomes unstable to modulation of the density along the mean velocity. The transition from the homogeneous CM state (i.e., the ordered state) to the disordered state proceeds via two first order transitions: one from homogeneous to banded, the next from banded to disordered.

Since this instability requires density variations, we reason that the instability might be eliminated by making the system incompressible. In this paper, we show that, indeed, the order-disorder transition is continuous in an incompressible system, and belongs to a new universality class. We emphasize that the incompressibility is crucial for the existence of this continuous transition; in a compressible system, as shown in the theoretical analyses listed above [12–18], the continuous transition is always preempted by the aforementioned banding instability. In an incompressible system, such bands cannot form, since the density cannot vary. We demonstrate this by finding, in a dynamical renormalization group (DRG) analysis of the hydrodynamic equations for an incompressible active fluid in d spatial dimensions, a novel stable fixed point that controls the transition. This calculation is done to order $\mathcal{O}(\epsilon )$ in an $\epsilon =4-d$ expansion; to the same order, we calculate the critical exponents of the transition. We also obtain two scaling laws relating these critical exponents which are valid to all orders in (i.e., exact).

Our results are testable in both experiments and simulations. Three potential realizations are:

• (i)  
  Systems with strong repulsive short-ranged interactions between the active particles. Incompressibility has, in fact, been assumed in, e.g., recent experimental studies on cell motility [19]. In such systems, the compressibility will be non-zero, but small. Hence, our incompressible results will apply out to very large length scales, or, equivalently, very close to the transition, but will ultimately crossover to the compressible behaviour (i.e., a small first order transition driven by the banding instability).
• (ii)  
  Systems with long-ranged repulsive interactions; here, true incompressibility is possible. Long-ranged interactions are quite reasonable in certain contexts: birds, for example, can often see all the way across a flock [20].
• (iii)  
  Motile colloidal systems in fluid-filled microfluidic channels. The forces exerted by the active particles are, of course, tiny compared to what would be needed to compress the background fluid, so that fluid is effectively incompressible. Since the active particles drag the background fluid with them, their motion is effectively incompressible as well. Indeed, experiments [21] show these systems do not exhibit the banding instability [12–18] found in all compressible active systems. This also suggests a numerical approach: simulating active particles moving through an incompressible fluid [22].

We formulate the most general hydrodynamic model for systems lacking both momentum conservation, and Galilean invariance, consistent with the symmetries of rotation and translation invariance, and the assumption of incompressibility. As the number density cannot fluctuate (by the assumption of incompressibility), the velocity field is the only hydrodynamic variable in the problem, which becomes soft as the transition is approached. Since the velocity is small near the transition, we can expand the equation of motion (EOM) in powers of the velocity. The symmetry constraints of translation and rotation invariance force the EOM valid at long wavelengths and times to take the form: [3–6]⁵

$$\begin{eqnarray}&&{{\partial }_{t}}{\bf v}+\lambda ({\bf v}\cdot \nabla ){\bf v}=-\nabla \mathcal{P}-(a+b|{\bf v}{{|}^{2}}){\bf v}+\mu {{\nabla }^{2}}{\bf v}+{\bf f},\end{eqnarray} \tag{ 1 }$$

where the pressure $\mathcal{P}$ enforces the incompressibility condition $\nabla \cdot {\bf v}=0$, ${\bf f}$ is a Gaussian 'white noise' with spatio-temporally Fourier transformed statistics:

$$\begin{eqnarray}&&\langle {{f}_{m}}({\bf k},\omega ){{f}_{n}}({\bf k}^{\prime} ,\omega ^{\prime} )\rangle =2D{{P}_{mn}}({\bf k})\delta ({\bf k}+{\bf k}^{\prime} )\delta (\omega +\omega ^{\prime} ),\end{eqnarray} \tag{ 2 }$$

and ${{P}_{mn}}({\bf k})\equiv {{\delta }_{mn}}-{{k}_{m}}{{k}_{n}}/{{k}^{2}}$ is the transverse projection operator. The compressible model [3–6] differs from this model only in the pressure $\mathcal{P}$ term, which is in this case a specified function of the density, and in that $\nabla \cdot {\bf v}\ne 0$ which allows a term proportional to ($\nabla$ · v)v as well. One need also then give an EOM for the density, which is just the usual continuity equation .

The EOM equation (1) reduces, when $a=0=b$, to the classic model of a fluid forced at zero wavenumber treated by [23] (their 'model B'). For $\lambda =0$, it reduces to a simple, time-dependent Ginzburg–Landau [24, 25] dynamical model for an isotropic ferromagnet with long ranged dipolar interactions [26] ⁶ .

Because our system lacks Galilean invariance, $\lambda$ need not (and in general will not) equal one, and the terms $-(a+b|{\bf v}{{|}^{2}}){\bf v}$ are allowed in the EOM. The latter terms are crucial as they make possible the existance of a polar ordered phase in an active system, which is not possible in a normal fluid.

At the mean field level, for $a\lt 0,b\gt 0$ the system is in the ordered phase⁷ with $|{\bf v}|=\sqrt{-a/b}$, and for $a\gt 0,b\gt 0$ it is in the disordered phase with $|{\bf v}|=0$. To go beyond this mean field description, we employ the DRG method [23] near the order-disorder transition. To do so, we spatio-temporally Fourier transform equation (1), and project orthogonal to wavevector ${\bf k}$; obtaining

$$\begin{eqnarray}\begin{array}{rcl} {{v}_{l}}(\tilde{{\bf k}}) & = & G(\tilde{{\bf k}})\left[ {{f}_{l}}(\tilde{{\bf k}})-\frac{{\rm i}\lambda }{2}{{P}_{lmn}}({\bf k}){{\int }_{\tilde{{\bf q}}}}{{v}_{m}}(\tilde{{\bf q}}){{v}_{n}}(\tilde{{\bf k}}-\tilde{{\bf q}}) \right.\left. -\frac{b}{3}{{Q}_{lmnp}}({\bf k}){{\int }_{\tilde{{\bf q}},\tilde{{\bf h}}}}{{v}_{m}}(\tilde{{\bf k}}-\tilde{{\bf q}}-\tilde{{\bf h}}){{v}_{n}}(\tilde{{\bf q}}){{v}_{p}}(\tilde{{\bf h}}) \right], \\ \end{array}\end{eqnarray} \tag{ 3 }$$

where we have adopted the reduced notations $\tilde{{\bf k}}\equiv ({\bf k},\omega )$ and ${{\int }_{\tilde{{\bf q}}}}={{\int }_{{\bf q},\Omega }}\equiv \int \frac{{{{\rm d}}^{d}}q}{{{(2\pi )}^{d}}}\frac{{\rm d}\Omega }{2\pi }$, and we have defined ${{P}_{lmn}}({\bf k})\equiv {{P}_{\operatorname{lm}}}({\bf k}){{k}_{n}}+{{P}_{{\rm ln} }}({\bf k}){{k}_{m}}$, ${{Q}_{lmnp}}({\bf k})\equiv {{P}_{\operatorname{lm}}}({\bf k}){{\delta }_{np}}+{{P}_{{\rm ln} }}({\bf k}){{\delta }_{mp}}+{{P}_{lp}}({\bf k}){{\delta }_{mn}}$, and the 'propagator' $G(\tilde{{\bf k}})\equiv {{(-{\rm i}\omega +\mu {{k}^{2}}+a)}^{-1}}$. This propagator is just the linearized response function giving the velocity induced, to linear order, by an external force.

Graphical representations of the various terms in equation (3) are shown in figure 1. For an excellent and detailed exposition on the meaning and use of graphical techniques in the dynamical perturbative renormalization group (DRG) approach, we refer the reader to [23].

Zoom In Zoom Out Reset image size

We now perform the standard DRG procedure [23], averaging over short wavelength degrees of freedom, and rescaling: ${\bf r}\to {{e}^{\ell }}{\bf r}$, $t\to {{e}^{z\ell }}t$ and ${\bf v}\to {{e}^{\chi \ell }}{\bf v}$. Here the dynamical exponent z and the roughness exponent $\chi$ respectively relate time and velocity scales in the renormalized system to those in the physical problem. The rescaling exponents z and $\chi$ are, of course, completely arbitrary. A very convenient choice of them, however, is one that keeps the scale of velocity fluctuations in the renormalized system fixed. The values of z and $\chi$ that we obtain by this criterion will, by construction, then give the scaling of the physical time scales and velocity fluctuations with length scale (hence the terms 'dynamical exponent' and 'roughness exponent').

Our procedure is identical to the calculation for model B in [23], except for a modified propagator, and some additional Feynmann graphs due to the extra $b|{\bf v}{{|}^{2}}{\bf v}$ nonlinearity in our problem. At the one loop level, the non-vanishing graphical contributions to the various coefficients in equation (3) are shown in figure 2. More details of the calculation are given in the supplemental materials (stacks.iop.org/NJP/17/042002/mmedia). We obtain, in d spatial dimensions, the following RG flow equations of the coefficients to one loop order and to linear order in $\epsilon \equiv 4-d$:

$$\begin{eqnarray}&&\frac{{\rm d}a}{{\rm d}\ell }=za-\frac{9{{g}_{2}}}{2}(a-\mu {{\Lambda }^{2}}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\frac{{\rm d}b}{{\rm d}\ell }=\left( 2\chi +z-\frac{17{{g}_{2}}}{2} \right)b,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\frac{{\rm d}\lambda }{{\rm d}\ell }=\left( \chi -1+z-\frac{5{{g}_{2}}}{3} \right)\lambda ,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\frac{{\rm d}\mu }{{\rm d}\ell }=\left( -2+z+\frac{{{g}_{1}}}{4} \right)\mu ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\frac{{\rm d}D}{{\rm d}\ell }=\left( -2\chi +z-d \right)D,\end{eqnarray} \tag{ 8 }$$

where we have defined dimensionless couplings:

$$\begin{eqnarray}&&{{g}_{1}}\equiv \frac{{{S}_{d}}D{{\lambda }^{2}}}{{{\left( 2\pi \right)}^{d}}{{\mu }^{3}}}{{\Lambda }^{-\epsilon }},\;{{g}_{2}}\equiv \frac{{{S}_{d}}Db}{{{\left( 2\pi \right)}^{d}}{{\mu }^{2}}}{{\Lambda }^{-\epsilon }},\end{eqnarray} \tag{ 9 }$$

and where ${{S}_{d}}\equiv 2{{\pi }^{d/2}}/\Gamma (d/2)$ is the surface area of a unit sphere in d dimensions, $\epsilon \equiv 4-d$, and Λ is the ultraviolet wavevector cutoff which is of order the inverse of some 'microscopic' length (e.g., the typical distance between active particles). Since our interest is in the transition, we have, in the last four recursion relations (5)–(8), set a = 0, and have worked to linear order in a in (4). It is straightforward to verify that higher order terms in a affect none of our results up to and including linear order in $\epsilon =4-d$.

Zoom In Zoom Out Reset image size

Now we derive coupled, closed recursion relations for ${{g}_{1,2}}$ from the recursion relations (4–8). We begin by taking the natural logarithm of our definitions equation (9) of ${{g}_{1,2}}$:

$$\begin{eqnarray}&&{\rm ln} \left( {{g}_{1}} \right)={\rm ln} \left( D \right)+2{\rm ln} \left( \lambda \right)-3{\rm ln} \left( \mu \right)+{\rm ln} \left( C \right),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rm ln} \left( {{g}_{2}} \right)={\rm ln} \left( D \right)+{\rm ln} \left( b \right)-2{\rm ln} \left( \mu \right)+{\rm ln} \left( C \right),\end{eqnarray} \tag{ 11 }$$

where we have defined the constant $C\equiv \frac{{{S}_{d}}}{{{\left( 2\pi \right)}^{d}}}{{\Lambda }^{-\epsilon }}$, which is 'constant' in the sense that it does not change upon renormalization. Differentiating both sides of (10) and (11) with respect to RG time ℓ therefore gives

$$\begin{eqnarray}&&\frac{1}{{{g}_{1}}}\frac{{\rm d}{{g}_{1}}}{{\rm d}l}=\frac{1}{D}\frac{{\rm d}D}{{\rm d}\ell }+\frac{2}{\lambda }\frac{{\rm d}\lambda }{{\rm d}\ell }-\frac{3}{\mu }\frac{{\rm d}\mu }{{\rm d}\ell },\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\frac{1}{{{g}_{2}}}\frac{{\rm d}{{g}_{2}}}{{\rm d}l}=\frac{1}{D}\frac{{\rm d}D}{{\rm d}\ell }+\frac{1}{b}\frac{{\rm d}b}{{\rm d}\ell }-\frac{2}{\mu }\frac{{\rm d}\mu }{{\rm d}\ell }.\end{eqnarray} \tag{ 13 }$$

Plugging the recursion relations (5)–(8) on the right-hand side of (12) and (13), gathering terms, and multiplying (12) by g₁ and (13) by g₂ gives two closed flow equations for ${{g}_{1,2}}$ for arbitrary $\chi$ and z:

$$\begin{eqnarray}&&\frac{{\rm d}{{g}_{1}}}{{\rm d}\ell }=\epsilon {{g}_{1}}-\frac{3}{4}g_{1}^{2}-\frac{10}{3}{{g}_{1}}{{g}_{2}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\frac{{\rm d}{{g}_{2}}}{{\rm d}\ell }=\epsilon {{g}_{2}}-\frac{1}{2}{{g}_{1}}{{g}_{2}}-\frac{17}{2}g_{2}^{2}.\end{eqnarray} \tag{ 15 }$$

Although not necessary, it is convenient to make a special choice of z and $\chi$ such that $\mu$ and D are kept fixed at their bare values (i.e, ${{\mu }_{0}}$ and D₀, respectively). The advantage of the (completely arbitrary) choice to keep $\mu$ fixed is that the correlation length $\xi$ is determined by the ratio $\frac{\mu }{a}$ (specifically, $\xi =\sqrt{\frac{\mu }{a}}$). Hence, by keeping $\mu$ fixed, we can determine the scaling of $\xi$ directly from the scaling of a, without having to worry about $\mu$ (since it is held fixed) at all. Likewise, keeping D fixed in addition keeps the fluctuations at the critical point unchanged (since those are controlled by the ratio $\frac{\mu }{D}$, which will, with this choice, be unchanged). This enables us to get the scaling of fluctuations at the critical point directly from the field rescaling $\chi$, which connects physical fluctuations to those in the renormalized system.

We will hereafter adopt this choice of z and $\chi$, which is

$$\begin{eqnarray}&&z=2-{{g}_{1}}/4+\mathcal{O}({{\epsilon }^{2}}),\quad \chi =\frac{z-d}{2}+\mathcal{O}({{\epsilon }^{2}}).\end{eqnarray} \tag{ 16 }$$

We will also hereafter use the subscript 0 to denote the bare (i.e., unrenormalized) values of the parameters.

Equation (4) now becomes

$$\begin{eqnarray}&&\frac{{\rm d}a}{{\rm d}\ell }=\left( 2-\frac{{{g}_{1}}}{4}-\frac{9}{2}{{g}_{2}} \right)a+\frac{9}{2}{{g}_{2}}\mu {{\Lambda }^{2}}.\end{eqnarray} \tag{ 17 }$$

Equations (14), (15), (17) have a non-Gaussian fixed point in $d\lt 4$:

$$\begin{eqnarray}\begin{array}{rcl} g_{1}^{*} & = & \frac{124}{113}\epsilon +\mathcal{O}({{\epsilon }^{2}}) \\ g_{2}^{*} & = & \frac{6}{113}\epsilon +\mathcal{O}({{\epsilon }^{2}}) \\ {{a}^{*}} & = & \left[ -\frac{27}{226}\epsilon +\mathcal{O}({{\epsilon }^{2}}) \right]\mu {{\Lambda }^{2}}, \\ \end{array}\end{eqnarray} \tag{ 18 }$$

which can be shown by analyzing the three recursion relations to be a stable attractor of all points on a two-dimensional surface (the 'critical surface' ) in the three-dimensional parameter space $({{g}_{1}},{{g}_{2}},a)$, but to be unstable with respect to displacements off this critical surface. The flows on the critical surface are illustrated in figure 3. This is exactly the topology of renormalization group flows that corresponds to a continuous phase transition with universal exponents controlled by the fixed point that's stable within the critical surface. Hence, we conclude that the order-disorder is generically continuous in incompressible active fluids. Furthermore, as far as we know, the fixed point we have obtained is novel, and thus the critical behaviour of incompressible active fluids belongs to a new universality class.

Zoom In Zoom Out Reset image size

As for equilibrium phase transitions [24, 25], there are a number of universal critical exponents associated with the order-disorder transition that we can obtain from this RG analysis. Readers interested in a more thorough discussion of these exponents, and their derivation from the renormalization group recursion relations, are referred to any one of a number of excellent textbooks (e.g., [24, 25]) that discuss this topic in far more detail than space allows us here. Suffice it to say that all of the analysis we present here is entirely standard in the study of critical phenomena.

All of these universal critical exponents are related to experimentally measurable correlation and response functions of the active fluid. For example, several of them can be obtained from the two point velocity correlation function $\langle {\bf v}({\bf r}+{\bf R},t+T)\cdot {\bf v}({\bf R},T)\rangle \equiv C({\bf r},t)$, which depends on bare parameters b₀, ${{\mu }_{0}}$, D₀, and ${{\lambda }_{0}}$ and, most importantly, the proximity to the phase transition $\delta {{a}_{0}}\equiv {{a}_{0}}-a_{0}^{c}$, where a₀^c is the value of a₀ at the transition. We will show in a moment that the RG predicts that this correlation function has a scaling form for large r near the transition:

$$\begin{eqnarray}&&C({\bf r},t)={{r}^{2-d-\eta }}{{Y}_{\pm }}\left( \frac{r}{\xi },\frac{t}{{{r}^{z}}} \right),\end{eqnarray} \tag{ 19 }$$

where the scaling functions ${{Y}_{\pm }}$ in equation (19) are different on the disordered (+) and ordered (-) sides of the transition, because the system is in different phases in the two cases. On the disordered side , we expect ${{Y}_{+}}(x,y)$ to decay exponentially with both x and y, while on the ordered side, ${{Y}_{-}}(x,y)$ has a more complicated scaling behaviour that we will discuss elsewhere [32]. Two of the universal critical exponents, $\eta$ and z, are displayed explicitly in (19). A third, the so called 'correlation length exponent' $\nu$, is implicit in the correlation length $\xi$ appearing in (19), which diverges as the transition is approached. In an experiment or simulation, one can approach the transition by tuning any one of many microscopic control parameters s (e.g., $s\;=$ mean density or noise strength). The RG then predicts that the correlation length $\xi$ diverges algebraically with this control parameter as the transition is approached:

$$\begin{eqnarray}&&\xi \propto |s-{{s}_{{\rm c}}}{{|}^{-\nu }},\end{eqnarray} \tag{ 20 }$$

where $\nu$ is a universal exponent and s_c is the value of the control parameter at which the transition occurs.

Another standard universal critical exponent is the 'order parameter exponent' $\beta$, which characterizes the growth of the order parameter (in our case, the mean velocity $\langle {\bf v}\rangle$ ):

$$\begin{eqnarray}&&|\langle {\bf v}\rangle |\propto |{{s}_{{\rm c}}}-s{{|}^{\beta }}.\end{eqnarray} \tag{ 21 }$$

The final two exponents characterize the response, near or at the transition, of the system to a weak external field ${\bf H}$; that is, simply adding a small constant vector ${\bf H}$ to the rhs of (1). On the disordered side of the transition, we expect linear response:

$$\begin{eqnarray}&&\langle {\bf v}\rangle ={{\chi }_{H}}(s){\bf H},\end{eqnarray} \tag{ 22 }$$

with the susceptibility ${{\chi }_{H}}$ (which should not be confused with the 'roughness' exponent $\chi$ introduced earlier) diverging with a universal exponent $\gamma$ as the transition is approached from the disordered side:

$$\begin{eqnarray}&&{{\chi }_{_{H}}}(s)\propto |s-{{s}_{{\rm c}}}{{|}^{-\gamma }}.\end{eqnarray} \tag{ 23 }$$

One can also ask about the response to this field right at the transition $s={{s}_{{\rm c}}}$. This defines the final universal exponent δ, via:

$$\begin{eqnarray}&&|\langle {\bf v}\rangle |\propto |{\bf H}{{|}^{\frac{1}{\delta }}},\end{eqnarray} \tag{ 24 }$$

where $\langle {\bf v}\rangle$ is in the same direction as the external field ${\bf H}$.

All of the scaling laws we have just specified are, of course, the leading order behaviour near the transition. Corrections to this themselves obey a scaling law, characterized by a 'correction to scaling exponent' ${{y}_{2}}\lt 0$; for example

$$\begin{eqnarray}&&{{\chi }_{_{H}}}(s)=A|s-{{s}_{{\rm c}}}{{|}^{-\gamma }}\left[ 1+{{\left( \xi \Lambda \right)}^{{{y}_{2}}}} \right],\end{eqnarray} \tag{ 25 }$$

where A is a non-universal constant.

We will now use the RG to derive all of these scaling laws, and calculate, to linear order in $\epsilon =4-d$, the values of the universal exponents.

We begin by noting that the exponential runaway from the critical surface in the unstable direction near the stable fixed point (18) grows like ${{{\rm e}}^{{{y}_{a}}\ell }}$, with the exponent

$$\begin{eqnarray}&&{{y}_{a}}=2-\frac{58}{113}\epsilon +\mathcal{O}({{\epsilon }^{2}}).\end{eqnarray} \tag{ 26 }$$

This eigenvalue determines the critical exponent $\nu$ governing the critical behaviour of the velocity correlation length $\xi$. In addition, the smaller (in magnitude) of the two negative eigenvalues of the critical fixed point gives the 'correction to scaling exponent' y₂ [24]; we find ${{y}_{2}}=-\frac{31}{113}\epsilon +\mathcal{O}({{\epsilon }^{2}}).$

To show how y_a determines $\nu$, and to derive the scaling law (19), we use the RG to connect the original $C({\bf r},t)$ to that of the rescaled system:

$$\begin{eqnarray}&&C({\bf r},t;\delta {{a}_{0}},{{b}_{0}},{{\lambda }_{0}})={{{\rm e}}^{2\chi \ell }}C\left( {\bf r}{{{\rm e}}^{-\ell }},t{{{\rm e}}^{-z\ell }};\delta a(\ell ),b(\ell ),\lambda (\ell ) \right),\end{eqnarray} \tag{ 27 }$$

where we have not displayed ${{\mu }_{0}}$ and D₀, since they are kept fixed in the RG. By choosing $\ell ={\rm ln} (\Lambda r)$, and using $\delta a(\ell )\approx \delta {{a}_{0}}{{{\rm e}}^{{{y}_{a}}\ell }}$, we can obtain from this a scaling form for large r:

$$\begin{eqnarray}&&C({\bf r},t)={{r}^{2-d-\eta }}{{Y}_{\pm }}\left( \frac{r}{\xi },\frac{t}{{{r}^{z}}} \right).\end{eqnarray} \tag{ 28 }$$

In particular, the equal time correlation function scales as ${{r}^{2-d-\eta }}$. In equation (28), the exponent $\eta$ is given by

$$\begin{eqnarray}&&\eta =2-d-2\chi =2-z+\mathcal{O}({{\epsilon }^{2}})=\frac{31}{113}\epsilon +\mathcal{O}({{\epsilon }^{2}}),\end{eqnarray} \tag{ 29 }$$

the scaling functions ${{Y}_{\pm }}$ by

$$\begin{eqnarray}&&{{Y}_{\pm }}(x,y)=C({{\Lambda }^{-1}},y;\pm {{x}^{{{y}_{a}}}},{{b}^{*}},{{\lambda }^{*}}),\end{eqnarray} \tag{ 30 }$$

where b* and ${{\lambda }^{*}}$ are the nonzero fixed values of $b(\ell )$ and $\lambda (\ell )$, respectively, and the diverging correlation length $\xi$ by $\xi \equiv {{\Lambda }^{-1}}|\delta {{a}_{0}}{{|}^{-\nu }}$, where the correlation length exponent

$$\begin{eqnarray}&&\nu =\frac{1}{{{y}_{a}}}=\frac{1}{2}+\frac{29}{226}\epsilon +\mathcal{O}({{\epsilon }^{2}}).\end{eqnarray} \tag{ 31 }$$

In our expression for the expansions for $\eta$, we have replaced $\chi$ and z by their values at the fixed point (18), which is valid given that r is large and the system is sufficiently close to the transition. We do so consistently when calculating other exponents as well. The first line of equation (29) is exact (i.e., independent of the -expansion), as it is simply the definition of $\eta$.

As discussed earlier, the order-disorder transition can be driven by tuning any one of many microscopic control parameters (e.g., density or noise strength). Whatever control parameter s is tuned, we expect ${{a}_{0}}-a_{0}^{c}\propto \left( s-{{s}_{{\rm c}}} \right)$ by analyticity near s_c, where s_c is the value of the control parameter s at the transition. As a result, the velocity correlation length $\xi$ just defined diverges as $\xi \propto |s-{{s}_{{\rm c}}}{{|}^{-\nu }}$ as any control parameter s is tuned.

Now we calculate the magnitude of the order parameter in the ordered state near the critical point. The RG connects the average velocity of the original system and that of the rescaled system with the relation

$$\begin{eqnarray}&&\langle {\bf v}\rangle (\delta {{a}_{0}},{{b}_{0}},{{\lambda }_{0}})={{{\rm e}}^{\chi \ell }}\langle {\bf v}\rangle \left( \delta {{a}_{0}}{{{\rm e}}^{{{y}_{a}}\ell }},b(\ell ),\lambda (\ell ) \right).\end{eqnarray} \tag{ 32 }$$

This relation holds for all ℓ. To make the best possible use of it, we will choose ℓ such that $\delta {{a}_{0}}{{{\rm e}}^{{{y}_{a}}\ell }}$ is of order 1. Therefore, ℓ is large since $\delta {{a}_{0}}$ is small near the critical point, and hence, both $b(\ell )$ and $\lambda (\ell )$ flow to their nonzero fixed values. Therefore, the quantity $\langle {\bf v}\rangle \left( \delta {{a}_{0}}{{{\rm e}}^{{{y}_{a}}\ell }},b(\ell ),\lambda (\ell ) \right)\approx \langle {\bf v}\rangle \left( \mathcal{O}(1),{{b}^{*}},{{\lambda }^{*}} \right)$ should be $\mathcal{O}(1)$, and is in any case independent of the bare parameters, in particular a₀. Hence, all the singular dependence on $({{s}_{{\rm c}}}-s)$ on the rhs of the equality (32) (indeed, all of the dependence on the bare parameters) is contained in the exponential (which depends in particular on $\delta {{a}_{0}}$ since the choice of ℓ that makes $\delta {{a}_{0}}{{{\rm e}}^{{{y}_{a}}\ell }}=\mathcal{O}(1)$ obviously does). This implies $|\langle {\bf v}\rangle |\sim |{{s}_{{\rm c}}}-s{{|}^{\beta }}$ with

$$\begin{eqnarray}&&\beta =-\nu \chi =\frac{1}{2}-\frac{6}{113}\epsilon +\mathcal{O}({{\epsilon }^{2}}).\end{eqnarray} \tag{ 33 }$$

The first equality in this expression, which is exact, can be rewritten in terms of $\eta$ using the definition of $\eta$ embodied in the first equality of (29), giving the exact hyperscaling relation

$$\begin{eqnarray}&&\beta =\frac{\nu }{2}(d-2+\eta ).\end{eqnarray} \tag{ 34 }$$

In this respect, our system is similar to equilibrium systems, in which (34) also holds [25].

We study next the response of the system to a weak external field ${\bf H}$; that is, simply adding a small constant vector ${\bf H}$ to the rhs of (1). In this case, the RG leads to the scaling relation

$$\begin{eqnarray}&&\langle {\bf v}\rangle (\delta {{a}_{0}},{{b}_{0}},{{\lambda }_{0}},H)={{{\rm e}}^{\chi \ell }}\langle {\bf v}\rangle \left( \delta {{a}_{0}}{{{\rm e}}^{{{y}_{a}}\ell }},b(\ell ),\lambda (\ell ),H{{{\rm e}}^{{{y}_{H}}\ell }} \right),\end{eqnarray} \tag{ 35 }$$

where $H\equiv |{\bf H}|$ and y_H is the RG eigenvalue of the external field H at the fixed point (18). As there are no one loop graphical corrections to the external field, we can obtain y_H to $\mathcal{O}(\epsilon )$ by simple power counting, which gives

$$\begin{eqnarray}&&{{y}_{H}}=z-\chi +\mathcal{O}({{\epsilon }^{2}}).\end{eqnarray} \tag{ 36 }$$

Again choosing ℓ such that $\delta {{a}_{0}}{{{\rm e}}^{{{y}_{a}}\ell }}$ is of order 1, we obtain

$$\begin{eqnarray}&&\langle {\bf v}\rangle (\delta {{a}_{0}},{{b}_{0}},{{\lambda }_{0}},H)={{\left( \Lambda \xi \right)}^{\chi }}\langle {\bf v}\rangle \left( 1,{{b}^{*}},{{\lambda }^{*}},H{{\left( \Lambda \xi \right)}^{{{y}_{H}}}} \right).\end{eqnarray} \tag{ 37 }$$

Since the expectation value on the right-hand side is evaluated in a system far from its critical region (since $\delta a=1$), we expect linear response to the external field on that side with an order one susceptibility. Hence

$$\begin{eqnarray}&&\langle {\bf v}\rangle (\delta {{a}_{0}},{{b}_{0}},{{\lambda }_{0}},H)\sim {{\left( \Lambda \xi \right)}^{\chi }}H{{\left( \Lambda \xi \right)}^{{{y}_{H}}\ell }}\propto {{\xi }^{\chi +{{y}_{H}}}}H,\end{eqnarray} \tag{ 38 }$$

which implies a linear susceptibility ${{\chi }_{H}}$ which diverges as $|s-{{s}_{{\rm c}}}{{|}^{-\gamma }}$ with

$$\begin{eqnarray}&&\gamma =\nu (\chi +{{y}_{H}})=\nu z+\mathcal{O}({{\epsilon }^{2}})=1+\frac{27}{226}\epsilon +\mathcal{O}({{\epsilon }^{2}}).\end{eqnarray} \tag{ 39 }$$

Note that equations (29), (39) seem to suggest that $\eta$ and $\gamma$ satisfy Fisher's scaling law $\gamma =(2-\eta )\nu$. However, since our system is out of equilibrium and thus the fluctuation dissipation theorem is not expected to hold, we do not expect Fisher's scaling law to hold; indeed, the $\mathcal{O}({{\epsilon }^{2}})$ terms probably violate it. The first line of equation (39) is exact, however, and can be used to derive another scaling law, as we will now show.

Turning on a small field right at the transition, we can again relate then the average velocities of the original and the rescaled systems using equation (35). However, $\delta a(\ell )$ now flows to 0 for large ℓ since the system is right at the critical point. Therefore, by choosing $\ell ={\rm ln} \left( 1/H \right)/{{y}_{H}}$, we obtain the H-dependence of ${\bf v}$:

$$\begin{eqnarray}&&\langle {\bf v}\rangle (\delta {{a}_{0}},{{b}_{0}},{{\lambda }_{0}},H)={{H}^{-\frac{\chi }{{{y}_{H}}}}}\langle {\bf v}\rangle \left( 0,{{b}^{*}},{{\lambda }^{*}},1 \right).\end{eqnarray} \tag{ 40 }$$

Again, all of the dependence on bare model parameters is now contained in the prefactor ${{H}^{-\frac{\chi }{{{y}_{H}}}}}$. Comparing (40) with our definition (24) of the critical exponent δ reveals that $\delta =-\frac{{{y}_{H}}}{\chi }$. Combining this with the exact first equalities in equations (33), (39), we obtain Widom's scaling relation

$$\begin{eqnarray}&&\gamma =\beta (\delta -1),\end{eqnarray} \tag{ 41 }$$

which is exact. Plugging the -expansions of $\gamma$ and $\beta$ into this relation, we find

$$\begin{eqnarray}&&\delta =3+\frac{51}{113}\epsilon +\mathcal{O}({{\epsilon }^{2}}).\end{eqnarray} \tag{ 42 }$$

We can estimate the numerical values of the exponents in spatial dimension d = 3 as follows: We first choose a scaling relation satisfied by any three exponents (e.g., equation (23) for $\eta$, $\beta$, and ν). We then determine numerical values for two of them (e.g., ν and $\beta$) by simply setting $\epsilon =1$ in the -expansion for them, and dropping the unknown $\mathcal{O}({{\epsilon }^{2}})$ terms. We now get the value of the third exponent (e.g., $\eta$) by requiring that the scaling law (i.e., equation (23) ) hold exactly in d = 3. In this example, this gives $\beta =1/2-6/113\approx 0.447$, $\nu =1/2+29/226\approx 0.628$, and $\eta =2\beta /\nu -1\approx 0.424$. Next, we take $\eta$ and $\nu$ to be given by their respective -expansions with $\epsilon =1$, and get $\beta$ from the exact scaling relation. This gives: $\nu =1/2+29/226\approx 0.628$, $\eta =31/113\approx 0.274$, and $\beta =\nu (1+\eta )/2\approx 0.400$. Finally, we take $\beta$ and $\eta$ from their -expansions, and get $\nu$ from the exact scaling relation, obtaining $\eta =31/113\approx 0.274$, $\beta =1/2-6/113\approx 0.447$, and $\nu =2\beta /(1+\eta )=0.702$.

Note that each exponent gets two possible values in this approach: one from directly setting $\epsilon =1$ in the -expansion, and another by obtaining the exponent from the exact scaling relation in d = 3.

Applying the same approach to Widom's exact scaling relation (i.e., equation (30)) and the three associated exponents $\gamma$, $\beta$, and δ gives the possible values: $\beta \approx 0.447$ or $\beta \approx 0.457$, $\delta \approx 3.451$ or $\delta \approx 3.503$, and $\gamma \approx 1.119$ or $\gamma \approx 1.096$.

So if we look at the range of values we have found for each of the exponents, we have $0.274\leqslant \eta \leqslant 0.424$, $0.628\leqslant \nu \leqslant 0.702$, $0.400\leqslant \beta \leqslant 0.457$, $3.451\leqslant \delta \leqslant 3.503$, and $1.096\leqslant \gamma \leqslant 1.119$. Assuming, as seems reasonable (and as is true for, e.g., the critical exponents for the equilibrium $\mathcal{O}(n)$ model [25]), that the correct values lie within the range spanned by the different approaches we have used here, we can conclude that, in spatial dimension d = 3, the critical exponents are as shown in the second column of table 1.

Comparing the critical exponents with the known values for the two equilibrium analogs of our system: the three-dimensional, three component Heisenberg model (i.e., the $\mathcal{O}(3)$ model) with and without dipolar interactions (third and fourth columns respectively in table 1), we see that $\nu$ and $\beta$ are very close in all three models. The situation is a little better for $\gamma$ and δ. The biggest difference, however, is clearly in $\eta$, which is much larger in the incompressible active fluid. Thus experiments to determine this exponent, which, as can be seen from equation (17), can be deduced from velocity correlations right at the critical point, will provide the clearest and most dramatic evidence for the non-equilibrium nature of this system, and the novelty of its universality class.

The values of the exponents in d = 2 obviously can not be reliably estimated quantitatively from the $4-\epsilon$expansion. We do note, however, that the ordered state is expected to exist and to have true long-ranged order. This is clear since true long-ranged order exists even in the compressible problem, which obviously has more fluctuations than the incompressible problem we have studied here. Hence, we do not expect this problem to be like the equilibrium $2d$ XY model, in which [28] the ordered state only has quasi-long-ranged order (i.e., algebraically decaying correlations). We therefore do not expect 2d incompressible active fluids to exhibit any of the singular behaviour of exponents found in the 2d equilibrium XY model; in particular, there is no reason to expect $\nu =\infty$. Beyond this, there is little we can say quantitatively about d = 2 beyond the expectation that the critical exponents $\beta$, $\nu$, $\eta$, δ, and $\gamma$ should be further from their mean field values $\beta =\nu =1/2$, $\eta =0$, $\delta =3$, and $\gamma =1$ than they are in d = 3. This implies that in d = 2, $\beta$ will be smaller, and the four other exponents will be bigger, than the values quoted above for d = 3. We also note that the exact scaling relations equations (34), (41) will hold in d = 2, and that all of the exponents will be universal (i.e., the same for all incompressible active fluids) in d = 2.

We have studied the order-disorder transition in incompressible active fluids using a dynamical $\epsilon =4-d$ expansion. This is the first study of the static to moving phase transition in active matter to go beyond mean-field theory, and include the effects of fluctuations on the transition. We find a stable non-Gaussian fixed point, which implies a continuous transition, whose critical exponents were calculated to $\mathcal{O}(\epsilon )$. This fixed point is new, and all of the critical exponents differ from those found for any previously known phase transition. Therefore, the universality class of this transition is new. In addition, we found that among the five critical exponents we calculated, there are two exact scaling relations which are the same as those in equilibrium ferromagnetic transitions. This is despite the fact that our system is fundamentally nonequilibrium, and that the universality class of its transition is new. Furthermore, we connected our model with two classic universality classes discovered in the early days of the renormalization group, namely the randomly stirred fluid model [23], and the dipolar ferromagnet model [26] (see figure 3). Specifically, we now know that the two fixed points associated to the two classic universality classes are unstable in the combined model we considered here.

Our predictions on the critical exponents can be tested in motile systems (both experimental and simulated) with strong repulsive interactions between the particles, e.g., in bacterial suspensions like those studied in [19, 29], provided that the system can be tuned to reach the critical point. The exponent $\beta$ can be determined by measuring the average velocity. The exponents $\eta$ and the correlation length (which determines $\nu$) can be obtained by the velocity correlation functions. The exponents $\gamma$ and δ can be obtained by measuring the response of the average velocity to an external perturbation. In magnetotactic bacteria; the perturbation can be a magnetic field [30], in chemotactic bacteria, it can be a nutrient gradient [31].

Future theoretical work [32] on this problem will include working out in quantitative detail the cutoff of the continuous transition by the banding instability in systems with a small, but non-zero, compressibility.

JT thanks Nicholas Guttenberg for explaining how to simulate systems with long-ranged interactions. He also thanks the Department of Physics, University of California, Berkeley, CA; the Max Planck Institute for the Physics of Complex Systems (MPI-PKS), Dresden, Germany; the Aspen Center for Physics, Aspen, CO; the Isaac Newton Institute, Cambridge, UK; the Kavli Institute for Theoretical Physics, Santa Barbara, CA; and the Department of Bioengineering, Imperial College, London, UK for their hospitality while this work was underway. He also thanks the US NSF for support by awards #EF-1137815 and #1006171; and the Simons Foundation for support by award #225579. LC acknowledges support by the National Science Foundation of China (under Grant No. 11474354) and the Fundamental Research Funds for the Central Universities (under Grant No. 2013XK04). CFL thanks the MPI-PKS where the early stage of this work was performed.