Léonie Canet

Hugues Chaté

Bertrand Delamotte

Nicolás Wschebor

nonperturbative renormalization group for the Kardar-Parisi-Zhang equation: general framework and first applications

Keywords: numbers: 64.60.Ht, 05.10.Cc, 68.35.Fx, 05.70.Ln, 05.40.-a

published or not. The documents may come

I. INTRODUCTION

In their seminal work [2], Kardar, Parisi and Zhang proposed a stochastic continuum equation to describe surface growth through ballistic deposition which reads ∂h(t, x) ∂t = ν ∇ 2 h(t, x) + λ 2 ∇h(t, x) 2 + η(t, x) (1) where η is an uncorrelated white noise of strength D, η(t, x)η(t ′ , x ′ ) = 2 D δ d ( x -x ′ ) δ(t -t ′ ), which models randomness in deposition. The term ν ∇ 2 h(t, x) provides the interface with a smoothening mechanism. The insightful feature of Eq. ( 1) is to take into account the nonlinearity of the growth velocity through the inclusion of the term λ(∇h(t, x)) 2 which plays an essential role in the large scale properties of the height profile h(t, x).

The KPZ equation ( 1) is maybe the simplest nonlinear Langevin equation showing non-trivial behavior [3], and as a consequence it arises in connection with an extremely large class of nonequilibrium or disordered systems [3] such as randomly stirred fluid (Burgers equation) [4], directed polymers in random media [5], dissipative transport [6,7] or magnetic flux lines in superconductors [8]. The KPZ equation has thus emerged as one of the fundamental theoretical models to investigate universality classes in nonequilibrium scaling phenomena and phase transitions [3]. It is only recently, though, that a definitively-convincing experimental realization has been brought out for a one-dimensional interface, confirming detailed theoretical predictions [9].

The profile of the stationary interface is usually char-acterized by the two-point correlation function

C(t, | x|) ≡ [h(t, x) -h(0, 0) -t ∂ t h ] 2 (2)
and, in particular, by its large scale properties. The KPZ growth leads to generic scaling. At long time τ and large distance L, C assumes the scaling form C(τ, L) ∝ τ 2χ/z g(L 2χ τ -2χ/z ) without fine-tuning any parameter of the model. The scaling function g is universal and has the asymptotics g(y) → const. as y → 0 and g(y) ∼ |y| as y → ∞. χ and z are the universal roughness and dynamical exponents respectively. In fact, these two exponents are not independent since the Galilean symmetry [4] -the invariance of Eq. ( 1) under an infinitesimal tilt of the interface -enforces the scaling relation z + χ = 2 as long as λ = 0. The KPZ equation encompasses two distinct scenarios depending on the dimension of the interface. Above two dimensions, there exist two different regimes separated by a critical value λ c of the nonlinear coefficient [2,4]. In the weak-coupling regime (λ < λ c ), the interface remains smooth, its properties are determined by the λ = 0 (Gaussian) fixed point -corresponding to the linear Edwards-Wilkinson equation [3] -with exponents χ = (2 -d)/d and z = 2. In the strong-coupling regime (λ > λ c ), the nonlinearity becomes relevant and the interface roughens. In this regime, the exponents are not known exactly, and some important issues are still controversial, such as the existence of an upper critical dimension (for a recent discussion, see e.g. [10]). This unsatisfactory situation has persisted up to very recently [11] because the strong coupling phase of the KPZ equation has remained out of reach of controlled analytical approaches. In particular, standard perturbation expansions [2,12,13] have been proved to fail at all order to find a strong coupling fixed point [14]. Some nonperturbative approaches have been devised, such as the modecoupling (MC) approximation [6,[15][16][17][18], or the weak noise scheme [19], but they are difficult to improve in practice. (Some specific comments regarding the MC approximation are put forward throughout this paper). Let us finally mention the self-consistent expansion, which is performed on the Fokker-Planck equation and can be systematically improved [20][21][22]. It has yielded important results for the KPZ equation in arbitrary dimensions, where it predicts an upper critical dimension [10]. Note that it has also been applied to different extensions of the KPZ model, such as the non local KPZ equation [20] or the KPZ equation with colored noise [21].

In two dimensions and below, the situation is different, both from a physical and a theoretical point of view. Physically, the interface always roughens, there is no transition. On the theoretical side, exact results are available in one dimension [1,[23][24][START_REF]During the completion of this work, the complete exact solution of the one-dimensional equation for curved initial conditions has[END_REF][26]. The critical exponents χ = 1/2 and z = 3/2 have been known for long since they are fixed by the existence of an incidental fluctuation-dissipation theorem in this dimension [2], but the scaling function g(y) (and other universal properties) have been computed exactly only very recently [1,23,[START_REF]During the completion of this work, the complete exact solution of the one-dimensional equation for curved initial conditions has[END_REF]26]. Note that there had been earlier attempts to determine the scaling function, in the framework of the mode-coupling approximation [15,27,28]. In Ref. [28], a refined ansatz for the scaling function in one dimension is devised to solve self-consistently the MC equations and it turns out that the result compares quite accurately with the exact solution. The one-dimensional scaling function has also been obtained within the selfconsistent expansion in one-dimension [22] and coincides in many respects with the MC result (see Section VI). Regarding the MC theory, it is an 'adhoc' approximation which consists of resumming one-loop diagrams while discarding vertex corrections. The quality of the results is all the more surprising that it was shown that the contribution of the neglected terms are of the same order as those kept [15]. The main drawbacks of the MC approach are that it strongly relies on the quality of the 'educated guess' for the ansatz -only available in one dimension up to now -and it cannot in practice be systematically improved by calculating higher orders.

Recently, we have proposed an analytical approach to the KPZ equation based on nonperturbative renormalization group (NPRG) techniques [11]. This early work has shown that the NPRG flow equations embedded a strong-coupling fixed point in all dimensions, and it has thus provided a complete, qualitatively correct phase diagram for the first time within a RG approach, as well as reasonable values for the strong-coupling fixed point exponents in physical dimensions [11].

In this paper we present a general and systematic framework for applying NPRG methods to the KPZ problem, strongly constrained by the symmetries of this model. We derive general Ward identities and introduce a 'covariantization' associated with the Galilean symmetry which, to the best of our knowledge, has never been reported before in the literature. Within this formalism, successive orders of approximation are easily made explicit. In the present contribution, we implement the minimal order in the response field of this approximation scheme. We review the (simplified) first account of it presented in [11] -postponing its complete revisited version -and derive new results related to the one-dimensional problem, in order to confront our approach to the available exact results. In particular, we compute the correlation function, show that it takes a scaling form at long time and large distance and extract the associated scaling function. The NPRG results are in even better agreement with the exact results of Prähofer and Spohn [1] than those obtained under the MC or the self-consistent approximations.

We stress the advantages of the NPRG formalism: (i) it is based on an exact flow equation, which, given a microscopic model, yields the macroscopic properties of the system. The only input is the bare action, that is, no a priori knowledge, other than the microscopic model and its symmetries, is required to compute the physical observables. In particular, one does not have to assume scaling at long time and large distance nor the form of the scaling function; (ii) the approximations, which are mandatory to solve the exact flow equation can be systematically improved and implemented in any dimension and, in contrast to most other nonperturbative approaches, the calculation of higher orders is achievable and have been performed in practice in other systems [29,30]; (iii) beyond scaling, a wealth of quantities can be calculated (for example, we compute here the correction to scaling exponent ω); (iv) it has already yielded many nontrivial and accurate results in systems at equilibrium such as frustrated magnets [31], the random bond and random field Ising model [32], membranes [33], bosonic systems [34], but also in systems out-of-equilibrium, where one can mention important advances in reaction-diffusion systems [35].

As the NPRG is a field theoretical method, our starting point is the field theory associated with Eq. ( 1), which follows from the standard procedure of Janssen-de Dominicis relying on the introduction of a response field h and sources (j, j) [START_REF] Janssen | C. de Dominicis[END_REF]. This procedure allows one to explicitly carry out the integration over the Gaussiandistributed noise η in Eq. ( 1) upon doubling the number of fields (see e.g. [38]). The generating functional reads:

Z[j, j]= D[h, i h] exp -S[h, h] + x jh + jh (3) S[h, h]= x h ∂ t h -ν ∇ 2 h - λ 2 (∇h) 2 -D h2 (4)
where x = (t, x).

The remainder of the paper is organized as follows. In Section II, we briefly review the nonperturbative renor-malization group formalism for out-of-equilibrium problems. In Section III, we analyze in detail the symmetries of the KPZ action (4) and derive Ward identities associated with the linearly realized ones. In Section IV, we build an approximation scheme based on a covariantization procedure rooted in the symmetries, and derive explicitly an ansatz at the minimal order in the response field. The determination of the complete phase diagram and of critical exponents in all dimensions, using a simplified version of this ansatz, is reported in Section V. Section VI is then devoted to the calculation (using the full ansatz) of the scaling function and of some other universal quantities in one dimension, which are compared with their exact counterparts. Technical details, such as some of the Ward identities, a discussion of the validity of our approximation scheme, the computation of vertex functions or the procedures for the numerical integration of the flow equations are reported in Appendices A, B, C and D respectively.

II. THE NONPERTURBATIVE RENORMALIZATION GROUP

The NPRG formalism relies on Wilson's RG idea, which consists of building a sequence of scale-dependent effective models such that fluctuations are smoothly averaged as the (momentum) scale κ is lowered from the microscopic scale Λ, where no fluctuations are yet included, to the macroscopic one κ = 0 where all fluctuations are summed over [29,37].

For out-of-equilibrium problems, one formally proceeds as in equilibrium, but with the presence of additional response fields and with special care required to deal with the consequences of Itō's discretization and with causality issues, as stressed in detail in [38] -from which conventions are taken throughout this paper. For future use, we define the Fourier conventions used in [38] and throughout this work:

f (ω, p) = d d x dt f (t, x) e -i p• x+iωt (5) f (t, x) = d d p (2π) d dω 2π f (ω, p) e i p. x-iωt (6)
≡ p f (p) e i p. x-iωt , (7)
where p = (ω, p).

To achieve the separation of fluctuation modes within the NPRG procedure, one adds to the original action S, a momentum-dependent mass-like term:

∆S κ = 1 2 q h i (-q) [R κ (q)] ij h j (q) ( 8
)
where the indices i, j ∈ {1, 2} label the field and response field respectively h 1 = h, h 2 = h, and summation over repeated indices is implicit. The matrix elements [R κ (q)] ij are proportional to a cutoff function r(q 2 /κ 2 ) (see Section III D), with q = q , which ensures the selection of fluctuation modes: r(x) is required to vanish as x 1 such that the fluctuation modes h i (q κ) are unaffected by ∆S κ , and to be large when x 1 such that the other modes (h i (q κ)) are essentially frozen. Since ∆S κ must preserve the symmetries of the model, we postpone the discussion of the precise structure of the matrix elements [R κ (q)] ij to Section III D after the analysis of these symmetries.

In presence of the mass term ∆S κ , the generating functional (3) becomes scale dependent (9) and the effective action Γ κ [ϕ, φ], where ϕ i = h i j, j are the expectation values of the fields h i in the presence of the external sources j and j, is given by the Legendre transform of W κ = log Z κ (up to a term proportional to R κ ) [29,38]:

Z κ [j, j]= D[h, i h] exp -S -∆S κ + x jh + jh
Γ κ [ϕ, φ] + log Z κ [j, j] = j i ϕ i - 1 2 q ϕ i [R κ ] ij ϕ j . ( 10
)
From Γ κ , one can derive two-point correlation and response functions,

[ Γ (2) κ ] i1i2 (x 1 , x 2 , ϕ, φ) = δ 2 Γ κ [ϕ, φ] δϕ i1 (x 1 )δϕ i2 (x 2 ) (11)
and more generally n-point correlation functions that we write here in a 2 × 2 matrix form (omitting the dependence on the running scale κ)

Γ (n) i3,...,in (x 1 , ..., x n , ϕ, φ) = δ n-2 Γ (2) (x 1 , x 2 , ϕ, φ) δϕ i3 (x 3 )...δϕ in (x n ) . (12)
The exact flow for Γ κ [ϕ, φ] is given by Wetterich's equation which reads (in Fourier space) [29]:

∂ κ Γ κ = 1 2 Tr q ∂ κ R κ • G κ with G κ = Γ (2) κ + R κ -1 (13)
the full renormalized propagator of the theory. When κ flows from Λ to 0, Γ κ interpolates between the microscopic model Γ κ=Λ = S and the full effective action Γ κ=0 that encompasses all the macroscopic properties of the system [38]. Differentiating Eq. ( 13) twice with respect to the fields and evaluating it in a uniform and stationary field configuration (since the model is analyzed in its long time and large distance regime where it is translationally invariant in space and time) one obtains the flow equation for the two-point functions:

∂ κ [ Γ (2) ] ij (p) = Tr q ∂ κ R(q) • G(q) • - 1 2 Γ (4) ij (p, -p, q) +Γ (3) i (p, q) • G(p + q) • Γ (3) j (-p, p + q) • G(q) (14)
where the κ and background field dependencies have been omitted, as well as the last argument of the Γ (n) which is determined by frequency and momentum conservation [38].

Solving Eq. ( 13) (or Eq. ( 14)) is in principle equivalent to solving the model. In practice this resolution cannot be performed exactly since (13) is a nonlinear integral partial differential functional equation. Hence one has to devise an approximation scheme. The main constraint on this approximation scheme is to preserve the symmetries of the problem. We thus now revisit the symmetries of the KPZ action.

III. SYMMETRIES OF THE KPZ ACTION

The KPZ action (4) possesses well-known symmetries, in addition to translation and rotation invariances: i) the Galilean symmetry and ii) the h-shift symmetry, which can be expressed as the invariance of the action (4) under the following transformations:

i) h ′ (t, x) = x • v + h(t, x + λ vt) h′ (t, x) = h(t, x + λ vt). ( 15
) ii) h ′ (t, x) = h(t, x) + c (16)
where v and c are arbitrary constant quantities.

In one dimension, the KPZ equation also satisfies a fluctuation-dissipation theorem that fixes the exponents exactly. This property roots in a time-reversal symmetry of the action which, as shown in [39], can be encoded in the transformation iii)

h ′ (t, x) = -h(-t, x) h′ (t, x) = h(-t, x) + ν D ∇ 2 h(-t, x). (17)
One can check that the invariance of the action (4) under the transformation (17) requires that the contribution ∇ 2 h(∇h) 2 vanishes, which is 'incidentally' true only in d = 1. The time-reversal symmetry thus only holds in this dimension. This set of symmetries entails Ward identities for the n-point vertex functions, as long as the mass term ∆S κ is appropriately chosen (see Section III D).

In fact, there exist even stronger symmetries of the KPZ action that, to our knowledge, were only pointed out in [40]. They consist of 'gauging' in time the transformations i) and ii) in the following way:

i') h ′ (t, x) = x • ∂ t v(t) + h(t, x + λ v(t)) h′ (t, x) = h(t, x + λ v(t)) (18)
ii')

h ′ (t, x) = h(t, x) + c(t) (19)
which we will refer to as Galilean-gauged and shiftgauged symmetries respectively. In these gauged versions, v(t) and c(t) are arbitrary infinitesimal functions of time. Note that the action (4) is not strictly invariant under the transformations (18) and ( 19) but the corresponding variations of the action are linear in the fields, and this behavior also yields useful Ward identities, which we derive in the following. These stronger forms of the symmetries will be thoroughly exploited.

Finally, an additional Z 2 symmetry, which is manifest on the Cole-Hopf version of the theory, is nonlinearly realized in the KPZ action (4). The Cole-Hopf field transformation writes as follows:

h(x) = 2ν λ log |w(x)| h(x) = w(x) w(x). ( 20
)
In terms of w and w -upon rescaling these fields and time -the KPZ action becomes:

S[w, w] = x w(∂ t w -∇ 2 w) - 1 4 g b (w w) 2 (21)
with the bare coupling constant

g b = λ 2 D ν 3 . ( 22
)
This action is invariant under the simple Z 2 transformation (iv) w(t, x) → w(-t, x), w(t, x) → w(-t, x). However, in terms of the original fields h and h, this transformation becomes iv)

h ′ (t, x) = -h(-t, x) + 2ν λ log | h(-t, x)| h′ (t, x) = h(-t, x). (23)
Its highly nonlinear form renders complicated the induced Ward identities among vertex functions. Consequently, they are not given here as they will not be exploited in the following.

A. Shift-gauged symmetry

The Ward identity associated with the shift-gauged symmetry can be derived by performing in the functional integral (9) the change of variables corresponding to the transformation (19). As this operation must leave the value of the integral unchanged, one obtains that

x c(t)j(x) -h(x) j, j ∂ t c(t) -∆S κ [c(t), h(x)] j, j (24)
must vanish. Here, [R κ (q)] 11 has been set to zero (as in Eq. ( 41)) for simplicity and for causality issues (see Section III D). As the mass term is quadratic in fields, ∆S κ [c(t), h(x)] is linear in h, and thus

∆S κ [c(t), h(x)] j, j = ∆S κ [c(t), h(x) = φ(x)].
The expression (10) of the Legendre transform then implies

x c(t) δΓ κ [ϕ, φ] δϕ(x) -φ(x)∂ t c(t) = 0. ( 25
)
After integrating by parts the second term, we conclude that the functional

Γ κ [ϕ, φ] - x φ(x)∂ t ϕ(x) (26)
is invariant under the transformation (19). In other words, the only non-invariant term φ∂ t ϕ of the bare action is not renormalized and the rest of the action is shift-gauged symmetric.

Let us express this property on the n-point vertex functions, that we denote from now on as Γ (l,m) κ (x 1 , . . . , x l+m ) which stands for the Γ (n=l+m) κ vertex involving l (respectively m) legs -derivatives of Γ κ with respect to ϕ (respectively φ) -with the l first frequencies and momenta referring to the ϕ fields and the m last to the φ fields. At this stage, it is convenient to work in Fourier space. The n-point vertex function are defined in Fourier space by

(2π) d+1 δ d+1 i p i Γ (l,m) κ (p 1 , . . . , p l+m-1 ) = x1•••x l+m Γ (l,m) κ (x 1 , . . . , x l+m )e i i ( xi• pi-tiωi) ( 27
)
where again, the last frequency and momentum, fixed by translational invariance in time and space, are implicit. The shift-gauged symmetry entails that the n-point vertex functions in Fourier space satisfy the following property:

Γ (m,n) κ (ω 1 , p 1 = 0, . . . , p m+n-1 ) = iω 1 δ m1 δ n1 (28)
which means that, apart from the contribution of φ∂ t ϕ to Γ

(1,1) κ

, the vertices vanish upon setting the momentum of one of the ϕ to zero. This is related to the fact that the field ϕ only appears in Γ κ with gradients. In particular, it roots the non-renormalization of the kinetic term φ∂ t ϕ of the KPZ action, which is well established in perturbation theory (see e.g. [15]).

B. Galilean symmetry

Global Galilean symmetry

Let us first review the Ward identities associated with the standard Galilean symmetry. As for the shift-gauged symmetry, one can prove that, as the transformation (15) encoding the Galilean symmetry of the bare ac-tion is affine, and provided the mass term ∆S κ is chosen Galilean-invariant (as in the form (41) -see Section III D), Γ κ also possesses this symmetry. The corresponding Ward identity is:

x ( x • v + λt v • ∇ϕ(x)) δΓ κ δϕ(x) +λt v • ∇ φ(x) δΓ κ δ φ(x) = 0, (29)
that is Γ κ is invariant under the same Galilean transformation (with the same parameter) as the bare action.

One can then proceed to the derivation of the Ward identities for the vertex functions, by taking functional derivatives of Eq. ( 29) with respect to the fields and then evaluating them for uniform and static fields. For instance, one gets the standard identity for the three-point function:

i ∂ ∂ p Γ (2,1) κ (ω = 0, p = 0; ω 1 , p 1 ) = λ p 1 ∂ ∂ω 1 Γ (1,1) κ (ω 1 , p 1 ). (30)
The Ward identity for a generic n-point function can be derived and one obtains:

i ∂ ∂ p Γ (m+1,n) κ (ω = 0, p = 0; p 1 ; . . . ; p m+n-1 ) = λ p 1 ∂ ∂ω 1 + • • • + p m+n-1 ∂ ∂ω m+n-1 Γ (m,n) κ (p 1 ; . . . ; p m+n-1 ). ( 31
)

Galilean-gauged symmetry

Let us now come to the gauged form (18) of the Galilean symmetry. As for the shift-gauged symmetry, the variation of the action under this transformation is linear in the fields and consequently, it entails a Ward identity which reads

d d x λ∇ϕ(x) δΓ κ δϕ(x) -x∂ t δΓ κ δϕ(x) +λ∇ φ(x) δΓ κ δ φ(x) -x∂ 2 t φ(x) = 0. ( 32
)
From this functional identity, one can again deduce identities for the vertex functions. They bare similar expressions as those for the global Galilean symmetry, but with a stronger content. For instance, the identity for the three-point function becomes:

iω ∂ ∂ p Γ (2,1) κ (ω, p = 0; ω 1 , p 1 ) = λ p 1 × Γ (1,1) κ (ω + ω 1 , p 1 ) -Γ (1,1) κ (ω 1 , p 1 ) , (33)
which in the limit ω → 0 coincides with (30). The gauged identity ( 33) is stronger as it constrains the whole frequency dependence and not only the zero-frequency sector. One could derive similar identities for generic npoint functions, but we now stress, instead, a more efficient way to exploit the Galilean symmetry, that will guide our construction of the approximation scheme.

Covariant time derivatives

The previous Ward identities ensuing from the Galilean symmetry do not clearly reflect the geometrical interpretation of this symmetry. In order to do so, one can analyze the Galilean invariance from another angle. One can build quantities which are manifestly scalar under the Galilean transformation (15), upon introducing an adequate covariant time derivative. Let us define a function f (x) as a scalar under the Galilean symmetry if its infinitesimal transform under ( 15) is given by

δf (x) = tλ v • ∇f (x). ( 34
)
With this definition, if f is a scalar then d d xf is invariant under the Galilean transformation and can be used to build an action possessing the Galilean symmetry. It follows from this definition that the response field h is a scalar, but that the field h is not unless one takes two successive space derivatives ∇ i ∇ j h. The gradient of a scalar remains a scalar, but not its time derivative since

δ(∂ t f (x)) = λ v • t∇(∂ t f (x)) + ∇f (x) . (35)
However, as in fluid mechanics, one can construct a covariant time derivative

Dt ≡ ∂ t -λ∇h(x) • ∇
which conserves the scalar property, i.e. if f is a scalar then so is Dt f . Note that the covariant time derivative of the field h itself, that we denote D t , bares a special form with a 1/2 factor:

D t h(x) ≡ ∂ t h(x) - λ 2 (∇h(x)) 2 (36)
since h is not a scalar on its own. These covariant derivatives will constitute the building blocks in the construction of our approximation scheme (see Section IV).

C. Time-reversal symmetry in d = 1

In d = 1, the action (4) with a mass term of the form (41) exhibits the additional time-reversal symmetry. Note that this is a discrete symmetry (i.e. there is no infinitesimal transformation corresponding to ( 17)). However, as it is linear in the fields, one can show using the same procedure as previously that Γ κ also possesses this symmetry. That is, it verifies:

Γ κ [ϕ(x), φ(x)] = Γ κ [-ϕ(-t, x), φ(-t, x)+ ν D ∇ 2 ϕ(-t, x)].
(37) Again, one can derive Ward identities for the n-point vertex functions by taking derivatives of (37) with respect to the fields and evaluating them at uniform and static field configurations. For the two-point functions, this yields, in Fourier space:

Γ (2,0) κ (ω, p) = Γ (2,0) κ (-ω, p) + ν D p 2 Γ (1,1) κ (-ω, p) + ν D p 2 Γ (1,1) κ (ω, -p) + ν D 2 p 4 Γ (0,2) κ (-ω, p), (38)
which, given that Γ

(2,0) κ (ω, p) = Γ
(2,0) κ (-ω, p) and spatial parity reduces to:

2ReΓ (1,1) κ (p) = - ν D p 2 Γ (0,2) κ (p). ( 39
)
For the three-point functions, one gets:

2ReΓ (1,2) κ (p 1 ; p 2 ) = - ν D p 2 1 Γ (0,3) κ (p 1 ; p 2 ). 2ImΓ (2,1) κ (p 1 ; p 2 ) = - ν D p 2 2 ImΓ (1,2) κ (p 1 ; p 2 ) - ν D p 2 1 ImΓ (1,2) κ (p 2 ; p 1 ). ( 40
)
Similar expressions can be derived for Γ

(3,0) κ and for the four-point functions, which are reported for completeness in Appendix A.

D. cutoff function

The mass term ∆S κ defined by ( 8) must be chosen such that the Ward identities are preserved all along the flow. This is not obvious a priori. In particular, a quadratic mass term cannot be invariant under the shift-gauged symmetry. However, as explained above, if the variation of ∆S κ under the transformation (19) is linear in the field, the Ward identity remains identical even in the presence of such a term. Similar comments hold for the Galilean symmetry, which is preserved as long as the cutoff matrix R κ does not depend on the frequency. We thus choose a frequency independent matrix R κ -which also preserves causality properties [38]. Moreover, the timereversal symmetry imposes in one dimension the following relation between the components of the cutoff matrix:

[R κ ( q)] 12 = -ν 2D q 2 [R κ ( q)] 22
, which we will impose in all dimensions.

Finally, we set [R κ ( q)] 11 = 0 since it is not necessary to consider a cutoff term proportional to ϕϕ and since the presence of such a term would spoil the property that Γ κ is proportional to φ, important for causality issues (see [38]) and also in order to preserve the shift symmetry. To summarize, the cutoff matrix is chosen with the following form:

R κ ( q) = r q 2 κ 2 0 ν κ q 2 ν κ q 2 -2D κ , (41)
where the running coefficients ν κ and D κ are introduced for convenience [11]. Here we use r(x) = α/(exp(x) -1) where α is a free parameter.

IV. APPROXIMATION SCHEME Two routes have been followed in the past to build approximation schemes to deal with exact NPRG flow equations such as Eqs. ( 13) and ( 14). Either one focuses on the long time and large distance properties, trying to describe as accurately as possible the zero-momentum and zero-frequency domain while approximating the other sectors. Or the aim is to compute the momentum and frequency dependence of two-point functions and then the approximation concerns the three-and four-point functions.

The usual implementation of the first route consists of performing a 'derivative expansion': Γ κ is expanded in powers of gradients and time derivatives of the fields [29,38]. This strategy has been undertaken to investigate the critical properties of many equilibrium [29], and also nonequilibrium systems such as reaction-diffusion processes [35], where it has yielded satisfactory results. However, for the KPZ problem, the signature nonlinear term involves a gradient, and this seems to preclude the use of the derivative expansion [39].

The second route is the one followed here. At equilibrium, one can compute the momentum dependence of the two-point function using the Blaizot-Méndez-Wschebor (BMW) approximation scheme which has proved to yield very accurate results for O(N ) models [41]. In the BMW framework, the momentum dependence of Γ κ is truncated. However, a direct implementation of the BMW scheme for the KPZ problem is hindered by the symmetries which impose strict identities between Γ (m,n) κ 's with different m and at different momenta, see e.g. Eq. ( 31), and render complicated the writing of a truncation.

To overcome this difficulty, we propose, instead of designing directly truncations on the three-and four-point functions, to construct an ansatz for Γ κ that manifestly preserves the shift-gauged and Galilean-gauged symmetries while allowing for arbitrary frequency and momentum dependencies in the two-point functions. This is easily achieved by combining together the basic Galilean scalars of the theory -φ, ∇ i ∇ j ϕ and D t ϕ, and arbitrary powers of their gradients and covariant time derivatives Dt -since, by construction, all calculations involving a functional of these quantities automatically satisfy the Galilean symmetry (the other symmetries, but the nonlinearly realized one, being easily enforced a posteriori).

On this basis, various truncations may be performed to obtain an ansatz for Γ κ . One could, for example, keep complete dependencies in φ and ∇ 2 ϕ, while treating only the zero-frequency sector. This 'super derivative expansion' option is left for future studies. An alternative choice, pursued here with the ultimate aim of confronting the obtained two-point correlation functions in one dimension with the exact ones, is to truncate the field dependence of Γ κ while preserving the complete dependence in both momentum and frequency of the two-point functions.

In the following, we hence proceed to a field truncation at the minimal order. It appears that the truncations of Γ κ at a given order in either ϕ or φ are not equivalent because the dependence in ϕ, contrary to that in φ, is constrained by the Galilean symmetry -ϕ enters in Dt . The minimal order in fields amounts to keeping terms, (i) at most quadratic in φ, (ii) linear in ∇ 2 ϕ or in D t ϕ, and (iii) to combine them with arbitrary powers of Laplacian ∇ 2 and covariant time derivatives Dt . This choice allows for an arbitrary dependence in momentum and frequency of the two-point functions. Notice that the dependence in ϕ is not restricted to be polynomial since arbitrary powers of this field are included through the covariant derivative Dt . To summarize, the ansatz for Γ κ considered in the following reads:

Γ κ [ϕ, φ] = x φf λ κ (-D2 t , -∇ 2 )D t ϕ -φf D κ (-D2 t , -∇ 2 ) φ - ν 2D ∇ 2 ϕf ν κ (-D2 t , -∇ 2 ) φ + φf ν κ (-D2 t , -∇ 2 )∇ 2 ϕ . (42)
Setting f λ κ = 1 and f ν κ = f D κ = D in the previous ansatz, one is left with the bare action (4) -written in terms of covariant time derivatives. Allowing non-constant (and scale-dependent) functions f X κ (X = D, ν, λ) of the Laplacian and of covariant derivatives is then the way to include an arbitrary dependence in momentum and frequency for the two-point functions that preserves all the symmetries, at the lowest (second) order in φ.

The operator Dt is squared in f D κ since this function is real. For f ν κ and f λ κ , to take the square of Dt is a choice that ensures a natural separation between real and imaginary parts of Γ (1,1) κ , without loss of generality for the two-point functions. Note that one needs to make precise the meaning of the expressions f X κ (-D2 t , -∇ 2 ) (X = D, ν, λ). Indeed, it is ambiguous as it stands since the Dt and ∇ operators do not commute. We take the convention that all Dt 's are on the left of the ∇'s. We also assume that these functions can be expanded in series of their arguments:

f X κ (-D2 t , -∇ 2 ) = ∞ m,n=0 a mn (-D2 t ) m (-∇ 2 ) n . ( 43
)
One can then consider the constraints stemming from the other symmetries. The shift-gauged symmetry imposes that f λ κ (ω 2 , p 2 = 0) ≡ 1. In dimension one, the time-reversal symmetry implies that the two functions f D κ and f ν κ become identical, and also that f λ κ (ω 2 , p 2 ) ≡ 1. In a generic dimension d, the ansatz (42) thus consists of three independent running functions of p 2 and ω 2 , and is reduced to a unique function of p 2 and ω 2 in d = 1.

To compute the flow equations of the functions f X κ (-D2 t , -∇ 2 ), X = λ, ν or D, one needs the expressions of the n-point functions up to n = 4. All the calculations are performed at vanishing fields. For the two-point functions, they are straightforward and yield

Γ (2,0) κ (ω, p) = 0 Γ (1,1) κ (ω, p) = iωf λ κ ω 2 , p 2 + ν D p 2 f ν κ (ω 2 , p 2 ) (44) Γ (0,2) κ (ω, p) = -2f D κ (ω 2 , p 2
). In this case, there is no ambiguity arising from the ordering of the gradient and covariant time derivative operators. The calculation and the explicit expressions of the three-and four-point functions, more lengthy, are detailed in Appendix C. One can check that all the Ward identities derived in Section III are satisfied by these functions.

We can be a priori confident in our approximation scheme, as it is close in spirit to the BMW approximation which can accurately capture the momentum dependence of two-point functions [41]. A more detailed discussion on the validity of our approximation scheme can be found in Appendix B.

V. SIMPLIFIED APPROXIMATION

The ansatz (42) in generic dimensions d can be further simplified by restraining the form of the functions f X κ (-D2 t , -∇ 2 ). The idea is to focus on a reliable description of the zero-frequency and zero-momentum sector of the theory while circumventing the limitations of the derivative expansion which is problematic here. For this, one can neglect all D2 t dependence in these functions:

f X κ (-D2 t , -∇ 2 ) → f X κ (-∇ 2
) for X = λ, ν and D. The ansatz proposed in [11] -which was derived before the formalization of the Galilean-gauged symmetry in terms of covariant derivatives -is actually endowed with an additional simplification since we imposed the one-dimensional constraint f λ κ = 1 in all dimensions (or stated otherwise we extended the condition f λ κ ( p = 0) = 1 holding for all d to the whole p sector). It reads

Γ κ [ϕ, φ] = x φ(x)D t ϕ(x) -φ(x)f D κ (-∇ 2 ) φ(x) - ν D φ(x)f ν κ (-∇ 2 )∇ 2 ϕ(x) . ( 45
)
As a matter of fact, this ansatz implies that all interaction vertices (Γ

(n)
κ with n > 2) are reduced to their bare form. Among these, the only non-vanishing one is Γ

(2,1) κ : Γ (2,1) κ (ω 1 , p 1 , ω 2 , p 2 ) = λ p 1 • p 2 . ( 46
)
With the ansatz ( 45), the two-point functions in generic dimensions simplify to

Γ (2,0) κ (ω, p) = 0 Γ (1,1) κ (ω, p) = iω + ν D p 2 f ν κ (p 2 ) (47) Γ (0,2) κ (ω, p) = -2f D κ (p 2
). With the aim of analyzing the fixed point structure, one introduces dimensionless and renormalized functions

f ν κ ≡ f ν κ /ν κ (respectively f D κ ≡ f D κ /D κ )
, where ν κ (respectively D κ ) are running coefficients (identifying at the microscopic scale κ = Λ with the bare parameters ν (respectively D) of the action (4)) that are related to anomalous dimensions κ∂ κ ν κ = -η ν (κ)ν κ (respectively κ∂ κ D κ = -η D (κ)D κ ). At a fixed point, these coefficients are expected to behave as power laws ν κ ∼ κ -η * ν and D κ ∼ κ -η * D and the scaling exponents are then expressed in terms of these anomalous dimensions as

z = 2-η * ν and χ = (2 -d + η * D -η * ν )/2.
The flow equations for the functions f ν κ and f D κ can be found in [11]. The Galilean invariance ensures that the flow of the dimensionless running coupling constant ĝκ ≡ κ d-2 λ 2 D κ /ν 3 κ is reduced to its dimensional part

∂ s ĝκ = ĝκ (d -2 + 3η ν (κ) -η D (κ)) (48)
with ∂ s = κ∂ κ , which thus enforces the identity z +χ = 2 at any non-Gaussian fixed point. The subleading exponent ω, governing the corrections to scaling, is independent of the leading critical exponents χ and z and can be [43,44]. See also Table I calculated from the flow behavior as s → -∞, using for instance that ĝκ ∼ ĝ * + ĝ1 exp(ωs) (49) in this limit. Despite its simplicity, this ansatz yields the correct phase diagram in all dimensions (Fig. 1(a)) -including the strong coupling fixed point, and reasonable exponent values in physical dimensions (Fig. 1(b)) [11]. More precisely, in all dimensions studied (i.e. up to d = 8), there exists besides the Edwards-Wilkinson fixed point F EW , a fully-attractive non-trivial strong-coupling fixed point F SC . In all d, generic scaling is found, i.e. the flow always reaches one of these fixed points. For d < 2, F SC is reached from any initial condition. For d > 2, F EW and F SC become locally fully attractive. More precisely, there exists a critical value of the (bare) coupling, denoted ĝc , that separates the basins of attraction of F SC and of F EW and for which the flow reaches the roughening transition fixed point F RT . In d = 2, F RT coincides with F EW , and becomes non-Gaussian for larger dimensions.

Because of Eq. ( 48) and since the fixed point coupling constant ĝ * SC is non zero, F SC is characterized by a single exponent. We thus only discuss below the values obtained for χ [11]. Table I and Fig. 1(b) contain our esti-mates for the roughness exponent χ and the subleading exponent ω. For d 4, χ decreases almost linearly with d, with the exact value 1 2 (respectively 2 3 ) recovered in d = 1 (respectively d = 0), and a reasonable but deteriorating agreement with numerical values for 2 ≤ d ≤ 4. In higher dimensions, χ increases with d, at odds with both numerical values and the scenario of d = 4 being an upper critical dimension beyond which χ = 0 [18,19,45]. The values for the critical exponents can be refined by computing the next order of our approximation scheme, which is work in progress.

Regarding F RT , we record negative values of χ for 2 < d < 5, which is reminiscent of perturbative results performed at fixed d [13] but in contradiction with exact results dictating χ = 0 [14,46]. This is to be imputed to our breaking of the Z 2 symmetry manifest in the Cole-Hopf formulation and non-linearly realized in the KPZ problem.

VI. ONE-DIMENSIONAL SCALING FUNCTION

We now consider the full ansatz (42) to calculate the momentum-and frequency-dependent two-point correlation function of the one-dimensional problem and to extract from it the universal scaling function. We recall that, due to the incidental time-reversal symmetry, the ansatz (42) simplifies in d = 1 to only one running function left:

f D κ = f ν κ ≡ f κ and f λ κ ≡ 1.
We also drop from now on the vector arrows and set ν = D = 1 since these two coefficients can be absorbed in the action (4) through an appropriate rescaling of the fields and of time and the change of coupling constant λ → √ g b = λD 1/2 /ν 3/2 .

A. Flow equations

The flow equation for the function f κ can be obtained either from the flow of Γ (0,2) κ or from the one of Γ (1,1) κ according to (44). These two flows are equal in d = 1 since the ansatz preserves all the symmetries including the time-reversal one. They can be computed using Eq. ( 14), where G κ is the propagator matrix

G κ = [Γ (2) κ + R κ ] -1
as defined in Eq. (13). Using the expressions (41) and (44)

(specialized to d = 1, i.e. setting ν κ = D κ , f D κ = f ν κ ≡ f κ and f λ κ ≡ 1)
for the matrix elements of R κ and Γ

(2) κ , one finds for the propagator

G κ (ω, q) = 1 P κ (ω 2 , q 2 ) 2k κ (ω 2 , q 2 ) Y κ (ω, q) Y * κ (ω, q) 0 ( 50
)
where k κ (ω 2 , q 2 ) = f κ (ω 2 , q 2 ) + D κ r(q 2 /κ 2 ), Y κ (ω, q) = iω + q 2 k κ (ω 2 , q 2 ) and P κ (ω 2 , q 2 ) = ω 2 + q 2 k κ 2 . The calculation of the flow equation for e.g. Γ (0,2) κ

, involving matrix products and trace according to Eq. ( 14), is then straightforward. This yields for the flow of f κ :

∂ κ f κ (̟, p) = - 1 2 q ∂ κ S κ (q) 1 P 2 κ (ω 2 , q 2 ) X κ (ω 2 , q 2 )Γ (2,2) κ (q, -q, p) + 1 2P 2 κ (ω 2 , q 2 )P κ ((̟ + ω) 2 , (p + q) 2 ) × 2X κ (ω 2 , q 2 )Y κ (̟ + ω, p + q)Γ (1,2) κ (-q, -p)Γ (2,1) κ (q, -p -q) -q 2 Y 2 κ (ω, q)Γ (1,2) κ (-p -q, p) × Y κ (̟ + ω, p + q)Γ (1,2) κ (-q, -p) + 2k κ ((̟ + ω) 2 , (p + q) 2 )Γ (2,1) κ (p + q, -q) -Y * κ (̟ + ω, p + q)Γ (1,2) κ (q, p) + 2k κ ((̟ + ω) 2 , (p + q) 2 )Γ (2,1) κ (q, -p -q) × q 2 (Y * κ ) 2 (ω, q)Γ (1,2) κ (p + q, -p) -2X κ (ω 2 , q 2 )Γ (2,1) κ (p + q, -q) ( 51
)
where q ≡ (ω, q), p ≡ (̟, p), X κ (ω 2 , q 2 ) = ω 2q 4 k 2 κ (ω 2 , q 2 ), S κ (q) = D κ r(q 2 /κ 2 ) and where the expressions for the three-and four-point vertex functions Γ (n,m) κ are given in Appendix C. The flow equation for f κ is hence an integral over the internal momentum q and frequency ω and depends on the external momentum p and frequency ̟.

As we intend to study the fixed point properties, we introduce dimensionless and renormalized quantities. Momenta are measured in units of κ, e.g. p = p/κ, and frequencies in units of D κ κ 2 , e.g. ̟ = ̟/(D κ κ 2 ). We define the dimensionless renormalized function fκ = f κ /D κ , as at the bare level, f κ=Λ = D = 1. The sole running anomalous dimension η κ is defined by κ∂ κ ln D κ = -η κ so that D κ ∼ κ -η * at the fixed point. The critical exponents in d = 1 are then given by z = 2 -η * and χ = η * . The absolute normalization of fκ and D κ is fixed by setting fκ (0, 0) = 1 for simplicity [41].

The flow equation for the dimensionless function reads:

∂ s fκ ( ̟, p) ≡ ∂ s 1 D κ f κ ̟ D κ κ 2 , p κ = η κ fκ ( ̟, p) + (2 -η κ ) ̟ ∂ ̟ fκ ( ̟, p) + p ∂ p fκ ( ̟, p) + 1 D κ κ∂ κ f κ (̟, p) (52)
where ∂ s = κ∂ κ . Once the substitutions for dimensionless quantities have been performed in (51), the last term 52) is dimensionless, and depends on the external dimensionless momentum p and frequency ̟. The dimensionless running coupling constant reads in one dimension ĝκ = κ -1 λ 2 /D 2 κ . Its flow equation is again reduced to its dimensional part due to Galilean invariance:

1 Dκ ∂ s f κ in (
∂ s ĝκ = ĝκ (2η κ -1) , (53)
and one finds as expected that χ = η * = 1 2 at any fixed point with ĝ * = 0. The subleading exponent ω is calculated according to (49) and we find values in agreement with the ones obtained with the simplified ansatz of Section V (see Tables I andII).

B. Numerical integration and error bars

The numerical integration of Eqs. ( 52) and (53) was performed using standard techniques and the numerical error stemming from the integration procedure was assessed by resorting to several system sizes, resolutions, etc. (see Appendix D). It turned out that the numerical error was negligible compared with the error coming from the approximation scheme itself, which constitutes the dominant source of inaccuracy and which we now discuss.

Ultimately, the error relative to the whole approximation scheme can be evaluated by comparing successive orders of approximation. However, when only a given order is available, one can assess the accuracy of this approximation by varying the parameter α of the cutoff function. Indeed, although physical quantities are obtained in the limit κ → 0 (where R κ=0 vanishes) from the exact NPRG flow equations ( 13) or ( 14) and thus do not depend in principles on the choice of the particular profile of R κ , any approximation introduces a spurious residual dependence on the regulator [30]. This sensitivity to the cutoff can be investigated and exploited to estimate the accuracy of the approximation implemented.

Here, α was varied between 1 and 60, and all the physical quantities computed and discussed in Section VI C exhibited a similar large plateau behavior. The value given in the following (text and Table II) for a quantity corresponds to the central value on the plateau, and the error bars reflect the dispersion of this quantity around its plateau value when α varies in the range [2,20]. We emphasize that this procedure only provides an intrinsic error at a given order of approximation related to an artificial regulator dependence. Accumulated experience in NPRG calculations seem to indicate that in general this procedure yields a reasonable estimate of the order of magnitude of the committed error [30]. However, in some cases, these error bars may represent only lower bounds on the uncertainties. Definite error bars can be inferred from the comparison with the next order of approximation.

Runs were started from various different initial conditions including the bare action ( fΛ ≡ 1) with ĝΛ typically between 1 and 10. In all cases, the function fκ smoothly deformed from its initial shape to eventually reach a fixed point where it stopped evolving, typically after s -20.

In other words, the NPRG flow equations entail generic scaling for the one-dimensional KPZ problem, since the long time and large distance behavior is always governed by this non-trivial fixed point, yielding a scaling regime (described in Section VI C). Fig. 2 shows a typical shape of the fixed point function f * ( ̟, p). It is fixed by the normalization condition to unity at vanishing momentum and frequency and decays as a power law when ̟ and/or p become large.

C. Scaling function

Extraction of the scaling function

We now turn to the description of the fully attractive fixed point. At this fixed point ∂ s fκ ( ̟, p) = 0 (in Eq. ( 52)) and g * = 0 such that η κ ≡ η * = 1 2 . Moreover, we verified numerically that the nonlinear term ∂ s f κ /D κ in (52) decouples (that is ∂ s f κ /D κ → 0) when ̟ ≫ 1 and/or p ≫ 1. As a consequence, at the fixed point, Eq. ( 52) reduces in the regime ̟ and/or p ≫ 1 to the homogeneous equation

1 2 f * ( ̟, p) + p ∂ p f * ( ̟, p) + 3 2 ̟ ∂ ̟ f * ( ̟, p) = 0. ( 54
)
Its general solution has the form

f * ( ̟, p) = 1 p1/2 ζ ̟ p3/2 (55)
where the function ζ cannot be determined from the homogeneous equation but can be extracted from the numerical solution of the full equation ( 52) -by tabulating the values p1/2 f * ( ̟, p) against the ratios ̟/p 3/2 . We observe that the fixed point function f * ( ̟, p) is regular for all ̟ and p (see Fig. 2), which allows us to deduce the limits of ζ. First ζ(0) has to be finite for the limit ̟ → 0 to exist and followingly f * ∼ ζ(0)p -1/2 for p ≫ 1 at fixed ̟. Second, ζ must behave as ζ(x) ∼ ζ∞ x -1/3 as x → ∞ for the limit p → 0 to exist and followingly f * ∼ ζ∞ ̟-1/3 for ̟ ≫ 1 at fixed p. We show below that the form of the solution (55) entails scaling for the (dimensionful) correlation function C(̟, p).

We now consider the fixed point dimensionful function

f (̟, p) = D κ f * ( ̟, p) = D κ f * ̟ D κ κ 2 , p κ . ( 56
)
The physical limit is obtained when the running scale κ tends to zero at fixed values of ̟ and p. This limit is precisely equivalent to ̟ ≫ 1 and/or p ≫ 1, which corresponds to the scaling regime described above where f * takes the form (55). Moreover, when κ → 0, D κ behaves as a power law D κ = D 0 κ -1/2 where D 0 is a non-universal constant. Hence, the physical dimensionful function f is expressed in terms of the scaling function ζ as

f (̟, p) = D 0 p 1/2 ζ 1 D 0 ̟ p 3/2 . ( 57
)
The correlation function is defined, upon inverting the matrix Γ (2) of the two-point vertex functions, by

C(̟, p) = - Γ (0,2) (̟, p) |Γ (1,1) (̟, p)| 2 = 2f (̟, p) ̟ 2 + p 4 f (̟, p) 2 (58)
where the ansatz (44), with f D κ = f ν κ ≡ f , f λ κ ≡ 1, and ν = D = 1, has been used in the second equality. Note that the definition (58) of the correlation function corresponds in real space to the connected mean value h(t, x)h(0, 0) c which differs by a factor -2 from the definition (2) in the co-moving frame. Restoring the ν and D coefficients also yields an additional factor, such that the correlation function defined by (58) finally relates to the Fourier transform of the correlation function defined by ( 2) via an overall multiplicative factor

C 0 = - 1 2 ν 2 D . ( 59
)
Replacing in (58) the dimensionful fixed point function f by its expression (57) yields multiplied by τ 7/3 to illustrate the power-law behavior of the tail F ∼ τ -7/3 (see text). Note that the small dispersion visible in the inset zooming the tail for τ 100 is hence magnified by a factor of order 100 7/3 .

C(̟, p) = 2 p 7/2 D 0 ζ 1 D0 ̟ p 3/2 ̟ 2 /p 3 + D 2 0 ζ2 1 D0 ̟ p 3/2 (60) = 2 p 7/2 1 D 0 ζ ̟ p3/2 ̟2 /p 3 + ζ2 ̟ p3/2 (61) ≡ 2 p 7/2 1 D 0 F 1 D 0 ̟ p 3/2 . ( 62
Eq. ( 60) shows that the correlation function takes a scaling form at large distance and long time. Let us emphasize once more that we did not assume scaling, the latter emerges from the fixed point solution of the flow equation starting from any reasonable microscopic initial condition. The scaling function F is hence universal, up to an absolute normalization and a rescaling of its argument by a multiplicative constant.

According to Eqs. ( 61) and ( 62), this function F can be computed from the numerical data for the dimensionless fixed function, selecting its values for arguments in the range ̟ and/or p ≫ 1 which corresponds to the scaling regime. We proceeded for various initial conditions and for different values of the parameter α. In all cases, we observed for each value of α the expected data collapse -generating the scaling function F -with a very high precision, as illustrated in Fig. 3.

Normalizations

Our aim is now to confront our scaling function F with the exact result obtained in Ref. [1]. The first step consists of fixing the normalizations. Indeed, two constants are involved in the standard definition (as recalled below Eq. ( 2)) of the universal scaling function g from the correlation function C C(τ, L) = const.τ 2χ/z g(const. ′ L 2χ τ -2χ/z ), that have to be fixed. The scaling function g(y) is normalized in [1] in the following way:

g(y) = lim t→∞ C 2λ 2 At 2 1/3 y, t 1 2 λA 2 t 2/3 (63)
with A ≡ D/ν. Then three functions f , f and f are introduced in [1], through the definitions

f (y) = 1 4 g ′′ (y) f (k) = 2 ∞ 0 dy cos(ky)f (y) (64) f (τ ) = 2 ∞ 0 dk cos(kτ ) f (k 2/3 )
and imposing the additional normalization condition

f (0) = 1. (65)
The last function f is proportional to C(ω, p)p 7/2 (see [1]), that is it corresponds to our function F reconstructed using (62) up to normalization factors. The analog of the functions f and f , that are computed in the following, will also be denoted with capital letters F and F respectively. The precise normalization constants between F and f can be established as follows. First one deduces from Eqs. ( 63) and (64) that:

f (τ ) = - p 7/2 2 5/3 λ 4/3 A 5/3 t 7/3 C p (2λ 2 At 2 ) 1/3 , τ p 3/2 t .
(66) Then, comparing this expression with our definition (62) and taking into account the multiplicative factor (59), we obtain the relation

f (τ ) = 2 2g * D νA F 2g * νA D τ (67)
where we used that the bare value λ is related to the fixed point coupling g * through g * = g b /D 2 0 = λ 2 D/(D 2 0 ν 3 ). Relating the Fourier transforms is then straightforward:

f (k) = 2 D νA F k 2g * νA D 1/3 ≡ 1 Fn F k k n . (68)
Finally, using the same normalization criterion (65) as in [1] gives the absolute vertical normalization Fn = νA 2D = F (0) and the absolute horizontal normalization k n = (4 Fn g * ) 1/3 .

Properties of the scaling function f

Let us first compare the scaling functions f and F . The function f is studied in detail in [1] and also in [18,22], where interesting features are highlighted. According to Ref. [1], the function f first decreases to vanish at k 0 ≃ 4.36236 . . . and then exhibits a negative dip of coordinates (k d ≃ 4.79079 . . . , fd ≃ -0.0012023 . . . ). After this dip, the function decays to zero with a stretched exponential tail, over which are superimposed tiny oscillations around zero, only apparent on a logarithmic scale. A heuristic fit of this behavior for k 15 is given in [1] f (k) ∼ 10.9k -9/4 sin

k 3/2 2 -1.937 e -1 2 k 3/2 . ( 69
)
We show below that we here recover qualitatively all these features, with reasonable estimates for the different parameters that characterize them. Following (64), the function F is defined from F by the integral

F (k) = ∞ 0 dτ π cos(τ k 3/2 ) F (τ ) (70)
which has to be computed numerically. The function F stems from the superposition of the numerous curves involved in the data collapse of Eq. ( 62). Although the data collapse is excellent, there subsists a high frequency and small amplitude noise due to the re-ordering of the points (apparent with a large zoom as presented in the insets of Fig. 3). As we intend to investigate the properties of the tail of F with high precision, we first need to eliminate this noise. For that, we devise an appropriate family of analytical fitting functions to smooth our data for F (τ ). The choice of the family of fitting functions is determined as follows.

First, they have to reproduce the large τ behavior of F . The latter can be inferred from the limits of the function ζ(x) established after Eq. (55). One obtains F (τ ) ∼ F∞ τ -7/3 as τ → ∞. For each value of the parameter α, the proportionality constant F∞ can be estimated graphically (up to 4-5 digits) from the curve τ 7/3 F (τ ). Moreover, the fitting functions have to be even, and finite at the origin. To satisfy these constraints, we build a family of fitting functions as an expansion in elementary rational polynomials of τ 2 raised to the adequate power to reproduce the large τ behavior:

Ffit (τ ) = a 00 + a 02 τ 2 1 + a 01 τ 2 + a 03 τ 4 + a 10 1 + a 11 τ 2 + a 20 1 + a 21 τ 2 + . . . 7/6 (71)
with the additional constraint (a 02 /a 03 + a 10 /a 11 + a 20 /a 21 + . . . ) 7/6 = F∞ [47]. The first fit is achieved using the three independent coefficients a 0i . The corresponding denominator introduces four non-analyticities in the points of the complex plane ±z ± 0 with FIG. 4. (Color online) Illustration of the fitting procedure described in the text, using the same data for F as in Fig. 3 but with unified colors (black curves: on the left F (τ ); on the right F (τ )τ 7/3 ). The four successive orders of the fit appear superimposed in the main graph; the two insets present large zooms (note the vertical scales) to differentiate them and highlight the convergence: the higher the order, the closer to the data (black-dotted curve) the fit lies. Note that for the right curve and inset F (τ )τ 7/3 the differences are again magnified by a factor τ 7/3 .

z ± 0 = 1 2a 03 -a 01 ± i -a 2 01 + 4a 03 (72)
whose coordinates were found to be very robust against the choice and the order of the fit as explained below.

Note that these four points appeared never to lie on the imaginary axis, nor on the real axis (-a 2 01 + 4a 03 > 0). We then add in turn simple monomials corresponding to purely imaginary poles up to nine independent coefficients (we observed that additional complex poles were systematically decomposed into two purely imaginary poles). This procedure turned out to converge rapidly for all α values, as illustrated in Fig. 4. In particular, the existence and the coordinates of the complex singularities ±z ± 0 appeared to be a robust feature of the data. Moreover, they are always found to lie the closest to the real axis when other poles are included according to (71) [48]. For the remaining calculations we work with the analytical expression Ffit which reproduces faithfully our data, and hence drop in the following the index fit .

We computed numerically the function F from F according to (70) and normalized it using (68), which yielded the following results. As displayed in Fig. 5, the overall agreement between the NPRG scaling function F and the exact one f of Ref. [1] is excellent. We find that F reproduces very accurately all the qualitative features of the exact function f , in particular the existence of the negative dip and the subsequent stretched exponential decay with the presence of oscillations in the tail. Indeed, we find for the position of the first zero k 0 ≃ 4.60 (6) and for the coordinates of the dip (k d ≃ 5.14 (6), Fd ≃ -0.0018(6)), which are close to the exact results, see Table II. As for the behavior of the tail of F , it can be inferred analytically from the pole structure of F . Let us define F (k) = F (k 2/3 ), which is hence the standard Fourier transform of F

F (k) = ∞ 0 dτ π cos(τ k) F (τ ). ( 73
)
As F (τ ) is C ∞ , the tail of F is dominated by the singularities of F (τ ) in the upper complex half plane lying the closest to the real axis, which are z + 0 and -z - 0 . Denoting z + 0 = a 0 + ib 0 with b 0 > 0, we obtain F (k) ∼ e -ikz + 0 + e ikz - 0 ∝ e -b0k cos(a 0 k) as k → ∞ and followingly

F (k) = F (k 3/2
) ∼ e -b0k 3/2 cos(a 0 k 3/2 ) as k → ∞.

(74) F hence decays following a stretched exponential with superimposed oscillations on the scale k 3/2 exactly as observed in the exact solution (see (69)). Regarding the tiny magnitude over which develop these features, this agreement is remarkable. Let us emphasize that within the MC approximation and the self-consistent expansion, the stretched exponential behavior and the oscillations were also highlighted, but not on the correct scale [1,18,22]. The NPRG method hence seems to provide more accurate results than other nonperturbative approaches, and with much less inputs (no assumed scaling nor a precise scaling form).

A consequence is that, conversely to the MC approach for instance, the coefficient of the exponential and the period of the oscillations can here be estimated -extracted analytically from the values of z 0 . We find b 0 ≃ 0.49 (1) for the coefficient of the exponential, and a 0 ≃ 0.28 (5) for the pulsation of the oscillations (checking that these values coincide with the same quantities estimated graphically from the tail of F computed numerically -the latter being determined with much less accuracy). These coefficients can be compared with the coefficients of the fit (69) of the tail of the exact function which are b 0 = a 0 = 1 2 . The NPRG coefficient b 0 essentially matches the exact one and the pulsation a 0 is of the same order as the exact one, which is already highly non-trivial. The discrepancy is only visible on a logarithmic scale, as illustrated in the inset of Fig. 5. Let us emphasize once more that the error bars attributed here to a given quantity reflect the weak variations of this quantity around plateau values when the α parameter of the cutoff function is varied. It represents an error estimate intrinsic to the order of approximation under study and does not imply that the value at the next order of approximation would necessarily fall within these error bars.

Scaling function in real space and universal amplitude ratio

We finally come to the real space scaling function f defined in (64). We computed it numerically by Fourier transforming F

F (y) = ∞ 0 dk π cos(ky) F (k). (75)
It is compared in Fig. 6 with the exact result. Again, the overall agreement is manifestly excellent. The NPRG function F reproduces very precisely the exact one f , though it is the less accurate of our three scaling functions since it stems from the two successive numerical (oscillating) integrations ( 70) and (75) of our raw data.

The tail is particularly sensitive to this loss of precision. The exact function f (y) is found to decrease as exp(-cy 3 ) when y → ∞ in [1], whereas the decay of the function F (y), though it starts with the correct behavior, rapidly crosses over to a simple exponential decay exp(-c ′ y).

From this function can be computed the universal amplitude ratio g 0 with the definition given in [1]

g 0 ≡ 4 ∞ 0 dy yf (y). ( 76
)
The exact result is the Baik-Rain constant [23] g 0 = 1.1503944783... (more digits can be found in [1]). Here, the universal constant g 0 can be estimated by performing the numerical integration corresponding to (76) using F . However, this amounts to achieving three successive numerical integrations from our raw function F and the resulting precision is low. Alternatively, one can compute part of the involved integrals analytically, using the definition given in [1]: g 0 ≡ g(0) where g(y) is the original scaling function proportional to the second derivative of f (y) (see (64)). Indeed, g(y) can be expressed integrating twice f through its Fourier transform f

g(y) -g(0) = 4 π y 0 dv v 0 du ∞ 0 dk cos(ku) f (k) = 4 π ∞ 0 dk (1 -cos(ky)) f (k) k 2 . ( 77
)
Taking the limit y → ∞ in this expression with the change of variables z = ky (and recalling that f (k = 0) = 1), one deduces that in this limit:

g(y) = 4 π ∞ 0 dk 1 -cos(ky) k 2 + O(1/y), (78)
and consequently

g(0) = lim y→∞ - 4 π ∞ 0 dk (1 -cos(ky)) f (k) -1 k 2 .
(79) Since the function ( f (k)-1)/k 2 is infinitely derivable, its Fourier transform vanishes in the limit y → ∞ and one is left with

g(0) = - 4 π ∞ 0 dk f (k) -1 k 2 (80) = 6 π 2 ∞ 0 dk √ k ∞ 0 dτ sin(τ k 3/2 )τ f (τ ) (81) = 2 π 2 Γ 1 3 ∞ 0 dτ τ 2/3 f (τ ), (82)
using the definitions (64). This universal quantity can thus be determined numerically with much higher precision from F by performing a unique, non-oscillating integral. We find g 0 ≃ 1.19 (1) which is in close agreement with the exact value. The defintion (76) suggests that the slight over-estimation is to be imputed to the contribution of the tail of the NPRG function F , which only decays exponentially, whereas it is suppressed as exp(-cy 3 ) in the exact function f . Note that the MC estimate g 0 ≃ 1.1137 [18] is of comparable accuracy with the NPRG result. However, the discrepancies between the MCT scaling function and the exact one appear much larger (see [18]), such that the estimate of g 0 probably benefits from some compensation between the bulk and the tail of the MCT function. All our estimates for the characteristic parameters of the different scaling functions obtained in this Section are summarized in Table II.

VII. CONCLUSION

In this work, we have presented a general framework to investigate the KPZ equation using the NPRG method. We have proposed a detailed and revisited analysis of the symmetries of the KPZ equation, both the standard symmetries -Galilean and shift -and their versions gauged in time, as well as the one-dimensional time-reversal symmetry and a nonlinearly realized Z 2 symmetry. We have derived general Ward identities associated with these symmetries (except for the Z 2 one) and proposed a convenient geometric interpretation of the Galilean symmetry in terms of covariant time derivatives.

We have then devised an approximation scheme based on an ansatz for Γ κ which, in contrast to the standard derivative expansion often used in the NPRG framework, preserves the momentum and frequency structure of the vertex functions. This ansatz, constrained by the symmetries, allows one to compute correlation functions in a nonperturbative, yet systematically improvable, way. We explicitly implemented the minimal order in the response field of this approximation scheme, specifically developing two applications: the determination of critical exponents in physical dimensions (with a simplified ansatz) and the computation of universal scaling func-tions in one dimension.

We found that, without any input other than the bare action and its symmetries (in particular without assuming scaling), the renormalization group flow is generically attracted toward a strong-coupling fixed point -which roots the existence of generic scaling -or to the EW fixed point for d > 2 and λ < λ c , yielding for the first time the full correct phase diagram of the KPZ equation within a RG approach.

The estimates of the critical exponents obtained in dimensions two and three compare reasonably with results from simulations, though the accuracy of these exponents decreases with the dimension. We are confident that these results, which were obtained using what remains, finally, a rather crude approximation, will get better as higher order approximations are studied -which is hopefully feasible in practice. Indeed, we know from experience that within the NPRG framework, the convergence is generally very fast, such that pushing it to even a slightly higher order can eventually provide excellent results [30]. Similarly, we are unable, at this current order of approximation, to settle on the existence of an upper critical dimension, but we should be able to solve this issue at next orders.

The scaling functions F , F and F obtained in one dimension, on the other hand, compare very accurately with the exact ones. In particular, we recover the stretched exponential with tiny superimposed oscillations of the tail of f (k) on the correct scale k 3/2 and with the correct coefficient b 0 ≃ 1/2 for the exponential decay. We also obtained an accurate estimate of the universal amplitude ratio g 0 ≃ 1.19 (1). Quantities beyond the leading scaling regime can also be computed, which we illustrated here with the determination of the correction to scaling exponent ω.

A further step will consist of refining the computation of critical exponents in all dimensions by improving the approximation, which is in progress. One could also compute the probability distributions of the h(t, x) field, which are again known exactly in one dimension but not in any other, and investigate the influence of boundary conditions.

where 'perm.' refers to the following: on the left-hand side, the permutation of the last (implicit) frequencymomentum argument (corresponding to the φ leg) together with the p 2 factor in front, from ω 1 over to ω 2 , ω 3 and -ω 1 -ω 2 -ω 3 ; on the right-hand side, the permutation of the first frequency-momentum argument (corresponding to the ϕ leg) together with the complementary product of p 2 in front, from -ω 1 -ω 2 -ω 3 over to ω 1 , ω 2 and ω 3 .

APPENDIX B: VALIDITY OF THE APPROXIMATION SCHEME

The approximation scheme presented in Section IV is closely inspired by the BMW scheme, which provides an accurate description of the momentum dependence of two-point functions at equilibrium [41]. The BMW scheme thus stands as the most adequate approximation scheme for the KPZ problem, both because we are intrinsically interested in the momentum and frequency structure of the two-point functions and because the derivative nature of the interaction vertex seems to require one to go beyond a simple derivative expansion [39]. However, the implementation of the BMW scheme for the KPZ problem is hindered by the very demanding symmetries of the KPZ action. As a consequence, it cannot be directly performed but has to be adapted regarding three main aspects.

First, the standard BMW approximation at equilibrium consists, at leading non trivial order, of performing truncations on the three-and four-point functions, keeping only the leading order in their internal momenta since high momentum contributions are suppressed by the cutoff function (see [41] for a full justification of the scheme). In the KPZ problem, as stressed in Section IV, the Galilean symmetry relates n-point functions with different n and different momenta in a complicated way such that a naive expansion of the three-and fourpoint functions in their internal momenta and frequencies would spoil these relations. To overcome this difficulty, we propose here instead to devise an ansatz for Γ κ that manifestly satisfies all the symmetries and such that they are hence automatically preserved in subsequent calculations (functional derivatives with respect to the fields). This ansatz imposes in turn the form of the three-and four-point functions. It appears that these functions retain some internal momentum and frequency dependencies that would be neglected in a direct implementation of the BMW approximation but which are here necessary to fulfill the symmetry constraints. Hopefully, these additional subleading terms are not expected to deteriorate the quality of the approximation.

Second, another specificity inherent to the KPZ problem resides in the absence of regulator on the frequency sector, which is prevented by the Galilean-gauged symmetry (see Section III D). As the whole BMW scheme relies on the presence of the cutoff term to perform expansions in internal momenta and frequencies, the absence of such a cutoff in the frequency sector would render an expansion in internal frequencies unjustified. Hopefully, the very same symmetry also cures this problem. Indeed, if vertices are expanded in momenta, then the Galileangauged symmetry fixes the frequency dependence associated with ϕ-legs (see e.g. Eq. ( 33)). Accordingly, no internal frequency expansion is needed for the ϕ-legs once an expansion in internal momenta is performed.

Finally, the previous discussion applies to the ϕ-legs. However, the φ-legs are treated in the proposed ansatz (42) in a much cruder way since Γ κ is simply expanded at quadratic order in that field. The quality of such an expansion has been analyzed in equilibrium statistical mechanics both within the derivative expansion (see, e.g. [29,30]) and within the BMW approximation scheme [42], and it seems to give reasonable results. Nevertheless, it probably constitutes the main source of inaccuracy of the present scheme. We leave for future work the improvement of the φ dependence of the approximation.

APPENDIX C: n-POINT FUNCTIONS AT THE MINIMAL ORDER OF THE APPROXIMATION SCHEME

To compute the three-and four-point functions at vanishing field ensuing from the ansatz (42) involves to take functional derivatives of the functions f X κ (-D2 t , -∇ 2 ) (X = D, ν, λ) with respect to the field ϕ and then to evaluate the result at ϕ = 0. This can be done using the expression of the series expansion (43):

δf X κ (-D2 t , -∇ 2 ) δϕ(t 1 , x 1 ) = ∞ m=1,n=0
ble. In all cases, a fixed point is reached after s -20.

The fixed point functions are recorded at s = -25. We studied separately the influence of the resolution (∆p and ∆ √ ̟), and of the mesh sizes (p max and ̟max ) on the precision level, and checked the convergence. The dif-ferences (in all physical quantities) between resolutions and/or domain sizes are in all cases smaller than typically 1% and these numerical errors are dominated by the residual variations observed when varying the cutoff parameter α.

[47] This is done to improve the quality of the fit, otherwise the relative weight of the tail is very small in the fitting procedure.

[48] The only exception to ±z ± 0 lying the closest to the real axis comes from spurious poles of the rational fraction inside the power law in (71) that are also zeros of the numerator and therefore do not contribute. Those artifacts appear ponctually and do not persist from an order of the fit to the other.

FIG. 1 . 1 d= 2

 112 FIG. 1. (Color online) (a) Flow diagram stemming from our (simplified) minimal order approximation in the (d, ĝv d ) plane (with v -1 d = 2 d+1 π d/2 Γ[ d 2 ] a normalization factor related to the integration volume). Red circles: renormalized value ĝ * SC at FSC. Dashed purple line: Gaussian fixed point FEW. Cyan squares: bare value ĝc, which separates the basins of attraction of FSC and FEW. Gray arrows symbolize flow lines. (b) Variation with d of χ = 2 -z for FSC (red circles: our results; orange squares: numerical values from[43,44]. See also TableI.)

FIG. 2 .

 2 FIG. 2. (Color online) Typical shape of the dimensionless fixed point function f * ( ̟, p) (recorded at s = -25 and for α = 1).

) 3 FIG. 3 .

 33 FIG. 3. (Color online) Leftmost curve: Scaling function F ̟ p3/2 corresponding to the data collapse (62), with the two insets zooming different parts to show the very high quality of the collapse (note the vertical scales); the various colors (gray shades) differentiate the contributions of distinct values of ̟. Rightmost curve: Scaling function F τ = ̟ p7/3

FIG. 5 .

 5 FIG. 5. (Color online)Comparison of the scaling function F (k) (red curve with dots) obtained in this work with the exact one f (k) (black curve with squares) from[1]. The inset shows the stretched exponential behavior of the tail with the superimposed oscillations, developing on the same scale k 3/2 . Note the vertical scale: this behavior develops with amplitudes below typically 10 -6 (see text for a detailed comment of the figure).

FIG. 6 .

 6 FIG. 6. (Color online) Comparison of the scaling function F (y) (red curve with dots) obtained in this work with the exact one f (y) (black curve with squares) from [1].

TABLE I .

 I .) Exponent values in integer dimensions. The average numerical values are extracted from[43]. To our knowledge, no estimates of ω are available in the literature.

	d	1	2	3	4
	χ [38]	0.50	0.33	0.17	0.075
	χ (numerics)	0.50	0.38	0.30	0.24
	ω [38]	0.817	0.70	0.63	0.54

  ) δϕ(t 1 , x 1 ) = λ ∇ x (δ(t -t 1 )δ (d) ( x -x 1 )) • ∇ x (∂ t -λ∇ x ϕ(x) • ∇ x ) +(∂ t -λ∇ x ϕ(x) • ∇ x )(∇ x (δ(t -t 1 )δ (d) ( x -x 1 )) • ∇ x ) . (87) , p 1 , ω 2 , p 2 ) = iλ [ p 1 • ( p 1 + p 2 )(ω 1 + 2ω 2 )] (ω 1 + ω 2 ) 2(m-k-1) (( p 1 + p 2 ) 2 ) n + (2 ↔ 3). , p 1 , ω 2 , p 2 ) = iλ ω 1 p 1 • ( p 1 + p 2 ) f D κ ((ω 1 + ω 2 ) 2 , ( p 1 + p 2 ) 2 ) -f D κ (ω 2 2 , ( p 1 + p 2 ) 2 ) -p 1 • p 2 f D κ (ω 2 2 , p 2 2 ) -f D κ ((ω 1 + ω 2 ) 2 , p 2 2 ) . (90)All the other functions can be computed the same way. The three-point functions read , ω 2 ;p 1 , p 2 ) = λ p 1 . p 2 f λ κ (ω 1 + ω 2 ) 2 , ( p 1 + p 2 ) 2 + λ p 1 . p 2 (ω 2 , ω 1 ; p 2 , p 1 ) + p 2 2 Γ (1,2)κ,ν (ω 1 , ω 2 ; p 1 , p 2 ) .

	with
	δ(-D2
	When Fourier transforming,
								δ(-D2 t ) δϕ(t 1 , x 1 )	→ -iλ p 1 • ( p 1 + p)(ω 1 + 2ω).	(88)
	One obtains for instance
	Γ (1,2) κ	(ω 1 ∞ m=1,n=0	a mn ×	m-1 k=0	ω 2k 2 (89)
	and thus:
	Γ (1,2) κ (ω 1 Γ (3,0) κ (ω 1 , ω 2 ; p 1 , p 2 ) = 0							(91)
	Γ (0,3)
		+ λ p 1 . p 2	ω 2 ω 1	ω 1 ω 2 2 , p 2 κ ω 2 2 -f λ κ (ω 1 + ω 2 ) 2 , p 2 f λ 2	f λ κ (ω 1 + ω 2 ) 2 , p 2 1 -f λ κ ω 2 1 , p 2 1
		-iλ	ν 2D		p 2 1	p 2 ω 2	.( p 2
		+ p 2 1		p 2 ω 2	. p 1 f ν κ (ω 1 + ω 2 )	2 , p 2 1 -f ν κ ω 2 1 , p 2 1	+ p 2 2	p 1 ω 1	. p 2 f ν κ (ω 1 + ω 2 )	2 , p 2 2 -f ν κ ω 2 2 , p 2 2
		+ p 2 2	p 1 ω 1	.( p 2	.	(93)
	Note that the combinations appearing in Γ	(2,1) κ	and Γ	(1,2) κ	are related. If one denotes Γ	(1,2) κ,D the expression (90) where
	the D index labels the function (f D κ ) appearing on the right-hand side, and similarly Γ where the function f D κ is replaced by f ν κ , one has	(1,2) κ,ν the same expression (90)
		iImΓ (2,1)						ν 2D	p 2 1 Γ (1,2) κ,ν
									a mn	m-1 k=0	(-D2 t ) k δ(-D2 t ) δϕ(t 1 , x 1 )	(-D2 t ) m-k-1 (-∇ 2 x ) n	(86)

t κ (ω 1 , ω 2 ; p 1 , p 2 ) = 0 (92) Γ (2,1) κ (ω 1 1 + p 2 ) f ν κ (ω 1 + ω 2 ) 2 , ( p 1 + p 2 ) 2 -f ν κ ω 2 1 , ( p 1 + p 2 ) 1 + p 2 ) f ν κ (ω 1 + ω 2 ) 2 , ( p 1 + p 2 ) 2 -f ν κ ω 2 2 , ( p 1 + p 2 ) κ

(ω 1 , ω 2 ; p 1 , p 2 ) = -

ACKNOWLEDGMENTS

L. Canet thanks the Universitad de la República for funding and hospitality during important stages of this work. The authors also wish to thank R. Blanch for his help to optimize the numerical codes and parallelize them, M. Prähofer and H. Spohn for providing us with their data for the scaling function f and its Fourier transform f , and M. Moore for discussions concerning the MC results. N. Wschebor acknowledges the support of PEDECIBA and ANII-FCE. The numerical parallel codes were run on the clusters FING and PEDECIBA (Montevideo) and on the cluster Healthphy (CIMENT, Grenoble).

APPENDIX A: TIME-REVERSAL WARD IDENTITIES FOR THE THREE-AND FOUR-POINT FUNCTIONS

All the Ward identities for the three-and four-point vertex functions are derived by taking functional derivatives of the identity (37) and evaluating them at uniform and static fields. For the remaining three-point function, one obtains in Fourier space:

For the four-point functions, we report the four additional independent identities obtained for uniform and static fields:

and finally

This type of relation is used in the following to shorten the expression of Γ

(3,1) κ

. The four-point functions are

where in the last line the D index refers to the involved function f D κ . We denote Γ(2,2) κ,D the two last lines of (95) with the same convention for the meaning of the index D. Finally, one finds, omitting the dependence in p i (which follows the same order as ω i ):

This appendix is devoted to presenting the details of the numerical procedure implemented to achieve the integration of flow equations, on the example of Eqs. ( 52) and (53). The momentum and frequency are discretized on a p × √ ̟ mesh of spacing ∆p and ∆ √ ̟ and sizes pmax and √ ̟max . The integrals over the internal momentum q and frequency ω in (51) are performed using Simpson's rule. The √ ̟ mesh is chosen because the integrand can have rather long tails in ω as there is no cutoff term (analogous to r(y)) suppressing the high frequencies. The values of f at ω ± ̟ which do not fall on √ ̟ mesh points are evaluated using cubic interpolations. The function f is extended outside the grid (for momenta p + q greater than pmax and/or frequencies ̟ + ω greater than ̟max ) using power law extrapolations.

Regarding the integration over q, as the integrand falls off exponentially due to the (derivative of the regulator) ∂ s S κ term in (51), the bounds ±∞ of the integral can be safely replaced by ±p max . This is not the case for the integral over the frequency ω as the decay of the integrand may be slow. We first compute the integral on [-̟max , ̟max ] and then evaluate the contribution of the integral on the boundaries [-∞, -̟max ] and [ ̟max , ∞] by performing in these regions the change of variable x = ω/ω max and using the extrapolated values of f . The precision on the momentum and frequency integral is of order 10 -4 for the typical resolutions ∆p = ∆ √ ̟ = 1/4 and mesh sizes ̟max = 15 2 and pmax from 20 to 45 increasing with the value of α (20 for α = 0.5 to 45 for α = 14). The derivative terms p∂ p and ̟∂ ̟ are computed using 5-point differences.

We use an explicit Euler time stepping with a typical time step ∆s = -2.10 -5 to integrate the flow equations on the renormalization time s, which turns out to be sta-