Critical Phenomenon of the Order-Disorder Transition in Incompressible Flocks

We study incompressible systems of motile particles with alignment interactions. Unlike their compressible counterparts, in which the order-disorder (i.e., moving to static) transition, tuned by either noise or number density, is discontinuous, in incompressible systems this transition can be continuous, and belongs to a new universality class. We calculate the critical exponents to $O(\epsilon)$in an $\epsilon=4-d$ expansion, and derive two exact scaling relations. This is the first analytic treatment of a phase transition in a new universality class in an active system.

Emergent properties of interacting non-equilibrium systems are of widespread and fundamental interest. One of the simplest, but most striking, of these is the selforganized phenomenon of "flocking"-that is, collective motion (CM) in large groups of motile organisms [1][2][3][4][5][6][7][8][9][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][13][14][15][16][17][18] of the hydrodynamic equations [3][4][5][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 flock 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, it can be eliminated by eliminating density fluctuations: that is, by making system incompressible. In this paper, we show that the order-disorder transition is continuous in an incompressible system. We demonstrate this by finding, in a dynamical renormalization group (DRG) analysis of the hydrodynamic equations for an incompressible flock in d spatial dimensions, a stable fixed point that controls * Electronic address: leimingyandongyu@gmail.com † Electronic address: jjt@uoregon.edu ‡ Electronic address: c.lee@imperial.ac.uk the transition. This calculation is done to order O( ) in an = 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: 1) 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 behavior (i.e., a small first order transition driven by the banding instability). 2) 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].
3) 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][13][14][15][16][17][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 fluctu-ate (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][4][5][6]23] where the pressure P enforces the incompressibility condition ∇ · v = 0, f is a "white noise" with spatiotemporally Fourier transformed statistics: and P mn (k) ≡ δ mn −k m k n /k 2 is the transverse projection operator. This EOM (Eq. (1)) reduces, when a = 0 = b, to the classic model of a fluid forced at zero wavenumber treated by [24] (their "model B"). With λ = 0, it reduces to a simple, time-dependent Ginzburg-Landau (TDGL) [1,25] dynamical model for an isotropic ferromagnet with long ranged dipolar interactions [2,28].
Because our system lacks Galilean invariance, λ need not (and in general will not) be one, and the terms −(a+ b|v| 2 )v are allowed in the EOM. The latter is crucial as it explains why there can be a polar ordered phase in an active system, which is not possible in a normal fluid.
At the mean field level, for a < 0, b > 0 the system is in the ordered phase with |v| = −a/b, and for a > 0, b > 0 it is in the disordered phase with |v| = 0. To go beyond this mean field description, we employ the DRG method [24] near the order-disorder transition. To do so, we spatio-temporally Fourier transform Eq. (1), and project orthogonal to wavevector k; obtaining where we have adopted the reduced notationsk ≡ (k, ω) and q = q,Ω ≡ , and we have defined P lmn (k) ≡ P lm (k)k n + P ln (k)k m , Q lmnp (k) ≡ P lm (k)δ np + P ln (k)δ mp + P lp (k)δ mn , and the "propagator" G(k) ≡ (−iω+µk 2 +a) −1 . Graphical representations of the various terms in Eq. (3) are shown in Fig. 1.
We now perform the standard DRG procedure [24], averaging over short wavelength degrees of freedom, and rescaling: r → re , t → te z and v → e χ v. Our procedure is identical to the calculation for model B in [24], except for a modified propagator, and some additional Feynmann graphs due to the extra b|v| 2 v non-linearity in our problem. At the one loop level, the non-vanishing graphical contributions to the various coefficients in Eq.
(3) are shown in Fig. 2. More details of the calculation are given in the supplemental materials. We obtain, in d spatial dimensions, the following RG flow equations of the coefficients to one loop order and to linear order in where we've defined dimensionless couplings: and where S d ≡ 2π d/2 /Γ(d/2) is the surface area of a unit sphere in d dimensions, ≡ 4 − d, and Λ is the ultraviolet wavevector cutoff. Since our interest is in the transition, we have, in the last four recursion relations (5)(6)(7)(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 = 4 − d.
From these RG flow equations, we can derive two closed flow equations for g 1,2 for arbitrary χ and z: Although not necessary, it is convenient to make a special choice of z and χ such that µ and D are kept fixed at their bare values (i.e, µ 0 and D 0 , respectively). We will hereafter adopt this choice of z and χ, which is We will also hereafter use the subscript 0 to denote the bare (i.e., unrenormalized) values of the parameters. Eq. (4) now becomes Eqs. (10,11,13) have a non-Gaussian fixed point in d < 4: which can be shown by analyzing the three recursion relations to be a stable attractor of all points on a twodimensional surface (the "critical surface" ) in the threedimensional 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 Fig.  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. The exponential runaway from the critical surface in the unstable direction near the stable fixed point (14) grows like e ya , with the exponent This eigenvalue determines the critical exponent ν governing the critical behavior of the velocity correlation length ξ. In addition, the smallest (in magnitude) of the two negative eigenvalues gives the "correction to scaling exponent" y 2 [25]; we find y 2 = − 31 113 + O( 2 ). A useful experimental probe of the transition is the velocity correlation function v(r + R, t + T ) · v(R, T ) ≡ C(r, t), which depends on bare parameters b 0 , µ 0 , D 0 , and λ 0 and, most importantly, the proximity to the phase transition δa 0 ≡ a 0 − a c 0 , where a c 0 is the value of a 0 at the transition. The RG connects the original C(r, t) to that of the rescaled system: where we have not displayed µ 0 and D 0 , since they are kept fixed in the RG. By choosing = ln(Λr), and using δa( ) ≈ δa 0 e ya , we can obtain from this a scaling form for large r: In particular, the equal time correlation function scales as r 2−d−η . In Eq. 17, the exponent η is given by the scaling functions Y ± by and the diverging correlation length ξ by ξ ≡ Λ −1 |δa 0 | −ν , where the correlation length exponent In our expression for the expansions for η, we have replaced χ and z by their values at the fixed point (14), Besides the unstable Gaussian fixed point (black diamond) and the stable fixed described in Eqs (14) (red square), there are two unstable fixed points: one at g * 1 = 0, g * 2 = 2 17 , which is the fixed point of an isotropic ferromagnet with long-ranged dipolar interactions [2] (purple circle), and one at g * 2 = 0, g * 1 = 4 3 , which is the fixed point of a fluid forced at zero wavevector (Model B of [24]) (blue triangle).
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 (18) is exact (i.e., independent of the -expansion), as it is simply the definition of η.
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 c 0 ∝ (s − s c ) by analyticity near s c , where s c is the values of the control parameter s at the transition. As a result, the velocity correlation length ξ just defined diverges as ξ ∝ |s − s c | −ν as any control parameter s is tuned.
The scaling functions Y ± in equation (17) 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 behavior that we'll discuss elsewhere [30]. 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 We choose such that δa 0 e ya is of order 1. Therefore, is large since δa 0 is small near the critical point, and hence, both b( ) and λ( ) flow to their nonzero fixed values. Then all the singular dependence on (s c − s) on the RHS of the equality (21) is included in the exponential. This The first equality in this expression, which is exact, can be rewritten in terms of η using the definition of η embodied in the first equality of (18), giving the exact hyperscaling relation In this respect, our system is similar to equilibrium systems, in which (23) also holds [1]. We study next the linear response of the system to a weak external field H; that is, simply adding a small constant vector H to the RHS of (1). In this case, the RG leads to the scaling relation where H ≡ |H| and y H is the RG eigenvalue of the external field H at the fixed point (14). As there are no one loop graphical corrections to the external field, we can obtain y H to O( ) by simple power counting, which gives Again choosing such that δa 0 e ya is of order 1, we obtain where b * and λ * are the nonzero fixed values of b( ) and λ( ), respectively. Since the expectation value on the right hand side is evaluated in a system far from its critical region (since δa = 1), we expect linear response to the external field on that side with an order one susceptibility. Hence, which implies a linear susceptibility χ H which diverges as |s − s c | −γ with Note that Eqs. (18,28) seem to suggest that η and γ satisfy Fisher's scaling law γ = (2 − η)ν. 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 O( 2 ) terms probably violate it. The first line of equation (28) is exact, however, and can be used to derive another scaling law, as we'll 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 Eq. (24). However,δa( ) now flows to 0 for large since the system is right at the critical point. Therefore, by choosing = ln This implies δ = − y H χ . Combining this with the exact first equalities in Eqs (22,28), we obtain Widom's scaling relation which is exact. Plugging the -expansions of γ and β into this relation, we find In summary, we have studied the order-disorder transition in incompressible flocks using a dynamical = 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 O( ). This fixed point is new, and so, therefore, is the universality class of this transition. In addition, we found that the critical exponents obey two exact scaling relations which are the same as those in equilibrium ferromagnetic transitions, despite the fact that our system is fundamentally nonequilibrium. We also presented predictions for the scaling behavior of the velocity correlation functions. Future theoretical work [30] 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.
J.T. thanks Nicholas Guttenberg for explaining how to simulate systems with long ranged interactions; and the Max Planck Institute for the Physics of Complex Sys-tems, Dresden; the Aspen Center for Physics, Aspen, Colorado; the Isaac Newton Institute, Cambridge, U.K., and the Kavli Institute for Theoretical Physics, Santa Barbara, California; 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. L.C. acknowledges support by the National Science Foundation of China (under Grant No. 11474354).

Supplmental Material
The various constituent elements of Feynman diagrams are illustrated in Fig. 1 in the main text. The one loop graphical corrections to various coefficients in the model Eq. (1) in the main text are illustrated in Fig. 2 in the main text. Since we are only interested in the RG flows near the critical point, we have evaluated all of these graphs at a = 0, except for the graph (a) for a, the adependence of which we need to determine the "thermal" eigenvalue y a . The corrections from graphs (a), (b), (c), and (d) are given respectively by: Combining these corrections (including the zero correction to D ) with the rescalings described in the main text, and setting dimension d = 4 in all of these expressions (which is sufficient to obtain results to first order in = 4 − d), leads to the recursion relations (4)(5)(6)(7)(8) in the main text, (although in (4) we have in addition expanded to linear order in a, which is sufficient to determine the exponents to O( )).

I. EVALUATION OF THE FEYNMAN DIAGRAMS
A. Graph (a) This graph represents an additional contribution δ (∂ t v l ) to ∂ t v l given by: where P mn (q) = δ mn − q m q n /q 2 is the transverse projection operator, and Q lmnp is defined in the main text after Eq. (3). In going from the second to the third line above, we have used the well-known identity P mn (q) q = 1 − 1 d δ mn , where q denotes the average over directionsq of q for fixed |q|. This identity is derived in part (II) of these Supplemental Materials. Clearly this correction to ∂ t v l (k, ω) is exactly what one would obtain by adding to the parameter a (hiding in the "propagator") in Eq. (3) in the main text a correction δa given by the coefficient of v l in (36); i.e., Eq. (32).

B. Graph (b)
This graph represents an additional contribution δ (∂ t v l ) to ∂ t v l given by: where the combinatoric prefactor of 18 arises because there are 3 ways to pick the leg with index o on the left. Then, once this choice has been made, there are 2 ways to pick the leg with index m on the left, and 3 ways to pick the one with index i on the right. To lowest order in the external momenta (k, p and h) and frequencies ω, ω p , and ω h , we can set all of these external momenta and frequencies equal to zero in the integrand of the integral overq in Eq. (37). This proves, as we shall see, to make this graph into a renormalization of the cubic non-linearity b. Making this simplification, Eq. (37) becomes, after using and dropping an integral of an odd function of Ω, Using the definition of Q nijs (q), and contracting indices, we obtain P mi (q)Q nijs (q) = P mn (q)δ js + P ms (q)P jn (q) + P jm (q)P ns (q) .
As we've done repeatedly throughout these Supplemental Materials, we'll replace this tensor with its average over all directions ofq. This is easily obtained from the known direction averages (70) and (84), and gives, after a little algebra, Inserting this into our earlier expression for δ (∂ t v l ) and performing a few tensor index contractions gives The trace in this expression can be evaluated as and the integral over frequency Ω and q is readily evaluated, and is given by Putting (43) and (44) into (42), and taking advantage of the complete symmetry of p, under interchanges of the indices j, o, and s to symmetrize the tensor prefactor gives This is readily recognized as a contribution to the − b 3 term in Eq. (3) in the main text; hence, the correction to b is given by (33).

C. The First Graph in (c)
This graph represents an additional contribution δ (∂ t v l ) to ∂ t v l given by: where P mn (k) is defined in the main text after Eq. (3), and the combinatoric prefactor of 6 arises because there are 2 ways to pick the leg with index m on the left, and 3 ways to pick the one with index i on the right. The piece of δ (∂ t v l ) linear in the external momentum k is immediately recognized as a contribution to the − i 2 λ term in Eq. (3) in the main text. Since there is already an implicit factor of k in the P lmn (k) in this expression, we can evaluate this graph to linear order in k by setting both k and the external frequency ω to zero in the integrand of the integral overq. Doing so, and in addition using (38) gives, after dropping an integral of an odd function of Ω that vanishes, where we've used the fact that Q nijs (q) is an even function of q. Using the definition of Q nijs (q), and performing the tensor index contractions, we can simplify the numerator of the integrand as follows: P mi (q)Q nijs (q) = P mn δ js + P jn P ms + P ns P mj , (48) where we've also used the fact that, e.g., P mi (q)P in (q) = P mn (q) (which is a consequence of the definition of P mi (q) as a projection operator). Using this result (48), and taking the angle average of the numerator of (47) (which is the only factor in the integral that depends on the direction of q) gives, In deriving this expression, we've made liberal use of the angle averages (84) and (70). The first term on the right hand side of (49) contributes nothing, since it contracts two of the indices on the prefactor P lmn (k) together, which gives zero, as can be seen from the definition of P lmn (k): Keeping only the second term gives where in the last step we have used the symmetry of P lmn (k) under interchange of its last two indices. This is immediately recognized as a contribution to the − i 2 λ term in Eq. (3) in the main text, which implies a correction to λ given by which is just the first term on the RHS of Eq. (34).

D. The Second Graph in (c)
This graph represents an additional contribution δ (∂ t v l ) to ∂ t v l given by where the combinatoric prefactor of 12 arises because there are 3 ways to pick the leg with index o on the left. Then, once this choice has been made, there are 2 ways to pick the leg with index m on the left, and 2 ways to pick the one with index i on the right. To extract from (53) the contribution to the − i 2 λ term in Eq. (3) in the main text, we need the pieces of δ (∂ t v l ) that are linear in either the external momenta k or h. We notice that the integral overq vanishes if we set h = 0 in its integrand. This implies the entire term is at least of order h. Thus, to obtain a correction to λ we can simply set the external frequency ω = 0 in the integrand; doing so, and integrating over Ω, we obtain This can be rewritten as jmn (h) + I where I (2) jmn (h) is given by equation (79) (with the obvious substitution k → h), and the other integrals are defined as: Immediately, by the properties of the projection operator we get We notice that both I jmn (h) and I (4) jmn (h) are already proportional to h, so we can simply set h = 0 inside the integral. Thus, we obtain Plugging the values (79), (60), (61), and (59) of the various integrals into Eq. (55), we obtain The third piece vanishes due to the incompressibility condition h j v j (h) = 0. The first and the second pieces can be grouped together since Q lmno is invariant under the interchange of m and n. Therefore, the above expression can be simplified as where in the second equality we have dropped the second and third pieces. The second piece can be dropped because it vanishes, as can be seen by simply changing variables of integration from h to k − h; this gives Adding the left and the right hand side of this equation, and dividing by 2, implies with the last equality following from P n (k)k n = 0. The third piece can be dropped due to the incompressibility condition h j v j (h) = 0. The remaining piece of (63) is readily recognized as a contribution to the − i 2 λ term in Eq. (3) in the main text. This implies a correction to λ given by the second piece on the RHS of Eq. (34).

E. Graph (d)
This graph represents an additional contribution δ (∂ t v l ) to ∂ t v l given by: The integral in this expression is readily seen to vanish when k → 0, since the integrand then becomes odd in q. Hence, the integral is at least of order k, so the entire term (include the implicit first power of k coming from the P lmn (k) in front), is O(k 2 ). Since we do not need to keep any terms in the equations of motion higher order in k and ω than O(k 2 ), this means that we can safely set ω = 0 inside the integral. Keeping just the O(k) piece of the integral then gives us a modification to the equation of motion of O(k 2 v), which is clearly a renormalization of the diffusion constant µ. So setting ω = 0 in the integrand for the reasons just discussed, and then writing P nij (k − q) using its definition as given in the main text after Eq. (3), gives for the integral in (65): The term proportional to k j in this expression can be dropped, since k j v j = 0 (this is just the incompressibility condition ∇ · v = 0 written in Fourier space). The term proportional to q i can be dropped since P mi (q)q i = 0 by the properties of the transverse projection operator P mi (q). This leaves two terms in the integral, which can be written as and Since I (1) jmn (k) already has an explicit factor of k in front, we can evaluate it to linear order in k by setting k = 0 inside the integral. Doing so gives The integral in this expression can now be evaluated by replacing the only piece that depends on the direction q of q, namely, the factor P mi (q)P jn (q), with its angle average. As shown in (II), this average is given by Inserting this into (69) gives Now let us expand I (2) jmn (k) to linear order in k. Changing variables of integration from q to a shifted variable p defined by: gives where we've defined and I jmn (k) to be the first and second terms in the expression for I (2) jmn (k). Since I (2.2) jmn (k) has an explicit factor of k in front, we can evaluate this term to linear order in k by setting k = 0 inside the integral. This leads to The calculation of I (2.1) jmn (k) requires more effort. Expanding the numerator we get dition k j v j = 0. Dropping them, and adding these two integrals I (1) jmn (k) and I (2) jmn (k) makes the entire correc-tion to the equation of motion coming from graph II become: Simplifying, and performing the tensor index contractions, gives where in the second step we have used the symmetry of P ljn (k) to write P ljn (k)k n = P lmj (k)k m . Now from the definition of P lmj (k), we have P lmj (k)k m = P lm (k)k j k m + P lj (k)k m k m . The first term in this ex-pression vanishes by the fundamental property of the transverse projection operator, while the second is just k 2 P lj (k). Thus, we finally obtain which is exactly what one would get by adding to µ (hiding in the "propagator") in Eq. (3) in the main text a correction δµ given by Eq. (35).

F. Vanishing diagrams
Besides the five non-vanishing one-loop diagrams, there are also two sets of one-loop diagrams that cancel exactly, giving zero contribution to the corrections (Fig. 4).

II. ANGULAR AVERAGES
In this section we derive the various angular averages used in the previous sections.
We begin by deriving the identity where q denotes the average over directionsq of q for fixed |q|. The identity (83) follows by symmetry: the average in question clearly must vanish when m = n, since then the quantity being averaged is odd in q. Furthermore, when m = n, the average must be independent of the value that m and n both equal. Hence, this average must be proportional to δ mn . The constant of proportionality in (83) is easily determined by noting that the trace of this average over mn is qmqm q 2 q = q 2 q 2 q = 1 q = 1. This forces the prefactor of 1 d in (83). From (83), it obviously follows that which is the identity we used in equation (36).
We now consider the average of two projection operators P mi (q)P jn (q) q that appears in (70). Using the definition of the projection operator, this can be written as follows: Exponents Incompressible active matter Heisenberg model [2] Heisenberg model with dipolar interactions [3] β 0.
The first two angular averages on the right hand side of this expression can be read off from (83). The last is new, and can be evaluated as follows: First, note that by symmetry, this average vanishes unless the four indices i, j, k, l, are equal in pairs. Furthermore, if they are equal in pairs, but the pairs are different (e.g., if i = j = x and m = n = z), then the average will have one value, independent of what the values of the two pairs of indices are (e.g., if i = j = y and m = n = x, the average would be the same as in the example just cited. The only other non-zero possibility is that all four indices are equal, in which case the average is the same no matter which index all four are equal to (i.e., the average when i = j = m = n = x is the same as that when i = j = m = n = z). Furthermore, this average must be completely symmetric under any interchange of its indices. This can be summarized by saying that the average must take the form: where Υ ijmn = 1 if and only if i = j = m = n, and is zero otherwise, and A and B are unknown, dimension (d) dependent constants that we'll now determine. We can derive one condition on A and B by taking the trace of (86) over any two indices (say, i and j). This gives q m q n q 2 q = (A + (d + 2)B)δ mn .
Comparing this with (83) gives A second condition can be derived by explicitly evaluating the angle average when all four indices are equal.
Since it doesn't matter what value they all equal, we'll chose it to be z. Defining θ to be the angle between the z-axis and q, we can obtain the needed average in ddimensions by integrating in hyperspherical coordinates: .
Comparing this with (86) evaluated for i = j = m = m gives .

III. NUMERICAL ESTIMATES OF THE CRITICAL EXPONENTS
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., Eq. (23) in the main text for η, β, and ν). We then determine numerical values for two of them (e.g., ν and β) by simply setting = 1 in the -expansion for them, and dropping the unknown O( 2 ) terms. We now get the value of the third exponent (e.g., η) by requiring that the scaling law (i.e., Eq. (23) in the main text ) hold exactly in d = 3.
Note that each exponent gets two possible values in this approach: one from directly setting = 1 in theexpansion, and another by obtaining the exponent from the exact scaling relation in d = 3.
So if we look at the range of values we've found for each of the exponents, we have .274 ≤ η ≤ .424, .628 ≤ ν ≤ .702, .400 ≤ β ≤ .457, 3.451 ≤ δ ≤ 3.503, and 1.096 ≤ γ ≤ 1.119. Assuming, as seems reasonable (and as is true for, e.g., the critical exponents for the equilibrium O(n) model [1]), that the correct values lie within the range spanned by the different approaches we've used here, we can conclude that, in spatial dimension d = 3, the critical exponents are as shown in the second column of Table I.
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 O(3) model) with and without dipolar interactions (third and fourth columns respectively in Table  I), we see that ν and β are very close in all three models.
The situation is a little better for γ and δ. The biggest difference, however, is clearly in η, which is much larger in the incompressible flock. Thus experiments to determine this exponent, which, as can be seen from equation (17) in the main text, 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 −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've studied here. Hence, we do not expect this problem to be like the equilibrium 2d XY model, in which [4] the ordered state only has quasi-long-ranged order (i.e., algebraically decaying correlations). We therefore do not expect 2d incompressible flocks to exhibit any of the singular behavior of exponents found in the 2d equilibrium XY model; in particular, there is no reason to expect ν = ∞. Beyond this, there is little we can say quantitatively about d = 2 beyond the expectation that the critical exponents β, ν, η, δ, and γ should be further from their mean field values β = ν = 1/2, η = 0, δ = 3, and γ = 1 than they are in d = 3. This implies that in d = 2, β 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 Eqs. (23,30) in the main text will hold in d = 2, and that all of the exponents will be universal (i.e., the same for all incompressible flocks) in d = 2.