Evidence of fluctuation-induced first-order phase transition in active matter

We study the properties of the Malthusian Toner-Tu theory in its near ordering phase. Because of the birth/death process, characteristic of this Malthusian model, density fluctuations are partially suppressed. We study this model using the perturbative renormalization group. At one loop we find that the renormalization group flow drives the system in an unstable region, suggesting a fluctuation-induced first-order phase transition.

Despite its simple definition, the Vicsek model [1,2] represents a challenging problem. Its rich phenomenology and complexity arise from the combined effect of the particles' ferromagnetic/alignment interaction and activity, that is particles' self-propulsion. The consequent rewiring of the interaction network implies feedback between density and velocity fluctuations, which turns an equilibrium second-order phase transition into an offequilibrium first-order phase transition [3,4]. The nature of this phase transition can be understood as a liquid-gas phase transition [5] and both numerical simulations and theoretical studies show that density fluctuations play a crucial role in determining the properties of the transitions in active matter [6][7][8][9][10].
A possible way to tackle this problem, and study the transition in the Vicsek model, consists of considering its corresponding hydrodynamic theory: the Toner-Tu (TT) theory [11]. However, precisely because of the velocitydensity coupling, deriving the properties of the TT theory in its near-critical regime is a highly complicated task. In the seminal paper [12], Chen et al. studied the critical dynamics of the TT theory in the incompressible limit (ITT) using the perturbative renormalization group (RG) -one of the first RG calculations showing that active matter systems can belong to a universality class of their own. They found that, if density fluctuations are completely suppressed, the phase transition is secondorder, and a novel active off-equilibrium fixed-point rules the critical dynamics [12,13]. If, on the one hand, this calculation shows the relevance of activity in the critical dynamics, on the other hand, it ignores completely density-velocity coupling; for this reason, it cannot explain why, and how, the order/disorder phase transition becomes first order.
A renormalization group calculation that explains the role of density fluctuations in turning the phase transition from second to first order is still missing. Here, we try to address this problem by studying the critical dynamics of a modification of the TT theory, namely the Malthusian Toner-Tu (MTT) theory [14]. In the MTT theory, density fluctuations are partially suppressed, and in this sense, MTT theory poses in between the ITT and the TT theories. This theory is easy enough to be studied with perturbative renormalization group techniques, but it is still rich enough to produce a behavior qualitatively different from the one of the incompressible theory.
The MTT theory, introduced by J. Toner in [14], is a model for self-propelled agents that reproduce and die. In its hydrodynamic description the equations of motion (EOM) for the velocity field v and the density field ρ, as introduced in [14], are: where P is a pressure term, f is a Gaussian white noise with variance 2D 0 , and (v · ∇)v is the standard material derivative [15]. The term g(ρ) accounts for the birth/death process and it tends to keep density near to a fixed value ρ 0 [14,16]. Near to ρ 0 the function g(ρ) acts as a harmonic restoring force, ∂ t ρ = g (ρ 0 )δρ, where δρ = ρ − ρ 0 , g(ρ 0 ) = 0, and g (ρ 0 ) < 0. At first glance the birth/death mechanism seems to further complicate the problem, but in fact it is a significant simplification. Because of the term g(ρ), density is not a conserved field [14], ∂ t ρ(k = 0, t) = 0, and in this sense ρ is not an hydrodynamic variable anymore [17]. Following [14,18,19], from the continuity equation (2) it is possible to derive a linear relation between density and velocity fluctuations: The force term F (v) = −dV /dv is responsible for the order/disorder phase transition, and it arises from the imitation/ferromagnetic interaction. The state of the system depends on the values of m 0 and J 0 [20]. While the coupling J 0 must be positive -otherwise the potential V would be unbounded -the sign of m 0 determines, at the mean-field level, whether or not the system is in a polarized phase v = 0. Precisely, if m 0 < 0 the average polarization of the system is non-zero, | v | = −m 0 /J 0 .
Before proceeding with the explicit derivation of the EOM, a remark is in order. In general, the TT theory breaks Galilean invariance [15,21], and the Malthusian TT theory is no exception; for this reason, we must include in the model all the possible terms with one derivative and two fields, namely v∇ · v and ∇v 2 [21]. The first of these two terms can be interpreted as a densitydependent alignment force [8]; assuming that the alignment force m 0 v depends on density, and using the relation (3), leads to: where m(ρ 0 ) is simply m 0 , and the second term reproduces the expected structure v∇ · v. The second non-Galilean-invariant term can be interpreted as a non-linear pressure term of the kind P ∼ −γ 3 v 2 . This would lead to the following pressure force, which reproduces the structure of the second non-Galilean-invariant term, ∇v 2 .
To derive the EOM of the model we need an explicit expression for the pressure force ∇P . Following the literature on active matter theories [20], the pressure term P can be expressed as a series in powers of density fluctuations: Using equation (3), the pressure force, at the leading order in the fluctuations, reduces to This leads to the following EOM for the velocity field, We remark that we obtained an EOM for the sole velocity field; however, a relic of the coupling with density is hidden in the parameter γ 2 , whose interpretation is that of a density-dependent ferromagnetic interaction. Equation (8) is the starting point of the perturbative renormalization group calculation.
Before starting to study the properties of the Malthusian theory near the transition, it is helpful to study its linear limit. The linearized EOM, in Fourier (k, t) space, is: It is evident that, even at the liner level, the EOM (10) is anisotropic. A part of the force is in the direction of the field v α , while another part of the force is in the direction of the wave-vector k α . For this reason, we have to distinguish between longitudinal and transverse fluctuations, namely fluctuations parallel and orthogonal to thek direction. We define the transverse P ⊥ αβ (k) and longitudinal P αβ (k) projection operators: We define also the longitudinal Γ 0 and transverse Γ ⊥ 0 kinetic coefficients, which in general are not equal. The responsible for their difference is the pressure force, indeed (Γ 0 − Γ ⊥ 0 ) is proportional to σ 1 . Intuitively, since longitudinal fluctuations are coupled to density fluctuations, as shown in equation (3), they relax faster than transverse fluctuations and Γ 0 > Γ ⊥ 0 . From the equation of motion (10) it is possible to derive the linear response and correlation where the functions G ,⊥ 0 (k, ω) and C ,⊥ 0 (k, ω) are defined as follows These expressions, (15) and (16), can be obtained with standard procedures [22], and a detailed derivation can be found in the SM. The ITT theory [12] can be formally recovered in the Γ 0 Γ ⊥ 0 limit; in this regime, both G 0 and C 0 are negligible, and the bare correlation and response functions are proportional to P ⊥ αβ (k), as in incompressible theories [12,23].
Renromalization group equations -The momentumshell RG approach [22] unfolds through two stages: i) integrating out small length-scale details, on the so-called The effect of step ii is to rescale each parameter by its corresponding naive dimension. The effect of step i is to provide each parameter with a perturbative scaling dimension, which is computed using the Feynman diagrams technique (the technical details can be found in the SM). Combining the perturbative and the naive scaling dimensions leads to the following set of recursion relations, which determine how the parameters of the model change when we observe the system at larger and larger length scales [22]: Here, δX represents the perturbative correction to the parameter X, and η is the perturbative anomalous dimension of the field. A set of 15 Feynman diagrams determines the perturbative contributions -their explicit expressions can be found in the SM. Iterating the RG transformation results in a flow in the space of theories: the renormalization group flow.
It is particularly useful to define a set of effective parameters and coupling constants, where r 0 and u 0 are the effective couplings relative to the force term m 0 v + J 0 v 2 v, and g i,0 are the effective couplings relative to the three non-linearities, v · ∇v, v∇ · v and ∇v 2 . The parameter µ is the ratio between the longitudinal and transverse kinetic coefficients, and it determines how the longitudinal fluctuations relax faster than the transverse fluctuations. The renormalization group flow of these parameters determines the large scale properties of the system. The recursion relations for the effective parameters can be derived by the relations (17)- (21): The naive scaling dimensions of the effective parameters are: where is the distance from the upper critical dimension = 4 − d. This means that in dimension d < 4 both u and the constants g i are relevant in the renormalization group sense. Incompressible case -In the Γ Γ ⊥ limit, the bare correlation functions of the Malthusian TT theory coincide with the ones of an incompressible theory [12,23]. Therefore, in this limit it is reasonable to expect that the transition is ruled by the exponents of incompressible active matter found in [12]. For Γ Γ ⊥ , which implies µ → ∞, the recursion relations (15)- (16) reduce to the ones of [12] (see SM). The large scale behavior for µ → ∞ is thus ruled by the critical exponents, which coincide with the ones of incompressible active matter [12]. However, for large but finite values of µ, namely in the neighborhood of the incompressible fixed point, the recursion relation of µ is revealing that the incompressible fixed point is unstable under a small deviation from incompressibility. This means that, for large but finite µ, the RG flow escapes this fixed point. Moreover, the larger µ, the longer the RG flow lingers near the incompressible active matter fixed point; this structure gives rise to a crossover [22].
To determine whether or not the critical dynamics is ruled by the incompressible fixed point we must consider the stopping condition of the RG flow. The RG flow stops when the correlation length, which trivially scales as ξ b = ξ/b, is of the order of the inverse cutoff Λ −1 , namely when ξ stop = Λ −1 , where l stop = log b (Λξ) is the number of RG iterations before stopping. If at the end of the RG flow µ stop is still large, the incompressible active matter fixed point rules the critical dynamics. For small enough correlation length, ξ µ φ 0 , the critical dynamics is ruled by the incompressible active matter fixed point (28) -where the crossover exponent [22] is In finite systems the correlation length is always bounded by the system size L. Hence, if L µ φ the critical dynamics is ruled by the incompressible fixed point, meaning that the results of Chen et al. [12] also hold in sufficiently small systems with mild density fluctuations µ 1. Conversely, if the system size is large enough, L µ φ , its properties at the transition point are not described by the incompressible active matter fixed point.
In that case we must study the recursion relations for finite µ. We remark that µ can be measured by looking at the system's transverse and longitudinal correlation functions, and, at least in principle, it could be measured both in numerical simulations and real experiments.
Fluctuation-induced first-order phase transition -The properties of the Malthusian TT theory in the vicinity of the transition are determined by the recursion relation (24)- (26), see SM. These recursion relations are complicated, and it is not possible to study them analytically. For this reason, we choose a set of initial values of the parameters and simulate the recursion relations (24)- (26) numerically. Usually, this procedure leads the system's parameters to an infrared-stable fixed point; however, things are more complicated in this model. In figure 1 we show the renormalization group flow obtained by numerical integration; we chose a large initial value of µ, meaning that the RG flow starts in a neighborhood of the incompressible active matter fixed point. As soon as the flow escapes this fixed point -which means that the longitudinal fluctuations, along with density fluctuations, become relevant -the RG flow enters an unstable region, where the ferromagnetic coupling u b becomes negative. Moreover, after that, the RG flow is characterized by run-away trajectories, and the couplings blow up to larger and larger values. The condition u > 0 ensures that the pseudo-potential V (v) is bounded, and as shown in figure 1 after the potential enters the u b < 0 region, the potential becomes unbounded.
To have an idea of what is happening, we consider how the pseudo-potential V (v) changes along the RG flow. At the beginning of the RG flow, the pseudo-potential V (v) takes a double-well shape, characteristic of λϕ 4 theories [22], guaranteeing that the order parameter is small. However, if the flow enters the u b < 0 region, the pseudopotential becomes unbounded. In this regime, the order parameter is clearly not fluctuating around zero anymore, suggesting a breakdown of the perturbative expansion. The situation may seem puzzling: all the theories along an RG trajectory describe the same system, which means that if the system is stable at the beginning of the flow, it must remain stable along the whole RG trajectory, but figure 1 clearly shows the contrary.
We can solve this apparent contradiction considering that the RG transformation may generate additional couplings that grant the stability of the theory also in the u b < 0 region. In figure 2 (right) we show two Feynman diagrams that generate a v 4 v term in the equation of motion, which may stabilize the system in the u b < 0 region. Even if the coefficient J b becomes negative [24], a force term −Kv 4 v could still guarantee the stability of the equation of motion, provided that K > 0. Adding this novel force term in the equation of motion is equivalent to modifying the pseudo-potential V new (v) as follows: If K 0 is zero, the potential V (v) develops a non-zero minimum only for m 0 < 0; moreover, these minima are arbitrarily close to zero provided that m 0 is small enough. This scenario gives rise to the phenomenology of 1; as soon as the flow escapes from the incompressible regime, the ferromagnetic coupling u b (green curve) becomes negative, and then the flow is characterized by run-away trajectories. In the u b < 0 regime the system is unstable. In the inserts we show qualitatively how the pseudo-potential V (v) changes before (blue dot) and after (red square) the flow enters the unstable region; when u b becomes negative the pseudo-potential is unbounded.
a second-order phase transition. Conversely, if K 0 > 0, the potential V may develop a minimum which is not close to zero when J 0 becomes negative. As shown in figure 2, when the constant J 0 becomes negative the potential may develop a global minimum far from v = 0, provided that both m 0 and K 0 remain positive; this would result in a discontinuous, first-order, phase transition.
In the literature, this scenario, where one or more couplings become negative and then run away, is often, but not always [25,26], related to a fluctuation-induced first-order transition. The Heisenberg model with cubic anisotropy and the scalar electrodynamics [27] are two paradigmatic examples in which run-away trajectories result in a first-order phase transition. Non-perturbative techniques are often needed to prove that a run-away RG trajectory results in a first-order phase transition [27,28], but this is beyond the purposes of this perturbative RG calculation.
The stopping condition of the RG equations is determined by the correlation length ξ, or by the system size L. Depending on ξ, the RG flow may enter or not the unstable region. If the correlation length, or the system size, is small the RG flow does not enter the unstable region, and the system displays the phenomenology of a secondorder phase transition, ruled by the incompressible active matter fixed point. Conversely, if the correlation length is large enough, the RG flow enters the u b < 0 region; when it happens because of the strong fluctuations, a phase transition, that should be second-order, becomes first-order [27]. In this sense, we call this phase transition a fluctuation-induced first-order phase transition [27].
Conclusions -We computed the renormalization group flow for the Malthusian Toner-Tu theory. In the incompressible limit this calculation reproduces the re-  sults found in incompressible active matter [12]. Thus, in agreement with the numerical evidences [13], the critical exponents found in [12] hold also for systems with mild density fluctuations, provided that the correlation length (or the system size) is not too large. Furthermore, we found a relation between two measurable parameters, µ and the correlation length ξ (or the size of the system L) that determines whether or not the critical dynamics is ruled by the incompressible active matter fixed point. If the system size is large enough, the RG flow escapes the incompressible active matter fixed point, and enters an unstable region. This suggests that, for large enough systems size L and correlation length ξ, a fluctuationinduced first-order phase transition may occur; corroborating the idea that fluctuations and renormalization of parameters play a crucial role in determining the order of the phase transition [8]. We believe that even though, in MTT theory, density fluctuations are partially suppressed, they are still strong enough to turn the secondorder phase transition into a first-order phase transition. Further studies are needed to clarify the situation; two strategies are including other couplings in the equation of motion or using non-perturbative techniques. This calculation provides one of the first RG evidence that the phase transition in active matter is first-order for large enough systems.