Tensorial Gross-Neveu models

We define and study various tensorial generalizations of the Gross-Neveu model in two dimensions, that is, models with four-fermion interactions and $G^3$ symmetry, where we take either $G=U(N)$ or $G=O(N)$. Such models can also be viewed as two-dimensional generalizations of the Sachdev-Ye-Kitaev model, or more precisely of its tensorial counterpart introduced by Klebanov and Tarnopolsky, which is in part our motivation for studying them. Using the Schwinger-Dyson equations at large-$N$, we discuss the phenomenon of dynamical mass generation and possible combinations of couplings to avoid it. For the case $G=U(N)$, we introduce an intermediate field representation and perform a stability analysis of the vacua. It turns out that the only apparently viable combination of couplings that avoids mass generation corresponds to an unstable vacuum. The stable vacuum breaks $U(N)^3$ invariance, in contradiction with the Coleman-Mermin-Wagner theorem, but this is an artifact of the large-$N$ expansion, similar to the breaking of continuous chiral symmetry in the chiral Gross-Neveu model.


Introduction
The large-N limit has a long tradition in field theory as it allows to restrict the perturbative expansion to a specific class of Feynman diagrams. In the case of vector fields with N components, the restricted class of diagrams is so simple (essentially just chains of bubbles diagrams) that it generally allows to solve many models, either diagrammatically or by saddle point methods [1]. For the next natural generalization in the sense of linear algebra, the case of N × N matrix fields, the large-N limit has led to numerous results in zero dimensions, in particular because of its connection to two-dimensional quantum gravity and string theory [2]. However, matrix field theories in two or higher dimensions turn out to be much more complicated than the vector case, as at leading order in the large N expansion one encounters all the planar diagrams. Somewhat surprisingly, it turns out that going one more step up in the rank, i.e. considering tensor fields, things simplify again, although not to the level of vector fields, thus making tensor fields the good candidates for models with a new and manageable, yet non-trivial, large-N limit.
Although tensor models in zero dimensions were introduced long ago as a generalization of matrix models for higher-dimensional quantum gravity [3,4], it is only more recently that the possibility of studying their 1/N expansion was discovered [5,6,7], and the class of leading order diagrams identified: the melons [8,9]. The latter are a class of diagrams with a tree-like combinatorial structure; they are therefore easier to sum than the planar diagrams, but are definitely richer than the chains of bubbles. Such combinatorial structures show up naturally in a variety of other contexts, and have for instance played a central role in renormalization group approaches to group field theories [10,11,12,13,14,15,16]. More recently, it is precisely the structure of the melonic graphs that has led to the observation that the Sachdev-Ye-Kitaev (SYK) model [17,18,19,20] and its tensorial generalizations [21,22,23,24,25,26,27] have a conformal IR limit and saturate the chaos bound [28], thus establishing them as interesting models of holography. Such models live in one (time) dimension, but have also been generalized to higher dimensions [22,29,30,31]. In close relation to such developments, a new expansion has recently been proposed to investigate the properties of matrix models of quantum black holes in large D dimensions [32,33,34,35].
In this paper, we explore two-dimensional models of fermions in a tensor representation of rank three, with quartic interactions. These are a natural generalization of the Gross-Neveu (GN) model [36], in which the fermions are in a vector representation. The GN model is asymptotically free in two dimensions, and although a mass term is forbidden by discrete chiral symmetry, dynamical symmetry breaking occurs and a mass is non-perturbatively generated. We expect a similar behavior in the tensorial generalization, thus in general we do not expect to find a conformal field theory in the IR, which is the main attractive feature of SYK-like models. However, the tensorial generalizations have several couplings and therefore it is in principle possible that for some specific values of the couplings a mass is not generated and a conformal theory is found. We will investigate this possibility by analyzing the large-N Schwinger-Dyson equations first, and then by the intermediate field method.
The two-point function of SYK models enjoys an emergent time reparametrization invariance (or conformal invariance) in the large N and strong coupling limit, apparent at the level of its Schwinger-Dyson equation. A crucial ingredient is that the self-energy dominates over the free propagator in the infrared regime. For massless fermions (with standard 1/ / p propagator), a naive dimensional analysis shows that this property can only hold if d(1 − 2/q) < 1, where d is the spacetime dimension and q is the order of the interactions. Hence d = 2 is the critical dimension at which this condition starts failing 1 : the quartic couplings being dimensionless, both the free propagator and the self-energy must be taken into account in the two-point Schwinger-Dyson equation. This is the most crucial difference between Gross-Neveu tensor models and SYK-like theories, leading to substantial qualitative differences in their large N behaviors.
The plan of the paper is as follows. In section 2 we introduce our two models -a U (N ) 3 model with Dirac fermions and an O(N ) 3 model with Majorana fermions -paying special attention to the parametrization of their action in terms of invariants. In section 3, we study the large-N (melonic) Schwinger-Dyson equations for the two-point functions of both models, and provide a first analysis of the mechanism underlying the spontaneous generation of mass. In the Dirac model, we find a specific combination of the coupling constants which seems to ensure that the theory remains massless. In the Majorana model, we find also a candidate massless theory, but with a marginally irrelevant coupling constant. Relying on the introduction of a matrix intermediate field for the Dirac model, we move on to an in-depth analysis of the would-be massless vacuum in section 4. The effective potential of the intermediate field is computed, and shown to develop a symmetry breaking pattern in the regime of interest. This results in the instability of the would-be massless vacuum, triggered by a non-zero expectation value of the intermediate field. We discuss in detail the apparent symmetry breaking of the U (N ) 3 symmetry, and derive the effective non-linear dynamics governing the associated massless modes. We recall how the latter should not be understood as Goldstone bosons since any apparently broken continuous symmetry must be restored in d = 2 [37,38], and we describe the mechanism of symmetry restoration at play in this model. We finally close with a summary and a few concluding remarks in section 5. Note added: During completion of this work a paper by Prakash and Sinha [39] has appeared, in which they also study tensorial fermionic models in more than one dimensions. They consider a model for complex fermions but with symmetry group U (N ) × O(N ) × U (N ) which allows one to write the equivalent of our I 2 (or tetrahedral) interaction. They do not consider the equivalent of our I 0 and I 1 interactions, and they focus on d = 2, searching for a behavior similar to the one in the SYK model, where the free propagator term in the Schwinger-Dyson equations can be discarded in the IR limit. Although the two works are thus in most part complementary, we make a connection between the two in section 3.2, where, relying on the RG analysis of appendix C, we provide a check of their conjecture that the models with tetrahedral interaction have a weakly interacting IR fixed point in 2 − dimensions.

The models
The Gross-Neveu (GN) model [36] is a model of N massless Dirac (or Majorana) fermions 2 arranged into a vector ψ i , with i = 1 . . . N , and invariant under U (N ) (or O(N )) 3 transformations ψ i → U ij ψ j , with U ∈ U (N ) (or O(N )) and summation on repeated indices The model is renormalizable and asymptotically free in d = 2 for g > 0, and it has a discrete chiral symmetry that (if unbroken) protects it from the generation of a mass term. A mass gap is however generated non-perturbatively, as most easily seen in the large-N limit. A large-N expansion is obtained by redefining g = λ/N and keeping λ finite. Then the beta function at leading order in 1/N is Since in the large-N limit such a beta function is valid also at large λ, one encounters a Landau pole in the running coupling. The latter is interpreted as the appearance of a tachyon due to an instability of the invariant vacuum. An analysis of the effective potential for the intermediate (or Hubbard-Stratonovich) field reveals that the stable vacuum breaks discrete chiral symmetry and leads to a dynamical mass generation.
Here we will study generalizations of the GN model, where we have N r fermions arranged into a rank-r tensor ψ i 1 ...ir , such that the action is invariant under U (N ) r (or O(N ) r for Majorana fermions) transformations There are many possible invariant actions; here we restrict to the case of d = 2 and quartic interactions, in order to have power-counting renormalizability. For simplicity we will restrict also to r = 3. We will refer to the invariance (3) as flavor symmetry, and we will also introduce a discrete permutation symmetry, which we will refer to as color symmetry [9,41]. Lastly, like it is done for the usual GN model, we will introduce also a discrete chiral symmetry as well as its continuous counterpart, which we define in appendix A, see (165) there. We split the action of our model into a free and an interacting part: where S free is bilinear in the spinor fields, while S int is quadrilinear. The construction of bilinear and quadrilinear terms which are invariant under Euclidean and chiral symmetries is recalled in appendix B. We will see below how to construct flavor-and color-invariant interactions. Since the choice of Majorana versus Dirac fermions-and hence O(N ) 3 versus U (N ) 3 symmetryleads to a slightly different set of interactions, as well as to a different normalization for the quadratic term, we shall discuss the two cases separately starting from Dirac fermions.

Dirac fermions
In the case of complex Dirac fermions we define the free part of the action as The discrete chiral symmetry (4) forbids a mass term in the action, and thus (6) is the only bilinear invariant under all our symmetries. In GN a mass is dynamically generated at the non-perturbative level, and we will address in some detail the question of dynamical symmetry breaking for our models in the upcoming sections. U (N ) 3 -invariant interactions are constructed by pairwise contraction of indices belonging to ψ fields with indices belonging toψ fields (which transform as the complex conjugate of (3), where U (i) ∈ U (N )). In order to represent the interaction terms in a compact way it is convenient and customary to introduce a graphical representation, as explained in figure 1, where in the left and central panels the following interaction terms are represented: Here, terms inside a parenthesis have all their spinorial indices mutually contracted, see (166) in appendix B. The matrices Γ X , with X = S, V, P , are defined as The letters S, V , and P stand for scalar, vector, and pseudoscalar, respectively, in reference to the transformation properties of the associated bilinears under the rotation group, see appendix B.
The complex nature of the fields (i.e., the fact that an index in ψ always needs to be contracted with an index inψ in order to form invariants) implies that the interaction graphs are bipartite; this means that the vertices can be divided in two sets, representing the ψ andψ fields, respectively, in such a way that two vertices within the same set are never adjacent. One can immediately verifies that the first two graphs in figure 1 are bipartite, while the third one is not.  (13) (right). Each vertex represents a tensor/spinor field, a solid line with label (color) i ∈ {1, 2, 3} joining two vertices represents the contraction of the i-th tensorial indices of the two fields at its end-vertices, and a dotted line with label X represents the contraction of the spinorial indices of the fields at its end-vertices, with the insertion of a matrix Γ X .
The interaction (8) is built upon one of three possible invariants. We label them by I X, 1 where X = S, V, P denotes the spinor structure Γ X , and = 1, 2, 3 describes the tensor structure. The interactions I X, 1 have the tensor index in the -th position contracted between different bilinears. In other words, I X,2 1 and I X,3 1 are obtained as the two cyclic permutations of the color indices in the central picture of figure 1 (the color being a code name for the index location associated to an edge). It is not difficult to verify that I X 0 and I X, 1 are the only possible flavor-invariant interactions compatible with the requirements of renormalizability and Euclidean invariance. The only non-trivial step is to show that the spinor contraction for the I X, 1 interactions can always be chosen as in figure 1, rather than for example parallel to the color-1 edge in that picture. This is done by using standard Fierz identities or, in other words, the completeness relation (154) of appendix A.
One can recognize that the I S 0 interaction is nothing but the usual GN interaction in disguise, and thus it is invariant under a larger symmetry group, namely U (N 3 ). 4 In fact the kinetic term has the same symmetry as well: it is the I 1 interactions that break such symmetry down to U (N ) 3 .
In order to reduce the number of invariants, we can demand that the action be invariant under simultaneous and identical permutations of all the tensor indices, known also as color symmetry. Such symmetry requires repackaging the I X, 1 interactions as 4 More general GN models, with also the V and P interactions have been studied in [42,43,44].
In this paper we will consider both models with and without color symmetry. 5 The most general renormalizable interacting action compatible with Euclidean, chiral, U (N ) 3 and color symmetries thus contains six independent couplings, and reads 6

Majorana fermions
In the case of real Majorana fermions we define the free part of the action as where the different normalization is chosen for later convenience, taking into account that the spinors are real valued.
In constructing the interacting part of the action one encounters a new invariant, besides those we discussed for the Dirac case above. This is because here graphs representing the interaction vertex no longer need to be bipartite-recall thatψ = ψ t γ 5 , and thereforeψ transforms like ψ with U (i) ∈ O(N ), see (3). The new graph is tetrahedral (complete, four vertex), depicted on the right in figure 1, and generalizes the Klebanov-Tarnopolsky-Carrozza-Tanasa interaction to two dimensions [47,22]. In order to reduce the number of possible interactions, we would be tempted to require again color symmetry. However, we must be a bit careful with how we average over colors for the I X, 2 interactions, because the symmetric structure of the complete graph associated to them implies crucial cancellations. It turns out that there is no I 2 interaction invariant under the full color permutation group; there is however (precisely) one invariant if we only require symmetry under the alternating subgroup of the color permutation group, as we will show.
Let us introduce two sets of I 2 interactions, labelled by a sign ±. The I 2,+ interactions are defined as and the I 2,− interactions as 5 Color symmetric models are the most commonly studied in zero dimensions [45,41] but non-symmetric ones have also been considered, e.g. in [46]. 6 The minus sign is for later convenience. We recall that at finite N , and finite UV cutoff, the fermionic functional integrals are analytic around the origin. The good sign, if any, of the couplings can only be determined by the large-N analysis of the renormalized theory.
These definitions have been chosen so that, for any X ∈ {S, V, P }, ∈ {1, 2, 3} and i ∈ {+, −}: • the orbit of I X, 2,i under the symmetric group S 3 (acting simultaneously on the three indices of all four tensors) is {I X, 2,j | j = ± , = 1, 2, 3}; • the orbit of I X, 2,i under the alternating group A 3 ⊂ S 3 is {I X, 2,i | = 1, 2, 3}. Now, the interactions I X, 2,+ and I X, 2,− are not independent. Indeed, for Majorana spinors one has that (ψ a ψ b ) = (ψ b ψ a ), (ψ a γ µ ψ b ) = −(ψ b γ µ ψ a ) and (ψ a γ 5 ψ b ) = −(ψ b γ 5 ψ a ), which implies: We therefore conclude that the fully color-invariant P and V sectors are trivial. It would look like we have one invariant in the S sector, i.e. 2 I S, 2,+ , but we will now show that this is also identically zero.
Let us determine the A 3 -invariant interactions, which we choose to parametrize in terms of the vertices I 2,+ . The completeness relation (154) can be used to prove that 7 for any color . By summing over , this immediately implies: where we have defined I X 2 := I X, 2,+ . Therefore the A 3 -invariant sector contains at most two independent coupling constants. Performing an odd permutation on the color indices in (18) we find that the left-hand side is invariant, while the right-hand side picks a minus sign; therefore, both sides of the equation are zero. We are thus left with no invariants under S 3 and at most one invariant under A 3 , which we could choose for example to be I P 2 . We are now going to show that this is in fact a non-trivial invariant, and that we can choose to write it in such a way that we do not have to perform a sum over colors.
To go further, we can invoke the completeness relations (155) and (156) to find: From now on, let us drop the ± index and simply denote I X, 2,+ by I X, 2 . For any given , the Fierz relations show that we are free to parametrize our action with the basis (I S, 2 , I V, 2 , I P, 2 ). Let us denote by g 2 = (g S, 2 , g V, 2 , g P, 2 ) t the coordinate vector in this basis, i.e. the interaction action in the sector 2 is S int,2 = − X g X, 2 I X, 2 . The Fierz identities can be summarized in matrix form by 7 The addition operation on color labels is to be understood modulo 3.
We can furthermore introduce the A 3 averaging operator This matrix turns out to have a two-dimensional null eigenspace generated by (1, 0, 0) t and (0, 1, 1) t , and a one-dimensional eigenspace with eigenvalue 1 generated by (0, −1, 2) t . We recover in particular that two of the a priori allowed invariants identically vanish: The A 3 -invariant space of I 2 interactions is therefore one-dimensional and generated by the -independent quantity (one can indeed check that (0, −1, 2) t is an eigenvector of F with eigenvalue 1): Summing over we find In summary, we conclude that: 1) there is no I 2 interaction invariant under the full permutation group S 3 ; 2) the space of I 2 interactions invariant under the subgroup of even permutations A 3 is one-dimensional, and generated by (24).
Remembering also that I P 0 = I V 0 = 0 for Majorana fermions, we write the A 3 -invariant interacting action for the Majorana case as

Large N : Schwinger-Dyson equation and mass generation
A powerful tool in the study of the large-N limit is provided by the Schwinger-Dyson (SD) equations. In zero dimensions they typically allow to solve tensor models and find their critical behavior. In one dimension, they play a crucial role for the SYK model and its tensorial generalizations, as they allow to identify a conformal IR regime and to solve the theory in that limit. In this section we will study the SD equations for our Dirac and Majorana twodimensional tensorial field theories, and in particular we will address the question of dynamical mass generation.

Dirac case
In two dimensions we do not expect any breaking of continuous global symmetries. This fact is known in the quantum field theory literature as Coleman theorem [38] and in the statistical mechanics literature as Mermin-Wagner theorem [37], and we will come back to it in section 4. Therefore, we should be allowed to assume U (N ) 3 -invariance of the theory, which implies that the propagator is proportional to the identity in tensor space: It is convenient to introduce a set of rescaled couplings defined as Then, the large-N SD equations are: where the factor 3 originates from the sum over colors in the I 1 interaction, see figure 2. Combining the two, in Fourier space we find Notice that this has the structureĜ −1 (p) = i / p − Σ 0 , where Σ 0 is a momentum-independent quantity (and a matrix in spinor space). Furthermore, the two-point function should be invariant under the rotation group (including parity transformations), so that Σ 0 cannot be proportional to γ µ or γ 5 . Therefore, Σ 0 must be proportional to the identity matrix, Σ 0 = −m1l; in the following we will omit writing explicitly the identity 1l unless we need to highlight its presence. This is consistent with (30) because γ µ and γ 5 are traceless, and it implies that the full prop- . Therefore, we can write a consistency equation for m, also known as the gap equation: Notice that the P and V interactions have dropped out of the equation. Since the integral on the right-hand side is UV divergent, we introduce a UV cut-off Λ, and obtain: By analogy with the GN model, we expect m = 0 to correspond to an unstable vacuum, while (33) should be the physical vacuum. We will verify this expectation in the following section, by studying the effective potential in the intermediate field formalism. The non-zero solution (33) is a non-perturbatively generated mass, which breaks the discrete chiral invariance of the bare theory. Requiring that this mass is invariant under a change of UV scale Λ yields the renormalization group equation: where β X i := Λ∂λ X i /∂Λ. Hence we must have: Moreover, by direct examination of the possible Feynman diagrams contributing to the two flows (see figure 3), one notices that β S 1 can only receive terms in (λ S 1 ) 2 , while β S 0 receives contributions in λ S 0 λ S 1 and (λ S 1 ) 2 . Hence we may conclude that the single beta functions are: Their associated flow is depicted in figure 4.
In conclusion, we have the same situation as in the standard GN model: asymptotic freedom (for both λ S 0 and λ S 1 ; we have not discussed the beta functions for the P and V couplings as they do not affect the 2-point function) and mass generation. Notice that (37) leads to (2) for λ S 1 = 0 and λ S 0 = λ, as expected.  (36) and (37). Arrows point towards the IR. The red line is the locus λ 0 + 3λ 1 = 0.
More interestingly, we see from (33) that we now have the opportunity to impose the massless condition: Indeed, this subspace of theories is preserved by the RG flow (35) (as seen also in figure 4), and the signs of the beta functions are consistent with asymptotic freedom: in such a way that where g > 0 is the value of 3λ S 1 at some reference scale Λ 0 . However, in view of equation (35), the massless condition (38) is IR unstable, and therefore it cannot be reached by a generic theory. Nevertheless, given that the condition λ S 0 + 3λ S 1 = 0 is RG invariant, one could define the theory to live in such a subspace. In this case we would seem to have a massless, although non-conformal theory: in fact there is still a non-trivial renormalization group running (defining a flow along the red curve in figure 4), hence the theory is not scale invariant. The only zero of the beta functions is at the origin, i.e. the free theory.
It is however premature to draw any conclusion on the existence of a massless theory. There are at least two possible pitfalls: the condition λ S 0 + 3λ S 1 = 0 could correspond to a trivial (non-interacting) theory in disguise; or it could correspond to a range of parameters for which the theory is unstable. In particular it is natural to worry about the second possibility, given that the massless condition forces the two couplings to have opposite sign. In order to check whether any such negative scenario is realized or not, we could look at the two-point function for the composite operatorψ a ψ a and check whether it decomposes trivially or whether it shows any tachyonic poles. An easier method to effectively do such computation is to introduce an intermediate field representation, as we will do in section 4.

Majorana case
In the case of real Majorana fermions we have a minor and a major change with respect to the complex Dirac case. The minor change is some combinatorial factors. The major change is the presence of I 2 in the action (25).
The large-N limit for real rank-3 tensors with O(N ) 3 invariance has been studied in [47], but here we have two main differences: the spacetime dependence, and the spinor indices. The latter lead to two types of contributions from the I 2 vertices: one without and one with a trace on spinorial indices, both depicted in figure 5. Defining the rescaled couplings as (in agreement with [47]) the large-N SD equations are: where Unlike the tadpole contributions proportional to λ S 0 and λ S 1 , the sunset ones proportional to λ 2 have a non-trivial dependence on the external momentum p, and therefore we no longer can assume a simple momentum dependence for the 2-point function. Assuming no parity breaking we can write the Fourier-transformed inverse 2-point function aŝ Locality imposes that A(p 2 ) and B(p 2 ) have no singularities for real p, henceĜ(0) −1 = B(0). Assuming also that A 2 (p 2 )p 2 +B 2 (p 2 ) is a real and monotonically increasing function of p, we see that we will avoid tachyonic poles in the 2-point function if B 2 (0) ≥ 0, and that for B(0) = 0 we have a massless theory. Unfortunately, solving the SD equations for A(p 2 ) and B(p 2 ) is beyond reach. In fact our goal is more modest, as we would like to understand if a conformal theory is a possible solution of the SD equations for some values of the couplings (i.e. if there is a non-trivial fixed point). To that end, we could in principle try to assume that the 2-point function takes the formĜ where we introduced two parameters b and α to be determined from the SD equation. Naïve dimensional analysis of the equation then would fix α = 1, i.e. as a consequence of the large-N SD equation a conformal theory necessarily has zero anomalous dimension. However, besides being UV divergent, in the massless case the loop integrals are also IR divergent, and we are forced to reintroduce a mass in order to regularize them. In practice we will use the following recipe to search for non-trivial fixed points: the leading contributions to the 4-point function have a colored structure as shown in figure 6: such diagrams scale as (1/N ) 2 × (1/N ) 3/2 × N = (1/N ) 5/2 , and are therefore suppressed by a factor 1/N with respect to the natural scaling of I 2 . Furthermore, the results of references [48] and [47] guarantee that such a suppression occurs at any order in the Feynman expansion. As a ℓ ℓ Figure 6: Structure of the leading contributions to the one-loop beta function of λ 2 : such diagrams are suppressed by a factor of 1/N with respect to the natural scaling of λ 2 .
consequence of the last observation, the requirement that λ 2 couplings have zero beta function will amount to imposing no wave-function renormalization. See appendix C for a derivation of the beta function of λ 2 at leading order in both 1/N and the coupling.
Let us now proceed to apply this recipe to the above SD equations. We start by plugging into equation (45) the ansatz forĜ(p) (46), further specifying A(p 2 ) = 1, B(p 2 ) = m and setting p = 0. After a lengthy but straightforward computation we obtain the following gap equation where J is the divergent integral encountered before, i.e.
As before, we still have a solution with m = 0. We would again expect it to correspond to an unstable vacuum, but verifying this expectation is more complicated in the present case because we do not have a useful intermediate field representation of I 2 . In any case, if tachyonic poles depend continuously on λ 2 , we would expect this vacuum to remain unstable, at least for small λ 2 . We also have two non-trivial solutions: where The solution (50) is non-zero whenever λ 2 = 0, which suggests that a non-perturbative mass is always generated in the A 3 -invariant model, irrespectively of the values of the coupling constants. This implies that such a theory can never flow towards a conformal fixed point. For λ S 0 +3λ S 1 > 0, the physically relevant solution is given by m − . In fact, in the limit λ 2 → 0 this yields back (33), On the other hand, m + diverges in the same limit, which suggests that the theory in this case does not admit a weak-coupling limit. One may wonder whether the conclusion of this analysis would still hold in a more general model without color symmetry. To settle this question, let us determine the gap equation associated to the most generic action with only I 2 interactions: where we work in the = 1 representation for definiteness. As we have determined in section 2, this covers all possible type-2 interactions, with or without invariance under permutation of the colors. The parametrization (53) has the computational advantage that all the Feynman diagrams it generates in the self-energy have a spinorial trace, as in the left panel of figure 5.
The self-energy now reads: Note that there is no term in λ V,1 2 λ S,1 2 or λ S,1 2 λ P,1 2 as these cancel by antisymmetry in q 1 ↔ q 2 . After a somewhat lengthy but straightforward calculation, one arrives at: where J was defined in (49), while I is an m-independent finite integral which is the subject of much literature (e.g. [49,50]), i.e.
Equation (48) is recovered for λ P,1 2 /2 = −λ V,1 2 = λ 2 and λ S,1 2 = 0, as it should. Interestingly, one realizes that the divergent term proportional to J 2 , which is the source of the spontaneous generation of mass by I 2 interactions, is cancelled in the gap equation whenever Following our recipe further, we should now impose that the condition (57), together with λ S 0 = λ S 1 = 0, is stable under radiative corrections. As explained before, since equation (57) is homogeneous of degree two in the couplings of type-2 interactions, it is stable under radiative corrections at leading order in 1/N , because the only possible source of flow of such couplings comes from the wave-function renormalization, which however is the same for all of them. Therefore, it only remains to look for additional conditions which ensure that no I S, 0 and I S, 1 interactions are generated by the RG flow at leading order in N , so that our SD equation with only type-2 interactions remains consistent. It is easy to check that type-0 interactions being topologically disconnected, they are never generated by the connected type-2 interactions. For the type-1 interactions, at first sight this looks doomed, since one would naively need to impose three conditions (one for each = 1, 2, 3, due to the color non-invariance) in a two-dimensional theory space (defined by (57)). As it turns out, we will see that these conditions are not independent and admit non-trivial solutions. The leading order 4-point functions of melonic theories are sums of ladder diagrams (see e.g. [48]). In the present theory, we have three types of elementary ladders, as shown in figure 7. The effective 4-point interactions generated by the ladders L are of type I X, 1 . Note that for = 1 one needs to use the Fierz identities to re-express the naturally generated interactions in our chosen basis.
In the following, we will assume thatĜ is of the form (46) with B(p 2 ) = 0. The kernel associated to the diagram L 3 is (we keep the conservation of external momentum implicit) Forgetting for a moment about IR divergences, at p 1 = p 2 , this takes the form while the latter is proportional to: In order to maintain λ S,3 1 = 0, the terms in δ IJ δ KL must cancel each other, yielding the condition: We can proceed similarly for the ladder L 2 , but by symmetry of the colors 2 and 3 in this problem one does not get any new condition: λ S,2 1 = 0 imposes no more than λ V,1 2 λ P,1 2 = λ V,1 2 λ S,1 2 . There is one extra condition to be derived from λ S,1 1 = 0. The relevant ladder is L 1 , with amplitude:  (154) and (155), the corrections to λ S,1 2 generated by these two terms compensate each other whenever which is the last consistency condition we can derive from the 4-point function.
In summary, we have learned from steps 2 and 3 of our recipe that we might obtain a massless theory provided one of the following two sets of conditions is satisfied: or (λ S,1 2 ) 2 = (λ P,1 2 ) 2 = 0 and λ V, Interestingly, such conditions are satisfied by models invariant under the continuous chiral symmetry (165).
Step 4 just amounts now to looking for possible conditions that avoid wave function renormalization. Unfortunately the previous steps left us with only one free coupling, thus either we are lucky and there is no need of wave function renormalization or such further requirement will lead us to a free theory. Assuming that one of the two set of conditions (65) and (66) hold, the free-energy Σ(p) evaluated at first order in the external momentum p is: which is of course both IR and UV divergent. In dimensional regularization one finds a simple pole divergence (see appendix C), or a logarithmic divergence in a momentum cutoff. This means that we have a non-trivial wave function renormalization, which can only be eliminated by the condition λ = 0, corresponding to a free theory. In conclusion, the conditions (65) and (66) do not correspond to a fixed point.
In appendix C we show that the coupling λ in equation (65) has, at leading order in 1/N and λ, the following beta function: As a consequence, unlike λ S 0 and λ S 1 , the coupling is not asymptotically free, it is IR free for both signs of the coupling. Such a result derives from the fact that at leading order in 1/N the beta function starts at two-loop order. The sign in (68) is thus compatible with results from the usual GN model, where the leading terms in 1/N of the two loop contributions to the beta functions have the same positive sign as our (68) (see for example [44]). The key difference is that for the GN model the two-loop part of the beta function is subleading in 1/N with respect to the one-loop part, while here the situation is reversed.
We conclude this section by noticing that in d = 2 − dimensions the beta function becomes thus showing a non-trivial IR fixed point of order √ , both at positive and negative coupling.
On the contrary, no UV fixed point is found in d = 2 + dimensions. The situation is thus reversed with respect to the GN model. Such result provides support for the conjecture in [39] that the theory in d = 2 − dimensions flows to a weakly interacting IR fixed point for small . 10

Intermediate field formalism for Dirac fermions
In the Dirac model, our study of the SD equation (30) suggests the existence of a subclass of theories in which the values of λ S 0 and λ S 1 conspire to give a vanishing non-perturbative mass (33). In the present section, we will investigate further the viability of the massless condition (38), by means of the intermediate field formalism (also known as Hubbard-Stratonovich transformation). As already apparent from the original GN literature [36], this formalism is particularly useful to elucidate the stability properties of the tentative solutions to the SD equations.
At first, it would seem that our tensorial version of GN necessitates the introduction of two types of intermediate fields: a) scalar fields associated to I X 0 interactions, as in the standard GN model; and b) matrix fields associated to I X 1 terms. However, it turns out that we can generate both I X 0 and I X 1 interactions with a single matrix intermediate field M X ij , provided we carefully adjust its free propagator.
In order to simplify the discussion we consider a reduced model, keeping only the interactions I S 0 and I S,3 1 in (11). This model is not color-invariant and it is not the most general one, but it 10 Although in appendix C we only discuss the models introduced in the present paper, it is straightforward to check that the calculation for the U (N ) × O(N ) × U (N ) model of [39] is essentially the same as for our model in (25). One finds the same anomalous dimension as in (185), with 3λ 2 2 replaced by the positive combination λ 2 = λ 2 1 + λ 2 2 − λ1λ2, with λ1 and λ2 being now the two couplings of the tetrahedral interactions of the model in [39]. Since the relation βi = 4ηλi for i = 1, 2 still applies in the large-N limit, both couplings are marginally irrelevant in two dimensions, and our conclusions for d = 2 − apply.
captures the essential qualitative properties of the invariant one: we have seen in the previous section that the P and V interactions do not affect the SD equations; furthermore, as we will see below, the difference between the color symmetric and the not color-symmetric case simply amounts to some factors of 3. Note that such a truncation is stable under the renormalization group flow at leading order in 1/N . In the language of [45], the 4-point melonic diagrams generated by I S,3 1 all have the same boundary graph: that of an I S,3 1 interaction. This ensures that no effective I S, 1 interaction with = 3 can be generated at leading order.
More explicitly, we consider in this section the reduced model: with where we have introduce the Hermitian matrix-valued bilinear We introduce the operator where is the projector on the trace part of an N × N matrix, and thus 1 − P is the projector on its traceless part. The eigenvalues of C are (1 − b) with multiplicity 1 (eigenvector proportional to the identity matrix), and 1 with multiplicity N 2 − 1 (with eigenvactors the basis of traceless Hermitian matrices). Thus C is always well defined, it is positive for b < 1 and it has a negative (zero) eigenvalue for b > 1 (respectively b = 1). It is thus invertible for b = 1, with inverse we can rewrite the interaction action as By the standard Hubbard-Stratonovich trick we can then rewrite the exponent of the interaction action as an integral over a Hermitian matrix-valued field M (the intermediate field); with obvious notation: Including also the fermions' free action, and making the notation completely transparent, the total Lagrangian density with the intermediate field reads This action is still invariant under discrete chiral transformations, which now act as (4), together with M → −M . As noticed above, for b = 1, C becomes degenerate, which manifests itself by the divergence of the coefficient in front of Tr[M ] 2 in the Lagrangian. 11 The condition b = 1 is equivalent to the massless condition derived in the previous section, and is therefore particularly relevant. It amounts to imposing the traceless conditions Tr[M ] = 0 in the path-integral, as can be seen by applying the Hubbard-Stratonovich trick once more, this time to the double trace. That is, by introducing a scalar φ and replacing We see that for b = 0 the field φ decouples from the rest of the action, while for b = 1 the quadratic part vanishes and the integral over φ gives a Dirac delta for the trace of M . We will use the representation with the field φ in appendix D, while for the remainder of this section we will stick to the double-trace representation.

Effective potential
The ψ fields can be integrated out from the partition function where and Tr includes a functional trace on top of the matrix trace.
To determine the vacuum configurations of the intermediate field, and hence check whether there is spontaneous symmetry breaking of the discrete chiral symmetry 12 , we want to compute the effective potential. By definition, the latter is the effective action (i.e. the generator of one-particle irreducible graphs) evaluated at constant field. In principle, in order to obtain the effective action we should introduce sources J coupled to M , perform the functional integral and lastly do a Legendre transform of W and we see that the action is of order N 3 (each trace being a sum over N elements). As a consequence, in the large-N limit we can bypass the integral over M and evaluate W [J] by saddle-point method, obtaining that it is simply given by the Legendre transform of S[M ]. 13 And since the Legendre transform is an involution, we find that the effective action in the large-N limit is simply given by S[M ] itself. Therefore, the effective potential we are looking for is obtained as From now on we focus on constant fields M . The effective potential in Fourier space is then: When expanding the logarithm in power series, the fact that the trace of an odd number of gamma matrices vanishes implies that only positive powers survive, so that we get: Including a UV regulator Λ as well as an IR regulator at intermediate steps since, even though their sum is not, the amplitudes are also IR divergent, we obtain: Since the effective potential has a global U (N ) invariance, it is convenient to express it in terms of the eigenvalues µ 1 , . . . , µ N of M . We obtain: where v(µ) := 1 2 is nothing but the effective potential of the standard GN model [36], up to a different normalization of the coupling constant.

Beta functions
The effective potential should satisfy as usual the renormalization group equation from which we obtain Using (76) we obtain also The beta functions so obtained for λ 0 and λ 1 are consistent with the ones obtained earlier for λ S 0 and λ S 1 in (37) and (36), up to extra factors of 3 due to the color permutations.

Stationary points and their stability
We can diagonalize them and rewrite where the variables µ i , for i = 1, . . . , N , are the eigenvalues of M . We seek a unitary invariant solution of the above equations of motion. The only U (N )invariant matrix is a matrix proportional to the identity, hence we look for a solution such that µ i = µ for every i. Equation (94) reduces to which admits the solutions µ = 0 and The latter breaks the discrete chiral symmetry, while the former is of course invariant. In order to investigate the stability of the solutions we need to look for the eigenvalues of the Hessian matrix ∂ 2 V ∂M ij ∂M kl . The Hessian is difficult if not impossible to write at a general field configuration, since in general M does not commute with its infinitesimal variation dM . 14 We didn't have this problem when writing ∂V ∂M ij because of the trace in the potential, but once we take one derivative the trace is gone, and for the second derivative we encounter this problem. However, the symmetric solution is proportional to the identity, so that it commutes with any matrix, and we can easily write the Hessian evaluated at the symmetric solution. It reads where and the operator P was defined in (74). Therefore, the Hessian has eigenvalues α(µ) with multiplicity N 2 − 1, and β(µ) with multiplicity 1. We can associate the latter with the variation of µ, since its eigenvector is proportional to the identity, while α(µ) is associated to infinitesimal SU (N ) transformations, its eigenvectors being given by the basis of traceless Hermitian matrices. When µ → 0 both α(µ) and β(µ) go to −∞, hence we conclude that the solution µ = 0 is always unstable, just as in the GN model.
For the non-zero solutionμ in equation (96) we find instead In this case, all the eigenvalues are positive for b < b c ≡ 2λ 1 π+2λ 1 , hence the solution is stable in this range. Therefore, we have spontaneous breaking of the discrete chiral symmetry, with the generation of a mass m for the fermions. The latter can be read off from (79) once we replace M ij → Nμδ ij : in agreement with (33), again up to a factor 3 due to the color permutations. However, since the N 2 − 1 degenerate eigenvalues α(μ) become negative for b > b c , we conclude that this vacuum becomes unstable before reaching the interesting value b = 1, corresponding to λ 0 + λ 1 = 0. We might be tempted to interpret the appearance of zero-modes at b = b c as a signal of a second-order phase transition, but there are two problems with such an interpretation. First, as we will see, at b = 0 many disconnected global minima are present, and the symmetric solution with µ i =μ ceases to be the global minimum at b > 0. Therefore, we have a first-order phase transition at b = 0. Second, and most important, any other minima break the U (N ) invariance of the model, but spontaneous symmetry breaking of a continuous symmetry is forbidden in two dimensions in virtue of the Coleman-Mermin-Wagner theorem. Therefore, we need further investigation in order to understand if the model at b > 0 is consistent, and in case it is, whether after restoration of the U (N ) invariance it corresponds to a different phase (e.g. with no breaking of the discrete chiral symmetry) or not. We will come back to such question in section 4.5.
The existence of other solutions is clear at b = 0: the equations of motion (94) decouple from each other, and reduce to Therefore, besides the usual zero solution, we have where the plus or minus sign can be chosen independently for each i = 1, . . . , N . And since at b = 0 the potential depends only on µ 2 i , for each choice of signs we obtain the same value of the potential. Hence we have N + 1 minima of the potential, one for each choice of signs (up to permutations), two of which are the b = 0 limit of (96), while all the others break the U (N ) invariance down to U (N + ) × U (N − N + ) (where N + is the number of plus signs in the given solution).
As we show in appendix D, for b = 0 all the N + 1 solutions above generalize to as many stationary points. However, the degeneracy gets lifted, and while for b < 0 the symmetric solution (96) is the global minimum, for b > 0 it actually becomes the stationary point with the largest value of the potential. The global minimum at b > 0 turns out to be the solution which maximizes the symmetry breaking, i.e. the one with an equal number of plus and minus signs (we assume from here on an even N ). 15 Such a solution exists at any b, and it is easily obtained by demanding the matrix M to be traceless. Then equations (94) reduce again to (102), with solution (103), plus the traceless constraint enforcing the number of plus and minus signs to be the same.
Evaluating the potential at the symmetric solution (equation (96)) we get while at the maximally broken solution (equation (103) At b = 0 we have v 1 = v 2 as discussed above, while v 1 < v 2 for b < 0 and v 1 > v 2 for b > 0, confirming that in the latter case the symmetric solution is no longer the global minimum.
To conclude the stability analysis, we can explicitly check that the maximally U (N )-breaking solution leads to the presence of 2N + (N − N + ) = N 2 /2 (would-be) Goldstone modes. Since such solution has µ 2 i =μ 2 independent from i, we know that at this configuration M 2 (but not 15 When N is odd, we can pick up N −1 2 positive (resp. negative) eigenvalues and N +1 2 negative (resp. positive) eigenvalues. The difference with respect to the N even case is negligible in the large-N limit.
M itself) commutes with dM , so we can expand the ln M 2 that appears in the first derivative 16 and obtain: whereα kl (μ) = We see thatα kl (μ) = 0 for k ≤ N/2 and l > N/2, as well as for l ≤ N/2 and k > N/2; hence we have N 2 /2 massless modes. These look like the Goldstone modes associated to the spontaneous symmetry breaking of U (N ) down to U (N/2) × U (N/2). However, spontaneous symmetry breaking of a continuous symmetry is forbidden in two dimensions [38,37], therefore we will refer to them as "would-be Goldstone modes", and we will discuss their role in section 4.5.

2-point function
The effective action for the intermediate field generates by definition the one-particle irreducible (1PI) n-point functions of the intermediate field. And since the latter is conjugated to a fermionic bilinear, one can say that the effective action for the intermediate field generates the 1PI n-point functions for that particular class of composite fields (see [36] for the precise relation in the GN model case). Here we compute the 2-point function first in the symmetric vacuum, and then in the maximally-broken one. For the former, we show that the instability at b > b c manifests itself in the presence of tachyonic poles. In the second case, we show that at small momentum p the behavior of the 2-point function is the standard 1/p 2 , corresponding to a logarithm in position space, the consequence of which we will discuss in the following subsection. The computation of the 2-point function is most easily performed using the action with the fermion fields not yet integrated out, i.e. from the Lagrangian (79), expanded either around the stationary point (96) or (103) with N + = N/2.
We define the intermediate field connected 2-point function (or propagator) as where the average is taken on the vacuum state. Its inverse, i.e. the second variation of the effective action at the field configuration corresponding to the vacuum state, satisfies the usual SD equation where Π is the self-energy. The latter is given in the large-N limit by a single diagram, with a single fermionic loop, shown in figure 8.
In the symmetric vacuum, the self-energy in momentum space is given by where we have used the massive fermionc propagator, with mass (101). The integral I(p) is UV divergent and thus requires a UV cutoff Λ. After a standard computation we obtain where we defined as in [36] B(p, m) ≡ Imposing the renormalization condition we obtain at last where Since B(p, m) is a monotonically increasing function of p, such that B(0, m) = 2 and at large momentum B(p, m) = ln p 2 m 2 + O( 1 p 2 ), we see that as soon as 2λ 1 π < b 1−b (i.e. for b > b c ) the function f (p) −1 has a pole at p > 0. Since we are in Euclidean signature, the latter is a tachyon.
Let us now consider the connected two-point function on the maximally broken vacuum. The calculation goes as before, except for the renormalization condition at zero momentum, which this time should reproduce the Hessian around the symmetry-breaking vacuum, and for the fact that we need to keep track of the different signs in the mass terms.

Effective non-linear sigma model and Coleman-Mermin-Wagner theorem
The intermediate field analysis of the Dirac model shows that the symmetric solution anticipated from the SD equation in section 3 is unstable in the limit b → 1 (i.e. λ 0 + λ 1 = 0). The true vacuum configuration of the intermediate field generates a non-zero mass and thereby breaks the discrete chiral symmetry of the bare action. More surprisingly, it also seems to break the continuous U (N ) symmetry of the model, which would be in obvious tension with the Coleman-Mermin-Wagner theorem [37,38]. This is reminiscent of a similar conundrum arising in the chiral version of the GN model [36]: in this theory, the discrete chiral symmetry is enhanced to a global U (1) symmetry, which looks naively broken by the spontaneously generated mass. However, it was shown by Witten that this impression follows from the incorrect assumption that the chiral angle admits a welldefined expectation value in the broken vacuum [51] (see also [52] for a recent discussion of this problem in the context of the AdS/CFT correspondence). Similarly, we will show that the U (N ) angles parametrizing the broken solutions of our model have a non-trivial dynamics governed by an effective non-linear sigma model, which suggests that a similar mechanism of symmetry restoration as in the chiral GN model is at play. In view of the involved dynamics of the massless modes we will uncover, we will not be able to provide a completely explicit check of symmetry restoration in our theory, but analogies with the chiral GN model will allow us to conclude with confidence that it does indeed occur.
Before moving to the specifics of our model, let us describe the main mechanism behind symmetry restoration in two dimensions. Consider a theory with a continuous global symmetry, and a set of covariant operators O A whose expectation values are zero by symmetry, but such that some linear combination of their 2-point functions On the contrary, if the limit is non-zero we have spontaneous symmetry breaking (with O A = ρ A = 0). In such case, by Goldstone's theorem there must be massless modes corresponding to the broken symmetries. In the deep IR we expect only the latter to be relevant, so it makes sense to freeze all the other (massive) modes and only look at the massless fluctuations around the vacuum. In dimension d > 2 these lead to consistently with the symmetry breaking picture (with the constant piece coming from the disconnected contribution). However, in d = 2 we obtain and since the logarithm diverges at large distances the expansion breaks down. What this means is that the massless modes fluctuate wildly at large distances and they cannot be treated perturbatively. Their nonperturbative treatment will generally lead to symmetry restoration (for example this happens if in (123) one has actually expanded |x| −α = e −α ln |x| to first order in α, which of course gives a wrong result at large |x|). We would like to explicitly check this scenario for our model at b > 0. We have seen in the previous subsection that indeed in the broken vacuum the would-be Goldstone modes lead to a logarithmic divergence for the (connected) 2-point functions at large x. However, this comes from having chosen one (simple) representative in the space of degenerate vacua and having expanded the unitary matrices around the identity. We want to keep them nonperturbative (or in other words average over the vacua) and compute whereM ab = sgn((N + 1)/2 − a)μδ ab , and θ(x) stands for the (x-dependent) N 2 − N angles parametrizing U (N )/U (1) ⊗N . In the second line we have frozen the N modes corresponding to shifts in the eigenvalues, which are massive and thus irrelevant in the deep IR. Furthermore, sinceM ab is invariant under the subgroup U (N/2) × U (N/2) (the stabilizer), we only need to consider unitary matrices in the coset space Gr(N, N/2) ≡ U (N )/(U (N/2) × U (N/2)), also known as the complex Grassmannian space. By symmetry arguments, the low energy effective action for the would-be Goldstone modes V ∈ Gr(N, N/2) must be given by a complex Grassmannian non-linear sigma model Let us investigate this point in more detail and determine the proportionality factor in equation (125). As previously argued, in order to capture the IR massless modes, we may restrict the pathintegral over M to the domain: Following the literature (see e.g. [53,54,55]), it is convenient to parametrize F in terms of orthogonal projectors P of fixed rank rk(P ) = N/2. This is simply given by the change of variables To explicitly check the nature of P , it suffices to note that When U is varied over U (N ), P spans the set of all orthogonal projectors of rank N/2. Hence the integration over the fields M (x) can be replaced by an integration over operators P (x) verifying: P 2 = P = P † and TrP = N/2. Moreover, the Jacobian of this transformation is independent of x and can therefore be ignored. Consider a fixed configuration P (x) = U (x)P U † (x) of the projector field. There clearly is some freedom in our choice of matrix field U : we can multiply it on the right by unitaries which commute withP , and on the left by unitaries which commute with P . By differentiation of the relation defining P , one finds that: The matrix P being Hermitian, the operatorP : A → [P, A] is itself Hermitian on the space of N × N matrices equipped with the inner product A, B = Tr(A † B). It has eigenvalues 0, 2 and −2, with multiplicities N 2 /2, N 2 /4 and N 2 /4 respectively. Hence the Maurer-Cartan form dÃ := U dU † decomposes uniquely into dÃ = dÃ 0 + dÃ + + dÃ − , where dÃ 0 is the orthogonal projection on the null eigenspace ofM and dÃ ± are the projections on the eigenspaces with eigenvalues ±2μ. Since dÃ 0 does not contribute to the flow of P (129), one can fix the gauge (up to a global transformation) by the condition: More explicitly, the null eigenspace ofP is immediately seen to be E 0 (P ) = P M N (C)P ⊕ (1 − P )M N (C)(1 − P ), with M N (C) being the space of complex N × N matrices, and the gauge-fixing condition is therefore equivalent to: In turn, these two equations are equivalent to: We therefore conclude that, in the chosen gauge, the Maurer-Cartan form dA := U † dU takes the form for some B. Another useful relation is obtained by conjugation of the flow equation (129): We can now proceed with the computation of the effective kinetic term of P , starting from the Lagrangian: We express P (x) = U (x)P U † (x) with U dU † ∈ E 0 (P ), and change variables in the path-integral over the fermions: ψ(x) → U (x)ψ(x) andψ(x) →ψ(x)U † (x). The path-integral measure is invariant under this local unitary transformation, but the kinetic term is not and generates an effective coupling between the Maurer-Cartan form and the the fermions: Importantly, our construction ensures that ∂ µ A := U † ∂ µ U is of the block-diagonal form (133).
Integrating out the fermions we generate as before an effective action, this time for the Maurer-Cartan form. Since ∂ µ A ij has dimension one, and since no terms in A ij without derivatives can be generated (because of the original global symmetry of the action), we are interested just in the part of the action which is quadratic in ∂ µ A ij (with no additional derivatives), higher order terms being irrelevant in the IR (and the linear part being zero). This is given by the same type of diagrams in figure 8 we encountered before for the self-energy of M (in the maximally broken vacuum). That is, we have and the fact that any trace of an odd number of gamma matrices is zero, we obtain Π µν ij;kl (0) = −N 2 1 2π Sincem is Λ-independent, we can take the limit Λ → +∞ at fixedm, and obtain Π µν ij;kl (0) = −N 2 1 2π sgn(µ k )sgn(µ l )δ µν δ ik δ jl .
Remembering that our gauge-fixing condition ensures a block-diagonal form of dA, see (133), we obtain: We can finally use (133) and (134) to conclude that: The key observation to be made is that the explicit dependence in our choice of unitary matrix U drops out, as it should. Remark also that (143) is mapped to (125) after identifying V (x) = 2P (x) − 1.
In order to explicitly check the restoration of the U (N ) symmetry at the level of the two-point function (124), we would need to determine the large |x| asymptotics of: Notice the N 2 factor in front of the action (143). As discussed before, this marks a difference with usual matrix-valued models, which would have just a factor N , and for which the Vandermonde determinant would thus play an important role. In our case, the large-N limit at leading order tells us simply that ∂ µ V = 0, i.e. V is constant. In such case, in evaluating (144) at leading order we find K(x) = Nμ 2 , and therefore we would seem to have spontaneous symmetry breaking of the U (N ) invariance. In the case of the apparent symmetry breaking of continuous chiral symmetry in the GN model, Witten determined that the two-point function of interest goes like |x| −1/N . This observation explains why one may find a non-zero constant in the limit N → ∞, while at the same time obeying the Coleman-Mermin-Wagner theorem. In particular, the U (1) chiral symmetry is preserved, and though there is an infrared massless mode in the theory (the so-called θ particle), it is not a Goldstone boson. By analogy, in the present model we might expect to find a power-law decay at finite N and large |x|, K(x) ∼ |x| −c/N α for some constant c and a positive exponent α, corresponding to quasi-long-range order. However, to the best of our knowledge there are no examples of this behavior for non-abelian symmetries, and on the basis of the generic features of non-linear sigma models on homogenous spaces, such as asymptotic freedom and dynamical mass generation, we expect to find more likely an exponential decay, K(x) ∼ e −m|x|/N α , corresponding to disorder. 17 Checking this explicitly in the complex Grassmannian sigma model (143) is however technically challenging, and certainly goes beyond the scope of the paper. We are not aware of any completely conclusive calculation in this respect, but we refer the reader to [56] for an early attempt based on a renormalization group analysis.

Conclusion and outlook
We have studied the large-N limit of a new class of fermionic models in two dimensions with quartic interaction and a tensorial index structure. In many respects the Dirac models with U (N ) 3 symmetry turn out to be still very similar to the usual Gross-Neveu model: they are asymptotically free and they have a non-perturbative mass gap. For such models we have introduced an intermediate matrix field representation, by means of which we have discovered a first order phase transition to a phase with apparent breaking of the U (N ) 3 symmetry, and with an effective action for the would-be Goldstone modes given by a complex Grassmannian non-linear sigma model. The U (N ) 3 symmetry is non-perturbatively restored in two dimensions, as we discussed in section 4.5, but it would be interesting to explore such new phase in higher dimensions.
The Majorana models with O(N ) 3 symmetry represent a major departure from GN-type physics. They are in fact the models which most closely resemble the SYK model, satisfying for example similar large-N Schwinger-Dyson equations; the latter imply in particular that their self-energy has a non-trivial momentum dependence, thus signaling a crucial departure from the standard GN model. However, they have a crucial difference also with respect to the SYK model: in two dimensions the coupling is marginal and the free-propagator term in the Schwinger-Dyson equation cannot be discarded neither in the IR nor in the UV. It turns out that at leading order in 1/N the coupling of their tetrahedral interaction is marginally irrelevant, i.e. the free theory is UV unstable. As a consequence, the theory develops a weakly interacting IR fixed point in d = 2 − dimensions for small , as recently conjectured in [39]. In precisely two dimensions we found no non-trivial conformal field theories within our space of theories.
One of the most appealing features of the Gross-Neveu models, together with their large-N properties, is that they are integrable, see e.g. [57,58]. This means that, as classical field theories, they possess an infinite number of conserved charges which constrain their dynamics [59,60]. Remarkably, integrability carries over to the quantum level, taking the form of Yangian symmetry, see for example [61]; this allows to fix the quantum scattering matrix of the theory exactly, even at finite N . From this point of view, tensorial Gross-Neveu model are an ideal playground for exploring integrability for 1 + 1-dimensional tensor models. In fact, for the simplest interaction which we considered, I S 0 , our model simply gives several copies of the usual GN model, and integrability should follow automatically, albeit somewhat trivially. Exploring more general interactions, however, might yield integrable deformations of the underlying Yangian algebra. It would be interesting to investigate this possibility, and in particular to elucidate whether there is any relation between integrability and color symmetry.
Another possible avenue for further applications of tensor models are large N gauge theories. In d = 1 models, gauging the G r symmetry is straightforward as it essentially amounts to restricting the set of observables to the singlet sector [22]. In contrast, the gauge potential acquires a non-trivial dynamics in higher dimensions, which might allow to construct interesting large-N tensorial gauge theories.
In two spacetime dimensions the Clifford algebra {γ µ , γ ν } = 2δ µν admits two-dimensional representations. In fact, it suffices to notice that the Pauli matrices σ x , σ y , and σ z satisfy the three-dimensional Euclidean Clifford algebra, hence one can select two of them as γ µ , and the third one as γ 5 (the subscript "5" really makes sense only in four dimensions, but we follow the very common convention of using the same notation also in two dimensions). We will mostly work in the following Majorana representation: In the case of Dirac fermions the following Weyl representation can also be useful: For both representations we have the identities: Furthermore, both γ µ and γ 5 are Hermitian matrices and they square to the identity matrix. We denote the spinorial indices by capital letters from the middle of the latin alphabet, e.g. I, J, K, L = 1, 2. We have the following useful relation: and the completeness relation By contracting the latter with either two γ 5 or two γ µ we find also For what concerns us, Majorana (Dirac) fermions ψ I a are real 18 (complex) Grassmann fields, with a spinorial index I labelling the two spinorial components. The boldface index a denotes a collective flavor index associated to a representation of the symmetry group G of the model (for example a ≡ i = 1, . . . , N in the O(N )-invariant GN model, or a ≡ abc, with a, b, c = 1, . . . , N for the O(N ) 3 -invariant tensorial version). We will usually omit one or the other type of indices, unless they are necessary.
The definition ofψ is a bit tricky in Euclidean signature, and it differs significantly in the literature. One option, followed by many, is to double the number of spinor fields in Euclidean space, with respect to Minkowski space, by treatingψ as an independent field not related to ψ by Hermitian conjugation. One negative aspect of such construction is that the resulting action is not Hermitian. This might not be a problem for fermions in general, but as we want to preserve the hermiticity property of the matrix-valued bilinear in equation (72) (true in Lorentzian signature, and needed for our intermediate field analysis), we follow here an alternative route, due mainly to Mehta [62] and van Nieuwenhuizen and Waldron [63], and which for our purposes boils down to the following observation. We can defineψ = ψ † Γ 0 for some matrix Γ 0 to be determined. Requiringψψ andψγ µ ψ to behave as a scalar and a vector, respectively, under rotations allows us to partially fix Γ 0 . In order to see that, we define the spinorial representation of the generator of rotations as It satisfies where is the vectorial representation of the rotation group, i.e. a rotation is written as where ω µν = α µν . Defining a finite rotation in spinorial representation as one can check that Note that in the Majorana representation is real (thus preserving the reality condition) but non-diagonal, while in the Weyl representation is diagonal but complex (as it should be, since there are no real one-dimensional representations of SO (2)). Therefore, forψψ andψγ µ ψ to transform appropriately we just need that if ψ → S(R)ψ, thenψ →ψS † (R). This is trivially achieved if Γ 0 is either the identity or γ 5 (or a linear combination of the two), because they both commute with S † (R). Either choice can be found in the literature, however, the requirement of a consistent Wick rotation from Lorentzian to Euclidean signature selects Γ 0 = γ 5 (see [62,63] 19 ). We can also define a continuous chiral transformation as with real θ. Under such a transformationψγ µ ψ is invariant whileψψ andψγ 5 ψ are not (here we are usingψ = ψ † γ 5 →ψe θγ 5 ). Furthermore, the combination (ψψ) 2 − (ψγ 5 ψ) 2 is invariant. However, for no real value of θ does such a continuous chiral symmetry reduces to the discrete chiral symmetry (4), which is instead obtained for θ = i π/2 (where one should remember that even for complex θ the transformation forψ isψ →ψe θγ 5 , and notψ →ψe θ * γ 5 ).

B Bilinears, quadrilinears, and symmetries
We denote a bilinear in the fermionic fields by means of parenthesis indicating full contraction of the spinorial indices, e.g.
where Γ stands for a generic product of γ-matrices. Euclidean space symmetry: In two dimensions we have only 4 independent bilinears, as opposed to 16 in four dimensions: which behave under the rotation group as scalars, vectors, and pseudoscalars, respectively. Besides S ab , we can form another scalar by inserting a derivative in V µ ab : where as usual we defined / ∂ µ = γ µ ∂ µ . We have only three independent scalars which we can construct out of four fields and with no derivatives: I P abcd = P ab P cd .
Quadrilinears with derivatives and higher order multilinears correspond to non-renormalizable interactions and are therefore not of interest for our purposes. Under the continuous chiral transformations (165), it's again (170) the only invariant bilinear, while among the quadrilinears the invariants are I V abcd and the combination I S abcd − I P abcd . Flavor symmetry: In order to complete the construction of the Lagrangian we need to contract the group indices in such a way to build invariants.
In the GN case, i.e. for G = O(N ) or G = U (N ), there are two possible ways to contract the four indices: one which we will call "disconnected contraction", a = b and c = d; and one which we will call "connected contraction", a = d and b = c. The reason for the terminology should be evident from figure 9. However, only one of the two types of contractions is independent. In fact, we can always use the completeness relation (154) to get rid for example of the connected contraction, in favor of a combination of invariants with disconnected contractions [42]. For Majorana fermions V µ aa = P aa = 0, hence we have simply (ψ a ψ b )(ψ b ψ a ) = − 1 2 (ψ a ψ a )(ψ b ψ b ), and thus one unique interaction. For Dirac fermions we have three independent interactions: I S aacc , I V aacc , and I P aacc . In the tensor case it is well known that there are many ways to construct invariants. 20 The invariants relevant to our case are discussed in detail in section 2. C Beta function of the λ 2 coupling at leading order In section 3.1 we have obtained the exact large-N beta functions for λ S 0 and λ S 1 by means of a Callan-Symanzik equation for the physical mass, and by a simple analysis of the possible diagrams at large N . The same method was not applicable to the λ 2 coupling in section 3.2, because it was not possible to solve the SD equations. In this appendix we are going to compute the beta functions for the coupling λ 2 of (25) and the coupling λ of (53) with conditions (65), and show that they are both IR free.
Let us start with the action (25). Note that, while λ 0 only receives quantum corrections from λ X 1 , there is no way to protect the latter from radiative corrections generated by λ 2 . However, the running of the latter is independent of λ 0 and λ X 1 at leading order in 1/N . In fact, as we explained in figure 6, at leading order in 1/N the vertex I 2 receives no quantum corrections at all. Its running can therefore only come from the wave-function renormalization Z, but since I 0 20 See for example [66] for an enumeration of invariants with n tensors. and I X 1 only contribute to the 2-point function with momentum independent tadpoles, Z must only depend on λ 2 . So for the purpose of calculating the running of λ 2 we can forget the other couplings.
Since in the presence of I 2 the 2-point function needs a multiplicative renormalization, we need to introduce a wave-function renormalization, Rewriting the action in terms of the renormalized fields, the effective coupling is Since at leading order in 1/N the part of the 4-point function with I 2 structure is exact at tree level, then λ 2,R is the renormalized coupling. Hence, its beta function is types of vertices). That is, we have Σ(p) = d d x e − i p· x Σ(x) = − i 12λ 2 c 3 2 π 8 1 + 2 ln(µ/p) γ µ δ νρ (δ µν p ρ + δ µρ p ν + δ νρ p µ ) + O( 0 ) = − i 6λ 2 c 3 2 π 1 + 2 ln(µ/p) / p + O( 0 ) , leading again to

D The other stationary points
We solve here the equations of motion for the Dirac model of section 4 in full generality. To that end, it is convenient to use the representation without the double trace, i.e. with potential The equations of motion read The latter is simply a constraint relating the new variable φ with the trace of M . The first equation can be simplified by a rescaling of the eigenvalues, defining We obtain which is solved (almost by definition) by where W (z) is the Lambert-W function (or product logarithm). Denoting by n + = N + /N the fraction of eigenvalues with the plus sign, and plugging the solution of (190) into (191), we arrive at a self-consistency equation determining z as a function of n + : Notice that the real branch of the Lambert-W function is defined for z ≥ −1/e, and since for 0 < n + < 1 we have both signs in its argument for the solutions above, we have the constraint |z| ≤ 1/e.
For n + = 1 we find (using e W (z) = z/W (z)) in agreement with the symmetric solution (96). However, in the range |z| ≤ 1/e, such equation can only be satisfied for hence the symmetric solution becomes disconnected from this family of solutions outside of this range of b (notice that for λ 1 > π/(2W (1/e)) the lower bound is replaced by b > −∞). Solving equation (195) for z is difficult, it is easier instead to solve it for n + . Since W (0) = 0, for the special case z = 0 (for which µ i reduces to (103)), we find n + = 1/2. That is, we find the traceless solution discussed in section 4.3 as a special case of the more general class of solutions. For z = 0, we obtain The computation of the Hessian for the general case is frustrated by the non-commutativity of both M and M 2 with the infinitesimal variation dM . We can however check that the value of the potential at these other solutions (with the constraint 0 ≤ n + ≤ 1), always lies between v 1 and v 2 , equation (104) and (105), respectively. Therefore, whether they are local minima, maxima or saddles, they do not affect the identification of the global minimum. Lastly, at b = 0 they all have the same value of the potential, i.e. they are the N + 1 degenerate solutions discussed below equation (103).