Functional renormalization group for $p=2$ like glassy matrices in the planar approximation: II. Ward identities method in the deep IR

This paper, as a continuation of our previous investigation [arXiv:2403.07577] aims to study the glassy random matrices with quenched Wigner disorder. In this previous work, we have constructed a renormalization group based on the effective deterministic kinetic spectrum emerging from large $N$ limit, and we extended approximate solutions using standard vertex expansion, at the leading order of the derivative expansion. Now in the following work, by introducing the non-trivial Ward identities which come from the $(\mathcal{U}(N))^{\times 2}$ symmetry broken of the effective kinetic action, we provide in the deep IR the explicit solution of the functional renormalization group for a model with quartic coupling by solving the Hierarchy to all orders in the local sector, which in particular imply the vanishing of the anomalous dimension. The numerical investigations confirm the first-order phase transition discovered in the vertex expansion framework, both in the active and passive schemes. Finally, we extend the discussion to hermitian matrices.


Introduction
In our recent investigation [1], we have considered a stochastic complex N × N matrix M , characterized by a quenched disorder, materialized by a Wigner Hermitian matrix of size N.In the large N limit, we showed that the equilibrium states look like a nonconventional Euclidean and non-local field theory, with an effective kinetics spectrum given by the Wigner law.Following the strategy developed in [2,3], we constructed a renormalization group based on a partial integration procedure of degrees of freedom, accordingly with the nontrivial notion of scale provided by the effective kinetics.There are three main difficulties with this effective field theory method in the point of view of the nonperturbative renormalization group techniques: 1.The intrinsic scale dependencies of the canonical dimension.
2. The non-local nature of the interactions.
3. The absence of an easily summable leading order sector.
The third point is specific to the matrix case, vectors and tensors being essentially branched polymers at the leading order [4,5].Note that the restriction to the leading order is enforced because we use the Wigner law for the effective kinetics spectrum.For random matrices, the leading order sector corresponds to the planar Feynman graphs, encoding dual triangulations of surfaces with zero genus [6].The first point is because the momenta distribution is a power law only asymptotically.To deal with these difficulties, we essentially focused on the deep IR, where the shape of the momenta distribution matches with a three-dimensional Euclidean field theory.Because of the second point, many popular methods usually considered in the functional renormalization group literature fail [7].In the previous work, we examined the relevance of the vertex and derivative expansion, whose relevance for nonlocal field theories has been established for tensorial field theories [8][9][10][11][12].In this framework we considered two different scaling, that we called respectively "active" and "passive" schemes, and provided some evidences in favor of a first-order phase transition between a dilute phase and a condensed phase for eigenvalues of the Hermitian matrix χ = M † M .
In this paper we are aiming to go beyond the vertex and derivative expansions, exploiting the non-trivial Ward identities arising because of the non-trivial kinetics, that break the global unitary invariance (U(N )) ×2 of the interactions.It is well known that Ward identities play a significant role in the renormalizability proofs and the construction of the renormalization group for Gauge theories [13,14].Ward identities have also played a significant role in the tensorial field theories context [15,16], characterized by similar non-localities as the model we consider in this paper.For these theories, the anomalous dimension receives corrections at one loop order, and the standard local vertex expansion fails to take into account the momenta dependency of the effective vertices, which, in contrast with ordinary field theories are significant.
The outline of this paper is the following.After a short recalling of definitions and main results obtained in the previous paper in sections 2 and 3, we derive and investigate the Ward identities in section 4 both for complex and hermitian matrices, in the deep IR.The main statement enforces the evidence in favor of the existence of a first-order phase transition, and we point out some opened issues in section 5. Some additional materials are given in Appendices A, B, and C to make the paper as self-consistent as possible.

The model and definitions
Let us briefly summarize the definition of the model.We consider a complex stochastic matrix M (t) ∈ C N×N , whose entries evolves accordingly with the stochastic process: where B(t) is a Gaussian white noise, with zero mean and variance, and H reads as: where both J and K are Hermitian Wigner matrices [17,18] with the same variance σ 2 .
In the large N limit, the spectrum of the Wigner matrices becomes deterministic, given by the well-known semi-circle law, and for any suitable test function f , we have the equality for N → ∞: dλ µ W (λ)f (λ) , (2.4) with: Then, at equilibrium, the large N partition function Z C eq [L, L] reads setting T = 2: where: and the dot product "•" is defined as: A • B := dt N i,j=1 A ij (t)B ij (t) .The resulting Boltzmann type partition function must be rewritten in a clever way introducing the two momenta, positives in the large N: together with the mass: m := a 1 − 4σ, such that the propagator becomes: C p 1 p 2 := E −1 (p 1 , p 2 ) , where E(p 1 , p 2 ) := p 1 + p 2 + m is the energy of the mode (p 1 , p 2 ).Then, in the large N limit, the equilibrium partition function looks like (but differs from) an ordinary Euclidean field theory, with nearly continuous momenta p.
Perturbation theory can be investigated using the ordinary Feynman graphs machinery.We materialize vertices of the theory as bipartite (bi)-colored completes graphs, called bubbles, such that black and white noises correspond respectively to fields M and M , hooked with colored edges materializing Kronecker delta between their indices.As an example: (2.9) A Feynman graph then looks like a (tri) colored bipartite regular graph, with dotted edges hooked between black and white nodes, corresponding to the Wick contractions (see Figure 4).It has to be noticed that a Feynman graph for such a theory is more than a set of vertices and edges, but shares an additional structure, the faces, for which we recall the definition: Definition 1 A face is an alternate cycle made of dashed and solid edges.It can be closed or open, and closed faces, corresponding to a complete trace, share a factor N.
Then, a typical Feynman amplitude A(G), depending on some Feynman graph G scale as where V k (G) denotes the number of vertices of type "2k", and F (G) the number of (closed) faces.
We considered two different strategies for constructing renormalization group flow, which we called respectively "passive scheme" or "scaling 1" and "active scheme" or "scaling 2".For scaling 1, we consider the Λ dependent family of models: where , and the propagator C (Λ) reads: where χ Λ (x) is the windows function: which is equal to 1 inside the interval (Λ, 4σ) and 0 outside.As explained in [1] it is suitable to view θ(x) as the limit of a sequence of smooth functions θ n (x), which converge weakly on the Schwartz space toward the step function θ(x) as n → ∞, and we assume to take implicitly the limit after the limit N → ∞ to replace sums by integrals.For instance: Then the effective average action: where R Λ := (C (Λ) ) −1 , interpolates between U int and: where Γ[Φ] is the full effective action (including all quantum corrections) and in these equations Φ p 1 p 2 is the classical field: (2.17) The flow equation can be deduced by taking the derivative with respect to Λ, and we get: where W Λ := ln Z eq [L].The local truncation, at the leading order of the derivative expansion, is then defined by the expansion: where from the initial condition z(4σ) = 0.The flow equations can be computed by taking the functional derivative with respect to the classical field.For µ 2 for instance, we get: or, introducing the renormalized dimensionless couplings: we get: In "scheme 2", we add in the classical Hamiltonian, the regulator such that, in order to interpolate with the classical Hamiltonian (for k → ∞) and the effective action Γ as k → 0. we impose to the regulator the following structure: where the function f (x) is inspired from the Litim's one: and Z is the wave function renormalization.Note that this regulator is something like a mixing between the standard Litim regulator and the Callan-Symanzik regulator [19].
The equation describing how the effective action Γ k changes as k changes is the so-called Wetterich-Morris equation [7,20], which reads for complex matrices: where Ẋ := k dX dk . (2.27) In this scheme, and at the leading order of the derivative expansion, it is suitable to make use of the following ansatz: From which, we easily deduce the flow equation for the dimensionless couplings: Explicitly, we get: (2.31) with: Remark 1 Note that, except for the anomalous dimension, which depends spuriously on the truncation, equation (2.30) is indeed true beyond local expansion, as the IR cut-off k is small enough [21].

Evidences for a discontinuous phase transition
In the rest of this paper, we essentially focus on scheme 2, and we will briefly summarize the evidence in favor of a first-order phase transition.For a sextic truncation, we get essentially two reliable fixed points in the deep IR, which we call f p1 and f p2, and which are respectively purely attractive and purely repulsive.The first order phase transition occurs along the solid red lines joining the two fixed points, at the point T with coordinate Y c ≈ 1.1035.We choose for order parameter the properties of the eigenvalue distribution of the matrix χ = M † M , assumed to be confined in the interval [a, b] (see Appendix B).Relevant solutions require a = 0, and the results, along the red curve, are summarized in Figure 3.The Figure shows that, as we move toward f p2, the eigenvalue distribution spreads out, until a critical value Y c , at which no finite solution for b exists.Hence, along this curve, b spontaneously jumps toward a finite value and plays the role of the (discontinuous) order parameter.

Beyond vertex expansion: closing hierarchy at equilibrium
In the previous section, we considered both derivative expansion and vertex expansion, by disregarding the role of symmetries and Ward identities.Taking into account the relations coming from Ward identities in the planar sector (because we implicitly consider the large N limit), we will: 1. Improve the computation of the anomalous dimension, taking into account the intrinsic external momenta dependency of the quartic vertex.
2. Close the hierarchy, expressing the effective local sextic in terms of the quartic and quadratic ones.
Note that this technique has been also considered in the context of tensorial field theory, in the leading sector (melonic), and beyond [15].

Symmetries and modified planar Ward identities
let's look back at the definition of the large N equilibrium partition function 4.96, with the Hamiltonian 4.97, that reads in terms of the p indices: where the Hamiltonian includes regularization, like for instance in 2.23, blinded in the kinetic kernel E k (p 1 , p 2 ) (the parameter k plays the role of IR cut-off as the UV cut-off and this parameter occurs in equation 2.23 previously).
Obviously, the Hamiltonian (4.1) is not invariant under U(N) 2 transformations: M → U M V , with U, V some N × N unitary matrices.However, the partition function 4.96 is invariant by construction, because we integrate over all matrices, and the Jacobian of the transformation is: We then have: Because the U(N) 2 transformations acting on the left and on the right are independent, we can consider infinitesimal transformations "on the left" or "on the right".For instance, we can consider the left transformation given by: where ϵ is an infinitesimal Hermitian matrix: ϵ = ϵ † .The interaction term is invariant by construction, and because the integration measure is also invariant, the variation of the partition function then comes to the kinetic and the source terms.Let us compute them separately.
Source term.The variation of the source is easy to compute from the definition, we have for instance at first order in ϵ Such that the variation of the source term is finally: Where we introduced the symbol "δ ϵ " as the variation at order 1 in ϵ.

Kinetic term.
In the same way and in first order, the variation of the kinetic term reads: (4.12)Then, the computation of the variation of the identity (4.3) at first order leads, because the cancellation has to be independent of ϵ:

Proposition 1 The partition function obeys the following differential equation (Ward identity):
where: The previous Ward identity can be easily expressed in terms of the connected correlation functions.Indeed, Z C eq =: exp +W C eq , and the connected correlations function are the n-th derivatives of the free energy W C eq with respect to the source field.To write down explicit expressions, we introduce the notation ⃗ p =: (p 1 , p 2 ) for the pair of indices shared by the source field L p 1 p 2 ≡ L ⃗ p .We furthermore make use of the norm: The connected correlation functions are then defined as: Note that odd correlation functions vanish by construction in the symmetric phase, where the classical field vanishes.More precisely, we define the symmetric phase as follows [15]: We call the symmetric phase the region of the phase space where an expansion of the classical action Γ k in the power of field is convergent enough.In this region, all the odd vertex functions vanish (including the classical field itself), and the two-point function The Ward identity (4.104) provides an infinite number of relations between observables, that we can obtain by taking successive derivatives with respect to the classical field Φ.
Note however that we have to impose the symmetric phase condition at the end, after the computation of all the derivatives.Finally, it is suitable to rewrite the Ward identities in terms of the 1PI vertex functions instead of connected generating functions, and we rewrite the Ward identity (4.104) as: where ⃗ p := (p 1 , p 2 ) and ⃗ p ′ := (p ′ 1 , p 2 ), and: where I is the identity matrix in the momenta space.In the rest of this section, we need only two Ward identities, obtained by taking respectively the second and fourth functional derivative of (4.105) with respect to the classical fields Φ and Φ. Taking the second derivative, and imposing at the end the symmetric phase conditions, we get after a straightforward computation (see for instance [15,22]): where ⃗ q, ⃗ q ′ are the external momenta i.e. the momenta of the variables with respect to which we take the derivative.Furthermore, contiguous repeated vector indices assumed are summed, and we made use of the identity: We focus on the planar approximation, in the symmetric phase, which imposes some conditions on the effective vertices Γ (2n) k .In the perturbative point of view, the effective vertices admit an expansion indexed by 1PI Feynman diagrams i.e. a set of vertices, edges, and faces (see definition 2), such that edges correspond to Wick contractions (see Figure 4).The edges can be contracted according to the following rule: In this section, we focus on the planar sector, which is the relevant one in the large N limit [6].The relevant Feynman diagram in this limit is the planar one, whose dual graphs are zero genus.Interestingly for our purpose, in the planar sector, all the boundary diagrams are connected, and we have the following statement: The proof can be done recursively and is given in Appendix A. We introduce the graphical notation: where "perm" indicates that we sum all the permutations of external momenta.This notation shows explicitly the structure of the internal paths {ℓ i } hooked to external nodes, and the function γ whereas previously the dotted edges materialized the effective propagator and the little gray circles materialized the operator δE k .It is easy to check that the first diagram creates a face (therefore a global factor N) whereas the second one does not.Hence, in the planar approximation, we keep only the first contribution.Finally, let us define the self energy Σ k (⃗ p ) as: Γ (2) Then, because of the definition of δE k , the previous equation (4.23) becomes finally at the leading order: .29) where the following constraint for external momenta is assumed: We will also need the Ward identity at the next order, relating the vertex functions at 6 and 4-points.To this end, we take the fourth derivative of (4.105) two times with respect to Φ and two times with respect to Φ.We denote respectively by (⃗ q 1 , ⃗ q 2 ) and (⃗ q ′ 1 , ⃗ q ′ 2 ) the external momenta shared by fields Φ and Φ, and we get: where external momenta are such that: assuming p 1 ̸ = p ′ 1 ̸ = q ′′ ; the configuration of momenta for the sextic effective vertex is shown in Figure 6.We call second modified Ward identity the equation (4.31).

Local approximation -Scheme 2
In this section, we investigate the Ward identities in the vicinity of the local potential approximation.This approximation usually focuses on the first order of the derivative expansion and neglects the momenta dependency of the effective vertices.In this case, the non-local structure of the interaction requires considering also the derivative of the effective vertex with respect to the external momenta [15].The Ward identities computed in the previous section allow computing the derivative of the vertices.More precisely, we will see that the first and second modified Ward identities introduced in the previous sections provide the full expression for the anomalous dimension in the local sector but also allow closing the hierarchy around the quartic sector.

First modified Ward identity.
First, let us consider the limit p ′ 1 → p 1 → 0 in equation (4.29).Then by setting q 2 = 0 we have: where we use the definition of Z(k), and we replace the discrete sum with a continuous integral involving the asymptotic distribution ρ(q).We then use the modified Litim's regulator, as defined in (2.24).This in particular implies that the windows of momenta allowed by the derivative with respect to p 1 in the previous equation are the same as the windows of momenta allowed by Ṙk in the flow equation, and the approximations used in the second case has then to be valid in the first case.Hence, because we focus on the local approximation, it makes sense to replace γ k,q000 by γ k,0000 : This approximation can be also motivated, without the local assumption, from the observation that, for k small enough, the allowed windows of momenta selected by the regulator are closed to q = 0.For the first term on the left-hand side, this approximation is difficult to motivate because far enough from the IR regime, the power counting can be very different, especially using scheme 1 for the scaling.Note however that the flow equation for Γ (2) Once again, if k is small enough, the ⃗ q dependency of Γ q can be neglected, setting ⃗ q = ⃗ 0 as the selected windows of momenta is closed to zero.Then: with the definition: NL n (k) := q G n (0, q) Ṙk (0, q) .(4.37) Hence, we have: To compute the sum over q of the last term, we can use the truncation to express G k , because the selected windows of momenta are the same as in the flow equations (it is not an additional assumption).The first term furthermore can be rewritten as: where s := ln(k).The UV contributions, unaffected by the regulator are essentially the same for G k (0, q) and G k+dk (0, q).Hence, we make the crucial assumption that, for k small enough, it makes sense to express G k with the IR truncation.Explicitly, assuming k is small enough (IR regime): where the functions F ± are defined as follows: The sign ± refers to the sign of the mass ũ2 , and the two functions are represented in Figure 7.For ũ2 > 0, we have two branches of solutions, which have the same limit for ũ2 → 0: Now, computing the derivative of the sum (4.40) with respect to s := ln(k), we get: Each term can be easily computed, and using the flow equation (2.30) for ũ2 , we get: , (4.46) where: Furthermore, L n (k) can be computed exactly: where J ν has been defined before in (2.32).Moreover, in the deep IR, we have explicitly: and we get finally, at the leading order in k: Finally, from the definition of R k : we have, again in the deep IR: (4.52) Furthermore, the truncation implies: and the first modified Ward identity reads finally: In the deep IR, assuming ν = O(1), we have σ → ∞, enforcing the condition: The fixed point solutions furthermore impose ũ4 ∼ σ3/2 , as we recalled in the previous sections 3. We then have: .
(4.56)This solution is the effective quartic coupling in the deep IR, and we can show its behavior on the left of Figure 8.The renormalization group flow is finally totally driven by the mass term: and its behavior is shown on the right of Figure 8.A direct inspection shows that it exists one interacting Wilson-Fisher-like fixed point for the value: Which is relevant, with critical exponents: The fixed point is indeed repulsive on both sides, and the flow reaches the positive mass region on the right, and the deep negative mass on the left, until the boundary fixed point at ũ2 = −0.5,where coupling becomes imaginary.We call FPA the fixed point on the right, near the positive mass region.We denote as FPB the fixed point on the left, on the negative region (which is also the boundary of the symmetric phase region, above this point the solution for the quartic coupling becomes imaginary).

Next to quartic order.
Local interactions of degrees higher than 4 can be expressed as well in terms of the effective renormalized mass ũ2 .Indeed, we can compute the u(W ) On the left: Behavior of the solution for the effective quartic coupling (rescaled by a factor σ3/2 ).On the right: Behavior of the beta function for ũ2 .
(4.56), and if we compare it with the flow equation (2.31), we deduce the expression of ũ6 , which can be used to deduce the expression for ũ8 from its flow equation, and so one.We then deduce in principle the expressions for ũ2n , even though concretely the expressions become rapidly intractable as n > 5.In Figure 9, we show the behavior of couplings from n = 1 to n = 5, and the behavior of the effective potential along the ũ2 axis, between fixed points FPA and FPB for a sextic truncation in Figure 10; In Figures 11-12 we can identified the corresponding eigenvalues distributions due to the confining, computed using formula (B.15) of appendix B. Starting from ũ2 ≈ −0.5, and moving toward positive masses, we show that the width of the eigenvalue distribution increases until ũ2 ≈ −0.22,where the distribution becomes unbounded.Indeed, up to this point, the potential becomes unbounded from below, as Figure 10 shows, and eigenvalues are repelled from the origin.Hence, up to this point, there are no finite values for the parameter b (see formula (B.15),where we make the assumption a = 0).Spontaneously, the value of b −1 takes a finite value b −1 ≈ 0.15 for ũ2 ≈ −0.22, a discontinuity reminiscent of a first order phase transition.The results are reminiscent of the phase transition from region I to region II recalled in section 3, once again we discover a discontinuous (first order) phase transition with respect to the parameter b, between a condensed to a dilute phase for the eigenvalues.
Claim 1 There exists a first-order phase transition between a high-temperature dilute phase and a condensed phase, with order parameter b −1 jumping from b −1 = 0 to b −1 ≈ 0.15 for a temperature scale ũ2 ≈ −0.22.
The values (ũ 2 , ũ4 , ũ6 ) ≈ (−0.22,0.77, −0.37) have to be compared with the corresponding values (ũ 2 , ũ4 ) ≈ (−0.22,0.92, −0.49) of the sextic truncation.Now it becomes necessary to try to respond if the above results are reliable.This question will be the subject at the end of this subsection and of the next one.Note that on some regions along the mass axis, the sign of the larger order coupling changes from n to n + 1, and the shape of the effective potential is poorly reliable for some points on the interval between fixed-points -see Figure 13 for instance.It has however been noticed that the region where the larger coupling is positive increases with the degree of the truncation, and moves, from FPB, toward the positive mass region.Figure 9 shows this explicitly,  and it can be checked that this phenomenon is repeated for higher-order truncations.Furthermore, in some points along the ũ2 axis, the rate of convergence is faster enough to make a reliable prediction.This is the case, for instance, for the potential near the fixed point FPB, as shown in Figure 14 up to order n = 5.Hence, we conclude that despite the instabilities observed at first orders of the vertex expansions, the results are reliable enough and make the prediction about first-order phase transition robust; higher-order couplings simply shift from right to the transition point, as Figure 10 shows.In Figure 15, we can see the eigenvalue distribution for truncations up to order 5 at FPB and the dependency of the transition point are summarized in Table 1.Note that the eigenvalue distribution around FPB remains close to the sextic prediction.Finally, note that the quartic order predictions are very poorly reliable because they predict only a continuous phase transition between a low-temperature regime where eigenvalues are spread around a non-zero vacuum and a high-temperature regime where eigenvalues are spread around the zero vacuum.
We conclude this part by making the following important remark.The effective couplings ũ2n σ3n/2 have the power count (3 − n)/2.Thus, by adding the couplings ũ8 and ũ10 , we have included irrelevant effects, in competition with other effects of the same order (relative to the Gaussian counting), derivative couplings of higher order for example.Let us remark that it seems difficult under these conditions to have complete confidence in our estimates of high-order couplings, in any case, this does not constitute a sufficient argument to doubt the satisfaction 1, which moreover seems to agree with the predictions of vertex expansion.Table 1: Summary of the properties of the transition point for truncations up to order n = 5.As expected, the transition remains first order since n > 0, and the transition point moves toward positive mass as n increases.

Second modified Ward identity.
In this section we then consider the sextic coupling ũ6 , which can be expressed using the second modified Ward identity (4.32).Setting p ′ 1 → p 1 , this identity reads: k,q,0,0,0,0,0 − 2(γ Furthermore, the expression of the derivative of γ (4) k with respect to the external momenta can be obtained from the flow equation, because of the constraint ν = 0. Indeed, coming back on the computation of ν, but taking into account the dependency of the effective vertex on the external momenta, we have: Ṙk (q, 0) Ṙk (q, p 1 ) or explicitly, because ν = 0: leading to the remarkably simple result: where we used the condition, coming from the chain rule: Now, let us return to the computation of the Ward identity (4.60).As explained before, the terms proportional to ∂ p 1 R k (0, q) can be computed using the local truncation, and we get: k,q,0,0,0,0,0 − 2(γ

.65)
The remaining terms can be computed using the same strategy as for the quartic case (see (4.38)).Then, we construct a contribution like q G 2 k (0, q) Γ( 2) k (q) from the first term (with γ (6) k ), which is exactly compensated by the same contribution coming from the term with two γ (4) k .Explicitly we get: k,q,0,0,0,0,0 − 2(γ and the Ward identity reads finally (remembering that Z = 1): k,p 1 ,0,0,0 | p 1 =0 =: ξ 1 , and we then deduce the second quantity d dp 1 γ (4) k,p 1 ,p 1 ,0,0 | p 1 =0 =: ξ 2 from equation (4.63).The behavior of ξ 1 and of the solution ũ4 /7 is shown in Figure 16, and are again continuous at ũ2 = 0. Obviously, the solutions do not match one with the other, and the Ward identity (4.2.3) indeed fails at the order of approximation we considered.This point requires a deeper analysis which we reserve for a forthcoming work.We suspect that the computation of sub-leading interactions like ξ 1 could require going beyond the derivative expansion, including for instance higher order derivative interactions, like for instance ∼ Z 2 (p 2 1 + p 2 2 ) in the 2-point function truncation (see [23]).In that way, the Ward identity will be used to fix the value of this interaction and so one; what could be improved also the computation of local interactions, for which we pointed out in the previous section the poverty of the predictions for high orders, which could be attributable to the fact of having ruled out the effect of derivative couplings.As mentioned before, we plan to address this question in the future.

Local approximation: Scheme 1
Here, we consider the Ward identities for scheme 1.In this case, the variation δE k ≡ δE Λ reads, in the limit p ′ 1 → p 1 : and we have: Due to the fact that χ Λ (p 1 ) do not vanish only on the interval p 1 ∈ [Λ, 4σ] (as n → ∞), we must have χ ′ Λ (p 1 ) ≃ 0 as soon as Λ ̸ = 0. Hence, the Ward identity (4.23) reads in that case, with the momenta configuration (4.30), but setting p ′ 1 → p 1 → 0: Λ,p 2 ,0,0,0 χ Λ (p 2 ) (Γ Because of the truncation (2.19), we have: To compute the sum of the left-hand side of equation ( 4.70), we use the same strategy as in the previous section.Because of the flow equation (2.20), and using the same arguments as before, we have, in the deep IR: Λ,p 2 ,0,0,0 L 2 (Λ) , ( where: then: Λ,p 2 ,0,0,0 χ Λ (p 2 ) (Γ which can be rewritten as: The second term can be computed from the same strategy as we used for the computation of the loops L n (Λ), we get: and, because of the definition (2.21) for μ2 : The first term can be computed from the same observation as before: Only the vicinity of p 2 = Λ is relevant.Hence, in the deep IR, we make use of the truncation, in a region where it is assumed to have good convergence properties, and we get after a tedious calculation: which can be expanded in power of Λ, and keeping only the leading orders, we have: then: Finally, because of the asymptotic expression for Λ ≪ 1: , the leading order Ward identity (of order Λ −1/2 ) reads finally, for m = 0 as: which implies that η = 0 , ( and then we can fix the asymptotic value for z, we choose: z = 1, and the next to leading order term of Ward identity (of order O(1)) is: which can be solved for μ2 as: Comparing with the flow equation (2.22), we deduce an explicit expression for μ4 : For m ̸ = 0, the Ward identities can be investigated in the deep IR again, assuming that | m| ≫ 1.At the leading order, we get, again η = 0, but the value for z is fixed by the term of order O( m0 , Λ 0 ), which reads:

.88)
The positive solution z + ≃ −0.38 ensures that 1 + z + ≈ 0.62 > 0, and we retain only that one solution.At the order O( m−1 ), we have: then: The Figure 18 summarizes the results for m ̸ = 0.

Short comment about Hermitian matrices
Let us discuss briefly our construction for Hermitian matrices.The stochastic quantization of Hermitian matrices uses the same kind of equation as (2.1), but includes the constraint: along the dynamics.In Appendix C we provided some details about the derivation of the corresponding Fokker-Planck equation, and as expected, the equilibrium theory is described by the quenched path integral: where here dM is the standard measure over N × N Hermitian matrices manifold, see Appendix C. In the large N limit, and working in the eigenbasis of the disordered matrix, the equilibrium partition function becomes (setting again T = 2): where explicitly Note that in this expression, and in contrast with the complex case, λ 1 and λ 2 are eigenvalues of the same Wigner hermitian matrix.We focus on scheme 2 for this section, and in replacement with the local truncation (2.28), we have: (4.98) Moreover, we keep the Litim regulator 2.24, but we add a factor 1/2 in front of the definition of ∆S k .With this parametrization of the theory space, the β-functions can be easily computed, and we get for instance, for the 2-points function (see [1] for more details): the definition of the dimensionless couplings being the same as for complex matrices.The equation is the same as for a complex matrix because of the normalization of the vertices.
The Ward identities can be deduced as for the complex case, by requiring the invariance of the partition function under some unitary transformation M → U † M U , where U ∈ U(N).Considering some infinitesimal transformation U = I + iϵ, with ϵ = ϵ † , the matrix field transforms as: Let us compute the variation of the kinetic action.We get, at the first order in ϵ: For the source term, we get in the same way: Hence, the first order Ward identity reads finally: The partition function obeys the following differential equation: Using the same notations as in section 4.1, we get: where ⃗ p := (p 1 , p 2 ) and ⃗ p T := (p 2 , p 1 ).Taking the second derivative with respect to the classical field Φ, ∂ 2 /∂Φ ⃗ q ∂Φ ⃗ q ′ , and imposing the symmetric phase conditions: This equation can be investigated as before, and the results are very similar.In particular, we get again that anomalous dimension have to vanish for the asymptotic flow in the deep IR regime.This condition has been found moreover in the tensorial field theory context, as a requirement for fixed point condition [24,25].

Conclusion and open issues
In this paper, we have presented a method going beyond the vertex expansion and exploiting Ward identities arising consequence of the symmetry breaking of the effective kinetics of the model.This second method imposes strong constraints on the renormalization group, notably imposing the cancellation of the anomalous dimension asymptotically towards the IR.For equilibrium theory this method confirms some predictions of vertex expansion, but not all.In particular, we enforce the evidence in favor of the existence of a first-order phase transition, resulting in the confinement of the eigenvalues seems to be inevitable.
However, this conclusion remains partial, and the proposed formalism must still be completed by subsequent work.In particular, the role of higher-order derivative interactions must be carefully evaluated and seems to play a determining role when we involve interactions of a higher order than 6.This open issue will be the subject of future work.Moreover, some aspects have been ignored in this formalism.In particular, the behavior of the flow in the IR but not deep IR, where the specific shape of the momenta distribution could play a relevant role, has been ignored.Furthermore, the role of 1/N corrections to the Wigner law, which could become significant in the deep IR as k, Λ ∼ 1/N, especially for sextic interactions, have been ignored as well and will be also the subject of forthcoming work.Finally, we only focused on first-order Ward identities, and it has been shown recently that next-to-leading-order Ward identities could be required for that kind of theory [26].B Computing eigenvalue spectra of M † M .
In this section, we shortly review the solution of the problem to determine the eigenvalue spectra of χ = M † M , in the large N limit, from the knowledge of the distribution law for the complex matrix M .We assume the partition is of the form: where V expands in power of χ: To perform the change of variable χ = M † M in the path integral, we use the clever expression of the identity: where the integral is performed on the space of Hermitian matrices.The Dirac delta function can be furthermore rewritten in the Fourier representation for the N (N + 1) independent components of χ: where the integration is performed on the hermitian matrix λ1 ; furthermore, it is suitable to exploit the unitary invariance of the integral over M to choose χ diagonal.Finally, the partition function reads: The Gaussian integral over M can be easily performed, , the Jacobian of the transformation is: and: Z ∝ dχ e −µ 1 Tr χ−Tr V (χ) .(B.11) In the large N limit, the integral is dominated by the saddle point.Using the polar decomposition of the integration measure over Hermitian matrices [6,17]: where Λ is a diagonal and positive definite matrix with eigenvalues λ i : dΛ ≡ i dλ i , ∆(Λ) := j>i (λ j − λ i ) is the Vandermonde determinant and dU is the Haar measure over the Unitary group.It is furthermore suitable to rescale χ → Nχ, such that the saddle point equation reads: This equation can be investigated in the continuum limit from Tricomi's formula [17,27,28].Indeed, denoting as µ(λ) the large N eigenvalues distribution, the previous equation reads: If the spectrum is bounded, the solution for µ(λ) as the integral form:  The formula (B.15) can be generalized in cases where we have only one hard wall or no hard wall.In the first case, assuming for instance we have a hard wall at x = b, the solution is [28]: and in the second case:

C Fokker-Planck equation for Hermitian matrices
In this appendix, we propose for the reader unfamiliar with the stochastic formalism to derive the Fokker-Planck equation for the Hermitian matrix (more detail about the general formalism can be found in the standard reference [13,14]).where the averaging is over the noise.It is easy to check that it is normalized to 1 with respect to the integration measure over the Hermitian matrix space: where dM := i dM ii j<i dReM ij dImM ij .We will prove that P [M, t] satisfy to the equation: Let q ∈ R N 2 the real vector with N 2 components, such that: the probability then can be rewritten in terms of the variable q as P (q, t) := ⟨δ(q − q(t))⟩.Now, because the system has no memory (the noise is randomly distributed with the same probability at each time), the probability as to be a Markovian process: P [q, t; q 0 , t 0 ] = d N 2 q ′ P [q, t; q ′ , t ′ ]P [q ′ , t ′ ; q 0 , t 0 ] , (C.5) where in the notation P [q, t; q 0 , t 0 ] we make appeared explicitly the initial condition M (t 0 ) = M 0 .Because the origin of time is irrelevant, let us consider the probability P [q, t 0 + ϵ; q 0 , t 0 ] to jump from M 0 to M in the (expected small) time interval ϵ.In the laps ϵ, the trajectory changes as: where f is the gradient force and b is the vector in R N 2 build from the noise matrix field B as the vector q, explicitly for instance: Now, we discretize the time interval into steps of size ϵ, and we define the discrete noise b i , for time t i = i × ϵ as: The discrete probability W (ϵ, q, q 0 ) := P [q, t 0 + ϵ; q 0 , t 0 ] to jump from q 0 to q in the laps ϵ is: Now, taking the Fourier transform with respect to the variable q, we get: (C.17) The linear term in ϵ is nothing but the matrix element ⟨p|H|q 0 ⟩, then, taking the inverse Fourier transform, we get the matrix element ⟨q|H|q 0 ⟩, where: is the unitary invariant Laplacian operator on Hermitian matrices [29].We then proved that: ∂P [q, t] ∂t = HP [q, t] , (C.21) with, returning to the matrix notation and using the explicit expression for the driving force (C.7), we get: where we make use of the identity

Figure 1 :
Figure 1: A typical Feynman graph involving four external edges, two quartic vertices, a sextic, and a quadratic vertex.

Figure 2 Figure 2 :
Figure 2: Qualitative illustration of the different phases.In the region in the interior of the thick black line, the denominator of ν is negative.

Figure 3 :
Figure 3: On the left: Evolution of the effective eigenvalue distribution along the red curve.On the right: Behavior of the inverse of the solution for the edge boundary with Ȳ , along the red curve.Ȳc ≈ 1.0135.

Figure 4 :Definition 3 Figure 5 :
Figure 4: A typical Feynman graph with 4 vertices, 4 external edges, and 6 internal edges.The red path with arrows materializes a typical external face.

Lemma 1
The Feynman expansion for Γ (2n) k can be indexed by planar diagrams having the same boundary diagram, corresponding to the single connected trace invariant vertices.

Figure 7 :
Figure 7: On the left: The two branches of solution, F ± .On the left: Continuity of the derivatives F ′ ± .

Figure 10 : 5 Figure 11 :Figure 12 :Figure 13 :
Figure 10: Evolution of the effective potential from FPB, along the mass axis.The abscissa x is a typical eigenvalue λ =: x 2 ≥ 0 of the positive Hermitian matrix χ := M † M .

Figure 14 :Figure 15 :
Figure14: The convergence of the potential in the vicinity of the two fixed points, for quartic (blue curve), sextic (purple curve), octic (yellow curve) and dectic (green curve).

Figure 16 :
Figure 16: On the left: The solution ξ1 (ũ 2 ) rescaled with a factor σ−3/2 , in the negative and positive region.On the right: Comparison with the solution ũ4 /7.

Figure 17 :
Figure 17: On the left: Behavior of the solution for the effective quartic coupling.On the right: Behavior of the beta function for μ2 (m = 0).

Figure 21 :
Figure 21: All the 1PI decompositions for the second type insertions: all the gray blobs denote 1PI components.

18 )
15) where C is adjusted from the boundary and normalization conditions.Consider for instance a purely quadratic confining potential, forf (λ) = µ 1 /2, we have, choosing C such that µ(b) = 0, µ(λ) = µ 1 λ(b − λ) 2πλ .(B.16)Note that it is suitable to set a = 0, and the boundary b is then fixed by the normalization condition:Not surprisingly, the distribution looks like a Marchenko-Pastur distribution and the typical shape is shown in Figure22.