Revisited functional renormalization group approach for random matrices in the large-$N$ limit

The nonperturbative renormalization group has been considered as a solid framework to investigate fixed point and critical exponents for matrix and tensor models, expected to correspond with the so-called double scaling limit. In this paper, we focus on matrix models and address the question of the compatibility between the approximations used to solve the exact renormalization group equation and the modified Ward identities coming from the regulator. We show in particular that standard local potential approximation strongly violates the Ward identities, especially in the vicinity of the interacting fixed point. Extending the theory space including derivative couplings, we recover an interacting fixed point with a critical exponent not so far from the exact result, but with a nonzero value for derivative couplings, evoking a strong dependence concerning the regulator. Finally, we consider a modified regulator, allowing to keep the flow not so far from the ultralocal region and recover the results of the literature up to a slight improvement.


I. INTRODUCTION
Random matrix models are specific statistical models describing (Euclidean) quantum fluctuations of a matrixlike field [1]. They appear as a framework for a very large kind of problems in physics and mathematics, from quantum gravity to biology (the list of references is very large, and we do not mention them here). In this paper, we essentially focus on quantum gravity interpretation, even if our results do not especially refer to this interpretation. The link between matrix models and twodimensional quantum gravity arises from the observation that the perturbation series of random matrix models can generate randomly arbitrary triangulated surfaces (see [2]- [11] and references therein); the precise relation between Feynman diagrams and elementary polygons being discussed on a concrete example in section II. Then, Feynman amplitudes of the perturbative expansion for such models are indexed by simplicial decomposition of the two-dimensional manifold; and as an important result (in particular for quantum gravity issues), the relative weight of two such a triangulation depends only on the genus of the corresponding manifold and the size N of the considered matrices [1].
In the large N limit, only planar surfaces survive, and the computation of the corresponding free energy shows the existence of a critical point, where infinitely refined * Electronic address: vincent.lahoche@cea.fr † Electronic address: dine.ousmanesamary@cipma.uac.bj triangulation starts to dominate the perturbative series; and interpreted as a continuum limit. Double scaling is a theoretical framework allowing keeping into account higher genus surfaces, taking the large N limit near the critical point in such a way that the relative weight of the different topological configurations are exactly compensated by their growth, fixed by a universal critical exponent. Renormalization group (RG) has been considered to be an alternative to the standard analytic method [12]- [19]. The argument is based on the interpretation of the correlation between coupling and N in the double scaling limit as a fixed point of the RG flow with N . In a Wilsonian perspective, integrating out matrix entries between N and N −δN generates effective actions, which drag the couplings so far from their initial values. This version of the RG flow and perturbative solutions has been investigated for twenty years [2]- [40], and reproduces semiquantitatively the exact results. More recently, a nonperturbative FRG framework has been considered to improve the perturbative results [12]. In this reference paper, the authors show convergence phenomena for the computed critical exponents toward the exact (i.e. analytic result) for double scaling.
In the following paper, we show that the naive approaches to solve the nonperturbative RG equations, especially based on a reduction of the theory space to the strictly local interactions or products of them are strongly incompatible with Ward's identities [41]- [49]. The origin of the incompatibility is traced to come from the regulator itself. Indeed, due to the presence of the regulator, the compatibility with Ward identities requires to enlarge the theory space to derivative couplings; which in turn seems to play a significant role in the fixed point structure, and finally introduce a spurious dependence on the regulator. To solve this issue, we introduce a modified regulator, parametrized in such a way that the contribution of derivative couplings in the Ward identities remains small in a significant domain of the RG flow, so that truncation involving only traces may be used without strong disagreements to approximate the exact solution of the RG equations. Note that we explicitly checked that the method used to derive the flow equations in the reference [12] (i.e. a systematic projection using a U (N )-invariant vacuum ansatz) is inconsistent, and the first part of this paper, voluntary pedagogical provide another derivation of the flow equations for truncation involving trace observables.
The outline is the following. Sections II reviews shortly on the matrix models, double scaling and functional renormalization group approach, including all the preliminaries required for the rest of the discussion. In section III we review the nonperturbative renormalization group flow following the reference [12], first in the local potential approximation and second including multitrace interactions. In section IV we show explicitly that the local potential approximation strongly violates Ward identities for some choice of natural regularization functions. We investigate the flow numerically and compare the numerical fixed point with the analytic solutions. In the last section (V) we provide some discussions and the conclusion of this work.

A. A short review on matrix models
To shortly reviewing matrix models, let us consider a concrete example for a trivalent model involving Hermitian N × N matrix φ, described from the partition function: where dM is the invariant measure on the N × N Hermitian matrices (for more detail see [2]- [4]). The classical action S[φ] := 1 2 Trφ 2 + g √ N Trφ 3 admits a natural U (N ) symmetry due to the global trace structure. Expanding the right hand side perturbatively in λ with the propagator we generate Feynman amplitudes labeled by connected ribbon graphs G, that is to say, a set of vertices, edges and faces. The interaction vertex has three external points, identifying the six strands pairwise. Propagator, vertex and their dual correspondence are pictured on Figure 1a, and an example of ribbon graph is given in Figure 1b. From Feynman rules, the partition function expands as a sum of Feynman amplitudes where, up to the rescaling φ → √ N φ, the amplitude A G depends on N and on the genius h of the dual representation ∆ G of G : A G = N 2−2h(∆ G ) . We stressed that matrix models are statistical models for triangulated surfaces, but we have not made contact yet with quantum gravity in dimension two. This correspondence can be heuristically traced as follows. Including cosmological constant Λ, classical gravity in dimension two is described by the action: where we used the Gauss-Bonnet theorem to compute the integral in terms of the Euler characteristic χ(M), and where we have denoted by A M the area of the surface M. Then, the theory only depends on two parameters, and we generally assume that only these two parameters are relevant to define the discretization. As a basic example, introducing an equilateral triangulation ∆ M of M, such that each triangle has a fixed area a, the action (4) can be discretized as and the quantum theory described by the partition function: matches the partition function (3), up to the identification: As a result, heuristically, the large N limit of matrix models (involving a lot of "microscopic" degrees of freedom) matches with the weak coupling regime of two dimensional topological gravity. Some formal results show the equivalence between matrix models and other quantum gravity approach. In particular, equivalence with Liouville theory at fixed topology has been stressed from agreement with KPZ relation, see [38]- [39].
The leading order contribution in the N → ∞ limit comes from the triangulations with zero genus, corresponding to a planar topology. Interestingly, a closed two-dimensional topological manifold is fully characterized by its genus and orientability. Note that Hermitian matrices only generate orientable surfaces, so that the genus fully determines the topology of the triangulation and allows to capture nonperturbative effects. Indeed, the perturbative expansion (3) can be rewritten as a topological expansion : where we have defined the sum over all triangulations with genus h as Z h (g). In the large N limit, the partition function for zero genus surfaces, Z 0 (g) has the following critical behavior [2]- [4]: where γ = −1/2. As a result, the free energy of the planar sector diverges around the critical point g = g c corresponding to the continuum limit, where Z 0 is dominated by arbitrary large triangulations. Going beyond the planar sector requires a double scaling limit, taking the two limits N → ∞ and g → g c in a correlated manner [1]. More precisely, we can show that Z h has the same critical point as Z 0 for any h: suggesting to take simultaneously the limits N → ∞ and g → g c in such a way that the ratio remains fixed such that all the topologies contribute to the free energy when we are close to the critical point : corresponding to the continuum limit, where the area diverge, like for the critical behavior of the naive N → ∞ limit. The double scaling limit may be analytically investigated, and exact results for g c and critical exponents have found, see [1] for details.

B. Flowing on the matrix theory space
In order to make contact with the reference papers [12]- [18] in view to compare our results with the ones, and in contrast with the model considered in the previous section, we focus on a quartic model, describing a random Hermitian matrix with the classical action: In addition to the U (N ) symmetry, this model as a discrete Z 2 symmetry φ → −φ, and generate squarulations rather than triangulations. This distinction is unimportant for the continuum limit that we will investigate, which does not depend on the choice of elementary discrete polygons used to build random surfaces. For this model, the critical value g c and the corresponding critical exponent θ in the continuum limit have been exactly computed [1]: The elementary intuition allowing to consider renormalization group approach to investigate the continuum limit for matrix comes from the constraint (11), freely interpreted as a fixed point for an appropriate scaling in N , at which any change as N → N − δN may be compensated by a change g → g+δg of the coupling, without change of the continuum physics [2]. This elementary observation suggests, in accordance with a Wilsonian point of view, to integrate out step by step the large N degrees of freedom on lines and rows, of the N × N matrices, reducing them to (N − 1) × (N − 1) matrices after a single step, (N − 2) × (N − 2) matrices after two steps, and so on. To each step, (N − i) × (N − i) matrices are described by effective action, which is a sum of two pieces: The classical action for (N − i) × (N − i) matrices, and the fluctuations term arising from integration of N −i+1 degrees of freedom. As a result, to each step, the couplings have the discrete change rule: where the notation suggests that we consider only the large N limit to define the β-function. Computing β(g i ) from a single step, we get the one-loop beta function [2]: which vanishes for g * = −1/6, in qualitative accordance with the exact result (14); the string susceptibility γ being related to the critical exponent −β (g * ) as [2]: .
The accuracy may be explicitly improved taking into account higher couplings and loop effects, an observation which strongly motivates a nonperturbative analysis, as suggested in [12]. In this reference, the authors introduced a FRG framework based on the Wetterich equation formalism. The new version of the Wilson RG procedure requires a splitting into modes, between UV scales (no fluctuations are integrated out) and IR scales (all the fluctuations are integrated out) dictating how the small distance fluctuations are integrated out. As we will see, different steps correspond to partial integration of modes between N and N −δN , and following the standard strategy in FRG formalism, we introduce a new term in the classical action: which behaves like a scale dependence mass term, the specific slicing in N depending on the shape of the regulator r N (a, b). Among the standard properties of the regulator, we recall the following (for more explanations, the reader may be consult the standard reviews [31]: 1. r N (a, b) has to have a nonvanishing "infrared" 3. r N (a, b) has to vanish in the limit N → 0, allowing to recover the original partition function. .
where the dot product is defined as A · B := mn A mn B mn . Due to the scale dependence of the regulator, the long distance physics effects ((m, n) N ) acquire a large mass and are frozen out, whereas the small distance effects ((m, n) > N ) are integrated out. The RG flow then relates Z N to Z N −δN . The transcription of this relation goes through a first order differential equation: which indicates how the average effective action Γ N change in the windows of scale [N, N − dN ] -the dot meaning derivative with respect to the RG parameter t := ln N :Ẋ = N d dN X. We recall that the average effective action is defined as slightly modified Legendre transform of the free energy W N := ln Z N : where Φ denotes the classical field: In the same way Γ (2) N in equation (20) denotes the second derivative of the average effective action : Even to close this section we have to add an important comment about the notion of canonical scaling. Scaling, that is to say the dependence of the quantities on the cutoff coming from their dimensions, plays generally an important role in renormalization. In standard quantum field theory for instance, dimensionality is closely related to renormalizability. For matrix models, the situation is quite different, because there are no referent space-time, no referent length and no canonical scaling coming from extra structure of the theory. However, the behavior of the RG flow with N in the vicinity of the Gaussian fixed point (i.e. keeping only the part of the scaling which is independent of the couplings), provides an intrinsic notion of dimension, that we call canonical dimension: Definition 1. For any trace observable g k Tr(φ k ) in the classical action, the canonical dimension of the coupling constant g k is defined in the vicinity of the Gaussian fixed point as the part of the scaling in N which is independent of g k and the other couplings.
We denote as d k the canonical dimension of g k , so that the intrinsic scaling writes as N d k +O(g1,g2, ··· ) . To find the explicit expression of d k , we then have to be investigate the behavior of the Feynman diagrams with N . This may be traced from the link between two-dimensional quantum gravity recalled in the previous section. Up to the rescaling φ → √ N φ we have stressed a relation between matrix coupling, N , Newton and cosmological constant. Keeping this relation implies that each Feynman diagrams scales exactly as N 2−2h ≡ N χ , where χ(∆) := V (∆) − E(∆) + F (∆) denote the Euler characteristic of the polygon decomposition ∆, having V vertices, E edges and F faces. It is not hard to see that this holds if and only if, up to the mentioned rescaling, the only N dependence of the classical action comes from a global N factor, enforcing the definition in agreement with formula (1). In this paper, we will consider also multitrace interactions at the level of the effective action, and we have to extend this formula for such interactions. In order to remain in accordance with the expected scaling N χ , we impose to cancel the additional N factors coming from the additional traces. As a result, for an observable of the form n j=1 Tr(φ k(j) ), one assigns the canonical dimension d For a double trace operator for instance Tr(φ k )Tr(φ l ) one gets d (2) kl = −(k + l)/2. As pointed out in [12], it is interesting to note that, even for a single trace operator, the canonical dimension is negative for k > 2, meaning that all non-Gaussian couplings are irrelevant. In this situation, the improvement of the scaling coming from radiative corrections plays an essential role in the fixed point structure.

C. Ward-Takahashi identity
Ward-Takahashi identities are a general feature of symmetry in quantum field theory and may be viewed as a quantum version of the Noether's theorem, resulting in the translation invariance of the Lebesgue integration measure in the path integral definition of the partition function (the reader could consult [48]- [49] for the first derivation of these identities in QED). Their interest in RG investigations has been largely discussed in the literature, and more specifically in the context of tensorial field theories in [41]- [49]. Our point of view is that Ward identities are nontrivial functional relations, depending on the regulator like flow equations, and with this respect have to take into account in the building of the RG approximations. This is what we will do in the next section. We will extend this discussion about the role of Ward identity in section III. In complement, the reader may consult [44]- [45].
Without a regulator term, only the source terms break the global U (Λ) invariance for some cutoff Λ. Infinitesimal variations provide the identity: to all orders of the perturbation. Note that for the rest of this paper we restrict our investigations into the symmetric phase where vanishing classical field Φ defined from equation (22) is expected to be a good vacuum and all the odd correlation functions vanish identically.
The regulator term 1 2 φ r k φ breaks explicitly the global U (Λ) invariance, and adds a new contribution to the asymptotic Ward identity (26). Let us consider an infinitesimal unitary transformation 1 + , being an infinitesimal anti-Hermitian operator. At the leading order, the transformation rule for the matrix field φ is: At the leading order in , the total variation of the generating functional Z N writes as Because S is a sum of traces, δS[φ] = 0. The variation of the source term is noting but: The variation of the regulation term can be deduced in the same way: where we assumed that [r N (a, b)] ab;cd = [r N (a, b)] cd;ab . This is exactly the same computation as for the source term, up to the replacement J ab → c,d [r N (a, b)] ab;cd φ cd , leading to: Moreover we assumed that r N (a, b) is a symmetric function with respect to a and b. Due to the translation invariance of the Lebesgue measure, Z N must be invariant up to a global reparametrization of the fields, therefore the variation of the left-hand side in (28) must be vanish δZ N = 0. From the identity: We finally deduce the following statement: In the symmetric phase, and along the path N = Λ to N = 0, the following relation holds: where we adopted the Einstein convention for repeated indices. Note that there are no summation over indices a and c.

III. SOLVING THE RG FLOW IN THE SYMMETRIC PHASE
The exact flow equation (20) cannot be solved exactly except for very special problems. Extracting information about the nonperturbative behavior of the RG flow then requires an appropriate scheme of approximation. In this section, we review a standard approach based on a crude truncation of the theory space. We focus on local interactions, or products of them, to remain closer to what we expect to be the theory space of the original matrix model, without regulator. As mentioned in the Introduction, this section is voluntarily pedagogical, due to strong disagreements with some results in the principal cited reference [12].
A. Local potential approximation i.) Local potential. The matrix action is nonlocal in the usual sense in field theory because all the interacting fields are not evaluated on the same point of the structure manifold. However, what allows to say that two objects interact locally is precisely the form of the interaction. The interaction then allows to define by themselves an appropriate locality principle, and we adopt the definition: Definition 2. Any global trace of the form Tr(φ k ) is said to be a local monomial interaction. In the same way, any functional of U [φ] which may be expanded as a sum of single traces is said to be a local functional.
Note that this locality principle reflects the proper invariance of the interactions concerning unitary transformations 1 .
The first parametrization of the theory space that we consider split the effective action Γ N (Φ) as a sum of two kind of terms: The last term U N (Φ) designates the purely local potential, expanding as a sum of single trace observables: (35) Following [12], we introduced a field strength renormalization Z N in front of the Gaussian term. The renormalized quantities are generally defined from a fixed coefficient in the Gaussian part of the original action. Rescaling the fields such that the mass term reduces to its free term 1 2 TrΦ 2 , we define the dimensionless and renormalized couplings u k,N as: As pointed out in the derivation of the Ward identity, the presence of the regulator breaks the U (Λ) invariance of the original action, and the RG flow has to generate noninvariante momentum dependent effective interactions such that, for instance: where the Taylor expansion of the function q starts at the order 1 in a/N and b/N . The terminology "momentum dependent" simply reflects the situation in ordinary quantum field theory, the indices of the matrix field playing the role of discrete momenta. Expanding q in power of a/N and b/N corresponds to the standard derivative expansion. As we will see from Ward identity, such a deviation from strict locality introduces relevant corrections at the leading order in 1/N , and must be kept in the large N limit. In particular, in the closure procedure around quartic interactions, the linear coupling: plays an important role in the fixed point structure, improving strongly the local potential approximation. In a first time, in order to compare them, we keep only the strong local part of the decomposition (34): The flow equations for couplings g n in the parametrization (39) can be deduced from the exact Wetterich equation deriving n time with respect to Φ and setting Φ = 0 (we recall that we work in the symmetric phase). Because Φ is a Hermitian matrix, Φ ab = Φ * ba , and : from which we get: where g ab,cd := δ ac δ bd is nothing but the "bare" propagators. For the regulator function, we chose a modified version of the Litim optimized regulator [25], allowing to make analytic computations: which obviously satisfy the requirements 1−4 given after equation (18). Taking the derivative with respect to the flow parameter t = ln N , we get straightforwardly: where we introduced the anomalous dimension Taking successive derivative with respect to Φ of the exact flow equation (20), and setting Φ = 0, we deduce the flow equations for all couplings involved in (39). For each step, all contributions involve some powers of the effective propagator G N = (Γ (2) N + r N ) −1 , evaluated for vanishing Φ, and for a + b ≤ 2N as: The one-loop sums that we will encounter in the derivation of the flow equations are all of the form: In the large N limit, the sum can be replaced by an integration up to 1/N corrections. Let us introduce the continuous variable 2N x := a + b, running from a/2N to 1: leading to: ii.) Truncated RG flow. Deriving the equation (20) twice with respect to the Φ fields, and setting Φ = 0, we get: N,ef,lm,ab,baG N,lm,cd , withG N,ab,cd := (G NṙN ) ab,cd and where once again we sum over repeated indices. To compute the sums, we have to take into account the symmetry structure of the external indices. From (45), we get for instance where we defined: For a fixed configuration of the external indices, there are two leading order contractions, both pictured on Figure 2, where in this graphical representation the dotted edges correspond to the contraction with the effective propagator G N,cd,efGN,lm,cd given by equation (50).
We now have to compute how many leading order contractions such as the one pictured on Figure 2 contribute. It is not hard to see that there are exactly 4 × 2 different Finally, computing the derivative of the left-hand side of equation for zero external momenta, we get : from which we deduce that: Divided by Z N , and from definitions (36), we then get finally; The computation of the beta function β 4 :=u 4,N follows the same strategy. Deriving once again twice with respect to the Φ fields, and setting Φ = 0 at the end of the computation, one gets, formally: The relevant diagrams corresponding to the two kinds of traces involved in these expressions, all including one internal face are pictured on Figure 3. Each of them may be easily translated into an equation like for the 2-point diagrams. For zero external momenta we get: and: Once again the numerical factors may be easily understood. For instance, for the diagram involving a 6-point vertex, there are six different ways to choose the first end point of the contracted edge, two to choose the second one, to make a leading order graph; and finally 4! ways to exchange the remaining external points. Because from definition Γ N,00,00,00,00 = 6g 4,N , it follows that: leading to: Remark 1. Neglecting the coupling u 6,N and expending the remaining right-hand side in power of u 4,N , up to order u 3 4,N , we do not reproduce the one-loop result (16). In particular, the numerical factor in front of u 4,N becomes 14/3. This cannot be viewed as a defect of the approach, the 1-loop beta function being nonuniversal for coupling with nonzero canonical dimension, as it can be easily checked.
Following the same procedure we can compute beta function for higher couplings, the flow equation for g k,N involving g k+2,N and so on, providing an infinite tower of hierarchical equations. The truncation method is the simpler approximation procedure, which truncates crudely in the full theory space, setting g k,N ≈ 0 for some k. This method has the advantage to be very tractable for (strict) nonlocal interactions, which is the case for matrix models. For k = 8, i.e. setting g 8,N ≈ 0 we find for the coupling g 6,N : where we neglected the term − 1 2 TrGΓ Computing each trace like for the two previous cases, we get two relevant diagrams at leading order: (61) Taking into account the permutation symmetries the two diagrams are respectively evaluated to 4! × 4g 6 g 4 I 0 ; and we obtain for β 6 : To summarize, we have the following statement: Proposition 1. In the large N limit, and in the local potential approximation, the truncated RG flow around φ 6 interactions is described by the following closed system: with: Note that the truncated RG flow becomes singular for u 4,N = −3, splitting the reduced phase space into disconnected regions. We call the perturbative region the region connected to the Gaussian fixed point. Moreover, we can remark that this singularity holds for arbitrary higher truncations. These equations allow to investigate the existence of nontrivial interacting fixed points for quartic and sextic truncations.
• For k = 6, the fixed point equation reduces to: The numerical plot of the beta function is given in Figure  4. We get two solutions: The first solution is under the singularity line u 4,N = −3, and therefore unconnected to the Gaussian fixed point. The second solution however is in the perturbative region, and corresponds to an UV-attractive fixed point.
As expected, the result seems to be improved when the order of the truncation is increased. The fixed point that we found is reminiscent to the standard Wilson-Fisher fixed point, with only one attractive and one repulsive direction in the UV (i.e. in the large N limit); the single positive critical exponents having to play the role of β (g * ) in equation (17).
The reliability of these results may be traced by investigating higher truncations. For k = 10, we have to add the contribution − 1 2 TrGΓ which becomes: Foru 8 , taking into account only the leading order contractions, we get: Solving numerically the flow equations, we get three fixed points, the first one being with critical exponents: and anomalous dimension η * ≈ 0.11.
Once again we find some results in qualitative accordance with the exact computation. We recover a Wilson-Fisher like fixed point having the expected characteristics, only one relevant direction with positive critical exponent. However, we do not observe significant improvement concerning k = 8 truncation. This seems to indicate that higher irrelevant operators do not contribute much to the accuracy of the critical exponents. Note that this result is in complete disagreement with the ones of [12], where a convergence phenomenon has been pointed out by the authors. We suspect that this disagreement is a consequence of the method used by the authors, which, setting a diagonal vacuum Φ ab = aδ ab to extract the flow equations, and therefore have selected more than strictly local interactions.
Interestingly, the numerical critical value for the coupling seems to be so far from the exact values than the ones obtained from k = 6 and k = 8 truncations. This value is not universal so that a disagreement with the exact value cannot be relevant for the reliability of the analysis. One expects that this is a defect of the LPA. Indeed, the fixed point arises essentially from the flow of irrelevant operators, which may be strongly coupled at the fixed point, where irrelevant interactions for the Gaussian counting can contribute significantly. Then, when we take into account higher interactions in LPA, we lost more and more information, coming especially from nonlocal and multitrace operators, as pointed out in [12]. To investigate the improvement coming from these operators, let us consider the k = 8 truncation, involving double and triple traces (to simplify the notation, we left the N index for couplings): The truncation for local interaction was based on the canonical dimension. For k = 8 for instance, one can say that we discarded interactions with a canonical dimension smaller than d = −3. If we think to build the same truncation including multitrace interactions, we could conclude that interactions such that Tr(Φ 2 )Tr(Φ 4 ), which have canonical dimension d = −3 must be discarded like Tr(Φ 8 ) interactions. However, the double trace increases the strength of the coupling. Then, in contrast to Tr(Φ 8 ), an interaction such that Tr(Φ 2 )Tr(Φ 4 ) contributes directly to the flow of g 4 at leading order, the tadpole contraction scaling as N 2 .
Starting with the computation ofŻ, we show that the contribution (54) holds, but has to be completed with double-trace diagrams. Then,at the leading order in N , we get, graphically: Computing the new diagram, and taking into account that we have 4 × 2 different permutations leading to the same diagram, we obtain: where we defined J (p) as: which can be approached by an integral as Finally, defining the dimensionless and renormalized couplings v i,j,k,··· as: we obtain, in replacement of the equation (54): leading to: In the same way, for β 4 , the previous computation has to be completed with the diagram: leading to: Finally, at this order for the truncation, the expression for β 6 is unaffected the multitrace interactions, except through the improvement of η. Now, we move on to the computation of the remaining beta functions, β 2,2 , β 4,2 and β 2,2,2 , respectively for couplings v 2,2 , v 4,2 and v 2,2,2 . For β 2,2 , we get two kinds of leading order contractions, nonvanishing ones being: For the first kind of diagram, we get the contribution: The second kind of diagrams comes from the term 3TrG Γ Finally, the third and last contribution involving a nontrivial loop arises from the term 1 2 TrG Γ Note that contributions involving loop without sum, such that the external indices fix the momentum along the loop vanish identically for zero external momenta. Indeed: Therefore, we obtain for β 22 : In the same way, for β 4,2 and β 2,2,2 , the nonvanishing typical diagrams are the following: leading to the following statement: Proposition 2. In the large N limit, the multitrace truncated RG flow around sixtic interactions interaction is described by the following set of equations: with: As for the single trace potential, we may investigate numerically the fixed point structure, increasing progressively the degree of the truncation.
• For k = 6, the set of equations that we have to solve is the following : Numerically, we get only one interesting fixed point, for v 22 = 0 and u 4 ≈ −0.20. This fixed point corresponds to the one discovered from the single trace k = 6 truncation, therefore, we expect that there are no improvements coming from multitrace at this order of approximation.
• For k = 8, we may distinguish two cases. In the first one we only consider the effect of double-trace interactions, neglecting the triple trace. In the second one, we include the triple-trace interaction. For the double-trace approximation, we get only one potential candidate fixed point, for the values: with anomalous dimension η * ≈ 0.14 and critical exponents: (θ 1 , θ 2 , θ 3 , θ 4 ) ≈ (1.04, −0.93, −1.03, −1.78) .
Our expectation about the role of the disconnected interactions seems to be disappointed. No significant improvement is observed, for the critical exponent the value is essentially the same as for LPA k = 8 and k = 10 truncations, and the value of the coupling u 4 * is the same as for the LPA. The only change is for the second critical exponent, whose value is slightly diminished from their purely local version. This, once again, is in complete disagreement with the results of [12], and seems to indicate a rapid convergence toward θ ≈ 1.04. This intuition is confirmed taking into account triple-trace interactions. We get once again the same fixed point, with v 222 = 0 and the same critical exponent θ ≈ 1.04. This seems to indicate that triple-trace interactions start to be relevant for k = 10 truncation, exactly as the double start to be relevant for k = 8.

IV. COMPATIBILITY WITH WARD IDENTITIES
As we explained in the previous section, there is no preferred notion of scale for the initial model (i.e. for the model without regulator). More precisely, there are canonical notions of deep UV and deep IR: the deep UV being related to the classical action S(φ), without integration over statistical fluctuations and the opposite deep IR scale, related to the effective action Γ[Φ], when all fluctuations are integrated out. However, there is no canonical way to reach the deep IR region from the deep UV one. All the fluctuations play the same roles, and are indistinguishable "UV" or "IR". All the ways that we think to cut through scales are a priori allowed and this difficulty is related to the triviality of the Gaussian term. For standard field theories, this is the spectrum of the kinetic operator which provides a canonical path from UV to IR, allowing to classify the fluctuations following their respective energy. But for matrix models, due to the U (Λ) invariance, all the eigenvalues of the kinetic operator are the same, and the fluctuations become indistinguishable.
This highlights the role of the regulator. The regulator that we introduced broke the global U (Λ)-invariance at the kinetic level, providing a preferred path from UV to IR and an ordering for partial integrations over quantum fluctuations.
The Ward identity that we derived in section II B is a consequence of this symmetry breaking. It arises from the nontrivial variation of the kinetic term under infinitesimal unitary transformation, like the flow equations arise from a nontrivial variation of the kinetic term under a change of the running scale N . Both are consequences of the symmetry breaking and have to be treated on the same footing, as nontrivial relations between effective vertices Γ (n+2) and Γ (n) . More precisely, and as it will be clearer in the rest of this section, one can say that the RG equation dictates how to move through increasing scales (from large to small N ) whereas the Ward identity dictates how to move in the momentum space. As we will see, because of the symmetry breaking, nonlocal derivative-like interactions such that (38) appear even in the strictly local sector, and play an important role in the behavior of the RG flow, especially around the UV fixed point.

A. Explicit Ward identities and enlarged theory space
i.) Explicit Ward identities. Taking the derivative with respect to ∂ 2 /∂J de ∂J d e of the Ward identity (33) and setting J = 0, we get: Setting d = a and d = e = e = c, and a = c, we get: where we used of the fact that G for any permutation π of four elements. In order to make contact with the parametrization (39), we decompose G where for convenience we introduced the reduced 2-point components g (2) ab , defined from the 2-point function G (2) ab,cd as: ab,cd =: g (2) ab δ ad δ bc .
(97) Now, return on equation (90), and consider the pair (ac). If a = c, we see that there are four different ways to choose the position of the pair, but for each choice, there remains only one way to fix the relative position of the three other pairs 3 . Therefore (90) may be rewritten as: ac g (2) cc γ (4) ac,cc,bc,ba g (2) bc g Dividing by g (2) ac g (2) cc , we deduce the following statement: Lemma 1. In the large N limit, the 4 and 2-point functions must satisfy the nontrivial relation : ac,cc,bc,ba g (2) bc g As explained in the section II B, the Ward identity dictate how to move into the momentum space from ultralocality, whereas the flow equations (20) dictates how to move through scales, from UV to IR. This is in this way that flow equations and Ward identities, both consequences of the U (Λ) symmetry breaking, cannot be considered separately.
Setting a = c + 1, and for sufficiently large N may be estimated from the same continuous approximation used to compute the sums in the previous section, that is to say: (99) Note that this approximation has to be used carefully, and for formal derivations, we may use derivative first as a notation. Computing the derivative for the Litim regulator, we get: In the same way, assuming that g (2) ab may be continued as an analytic function g (2) (x, y) for the continuous variables x, y := a/N, b/N , we get: (100), it is clear that the windows of momenta allowed in the sum over b from df /dx := f is the same as the one allowed byṙ N in the flow equation (20). Therefore, the same approximations used to solve the RG equations may be used for g (2) bc , g (2) ba and γ (4) ac,cc,bc,ba . The same situation has been observed for tensor field theory (see [44]), for several choices of regulator function. Then, one expects that this is not a well consequence of the Litim regulator, but a general feature that the allowed windows of momenta forṙ N cover the one of f . Moreover, equation (101) points out the existence of a strong relation between 4-point functions and the momenta variations of the 2-point functions along the path from the deep UV sector to the IR sector. Therefore, and as we will see explicitly, even in the large N limit, nonlocal interactions such that (38) survive at the leading order in 1/N and cannot be discarded from any relevant parametrization of the phase space. This argument shows that strictly local potential approximation have to be enlarged with derivative-like interaction to become compatible with Ward identity. As a first improvement, we can consider the following minimal enlargement : where U N [Φ] expands as a single trace like in equation (35). We call improved LPA this parametrization allowing a small deviation from the crude LPA. From this approximation, and γ (4) ac,cc,bc,ba → Inserting these relations into the lemma 1, we see that the second term on the left-hand side is exactly compensated with the same term on the right-hand side. Then, setting c = 0, we get the following statement: Proposition 3. Up to 1/N corrections, and in the improved local potential approximation, the 2-point derivative coupling γ and the local 4-point coupling u 4 satisfy: where Z Nγ =: γ and L N = (Z N ) −2 NL N bc ) 2 f (c, b) c=0 .
With the truncation (102),L N depends only onγ. Therefore, deriving equation (105) with respect to t leads to This equation that we call Ward constraint relies on two beta functions along the history of the RG flow, since N remain large. As an important consequence: Corollary 1. In the large N limit, any fixed point of the flow equations satisfies the Ward constraint (105).
The flow of the nonlocal kinetic coupling γ receives two kinds of contributions. A first contribution arises from the derivative with respect to one external momentum of the loop integrals, but a direct computation shows that these variations vanish identically. A second contribution arises from the derivative of the effective vertex themselves. In the local potential approximation, the vertex does not depend on the external momenta. But from the Ward identity, it follows that the ultralocal information determines completely the first derivative with respect to the external momenta, like ultralocal 4-point coupling u 4 determines γ in lemma 1. In order to obtain the first derivative of the 4-point function, we need to the Ward identity involving 6-point functions (i.e. derived from (33) deriving four time with respect to the source J.).
It is more convenient to write the original Ward identity (33) as: Taking the derivative four time with respect to the classical field ∂ 4 /∂Φ de ∂Φ d e ∂Φ pq ∂Φ p q of the equation (108), we get: ba g (2) cb Γ (6) ab,bc,de,d e ,pq,p q − 6g ap ,de,d e ,pq δ q c − Γ ap,de,d e ,p q δ qc − Γ ad ,de,pq,p q δ e c − Γ (4) ad,d e ,pq,p q δ ce + Γ (4) q c,de,d e ,pq δ ap + Γ (4) qc,de,d e ,p q δ ap + Γ (4) e c,de,pq,p q δ ad + Γ (4) ec,d e ,pq,p q δ ad , where the 6 on the second term in the left-hand side is a short notation for the 3 × 2 terms corresponding to the different pairing of the derived variables. Setting d = a and d = e = e = p = q = q = p = c for c = a; and keeping only the leading order contractions in the large N limit, the previous relation reduces to the following lemma: Lemma 2. At the leading order in the 1/N expansion, the 6, 4 and 2-point vertex functions must satisfy the following relation: ba g (2) cb Γ (6) ab,bc,cc,ca,cc,cc − 6g (2) ba g (2) cb g (2) c b Γ (4) c b ,bc,cc,cc Γ (4) b c ,ab,cc,ca = − 3Γ (4) ca,ac,cc,cc − Γ (4) cc,cc,cc,cc .
Setting a = c + 1, and keeping only the first term in the 1/N expansion of the difference f (a/N, b/N ) − f (c/N, b/N ), the argument used for the previous explicit Ward identity holds : the windows of momenta allowed by the distribution f (c/N, b/N ) are the same as forṙ N involved in the flow equation (20); and to make sense, the same approximations used to solve this one have to be used in the computation of the Ward identities. Like for the 4-point vertices we introduce γ (6) , with fixed boundaries, such that: the sum over π running through the permutation of six elements, and on the left-hand side, into the sum over b, we replace f c,c,c,b by g 4 /4 and γ (6) cc,cc,cc,cc,bc,cb → g 6 6 . (110) On the right-hand side, from definition (97), Γ c,c,c,c = 4!f c,c,c,c . For Γ (4) ac,ca,cc,cc however, there are only 4 × 2 different configurations for the external indices providing a nonzero contribution. As a result: 3Γ (4) ca,ac,cc,cc − Γ (4) cc,cc,cc,cc = 4!(f a,c,c,c − f c,c,c,c ) . (111) As for the 2-point function, we assume that in the large N limit f c,c,a,c behaves like a continuous functioñ f (x, c/N ) for the continuous variable x = a/N , such that f (a/N, c/N ) ≡ f c,c,a,c . We then define, at leading order in 1/N : and at the first order in 1/N , the lemma 2 becomes, setting c = 0 and simplifying the global factor 1/N : where we defined: From the definition u 6 = (Z N ) −3 N 2 g 6 , we finally deduce the following statement, between renormalized quantities: Proposition 4. Up to 1/N corrections, and in the improved local potential approximation, the 4-point derivative coupling and the local 4 and 6-point renormalized couplings u 4 and u 6 are related as: where we defined: Equations (105) and (115) show explicitly that the strictly local flow strongly violates the Ward identity. This is especially true at the fixed point, where, from equation (105) we see thatγ and u 4 have to be of the same order, indicating that the regulator scheme strongly influences the nature of the theory space. As we will see in the next section, a systematic analysis, including the flow of the derivative couplings seems to confirm this pessimistic forecast, despite the accordance of the resulting critical exponent with the expected value.

B. Strong deviation with local fixed point
Now, we move onto derivation of the flow equations in the parametrization (102). A first change concerns the effective propagator (45), which becomes: such that the integral I (p) a , equation (47) becomes: To simplify the discussion, we introduce the sequence ι p,q (y) for the continuous variable y = a/2N such that: and: (120) In addition we defined the renormalized loopĪ a . Equation (54) is then transformed as: solved as: where we used the concise notation ι p,q ≡ ι p,q (0). In the same way, we get for β 4 , in replacement of (60): The flow equation for γ can be deduced from (49), like η N . From definition: we get (β γ ≡γ): It is easy to check that the involved derivatives ι p,q ≡ ι p,q (0) vanish identically such that the equation forγ reduces to: where we took into account that what we calledL N and U N may be expressed in terms of the sequences ι p,q , Note that the origin of the factor 12 ≡ 4 × 3 in front of Ξ in equations (128) counts the different localizations for the derivativef (see equation (112)) on the vertex itself, as pictured on Figure 5 below.
Another expression forγ comes from the Ward identity, equation (107), namely Obviously ι 2,0 = −4ι 3,2 (where the prime means the derivative with respect toγ), so that the equation foṙ γ reduces to: From β 4 given by equation (124), and equation (128), we can deduce u 6 in terms of u 4 andγ dynamically along the RG flow, , which, from the Ward constraint (3) can be translated as a function onγ only. At this level of approximation, the problem is then completely closed. The two parametersγ and u 4 fix u 6 , which fix u 8 and so one. This conclusion highlights two points. First, the role played by the derivative couplings, second that the improved local potential parametrization (102), which involve an infinite number of couplings can be, in fact, reduced to a one-dimensional manifold. Obviously, enlarging the theory space with more derivative and/or disconnected interactions, we lost this property. Moreover, note that we neglected the flow of the derivative couplingΞ ≈ 0.
We now move on the essential motivation to build the improvement discussed in the previous paragraph: the investigation of the global fixed point solutions of the flow equations. Our strategy is the following. Setting β 4 = 0, from the linearity of the equation in the sixtic coupling, we fix u 6 uniquely in terms of u 4 andγ through a relation of the form u 6 = F (u 4 ,γ). Moreover, from the first Ward identity given by Proposition (3), u 4 andγ are not independent, u 4 =γ/ι 0,2 , therefore : Explicitly: .
(133) Inserting these relation into equation (128), and settinġ γ = 0, we deduce the following: Proposition 5. In the large N limit, all the fixed points of the improved LPA have to be solution of the following equation: This equation can be solved numerically. One may expect that the dynamical definition of u 6 , u 6 = f (u 4 ,γ) breaks down at the fixed point, because both β γ and β 4 vanish at this point. It is not hard however to show that: Lemma 3. The effective RG flow, described by the function f (u 4 ,γ) satisfies: ensuring continuity at the fixed point. Proof. The proof is elementary. Let us rewrite our set of flow equations as: and the relation between them coming from Ward identity as β γ = Aβ 4 . Using the last one, we get the explicit expression for f : Moreover, setting β γ = 0 on one hand, and β 4 = 0 on the second hand; we get respectively the two solutions: Inserting these two solutions into the expression of the dynamical coupling u 6 , we get: For a global fixed point u * 6 = u * * 6 . Therefore, without singularity of the involved coefficients, we get u 6 ≡ u * 6 .
The numerical plot of the β-function β γ is provided in Figure 6, showing three zeros. The first one, forγ = 0 corresponds to the Gaussian fixed point u 4 = 0; the other, forγ ≈ −0.19 is UV-attractive, and seems to be in qualitative agreement with the UV attractive fixed point relevant for double scaling limit. Computing u * 4 from Ward identity, we get: The critical exponent can be computed from (131), computing the derivative of β γ with respect toγ at the fixed point. Numerically, we get: We discovered a relevant direction, but in strong disagreements with the theoretical predictions for the local matrix model. This is not necessarily surprising. Indeed, the value ofγ at the fixed point is relatively large, indicating a large deviation from the local model, likely to generate important qualitative differences. Thus, taking into account the Ward identities with an arbitrary regulator generates non negligible effects, making the predictions of the RG for the local model unreliable. However, the approach has the merit of indicating a clear criterion: we expect the predictions to be closer to those of the local model the smaller the value ofγ at the fixed point is. Thus, any regularization minimizing its effects will gain in reliability with respect to investigations in the theoretical space of local models and their critical properties. Let us finally note that the flow ofγ being an effect of the regulator, it also gives an indication concerning the dependence of the results on the choice of the latter.
We considered only the minimal crude truncation in the space of derivative couplings, showing the instability of the ultralocal sector due to the Ward identity. But morally, all the derivative couplings have to contribute on the left-hand side of the Ward identities, and equation (108) can be viewed for instance as the minimal truncation of a complete equation, involving an infinite set of couplings. To be more precise, let us introduce the following graphical notation. For each derivative operator like : we adopt a graphical representation as: the dots counting the "derivative insertions". Keeping only connected interactions, the enlarged theory space then can include all the possible "dotted" interactions, Such that with this short notation, the Ward identities (108) and (115) rewrite as: and: and so one for higher derivative couplings; the dotted edge meaning sums like L N and U N , or higher momenta when dots appear along the resulting closed face. We have then to deal with a proliferating number of derivative couplings. A crude truncation over the theory space, like we considered in this section improves the result with respect to a naive ultralocal truncation. But the result seems not to be satisfactory, because of the strong dependence of the derivative couplings. A way to solve this difficulty could be to investigate the dependence of the fixed point on the choice of the regulator. If the resulting fixed point and its associated critical exponents depend slightly on the choice of the regulator, for a large range of them, this can be a strong argument in favor of the reliability of our result. This is the strategy that we will discuss in the next section. Before starting the next section, and without going into the technical details concerning the optimization of the regulator, let us provide here some important comments on the choice of the regulator and its optimization condition. For the well known exact models (we denote by exact model the solvable model with the exact solution on the flow), a regulator of the FRG analysis is said to be optimal if the corresponding fixed points and critical exponents are very close to the exact results. Note also that all the regulators must carefully check the limits Γ N →Λ = S (the microscopic action) and Γ N →0 = Γ (the effective action). In the case of the nonsolvable model, we do not have rigorous criteria to fix the choice of an optimal regulator because no comparison of the results coming from the FRG can be made with the exact results. In the case of this paper the exact result is well known in the literature [1] and therefore the comparison is possible as well as the choice of the optimal regulator. Our approach to optimization is different from the general field-theoretical approach to optimization discussed in [25], [27]- [29]. The choice of a truncation and a optimal regulator will be determined by their agreement with the exact results. Thus, we will deform the Litim regulator in the hope that the deformation parameters will be chosen in order to properly approach the exact results well known in the literature. Finally the Ward identity is also a constraint that can chart the path in choosing such an optimal regulator.

C. Dependency on the regulator: a first look
In order evaluate qualitatively the dependence of our results on the choice of the regulator, a simple way is to introduce a parametrization depending on a small number of parameters. Following [30], we consider the following parametrization: We get the effective propagator into the range a+b ≤ 2N : such that, in the integral approximation, I and: the regulator; that we expect to be a consequence of the role played in the computation of the critical exponents by the derivative coupling γ. This dependence seems to be very strong, except in two regions, in the vicinity of α ≈ 1 and α ≈ 2.5, where critical exponents reach values θ ≈ 1.5 and θ ≈ 1.2 respectively. Around these two extrema, the dependency with respect to the regulator (in the considered parametrization) is small. Note that, it is possible to fine tune the value of α to reproduce exactly the expected critical exponent. On find θ ≈ 0.8 for α ≈ 0.563. However, in this region of the parameter space, θ depends strongly on α, and the reliability is poor.
Another popular choice for the regulator is the exponential one: and numerical investigations leads to the same conclusions as for the Litim's choice: A strong dependency of the relevant quantities on the regulator.

V. DISCUSSIONS AND CONCLUSION
The Ward identities, like (108) and (115) highlight the role played by the regulator in the emergence of the derivative couplings. We can stress a parallel between flow evolution and divergence of the flow toward the derivative sector. In both cases, this is the variation of the propagator -for N or a/N ≡ x that generates the moving into the theory space, in "scales" or "momenta" directions respectively. Moreover, the two transformations are not generally independents. For a regulator of the form: we get: r N (x, y) = η N r N (x, y) − Z N x ∂f ∂x + y ∂f ∂y . (153) The first term is intrinsically associated with the RG flow; however, the second part involves derivative for the momenta, which are the generators of the momentum displacements in Ward identities. Heuristically, one may picture the global dynamics as follows (see Figure 8). Starting from a point at scale N in the full theory space, the Ward operator allows moving horizontally, at fixed N toward the derivative interactions world. In the same way, the RG map allows to move vertically, from the scale N to the scale N + δN , but due to the second terms of the right-hand side of (153), the RG transformation generates as well a horizontal displacement. This is another way to understand the instability of the local phase space, pointed out in the previous section. Therefore, and despite the accordance of our results with the expected ones, especially about the value of the critical exponent, and the apparent qualitatively small dependence on the regulator in a small range of values around α = 1; one cannot conclude that our results have anything to do with the original model, the explored region of the theory space being very far from one of the original ultralocal ones.
From the last picture, a question remains open. Can you build a RG map which is the most vertical as possible, at least for N sufficiently large, in such a way that L N , U N and their higher momenta remains small enough, such that derivative couplings can be discarded from the RG flow? This question seems to be very difficult in regard to the complex hierarchical structure of the flow equations that we discussed in this paper. A heuristic attempt to solve this problem, or at least to build a flow which remains vertical for a long time is to choose a regulator such that L N , U N vanish or become small for vanishingγ. This can be achieved for instance with a regulator of the form: where the parameters α have to be fine-tuned such that L N vanish and , u 2 4 U N becomes small forγ = 0. By solving L N = 0 we get the solution α = 2. It is easy to check that this regulator satisfies the four requirements enumerated above equation (18). The corresponding flow equations can be easily deduced from our previous analysis. The condition L N | γ=0 = 0 and |U N | γ=0 | ≈ 1 allows to keep Ξ =γ = 0 with a very good approximation along the flow for a long time, in regards to the rapidity of the convergence of the truncated expansions. Settingγ = 0, the flow equations become the following: Investigating numerically the successive truncations, for k = 6, k = 8 and k = 10 like in the section III, we get only one fixed point with one relevant direction, the details being summarized in Table 9 below. Interestingly, the convergence of the truncations seems to be improved with respect to the ones considered in section III. We get only one relevant direction, with a critical exponent matching with the perturbative result. Note that no significant improvement arises from the nonperturbative effects. Moreover, the value of the relevant critical exponent seems to be insensitive to the level of the truncation. However, the value of the corresponding coupling is in strong discordance with the expected one. We get a positive and very small value for u 4 , u 4 ≈ 0.016 for k = 6, u 4 ≈ 0.01 for k = 8 and u 4 ≈ 0.008 for k = 10; the values of the other couplings being of very small magnitude with respect to these values. This observation, as mentioned, is not presented as a rigorous way to build a solution for nonperturbative RG equation, but as a qualitative illustration of how we can deal with derivative couplings to keep the flow in the purely local sector. A more complete investigation has to be carried out on this subject. Other more sophistical methods, using, for instance, background fields to constrain the flow along the vertical direction are expected to be helpful to realize such an RG map. We keep these investigations to a forthcoming work.