Novel physics arising from phase transitions in biology

Phase transitions, such as the freezing of water and the magnetisation of a ferromagnet upon lowering the ambient temperature, are familiar physical phenomena. Interestingly, such a collective change of behaviour at a phase transition is also of importance to living systems. From cytoplasmic organisation inside a cell to the collective migration of cell tissue during organismal development and wound healing, phase transitions have emerged as key mechanisms underlying many crucial biological processes. However, a living system is fundamentally different from a thermal system, with driven chemical reactions (e.g., metabolism) and motility being two hallmarks of its nonequilibrium nature. In this review, we will discuss how driven chemical reactions can arrest universal coarsening kinetics expected from thermal phase separation, and how motility leads to the emergence of a novel universality class when the rotational symmetry is spontaneously broken in an incompressible fluid.


Introduction
Collective phenomena are intimately linked to the phenomenon of phase transitions in physics. At a typical phase transition, a many-body system with constituents that interact only locally with their neighbours, be they molecules or living organisms, can collectively change their behaviour upon a subtle change of a single parameter, to the extent that the qualitative behaviour of the whole system is modified. Phase transitions encompass many everyday phenomena such as oil drop formation in a salad dressing and magnetisation in some metals. The study of phase transitions is of fundamental interest to physicists because of the emergence of universal behaviours at a phase transition. By universal behaviour, we mean that certain properties of the system are highly independent of the system's microscopic details. In the salad dressing example, such property can be the power law exponent that governs how the average size of oil drops changes with time; in the example of magnetisation, it can be the power law exponent that governs how the correlation function of two atomic spins decays with respect to their distance. Recently, phase transitions in living systems have also been under intense attention. Indeed, the generic non-equilibrium nature of biological systems have given rise to novel universal behaviours not seen before. In this review, we will focus on two such examples: phase separation with driven chemical reactions, motivated by the mechanism underlying the formation of some membraneless organelles in cells [1,2], and spontaneous symmetry breaking in incompressible active matter, motivated by its relevance to biological tissues [3,4,5] (Fig. 1).
In Sect. 2, we will first describe the relevance of phase separation in cytoplasmic organisation and then review the latest findings on how driven chemical reactions (e.g., adenosine triphosphate (ATP)-driven phosphorylation) can lead to co-existing phase-separated protein drops in the cytoplasm, contrary to the universal coarsening behaviour expected in its equilibrium counter part. In Sect. 3, motivated by the collective behaviour found in motile organisms, we will introduce a generic model of incompressible active fluids from a symmetry consideration. We will then elucidate how a novel critical behaviour emerges at the onset of collective motion, and discuss the universal behaviour of a two dimensional incompressible active fluid in the ordered phase. Finally, we will end with Conclusion & Outlook.
2. Non-equilibrium phase separation: a mechanism for cytoplasmic organisation

Membrane-less organelles
Biological cells organise their contents in distinct compartments called organelles, typically enclosed by a lipid membrane that forms a physical barrier and controls molecular exchanges with the surrounding cytosol. Recently an intriguing class of organelles lacking a membrane is being studied intensely [8]. Membrane-less organelles have attracted an intense interest from the biology community as they are present in many organisms from yeast to mammal cells, and are critical for multiple biological functions. For example P granules are involved in the asymmetric division of the Caenorhabditis elegans embryo [9], and stress granules assemble during environmental stress and protect cytoplasmic RNA from degradation [10] (Fig. 1 a)). Membraneless organelles are generally spherical, fuse together upon contact [11,12], and their components quickly shuttle in and out [13,14], thus resembling liquid drops. Indeed, Figure 1. Cytoplasmic phase separation and tissue dynamics as active matter. a) In many distinct types of cells, certain proteins can phase separate from the cytosol to assemble membrane-less organelles, such as the stress granules (yellow drops) shown here in human epithelial cells (HeLa) [6], akin to oil drop formation in an oil-water mixture (b). c) In a monolayer of Madin-Darby Canine Kidney (MDCK) cells, the cells in the tissue can undergo dynamical rearrangement as shown by the snapshot of the velocity field shown in (d) [7]. strong experimental evidence indicates that membrane-less organelles are assembled via liquid-liquid phase separation [2,15,16], a common phenomenon in every day life responsible for example for oil drop formation in water ( Fig. 1 b). Under the equilibrium condition phase separation is well understood [17]. However cells are driven away from equilibrium by multiple energy-consuming processes such as ATPdriven protein phosphorylation [18], which can potentially affect the phase-separating behavior of membrane-less constituents. For example P granules do not distribute homogeneously in the cytoplasm but preferentially to the posterior side of the cell [19], and stress granules form and dissolve according to environmental cues [20]. The fascinating physics associated to membrane-less organelles are only beginning to be investigated [12,21,19,11,22,23,24,25].
In this section, we will start with a brief summary of relevant principles of equilibrium phase separation in Sec. 2.2. We will then review the latest progress on phase separation driven out of equilibrium by energy-driven chemical reactions in Sec. 2.3. Specifically we will focus on a ternary fluid model of the cell cytoplasm where chemical reactions can convert phase-separating molecules into soluble molecules and vice versa. We will show how such reactions can control drops assembly and size, and suppress Ostwald ripening, allowing a collection of organelles to coexist in the Equilibrium phase diagram of a ternary mixture composed of molecules P (red disks), S (blue disks) and C (not shown). Outside the phase boundary (green curve) the system is homogeneous (" " symbol). Inside the phase boundary ("♦" symbol) the system phase separates into two phases "in" and "out" of distinct concentrations. The coexistence concentrationŝ P in,out ,Ŝ in,out are given by the intersections between the tie-lines (straight lines) and the phase boundary. cytoplasm.

Equilibrium phase separation
Interactions between molecules can cause a homogeneous system to undergo a phase separation, i.e. the spontaneous partitioning of a system into multiple phases of distinct properties such as concentration [17]. The transition from the homogeneous state to the phase-separated state is controlled by parameters such as temperature, pressure or concentrations. The set of parameters leading to phase separation are represented in a phase diagram, as shown in Fig. 2 for a ternary mixture composed of molecules P (red disks), S (blue disks) and C (not shown). The molecular concentrations are labelled by the same symbols P, S, C. We assume incompressibility and that all three types of molecules occupy the same volume, so the combined concentration ψ ≡ P + S + C is homogeneous. The concentration C at any point in the phase diagram is therefore given by ψ − P − S. Outside the phase boundary (green curve) the system is homogeneous (" " symbol). Inside the phase boundary ("♦" symbol) the system phase separates into two phases ("in" and "out") of distinct concentrations (P in,out ,Ŝ in,out ), given by the intersections between the tielines (straight lines) and the phase boundary.
At the equilibrium condition a multi-drop system is unstable due to Ostwald ripening that causes large drops to grow and small drops to evaporate [26,27], and coalescence caused by the fusion of drops upon contact [28]. Eventually a unique drop remains in a finite system. Since the crowded environment of the cytoplasm inhibits the diffusion of macromolecular aggregates [29] we will ignore drop coalescence in this review and focus on Ostwald ripening. Ostwald ripening is caused by two ingredients. One is the Gibbs-Thomson relation that relates the coexistence concentration to the drop radius. For example for the P concentration we have: were l c is a capillary length [25] andP in,out are the coexistence concentrations for a flat interface (R → ∞, Fig. 2). The smaller the drop, the larger the concentration outside which is a consequence of the Laplace pressure [17]. The second ingredient driving Ostwald ripening is the existence of a diffusive concentration profile between drops, which can be approximated by an ideal gas diffusion profile in the case of small concentration outside drops [27]: where P out (r, t) is the P profile outside a drop at time t and distance r from the drop centre, and we have assumed spherical symmetry centred on the drop. D is the molecular diffusion coefficient. We use a crucial assumption known as the quasi-static approximation: the dynamics of drop radii is much slower than the equilibration of the concentration profiles. Therefore the profiles can be assumed to be always at steady state, and imposing the Gibbs-Thomson relation at the interface (Eq. (2)) we find: where P out (∞) is the concentration far from drops and ∆ ≡ P out (∞) −P out is commonly referred to as the supersaturation. We have also assumed small drop density so that drops are on average far from each other and the concentration far from drops P out (∞) is homogeneous. In other words, drops interact with each other only via this common far-field. The diffusive profile leads to a flux J out→in = D∇P out | r=R of molecules P entering the drop at the interface [27]: When the flux J out→in is positive molecules P accumulate at the interface leading to drop growth, while the drop shrinks when J out→in < 0. We show in Fig. 3 a) the flux J out→in for varying drop radius R, assuming fixed supersaturation ∆. There exists a steady state radius R * (J out→in = 0, purple disk) that is unstable, called nucleus radius: smaller drops evaporate (J out→in < 0, left arrows) and larger drops grow (J out→in > 0, right arrows). This is due to the Gibbs-Thomson relation (Eq. (2)) dictating that the concentration near a small drop is larger than that near a large drop. As a result a diffusive flux is directed from small to large drops ( Fig. 3 b), red arrows). A multi-drop system is therefore unstable against Ostwald ripening, i.e large Ostwald ripening in equilibrium systems. a) The influx J out→in of molecules P at a drop interface for varying drop radius R and constant supersaturation ∆ (Eq. (5)). The steady state R * (J = 0, purple disk) is unstable: drops larger than R * grow while smaller ones shrink. Insert: as ∆ decreases, the critical radius R * increases. b) Schematic of P concentration profiles along an axis connecting the centres of two drops of different radii R 1 > R 2 . The Gibbs-Thomson relation (Eq. (2)) dictates that the solute concentration is lower close to the large drop (Pout(R 1 )) than close to the small drop (Pout(R 2 )). This causes a diffusive flux from small drops to large drops (red arrows). c) A multi-drop system is therefore unstable against Ostwald ripening: small drops evaporate and large drops grow. As fewer drops survive the supersaturation ∆ decreases causing more drops to dissolve (see insert in a)). Eventually a unique drop remains in a finite system. drops grow at the expense of small ones (Fig. 3 c)). As drops disappear, the average drop radius increases and the P concentration near drops decreases according to the Gibbs-Thomson relation (Eq. (2)). Hence the supersaturation ∆ decreases because the total number of molecules P in the system is fixed. This causes the critical radius R * to increase, as shown in the insert of figure Fig. 3 a). Therefore Ostwald ripening occurs until, in a finite system, a single drop survive (Fig. 3 c)) [27].
In the cell cytoplasm Ostwald ripening is not a desirable feature since it forbids the stability of multiple membrane-less organelles. In Sec. 2.3 we show how nonequilibrium chemical reactions in our ternary mixture can arrest Ostwald ripening.

Phase separation in presence of non-equilibrium chemical reactions
The presence of non-equilibrium chemical reactions have been proposed recently to explain multi-drop stability in the cytoplasm, as well as being a mechanism to control the formation, dissolution and size of membrane-less organelles [22,23,30,31]. We investigate in this section the physical mechanisms involved and recover the results from [30] by using a different and more intuitive approach.
We consider the ternary mixture discussed in Sec. 2.2 and now assume that P and S can be inter-converted by the chemical reactions: where k, h denote the forward and backward reaction rate constants (Fig. 4). In an equilibrium system k, h should depend on the local concentrations. Indeed, since the mixture phase separates, the interactions between molecules P , S or C must be distinct, leading to concentration-dependent reaction rate constants, as detailed in Fig. 5. In our non-equilibrium system however, we assume that energy consuming processes such as ATP-dependent phosphorylations drive the chemical reactions [17] such that the rate constants k and h are concentration-independent. In Fig. 6 we show the results of Monte Carlo simulations of our non-equilibrium ternary mixture on a two-dimensional lattice [30]. Molecules P are represented by red dots, molecules S by blue dots and molecules C are not shown. Molecules P form energetic bounds with neighbouring molecules P to induce phase separation, and P and Non-driven and driven chemical reactions in phase separating systems. We consider here for clarity a particular example from our ternary mixture, where P molecules (red disks) form bonds (black lines) with neighbouring P molecules while molecules S (blue disks) do not interact. We also concentrate on the forward reaction P → S. The chemical conversion of a given P into a S requires an activation energy ∆U to break these bounds. Since drops are enriched in P , many bonds need to be broken and ∆U is high (upper graph). Conversely the cytoplasm is poor in P and ∆U is thereby small (lower graph). At thermal equilibrium the chemical conversions are non-driven and the energy required to overcome the barrier ∆U is provided by thermal fluctuations alone, and therefore the reaction rate constant k decreases exponentially with ∆U [32]. As a result, k is thus concentration-dependent: k is small inside drops where ∆U is high, and large in the cytoplasm. In the case of driven chemical reactions, such as ATPdependent phosphorylation, an external source of energy is provided [18]  S randomly convert into each other with fixed transition probabilities according to Eq. (6). In the initial state all three types of molecules are distributed homogeneously (leftmost snapshot). Phase separation spontaneously occurs, P -rich drops appear and grow (middle snapshot), and the system undergoes Ostwald ripening leading to an increase of the average drop radius. Eventually a steady state is reached where the multi-drop system is stable and drops have similar radii (rightmost snapshot). Therefore the presence of the non-equilibrium chemical reactions arrest Ostwald ripening in our system. Interestingly drops tend to be of almost equal distance from each other, thus forming a close-packing lattice. In the next sections we will provide a theoretical analysis of the observed behaviour.
2.3.1. Introduction of a new length scale (ξ). We assume that our system remains close to equilibrium to the extent that local equilibrium applies. The interfacial Figure 6. Monte Carlo simulation of phase separation in a ternary mixture with non-equilibrium chemical reactions [30]. Molecules P (red dots) form energetic bounds with neighbouring P , while molecules S (blue dots) and C (not shown) do not form bounds. Chemical reactions transform P into S and vice versa with fixed reaction rate constants (Eq. (6)). In the initial state molecules are randomly distributed (leftmost snapshot). As time progresses drops form and undergo Ostwald ripening leading to the increase of the average drops radius (middle snapshot). Eventually a multi-drop steady state is reached (rightmost snapshot). At the steady state drops have roughly the same radius and are evenly distributed. Simulation details can be found in [30]. region therefore remains governed by equilibrium principles and the coexistence concentrations at the interface P in,out (R) are given by the Gibbs-Thomson relations (Eqs (1) and (2)). For simplicity we will assume that S is inert to phase separation in the sense that its concentration inside and outside drops are identical (S in (R) = S out (R)). The more general scenario where S can segregate inside or outside drops is treated in Ref. [30]. Since only P phase separates we will refer to P as the solute. For simplicity we assume the same diffusion coefficient D for both species P and S, and moreover that D and the reaction rate constants k and h are identical both inside and outside drops. In the context of membrane-less organelles this assumption is justified by the fact that drops are not highly packed but porous, and components rapidly shuttle in and out [13,33,14]. The concentrations P and S therefore obey the following reaction-diffusion equations: where we used again the quasi-static assumption (∂ t P = 0, ∂ t S = 0). In a homogeneous system ∇ 2 P = ∇ 2 S = 0 and the concentrations are at their chemical equilibrium: S in,out = P in,out k/h. Adding Eq. (7) and (8) and imposing no-flux boundary conditions in drop centres and at the system boundary [30] we find that the combined concentration P + S is homogeneous inside and outside drops: Non-equilibrium concentration profiles in a multi-drop system at steady state. The profiles of P (red curve) and S (blue curve) (Eq. (10)) are shown along an axis connecting the centres of two drops of radius R. The interfaces are at local thermal equilibrium so the coexistence solute concentrations (P in,out (R)) are given by the Gibbs-Thomson relations (Eqs (1) and (2)). The concentration S is continuous at the interface (S in (R) = Sout(R)). The non-equilibrium chemical reactions coupled with diffusion and phase separation create a circulating flux of molecules P and S between drops and cytoplasm. Drops are rich in P so the chemical conversion P → k S dominates, leading to an accumulation of S molecules inside drops (downward red/blue arrows). The excess of S is then transported by diffusion toward the cytoplasm (blue arrows). In the cytoplasm the reverse reaction dominates, leading to creation and accumulation of P molecules between drops (upward blue/red arrow), which diffuse toward the drops (red arrows).
where φ in,out are independent of r. The reaction-diffusion equations (Eqs (7) and (8)) therefore decouple and we obtain the solute concentration profiles: in,out e r/ξ + U with U (n) in,out independent of r and is a new length scale introduced in the system by the chemical reactions, which is the length scale of the reaction-induced concentration gradients. We show in Fig. 7 a schematic of the concentration profiles of P (red curve) and S (blue curve) in a multi-drop system at steady state, along an axis connecting two drop centres. The chemically-induced concentration gradients and chemical reactions create a circulating flux of molecules P and S from drops to cytoplasm and vice versa. Inside drops, the concentration of P is high and the chemical conversion P → S dominates. Therefore the concentration of S increases inside drops (downward red→blue arrows). Molecules S are then transported outside drops by diffusive fluxes (blue arrows). In the cytoplasm on the contrary, the reverse reaction S → P dominates so the concentration of P increases between drops (upward blue→red arrow). Molecules P are then transported by diffusive fluxes toward drops (red arrows). Molecular fluxes at a drop interface for varying drop radius R and fixed supersaturation ∆. The efflux J in→out (red curve, Eq. (13)) driven by chemical reactions competes with the influx (brown curve) which is identical to that at the equilibrium condition (i.e without chemical reactions, see Eq. (5) and Fig. 3 a)). Drops grow if J out→in > J in→out and shrink otherwise. Two steady states radii R * − < R * + exist. R * − (purple circle), reminiscent from the equilibrium case, is unstable against Ostwald ripening: smaller drops shrink (left arrow) and larger drops grow (right arrows). At large radius the efflux J in→out dominates so drops shrink. There exists therefore a second steady state R * + (green disk) that is stable against Ostwald ripening: larger drops shrink and smaller drops grow. The insert shows the effect of decreasing the supersaturation ∆. The smallest unstable steady state is labelled Rc, the largest stable steady state is labelled Ru, and the critical radius R l (black square) is the boundary between the unstable and the stable region (Eq. (14)).
We will now study how these fluxes can arrest Ostwald ripening, in two limiting regimes based on the drop size relative to the gradient length scale ξ (Eq. (11)).

Small drop regime (R ξ).
We start by considering the regime of small drop radius R compared to the gradient length scale ξ (Eq. (11)). In a homogeneous system (not phase-separated) the concentrations are at their chemical equilibrium (∇ 2 P = ∇ 2 S = 0 in Eqs (7) and (8)). If the initial supersaturation is small, the cytoplasmic concentrations change only slightly during the phase separation and thus remains close to chemical equilibrium. We will therefore neglect chemical reactions in the cytoplasm and assume purely diffusive profiles outside the drops (Eq. (4)). The influx J out→in of P molecules at a drop interface is therefore identical to that in the equilibrium case (i.e. without chemical reactions, Eq. (5)). A quantitative analysis supporting this approximation is given in Appendix A.
Inside drops, since R ξ, the exponential terms in Eq. (10) are close to one. Imposing the no-flux boundary condition in the drop centre and enforcing the Gibbs-Thomson relation at the interface Eq. (1) we get: We find that the solute profile remains flat inside drops, unaffected by the chemical reactions. This can be interpreted as follow. The dominant reaction inside drops, P → k S, creates an excess of S molecules. Since drops are small compared to the gradient length scale ξ (Eq. (11)) the diffusion coefficient D is large enough so that the excess of S is quickly evacuated outside drops by diffusion, leaving the drop concentrations unperturbed. The subsequent efflux J in→out of molecules S at the drop interface is therefore simply proportional to the degradation rate of P molecules (kP in ) multiplied by the drop volume and divided by the drop surface: We show in Fig. 8 the fluxes J in→out (brown curve) and J out→in (red curve) at the drop interface, for varying drop radius R and at fixed supersaturation ∆. The only difference with the equilibrium case (without chemical reactions, Fig. 3 a)) is the existence of the chemical reaction-driven efflux J in→out , which competes with the influx J out→in . At small drop radius J out→in is negligible, and the steady state R * − (purple disk) is reminiscent from the equilibrium case ( Fig. 3): the system is unstable against Ostwald ripening, i.e. smaller drops dissolve (J in→out > J out→in ) while larger drops grow (J in→out < J out→in ). At large radius however the efflux J in→out dominates causing drop shrinkage. This introduces a new steady state R * + (green disk) that is stable against Ostwald ripening: larger drops shrink and smaller ones grow.
As ∆ gets smaller, R * − increases while R * + decreases, as shown in the insert of Fig.  8, and we can graphically identify three critical radii for which we provide approximate expressions [30]: where φ is the average concentration of P + S in the entire system. Due to the incompressibility condition the overall solvent concentration C in the entire system is given by ψ − φ, where we recall that ψ is the combined concentration of all species P + S + C (Sec. 2.2). R c is the smallest unstable steady state R * − , R u the largest stable steady state R * + , and R l is the boundary between the unstable and the stable regimes (black square symbol). We note that similar scaling laws to R l and R u have been previously found in binary mixtures i.e. without the solvent C [23]. There exists also a critical forward rate constant k c above which the slope of J in→out is so large that all drops dissolve (J in→out > J out→in for all R). k c is maximally bounded as follow [30]: When k > (φ−P out )/P out h the conversion P → k S is so fast that the system is outside the phase-separating region ("♦" symbol in Fig. 4).
We show the stability diagram in Fig. 9 (lower region, small drop regime) for varying forward rate constant k and drop radius R. A multi-drop system at steady state R is unstable against Ostwald ripening if R c < R < R l (upward arrows), and stable if R l < R < R u (grey region). Figure 9. Stability diagram of a multi-drop system for varying forward rate constant k and drop radius R. A steady state exists within the continuous curve. Outside this region all drops dissolve (downward arrows). Outside the grey region but within the continuous curve, drops are unstable against Ostwald ripening and coarsen (upward arrows). In the grey region drops are stable against Ostwald ripening. The stability-instability boundary is shown by a dashed curve. The analytical results for the critical radii in the small drop regime are showed by the upper dotted line (Ru, Eq. (14)) and lower dotted line (R l ). The analytical results for the critical rates kc, k l and ku (Eqs (15), (20) and (21)) are showed by arrows. Large drops dissolve beyond ku, all drops dissolve beyond kc, and large drops are stable against Ostwald ripening beyond k l . Parameters: where ν is the molecular volume of P and S and can be chosen arbitrarily.
As the rate constant k decreases the drop radii tend to increase (R u , R l ). When the radii get comparable or larger than the gradient length scale ξ, the system exits the small drop regime and enters the large drop regime, which we will now discuss.

Large drop regime.
We now concentrate on the regime where drops are large compared to the gradient length scale (R ξ, Eq. (11)). The quantity ξ/r in the profile expression Eq. (10) is therefore small in the cytoplasm (r > R). Then, enforcing the no-flux boundary condition in drop centre and the Gibbs-Thomson relation at the interface (Eqs (2)), the profiles become [23,30]: and S in,out (r) = φ in,out − P in,out (r), with φ in being independent of r (Eq. (9)). We show in Fig. 10 a schematic of the concentration profiles of P and S along an axis passing through a drop centre. Far from the interface the profiles are flat, implying that P and S are in chemical equilibrium (∇ 2 P = ∇ 2 S = 0 in Eqs (7) and (8)). Concentration gradients and fluxes exist only in the interfacial region. Interestingly this large drop regime cannot exist in the absence of the solvent  11)) are localised at the interfacial region only. Hence the reaction-diffusiondriven circulation of molecules P and S between drops and cytoplasm (blue↔red circular arrow) is also localised at the interfacial region only. See Fig. 7 for details about this circulation.
(C), i.e. in a binary mixture. Indeed, imposing in this case the incompressibility condition P in,out (r) + S in,out (r) = ψ with ψ being independent of r, and the chemical equilibrium condition S in,out (r) = P in,out (r)k/h, there is a unique solution for the solute concentration: the solute concentration, both inside and outside drops away from the interfaces, equals the overall solute concentration in the whole system. Since the solute concentration inside and outside drops is supersaturated ("♦" symbol in Fig. 4), new drops will be created through further phase separation, ultimately leading the system to the small drop regime (Sec. 2.3.2). In our ternary mixture however the incompressibility condition is P in,out (r) + S in,out (r) + C in,out = ψ with C in,out independent of r due to Eq. (9). Since C can distribute differentially inside and outside drops (C in = C out ), and adding now the chemical equilibrium conditions S in,out = P in,out k/h we find two distinct solutions for the solute concentration.
In other words the added degree of freedom from the solvent allows for distinct chemical equilibrium concentrations inside and outside drops, as shown in Fig. 10 (P in (0), P out (∞)).
Using the solute concentration profiles Eqs (16) and (17) we obtain the influx of solute P (J out→in = DdP out /dr| R ) and the efflux of molecules S (J in→out = −DdS in /dr| R = DdP in /dr| R , Eq. (9)): Neglecting for the time being the terms ξ/R which are small in this regime, we find that both fluxes are simply proportional to the difference of concentrations between the interface and far from the interface. As the drop radius R increases, the solute cytoplasmic concentration P out (R) close to the interface decreases according to the Gibbs-Thomson relation (Eq. (2)). This causes an increase of the solute influx J out→in into the drop, and therefore contributes to further drop growth. Hence, as in equilibrium systems (i.e without chemical reactions, Sec. 2.2), the Gibbs-Thomson relation supports Ostwald ripening.
Let us now study the effect of the term ξ/R in the flux expressions Eqs (18) and (19). The term ξ/R, which vanish at large radius R, captures the concentration profile asymmetry inside and outside drops (Eqs (16) and (17)) arising from the spherical drop shape. As a result drop expansion tends to increase the influx J out→in while decreasing the efflux J in→out , acting against further expansion. The chemical reaction-induced term ξ/R therefore tends to stabilise a multi-drop system against Ostwald ripening.
In summary in this large drop regime we found that the Gibbs-Thomson relations tend to destabilise a multi-drop system causing Ostwald ripening, and chemical reactions on the contrary have a stabilisation effect, as in the small drop regime (Sec. 2.3.2). In Appendix B we provide a quantitative analysis of Eqs (18) and (19) and find the critical forward rate k l beyond which Ostwald ripening is arrested: and the maximal forward rate constant k u above which drops dissolve: and we have recovered the results from [30]. Note that for k > k u drops may still exist in the small drop regime (Sec. 2.3.2). We show in the stability diagram in Fig. 9 [30] the stable and unstable regions depending on the forward rate constant k and drop radius R We also note that there is no upper-bound on the drop radius R contrary to the small drop regime (Sec. 2.3.2).

Spatial organisation
Another interesting phenomenon resulting from this type of non-equilibrium phase separation is the potential spontaneous spatial organisation of drops on a lattice, as observed in Monte Carlo simulations shown in Fig. 11a [30]. In this section we provide a simple intuitive argument that accounts for the observed lattice organisation. Note that these simulations do not capture momentum transfer between molecules, hence drop diffusion is not reproduced accurately. Whether this lattice structure survives in the presence of drop thermal diffusion remains to be investigated. At steady state we have seen that a flux balance exists at the interface between the influx of solute P and the efflux of S (Fig. 7). We now consider a system where drops are regularly distributed on a lattice initially, and study what happens when a drop moves away from its initial position (Fig. 11b)). Since the reaction S → P dominates in the dilute phase between drops, this leads to an accumulation of molecules P in the cytoplasm, and a subsequent diffusive flux of P toward drops (Fig. 7). As the time required for P molecules to reach drops increases with the inter-drop distance, P accumulates predominantly where the inter-drop distance is large. A similar intuitive argument was also provided to explain the concentration profiles inside drops (Eq. Figure 11. Drop spatial organisation. a) Snapshot of Monte Carlo simulations of a multi-drop system at steady state (no momentum transfer between molecules, see main text for details). Drops spontaneously distribute homogeneously in space. b) Schematic of the solute profile across three adjacent drops (red curves). The influx of molecules P is shown by red arrows and the efflux of molecules S by blue arrows. The solute profile gradient between drops far from each other is stronger than between close drops. The solute influx is therefore larger on one side of the drop (thick red arrow) and weaker on the opposite side (thin red arrow), causing the drop to grow on one side and dissolve on the other side. This leads to effective drop displacement tending to align drops on a lattice. (12)). In addition, due to the concentration P at a drop interface being fixed by the Gibbs-Thomson relation (Eq. (2)), the concentration gradient of P is increased on the side of the drop where the inter-drop distance is the largest. The influx of solute on this side of the drop is therefore strengthened (thick red arrow). On the contrary on the side of decreased inter-drop distance, a weaker solute gradient leads to a smaller solute influx (thin red arrow). The imbalance of the solute influx on each side of the drop causes effective drop displacement toward the larger influx, i.e., toward the drop initial position. As a result, chemical reactions in our multi-drop system tend to distribute drops on a lattice structure.

Active matter: motile organisms in the incompressible limit
Active matter refers to physical systems in which some or all constituents of the system can exert forces continuously on their surrounding environment [34]. For instance, in the case of a bird flock, the birds fly by flapping their wings to move the air around them; in the case of a cell tissue on a substrate, the cells move via coordinated and ATP-driven remodelling of biopolymers beneath their cell membranes [35]. Active matter constitutes a non-equilibrium system and the energy is provided either through a continuous supply of fuel or by energy already stored in the system.
Here, we will focus exclusively on active matter in the condensed state, to the extent that the system can be viewed as incompressible. Biological examples include a dense collection of motile bacteria [36], or a cell tissue in which the cells are undergoing dynamic rearrangement [37,38]. In the hydrodynamic limit (the limits of long time and long wavelength), an active matter can usually be described by equations of motion (EOM) of field variables that correspond to some coarse-grained properties of the system, such as the local density and the local velocity [39,40,34]. Such an EOM can generically be written down based on symmetry consideration alone and the associated universal behaviour of the system can then be analysed using analyical methods such as dynamical renormalisation group (DRG) methods [41,42], or numerical methods. In this review, we will focus on the former approach and discuss how it enables us to elucidate the universal behaviour of an incompressible active fluid at criticality and in the ordered phase in certain spatial dimensions.

Hydrodynamic theory of incompressible passive fluids -Navier-Stokes equation
For an equilibrium system, symmetry constrains the allowable form of the Hamiltonian of the system [43]. For a non-equilibrium system, although a Hamiltonian may no longer be relevant, we can still use symmetry to deduce the form of EOM [44,45]. To illustrate this approach, we will now review how such a symmetry consideration can help us derive the incompressible Navier-Stokes equation.
In an incompressible fluid, the obvious field variable is the velocity field v(r, t), whose dynamics can be written as: where ρ is the density field and F corresponds to the local force density. Since the system is incompressible, ρ is constant everywhere and we will ignore this constant factor from now on. We now impose the following symmetries: (i) Temporal invariance: F does not depend on time t explicitly, hence forbidding terms like tv. This symmetry means that experimental results on the fluid motion do not depend on the day of the week on which the experiments are done. (ii) Translational invariance: F does not depend on the spatial location r explicitly, hence forbidding terms like r. This symmetry means that experimental results do not depend on the location where the experiments are done (iii) Rotational invariance: the EOM is invariant if the reference frame is rotated, hence forbidding terms like w for some constant vector w since a term like this will pick out a particular direction in the system. This symmetry means that experimental results do not depend on which direction the experimental apparatus are positioned towards. (iv) Parity invariance: the EOM is invariant under spatial inversion, hence forbidding terms like ∇ × v. This symmetry means that the physical system has no chirality, i.e., the physics of fluid motion has no handedness. Imposing these symmetries, and expanding F in powers of v and of the spatial derivatives ∇, we arrive at the generic EOM: where v ≡ |v| and ". . ." refer to higher order terms permissible in F that are not shown. Note that the first term on the R.H.S. of the above equation, − κ, is a vectorial Lagrange multiplier there to enforce the incompressibility condition ∇ · v = 0. By the Helmholtz decomposition, we can write F as ∇p + ∇ × A where p is a scalar field and A is a vector field [46]. Since we want to subtract off part of F that is not divergence-free, we have κ = ∇p.
Our EOM so far does not look like the Navier-Stokes equation yet as we are still missing one crucial symmetry: the Galilean invariance.
(v) Galilean invariance: when no external forces are acting on the system, the EOM is invariant if the reference frame is boosted to another reference frame that is travelling at a constant speed in an arbitrary direction.
Under this additional symmetry, the EOM remains invariant if we perform the following simultaneous transformations: r → r − wt and v(r, t) → v(r − wt, t) + w, for some arbitrary vector w. Imposing this constraint, the EOM, to order O(∇ 4 ), is which is exactly the incompressible Navier-Stokes equation, with p interpreted as the pressure divided by the density. If we are only interested in the coarse-grained (long wavelength) behaviour of the system, we can argue that higher order terms, such as ∇ 4 v, are unimportant compared to ∇ 2 v, and thus Eq. (24) can be viewed as the hydrodynamic equation of incompressible fluids. In a physical system, fluctuations, e.g., thermal in origin, are inevitable. Using the fluctuation-dissipation relation, fluctuations can be added to the above Navier-Stokes equation that renders it suitable to describe incompressible thermal fluids [47]. Analytical treatment of the resulting stochastic partial differential equation can then be used to elucidate the universal behaviour of the system. For instance, the existence of long-time tail of various correlation functions of thermal fluids, first discovered via simulations [48,49], have been confirmed using diverse analytical methods such as kinetic theory [50,51] and DRG analysis [42]. We will not review these well known results here, instead we will now turn our attention to incompressible active fluids.

Incompressible active fluids
We will focus exclusively on the so-called "dry" active matter [39,34], in the sense that there exists a fixed background in the system for the active constituents to exert forces on. Experimentally, the active constituents can be motile cells and the fixed background can be a gel substrate that the cells crawl on. In contrast, wet active matter describes motile organisms in a fluid medium in which organisms move by exchanging momentum with the surrounding fluid, and the resulting fluid flow can in turn affect the motion of the organisms [52,53].
In dry active matter, due to the ability of each active volume element to generate forces against a fixed background, the Galilean invariance no longer applies. Omitting this symmetry, the general EOM of a generic incompressible active fluids is of the form of Eq. (23), which is in fact exactly the incompressible version of the Toner-Tu equation devised to describe the flocking behaviour [44,45,54].
Ignoring the blue terms in Eq. (23) for the time being (whose omissions will be justified later), and focusing on spatially homogeneous states (so all terms involving Figure 12. In a generic incompressible active fluid, two distinct phases are possible. At the mean-field level, a disordered phase exists when the parameter a is negative, and an ordered phase (characterised by a non-zero mean speed of the system) emerges when a is positive. The transition between these two phases is continuous, or critical (region depicted in red). The surface plots depicts the "potential energy landscape" at the mean-field level for an active fluid in two dimensions. In the disordered phase, the energy landscape is like a parabolic bowl, while at the transition, the global minimum of the landscape becomes very flat. The landscape transitions further into the shape of a Mexican hat in the ordered phase.
∇ become zero), the simplified EOM can be written as where H(v) = −av 2 /2 + bv 4 /4. H can be viewed as a "potential energy" term, whose forms, depending on the parameter a, are depicted in Fig. 12 for a two dimensional system. When a is negative, H has only one minimum at v = 0, which suggests that the only steady-state solution is the v = 0 homogeneous state. We call this the disordered phase. As a increases beyond zero, a continuum of minima emerges and all of these will have a non-zero mean speed. This corresponds to the ordered phase, or the collective motion phase. The transition between these two phases is continuous and thus constitutes a critical transition. From equilibrium statistical mechanics, we know that when spatial heterogeneity and fluctuations are restored, the system can possess scale-invariant features at criticality [55,56,43]. This is also what happens in our active fluid system, as we shall show next.

Universal behaviour at the critical point
To understand the emergence of scale-invariant structures at the critical point when the system transitions from the disordered phase to the ordered phase, we will first analyse the EOM at the linear level and then incorporate the nonlinear effects using DRG methods.

Linear theory.
To arrive at the linear equation, we tune all the coefficients in the EOM to zero except for the terms below: where we have added the Gaussian noise term f . Since we are interested in an incompressible system, we would like the noise to be incompressible as well. In Fourier space, f (q, t) = d d re −ir·q f (r, t), the incompressibility conditions implies q · f = 0. Now, given any Gaussian noise termf with statistics we can use the transverse projection operator P ij (q) ≡ δ ij − q i q j /|q| 2 to define an incompressible noise term as f i = P ijfj . Since q · f = q i P ijfj = 0, f is incompressible as desired. In the Fourier transformed space, f has the statistics Note that the form of the noise term f also respects all of the symmetries (symmetries (i)-(iv)) imposed on our system. Furthermore, we no longer need the Lagrange multiplier in the linear EOM (26) as it is intrinsically incompressible.
To investigate the scale invariant properties of our linear model, we now perform the following re-scaling for some dimensionless number that describes how the spatial length scale is modified. The field variable v and time t are also re-scaled, albeit with distinct exponents: the roughness exponent χ and the dynamic exponent z, respectively. The numerical values of these two exponents are yet to be determined.
Applying the re-scaling to Eq. (26), we find where d is the spatial dimension. The prefactor in front of the noise term originates from the form of the noise term (28) and the fact that the delta function scales inversely to its argument, e.g., δ(t) → e −z δ(t).
Re-writing Eq. (30) as we see that the transformed equation is exactly of the form of the original EOM (26) except that the coefficients µ ≡ e (z−2) µ and D = e (z−2χ−d) D have acquired a dependency on . What it means is that if we re-scale the spatial coordinate, then the coefficients in the EOM will generically be modified. We can express the coefficients' dependencies on in the form of differential equations: We shall call the above the flow equations of the coefficients.
If we now pick z to be 2 and χ to be (2 − d)/2, then µ and D remain unchanged as changes. In other words, given this choice of the exponents, the coefficients in the linear EOM are invariant under re-scaling. The beauty of this invariance is that it enables us to obtain the power-law behaviour of the temporal and spatial correlation functions of the system [57]. For instance, we can relate the equal-time correlation function at different distance because where we have picked such that = ln r, and the second equality follows from the fact that the re-scaling of r can be absorbed by re-scaling the field variable v according to v → e χ v.
What we have seen is that in the linear theory, by suitably re-scaling the field variable and time, the coefficients in the EOM will remain invariant under spatial rescaling, which leads to a power-law behaviour of the correlation function. Importantly, the power law exponents follow purely from the structure of the equation, and are independent of the actual coefficients in the EOM. We will now look at how the incorporation of other terms in the EOM affects this conclusion. (23) is a stochastic, nonlinear partial differential equation, as such, it is notoriously difficult to analyse. Here, we will employ a DRG method to treat this problem analytically [41,42]. The strategy is to use the results from our linear theory as our reference point and then to incorporate the nonlinear effects perturbatively. To proceed, we will first employ the scaling exponents from our linear theory to gauge the importance of the additional terms in our full EOM. For example, by absorbing the scaling transformation on the term v 4 v into the coefficient c, we have c = e (4χ lin +z lin ) c = e (6−3d) c. If we now take to be big, then this term becomes small as long as the spatial dimension is greater 2. In fact, when nonlinearities are taken into account (42), this term can be shown to be negligible in the hydrodynamic limit even at d = 2. Namely, if we focus on the large distance properties of the system, this particular term will become negligible asymptotically as → ∞. This is also true for the (∇ 2 ) 2 v term since µ = e (−4+z lin ) µ = e −2 µ . In fact, all the blue terms in the EOM (23) can be shown to be asymptotically negligible if 2 < d < 4 according to the scaling exponents from the linear theory. We call these terms irrelevant.

Nonlinear effects. The full EOM
How about the advective term (v · ∇)v? One can readily see that λ = e (χ lin −1+z lin ) λ = e (4−d) /2 λ, which means that this term becomes ever more important if d < 4 as → ∞. The same applies to the term v 2 v as b = e (2χ lin +z lin ) b = e (4−d) b. These two terms are therefore relevant for d < 4.
What we have seen so far is that if d is below 4, d = 3, say, then in the hydrodynamic limit ( → ∞), the full EOM can be reduced to Note that in contrast to the linear theory, in which the scale-invariant properties hold for any , positive or negative, the full EOM only gets simpler as → ∞. In Eq. (34), we have also re-instated the linear term av as it is needed in order to fine-tune the system to criticality. The role of a is similar to the role of temperature in the Ising model, which need to be fine tuned to lead the system to the critical point [55,56,43].
As mentioned before, using the exponents from the linear theory, the coefficients in Eq. (34) only remain invariant upon re-scaling at d = 4. For d < 4, both the λ and b terms diverge upon zooming out spatially ( → ∞). Looking back at Eqs (32), we have implicitly assumed that the two flow equations are uncoupled from each other, and to find the fixed point for these two equations, it was enough to use the two free variables, χ and z, to make the R.H.S. of the two equations zero. However, now that we have four coefficients in our reduced EOM (λ , µ , b and D ), two variables will generically not be enough to set the four flow equations to zero. The problem here lies in the fact that how the coefficients vary upon increasing are not independent. From the EOM, we have Note that we have left out the coefficient a associated to the term v in the EOM. This is because a is our fine-tuning parameter to take the system to criticality, and as such, the flow of a does not couple with the flows of other coefficients at the critical point [58]. At d = 4, there are only two dimensionless quantities that can be constructed out of these four coefficients. Two particular choices are: If coupling between the flow equations were to occur, then the coupling terms can only be functions of the above two coupling coefficients.
Going below four dimensions, the above two quantities are no longer dimensionless (because the dimension of D depends on d), so we have to find other dimensionless coupling coefficients. In addition, we will have to write down exactly how these coupling coefficients enter the flow equations of the coefficients. These two tasks are dealt with using a DRG transformation together with the -expansion method [41,42]. Under a DRG transform, fluctuations associated with the short distance behaviour of the system are averaged over and the effects of the averaging are then incorporated back into the EOM. Since we have started with the most general EOM possible allowed by symmetry, the form of the resulting EOM must remain the same. Hence, the effects of a DRG transformation can only lead to the modifications of the coefficients in the EOM, which will be encoded in the flow equations of the coefficients. The -expansion method refers to the fact that we are calculating the nonlinear effects of the system perturbatively, where the perturbation is from the linear theory (which is applicable when d = 4), and the small parameter is the deviation from the dimension below which the linear theory breaks down, i.e., = 4 − d in our case. A detailed discussion of the DRG calculations for our system is beyond the scope of this review and we will only quote the results here [58]: where and S d = 2π d/2 /Γ(d/2) is the surface area of a unit sphere in d dimensions, and Λ is some fixed short wavelength cutoff of dimension [L] −1 , whose inverse corresponds roughly to the size of the motile organism under consideration, e.g., it can correspond to the cell diameter in a tissue. Note that the form of the coupling coefficients g in Eq. (38), obtained here from the DRG calculation, are the same as those in Eq. (36), except for the introduction of Λ to keep the g's dimensionless.
Since these two dimensionless quantities themselves vary with , we can study their own flow equations. By their definitions (38), we have Using the expressions in the flow equations (37), we can re-write the R.H.S. above as Eqs (40) now define two coupled ordinary differential equations, and their dependencies on are shown in Fig. 13. As → ∞, the coupling coefficients flow to three fixed points. Two are unstable (depicted by the purple circle and the blue triangle in Fig. 13) and were already known to physicists in the context of ferromagnetism with dipolar interactions [59,60,61] and randomly stirred fluids (model B in [42]). The stable fixed point (red square), in the sense that all flow lines around the fixed point converge to it as increases, was a novel addition whose discovery was solely motivated by the study of active matter [58]. We will now focus on this fixed point.
At the stable fixed point, and the R.H.S. of Eqs (40) are zero and thus enforce two algebraic constraints on the flow equations of the coefficients (37) via Eqs (39). If we now pick z and χ such that two of the flow equations are zero, those of µ and D , say, then together with the two previous constraints, they imply the invariance of all coefficients as → ∞. In this asymptotic regime, the system will again exhibit scale invariant properties as in the linear theory, except now, due to the nonlinearities in the system, the roughness and dynamics exponent are modified to In particular, in contrast to the scaling predicted by the linear theory in Eq. (33), when the nonlinear effects are taken into account, the equal-time correlation of the system actually follows the power-law: Note that we focus here on the system's behaviour right at the critical point so that the correlation length of the system is infinite. Away from criticality, there is another exponent that quantifies how the correlation length varies as one gets close to the critical point. We refer the interested readers to [58] for further details.

Universality class.
We have seen that by a symmetry consideration alone, we can derive a model EOM that generically describes incompressible active fluids. By focusing on the large-distance behaviour (the limit → ∞), we have managed to calculate the scaling behaviour of the system in spatial dimensions 4 − . Importantly, at no point did we need the inputs of the actual model parameters. Therefore, the exponents governing the scaling behaviour depend only on the symmetry of the system, and are otherwise oblivious to the quantitative parameters that describe the actual systems. This suggests that a large class of dynamical systems having the same form of EOM, but with distinct model parameters, will exhibit identical scale invariant behaviour. These distinct systems are said to belong to the same universality class.

Ordered phase in two dimensions
We have seen that at the critical transition, the scaling behaviour of a generic incompressible active fluid constitutes a novel universality class in non-equilibrium physics. Here, we will describe how in two dimensions, the ordered phase in incompressible active fluids also exhibits universal behaviour, albeit with scaling behaviour that belongs to a well known universality class: the Kardar-Parisi-Zhang (KPZ) universality class that originated from modelling surface growth in the nonequilibrium regime [62].
3.4.1. Linear theory. As before, the symmetry consideration alone has fixed the structure of the EOM. The difference from the previous section is that while at criticality, v is small and thus an expansion with respect to v makes sense, this is no longer true in the ordered phase since v can be large. As a result, the generic EOM in the ordered phase is different from the EOM in Eq. (34) [63]. We will however continue to use the restricted form (34) for simplicity and our conclusion remains the same even for the generic case [63]. Without loss of generality, we can assume here that the collective motion is along the x-direction and re-write the field variables v in terms of the deviation field u from the mean velocity and we assume u ≡ |u| is much smaller than v 0 . Re-writing the EOM in terms of u and then keeping only the leading order linear terms, we have where v 0 = a/b, and p(r, t) again serves to impose the incompressibility condition ∇ · u = 0. Specifically, we have In the Fourier-transformed space, where q ≡ |q|.
In Eq. (45), we have also allowed for distinct "viscosity" coefficients µ x,y for the two directions since in the ordered phase, the rotational symmetry is broken and there is no reason to expect that the scaling behaviour will be identical in both x and y directions. For the same reason, we will allow for two distinct direction-dependent roughness exponents χ x and χ y , as well as an anisotropic exponent ζ in our re-scaling scheme: x → e x , y → e ζ y , u x → e χx u x , u y → e χy u y , t → e z t .
Although there are now seemingly five distinct exponents, they are not all independent as some of the exponents are related via the incompressibility condition: q x u x = −q y u y → e (χx−1) q x u x = −e (χy−ζ) q y u y , which implies that in the hydrodynamic limit, In two dimensions, the "potential energy" term H is depicted in Fig. 12, whose functional derivative with respect to v gives rise to the terms (a + b 2 v 2 )v in the EOM. When the symmetry is broken, we know from equilibrium physics that fluctuations in v y dominate over those in v x because moving along the trough of H costs much less energy than moving up the valley in H. We assume the same principle applies to our non-equilibrium problem here, which we will verify a posteriori.
Focusing therefore on u y , we first go to the boost frame u(x, y, t) → u(x − λv 0 t, y, t) to eliminate the second term on the L.H.S. of Eq. (45). Without this term, we have where (50b) follows from (47). We now use the incompressibility condition q x u x + q y u y = 0 to express u x as −q y u y /q x . Eq. (50b) can thus be written as Upon performing the re-scaling as in Eq. (31), we have We have argued previously that fluctuations in u y dominate over those of u x , we thus expect that χ y > χ x , this in turn leads to the inequality ζ > 1 via Eq. (49). Therefore, q 2 → e −2 q 2 x + e −2ζ q 2 y ∼ e −2 q 2 x . This explains the prefactor in the second term on the R.H.S. of Eq. (52). Furthermore, we can conclude that µ x q 2 y u y , the first term in the R.H.S. of Eq. (52), is not as important as the third term upon coarsegraining. Ignoring the first term, we can make the remaining three terms invariant upon re-scaling by choosing the following exponents: In addition, Eq. (49) implies χ x,lin = −3/2.

Nonlinear effects.
Going back to the simplifed EOM (34), we can now gauge the relevance of the two distinct nonlinear terms. We consider first the advective term: λ(v · ∇)v. Focusing again on the EOM of u y , the nonlinearities coming from this term are λu x ∂ x u y and λu y ∂ y u y . Upon re-scaling, these terms become e (z+χx−1) λu x ∂ x u y and e (z+χy−ζ) λu y ∂ y u y . Using the exponents from the linear theory, we can see that both terms are irrelevant as → ∞. On the other hand, the nonlinearities associated to the coefficient a are of the form which, upon re-scaling, becomes In this case, both the first and third terms are relevant and we will thus keep these two terms.
In summary, for an incompressible active fluid in two dimensions, the governing EOM (34) can be further simplified, in the hydrodynamic limit, to where we have omitted the subscript y in µ above to ease notation. Note that we have not considered the role of the pressure term p in detail in the above analysis, however our conclusion remains the same even if we do so [63].

3.4.3.
Mapping active fluids in two dimensions onto an equilibrium system. To make analytical progress with the EOM (56), we will first map the class of systems described by (56) onto an equilibrium system. As the advective term λ(v · ∇)v is irrelevant in the ordered phase, its omission enables us to re-write the hydrodynamic EOM Eq. (56) as where Viewing H as the Hamiltonian of the system, and the noise term f as thermal fluctuations, we can analyse the static, or equal-time, properties of the system by studying the corresponding partition function: where β = 1/D (27) and D 2 u ≡ lim N →∞ N k du x (r k , t k )du y (r k , t k ), with the index k enumerating the discretised but infinitesimal elements of both space and time. Note that the incompressibility constraint ∇ · u = 0 is enforced by the delta function δ(.) in the partition function above. Since we are in two dimensions, we can further eliminate the incompressibility constraint by using the stream function h(x, y) defined as In terms of h, the Hamiltonian H becomes which is exactly the Hamiltonian that describes a dislocation-free smectic liquid crystal in two dimensions (Fig. 14), with B ≡ 2av 2 0 and K ≡ µv 2 0 being the compression modulus and the bending modulus, respectively [64]. In a smectic liquid crystal, the liquid crystals (depicted as red ellipsoids in Fig. 14) formed a layered structure in which the layers are parallel to x-axis on average and h(x, y) describes the height deviation of the layers from the expected location. The word smectic comes from the Greek word for soap, whose slippery surface upon lubrication consists of layered lipid bilayers which can slide across each other easily.
Note that since H S has to be invariant under a rotation with respect to the axis perpendicular to the xy-plane, the compression term should also include the term −(∂ y h) 2 /2 in the curly brackets. However, such a term is irrelevant in the hydrodynamic limit and thus omitted.
Given the Hamiltonian H S , at thermal equilibrium, the probability of having a particular profile h(x, y) is given by where the partition function is 3.4.4. Kardar-Parisi-Zhang universality class. Besides the surprising connection between incompressible active fluids and smectics in two dimensions, we will now make one further connection to a well known physical model: the KPZ surface growth model, based on a mapping devised in [65,66]. To do so, we first add the two boundary terms below (hence immaterial due to the assumed periodic boundary condition) to the smectics Hamiltonian: where the last equality follows from an integration by parts on the first term. Adding these immaterial boundary terms to the smectics Hamiltonian enables us to "complete the square" in the integrand and to re-write H S as Employing the delta function δ[.], we can now re-write the partition function (63) as where η may be viewed as an auxilliary field to be integrated over. However, another way to interpret η(r, t) is to view it as a random variable for any particular r and t, such that its probability density is proportional to e −βη(r,t) 2 . From this perspective, we can generate the same statistics of height profile h as in Eq. (62) by using the Langevin equation: where η(x, y) = 0 , η(x, y)η(x , y ) = 2Dδ(x − x )δ(y − y ) .
Eq. (67) is exactly the (1+1)d Kardar-Parisi-Zhang (KPZ) surface growth model if one interprets the symbol y as time and x as the one spatial dimension in the KPZ equation.
The roughness and dynamic exponents of the KPZ model are known exactly: z KPZ = 3/2 and χ KPZ = 1/2 [62]. Translating this back to our active fluid model, the KPZ dynamic exponent becomes the anisotropic exponent: ζ = 3/2. Furthermore, via the relationships between u and h in Eq. (60), we have Compared to the exponents obtained from our linear theory (53), we see that the anisotropic exponent ζ, and as a result also the roughness exponent along the x direction χ x , are modified due to the nonlinear terms in the EOM. Furthermore, since χ y > χ x , we can now verify the previous assertion that u y dominates over u x in the hydrodynamic limit (Fig. 15). Looking at the re-scaling scheme in (48), the only exponent that remains unknown is the dynamic exponent z. As we have mapped the active fluids onto an equilibrium system and then focused on the partition function, the connection established actually does not allow us to determine the dynamic exponent z. Hence, what the value of z is remains an interesting open question.

Conclusion & Outlook
The behaviour of a system undergoing a phase transitions is an archetype of collective phenomena, and the emergence of universality at a phase transition makes it a fascinating subject for physicists. Motivated by recent studies focused on phase transitions in biological systems, we have discussed how novel physics can arise from the generic non-equilibrium nature of living matter, be it driven chemical reactions or self-generated mechanical forces. For the former, we have seen in Sect. 2 how driven chemical reactions in the cell cytoplasm can stabilise protein drops that form by phase separating out of the cytosol, contrary to the universal coarsening behaviour found in phase separating systems at thermal equilibrium. In Sect. 3, we focused on incompressible active matter and we have shown that a novel universality class emerges at the critical order-disorder transition due to the activity of the system. Furthermore, in the ordered phase in two dimensions, we have discussed the surprising connections between active fluids, smectics, and the Kardar-Parisi-Zhang surface growth model.
Besides the two particular systems discussed in this review, there are other examples of phase transition phenomena in biology that have led to the discovery of novel physics. Here, we will mention two particular examples that are receiving intense attention from physicists.
(i) Motility-induced phase separations. Motivated by the run-and-tumble motion of bacteria [68], the study of active particles that interact solely via volume exclusion led to the discovery of liquid-gas phase separation driven purely by motility [69,70,71,72]. The condensed drops in motility-induced phase separations (MIPS) show interesting interfacial properties, such as the existence of a stable interface with a negative surface tension [73]. Away from the critical point and close to the phase boundary, the coarsening kinetics of MIPS has been argued to be identical to the Liftshitz-Slyozov scaling law found in equilibrium system [74]. However, what the universal behaviour of MIPS is at criticality remains an interesting open question. (ii) Active polymer networks. Throughout the review, the non-equilibrium processes, be they active motion or chemical reactions, occur in a homogeneous environment that corresponds to an Euclidean space. What happens if active forces are now transmitted through a network of irregular structures instead? A recent discovery found that a biologically relevant active polymer network under fragmentation can self-organise itself to exhibit a scale-invariant signature of a critical system [75,76]. While the exponents observed are close to that of the static percolation universality class [77,78], the question of whether the critical phenomenon in active networks actually belongs to the static percolation universality class remains unsettled [79,80,81,82].
In terms of outlook, we believe the following future directions will expand the horizon of both biology and physics.
(i) In Sec. 2 we have studied how driven chemical reactions can stabilise a multidrop, ternary system. As the cell cytoplasm is a complex mixture of thousands of different molecules [83,84] it will be interesting to see how these results may be modified in many-component mixtures. (ii) In Sec. 2.4 we have provided intuitive arguments to explain the appearance of a lattice structure of phase-separated drops in our Monte Carlo simulations. Such a structure naturally suggests a kind of repulsive interactions between drops, which may serve to stabilise a multi-drop system against coarsening via coalescence due to drop diffusion. Whether this is indeed the case remains to be investigated.
(iii) In Sect. 3, we have studied the simplest kind of symmetry: the rotational symmetry and the associated universal behaviour when the symmetry breaks spontaneously in an active system. But what are the other relevant symmetries in biology, and will they also lead to novel universal behaviours?
Appendix A. Non-equilibrium phase separation in the small drop regime: solute concentration profile in the cytoplasm We analyse the cytoplasmic solute concentration P out (r) (Eq. (10)) in the small drop regime (R ξ, Eq. (11)), in two limiting cases. First, we assume the system contains few drops so that the inter-drop distance is large compared to the gradient length scale ξ. The term U (1) in is zero to avoid a diverging concentration far from drops (r ξ), and we fix the interface concentration (r = R) according to the Gibbs-Thomson relation Eq. (2): where P out (∞) is the concentration far from drops. We now assume that many drops are present so that the inter-drop distance is small compared to ξ. The quantity r/ξ is then always small and the exponential terms in Eq. (10) are close to one. Imposing again the Gibbs-Thomson relation at the interface we find P out (r) = P out (∞) + (P out (R) − P out (∞))R/r, which is equivalent to Eq. (A.1) for vanishing r/ξ. Therefore Eq. (A.1) is a good approximation for the cytoplasmic profile P out (r) in both regimes of small and large drop number.
The cytoplasmic profile P out (r) is closely similar to that at equilibrium (i.e without chemical reactions, Sec. 2.2), for identical supersaturation ∆. This can be seen by comparing Eqs (4) and (A.1) in the vicinity of and far from the drop interface (r ≈ R and r R). The influx of solute entering the drop, J out→in ≡ DdP out /dr| r=R is Since R ξ is small by definition in the small drop regime we recover the solute flux in absence of chemical reactions Eq. (5).

Appendix B. Non-equilibrium phase separation in the large drop regime
In the large drop regime the drop radii R are large compared to the gradient length scale ξ (Eq. (11)). The concentration profiles (Eqs (16) and (17)) are shown in Fig.  10 and the molecular fluxes at the drop interfaces are given by Eqs (18) and (19). We will first study the system at the steady state, all drops having the same radius R * , and then we will analyse its stability against Ostwald ripening.

Appendix B.1. Steady state
Besides our assumption that S does not phase separate i.e. S in (R) = S out (R) (see Ref [30] for the more general analysis without this assumption), we enumerate the other constraints on the system at steady state: (i) The Gibbs-Thomson relations dictate the concentrations at the interfaces (Eqs (1) and (2)), which follow from the assumption that the system is close to equilibrium to the extent that local thermal equilibrium is true (see main text). (ii) Drops and cytoplasm are at chemical equilibrium: in this regime drops are large compared to the gradient length scale ξ, and since we only focus on low drop density systems the same is true for the cytoplasm. Therefore in both phases and far away from interfaces the concentration profiles must be flat. Taking ∇ 2 P = ∇ 2 S = 0 in the reaction-diffusion equations Eqs (7) and (8) it follows that the concentrations are at chemical equilibrium. This imposes the following constraints on the solute concentration in the drop centre (P in (0)) and far from drops (P out (∞)): The solute mass conservation imposes the relation: where N in the number of drops, V the total volume of the system and P tot = φh/(k + h) is the global concentration of P in the entire system [30]. Note that we have neglected the concentration gradients near the interface due to the relative small size of the interfacial region (∼ ξ) compared to the size of the drops and the cytoplasm in this regime. (iv) The combined concentration P + S is homogeneous in drops and in the cytoplasm, as already shown in the main text, Eq. (9). This allows to express the concentrations P in (0), S in (0) in the drop centre and P out (∞), S out (∞) far from drops as functions of the interface concentrations: (v) The flux balance condition: the fluxes at the drop interfaces J out→in and J in→out (Eqs (18) and (19)) must balance each other at steady state so that drops neither grow nor shrink: where we have assumed infinite drop radius (ξ/R → 0).
Note that the conditions from (i) to (iv) also apply when drop radii are not at steady state because we assumed that the concentration profiles are always at the steady state (∂P/∂t = ∂S/∂t = 0, Sec. 2.2). Imposing all constraints above we can determine the steady state drop radius R * [30]: Contrary to the small drop regime (Sec. 2.3.2) the drop radius R * scales with the system size V , as in equilibrium systems i.e., systems without chemical reactions. This result shows that there exists a maximal forward rate constant k u above which drops dissolve (R * = 0) [30]: For k > k u drops cannot exist in the large drop regime but may still exist in the small drop regime (Sec. 2.3.2).

Appendix B.2. Stability of the steady state against Ostwald ripening
Let us now analyse the stability of our system against Ostwald ripening. We consider a multi-drop system initially at steady state, all drops radii being R * (Appendix B.1).
The molecular fluxes at a drop interface are given by Eqs (18) and (19). We now perturb the system by randomly increasing or decreasing each drop radius R by a fixed quantity , much smaller than R * . Let us concentrate on a specific drop i whose radius R i is increased: R i = R * + . The new interface fluxes for this drop are (Eqs (18) and (19)) Note that the values of the solute concentrations appearing in these expressions may also change due to the perturbation, according to the constraints in the system. If J out→in > J in→out the drop further expand, accentuating the initial perturbation and the system is therefore unstable against Ostwald ripening. On the contrary if J out→in < J in→out the drop shrinks back, implying that the system is stable. We will solve this problem by performing a linear stability analysis in the small parameter /R * . We first derive the dependency of the concentrations appearing in the flux expressions (Eqs (B.8) and (B.9)). Since the Gibbs-Thomson relations dictate the interface concentrations of the solute (constraint (i) in Appendix B.1) we find: P in (R i ) = P in (R * ) (B.10) P out (R i ) = P out (R * ) −P out l c (R * ) 2 + O(( /R * ) 2 ) . (B.11) Then, enforcing the constraints of chemical equilibrium far from the interface, and of homogeneous combined concentrations P + S both inside and outside the drop (constraints (ii) and (iv) respectively, Appendix B.1), we find P in (0) = P * in (0) +P out l c (R * ) 2 + O(( /R * ) 2 ) , (B.12) where P * in (0) is the solute concentration in the drop centre at the steady state (R i = R * ). To obtain this result we used the fact that the far-field concentrations P out (∞), S out (∞) are not affected by the perturbation.
Indeed, the random perturbation of the drop radii leaves the total drop volume unchanged. Then, the interface solute concentrations are perturbed according to the Gibbs-Thomson relations (Eqs (1) and (2)), thereby modifying the profiles near drops. However since the solute interface concentrations are randomly perturbed (up and down), we expect that the far-field profiles remain unperturbed up to the first order in the perturbation. We also used k h since k must be smaller than k u (Eq. (B.7)) in the large drop regime.
We can now express the fluxes as functions of the perturbation : To obtain this expression we used the results from the steady state (Appendix B.1) that lead to P out (∞) − P out (R * ) ≈ P in (R * ) − P * in (0) ≈ kP in /2. At small forward rate constant k the termP out l c originating from the Gibbs-Thomson relation (Eq. (2)) dominates over the chemical reaction induced-term kP in ξ/(2h). This implies the influx J out→in being larger than the efflux J in→out , and therefore drop growth. In this scenario the perturbation is amplified, and the system is therefore unstable against Ostwald ripening. On the contrary when k is large enough we have J out→in < J in→out , leading to drop shrinkage and the multi-drop system is stable. This confirms our intuitive arguments in Sec. 2.3.3 that chemical reactions tend to stabilise a multi-drop system. The critical rate k l of the transition between the unstable and stable regime is the solution ofP out l c = k lPin ξ/(2h), hence Eq. (20).