Normal form for renormalization groups

The results of the renormalization group are commonly advertised as the existence of power law singularities near critical points. The classic predictions are often violated and logarithmic and exponential corrections are treated on a case-by-case basis. We use the mathematics of normal form theory to systematically group these into universality families of seemingly unrelated systems united by common scaling variables. We recover and explain the existing literature and predict the nonlinear generalization for the universal homogeneous scaling functions. We show that this procedure leads to a better handling of the singularity even in classic cases and elaborate our framework using several examples.

The results of the renormalization group are commonly advertised as the existence of power law singularities near critical points. The classic predictions are often violated and logarithmic and exponential corrections are treated on a case-by-case basis. We use the mathematics of normal form theory to systematically group these into universality families of seemingly unrelated systems united by common scaling variables. We recover and explain the existing literature and predict the nonlinear generalization for the universal homogeneous scaling functions. We show that this procedure leads to a better handling of the singularity even in classic cases and elaborate our framework using several examples.

I. INTRODUCTION
Emergent scale invariance is a key to many of our current scientific and engineering challenges, including cell membranes [1], turbulence [2], fracture and plasticity [3,4], and also the more traditional continuous thermodynamic phase transitions. The current formulation of the field has an elegant framework which can explain observables that scale as power laws times homogeneous functions. However, the literature on corrections to this result, including logarithms and exponentially diverging quantities, is much more scattered and does not have a similarly systematic framework.
The renormalization group (RG) is our tool for understanding emergent scale invariance. At root, despite challenges of implementation, the renormalization group (RG) coarse grains and rescales the system to generate ordinary differential equations (ODEs) for model parameters as a function of the observed log length scale . A fixed point of these flows represents a system which looks the same at different length scales; systems near criticality flow near to this fixed point. In cases where the flow can be linearized around the fixed point, the RG implies that observables near criticality are given by a power law times a universal function of an invariant combinations of variables; e.g. the Ising model has magnetization m ∼ t β M(Lt ν ) where L is the system size and t = (T − T c )/T c is the deviation of the temperature T from the critical temperature T c .
Here we use the fact that the predictions of the RG can be written down as a set of differential equations in the abstract space of Hamiltonians. This allows us to apply a branch of dynamical systems theory, normal form theory [20,21] to provide a unified description applicable to all of these systems. We arrange these systems into universality families of theories, each defined by its normal form. Each family has universal terms (linear and nonlinear), whose values determine a system's universality class within the family. Finally, each family's normal form predicts the natural invariant scaling combinations governing universal scaling functions.
The perspective we present here is transformative: unifying, simplifying, and systematizing a previously technical subject and promising new developments in the field. Our best analogy is to the introduction of homotopy theory in the 1970's [22][23][24][25] as a systematic method that unified the treatment of some of the many defect structures studied in materials and field theories. Just as there have been several previous works that correctly identified the universal effects of nonlinear terms for phase transitions where an analytic RG approach is available [16,[26][27][28][29][30][31][32][33], the Burger's vector, winding number, and wrapping number of dislocations in crystals, disclinations in liquid crystals, and Skyrmions in nuclei were understood individually before the mathematics of homotopy theory was seen as the natural mathematical framework. Just as homotopy theory facilitated the study of defects in more complex systems (metallic glasses, cosmic strings, quasicrystals), so our normal form methods are allowing the correct identification and characterization of the singularity in systems in experimental and numerical explorations where analytic RG calculations do not yet exist [34]. Finally, homotopy theory quickly uncovered the fascinating entanglement and transformation properties of non-abelian defects [24], with early speculative applications in glass physics [35] and eventually inspiring the closely related nonabelian braiding being developed for topological quantum computing. Similarly, we demonstrate here that our methods allow, for what appears to be the first time, the use of the correct, remarkably rich, invariant scaling variables in the universal scaling functions for systems where universal nonlinear RG terms are needed, and we have discussed elsewhere [36] how our methods can be powerful tools for systematically incorporating corrections to scaling near critical points even when universal nonlinear terms are not needed. For example, in the future the normal-form change of variables we introduce here could become an inner expansion matched to series and virial expansions at extremes of the phase diagram; this would allow rapid and accurate convergent characterizations of materials systems close to and far from criticality.
Our machinery provides a straightforward method to determine the complete form of the critical singularity in these challenging cases. Our initial results are complex and interesting; they pose challenges which we propose to address in future work. The coordinate transformation to the normal form embodies analytic corrections to scaling, which allow us to address experimental systems as they vary farther from the critical point. Finally, bifurcation theory is designed to analyze low-dimensional dynamical systems without detailed understanding of the underlying equations; our methods should improve scaling collapses in critical phenomena like 2-d jamming [37] where there is numerical evidence for logarithms but no RG framework is available.
We begin by distinguishing our work from previous literature connecting the RG to normal form theory. The previous approach [38][39][40] compared the application of RG-like methods and normal form theory to solving nonlinear differential equations using perturbation theory. The connection we are making is different. We are applying normal form theory to the RG flow equations. Hence, our approach is to apply normal form theory to make predictions about the general structure of the flows given the topology (nature and number of fixed points), rather than to apply it to the model that produces these flows.
We give an introduction to normal form theory in Section II. We give a survey of the previous literature on nonlinear scaling in the RG in section III. We show how the application of normal form theory allows us to define universality families of fixed points in Section IV. We present several worked out examples starting with the 4d Ising model in Section V A and the Random Field Ising model in Section V B. We then work out the application of normal form theory to the Ising model in dimensions 1, 2 and 3 in Sections V D-V F.

II. NORMAL FORM THEORY
Normal form theory [21] is a technique to reduce differential equations to a 'normal form' by change of coordinates, often the simplest possible form. This is achieved by making near-identity coordinate transformations to get rid of as many terms as possible from the equation. It was developed initially by Poincaré to integrate nonlinear systems [41,42]. The physical behavior should be invariant under analytic changes of coordinates, and the length (or time) parameter should stay the same, which the mathematical literature addresses by perturbative polynomial changes of coordinates (attempting removal of nth order nonlinearities in the flow by using nth order or lower terms in the change of variables). To any finite order this gives an analytic change of coordinates, but it is not in general guaranteed to converge to an analytic transformation; we will thus call it a polynomial change of coordinates.
We give a brief introduction to normal form theory here for completeness. A more detailed treatment can be found in Ref. [21]. Typically one starts with a set of differential equations of the form where is some parameter, θ = {θ i } is the vector of state variables and the vector field g defines the flow. In the context of statistical mechanics and renormalization group flows, θ i 's are parameters or coupling constants that enter into the free energy and is the difference in dimension from the lower or upper critical dimensions. Let us first work with the case where does not enter into the equations. The first step is to find the fixed point of the equation and use translations to set the fixed point of each θ i is at 0. The next step is to linearize about the fixed point and reduce the linear part to the simplest possible form. In general, this is the Jordan canonical form but it is often just the eigenbasis. Then, the equation is where J is the linearized matrix of the flow and the remaining terms are in the vector field f (θ) ∼ O(θ 2 ). Terms of order k are defined to be made up of homogeneous polynomials of order k. So for θ = (θ 1 , θ 2 , θ 3 ), θ 2 1 θ 2 θ 3 ∼ O(θ 4 ). We will denote terms of order k by a lower index. So Note the index is giving the order of the polynomial and not enumerating the components of the vector field. Let the lowest non-zero term be at some order k ≥ 2 (usually 2). Then we can write The idea is to try and remove higher order terms by making coordinate changes. To remove the term f k , we try to do a coordinate change of order k, where h k ((θ)) is a polynomial inθ. This construction is similar to nonlinear scaling fields [26,43] which try to linearize the RG flow equations with a subtle difference that we will remark on later. The higher order terms which we can remove by coordinate changes correspond to analytic corrections to scaling. Then, to find the equations in the new variables.
Dh k is the matrix of partial derivatives of the vector field h k with respect to the parameters θ. Now, substituting this into the equation which upon simplification gives For the last line, notice that the matrix J is the same as DJθ (i.e. the same as the matrix of partial derivatives with respect to parametersθ of the vector Jθ).
Two of the terms can be written as the Lie bracket (a commutator for vector fields) defined as [h k , Jθ] = −((Dh k )Jθ − (DJθ)h k ) to give the final equation So, if we want to remove the term f k , we need to solve the equation [h k , Jθ] = −f k for h k . It's important to note that whether this equation can be solved or not depends only on the linear part of the equation given by the matrix J. That is, within the space of transformations that we are considering, the linear part of the equation completely determines how much the equation can be simplified and how many terms can be removed. This is not true if there are zero eigenvalues and one then has to consider a broader space of transformations which we will consider later. To see when the equation can be solved, we first note that the space of homogeneous polynomials is a vector space with a basis constructed in the obvious way θ α1 1 ...θ αn n . Any term at order k can be written as a sum of such terms for which i α i = k. The Lie bracket can be thought of as a linear operator on this space. To find the set of possible solutions is to find the range of this linear operator. Let us take the case where the linear part is diagonalizable and so just consists of the eigenvalues λ i . Let us say for simplicity that the j component of the vector f k (f k ) j = cθ α1 1 ...θ αn n for some set of {α i }. Then, the jth component of the matrix equation reduces to This can be solved easily by choosing (h k ) j = aθ α1 1 ...θ αn n and a = c When all nonlinear terms can be removed by such a coordinate transformation, then the usual case of power law scaling is obtained. The fixed point, in this case, is called hyperbolic. Alternatively, if we have a term with λ j = i λ i α i (a linear combination of other eigenvalues λ i with positive integer coefficients α i ), this term is called a resonance and cannot be removed from the equation for dθ j /d . This contributes to the singularity at the fixed point which is no longer given by power law combinations.

A. Bifurcations
Notice a special case of these equations when for some k, a particular λ k = 0. In this case, it is possible to get an infinite number of resonances because the equation λ i = λ i + α k λ k is also true for all α k and λ i . This case, when one of the eigenvalues goes to 0 depending on some parameter is called a bifurcation. If all linear eigenvalues λ i of the flows are distinct and non-zero, which terms can be removed using polynomial coordinate changes depends only on these λ i . As we saw, this approach can be formulated elegantly as a linear algebra problem of the Lie bracket on the space of homogeneous polynomials. For more general cases-including bifurcations-'hypernormal form' [44][45][46] theory develops a systematic but somewhat more brute-force machinery to identify which terms can and cannot be removed perturbatively by polynomial changes of coordinates. Classic bifurcations include the pitchfork bifurcation, the transcritical bifurcation, the saddle node and the Hopf bifurcation [47].
Confusingly, bifurcation theory separately has its own 'normal form' of bifurcations. These normal forms are derived in a very different way using the implicit function theorem. The basic idea is to ask for the smallest number of terms in the equation which will preserve the qualitative behavior of the fixed points (e.g. exchange of stability of fixed points), and then map any other equation on to this simple equation using the implicit function theorem. This mapping allows for too broad a class of transformations to be useful for our purposes. An important feature of the flows that we want to preserve is their analyticity, we therefore only consider polynomial changes of coordinates.
An explicit example is given by the 4-d Ising model. It is known that the magnetization M ∼ t 1/2 (log t) 1/3 with log log corrections. The quartic coupling u and the temperature t have flow equations which traditional bi-furcation theory would simplify to Calculating the magnetization with this set of flow equations leads to the wrong power of logarithmic corrections. By allowing too broad a class of coordinate transformations, bifurcation theory hides the true singularity in the non-analytic coordinate change. We will show that normal form theory instead predicts which does predict the correct behavior. We will present the explicit solution of this equation in Section V A. Here, we just note that the traditional log and log log terms follow from the solution's asymptotic behavior. To get these equations, we will remove higher order terms in u by using a coordinate change that is lower in order (broadening the formalism we considered in Section II). Using lower order terms to remove higher order terms is part of hypernormal form theory. For our purposes, the distinction is somewhat artificial and here we simply use normal form theory to denote any procedure that uses only polynomial changes of coordinates to change terms in flow equations. In Sections V A and V B, we will explicitly work out the case of a single variable undergoing a bifurcation for the 4d Ising model and the 2d Random Field Ising model and show how there are only a finite number of terms which cannot be changed or removed. It is worth mentioning here that there can be cases in which two variables simultaneously have 0 eigenvalues. The XY model [48] offers an example where this happens. The dATG transition in 6 dimensions has two variables that simultaneously go through a transcritical bifurcation [49,50]. Polynomial changes of coordinates in both variables can be used here too, but because there are generically more terms at higher order than at lower order (there are many more ways to combine two variables into a sixth order polynomial than there are to combine them into a third order polynomial), we usually do not have enough freedom to remove all terms. Therefore, simultaneous bifurcations in more than one variable often have an infinite number of terms in their flow equations that cannot be removed.

III. EARLIER WORK
The approach we take is inspired by Wegner's early work [8,26], subsequent developments by Aharony and Fisher [51,52], and by studies of Barma and Fisher on logarithmic corrections to scaling [31,32]. The approach of Salas and Sokal on the 2-d Potts model [16], and of Hasenbusch [33] et al. on the 2d Edwards-Anderson model is similar in spirit to ours.
Wegner [26] first constructed nonlinear scaling fields which transform linearly under an arbitrary renormalization group. His construction is very similar to the coordinate changes we considered above for normal form theory. The one difference is that Wegner explicitly allows the new coordinates to depend on the coarse graining length . We will not allow this explicit dependence on in our change of coordinates, as it doesn't seem to offer any advantage over regular normal form theory.
Eventually, the goal of using normal form theory to understand the differential equations that describe RG flow is to simplify and systematize scaling collapses. This requires a systematic way of dealing with corrections to scaling beyond the usual power laws. There are three different types of corrections to scaling that have appeared in the literature. These include logarithmic, singular and analytic corrections to scaling. Logarithms in the scaling behavior typically occur at an upper critical dimension or in the presence of a resonance. Wegner and Riedel [8] considered the case of a zero eigenvalues which occurs at the upper critical dimension of Ising and tricritical Ising models. They derived the form of the scaling in terms of logarithmic corrections to scaling. However, they used perturbation theory to ignore higher order terms in the flow equations rather than only keeping those terms which cannot be removed by an application of normal form theory. Here, we will solve the full flow equations and see that the logarithmic corrections to scaling are better incorporated as part of the true singularity using normal form theory.
Analytic corrections to scaling were explored by Aharony and Fisher [52] who gave a physical interpretation of the nonlinear scaling fields (see below Eq.( 5)) in terms of analytic corrections to scaling in the Ising model. Analytic corrections to scaling capture the difference between the physical variable T and H (that your thermometer or gaussmeter measures) and the symbolst andh in the theory of the magnet. The liquid gas transition is in the Ising universality class but a theory of the liquid gas transition has to include analytic corrections to scaling to match with the universal predictions of the Ising model. Moreover, such corrections are also needed to explain the non-universal behavior away from the fixed point. Analytic corrections to scaling will correspond to terms in the differential equations that can be removed by coordinate changes.
The singular corrections to scaling are also incorporated as part of the true singularity with the addition of irrelevant variables. Finally, the ability to change the renormalization scheme leads to what are called redundant variables. In related work [36], we argue that these variables can be seen as a gauge choice which contributes to the corrections of scaling. In forthcoming work [53], we will explore the consequence of this distinction between gauge corrections and genuine singular corrections to scaling further. Despite similar inclinations, none of these works make the complete connection to normal form theory. One advantage of our approach is precisely that it brings together this disparate literature into a unified framework. The analysis presented here is general and applicable to all kinds of situations, ranging from old problems like the nonequilibrium random field Ising model (NERFIM) [54], to newer research problems like jamming [37].

IV. UNIVERSALITY FAMILIES
Traditionally, the RG contains the concept of a universality class. The universality class is essentially determined by the critical exponents which explain the scaling behavior of a model, i.e. by linearized RG eigenvalues. Normal form theory suggests another possible classification. Each fixed point can be classified by the bifurcation or resonance that it is at. The simplest case, which is also the traditional one, is the hyperbolic universality family. In the hyperbolic case, it is possible to remove all nonlinear terms in the flow equations by changes of coordinates. Hence, the RG can be written as a linear flow to all orders in perturbation theory. Different values for the linear eigenvalues correspond to different universality classes. While traditionally this is a statement about the linearization of the RG, here it is a statement about the only terms in the flow equations that are universal in the sense that they can not be removed by a coordinate change.
The need for this generalization becomes clear when we examine cases which are not traditional. In Table I we present common universality families and well-studied statistical mechanics systems governed by each. The pitchfork bifurcation shows up in the 2-d Random Field Ising model; it has a cubic term in the equations for w, the ratio of the disorder to the coupling [55]. We have derived that the correct equations require an additional w 5 term [34], which was not included in previous work. The 2-d Ising model has a well known logarithmic correction to the specific heat, which Wegner associated with a t 2 resonance term in the flow equation [26]. The 1-d and 4-d Ising models have transcritical bifurcations. The 1-d Ising case is somewhat special and we will cover it later in Section V E. These cover all the important bifurcations with one variable [56].
When more than one variable is undergoing a bifurcation, or if more than one variable has an inherently non-linear flow, the analysis becomes considerably more complicated. This is evidenced in the the 2-d XY model at the Kosterlitz-Thouless (KT) transition [57]. It has been shown that the simplest normal form of its flow equations (in the inverse-temperature-like variable x ∼ 1/T − 1/T c and the fugacity y) has an infinite number of universal terms, which can be rearranged into an analytic function f [30] (Table I). We conjecture that the very similar transition observed in randomly grown networks [58,59] is not in the KT-universality class, but rather is in the same universality family. It is not to be expected that a percolation transition for infinite-dimensional networks should flow to the same fixed point as a 2-d magnetic system, but it is entirely plausible that they share the same normal form with a different universal function f . Different universality classes within the same universality family, such as those of the 4-d Ising model and the 2-d NERFIM have different power laws and scaling functions. However, as shown in Table I, because they both have a transcritical bifurcation the two classes have the same complicated invariant scaling combinations [60]. This hidden connection is made apparent in the shared normal form, where the quartic coupling and temperature (u, T ) in the first class are associated with the (negative of) disorder strength and avalanche size (−w, S) in the second [61].
Indeed, the normal form not only unites these universality classes, but allows a more precise handling of their singularity. It is usually stated that the magnetization M ∼ t 1/2 (log t) 1/3 , the specific heat C ∼ (log t) 1/3 and the susceptibility χ ∼ (log t) 1/3 /t with log log corrections [8]. We show in the supplementary material that the true singularity of the magnetization at the critical point is M ∼ t 1/2 W (xt −27/25 ) 1/3 , where W is the Lambert-W function defined by W (z)e W (z) = z, and x(u) is a complicated but explicit function of the irrelevant variable u. (The traditional log and log-log terms follow from the asymptotic behaviors of W (x) at large and small x. The universal power 27/25 becomes manifest in the complete singularity, but is disguised into a constant factor up to leading logs.) We now show how to apply normal form theory to specific examples.

V. APPLICATION TO SPECIFIC SYSTEMS
In the sections below, we derive in detail the scaling form for the entries shown in Table I [50] .

A. 4-d Ising
The study of critical points using the renormalization group was turned into a dynamical system problem by Wilson [62]. These RG calculations are done by first expressing the Ising model as a field theory with a quartic potential uφ 4 . Then by coarse-graining in momentum space and rescaling, one obtains the flow equations dt/d =2t −Ātu +Ctu 2 +Ētu 3 +Ḡtu 4 +Ītu 5 +Ktu 6 ..., du/d = u −Bu 2 +Du 3 +F u 4 +Hu 5 +Ju 6 +Lu 6 ..., where t is the temperature, f is the free energy and u is the leading irrelevant variable. This is the highest order This gives the equations up to order u 4 , dt/d =2t −Ātũ + (−Āb 1 + a 1B +C)tu 2 + ..., (21) dũ/d = −Bũ 2 +Dũ 3 Note that any term of the form u m t in the equations for dt/d and any term of the form u m in the equations for du/d is a resonance. Hence, the coefficientsĀ,B andD remain unchanged with this change of variables. However, the coefficientsC andĒ are changed (though the change is independent of a 2 and b 3 because they are resonances) and in particular, can be set to 0 by an appropriate choice of coefficients. This creates a general procedure for reducing this flow to its simplest possible form. First, all terms that are not resonances are removed in the usual way by solving Eq. ( 11). Then, we perturbatively remove most of the resonances using the following procedure. First consider the u flow. Suppose the lowest order term in the flow after the u 3 term is u n , i.e.
with n > 3. Consider a change of variables of the form u =ũ + b n−2ũ n−1 . Then Evidently, the coefficient of theũ n term can be set to 0 with an appropriate choice b n−2 = N n /(B(3 − n)). So all terms of the form u n with n > 3 can be removed by a change of coordinates. Incidentally, this derivation also shows why it is not possible to remove the u 3 term. Now consider the t equation with We consider a change of coordinates t =t + a n−2tũ n−2 .
Any term which is not of this form can be removed in the usual way by solving Eq. 11. Finally, we have another degree of freedom that we have used. We can rescale u and t to set some of the nonlinear coefficients to 1. This reflects the fact that the original coefficientsĀ,D depend on an arbitrary scale of u and t that we have chosen. By choosing the scale soB = 1 and the coefficient of the resonance is -1, defines D =D/B 2 and A =Ā/B. Hence, by considering all such polynomial change of coordinates, we can reduce this set of equations to their normal form The resultant equations then have 2 parameters A and D which are universal, in a way that is similar to the eigenvalues of the RG flows as in Table I. The normal form variablest,ũ,f are equal to the physical variables t, u and f to linear order (up to a rescaling). Corrections to these are analytic corrections to scaling. Hence, we will henceforth simply refer to the normal form variables as t, u and f . It is important to note that we are making a particular choice for the analytic corrections to scaling by setting them equal to zero. It is possible to make a different choice for the higher order coefficients. In particular, the equation for du/d goes to ∞ at finite if u starts at a large enough value. Hence, it may be more useful to make a different choice for the higher order coefficients. All of these choices will agree close to the critical point but will have different behavior away form the critical point. Later, we will consider a different choice for the higher order terms.
The 4-d Ising model has both a bifurcation and a resonance. The u 2 , u 3 and Aut terms come from the bifurcation and cannot be removed by an analytic change of coordinates. The t 2 term is a consequence of an integer resonance between the temperature and free energy eigenvalue, λ t = 1/ν = 2, λ f = d = 4.
Before examining the full solution of Eqs. (30 -32), we will first study the effect of each part of the RG flows. First, considering only the linear terms and coarsegraining until t( * ) = 1, the free energy is given by f ∼ t 2 . This is the mean-field result and also the traditional scaling form that RG results take in the absence of nonlinear terms in the flow equations. Second, we include the resonance between the temperature and free energy eigenvalue, which leads to an irremovable t 2 term in the flow equation for the free energy. This term cannot be removed by analytic coordinate changes, and yields a log correction to the specific heat. Third, the irrelevant variable u undergoes a transcritical bifurcation. Results in the hyper-normal form theory literature, as well as some articles in the high-energy theory literature [28,29] recognize that the simplest form that the equation can be brought into is Eq. 31. The solutions  Let us use this to derive the finite-size scaling form of the free energy. Early finite-size scaling work [65][66][67] attempted scaling collapses with logs; recent work does not attempt collapses at all [68]. Finite-size scaling requires an equation for the magnetic field, h, given by dh/d = 3h. Explicit calculations show that the coefficient of the hu term is zero (see below). The free energy is then a function of three scaling variables, u( ), t( ) and h( ). It is given by To get a finite size scaling form, we coarse grain until = log L, the system size. Note that u(L) cannot just be ignored because it is a dangerous irrelevant variable. However we can account for it by taking the combination t(L)/(u(L)) 1/2 and h(L)/(u(L)) 1/4 as our scaling variables [69]. The scaling form of the free energy then depends on u 0 which we do not have a way to change or set in the simulation. Instead, we treat u 0 as a fit parameter in the scaling form of the susceptibility: (34) At the critical point t = 0, the function Φ must be analytic for finite L (since non-analyticity requires an infinite system size). Φ(0) is therefore a constant independent of L and u 0 at t = 0. Using this, u 0 may be estimated from χ at different values of L by fitting to its predicted dependence χ ∝ L 2 (W (y[u 0 ]L 1/D ) + 1) 1/2 where y[u 0 ] is defined above. Figures 1-2 shows the scaling collapse of the magnetization and susceptibility. The magnetization is collapsed using the best-fit value of u 0 = 0.4. Though our collapses are not significantly better than the traditional logarithmic forms, the correct form of the singularity will be more apparent at larger values of u 0 . This is because the log log term which is the second term in the asymptotic expansion of the W function is very small compared to the log except at large u 0 and small L. Changing the value of u 0 will require a model different from the nearest neighbor square lattice Ising model.
So far, we have been considering the effects of changing coordinates in the control variables on the predictions of the theory. Wegner [70] had also considered changing coordinates in the degrees of freedom of the theory. These changes lead to 'redundant' variables, the corrections from which can be removed by coordinate changes. We discuss them in separate work. Here we merely note that they can be used to explain some features of the scaling, like the fact that the coefficient of the hu term is zero.

Choice of normal form
There are certain choices we have made in our application of normal form theory. One is to keep the flow parameter unchanged. Some of the dynamical systems literature considers changing to depend on other parameters. This would be unusual since the coarse graining length would depend on the physical parameters but does not seem to be disallowed. We show in the supplementary material that this does not change the predictions for the 4-d Ising model.
Normal form theory makes a particular choice for what to do with the coefficients that can be changed by coordinate changes: it sets them equal to zero. In general, however, it is not clear that the best choice to make is to set them equal to zero. Consider the equation which, as we saw, has the solution u( ) = 1/(D(1 + W (ye /D ))). Here, y = (1/(Du 0 ) − 1) exp(1/(Du 0 ) − 1). Note that u 0 > 0 as a requirement for the stability of the free energy. If u 0 < 1/D, then y > 0, and if u 0 ≤ 1/D, y ≤ 0. Hence, the domain of attraction of the fixed point at u 0 = 0 has a length 1/D. If we have a system where u 0 > 1/D, then this will lead to u( ) → ∞ in a finite coarse graining length. This is reflected in the branch cut of the W function at −1/e. In the context of high energy physics, some have tried to find deep meaning in this pole [29]. However, for scaling purposes, we generally prefer a choice of coordinates for which there is no such unphysical behavior. One natural choice is to instead use the equations For small u, this has the same behavior as Eq. 35. However, the behavior at large u is now well behaved. The solution of this equation is which is in fact somewhat simpler that the solution to Eq. ( 35). Scaling collapses with this choice of normal form for the susceptibility are shown in the case of the Ising model in Figure 3. Better numerics are needed to tell if this choice of normal form is really useful. It turns out that this form has been implicitly used before in the Random Field Ising model as we show explicitly in the next section.

B. Random Field Ising model
Finding critical exponents for the random field Ising model has been a longstanding challenge in physics. Some initial results used supersymmetry to prove an equivalence of the Random Field Ising model in dimensions d + 2 with the Ising model in dimensions d [71,72]. It was later shown that the lower critical dimension of the Random Field Ising model is not 3 (as would be expected from such a correspondence) but rather 2 [73]. The upper critical dimension is 6. Here, we will look at the scaling behavior of the Random Field Ising model at its lower critical dimension, d = 2. Consider a spin system with a random field.
where, J is the nearest neighbor coupling and h i is a random field chosen from a Gaussian distribution with width r. A phenomenological theory for the RG was formulated by Bray and Moore [55]. It turns out to be useful to define a quantity w = r/J. Then, using heuristic arguments on the stability of domain walls, they derive By considering the cost of roughening the domain wall because of the presence of random fields, which goes as r 2 , they are able to derive the next term in the equation for J which is now For the random field r, the energy of a region of size b d is proportional to b d/2 . Any corrections requires forming a domain of 'wrong spins' which, being akin to a barrier crossing problem, is exponentially suppressed. Hence the equation for r is given by with exponentially small corrections. These two equations together can be used to derive Eq.( 39). Bray and Moore conjecture that Eq.( 41( holds exactly to all orders in w (up to exponential corrections). However, it is possible for Eq. (40) to have higher order terms in w and thus Eq.( 39) is only correct to order w 5 . Integrating Eq.( 39), we get ∼ −1/(2Aw 2 ) + 1/(2Aw 2 0 ). This implies that the correlation length is For finite size systems, the system size L ∼ exp(1/(2Aw 2 0 ). Meinke and Middleton [27] showed that their finite size data was much better fit by a function of the form w −2y 0 exp(C/w 2 0 ) where C is a constant they fit to (C = ∆ 0 in their notation) and y = 1.07. We will show that this prediction is consistent with the results of normal form theory.
As we have already argued, there is no reason Eq.( 39) is true to all orders in w. Indeed the, normal form prediction for the flow equations can be derived in a straightforward way. Consider adding a term A n w n to Eq.( 39) at = 0. This is a resonance and cannot be removed usually under normal form theory. Suppose we make a change of coordinates w =w + a nw n−2 . Then, to order O(w n ), we get dw d = Aw 3 + (3Aa n − Aa n (n − 2) + A n )w n + O(w n+1 ).
(43) We can set the coefficient ofw n = 0 if we use a n = A n /((n − 5)A). This procedure fails for n = 5 but works for all n > 5.
[74] Hence, the normal form of the equilibrium RFIM is given by As before, we have used the freedom to rescale w to set the coefficient of thew 3 term to 1. The solution of this equation gives us an expression for the correlation length This scaling form could explain the data in Meinke and Middleton with D as a fit parameter. Notice that for this to work, D must be positive. However, this solution has the strange property that the correlation length goes to 0 for w 2 = 1/D. If w 2 > 1/D, w 2 (l) decreases till it reaches 1/D. If w 2 < 1/D, it increases till it reaches 1/D. As in the 4-d Ising model, it may be more useful to consider instead the flow equation This gives the scaling form This is exactly consistent with the scaling form Meinke and Middleton use to collapse their data. Their data would predict the universal value for D = 2.14 [75]. Any system in the same universality class should see a value of D consistent with this value. However, different values of D would correspond to different universality families within the same class. We now turn to discussing the XY model before returning to the Ising model in dimensions 3, 2, and 1.

C. 2-D XY Model
The 2-d XY model is a remarkable system for several reasons. It was the site of recently celebrated insight into the connection between ground-state topology and phase transitions [76]. Thermodynamic quantities have essential singularities at its phase transition, not ordinary power laws, and their derivatives remain continuous to arbitrary order, making its phase transition infinite order [48,77,78]. This is related the fact that its RG flow equations are inherently nonlinear: they have no relevant and two marginal state variables and the procedure laid out by (11) for removing higher order terms from the flow equations contributes nothing to their simplification.
The XY model is usually posed as ferromagnetically interacting planar spins. Its partition function is exactly equivalent to the product of a trivial Gaussian modelcorresponding to spin wave degrees of freedom-with a neutral Coulomb gas-corresponding to the interaction of spin vortices [79]. The latter component contains the interesting critical behavior, which is characterized by these vortices going through an unbinding transition. The flow equations for a Coulomb gas in dimension d are given by where K ∼ T −1 and y is the fugacity of the vortices [57], which for an XY model is a function of temperature and cannot be tuned independently but is a free parameter in other equivalent models, e.g., the Coulomb gas itself. For d > 2 there is no phase transition in this system, and for d < 2 a nontrivial unstable fixed point appears and there is a phase transition in the hyperbolic universality family. It is worth noting that these flow equations do not describe the XY model for any dimension besides d = 2; 2 is the upper critical dimension of the Coulomb gas and these flow equations, while it is the lower critical dimension for the XY model. At d = 2 the flow equations undergo a novel bifurcation: there appears a line of stable fixed points at y = 0 for all K > 2, terminating at K = 2. This termination is the Berezinskii-Kosterlitz-Thouless (BKT) critical point. The flow equation near this point with x = K − 2 is These flow equations are zero to linear order and have zero Jacobian at the fixed point.
In principle arbitrary higher-order terms in these equations exist, but there are several constraints on their form. There is a symmetry y → −y in the partition function arising from the neutrality condition-y enters the partition function in factors of y − r n 2 r for r n r = 0which implies that dx/dl be even in y and dy/dl be odd. In addition, when the fugacity is zero the model is trivial and x cannot flow, meaning that dx/dl must only have terms proportional to y. Having applied these constraints, the simplest normal form has been proven by induction in polynomial order (Appendix A of [30]) to take the form For the BKT point in the sine-Gordon model, which is thought to display to the same universality as the XY model, it is known that b 0 = 3/2 [30,80]. An infinite number of coefficients remain, represented here in the form of the Taylor coefficients of an analytic function f (x 2 ). These numbers are universal in the sense that there is no redefinition ofx andỹ such that the flow equations take on the form above and contain different coefficient values. Unlike those in the previous sections, this bifurcation does not have a named classification as far as the authors know. A constant of the RG flow can be found by integrating these forms. First, dividing the equations (54) by (53) (and dropping the tildes), we find which separates into Integrating both sides and choosing l 0 such that x(l 0 ) = 0, we find dx. (58) It follows that is a constant of the flow. The expansion of the integral can be taken to arbitrary order with ordinary computer algebra software. The finite-size behavior of the flow is rather complicated and doesn't yield closed form results, details can be found in [30]. The XY model and other infinite-order transitions are usually characterized by the anomalous exponent σ parametrizing the essential singularity in the correlation length, which for the BKT transition is σ = 1/2 [48]. Conformal field theory predicts the presence of infinitely many models with this anomalous exponent [81]. The value of σ been shown to be fixed by the quadratic-order truncation of the system's flow equation, independent of any higherorder terms [82]. There are six possible quadratic-order terms in flow equations with two variables. Of these, two can be removed by linear transformations of the two variables. Two more can be set to 1 by rescaling the variables. Hence, there are two parameters at quadratic order which determine the universality family that the system belongs to, and infinite number of subsequent terms which determine the universality class. Giving a full classification of possibilities is beyond the scope of this paper but we give some examples below. For instance, when the requirement of symmetry under y → −y is lifted, the flow equations can no longer be brought to the form Eqs. (53) and (54); though the simplest form that results isn't yet known, it is certainly different from the symmetric case, a fact that can be found by simply trying to eliminate the nonsymmetric cubic terms. In such a case the codimension of the bifurcation would likely be different, corresponding to the fact that no a priori reason exists for the vanishing of the term linear in y in dx/dl. Such linear terms would change the universality family. Among the infinite collection of BKT-like conformal theories-including the many physical models identified as having a BKT-like transition because their behavior resembles Eq. (62), like percolation in grown networks [58,59]. These may already be examples of models with σ = 1/2 but belonging to another universality class. It could also be the case that all BKT-like transitions are in fact members of the same universality family and class.
Other universality classes and families definitely do exist, characterized by novel values for σ. The level-1 SU(N ) Wess-Zumino-Witten model has been found to be characterized by σ = N/(N + 2) [83]. Dislocatedmediated melting alone has produced a melange of anomalous exponents, with σ = 1/2, σ = 2/5, and σ = 0.369 63 . . . depending on precise specification of the model and the lattice geometry [84,85]. Topological transitions in systems whose vortices are non-Abelian produce several series of σ values dependent on particular symmetry [86]. Each value of σ indicates either a different universality family or merely a different class within the same family depending on how it affects the terms at quadratic order. A classification of possible bifurcations and corresponding simplest normal forms is in order for flow equations whose leading order is quadratic, and whose expansions are constrained or not by various symmetries. This would be the first step in developing techniques for distinguishing between universality classes and families of this type using experimental or simulation data.

D. 3-d Ising model
There is a sense in which the Ising model is simplest in 3 dimensions because it is part of the hyperbolic universality family. It is also the first natural application of the expansion. The transcritical bifurcation at 4 dimensions leads to an exchange of stabilities of the Gaussian fixed point and the Wilson-Fisher fixed point at a nonzero value of u = u * . About this Wilson-Fisher fixed point, the flow equations of the 3-d Ising model are in the hyperbolic universality class with linear coefficients which define the Ising universality class.
However, another approach is to consider the scaling form as a function of the dimension in a way that is well defined even at = 0. Doing this, naturally requires us to keep nonlinear terms in the equation because we already know that the 4-d Ising model has nonlinear terms in its flow equations.
We want to write the flow equations about the 3-d fixed point but keep the nonlinear terms required for the scaling form to have the correct limiting behavior in 2-d and in 4-d. We can write the normal form of the flow equations as We have included the nonlinear terms in u required for the correct scaling behavior and the resonance between the temperature and the free energy. As usual, we switch notation to t, h and u with the understanding that they are different from the normal form variables by analytic corrections. Let us look at the scaling variable formed with t and u which can be obtained by solving The solution of this equation gives the scaling variable where u 1 and u 2 are the two non-zero roots of the denominator on the r.h.s of Eq.67 which to first order in λ u are given by u 1 = λ u and u 2 = 1/D − λ u . The form of the scaling variable is interesting, it is essentially given by a product of the linearized scaling variables at the three fixed points that the equation has. Taking the limit → 0, we get which is the right scaling variable in 4-d. We have not yet been able to obtain an analytical form for the scaling variable involving t and h. This is because the equation for u(l) does not seem to have a closed-form solution here (unlike the 4-d case). Nevertheless, we are motivated by an attempt to create scaling variables which interpolate between different dimensions and have the correct scaling behavior in many dimensions going down from 4 to 1.
Once the full scaling variables are written down, a first test would be to see if these scaling variables do better collapsing the numerical data in 3-d.

E. 1-d Ising model
The 1-d Ising model is somewhat different because it is the lower critical dimension and does not have a phase transition. The 1-d Ising model has an exact solution which can be obtained by using transfer matrices. The partition function can be written as the trace of a transfer matrix T N where N is the number of spins in the system. The matrix T ij = e −βH(sisj ) . Coarse graining here can be done by a well-defined procedure, the coarse grained transfer matrix is defined asT = T b where b is the coarse graining length scale. Defining = log b and expanding for b close to 1, we can get flow equations for the temperature T This is different from the flow equations we have considered so far because of the presence of non-analytic terms in the flow. The non-analytic term which multiplies the T 2 term is = −1 at T = 0. So, this equation corresponds to a transcritical bifurcation where the additional terms are non-analytic at T = 0. This can be used to derive a correlation length χ ∼ exp(2/T ). To interpret the flow further, consider the change of coordinates κ = exp(−2/T ). In these variables, the flow is Evidently, the flow is analytic in this variable. Solving the full flow Eq.( 70) gives χ ∼ −1/(log tanh(1/T )). For non-zero , this argument is usually extended in what is called a Migdal-Kadanoff procedure for doing RG [87,88]. The flow equations are identical except for the presence of a − T term which serves as the bifurcation parameter. The 1+ expansion can be summed completely because the flow equation is known to all orders. It does not yield very accurate critical exponents though it gives the exact value of the critical temperature in 2-d (because it respects duality symmetry). Several people have improved the expansion [89,90].
The presence of non-analytic terms in the flow equations complicates the application of normal form theory.
We will come back to it when discussing Legendre transform of flow equations.

F. 2-d Ising model
The 2-d Ising model is a particularly nice example because it has an exact solution in the absence of a magnetic field. All predictions then can be compared to the exact solution. Surprisingly, despite the known exact solution, the scaling behavior of the 2-d Ising model is still not completely understood. A full discussion of the 2-d Ising model will be given in separate work [53]. Here, we give a brief summary of the issues involved.
The only variable required to describe the 2-d Ising model in the absence of a field is the temperature t. The linear eigenvalues of the free energy and the temperature are 2 and 1 respectively. The normal form of the flow equations can be written as We have used the fact that the only term which cannot be removed by traditional normal form analysis is the resonance t 2 . In fact, it cannot be removed by any analytic change of variables. We have also used the freedom to rescale t to set the coefficient of the resonance equal to -1 [91]. The solution to this can be written ast =t 0 e and the free energỹ f (t 0 , ) = e −2 f (t 0 e ) −t 2 0 .
Now, the normal form variablet 0 is some analytic function of the physical variable t 0 . It is linear to first order in t 0 . Hence, we can write it ast 0 = t 0 (1 + c(t 0 )) where c is some analytic function. Then, we can expand f (t 0 ) = t 2 0 (1 + c(t 0 )) 2 f (1) + t 2 0 (1 + c(t 0 )) 2 log t 0 (1 + c(t 0 )), = t 2 0 (1 + c(t 0 )) 2 f (1) + t 2 0 (1 + c(t 0 )) 2 log(1 + c(t 0 )) + t 2 0 (1 + c(t 0 )) 2 log t 0 , where both a(t 0 ) and b(t 0 ) are some analytic functions of t 0 . Meanwhile any change of coordinates which adds an analytic function of t 0 tof can be absorbed in the definition of a(t 0 ). Hence, we can write the final most general form of the free energy of the 2-d Ising model as f = a(t 0 ) + b(t 0 ) log t 0 . Indeed, the exact solution of the 2-d Ising model can be written in this form [92].
While the basic solution of the 2-d Ising model is sim-ple, some challenges still remain. The scaling form in the presence of other variables (like the magnetic field and other irrelevant variables) which has so far only been conjectured [52,92] naturally follows from an application of normal form theory. It is given simply by including other variables in the argument of the free energy in Eq. ( 75) before coarse graining till t( ) = 1. Irrelevant variables are the source of singular corrections to scaling. An interesting unresolved issue is the presence of higher powers of logarithms in the susceptibility which are not found in the free energy [93,94]. This is usually attributed to the presence of irrelevant variables. Here it is possible to show that the irrelevant variables which are derived from conformal field theory [92] would in fact lead to higher powers of logarithms in the free energy which are not observed. Hence, they cannot explain the higher powers of logarithms in the susceptibility. It is possible that there are other irrelevant variables in the 2-d square lattice nearest neighbor Ising model with a field which are not predicted by conformal field theory but can capture the higher powers of logarithms in the susceptibility, as they turn on with a field. The logarithm due to the resonance in the 2-d Ising model is most apparent in the specific heat. It is easy to derive the flow equations for the inverse specific heat which have the form and has a transcritical bifurcation in two dimensions. This raises a question, is it legitimate to talk about a bifurcation in two dimensions for the Ising model if it happens in the space of results rather than the space of control variables? Intriguingly, though perhaps unrelated, a bifurcation has been observed in 2 dimensions using methods of conformal bootstrap [95,96]. In thermodynamics, a natural framework which interchanges between results and control parameters is given by Legendre transforms. However, the flow equations for the Legendre transformed coordinates generically have nonanalyticities in them. We suspect that the variable t (and h etc.) is uniquely specified as the correct variable for RG. It is possible that it is more natural to consider removing degrees of freedom in the canonical ensemble (t and f ), then in a microcanonical one (E and S) [97]. A fuller discussion will be given in forthcoming work [53].

VI. CONCLUSION
We have shown how normal form theory leads to systematic procedure for handling the singularity in RG flows. The concept of universality families broadens the notion of a universality class and we have elucidated it with several different examples. We have focused on getting a precise handle on the singularity at the critical point. However, normal form theory also gives an elegant way to fit corrections to scaling. Interestingly, even the scaling of the 2-d Ising model which has an exact solution has some unresolved mysteries which we are exploring. It is possible that interpolating between dimensions in a way that captures the correct singularities can improve scaling collapses in all dimensions. Finally, we are exploring the application of our methods to systems like jamming in 2-d [37], where logarithmic corrections are observed but no renormalization-group theory is available. In general, we expect this fruitful confluence of dynamical systems theory and the renormalization group will not only clarify and illuminate previously known technical calculations, but will also facilitate quantitative analysis of experimental and theoretical systems farther from their critical points and before the underlying field theory is well understood.