{'fraction_non_alphanumeric': 0.039321433938321285, 'fraction_numerical': 0.01218880967090214, 'mean_word_length': 4.495183044315993, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simultaneous Critical Order-Disorder and Critical Phase Separation in Compressible Polar Active Fluids: Dynamic and Functional Renormalization Group Analyses

The collective phenomenon of swarming behavior can, in the hydrodynamic limit, be described as a polar active fluid (PAF) - the same theory that governs the behavior of a collection of self-propelling particles with alignment interactions. Although a PAF is fundamentally out of equilibrium at the microscopic level, it is unclear whether the macroscopic behavior reflects this fact, e.g., through the violation of the fluctuation-dissipation relation valid in equilibrium systems. In this work, we investigate the multicritical point (MCP) of compressible PAFs, where the critical order-disorder transition coincides with critical phase separation. We first study the critical phenomenon using a dynamic renormalization group analysis and find that it is insufficient since two-loop effects are important to obtain a nontrivial correction to the scaling exponents. We then remedy this defect by using a functional renormalization group analysis. We find three novel universality classes and obtain their critical exponents, which we then use to show that at least two of these universality classes are out of equilibrium because they violate the fluctuation-dissipation relation.


I. INTRODUCTION
Active matter refers to many-body systems in which the microscopic constituents can exert forces or stresses on their surroundings and, as such, detailed balance is broken at the microscopic level [1,2]. However, even if the microscopic dynamics are fundamentally different from more traditional systems considered in physics, it remains unclear whether novel behavior will emerge in the hydrodynamic limits (i.e., the long time and large distance limits [3]). One unambiguous way to settle this question is to identify whether the system's dynamical and temporal statistics are governed by a new universality class (UC), typically characterized by a set of scaling exponents [4][5][6]. These exponents can in principle be determined using either simulation or renormalization group (RG) methods. However, simulation studies can be severely plagued by finite-size effects (e.g., two recent controversies concern the scaling behavior of active polymer networks [7,8] and critical motility-induced phase separation [9][10][11]). Therefore, RG analyses remain as of today the gold standard in the categorization of dynamical systems into distinct UCs. Indeed, for polar active fluids (PAFs) [12][13][14], an archetype of active matter systems, the use of dynamic renormalization group (DRG) [15] analyses have led to, on one hand, surprising realizations that certain types of PAFs are no different from thermal systems in the hydrodynamic limit [16,17], and on the other hand discoveries of diverse novel phases [13,14,[18][19][20][21][22][23][24][25] and critical phenomena [26,27]. However, due to the inherent technical difficulties in DRG methods, all of these studies have been restricted to PAFs in the incompressible or infinitely compressible (i.e., Malthusian) limits except for rare exceptions [19,20]. Further, when a DRG analysis was used in conjunction with the -expansion method, which was typically the case, it has always been restricted to the one-loop level.
In this work, we apply for the first time functional RG (FRG) methods on compressible PAFs and overcome some of these technical challenges. Specifically, we investigate a multicritical region of dry compressible PAFs. Although experimentally less accessible than simple critical points, multicritical points (MCPs) can offer surprising new physics, even in models that are thought to be well understood. For instance, nonperturbative fixed points have been discovered in the extensively studied O(N ) model [28], and in systems where two order parameters compete, whose individual critical points belong to equilibrium universality classes, the multicritical region where both critical points coincide can be manifestly out of equilibrium and demonstrate very interesting, spiral phase diagrams [29].
We will first apply a traditional one-loop DRG approach to the MCP of our interest, demonstrating how it is insufficient to capture its universal physics and then, for the first time for PAFs, apply a FRG [30][31][32][33][34][35][36][37] analysis that goes beyond the equivalent perturbative one-loop level.
The outline of this paper is as follows. In Sec. II, we introduce the hydrodynamic theory of compressible polar active matter and discuss salient features in its phase diagram, which enables us to define the multicritical point of interest. We then show how general scaling invariance of the equations of motion leads to powerlaw behavior in the correlation functions in Sec. III. For the multicritical point, we show this first in the linear regime in Sec. IV and then turn to the nonlinear regime in Sec. V. In Sec. VI, we perform the one-loop DRG calculation and argue why it is not sufficient to take into account the nonlinearities and then present the FRG approach in Sec. VII that we use instead. Using this method we find three RG fixed points which represent three novel nonequilibrium UCs. We discuss them and their scaling behavior in Sec. VIII. Finally, we summarize our findings and give an outlook on future work in Sec. IX.

II. COMPRESSIBLE POLAR ACTIVE FLUIDS
A. Equations of motion from symmetry and conservation laws Polar active matter aims to describe the collective behavior of swarming animals, e.g., flocks of birds, schools of fish or bacterial swarms. In the fluid state, as in passive fluids that are describable by the Navier-Stokes equations, the relevant dynamical variables are the momentum density and density fields, denoted by g and ρ respectively. Without any assumptions about the microscopic realization, one can then, based on symmetry and conservation laws, construct a generic set of hydrodynamic equations of motion (EOM) for these variables. Here, we assume that the particle number is conserved (as opposed to, e.g., a Malthusian system in which birth and death of particles can occur [18,22,23]). We thus arrive at a continuity equation as the EOM of the density field ρ: (1) For the momentum density field, we assume temporal, translational, rotational and chiral invariance. In addition, we focused on active systems in which the constituents move on a fixed frictional substrate, i.e., dry active matter systems (as opposed to wet active systems such as active suspensions) [2]. This symmetry consideration leads to the following generic hydrodynamic EOM of g [13,14,56]: where the ellipsis represents the omitted higher-order terms (i.e. terms of higher order in spatial derivatives and the momentum density field). In the EOM above, all coefficients are functions of ρ, and the noise term f (r, t) is a zero mean Gaussian white noise of the form The above hydrodynamic EOM (1,2) are termed the Toner-Tu EOM [13,14]. However, in contrast to the original Toner-Tu formulation, we have chosen to use the momentum field as the hydrodynamic variable instead of the velocity field so that there is a linear relationship between ρ and g, which facilitates our discussion later.

B. Mean-field theory and homogeneous phases
Given the hydrodynamic EOM, one of the first, and simplest, question to ask is: what are the mean-field solutions to the EOM? Answering this question amounts to focusing on temporally invariant, spatially homogeneous, and noise-free solutions to the EOM. These solutions are readily seen to be where β is taken to be always positive for reasons of stability. Since g can point in any direction, the |g| > 0 state implies a spontaneous symmetry breaking of the rotational symmetry, and corresponds to the homogeneous ordered phase of polar active matter where collective motion emerges. In contrast, the |g| = 0 state corresponds to the homogeneous disordered phase, i.e., there is no collective motion.
C. Phase diagram: phase separations, critical and multicritical phenomena The homogeneous phases discerned from the previous mean-field analysis are, however, not always stable, even in the absence of the noise term f (3). The standard way to ascertain the in/stability of the homogeneous phases in this noiseless regime is to use a linear stability analysis. Here, the temporal evolution of an initially small perturbation to a homogeneous solution is studied and the growth or decay of its amplitude signifies whether the homogeneous state is unstable or stable, respectively. Using ter, e.g., α0, and the average density ρ0, the system can be in the homogeneous disordered phase (white region denoted by D) or the polar ordered phase (yellow region denoted by O); It can also phase separate into two disordered phases with different densities (blue region), or into one ordered phase and one disordered phase, again with different densities (green region flanking the homogeneous ordered phase). The critical behavior associated with the first type of phase separation is generically described by the Ising UC (blue triangle) [10], and that associated with the second type is described by a putatively novel UC yet to be described (yellow inverted triangle) [58]. (b) Upon further fine-tuning, these two critical points can coincide (red circle) [59], and the resulting critical point is described by a novel UC uncovered in the present work.
this method, complex phase diagrams of polar active fluids have been uncovered [57][58][59]. In particular, distinct types of phase separations, i.e., an inhomogeneous state where two different phases co-exist, have been demonstrated, which cannot be captured by mean field analysis.
In particular, expressing α and κ in Eq. (2) as where δρ = ρ − ρ 0 with ρ 0 being the average particle density in the system, two disordered phases (with distinct densities) can co-exist if α 0 > 0 and κ 0 < 0 (blue region in Fig. 1(a)) [10], while an ordered phase can coexist with a disordered phase if κ 0 > 0 and α 0 < 0 (green region) [58]. Further, the system can become critical upon fine-tuning: if α 0 > 0 and κ 0 = κ 1 = 0, the resulting critical behavior belongs to the Ising universality class (UC) (blue triangle) [10], while if α 0 = α 1 = 0 and κ 0 > 0, the associate critical behavior corresponds to a yet to be characterized UC (yellow inverted triangle) [58]. Recently, a third type of critical behavior was identified [59], which corresponds to the merging of these two distinct critical points by simultaneously fine-tuning α 0 , α 1 , κ 0 and κ 1 to zero (red circle in Fig. 1(b)). The universal behavior of this new multicritical point is the focus of this paper.

III. SCALE INVARIANT EQUATIONS OF MOTION
It is generally expected that the EOM of general systems at critical points become invariant under rescaling of lengths, time and fields [4][5][6]. Hence, at the multicritical point (MCP), we expect that the EOM (1,2) are invariant under the rescaling r → re , t → te z , ρ → ρe χρ , g → ge χg , (6) for some exponents z, χ ρ and χ g , that are a priori not known. If this is the case, however, this defines a rescaling symmetry of the theory which the correlation functions have to obey as well. Take for example the densitydensity correlation function: (7) Choosing = − ln r, r being measured against some reference scale, we see immediately that where S ρρ (.) is a scaling function that only depends on the ratio t/r z which is invariant under rescaling (6). So we immediately see that, if the EOM are invariant under a rescaling transformation, the correlation functions will generally express powerlaw behavior. Likewise, this argument can be applied to the momentum-momentum correlation function C g (r, t) = g(r, t)g(0, 0) = r 2χg S gg t r z , (9) where S gg is again a scaling function with similar properties as S ρρ . Ultimately, we will demonstrate that the EOM does become scale invariant and determine the scaling exponents using a FRG analysis, but first, we will illustrate the scale invariance discussed here using the simple, but quantitatively incorrect, linear theory.

IV. LINEAR REGIME
In the linear regime, i.e., when the non-linear terms in Eq. (2) are neglected, the scaling behavior discussed above can readily be seen. Around the MCP at which critical disordered phase separation (blue triangle in Fig. 1a)) merges with critical disorder-order of a generic compressible PAF (yellow inverted triangle), |g| ≈ 0, such that the linearized EOM are where we have introduced the term characterized by ζ since, when κ 0 is fine-tuned to zero, this term is now the leading order term linear in ρ. In Eq. (10), we have redefined ρ to be δρ to ease notation, and we will continue to do so from now on.
A. Scaling exponents

Correlation functions
Upon rescaling time, lengths, and fields according to Eq. (6), the linearized EOM (10) become They thus remain unchanged if At the linear level we can therefore directly conclude that using the argument from Sec. III. Since the linearized EOM (10) are solvable analytically by performing a spatio-temporal Fourier transform, the scaling behavior of the correlation functions can in fact be demonstrated explicitly. This has the added advantage that the expressions of the aforementioned scaling functions, S ρρ and S gg , can be obtained in the form of integrals.
Specifically, by performing a spatiotemporal Fourier transform, the linear EOM can be written as where g (q) = g(q) ·q, g ⊥ = g − g q, withq being the unit vector in the direction of q, q = |q|,q = (q, ω), and the G's in (14), or the "propagators", are: where µ ≡ µ 1 + µ 2 . Given the above expressions, the correlation functions are now obtained straightforwardly: where q ≡ d d qdω/(2π) (d+1) ,q ·r = q · r − ωt, and P ij (q) ≡ q i q j /q 2 and P ⊥ ij (q) ≡ δ ij − q i q j /q 2 are the projectors parallel and transverse to q respectively.
Focusing on C ρ (r, t) (16a) as an example, the substitutions ω = Ω/r 2 and q = Q/r in the integral lead to which demonstrates the scaling form (8) with the scaling exponents from the linear theory (12), and As we will show later, all the scaling exponents from the linear theory (12) are in fact incorrect for describing the hydrodynamic behavior around the MCP due to the nonlinearities in the EOM.

Divergence of correlation length
Besides the scaling exponents in the correlation functions right at the MCP, the divergence of the correlation length, as one approaches the MCP, is also governed by another set of scaling exponents. For the Ising model, this is the temperature. For the present MCP however, this divergence is associated to two parameters α 0 , κ 0 . The other two relevant parameters, α 1 and κ 1 , take a role akin to the magnetic field in the Ising model.
As the divergence of the correlation length is described by two parameters, α 0 and κ 0 , there are also two exponents which we call y 1 and y 2 . At the linear level, these exponents correspond expectedly to their mean-field values: We note that these scaling exponents are again expected to be modified by the nonlinearities, as we shall see in the next section. Experimentally, these exponents define a relationship between the correlation length ξ and the distances t 1 and t 2 in the phase diagram (see the inset of Fig. 2) from the critical points of disordered phase separation (blue upwards triangle in Fig. 1 and blue line in the inset of Fig. 2) and the critical order-disorder transition (yellow downwards triangle in Fig. 1 and yellow line in in the inset of Fig. 2), This means that, upon approaching the MCP, one has to enforce the relationship given by the right proportionality in Eq. (20) to see a clear scaling behavior in the correlation length. This relationship is also visualized in  shows the critical scaling behavior characterized by the scaling exponent χρ (slope triangle). In reality though, the system will never be exactly at the critical point. This can be characterized by the distances t1 and t2 (gray lines in inset) from the critical point of disordered phase separation (blue line in inset) and from the critical order-disorder transition (yellow line in inset) respectively. This manifests in a finite correlation length ξ above which the scaling behavior breaks down (gray lines in main figure). As one approaches the fixed point, by decreasing t1 and t2, the correlation length diverges. If this is done carefully, such that the second relation in Eq. (20) remains unchanged, for example by rescaling t1 and t2 by a factor s y 1 and s y 2 respectively (gray arrow in inset), the correlation length rescales according to Eq. (20), i.e. by a factor s −1 (gray arrow in main figure).

V. NONLINEAR REGIME
While the scaling behavior described in Sec. IV is qualitatively expected generally, the critical exponents (12,19) obtained in the linear theory are only expected to be exact when the spatial dimension d is high enough. Below a certain upper critical dimension d c , nonlinear terms become important which modify the scaling and correlation length exponents. The linear exponents (12) can however be used to gauge the importance of various nonlinearities in the EOM (2) as d is lowered.
We now turn to the full EOM of g (2), and perform a rescaling (6) with the linear exponents (12). If the spatial dimension d is large enough, all nonlinear terms are irrelevant, i.e., they vanish as → ∞. As d decreases from, say infinity, the nonlinear terms that first become relevant (and are not fine-tuned to zero), i.e., terms that diverge as → ∞, are α 2 ρ 2 g and κ 2 ρ 2 ∇ρ , which happens at the upper critical dimension d c = 6. These non-linear terms, together with the linear terms, support the following symmetry: which emerges around the MCP. One can therefore simplify the consideration by restricting to the subspace of EOM compatible with this symmetry, in which α 1 and κ 1 are vanishing and are not generated under RG transformations.
A physical interpretation of this emergent symmetry (22) corresponds to the equivalence between a highdensity band traveling in directionn and the corresponding low-density band (with the same profile but inverted) traveling in the −n direction.
Just below six dimensions, the universal hydrodynamic EOM (2) is therefore where γ is a dimensionless coefficient introduced to allow for renormalization of the temporal derivative term associated to g. Note that such a RG correction to the temporal derivative is also present in a recent study of incompressible active fluids with quenched disorder [25]. Further, we note that the signs of the nonlinear terms (with α 2 , κ 2 > 0) are chosen for the sake of stability. By the same token, the term ζ∇ 2 ∇ρ stabilizes the system in the case of κ 0 < 0. In fact, this term is marginal according to our linear theory and is therefore required in our discussion. It can be interpreted as an effective "pressure" or compressibility term for the momentum density field which, in the limit of small ζ manifests as an effective diffusion for the density mode. Its dispersion relation is given as where the last approximation is valid in the limit of ζ µ 2 /(4γ) such that the diffusion constant of the density mode is given by ζ/γ.

VI. DRG ANALYSIS
Traditionally, a DRG analysis together with theexpansion method is now applied. As we will demonstrate now however, a one-loop calculation as is usually performed, will not be sufficient.
The DRG [15] is usually performed by first transforming the EOM (1,23) to Fourier-space and then splitting the fields into small and large scale modes, arbitrarily split at an inverse length-scale Λ = Λe − , which is a fraction of the physical cutoff scale Λ that defines the smallest length scale of the system, e.g., the average distance between individual particles: One can then eliminate the small scale modes ρ > and g > from the EOM by recursively reinserting the formal solution for the small scale modes in terms of the large scale modes provided by the EOM. This generates a hierarchy of terms, which can only be truncated by assuming the interaction terms are small, i.e., by going to the perturbative limit. Averaging this expression over the small scale noise terms f > that remained so far in the equation, then generates effective contributions to the coefficients of the EOM (23) that depend on the coarse-graining scale Λ . These terms can be represented diagrammatically through Feynman diagrams which can be classified by their number of loops, i.e., number of integrals one needs to solve to determine the correction.
The conventions we will be using in this paper arẽ where id is the identity matrix and G is the propagator of the momentum density fields In general, each diagram expression is of the tensorial rank equal to the number of open unbroken lines. These, so-called, graphical corrections are where outgoing lines without a wave vector imply that the corresponding wavevector has been set to zero. The loop wavevector has not been written explicitly and its in-tegral is implied. Whether the external wavenumber q is routed through the upper or the lower part of the loop is irrelevant up to order . Also, the tildes are omitted since the external frequency ω q is set to 0 for every diagram in our approximation. Only internal lines correspond to a propagator according to Eq. (26). External lines only mark open vector indices and external wavevectors.
We further note that the different prefactors in Eq. (29) stem from the fact that the right-hand side has to be expanded to linear order in q. For diagrams whose first vertex has three outgoing density propagators (dotted lines), i.e. the vertex in Eq. (26f), it is immaterial at this order whether the external wave vector leaves before or after the loop part, e.g., For diagrams where the first vertex has two outgoing density fields and one outgoing momentum field (two dotted and one unbroken line), i.e. the vertex in Eq. (26e), contributions where the external wave vector leaves before passing through the loop part vanish identically since they are wave vector independent and the wave vector independent part vanishes due to antisymmetry of the integrand, e.g., This attributes the relative factor 2/3 between these two types of diagrams since all permutations of the outgoing wavevector need to be considered.
At the 1-loop level, we can set κ 0 and α 0 to zero on the right-hand-side of Eqns. (28) and (29). For Eq. (30), the expansion must be carried out explicitly to second order in κ 0 and α 0 however, which is most conveniently done after the frequency integration. Further details of the evaluation of the diagrams can be found in App. C.
As an example, for one of the specific 1-loop diagrams from Eq. (28) that provides a RG correction to α 2 the analytical expression is, We note here that the ζ term is indeed crucial to regularize the RG calculation, as without it the above 1-loop frequency integral is clearly divergent.
Now together with these corrections, stemming from the nonlinear terms, we can reattempt the rescaling in Eq. (6) which we include as an effective rescaling of the couplings. In addition, the EOM (23) is divided by γ, to fix the time-derivative coefficient to unity: If we now consider that is an infinitesimal number → d , defining the so-called Wilsonian momentum shell, we can write down the DRG flow equations: This means that there are three parts that contribute to the flow equation of each coupling: the rescaling of fields, lengths and time, the graphical correction of each term, which have been rescaled by the same factor as the respective coupling, and finally the graphical correction of γ from dividing the total EOM by this factor. The coupling γ itself is not rescaled and will in fact not approach a fixed point. But this does not matter since the rescaled EOM no longer depend on γ.
If the flow equations (34) are vanishing, this means that the EOM are invariant under this rescaling transformation in the presence of nonlinearities, which implies power-law correlations as discussed in Sec. III.
To facilitate the comparison between our DRG calculation and our FRG analysis to be presented later, we further define the dimensionless couplings through which all the flow equations can be expressed, where the geometric factor S d /(2π) d , with the surface area of a d-dimensional unit sphere S d = 2π d/2 /Γ(d/2) and the Euler gamma function Γ, was introduced for convenience. In this way, the rescaling introduced earlier is removed again from the flow equations, which might seem surprising. However, engineering and scaling dimension are actually very closely related. Any anomalous scaling exponent is generated by the renormalization of couplings that relate different units. For example, the diffusion constant µ /γ relates time and length scales. The four couplings carrying that information in this system are µ , ζ, γ and D which are all removed from appearing explicitly in the flow equations due to making them dimensionless (36). What remains in the flow equations, however, is their graphical corrections The dimensionless flow equations can then be written as The dimensionless graphical correctionsḠr DRG are defined analogously to their respective couplings. Now, since all nonlinearities in the EOM are cubic in nature (the α 2 and κ 2 terms), the only coefficients that receive graphical corrections are α 0,2 and κ 0,2 at the 1loop level, which are shown in Eq. (28)(29)(30). The other graphical corrections are all zero, i.e., and as a result, Therefore, we can directly infer that the scaling exponents are unchanged from the linear theory The remaining flow equations, perturbatively expanded to second order inᾱ 0 ,κ 0 ,ᾱ 2 andκ 2 are

VII. FRG ANALYSIS
Our FRG analysis is based on the so-called Wetterich equation [30][31][32]: where Γ k is the k-dependent effective average action, with k being the inverse length scale up to which fluctuations have been averaged out (i.e., k plays the equivalent role of Λ = Λe − in the previous DRG analysis.) The functional Γ k interpolates from the microscopic action Γ Λ to the macroscopic effective average action Γ 0 . As the Legendre transform of the logarithm of the partition function, it contains all statistics of the many-body problem and can therefore be regarded as its full solution. It can further be regarded as the classical action for the average of the fields. The gradual incorporation of fluctuations as k → 0 is facilitated by the "regulator" R k , which serves to suppress fluctuations of length scales greater than k −1 . The regulator can be chosen arbitrarily as long as R Λ ≈ ∞ and R 0 = 0 to ensure the correct boundary conditions for Γ k . Further, Γ k in Eq. (43) denotes the field dependent matrix of the second functional derivatives of Γ k (i.e., entries are of the form δ 2 Γ k /(δgδρ), etc), and Tr stands for the matrix trace over internal indices and integration over the internal wave vector and frequency.

A. Choice of functional in the Wetterich equation
While the Wetterich equation (43) is in principle exact, the actual implementation of the RG flow relies on restricting the functional Γ k to a manageable form. Here, we will take Γ k to be the functional obtained from the EOM (1,23) via the Martin-Siggia-Rose-de Dominicis-Janssen formalism [37,[61][62][63]: where r ≡ d d rdt, and all coefficients above (µ 1 , µ 2 , α 0 , etc.) are now k dependent. The response fields introduced by the formalism are denoted byḡ andρ. Note the absence of any coefficients in the density 'sector' of Γ k , i.e., terms proportional toρ in (44). This is due to the fact that it does not renormalize because of an extended symmetry as defined in Ref. [64]. Specifically, this extended symmetry stems from the linear nature of the continuity equation (1), which implies that, under the transformation ρ(r, t) →ρ(r, t) + ε(r, t), with an arbitrary field ε, the microscopic action Γ Λ transforms linearly in the fields, Since the transformation (45) is scale-independent it commutes with the scale derivative ∂ k and we can use the Wetterich equation (43) to see how this relation changes under RG transformations in the case of infinitesimal ε: As δΓ k is linear initially at the scale k = Λ, δΓ k vanishes, so δΓ k remains unchanged at the infinitesimally larger RG scale k = Λ + dk. This argument can be repeated at this scale and so on, showing that ∂ k δΓ k = 0 and δΓ Λ = δΓ k = δΓ 0 at all scales. Thus, the density sector does not renormalize, which is why, in Eq. (44), we have set the coefficients characterizing it to unity, fixing the engineering dimensions of the fields. This nonrenormalization further implies the hyperscaling relationship between the density and momentum density field, For the momentum density 'sector' of Γ k , i.e., terms proportional toḡ in (44), we know from our linear theory that this form of Γ k is sufficient only around the critical dimension d c = 6. As a result, we expect that the validity of our quantitative predictions is limited to around d c . Therefore, we will express our results as corrections to the linear theory in terms of = d c − d. In particular, our results for universal exponents will coincide with the perturbative DRG results to order , if the fixed point values forμ andζ are put into the DRG calculation by hand. At the same time, corrections of order O( 2 ) are expected to differ from the DRG results at the same order, since the FRG analysis is nonperturbative in nature. Even though our approach, therefore, becomes perturbative in the couplings α 0 , κ 0 , α 2 and κ 2 , since their FP values are controlled by , our approach is not fully perturbative, since we take the full dependence of the flow equations onμ andζ into account, whose FP values are not controlled by .

B. Regulator
Besides the form of the average action, the regulator R k needs to be specified, which we choose to be, in spatiotemporally Fourier transformed space, where the ordering of the matrix entries is: (ḡ, g,ρ, ρ).
The choice of a time-independent regulator is common for dynamical systems [37]. Also, this matrix form does not regulate the density sector directly but rather introduces a k dependent "pressure term" in the momentum field sector. This regularization sufficiently cuts off large and small scale fluctuations (for appropriate choices of A k and B k ) while following the overall structure of the EOM. In particular, it leaves the extended symmetry unmodified, implying that density remains conserved, even in the regulated theory, which we believe to be crucial to obtain the correct scaling behavior (e.g., the value of the dynamic exponents in dynamic Ising models depends on whether the dynamics are conservative or not [4]).
We also define the following in Eq. (49): where we write the k-dependence of the couplings explicitly and a is an arbitrary positive constant. In principle, all results obtained should be independent of the regulator choice, however, truncating the form of Γ k usually introduces some form of regulator dependence. This dependence can be judged by the a-dependence of the critical exponents. It turns out that for an algebraic regulator as in Eq. (50b), the critical exponents are independent of a. This is shown numerically below, but an analytical argument has been given in Ref. [65] as well.
We further verify our results using also another class of regulator, a generalization of the Litim regulator [66] m(y) = a(1 − y) 4 Θ(1 − y) .
The fourth order is required to ensure continuous integrands in the RG flow equations, as derivatives of m(y) up to fourth order appear.

C. FRG flow equations
With the forms of Γ k and R k defined, one can then use the Wetterich equation (43) to project a set of coupled ordinary differential equations (ODEs), one for each coefficient in the functional (44). For instance, since where V T = (2π) d+1 δ d+1 (0) is the spatio-temporal volume, we obtain from the Wetterich equation (43) that where ≡ − ln(k/Λ). The full set of such FRG flow equations correspond to the graphical corrections to the RG flow equations in our previous DRG analysis (28)(29)(30). At the same time, as in our DRG analysis (36), it is convenient to introduce dimensionless couplings, for the determination of potential RG fixed points (FPs).
In the FRG formalism, the non-dimensionalization and rescaling are performed in a single step, where the inverse scale k takes the role of e − Λ. This essentially skips the step where the flow equations are written as in Eq. (34). Specifically, we define the following: The rescaling is mostly prescribed by the dimensionality of the couplings. However, we have taken the liberty to rescaleκ 0 andκ 2 with a factor that contains the dimensionless couplingζ. This particular choice ensures that the flow equations remain regular in the case thatζ → 0.
Like in the DRG calculation, we can now again define the anomalous scaling dimensions Note that the flow equation ofζ admits two different kinds of FPs: either η ρ = 0 andζ = 0 or η ρ = 0 andζ = 0. This motivates our choice of dimensionless couplings.
The nondimensional graphical corrections to the FRG flow equations can be written in terms of the second-order functional derivatives: evaluated at vanishing fieldsḡ = g =ρ = 0, except for the density which is set to a value ρ unif , uniform in space and time, or their dimensionful equivalents, F g = µ k 2F g , F ρ = ζk 2F ρ and F D = DF D . Akin to the DRG, the F terms can be represented diagrammatically, hence we call them graphical corrections too. Together with the detailed analytical expressions, they are derived in App. A.
Before we proceed to detail what these graphical corrections are exactly, we will first discuss how FRG enables us to go beyond the 1-loop calculation in the previous DRG calculation. The strategy that we use here follows from Ref. [35], which is to evaluate the aboveF's at a non-vanishing (i.e., off-critical) density. This procedure is supported by the following physical argument: The effective average action at a nonzero scale k, Γ k , serves as an effective theory that describes subsystems of size k −1 . Within this subsystem, the mean density background can be different from the total density background which is vanishing at the MCP, ρ k = ρ 0 = 0, as mass can be exchanged between the subsystems. This is the case when a linear stability analysis on Γ k , which includes nonlinear effects on scales smaller than k −1 , predicts that the homogeneous state ρ = 0 is unstable, which happens when α 0,k < 0 and/or κ 0,k < 0 (the scale dependence has been made explicit here to emphasize that the instability is scale dependent). Then, induced by fluctuations, the system will locally phase separate and spontaneously select a new local density until local stability is reached again. The new local equilibrium is reached when α(ρ) = α 0 + α 2 ρ 2 = 0 and κ(ρ) = κ 0 + κ 2 ρ 2 > 0 or κ(ρ) = 0 and α(ρ) > 0, whichever happens first. Since this is the physical state of the system, we choose it as the constant background field value for ρ unif when evaluating the second-order functional derivatives at scale k (57) for the flow equations for µ 1 , µ and ζ, i.e., For this definition to work, the flow equations for α 0 , κ 0 , α 2 and κ 2 must be evaluated at ρ unif = 0 though, e.g., as Tr by definition.
Note however, distinct from the procedure applied to the Ising model in Ref. [35], we do not make the change of variables: α 0 → ρ unif , in the flow equations. Instead, we treat ρ unif as an auxiliary variable that is determined from Eq. (58) at each RG step. Our procedure is more advantageous here because of the following: since α 0 , κ 0 , α 2 and κ 2 are always changing smoothly in k, the flow equations never become singular, whereas, since ρ unif can have cusps, e.g., if α 0 flips its sign while κ 0 > 0, the flow equation for ρ unif has singular behavior at ρ unif = 0.
Finally, we can then write the graphical corrections to Eq. (56) as where theF's are given in Eq. (57), and we have introduced the dimensionless wavevectorq = q/k,q = |q|, and frequencyω = ωµ /(ζk 2 ). As described above, the couplings α 0 , κ 0 , α 2 and κ 2 are evaluated at vanishing background density fluctuation, while the remaining couplings and anomalous dimensions are evaluated at the value shown in Eq. (58).
We now impose the small expansion by neglecting all but the terms from leading order in by usinḡ which we know already from our DRG analysis. The flow equations, akin to Eq. (42) in the DRG analysis, can then be written as where the anomalous dimensions are The evaluation of theF's at a nonvanishing density can now be seen to serve two purposes: First, it enables the flow equations for the nonlinear couplingsᾱ 2 and κ 2 to be projected fromF g andF ρ by taking a secondorder derivative with respect to ρ unif . Secondly, if we were to set ρ = 0 for the evaluation of the flow equations ofμ andζ as well as the anomalous dimensions, they would be vanishing. Then the flow equations would be equivalent to the one-loop DRG result. (In fact, we will show in App. C that one can obtain the one-loop DRG equations exactly from the FRG formalism by using a specific "sharp" regulator, as in Ref. [67].) Choosing a nonvanishing ρ unif , therefore, allows us to incorporate effects that go beyond the one-loop level.

VIII. NOVEL RG FIXED POINTS
The fixed points of the FRG flow equations (63,64) determine the universality classes of the system and their associated scaling behavior. While the flow equations (63,64) can in principle be expressed analytically, the number of terms involved renders them unilluminating. Further, we are unable to analytically solve some of the integrals buried in the definitions of the F's (A10). We, therefore, use a combination of computer algebra and numerical methods to solve the FRG equations and thus discern the flow of these couplings upon decreasing the inverse length scale k. Details of the implementation are given in App. B.
In a typical perturbative DRG calculation to one-loop order, one would find that the flow equations for the nonlinear couplings κ 2 and α 2 decouple from the relevant couplings κ 0 and α 0 . Their FP values can therefore be easily obtained even in a numerical calculation since usually at least one FP in this subspace is attractive. In our FRG approach, however, the flow equations for the amplitude ratiosμ andζ, which the non-linear couplings depend on, are directly proportional to the relevant couplingsᾱ 0 andκ 0 through Eq. (58). Therefore, one has to solve all flow equations simultaneously. This is problematic since the relevant couplings diverge from the FP. To tackle this problem in an FRG calculation, one typically invokes the shooting method [33,35] to fine-tune the relevant parameters, which however becomes difficult when there are many parameters to fine-tune. Here we have developed the following simple method to tackle this problem.

A. Fine-tuning by reversing RG flows
To steer the couplings towards the fixed points, we invert the sign of the relevant flow equations. This operation manifestly leaves the locations of the FPs invariant, but changes their stability. The flow equations, therefore, fine-tune themselves. Once the fixed point solution is found, the original signs can be restored to obtain the critical exponents. This method can also be extended to explore other unstable FPs by inverting additional flow equations.
With the help of this simple trick, we find a total of four FPs (Fig. 4 and Tab. I). One, FP3, is stable and therefore governs generically the universal critical behavior of the MCP under consideration. It is denoted by the red circle in Fig. 4 and reached by performing the following inversions Here, "stability" refers to the stability within the "critical manifold". We also obtain two other unstable nontrivial FPs: FP2, the green square in Fig. 4 reached by, additionally to (65), inverting and FP4, the blue diamond in Fig. 4 reached by, additionally to (65), inverting Finally, there is the trivial Gaussian FP, FP1, denoted by the yellow pentagon.
To the best of our knowledge, the universality classes associated to all FPs are novel, except for the Gaussian FP (yellow pentagon).

B. Genuine nonequilibrium UCs
In equilibrium the fluctuation dissipation theorem implies η γ = η D . Since this is clearly broken for FP3 and FP4 (see Tab. I), we can conclude that FP3 and FP4 are novel nonequilibrium universality classes. While the fluctuation dissipation theorem does not seem to be broken for FP2, this does not necessarily imply that FP2 describes the critical phenomenon of an equilibrium system. We discuss this further in Sec. VIII E.

C. Nonlinear scaling
Now that we have actually found fixed points through our RG analysis, we will revisit the scaling behavior of the theory as well as the correlation functions.
At the FPs, the anomalous dimensions, i.e., the η's (55), take on universal FP values. Therefore, under a RG transformation from a reference scale k , where the system is already sufficiently close to the FP, to the scale k, the dimensionful EOM (23) in our truncation transform as where = log(k /k), primed couplings denote couplings at the reference scale k and starred couplings denote the FP value of the couplings. This EOM together with the continuity equation (1) is scale-invariant if we rescale lengths, time and fields, with the nonlinear scaling exponents: Note that this only works if η ρ = 0, which, for the nontrivial FPs, is only the case for FP2 and FP3, but not FP4. We will discuss the case η ρ = 0 for FP4 in Sec. VIII F. If η ρ = 0, we have found a rescaling transformation under which the EOM are invariant, and can, therefore, apply the argument from Sec. III to deduce the scaling of the correlation functions (8,9). But why is it even possible to extract the scaling behavior from the coefficients at finite k, i.e. in a regulated theory with a finite IR cutoff? To see this, consider the following argument (compare also to [23,68]). Suppose, we have the inverse propagator in Fourier-space sufficiently close to the FP, i.e. Γ (1,1,0,0) k (ω, q) scales homogeneously under an RG transformation, k → sk, and simultaneous rescaling of q → sq and ω → s z ω (the scaling behavior of Γ (1,1,0,0) is inverse to that of G), i.e.
for some χ. Now we can consider two equivalent cases. First we set s = k /q with a constant scale k where we can now safely take the limit k → 0, showing that the propagator of the effective action, which is the full solution to the many-body problem, all nonlinear fluctuations included, follows a powerlaw in q with the exponent χ. Secondly we set s = k /k where in the last line, we developed Γ (1,1,0,0) k to zeroth order in q and ω, corresponding to our truncation which, by neglecting higher order derivative terms, also assumes that q k. This last result shows, that the constant part of the inverse propagator scales in k with the same exponent χ. And this exponent has already been obtained through our FRG analysis, i.e., χ = η µ − 2.
A similar argument can be made for the rest of the entries is Γ (2) k , showing the scaling behavior (8,9) for the realspace propagators, obtained from inverting and then Fourier-transforming Γ (2) k .

D. Universal critical exponents and amplitude ratios
The resulting values for the critical exponents at these FPs, depending on the regulator parameter a are plotted in Fig. 5 for both regulators, Eq. (50b) and (51), exemplarily at = 0.1. This clearly shows that, for the algebraic regulator (50b), the critical exponents are independent of the parameter a. While we have shown this numerically, an analytical argument for this is given in [65]. In contrast, the results using the Litim regulator (51) do depend on a, however, the estimated values for the critical exponents are compatible with those of the algebraic regulator, i.e. the deviations are all smaller than 3 . Due to the expansion accuracy can anyways only be expected to this degree. By virtue of the principle of minimal sensitivity [69], we, therefore, chose the result of the algebraic regulator as our main results which, expressed in terms of = (d c − d), are shown in Table II. In addition to these critical exponents, we can also provide quantitative predictions on the two universal amplitude ratios:μ andζ, shown in Tab. I. These could in principle be measured experimentally from diffusion constants and the wavelength of density waves.
Finally, we can also determine the exponents y 1 and y 2 , describing the divergence of the correlation length according to Eq. (20), which are also shown in Tab. II.
The dependence of the amplitude ratiosμ andζ, as well as the correlation length exponents y 1 and y 2 on the regulator parameter a is similar to the other exponents in Fig. 5, i.e. they are independent of a for the algebraic regulator (50b) and the results obtained with the Litim regulator (51), while varying in a, are compatible with those of the algebraic regulator up to a correction of nextto-leading order in . For the second fixed point, we observe that the inverse correlation length exponent y 1 agrees with the correlation length exponent of the Ising universality class in 4 − dimensions. Though the difference in upper critical dimension between the two universality classes implies their distinctness, it is nevertheless interesting to explore how this relationship arises, which we will do in this subsection. At FP2, the FP values of all α-couplings vanish, which implies that the EOM reduce to which is linear in g. We can immediately see, that the transverse component decouples from both the density and longitudinal momentum field and that its dynamics are given by the mean field critical O(d − 1) EOM In the parallel sector, we can eliminate g from the EOM to obtain This equation is again reminiscent of Model B dynamics, except that the time-derivative term is heavily modified. The linear mode of the momentum field, though eliminated from the equation, manifests now in the second order time-derivative. Here, all anomalous dimensions are zero. As a result, the fluctuation-dissipation relation is not explicitly broken and it remains to be seen whether the model equation (76) corresponds to an equilibrium system or not.
This EOM can therefore be seen as an Ising model with exotic two-mode dynamics, that rises the scaling dimension of the field, and therefore also the upper critical dimension.

F. FP4: Emergence of two time-scales
So far, we have discussed the FPs, where η ρ = 0. Now we turn to the case η ρ = 0, which is the case for FP4. In this case, it is impossible to choose rescaling exponents z, χ g and χ ρ such that all the terms in the EOM (68) rescale homogeneously. For instance, if we were to choose the same exponents as for FP2 and 3 (70), in the large , i.e. hydrodynamic limit, the continuity equation (1) would reduce to a statement of staticality, and all the "pressure-terms", proportional to ∇ρ in (68), would vanish, decoupling the momentum density field from the density field at the linear level. Since the continuity equation is modified, this is the only FP where the hyperscaling relation χ ρ − χ g − 1 = z − 2, enforced by the extended symmetry of the continuity equation, is broken.
If we interpret this as the momentum-densitycouplings becoming irrelevant at this fixed point, we can simply omit these and obtain a scale invariance on the remaining terms. For the momentum density correlation function, we can therefore conclude that, At the linear level, they would look just like in Eq. (16d) with ζ = 0. In this decoupled limit however, the density correlation function cannot be determined since the density field seemingly decouples from the noise term f . This suggests, that this choice of rescaling is not the correct one when looking at density-density correlations. If we were to choose instead a different time and momentum field rescaling, the continuity equation remains scale invariant and the "pressure-terms" stay relevant. However, now the timederivative term on the left-hand-side of Eq. (68) vanishes in the large , i.e., hydrodynamic, limit. Now, again regarding this term as irrelevant, the remaining terms support a scale invariance, from which we can conclude that the density correlation function scales as At the linear level, with this rescaling, the correlation functions would be C g (r, t) = ζ q e iq·r 2Dω 2 P (q) TABLE III. Fixed point values and universal critical exponents for the two additional fixed points in theζ = 0 plane, expressed as an -expansion from the upper dimension dc = 6. The fixed point values of α0 and κ0 are again normalized such that they are independent of the regulator parameter a. When no value forμ is given, it can take any arbitrary value and is not universal in this case.
The transverse components of the momentum density remain unmodified. This clearly shows that at FP4 two separate time-scales are emerging; one fast time-scale at which the momentum density can react quickly to perturbations in the presence of a density field frozen in time, and a slow time-scale over which the density relaxes once momentum fluctuations have already dissipated long ago. In other words, the momentum density ceases to be a hydrodynamic variable at this fixed point, since it becomes a "fast" mode.
Even though our analysis has shown that in the vicinity of FP4 the "pressure-terms" proportional to ζ are irrelevant, they cannot be neglected in the RG analysis. First, as discussed in Sec. VI if ζ is naively set to zero, the frequency integrals are clearly divergent. Secondly, the FP value of κ 2 plays an important role in determining the critical exponents of FP4. This coupling must therefore be included necessarily in the discussion.
The right coupling to neglect instead is therefore the time-derivative term of the momentum density. Then, its EOM becomes an exact, time-invariant identity enslaving the momentum field to the density field. It can be used to eliminate the momentum field from the continuity equation to obtain, for simplicity at the linear level, that which, except for the µ term, is exactly the general equation of model B dynamics one would write down for a conserved, scalar quantity [4]. Fine-tuning κ 0 → 0 alone yields the critical model B, Ising universality class. We can therefore interpret FP4 as a genuine nonequilibrium multicritical point of model B, where, in addition to κ 0 , also the time-derivative term, characterized by α 0 , is finetuned to zero. For reasons of stability, the higher-order derivative term characterized by µ , is needed.
The ζ = 0 manifold, therefore, describes the model B theory subspace. Within it, in addition to FP4, we find two additional FPs which are related to FP2 and FP3, so we name them FP2' and FP3' accordingly. Since they are unstable in the ζ direction and already described by the model B EOM, we will not discuss them in detail.
Instead, we will simply report here the FP values and critical exponents in Tab. III.

IX. SUMMARY & OUTLOOK
We have demonstrated in this work the whole process from first formulating the large scale hydrodynamic equations to discovering novel universality classes on the particular problem of the multicritical point in the phase diagram of compressible active matter, where the critical point of the flocking transition coincides with the critical point of disordered phase-separation.
We started with the Toner-Tu equations, derived purely from symmetry arguments and conservation laws, which therefore describe general compressible polar active fluids, i.e. systems of self-propelling and aligning particles at the hydrodynamic level. Using mean-fieldtheory and linear stability analysis one can show, first, the existence of two homogeneous phases and, secondly, the existence of regions where the homogeneous state is unstable indicating phase-separation, revealing the phase diagram of compressible polar active matter. Linear stability analysis further reveals the existence of two critical points, the critical order-disorder transition and the critical point of disordered phase-separation. Further, there is a multicritical point in the phase diagram where both these critical points overlap.
Analyzing the linearized equations of motion around this multicritical point, then revealed the scale invariance of the theory at the linear level, which results in powerlaw scaling of the correlation functions. We elucidated how these functions, and thereby the scaling exponents as well as the scaling behavior of the correlation length, could be measured in principle.
The scaling behavior of the linear theory further informed us about the relevance of the possible nonlinear terms and the critical dimension, above which the linear theory is expected to be exact. While in the nonlinear regime, the correlation functions could no longer be calculated explicitly, scale invariance will in general still lead to powerlaw correlation functions, though their scaling exponents are no longer trivially determined.
To determine these nontrivial exponents, we attempted a one-loop dynamic renormalization group (DRG) approach within the -expansion, which ultimately failed, since two-loop effects are necessary to capture the universal physics of this problem.
To take these two-loop effects into account, we set up a functional renormalization group (FRG) ansatz which describes effective theories at intermediate scales and whose scale dependence is described by the Wetterich equation. A key ingredient to this approach was the physical insight that subsystems at an intermediate scale are locally in a homogeneous state even if the total system is phase separating. This enabled our FRG approach to go beyond the perturbative DRG approach.
Using various computer algebra and numerical methods, the renormalization group flow equations could then be evaluated. Further, by inverting the sign of the flow equation for the relevant couplings, we found three renormalization group fixed points. Finding them proves the existence of three novel universality classes, at least two of them being genuinely out of equilibrium, since they allow for nontrivial scaling of the correlation functions which we are able to predict quantitatively. This is the first time that novel universality classes have been shown to exist for compressible polar active matter and it is the first time that the FRG has been used in the field of active matter to uncover new critical phenomena beyond the one-loop level of the perturbative DRG.
Possible future directions are using the FRG to explore open questions in compressible active matter such as: what are the universality classes of the critical orderdisorder transition [58] and of the ordered phase [56,70].

Appendix A: FRG Flow equations
In the following we give details on how the right-hand-side of Eqns. (57) can be obtained. As a first step, to reduce some of the complexity of the following calculation, we use a trick commonly used in FRG. We reformulate the Wetterich equation (43) to the form where the derivative ∂ = k ∂ k on the right-hand-side, only acts on the k-dependence of the regulator. This reduces the number of propagators (Γ (2) k + R k ) −1 in the following expressions by one. As the projections in Eqns. (57) involve two functional derivatives of Eq. (A1), the right-hand-side of Eqns. (57) after derivation will involve, in addition to Γ (2) k also the third and fourth order functional derivatives of Γ k : where we didn't write explicitly the dependence of Γ k and Γ (4) k , on the external wavevectors and frequencies that come from the functional derivatives. As written in Eqns. (57), the expressions are then to be evaluated at the constant field ρ(r, t) = ρ unif . The expressions for Γ (2) k , Γ k and Γ (4) k with this background field are straightforwardly determined from our ansatz, Eq. (44): where we introduced the short-hand notationδ qp... = (2π) d+1 δ d+1 (q +p + . . . ). All other functional derivatives of Γ k , up to fourth order, evaluated at these background fields and not mentioned here are vanishing.
The propagator G is essentially the same as in the DRG calculation (27), except that it now contains the terms introduced by the regulator, which make the flow equations IR convergent, where in the DRG calculation this is taken care of by the integral boundaries.
If we regard the individual entries of Eq. (A2), especially those corresponding to the F's in Eq. (57), we can arrange the vertices, Eqns. (A3f-A3i), in matrices similar to that in Eq. (A4), such that the matrix products in Eq. (A2) can be carried out and traced over. Then the F's can be decomposed into individual terms, which can then be represented diagrammatically, In dimensionless units,Ḡ   (63) is then solved using a fourth order adaptative Runge-Kutta-Fehlberg (4,5) algorithm provided by the GNU Scientific Library [71], where at each RG "time-step", the wave number integral (B6) in (A10) is solved numerically using an adaptative quadrature routine, with the 15 point Gauss-Kronrod rule for infinite boundary integrals in case of the algebraic regulator and the 61 point Gauss-Kronrod rule for finite boundary integrals in case of the Litim regulator, again provided by the GNU Scientific Library [71]. For the adaptative ODE-solver we set a maximum relative and absolute error of e r = 10 −11 and e a = 10 −13 respectively. For the integration they are set to e r = 10 −12 and e a = 10 −10 . Since, we are modifying the flow equations such that each fixed point we would like to investigate is attractive, the initial conditions do not matter too much (as long as they are within the attractive basin). In our analysis we set them close to the Gaussian FP. Since we obtained all our results in Tab. II at = 0.1, the following initial conditions are sufficiently close to the Gaussian FP:μ =ζ = 1,ᾱ 2 =κ 2 = 10 −3 andᾱ 0 =κ 0 = 0. We use these initial conditions to analyze all FPs, except for the following modifications: if we are interested in FP2, we setᾱ 2 = 0 (with this condition it is not necessary to invert the flow equation forᾱ 2 , as it is not being generated if it is vanishing initially), and when we investigate FP4, we setκ 2 = −10 −3 .
When solving the ODEs we set an initial "time-step" of d = 5 × 10 −4 initially, but it is being quickly changed by the adaptative algorithm. The adaptative nature of the algorithm is key to solve this problem, since there are three vastly different RG "time-scales", given by the critical exponents, to resolve. First, the relevant couplings α 0 and κ 0 converge very quickly with a time-scale of ∼ 0.5 to the critical surface. The next time-scale is that of nonlinear couplings α 2 and κ 2 and of order ∼ 1/ . Then the final and longest time-scale sets in, wherein the amplitude ratiosμ andζ converge. This time-scale is approximately of the size of the anomalous dimensions, i.e., of order ∼ 1/ 2 . However, since the prefactor of the O( 2 ) critical exponents is typically a lot smaller than unity, we find total convergence of the fixed points only at about ∼ 5 × 10 4 , though FP4 converges somewhat quicker, since the scaling exponents are generally larger, compare Fig. 5 and Tab. II.
Finally, once the ODEs have converged to the desired FP, the correlation length exponents can be obtained by first restoring the original signs of the flow equations, and then taking discrete derivatives with a finite difference of dc = 10 −7 with respect to the couplings. The largest two eigenvalues of the so obtained matrix characterizing the linearized flow equations around the FP give the correlation length exponents y 1 and y 2 .
The complete algorithm is summarized again in Alg. 1.