Critical line of the $\Phi^4$ scalar field theory on a 4D cubic lattice in the local potential approximation

We establish the critical line of the one-component $\Phi^4$ (or Landau-Ginzburg) model on a simple four dimensional cubic lattice. Our study is performed in the framework of the non-perturbative renormalization group in the local potential approximation with a soft infra-red regulator. The transition is found to be of second order even in the Gaussian limit where first order would be expected according to some recent theoretical predictions.


Introduction
It is a real pleasure and a great honor for the author to contribute, with this paper, to the festschrift dedicated to Professor Myroslav Holovko at the occasion of his 70t h birthday. Myroslav is an expert of the collective variables (CV) method introduced by the Ukrainian school in the framework of which Wilson's ideas on the renormalization group (RG) [1] can be implemented with great effect [2]. Here we expose recent post-Wilsonian advances on the RG in the framework of statistical field theory. Obviously, many of the ideas exposed here could easily be transposed to the CV "world" by the readers of references [3,4] where the links between the CV method and standard statistical field theory are established.
These last past years, Wilson's approach to the RG [1,5] has been the subject of a revival in both statistical physics and quantum field theory. Since the seminal work of Wilson, two main formulations of the non-perturbative renormalization group (NPRG) have been developed in parallel. Very similar to the works of the Ukrainian school on the of CV formalism we have the approaches initiated independently and in parallel by Wetterich et al. [6][7][8][9] on the one hand and Parola et al. in the other hand [10][11][12]. In this corpus of works one is interested to establish and to solve the flow equations of the Gibb's free energy by means of non-perturbative methods. In an alternative formulation, Polchinski and his followers consider rather the flow of the Wilsonian action [13,14], instead of that of the free energy, which makes the method more abstract and less predictive than that of Wetterich, although more in accord with Wilson's ideas. The link between these two formulations can however be established, see for instance references [15,16]. Other non-perturbative methods based either on the CV or Monte Carlo methods are also the subject of active studies and are discussed, for instance, in reference [17] and references quoted herein.
The NPRG has proved its ability to describe both universal and non universal quantities for various models of statistical and condensed matter physics near or even far from criticality. Recently it has been extended to models defined on a lattice [18]. Successful applications to the three-dimensional (3D) Ising, XY, Heisenberg models [19] and Φ 4 model [20] are noteworthy. Here we extend the study of reference [20] on the Φ 4 model in three dimensions of space to the case D = 4; it was made possible by the recent publication by Loh of a novel numerical method to compute the lattice Green's functions [21]. The D = 4 find available Monte Carlo simulations to which to compare our data. We stress that, in the wide range of parameters considered in our study (see table 1), we exclude the occurrence of a first order transition. This conclusion seems in agreement with a general analysis of the criticality of the model made in reference [29,30].
Our paper is organized a follows : In section 2 we review briefly the basic definitions and results concerning the statistical mechanics of scalar fields on a lattice. Section 3 is devoted to theoretical and technical aspects of the NPRG on the lattice. We then present our numerical experiments and discussed the results in section 4. We conclude in section 5

Model
Let us consider some arbitrary field theory defined on a 4D hyper-cubic lattice Λ = aZ 4 = {r|r µ /a ∈ Z; µ = 1, . . . , 4} where a is the lattice constant. The real, scalar field ϕ r is defined on each point of the lattice. It is convenient to start with a finite hyper-cubic subset of points {r} ⊂ Λ and to assume periodic boundary conditions (PBC) for the ϕ r before taking the infinite volume limit, although no difficulties are expected to arise from this operation.
In the case of short-range interactions between the fields, the action of the theory can quite generally be written as [22] where ϕ is a shortcut notation for {ϕ r } and is the Fourier transform of the field and the N momenta q are restricted to the first Brillouin zone B = [−π/a, π/a] ⊗4 of the reciprocal lattice. The inverse transformation reads :  Obviously one has ǫ 0 q ∼ q 2 for q → 0 and max q ǫ 0 q = ǫ max 0 = 16/a 2 . We will also define for convenience k max ≡ 4/a by ǫ max Note that in a system of units where the dimension of wave-vector q µ is [q µ ] = +1, the dimension of the fields are [ϕ r ] = 1 and [ ϕ q ] = −3 so that the kinematic part of the action S ϕ is dimensionless. Henceforth we shall only consider the Landau-Ginzburg polynomial form U (ϕ) = (r /2) ϕ 2 + (g /4!) ϕ 4 . Since [ϕ r ] = 1 and [a 4 U (ϕ)] = 0 it follows that [r ] = 2 and [g ] = 0. Therefore, in the thermodynamic limit, the physics of the model depends only upon the two dimensionless parameters r = r a 2 and the dimensionless (only in D = 4 ) g = g .
Another way of writing the action (2.2), which is useful for numerical investigations, is [22,28] where the 4 unit vectors e µ constitute an orthogonal basis set for R 4 . The field ψ and the parameters (κ, λ) are all dimensionless and they are related to the bare field ϕ and dimensionless parameters (r , g ) through the relations (2.7c)

Thermodynamic and correlation functions
The thermodynamic and structural properties of the model are coded in the partition function [31] Z [h] = Dϕ exp −S ϕ + h|ϕ , (2.8) where the dimensionless functional measure is given by Dϕ = n dψ n , (2.9) where r = an, the dimensionless ψ n is defined at equation (2.7a), h is an external lattice field, and the dimensionless scalar product in (2.8) is defined as h|ϕ = a 4 r h r ϕ r . (2.10) The order parameter is given by ∂h r , (2.11) where the brackets 〈· · · 〉 denote statistical ensemble averages and the Helmholtz free energy . Note that in the continuous limit, i.e. L = N a fixed, a → 0, the partial derivatives tend to functional derivatives, i.e. a −4 ∂ · · · /∂h r → δ · · · /δh(r). It follows from first principles that W [h] is a convex function of the N variables {h r }; it is also the generator of the connected correlation functions G (n) (r 1 . . . r n ) = a −4n ∂ n W [h] /∂h r 1 · · · ∂h r n , where ∂ · · · /∂h r denotes a partial derivative with respect to one of the N variables h r .

Lattice NPRG
An elegant procedure to implement the lattice NPRG was given by Dupuis et al. in references [18,19]; it extends to the lattice the ideas of Wetterich [6,7] for the continuum, i. e. the limit a → 0 of the model; it is very similar to the Reatto and Parola hierarchical reference theory of liquids [10][11][12]. We add a quadratic term to the action (2.2) where R k q is positive-definite, has the dimension [ R k ] = 2 and acts as a q dependent mass term. The regulator R k q is chosen in such a way that it acts as an infra-red (IR) cut-off which leaves the highmomentum modes unaffected and gives a mass to the low-energy ones. Roughly R k q ∼ 0 for ||q|| > k and R k q ∼ Z k k 2 for ||q|| < k. The scale k in momentum space varies from Λ ∼ a −1 , some undefined microscopic scale of the model yet to be defined precisely, to k = 0 the macroscopic scale. To each scale "k" corresponds a k-system defined by its microscopic action S k ϕ = S ϕ + ∆S k ϕ . We denote its partition function by Z k [h], its Gibbs free energy byΓ k φ , etc. The generalization of equation (2.13) is where the so-called average effective action Γ k φ , which was introduced by Wetterich in the first stages of the NPRG, is defined as a modified Legendre transform of W k [h] which includes the explicit subtraction of ∆S k φ [6,7], i. e.
Note that the functional Γ k φ is not necessarily a convex functional of the classical field φ by contrast withΓ φ which is the true Gibbs free energy of the k-system. The choice of the regulator R k q would not affect exact results but matters as soon as approximations are introduced. We have retained the Litim-Dupuis-Machado (LMD) regulator introduced by Dupuis and Machado [18,19] for the lattice as an extension of Litim's regulator widely used for off-lattice field theories [23]. Sharp cut-off regulators often yield unphysical behaviors, notably in the local potential approximation, and should be avoided, see e. g. [20,25]. The LMD regulator reads where ǫ k = k 2 and Θ is the Heavyside's step function. At scale "k", the effective spectrum of the k−model We note that for ǫ 0 q > ǫ k the regulator R k q vanishes in agreement with the fact that the high energy modes are of affected, i. e. one has ǫ eff. k (q) = ǫ 0 q . Conversely, for ǫ 0 q < ǫ k , a constant massive contribution is associated to the low-energy modes, with a tendency to a freezing of their fluctuations, i. e. one has ǫ eff.
It is easy to show the the average effective action satisfies the exact flow equation [6-8, 18, 19] where G (2) k is the Fourier transform of connected pair correlation function of the k−system defined as For an homogeneous configuration of the field φ r = φ we have, on the one hand, Γ k φ = N a 4 U k (φ) where the potential U k (φ) is a simple function of the field φ and, on the other hand, the conservation of momentum at each vertex which implies, with the usual abusive notation, G (2) k (q, −q) = N a 4 G (2) k (q); from these remarks it follows that : where the second line (3.8b) is valid in the thermodynamic limit (a fixed, N → ∞). Note that in order to establish equation (3.8) we also took into account of equation (2.14) in Fourier space for the k−system, , for an homogeneous system. The reader will agree that equation. (3.8), which is exact, is an extremely complicated equation since the vertex function Γ (2) k (q, −q), which is the Fourier transform of the second-order functional derivative of Γ φ with respect to the classical field φ, depends functionally upon φ.
The implicit solution (3.2) of (3.8) allows us to establish precisely the initial conditions. The initial Physically it means that all fluctuations are frozen and the mean-field theory becomes exact. When the running momentum goes from k = Λ to k = 0 all the modes ϕ q are integrated out progressively and the effective average action evolves from its microscopic limit Γ Λ

Local models and the initial condition of the flow
Some members of our family of k-systems are nice fellows. It follows from (3.5) that, for Λ > k > k max , or equivalently ǫ k > ǫ 0 q for all vectors q of the first Brillouin zone, we have ǫ eff. k ≡ ǫ k which means that the action S k ϕ of the k-system is local and reads S k ϕ = a 4 {r} U (ϕ r ) + (1/2) ǫ k ϕ 2 r . Therefore, at scale "k", we have a theory of independent fields on a lattice, which is trivial.
where we have introduced the dimensionless variables ϕ = aϕ, h = a 3 h, and ǫ k = a 2 ǫ k . Note that U (ϕ) = a 4 U (ϕ). The Helmholtz free energy and Wetterich effective action can be written as lattice sums where the convex functions ln z k (h) and γ k (φ) are related by a Legendre transform γ k (φ)+ln z k (h) = φ h, with, for instance φ = d ln z k (h)/d h. In general, the quantities ln z k (h) and γ k (φ) cannot be computed analytically but can easily be evaluated numerically for any value of Λ > k > k max .
It is interesting to note that the implicit equation (3.2) now reads which leads us to two remarks. First, the choice Λ = ∞ implies γ Λ = U since we can replace the Gaussian Our second remark is that one can derive from the equation (3.11), i.e. from his solution!, the flow Noting that, for a homogeneous system, with ∂ t = −k∂ k . Clearly, equation (3.13) can also be obtained directly from (3.8) in the range Λ > k > k max .
We are now in position to explicate the initial conditions which can be used to solve the flow equation (3.8) for the local potential • either Λ = ∞ and U Λ = U (Mean field theory like initial conditions). In this case the flow equation (3.13) must be solved numerically for Λ > k > k max . Note that U Λ (φ) can be non-convex.
In our numerical experiments we retained the second term of the alternative.

The general case
A non-perturbative, but intuitive approximation to solve the flow equation. (3.8) is to make an ansatz on the functional form of Γ k [φ]. In the local potential approximation (LPA) one neglects the renormalization of the spectrum and assume that [18,19] (LPA ansatz) (3.14) For a uniform configuration of the classical field φ r = φ and, in the thermodynamic limit, the flow equation (3.8b) becomes :

The LMD regulator
With the LMD regulator (3.4) the loop-integral in the r.h.s. of equation (3.15) can be worked out analytically which leaves us with a much simplified flow equation for the potential is the threshold function [7] which takes a very simple expression with the LMD regulator and finally N (ǫ) = q∈B Θ(ǫ − ǫ 0 (q)) (3.18) denotes the (normalized) number of states (note that we set a = 1 to simplify the algebra). It proves convenient to introduce also the density of states D(ǫ) = q∈B δ(ǫ − ǫ 0 (q)) , (3.19) so that (3.20) The two functions D(ǫ) and N (ǫ) are obviously related to the lattice Green function which, for a SC lattice, reads [21,32,33] G(τ) = 1 π 4 π 0 d q 1 . . . π ImG(τ) , (3.22) with τ = 4 − a 2 /2ǫ. Note that the interval of the spectrum 0 ≤ ǫ k ≤ ǫ max 0 corresponds to the interval −4 ≤ τ ≤ 4 for the auxiliary variable τ. Recently, in reference [21], Loh has obtained a novel integral representation of the Green's function of simple hyper cubic lattices. The resulting one-dimensional integral obtained for G(τ) involves non-oscillating, well behaved functions and it can thus be computed precisely by means of a Gauss quadrature. From the results of reference [21] we obtained where I e (x) = I 0 (x) exp(−x), K e (x) = K 0 (x) exp(x), I 0 (x) and K 0 (x) being the modified Bessel Functions of first and second class respectively.

Various limits
We first note that, for ∞ > k > k max one has the trivial identity N (ǫ k ) = a −4 . Therefore the LMD flow equation (3.16) is identical to the exact NPRG equation (3.13) for the local potential. LMD approximation is thus exact for local theories [20].
Secondly we consider the scaling limit k → 0. We have where v 4 = 1/(32π 2 ) is a geometrical factor, then, the flow equation (3.16) reduces to for arbitrary dimension D and regulator L , can be found in reference [25] while the case of a sharp cut-off was discussed for the first time in the inspiring paper of Hasenfratz-Hazenfratz [38].
?????-8 4D scalar field theory on a lattice Fixed point solutions make sense only for an equation involving but dimensionless functions and variables and emerge in general the limit k → 0. We introduce the dimensionless field x = k −1 φ and potential u k (x) = k −4 U k (φ). The adimensioned flow equation can thus be written = 0 is obviously a special solution. By integration it gives u ⋆ ′ (x) = 0 (Z 2 symmetry) and u ⋆ (x) = v 4 /4, this is the Gaussian fixed point. In order to study the stability of the fixed point we linearize (3.27). Let us define (3.28) and expand equation (3.27) in powers of h, it yields Let us start the analysis with the linearized RG equation    where χ p (x) is a convenient redefinition of Hermite's polynomial H 2p such that its coefficient of degree 2p is one. We have χ 0 (x) = 1, Clearly for p = 0 we have a trivial constant solution. p = 1 corresponds to λ 1 = 2 thus χ 1 (x) is a relevant field. The case p = 2 corresponds to λ 2 = 0 and χ 2 (x) is a marginal field. For all p ≥ 3 the eigenvalue λ p < 0 (for instance λ 3 = −2) correspond to irrelevant solutions χ p (x). The stability of the marginal field χ 2 (x) can be obtained by finding a solution of equation (3.29a) equal to χ 2 at the dominant order. An analysis similar to that of reference [38] reveal that χ 2 is in fact irrelevant beyond the linear approximation. The picture of the scaling fields χ p (x) in D = 4 is thus consistent with a critical point [31]. The usual analysis [31] then yields for the critical exponent ν the classical value ν = 1/λ 1 = 0.5. Since Fisher's exponent η = 0 in the LPA all other (classical) exponents are deduced from scaling relations.
It is generally admitted, and was confirmed by the recent numerical studies of Codello [40] that there is no other fixed point than the Gaussian fixed point in D = 4. We have just shown that the LPA/LMD theory, albeit approximate, supports the existence of this fixed point.

A change of variables
We pointed out in section 3.4 that in the asymptotic limit k → 0 the lattice and off-lattice LPA flow simple pole ω = −1 in the threshold function L (ω) (see equation (3.17)). This point has been studied at length in references [24,25]. Specializing this discussion to the case D = 4 we note that in the limit is a precursor of the spontaneous magnetization φ 0 = lim k→0 φ 0 (k). It follows that the threshold function L diverge in this interval as k −2 . This yields a universal behavior L (φ)/L (φ = 0) = 1 − φ 2 /φ 2 0 . Moreover, as a consequence, U k (φ) becomes convex as k → 0, in particular it becomes constant for −φ 0 < φ < +φ 0 .
The divergence of the threshold function makes impossible to obtain numerical solution of the nonlinear PDE (3.16) in the ordered phase, we really deal with stiff equations. In order to remove stiffness, In contradistinction with equations (3.16) the quasi-linear parabolic PDE (4.1) can easily be integrated out. As in references [10,20,24,25] we made use of the fully implicit predictor-corrector algorithm of Douglas-Jones [27]. This algorithm is unconditionally stable and convergent and introduces an error of O((∆t ) 2 ) + O((∆φ) 2 ) (∆t and ∆φ discrete RG time and field steps respectively) and can be used below and above the critical point as well. In the ordered phase we note that [25]  The initial conditions on the local potential U k at k = Λ are easily transposed to the field L k . It follows from the discussion at the end of section 3.2 that the simplest choice is to choose Λ = k max = 4/a and L Λ = L (a −4 γ ′′ k max (φ)) where γ k max is the local Wetterich function and φ = aφ for all values of the order parameter φ.
Of course, in practice, a cut-off must be imposed on φ and boundary conditions must then be introduced such that the PDE is solved only on the interval −φ max < φ < φ max for all k with some specifications on the boundaries. We made the consistent choice [20,25] L k (±φ max ) = a −4 L (a −4 γ ′′ k (φ max )). Here Wetterich effective function γ k (φ max ) is evaluated in the first approximation of the hopping parameter expansion (see e.g. reference [22]) by assuming the validity of the local approximation. Table 1. Critical parameters of the Φ 4 scalar field theory on a 4D simple cubic lattice in the LPA approximation using the LMD regulator (3.4). From left to right : g , r c (g ). The data were obtained by fixing g and determining r c (g ) by dichotomy. An uncertainty of at most ±1 affects the last digit. g r c (g ) g r c (g ) 0. 10

Solving the flow equations
We solved equation (4.1) with the Douglas-Jones algorithm [27]. We used for most our numerical experiments ∆t = 10 −4 , a maximum of N t = 3 10 5 time steps, ∆φ = 10 −4 and N φ = 30000 field steps (i. e. φ max = 3.). Note that the functions N (ǫ) and D(ǫ) can be computed once for all with the desired precision. In order to determine the critical point r c (g ) one proceeds by dichotomy, g is fixed and one varies r . An illustration of the method is given in figure (2) in the case g = 1000. The renormalized coupling constant u (2) k ≡ U ′′ k (M = 0)/ǫ k , with ǫ k = a 2 k 2 , discriminates the state of the system by its behavior in the limit k → 0.
Of course the Gaussian fixed point, characterized by u (2) k = 0, is never reached but approached only asymptotically for r = r c (g ). As soon as r r c (g ) the flow deviates from the fixed point due to the relevant fields. For r < r c (g ) the coupling constant u (2) k → −1 as t increases; this is the expected behavior in the ordered phase. For r > r c (g ), u (2) k → +∞ when k → 0 (and thus ǫ k → 0) since the compressibility U ′′ k (φ) remains finite for all values of the order parameter φ; the curves escape to +∞ as can be seen on the

figure.
A few dichotomies of r thus yield a very precise estimate of r c (g ). We checked that our values for the parameters ∆t , ∆φ, etc give at least 8 stable figures for r c (g ). We report only 7 figures in the table 1 with the last figure rounded-up. Precision could be enhanced with codes in quadruple precision, unfortunately no such public domain FORTRAN code exists for the calculation of Bessel functions. We explored a wide range of values of parameters with g varying in the range g = 10 −5 (the Gaussian limit) up to g = 100000  Figure 4 displays the inverse compressibility U ′′ k (φ = 0) in the limit k → 0 for g = 0.00001. The fixed point is attained and the expected linear classical behavior of U ′′ k (φ = 0) ⇐ (δr ) is eventually obtained. A linear regression of the right part of the curve gives an exponent of γ −1 = 0.99985 in agreement with the classical value of the compressibility exponent γ = 1. A weak first order transition would yield a discontinuity at some value of r which is never observed for g ≥ 10 −5 . Numerically it proved very difficult to consider smaller values of g smaller than 10 −5 and a code written in quadruple precision should be -5e-10 0 5e-10 δr  necessary to investigate further this question.

Conclusion
In this paper we have computed the critical line of the Φ 4 one-component model on the simple cubic lattice in four dimensions of space in the framework of the NPRG within the LPA approximation. We made use only the smooth LMD regulator which is expected to give the better results. The flow equations have been solved for the threshold functions rather than for the potential. This trick allows to obtain numerical solutions in the ordered phase where the PDE for the potential are stiff and fail to converge. A dichotomy process based on the generically different asymptotic behaviors of the dimensioned inverse susceptibility U ′′ k (φ = 0)/k 2 in zero field, below and above the critical point, yields a very precise determination of the critical line r c (g ). The model is trivial in the sense that all the solutions belong to the basin of attraction of the Gaussian fixed point for all considered values of g . We did not observe a weak first order transition in the Gaussian limit g → 0, at least, numerically, for g > 10 −5 . A numerical exploration of still lower values of parameter g would require a quadruple precision code which is out of reach for the moment.
In reference [20] we obtained an excellent agreement between our estimates of the critical line of the 3D Φ 4 model on a simple three dimensional lattice and that of Monte Carlo simulations of Hasenbush [28]. In D = 3 the LPA approximation does not yield the exact critical exponents contrary to the case D = 4 where the classical exponents are found. One can thus a fortiori expect an excellent agreement for the critical line between the theory and the simulations in 4D. Unfortunately, we were unable to find estimates of the critical line of the 4D version of the model by means of Monte Carlo simulations in the literature.