Uncovering novel phase transitions in dense dry polar active fluids using a lattice Boltzmann method

The dynamics of dry active matter have implications for a diverse collection of biological phenomena spanning a range of length and time scales, such as animal flocking, cell tissue dynamics, and swarming of inserts and bacteria. Uniting these systems are a common set of symmetries and conservation laws, defining dry active fluids as a class of physical system. Many interesting behaviors have been observed at high densities, which remain difficult to simulate due to the computational demand. Here, we devise a new method to study dry active fluids in a dense regime using a simple modification of the lattice Boltzmann method. We apply our method to an active model with contact inhibition of locomotion, which has relevance to collective cell migration, and uncover multiple novel phase transitions: two first-order and one potentially critical. We further support our simulation results with an analytical treatment of the hydrodynamic equations.

The dynamics of dry active matter have implications for a diverse collection of biological phenomena spanning a range of length and time scales, such as animal flocking, cell tissue dynamics, and swarming of inserts and bacteria. Uniting these systems are a common set of symmetries and conservation laws, defining dry active fluids as a class of physical system. Many interesting behaviors have been observed at high densities, which remain difficult to simulate due to the computational demand. Here, we devise a new method to study dry active fluids in a dense regime using a simple modification of the lattice Boltzmann method. We apply our method to an active model with contact inhibition of locomotion, which has relevance to collective cell migration, and uncover multiple novel phase transitions: two first-order and one potentially critical. We further support our simulation results with an analytical treatment of the hydrodynamic equations.

PACS numbers:
Symmetry serves a foundational role in all areas of physics today. By identifying the underlying symmetries of a classical many-body system, one can derive the hydrodynamic equations of motion (EOM) that govern the macroscopic dynamics of that system, and, crucially, any other system that respects the same symmetries [1]. In the case of simple (passive) fluids, temporal, rotational, translational, chiral, and Galilean symmetries lead to the Navier-Stokes equations in the hydrodynamic limit [2]. Removing the Galilean symmetry and breaking the conservation of momentum lead instead to the Toner-Tu EOM [3][4][5], which generically describe dry polar active fluids [6] -systems typically composed of particles that generate non-momentum conserving propulsion through interaction with a fixed substrate.
Analysis of a hydrodynamic theory can elucidate the universal behavior exhibited by all generic systems respecting the prescribed set of symmetries; conversely, any particular many-body system defined by certain microscopic rules that respect the same set of symmetries can also be used to study the associated universal behavior in the hydrodynamic limit. An example of the former is the use of dynamical renormalization group methods to confirm the universal long-time tail phenomenon [7] in thermal fluids first discovered by simulations [8,9]; and an example of the latter is the use of lattice gas cellular automata to study the Navier-Stokes equations [10]. Indeed, it came as a surprise in the 1980s that a simple lattice gas model, which resembles nothing like a real fluid microscopically, exhibits hydrodynamic behavior identical to the Navier-Stokes equations. In dry active fluids, similar lattice gas models have helped to clarify the nature of the order-disorder transition in polar active fluids [11], and the critical behavior of motility-induced phase separation [12]. Superseding the lattice gas cellular automata is the celebrated lattice Boltzmann method [13][14][15][16][17], which led to a drastic improvement in the computational efficiency of the simulation requirement.
With their suitability to the almost incompressible regime, lattice Boltzmann methods (LBM) have unsurprisingly been employed to simulate the fluid part of suspensions of liquid crystals and active swimmers [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. These examples all describe wet active matter with momentum conserving forces [6], however the development of the LBM for dry polar active fluids is still missing, which is what we accomplish here. In doing so we are able to study effectively the dynamics of condensed active matter, found in many biological systems. In particular, we apply our simulation method to explore the rich phase behavior of a two dimensional active fluid model with contact inhibition -a model that is of direct relevance to migrating cell layers [35,36], and uncover two novel first-order and a potentially novel critical transition in our system.
Dry polar active fluids. Our system involves the evolution of a distribution function f i (t, r), i ∈ {0, 1, ..., 6}, defined on a two-dimensional triangular lattice (D2Q7 in standard LBM notation [37]) with corresponding lattice vectors e i given by e 0 = 0, and e i = cos[(i − 1)π/6]x + sin[(i − 1)π/6]ŷ for i > 0. The distribution function f i is related to the local mass density and local velocity of the system as follows: where c is the ratio of the lattice spacing ∆x to the time increment ∆t. Typically, scales are chosen so that c = 1, and |u| c is required for stable solutions [38]. We then evolve the discretized Boltzmann equation using the Bhatnagar-Gross-Krook collision operator [39], albeit with fluctuations included, as follows: where τ is a relaxation parameter, f ST i is the steady-state distribution, and whereη i are temporally and spatially uncorrelated random variables uniformly distributed in the interval [−σ, σ]. This particular choice of noise alters the local velocity, but Eq. (3) ensures the density is conserved. Note that the fluctuations do not conserve momentum as there is no reason to do so when considering dry active fluids.
In standard LBM, f ST i is constructed so as to preserve the moments ρ, u at each site, accomplished in a D2Q7 model by with u * = u, where w i are lattice specific weights: w 0 = 1/2 and w i =0 = 1/12. However, since the conservation of local momentum no longer holds in a dry active fluid, the restriction on the form of the steady-state distribution f ST i is drastically relaxed. Here, we will make a simple modification to f ST i , to enable collective motion in our model. Specifically, we maintain the form of f ST i given in Eq. (14), however we take u * = U 0 (ρ)u/|u|, where U 0 is a predetermined local mesoscopic speed, which may be density-dependent. This update rule conserves mass density, but does not conserve momentum nor energy, thus rendering the system in the realm of the hydrodynamic theory of Toner and Tu, rather than Navier-Stokes [38]. Indeed, we will see that such a minimal modification of the update rule is sufficient to elucidate the known phase behaviors of dry active fluids, as well as enabling us to explore novel phase transitions.
For the remainder of this letter we will take U 0 to be of the form where A, B are positive constants. We expect this to reproduce the dynamics of active polar matter because it mimics the key behavior that local particles' movements tend to align. In addition, we assume here that collective motility can be impeded when the density gets too high or even stop altogether, which can be viewed as a form of contact inhibition of locomotion, as observed in cell monolayers [40][41][42] and traffic jams [43,44]. We simulate our active fluid model in a rectangular geometry with periodic boundary conditions ( Fig. 2(d)). This ensures that any ensuing collective motion aligns itself along the x-axis, which helps to reduce noise in our data collection. We use the average speed of the whole system as an order parameter ϕ, given by Necessarily 0 ≤ ϕ ≤ B. Fig. 2 shows how ϕ behaves with varying average densityρ. Two curves have been traced by slowly increasing (blue) and decreasing (red) the total density. The two curves coincide with the exception of four regions, in which there exist hysteresis loops indicating the existence of meta-stable states. We will use these hysteresis loops as proxies for the locations of first-order phase transitions that separate distinct phases. At low densities noise dominates, and the system is in the homogeneous and disordered (HD) phase. As the density increases, alignment effects take over and dense locally aligned bands spontaneously develop, which propagate through a low density disordered region. We refer to this as the banding (B) regime, which corresponds to the coexistence of the HD phase in the dilute regions and a homogeneous, ordered (HO) phase in the dense regions. Increasing the density further causes the dilute HD phase to vanish, resulting in a single HO phase. The first HD-B transition is easily recognised from previous flocking studies [11,[45][46][47] which found it to be first order. The second B-HO transition is less studied, likely because higher density molecular dynamics simulations are more difficult to compute, but has recently been shown to also be a first order transition [47]. Our results are in perfect agreement with both of these as evidenced by the aforementioned hysteresis loops. This lends credence to the conclusion that the two remaining hysteresis loops at higher densities also indicate firstorder phase transitions, which to our knowledge have not been reported previously. These transitions separate the HO region from another HD region, again through a coexistence regime similar to the B regime ( Fig. 2(d)), but with the properties inverted -the high density regions  [38]. The plots show a dot for each of our 180 × 60 triangular lattice sites, with the color denoting the direction of u and the size of the dot depicting the density ρ. White space implies a region of low density. Color wheel is the same as in Fig. 1. Colored symbols are identifiers for the data points in Fig. 3(a). The parameters are: c = 1, τ = 1, and σ = 1/70. See [38] for details on simulation procedure.
are disordered, whilst the low density regions are aligned and moving. In the steady state, as material flows into a dense band it comes to a halt, extending the band on that side, whilst on the other side the band decays as the material flows away. The net result is that the dense band, although locally stationary, appears to propagate against the direction of flow, hence we call it reverse banding (RB). In other words, the dense band moves due to condensation on one side and evaporation on the other. In active matter, this phenomenon was first reported in a collective cell migration model [44], and likened to density waves observed in traffic flow [43]. Interestingly, the same mechanism also underlies the directed motion of condensed drops under the influence of a chemical gradient [48], as well as the formation of a drop lattice in a chemical reaction-controlled phase separated drop system [49].
Phase behavior. To further elucidate the phase behavior of the system, we now vary the contact inhibition parameter A as well asρ ( Fig. 3(a)). When A = 0, we recover the expected phenomenology of a typical polar active model, such as the Vicsek model [50], in that the system transitions from the HD phase to the B regime followed by the HO phase as density increases. However, as A grows, we see that the HO phase will transition back to the HD phase via the RB regime as density increases further. Most interestingly, as A increases, all four firstorder phase transition lines converge to a single point, indicated by a star in Fig. 3(a). If the phase boundary is now traversed by varying A at this fine-tuned density, the system will go directly from the HD phase to the HO phase, without any spatial heterogeneity as displayed in the B or the RB regimes. The order parameter ϕ also rises from zero in the HD phase to a positive value in the HO phase in a continuous manner. We therefore conclude that such a transition is critical and the star thus indicates a critical point of this active fluid model. A snapshot of this critical behavior is shown in Fig. 3(b). We will discuss this critical point further in our hydrodynamic analysis. Finally, we note that in addition to the four regimes already discussed we recognise a fifth somewhat trivial phase of no motion, where |u| = 0 since (a) Phase diagram comparing average densityρ with coefficient A (5). Each sample point was identified as belonging to one of five regimes, categorised by the coloured symbols (see Fig. 2b)) using a method detailed in [38]. The red asterisk region has zero motion because |u| = B −Aρ 2 ≤ 0 with equality depicted by the red line. Black lines depict approximate locations of phase transitions (drawn by eye) which meet at a critical point (star). (b) Snapshot of behavior close to the critical point atρc ≈ 0.195 and Ac ≈ 2.55, using same display scheme as in Fig. 2(b). Simulation details are in [38].
Hydrodynamic description of the transitions. In the hydrodynamic limit, our model is necessarily described by the Toner-Tu equations due to our imposed symmetries, and we confirm this analytically using the Chapman-Enskog expansion method in [38]. We now describe our simulation findings analytically by performing a linear stability analysis on the following reduced Toner-Tu EOM as they suffice for our discussion: where g is the momentum density and g ≡ |g|.
Since the bands in the B and RB regimes are both perpendicular to the direction of collective motion in the dense and dilute regions, respectively, the most unstable wave vectors in the HO state are expected to be parallel to the direction of motion. Therefore, we need only consider a 1D system, where the coordinate x is the direction of the collective motion. Expanding about a homogeneous collective motion state with , we have to linear order: Note that for a stable system β > 0 necessarily, thus α 0 > 0 implies collective motion, with α 0 ≤ 0 giving a stationary state. Solving for s and focusing on the real part, which determines the stability of the system, we have The first eigenvalue 2α 0 corresponds to the fast relaxation when the momentum deviates from the mean field value α 0 /β in the absence of spatial variations. The second eignvalue quantifies when the instability sets in.
In particular, if α 0 → 0 and α 1 = 0, the system is always unstable. In the previously mentioned studies, [11,46,47,50], α 1 is typically assumed to be positive, so that alignment interactions increase with density, and as α goes from negative to positive the system will necessarily become unstable, resulting in the banding instability. However, it is also clear that the instability condition does not depend on the sign of α 1 . Hence, if α goes from positive to negative due to contact inhibition, the system will again be unstable. This is exactly what occurs in our simulation and corresponds to the Reverse-Banding (RB) regime in our phase diagram. In our simulation we found that all four first-order phase transition lines converge at a critical point. With our hydrodynamic argument, we can see that this critical point is reached by setting both α 0 and α 1 to zero, which is achieved by fine-tuningρ and A in our simulation [38].
We believe that this critical point is distinct from any critical behavior previously discussed in the context of active matter for the following reasons. We are aware of critical phenomena analyzed, either analytically or computationally, in the following active matter systems: 1) self-propelled particles with long-ranged metric alignment interactions [51,52], 2) incompressible active fluids [53], 3) critical motility-induced phase separation [12,54,55], and 4) self-propelled particles with velocity reversals and alignment interactions [56]. Our critical system is different from system 1) because ours does not have long-range interactions, from 2) because our velocity field is not divergence free, from 3) because critical MIPS happens in the disordered phase while criticality in our system occurs at the onset of the ordered phase, from 4) because our system exhibits collective motion while system 4) exhibits only orientational order, but no collective motion. Besides the above comparison, as our model spontaneously breaks rotational invariance at the critical point, it is also distinct from critical active matter in which the direction of motion is determined by an external field [57,58].

Conclusion & Outlook.
Through symmetry considerations, we have generalized LBM to simulate a dry active fluid system with contact inhibition, and uncovered novel phase transitions in the dense limits of the system. Our work is of direct relevance to the fields of active matter, cell biophysics, and non-equilibrium phase transitions. Looking ahead, we believe that the drastic improvement on computational efficiency will help solve interesting open questions in dry active fluids, such as the dynamic scaling in the incompressible limit [59], and universal behavior in chiral active fluids [60][61][62]. Furthermore, the fact that LBM are adaptable to complex geometries renders them suitable for investigating interfacial instabilities in a multiphase active fluid system, such as in the context of tissue regeneration [63][64][65]. Finally, we believe our method can be straightforwardly extended to higher dimensions. D.N. was supported by the EPSRC Centre for Doctoral Training in Fluid Dynamics Across Scales (grant EP/L016230/1).

SUPPLEMENTAL MATERIAL CONVENTIONAL LATTICE BOLTZMANN METHODS
The lattice Boltzmann methods (LBM) refer to solving the Boltzmann equation on a lattice. The Boltzmann equation itself concerns the temporal evolution of a distribution function f (t, r, ξ), which corresponds to the mass density of matter moving with velocity ξ at location r and time t.
The typical variables of interest are where D is the spatial dimension, ρ is the mass density field, and u is the velocity field. Discretizing ξ into a finite number of velocity vectors directed towards neighbouring spatial nodes naturally constructs a lattice in the domain. There are many suitable choices of lattice, but they must contain the necessary number of symmetries to preserve enough higher moments of the system in question according to the desired level of accuracy. [1,2]. The DQ notation is conventional for describing the lattice structure of the domain, where D is again the dimension, and Q is the number of nodes each is connected to, including itself. Typically a 'zero' vector is included in the discrete directional space to help balance higher moments. Here, we will focus on a D2Q7 lattice, which is a triangular lattice in two dimensions with a discretized distribution function f i (t, r), where i = 0, 1, ..., 6. The corresponding vectors are e 0 = 0, and e i = cos[(i − 1)π/6]x + sin[(i − 1)π/6]ŷ for i > 0. The local mass density and velocity are thus given by where c is the ratio of the lattice spacing ∆x to the time increment ∆t. For accurate solutions velocities must obey |u| c, which is known as the low Mach number assumption. Typically scales in the LBM scheme are chosen such that c = 1.
The discretized Boltzmann equation can then be evolved using the Bhatnagar-Gross-Krook collision operator [3] where f eq i (t, r) is a lattice dependent steady-state (equilibrium) distribution and τ is a relaxation parameter. When simulating the Navier Stokes equations, τ is linearly proportional to the viscosity. The evolution of this equation is typically split into two processes-the streaming step and the collision step, the latter of which is effectively a local redistribution. The design of the steady-state distribution is highly flexible and the only requirement is that it respects the symmetries and conservation laws imposed on the system. For the Navier-Stokes equation, f eq i is obtained by discretizing the Boltzmann distribution for velocities at thermal equilibrium [4] and in a D2Q7 framework is given by where w i are lattice specific weights: w 0 = 1/2 and w i =0 = 1/12. This discretization is possible through the low Mach number assumption |u| c.

SUMMARY OF MODIFIED LATTICE BOLTZMANN MODEL ALGORITHM
A single time step is partitioned into the following steps

NEGATIVE CORRECTION
In our modified lattice Boltzmann method, we have included stochastic noise as detailed in the main text. It is possible for one or more of the f i values to become negative in this process. This is not permitted, hence we correct this using the following procedure at each lattice site.
2. If f i < 0 for i > 0, then take j, such that e j is the reverse direction of e i ( 1 ←→ 4, 2 ←→ 5, 3 ←→ 6), and add |f i | to f j , then set f i = 0.
3. By now f i ≥ 0 for all i, however it is possible that the density has changed. This is corrected by rescaling Also note that due to both the noise and the rescaling of u in f ST i , it is possible that when τ < 1, even if f i > 0 for all i. Therefore to ensure that no values become negative, we only use τ ≥ 1. Different values of τ were tested, yielding results that were qualitatively similar to those described in the letter.

General notes
For all of the results in the letter we used a 180 × 60 node triangular lattice with periodic boundary conditions. Every simulation was initialized by choosing 2 uniformly distributed random numbers for each lattice site -the first deciding the initial direction for u with each direction in [0, 2π] equally likely, and the other choosing the initial density such that it falls within ±10% ofρ, the prescribed average density. The initial magnitude of u is determined by Eq. (5) in the MT, and the initial distribution given by f ST i (14). The system was then evolved 10000 steps before a "run" would formally begin. This marks the end of the initialization process.
As detailed below, during the course of each run the mean density would be deliberately changed by a small amount δρ. Each time this occurred, this was achieved by adding δρ/7 to every f i across the lattice -the lattice was not reinitialized. For measurements involving averaging over multiple runs the initialization process was repeated between each run.
With the exception of A andρ, the parameters were held constant with ∆x = 1, ∆t = 1, c = 1, B = 0.3, τ = 1, and σ = 1/70. Simulation procedure in Fig. 2 For Fig. 2(a) A = 2, and onlyρ is varied. The goal was to detect any hysteresis effect when adiabatically increasing or decreasing the density. For the increasing density run, we start withρ = 0.1 and increase it by |δρ| = 0.0025 every t sample steps. A sample was taken just before each density change. For higher densities the system would reach equilibrium faster, hence in order to be able to detect the hysteresis effect between HO-RB and RB-HD, a lower value of t sample = 500 was used aboveρ = 0.19, with t sample = 1000 otherwise. This is because in a finite system, running the system long enough would eliminate all hysteresis effects.
For the decreasing density run, we start withρ = 0.34 and again decrease it by |δρ| = 0.0025 every t sample steps, which are again 500 aboveρ = 0.19 and 1000 otherwise. The system was reinitialized between increasing and decreasing runs. As a result, they should be treated as independent.
The whole increasing and decreasing density cycle was repeated over 90 runs. The order parameter ϕ was evaluated for each sample, and the values for samples at equal densitiesρ taken in the same direction of the cycle were collated and averaged. The error bars are the standard deviation ofφ, that is, the estimate of the standard deviation of the samples σ ϕ scaled by the square root of the sample size, σ ϕ / √ 90. Simulation procedure in Fig. 3 For Fig. 3(a) A was varied as well asρ. Results were obtained by fixing A then slowly varyingρ, much as for Fig. 2(a) except: there were no descending runs, ∆ρ = 0.01, and t sample = 3000, which is enough for the system to reach its steady state. Each run was repeated 10 times. The 10 samples for each point (ρ, A) were used to calculate the parameters ϕ, ∆ρ, ζ (see below) which were then averaged.

REPRESENTATION OF STATIC CONFIGURATIONS
Each of the plots in Fig. 2(b) and Fig. 3(b) are scatter plots of 180 × 60 circles, arranged in the triangular lattice configuration. They represent the density and velocity fields for the final time step before the density was changed. In Fig. 2 The colour of each dot corresponds to the direction of the macroscopic velocity u measured at that site, according to the colour wheel given in Fig. 4. As our system is two dimensional a single angle (color) is sufficient to determine the direction. The area of each dot is assigned by the following formulae where ρ max , (ρ min ) is the maximum (minimum) value of ρ over all of the lattice sites in that frame, a 0 is an appropriately chosen value for the size of the biggest circle, which depends on the spacing of the particular plot. is a small value to ensure no values are exactly zero.

MOVIES
Five movies are provided to demonstrate the behaviour of the system as shown in Fig. 2(d) and Fig. 3(b), using the same visual representation method. For each of the first four movies the system was initialized with A = 2 andρ again given by A frame was captured every 5 time steps. Each movie is constructed of 200 frames displayed at 12 frames per second.

REGIME CATEGORISATION IN THE PHASE DIAGRAM
For the phase diagram in Fig. 3(a) we categorize the system's behavior using the following systematic method. Firstly, anything above the red line B − Aρ 2 = 0 was marked as motionless.
Since we have a finite system, even in the most disordered phase ϕ is never exactly zero, and as it is strictly positive it cannot average to zero either. Therefore we had to impose a cutoff to distinguish between a completely disordered state and one with an ordered phase (including the B and RB phases). The cutoff ϕ threshold is 0.025, with ϕ < ϕ threshold marked as HD.
The other key statistic we used to classify the behaviour was large scale density variation. Due to the aspect ratio of the system and the periodic boundary conditions, any travelling wave fronts in the system will naturally tend to align themselves with the y-axis, hence we averaged the density over the y-direction in order to reduce the effect of the stochastic fluctuations. We then quantified the density variation by a single parameter ∆ρ, calculated by integrating over the power spectrum of the density deviation from the mean in the x-direction. Explicitly where ρ y (x) = 1

Ly
Ly 0 dyρ(x, y) and F is the Fourier transform. As with ϕ we had to impose a cutoff to distinguish between phase separated regions and local noise. We chose this to be ∆ρ threshold = 12, with ∆ρ > ∆ρ threshold considered to be either B or RB. Fig. 5 shows a 3D representation of ∆ρ. The two high ridges clearly correspond to the B and RB regions. These narrow in thickness and the peak decreases as they approach each other and seem to vanish at the critical point. The two aforementioned quantities, ϕ and ∆ρ, do not enable us to distinguish between the B and RB phases. Therefore we also measured the following parameter: where ζ is large if high density regions are stationary, as in RB. Hence we identified RB values with ζ > ζ threshold = 0.1. Visual inspection of the samples at several distinct points in the phase diagram shows the categorization to be largely accurate.

LATTICE BOLTZMANN TO TONER-TU EQUATIONS VIA CHAPMAN-ENSKOG EXPANSION
In this section we will demonstrate that our lattice Boltzmann model simulates a version of the Toner-Tu equations. By expressing our active lattice Boltzmann model in the form of a conventional lattice Boltzmann scheme with a body force term, we can use the results of [5] to derive the macroscopic equations.
Let's introduce the notation Hence in the conventional LBM f eq i = ξ i (ρ, u) and in our active model f ST i = ξ i (ρ, u * ). The evolution equation becomes where Note that As in [5], we define the 'actual velocity' according to After expressing Eq. 21 in terms of the actual velocity, we can now use the results of [5] to determine our macroscopic equations. They are the continuity equation and the momentum equation whereΠ = ρν ∇v + (∇v) T .
Note, that in deriving this several O(u 3 ) terms have already been ignored, consistent with a small Mach number assumption. The final term is given by τ i e i e iFi = ρu * u * − ρvv .
Recall that u * = U 0 u |u| . We can express both u and u * in terms of v and U 0 by and Thus This is the key term of interest. It corresponds to a linearised version of α(ρ)g − βg 2 g in the Toner-Tu equations quoted in the main text (g = ρv), where the linearisation has occurred about the stationary points |g| = α/β. Matching terms gives

FLUCTUATION-INDUCED RENORMALIZATION OF α
The above analysis did not incorporate noise into the picture, we will do so now and demonstrate how, when the noise is included, the renormalized coefficient α will naturally depend on the density ρ as well. Here, we aim mainly to connect our hydrodynamic EOM derived from our Chapman-Enskog expansion to the linear stability analysis in the main text, and so will not perform a full renormalization procedure on the EOM. Specifically, we will show how the β term, together with the noise term with strength σ, renormalize α around the critical region when α ≈ 0.
The calculation mirrors almost exactly the corresponding calculation in the study of incompressible active polar fluids at criticality [6], except that the projection operator is absent here since our system is not incompressible. Specifically, the shift in α is given by where S D ≡ 2π D/2 /Γ(D/2) is the surface area of a unit sphere in D dimensions, Λ −1 corresponds to the small length scale cutoff, and Λd corresponds to the width of the momentum shell being integrated over in the above coarse-graining process.
Adding this shift to α modifies it to, to O(α): where we have used the expression for β in 34.
The above expression indicates that the critical transition when U 0 and ρ are fine tuned such that α R = 0 and ∂α R /∂ρ = 0, which can be achieved by fine tuningρ and A in the expression of U 0 . These are given bȳ In our simulations we foundρ c ≈ 0.195, A c ≈ 2.55. These values give, B ≈ 0.29, C 1 ≈ 0.0014, which in particular is consistent with our actual value of B = 0.3.