Conformal Symmetry and Composite Operators in the $O(N)^3$ Tensor Field Theory

We continue the study of the bosonic $O(N)^3$ model with quartic interactions and long-range propagator. The symmetry group allows for three distinct invariant $\phi^4$ composite operators, known as tetrahedron, pillow and double-trace. As shown in arXiv:1903.03578 and arXiv:1909.07767, the tetrahedron operator is exactly marginal in the large-$N$ limit, and for a purely imaginary tetrahedron coupling a line of real infrared fixed points (parametrized by the absolute value of the tetrahedron coupling) is found for the other two couplings. These fixed points have real critical exponents and a real spectrum of bilinear operators, satisfying unitarity constraints. This raises the question whether at large-$N$ the model is unitary, despite the tetrahedron coupling being imaginary. In this paper, we first rederive the above results by a different regularization and renormalization scheme. We then discuss the operator mixing for composite operators, and we give a perturbative proof of conformal invariance of the model at the infrared fixed points, by adapting a similar proof from the long-range Ising model. At last, we compute the two- and three-point functions of $\phi^4$ and $\phi^2$ composite operators. In the correlator of one tetrahedron and two quadratic operators, we find a logarithmic dependence on the external positions that is precisely of the type encountered in logarithmic conformal field theories, which implies that the model is actually non-unitary. We highlight how this and other aspects of the model show similarities to the bi-scalar conformal fishnet theory of G\"urdogan and Kazakov.


Introduction
Finding and studying non-supersymmetric interacting conformal field theories (CFT) in d > 2 dimensions is a very challenging task, with a long history. A young and successful approach is the conformal bootstrap [3], working completely within the framework of CFT and seeking to identify CFT's by imposing consistency conditions. The most standard and historical approach is however based on the renormalization group: here the challenge is to find non-trivial (interacting) fixed points. They correspond by construction to scale-invariant theories, but very commonly invariance under the full conformal group arises as well [4]. The great hurdle faced by the renormalization group is that our main tool of investigation, perturbation theory, demands that we approach the interacting fixed point by perturbing the free theory. Clearly, the task becomes more daunting the further the two theories are, and therefore it is extremely important from the theoretical point of view to have adjustable parameters that allow us to bring the interacting fixed point closer to the free theory or to tame the perturbative series. There are mainly two widely exploited parameters of this sort [5]: one is the analytically continued spacetime dimension, as typically the interacting fixed point collapses into the non-interacting one at some critical dimension; the other is the number of field components N , when their interactions are constrained by a (global or local) symmetry group, as the large-N limit can lead to drastic simplifications of the perturbative expansion.
Tensor models, in which N = N r with r the rank of the tensor, are a recent entry in the menu of field theories admitting an interesting large-N limit. They typically admit a large N limit dominated by melonic diagrams [6][7][8][9][10] different from both the large N limit of vector models (r = 1, dominated by bubble diagrams [11]) and the one of matrix models (r = 2, dominated by planar diagrams [12]). The melonic dominance translates into a perturbative expansion which is richer than the one of vectors and more manageable than that of matrices. It is therefore interesting to construct models of tensor field theories in d dimensions and look for their fixed points at large N [1,2,[13][14][15][16][17][18]. 1 We call the resulting conformal field theories melonic.
The main result of Ref. [1] is that in the large-N limit the tetrahedron coupling g is exactly marginal (its beta function is identically zero) while the beta functions of the other two couplings are quadratic in the two couplings themselves, with g-dependent coefficients. The model has four g-dependent fixed points, which are real for g purely imaginary and below some critical value. One of them is IR attractive for both the pillow and the double-trace couplings. Notice that from the perspective of the long-range Ising model (N = 1, only one quartic interaction) this infrared fixed point is surprising. In fact, our kinetic term corresponds to the transition point between the long-range behavior and the mean field theory one, in which the infrared fixed point disappears in the Ising case. The existence of a non-trivial IR fixed point in our O(N ) 3 model is exclusively due to thhe tensor structure and the large-N limit. The critical exponents at the IR fixed point, that is, the scaling dimensions of the pillow and double-trace operators, are real and above the unitarity bounds. The spectrum of dimensions of bilinear operators with arbitrary spin, as well as their OPE coefficients with two fundamental fields, has been computed in Ref. [2], where it was found again to be real and above the unitarity bounds.
Given the fact that the model has an imaginary coupling, it is to be expected that it is non-unitary, but the results of [1,2] raise the tantalizing possibility that in the large-N limit we could find a unitary theory. 2 For example, non-unitarity could manifest itself in some dimensions or OPE coefficients having imaginary parts which are suppressed in 1/N .
In this paper we address the following two questions: (i) does the large-N infrared fixed point of Ref.
[1] define a conformal field theory? (ii) Does it define a unitary theory?
In order to tackle the first question, we will adapt to our model the methods of Ref. [52], which gives a proof of conformal invariance to all orders in perturbation theory for the infrared fixed-point of the long-range Ising model with propagator C(p) = 1/p (d+ )/2 . Most of that proof is built on standard ideas (e.g. from Ref. [56]), except that the non-local propagator of the long-range model implies the absence of a local energy-momentum tensor. The main point is to use the Caffarelli-Silvestre trick [57] of localizing the kinetic term by means of an embedding of the theory in d + p dimensions, with p = 2 − (d + )/2.
The main differences in our case are that: first we are interested in = 0, and second we must deal with multiple quartic interactions which mix under renormalization. We thus need to revisit the results of Ref. [1] and [2] with an analytic regularization (rather than using a momentum cutoff), and construct renormalized composite operators. Once this is done we conclude along the lines of [52] that the fixed-point theory of our O(N ) 3 model is indeed conformally invariant.
We then compute two and three-point functions among the renormalized composite operators, which we use to further test conformal invariance and to address our second question. Conformal invariance greatly constraints these correlators [58]: the two-point functions between operators of different scaling dimension are zero, the three-point functions among three scalar operators are completely fixed by their dimensions (up to an overall OPE coefficient), and so on. Such constraints are respected by all the correlators we have computed, with a caveat in one case on which we comment further below.
For special values of the dimensions of operators some subtle issues appear [59]. For instance one should be careful to distinguish between the d-dimensional Dirac delta δ(x), which is a homogeneous distribution under conformal transformations, and |x| −d which has a singular Fourier transform and upon regularization does not transform homogeneously [60]. We need to deal with this issue as in our model we have an exactly marginal operator, the tetrahedron, which by definition has dimension d. We find contact terms (i.e. terms including a delta function) in its three-point function with itself Eq. (6.25), as expected [61,62]. These contact terms do not lead to an anomaly as their coefficients are finite. However, an anomaly can arise from its two-point function [63], which has the functional form 1/|x| 2d , and which in two dimensions has a singular distributional limit: 1/|x| 2d− ∼ 1 (∂ 2 )δ(x). We therefore expect a conformal anomaly in d = 2. Such a conformal anomaly in two dimensions would not be a big surprise, but in the absence of a local energy-momentum tensor it is not obvious how it should be interpreted. It should also be noted that in Ref. [2] a puzzling discontinuity was found in the spectrum of bilinear operators: the computation at d = 2 differs from that at d = 2 + in the limit of vanishing . We hope to come back in future work to the two-dimensional case to clarify these issues.
The computation of two-and three-point functions allows us to answer also the question about unitarity. Unitarity constrains the correlators [3]: we can check whether our two-point functions satisfy reflection positivity, whether the OPE coefficients appearing in our three-point functions are real, and so on. It turns out that these constraints are satisfied by most of the correlators we have computed, but we also found a surprise.
For the three-point function of one tetrahedron operator with two φ 2 operators, Eq. (6.34), we find a logarithmic dependence on positions of the type encountered in logarithmic CFT [64], rather than a standard power law. This strongly suggests that our large-N fixed point defines a logarithmic CFT and that the tetrahedron operator belongs to a logarithmic multiplet of exactly marginal operators. If correct, this in particular means that our large-N CFT is not unitary (hence quashing our original conjecture) as in a logarithmic multiplet the logarithmic dependence shows up already in the two-point functions, which thus cannot satisfy reflection positivity. Unfortunately, we have not been able to identify the other members of such a multiplet, hence we cannot conclude with certainty that we have a logarithmic CFT. However, given our proof of conformal invariance, and the fact that the logarithmic bahavior can only be compatible with it in the presence of a logarithmic multiplet, we see no alternative.
Melons vs fishnets. We conclude this overview of results by a remark. The appearance of logarithmic correlators in our computation has brought to our attention a number of similarities between our model and the conformal fishnet theory introduced in Ref. [65], or more precisely with its generalization to dimension d < 4 [66], which requires a long-range propagator like ours. The conformal fishnet theory is a model of two complex matrices, with a single-trace chiral quartic interaction, without its hermitian conjugate. The interaction is therefore complex, as in our model. Furthermore, such interaction is exactly marginal in the large-N limit, also as in our model. Moreover, renormalizability requires the introduction of doubletrace interactions, and the four-point function renormalizing them is built out of ladders and bubbles in a similar fashion to what was found in Ref. [1] for our model (see also Sec. 3.2.1 below). The resulting beta functions for the running couplings are therefore quadratic in both models, with coefficients parametrically depending on the exactly marginal coupling (compare Eq.(13) of Ref. [67], with our beta functions (3.42) below). Lastly, the conformal fishnet theory is a logarithmic CFT [68], which is what called our attention to the similarities between the two models. In the fishnet theory the logarithms appear due to the fact the large-N limit suppresses the tree-level contribution to certain two-point functions, leaving behind some logarithmic terms from the quantum correction, which then cannot be interpreted as a correction to the scaling dimension. We found exactly the same type of mechanism in our model, but for a three-point function.
There are of course also important differences. The different names that have been attached to the two conformal theories are not an accident: whereas at large-N our model is dominated by melonic diagrams, the model of Gürdogan and Kazakov is dominated by fishnet diagrams. What is accidental is the fact that for the four-point functions of fundamental fields both type of diagrams reduce to ladders: indeed a ladder can be thought either as melonic graph which has been open on two edges (and with resummed propagators), or as a fishnet with periodicity of length two in one direction. However, all the other n-point functions of fundamental fields are different. In particular, while the two-point function in a melonic CFT is given by a sum over melonic two-point diagrams, in the fishnet CFT there is no correction to the bare propagator (and hence no mass or wave function renormalization) at leading order. Such a difference is not very important for our long-range model in d < 4, which has no wave function renormalization anyway, but it becomes relevant for models with a standard short-range propagator, including ours at d = 4.
Plan of the paper. In section 2, we introduce and review the model. In section 3, we discuss renormalization and fixed points of the model, both in the Wilsonian picture and in the minimal subtraction scheme. Then, in section 4, we discuss the mixing of the φ 4 composite operators under the renormalization group flow. In section 5, we give a proof of the conformal symmetry at the infrared fixed point of a class of correlations, based on the D = d + p dimensional embedding method of Ref. [52]. Lastly, in section 6, we use the perturbative expansion at the fixed-point (our small parameter being the exactly marginal tetrahedron coupling) in order to compute the two-and three-point functions (and hence the OPE coefficients) among the φ 4 and φ 2 composite operators.
In the appendices we include some detailed computations and additional remarks. Appendix A contains detailed computations for several integrals we use in the main text. The beta functions for the higher spin bilinear operators are presented in appendix B. In appendix C, we discuss the large-N scaling of the maximally single-trace (MST) and maximally multi-trace (MMT) operators, and in appendix D, we discuss correlators containing the pillow operator, which is neither of MST nor MMT. Finally, in appendix E, we give a comparison with the long-range Ising model, that loosely corresponds to N = 1.

Overview of the Model
We study the tensor model of [1,2], that is, the O(N ) 3 tensor model of Klebanov and Tarnopolsky [26] and Carrozza and Tanasa [44] (CTKT model) with a long range covariance. The fundamental field is a real tensor field of rank 3, φ a 1 a 2 a 3 (x), transforming under O(N ) 3 with indices distinguished by the position, and we denote a = (a 1 , a 2 , a 3 ). The action of the model is: ab;cd + λ 2P where repeated tensor indices are summed over a i = 1, · · · , N and we introduced the projectors: ab;cd = 3(δ p ab;cd −δ d ab;cd ) ,P (2) ab;cd =δ d ab;cd . (2.2) and the rescaled operators: Here t stands for tetrahedron, p for pillow, and d for double-trace. Such names refer to the graphical representation of the respective pattern of contraction of indices, as recalled below. We use the following shorthand notations for the quadratic invariant: and for the quartic invariants: The difference between this model and the CTKT model is that the Laplacian is allowed to have a non integer power 0 < ζ ≤ 1. This modification preserves the reflection positivity of the propagator: the free theory is unitary for any ζ ≤ 1. The choice ζ = d/4 renders the quartic invariants marginal in any d [1]. It is this value of ζ that interests us in this paper. Moreover, we will restrict to d < 4 in the following, in order to avoid a wave function renormalization (see Sec. 3.1, and footnote 3 in particular).
We have not assigned any subscript to the coupling of the tetrahedral invariant, as it plays a special role in the model. Observe also that the infrared fixed point found in [1,2], and that we aim to study, occurs for a purely imaginary tetrahedral coupling, hence we have chosen here to make that explicit from the onset, by writing the coupling as i λ, with λ ∈ R.
As usual, it is convenient to introduce a graphical representation of the O(N ) 3 invariants, which also justifies their names. We represent every tensor (φ a , φ b and so on) as a three-valent node and every contraction of two indices (a i and b i for instance) as an edge with a color i = 1, 2, or 3 (red, green, or blue) corresponding to the position i of the indices. As a result, O(N ) 3 invariants are represented by 3-colored graphs. The graphs corresponding to the quartic invariants of Eq. (2.1) are depicted in Fig. 1. We can expand the free energy and the connected n-point functions perturbatively around the Gaussian theory. We introduce two graphical representations for the terms in the perturbative expansion. In the first, each interaction invariant is a 3-colored graph, and the Feynman propagators are represented as edges with a new color, which we call 0 (pictured in black), connecting the tensors. This leads to a representation of the perturbative expansion in terms of 4-colored graphs, as for example in Fig. 2. In the second representation, we simplify the graphs by shrinking each interaction invariant (i.e. all its edges with colors from 1 to 3) to a point. We call the resulting object Feynman diagrams, as they represent in a more straightforward way the spacetime integrals associated to the amplitude. An important class of graphs and diagrams are the melonic ones, which however have very different features depending on whether it is the 4-colored graph which is melonic, or the (single color) diagram. The pillow and doubletrace interactions are examples of melonic 3-colored graphs, and models based on such type of interactions are known to be dominated by melonic 4-colored graphs at leading-order in 1/N [45]. The corresponding Feynman diagrams are cactus diagrams, as in vector models. On the contrary, the tetrahedron interaction is not a melonic 3-colored graph, but it leads to melonic Feynman diagrams in the large-N limit [26,44].
As a result of the combination of pillow, double-trace, and tetrahedron interactions, our model has a 1/N expansion dominated by melon-tadpole diagrams [1] with melons based on couples of tetrahedral vertices and tadpoles based on either pillow or double-trace vertices (see Fig.4).

Renormalization and Fixed Points
We consider the d-dimensional theory in Euclidean signature, and our aim is to study the case ζ = d/4.
We introduce an infrared regulator µ by modifying the free covariance of the theory to: and we regulate the UV divergences by setting: This implies that the ultraviolet dimension of the field is The reader should keep in mind that for us is just a regulator and we always intend to send → 0 in the end. In particular, we are not looking for Wilson-Fisher type of fixed points at small but finite . We will need the following Fourier transform, which holds for d/2 > γ > 0: .

The two-point function
We now discuss the bare and the full two-point functions of the model.
The bare propagator. The bare propagator in the direct space, with infrared regulator set to zero, is obtained from Eq. (3.4) by simply setting γ = ζ: (3.5) We will encounter below the convolution of the cube of the bare propagators with another propagator. Using repeatedly the Fourier transform in Eq. (3.4), we obtain the formal result: The problem with this formula is that we used Eq. (3.4) for γ = 3ζ − d = −d/4 + 3 /4 < 0, hence the result is only formal. In fact this convolution hides an ultraviolet divergence which needs to be subtracted.
Taking this into account, we obtain: 3 where this time the integral is convergent, and we have used the distributional limit lim →0 /|x − y| d−2 = δ(x − y) (alternatively one can take the limit in Fourier space to reach the same conclusion).
The full two-point function. We now discuss the full two-point function of the model. In the absence of spontaneous symmetry breaking, 4 the full two-point function is diagonal in the tensor indices. That is, it can be written as: Due to the infrared regulator, the full two-point function also aquires a µ-dependence, hence we write its regulated diagonal component as G µ (x − y).
Following [1], we observe that G µ (p) at leading order in N respects the melonic Schwinger-Dyson (SD) equation: where the sign of the last term is changed with respect to [1] due to the explicit i factor in our action. Similar to Eq. (3.6), this equation exhibits a power divergence in the ultraviolet. The critical theory is obtained by tuning the bare mass to exactly cancel this divergence such that the renormalized mass is zero. The melonic SD equation in the limit µ → 0 for the critical theory is: where G(p) denotes the full two-point function with no cutoff. It is solved self consistently at = 0, for any d < 4 by the ansatz G(p) −1 = Zp 2ζ : provided that the constant Z is chosen as: In Fourier space the subtracted melon contribution is: where the integral over the ai parameters is convergent for d < 4 and < d/3. We computed it in Appendix A. Multiplying Eq. (3.7) by p −2ζ , Fourier transforming back to the direct space (which is allowed for ≥ 0) and using the analytic continuation of the Γ function we obtain Eq. (3.7) 4 Spontaneous symmetry breaking in tensor field theories has been so far not much explored, but see Ref. [17,69,70].
This in turn fixes the bare mass to m 2ζ = m 2ζ c , where: Summarizing, in the critical theory, the net effect of all the Feynman diagrams on the two-point function is a just multiplicative factor: G(p) = Z −1 C(p). The constant Z can be seen as a finite wave function renormalization which resums all the melonic insertions. By resumming an infinite series of diagrams, Z carries non-perturbative information on the radius of convergence of the perturbative series [1] at large N .
The fact that we can perform such a resummation is an important property of our tensor model at = 0. Nevertheless, in order to regulate the UV divergences, below we will consider > 0. In this case, the resummation of the leading-order two-point function can not be performed explicitly. Furthermore, for theories with a non-local propagator, due to the locality of counterterms, there is no need to renormalize the fields. Therefore, it is common to set the the wave function renormalization to 1 (e.g. [52]), and we will do the same here. We thus seem to have a discontinuous Z at = 0, but in fact this is not the case. It is all only a matter of repackaging of diagrams: the two-point melonic insertions contribute also at > 0, but as we let → 0, we can conveniently resum them by working with melon-free diagrams and propagators divided by Z. Differences appear if we decide to rescale or not the fields by √ Z: rescaling is a natural choice at = 0, as it brings the two-point back function to its original form, but it is not justified at > 0 as G(p) and C(p) have different functional forms. In any case, rescaling by a finite constant is just a choice of renormalization scheme like others, and as such it will not affect the universal part of the beta functions. For = 0 the melonic SD equation in direct space is: which can also be written as: The SD equation simplifies if the integral over z is understood in the sense of dimensional regularization. In this case the local part of the melonic correction and the tadpoles are set to zero. In the Wilsonian picture the melon integral combines with the explicit mass counterterms which provide the subtraction. This works as long as the position y is an internal position. If y is an external argument, that is when this is a contribution to a correlator . . . φ 4 t (y) one needs to replace the tetrahedral operator by the renormalized tetrahedral operator:

The beta functions
We review the β functions of the theory defined by the action (2.1) tuned to criticality.

Quartic couplings
We denote i Γ R t , Γ R 1 , and Γ R 2 the appropriately normalized one-particle irreducible four-point functions at zero external momentum. They are computed using the bare expansion in terms of connected amputated one-particle irreducible four-point diagrams G with amplitude: where n(G) denotes the number of vertices of G, e ∈ G denotes the edges of G, and T runs over the spanning trees in G (each having n(G) − 1 edges). The tetrahedral four-point function is trivial, as it receives no radiative corrections at large N : At leading order in 1/N the remaining four-point functions Γ R 1 and Γ R 2 are identical up to replacing λ 2 by 3λ 2 andλ 1 byλ 2 . We discuss Γ R 1 . Only chain diagrams [1] contribute to Γ R 1 at leading order in 1/N . A chain diagram G is a sequence of irreducible pieces connected by pairs of parallel horizontal edges. The irreducible pieces are either vertical ladder rungs with two tetrahedral couplings, or bare vertices λ 1 . There are 2 n chain diagrams with n irreducible parts (that is vertical rungs or bare vertices). The edges of a chain diagram are decorated by arbitrary melonic insertions, but we do not include tadpoles, as they are assumed to have been taken care of (by dimensional regularization or mass subtraction).
The chain diagram consisting in a bare vertex has amplitude 1. We denote n t (G) the number of tetrahedral vertices of G (which is always even), n 1 (G) the numbers of vertices λ 1 of G, and G the set of connected chain diagrams with at least two internal vertices. We have: The case = 0 is special. The series in (3.20) can be further simplified: all the melonic insertions can be analytically resummed at the price of multiplying the bare propagator by Z −1 from Eq. (3.12). Furthermore, in this case the four-point function itself should be divided by Z 2 . The overall effect is that Eq. (3.20) can be rewritten at = 0 by dividing all the couplings by Z 2 and reducing the sum to chain diagrams with no melonic insertions.
The chain diagrams can be analyzed in terms of their one-vertex irreducible components. Adapting the notation of [1] to the case > 0, we denote: •Û r the sum of dimensionless amplitudes of the ladders with r ≥ 1 rungs, and with melonic insertions; we include inÛ r theλ-dependence due to melonic insertions, but not that due to the pairs of vertices in a rung. Therefore, we write the sum over the ladders of arbitrary length as: is not the generating function of the amplitudeŝ U r (x), due to the x-dependence of the latter. We also define U r ≡Û r (0), for the amplitudes without melonic insertions (Fig. 5). •Ŝ r the sum of dimensionless amplitudes of the caps with r ≥ 1 rungs, i.e. ladders with r rungs closed on a λ 1 vertex on one side, with melonic insertions; we write the sum over caps as: with x = µ − λ . We also define S r ≡Ŝ r (0), for the amplitudes without melonic insertions (Fig. 6). Figure 6: The caps S 1 , S 2 , and S 3 . The blue vertex represents λ 1 .
•T r the sum of dimensionless amplitudes of the double-caps with r ≥ 0 rungs, i.e. ladders with r rungs closed on a λ 1 vertex on each side, with melonic insertions; we write the sum over double-caps as: We also define T r ≡T r (0), for the amplitudes without melonic insertions (Fig. 7). Observe that for the caps and double-caps exactly one and two vertices, respectively, correspond to couplings λ 1 . The leading 1/ behavior of these amplitudes is U r ∼ −(2r−1) , S r ∼ −2r , T r ∼ −(2r+1) . In terms of the resummed amplitudes, the bare expansion writes [1]: which exhibits an ascending series of poles in 1/ . In the limit → 0, as discussed above, we can restrict to ladders, caps, and double-caps with no melonic insertions (that is, we can remove hats fromÛ r ,Ŝ r , and T r ), provided that the couplings and the amplitude are rescaled by Z −2 .
Wilsonian beta functions. In the Wilsonian picture we identify the four-point functions with the running couplings: and for the theory at = 0 we would further rescale the right-hand side by Z −2 . The beta functions are the scale derivative of the running couplings at fixed bare couplingsλ,λ 1 ,λ 2 . The beta functions ofg 1 and g 2 decouple. In order to compute them, we invert the bare series: .
We first observe that the β function of the tetrahedral coupling is trivial: In particular β t is identically zero at = 0 indicating a line of fixed points in that case [1]. For the other couplings, the beta functions are obtained by taking µ d dµ in equation (3.26), leading to: where U = −g 2Û 1 + . . . , S = −g 2Ŝ 1 + . . . and T =T 0 −g 2T 1 + . . . and their derivatives are evaluated at the effective tetrahedral couplingg. Comparing with the results of [1], we see that the beta functions remain a quadratic polynomial ing 1 andg 2 also at > 0: where the coefficients are: (3.31) Up to cubic order in the couplings (that is, two loops), and using the integrals D ≡ T 0 (0) = U 1 (0) and S 1 (0) discussed in Appendix A, we find: and at first order in the coefficients β (0) , β (1) and β (2) are: which reproduce the results of [1] upon dividing all the couplings by Z 2 .
Minimal subtraction. The minimal subtraction consists in fixing a series of counterterms for each coupling, having ascending series of poles in 1/ , with residues chosen such that the four-point functions expressed in terms of the renormalized couplingsg,g 1 ,g 2 do not have any poles in 1/ . The tetrahedron coupling is still trivial, as before. For the others we have: where the B (k) (g 1 ,g) are -independent. A standard manipulation [71] then leads to: The difference between minimal subtraction and Wilsonian scheme amounts to a mapping between renormalized couplings. In fact, in minimal subtraction we find: where we have used: and similar expansions for S (0) (g) and T (0) (g). The superscript here denotes the order of the pole in , and each amplitude U r and so on has an expansion of the form (3.37). Therefore, we have: We could in principle plug this transformation in Eq. (3.26), find the expansion in Eq. (3.34), and write the beta function in terms of the one-vertex irreducible amplitudes. The result is rather cumbersome. Let us instead compute directly the β functions in minimal subtraction up to two loops (order λ 3 ). At this order, the beta functions of marginal couplings are scheme independent up to O( ) terms because F (g 1 ,g) =g 1 + O(g 2 ). The bare series of Γ R 1 at cubic order is: They are cancelled by choosing: which leads to: The reader can check that these functions coincide with the Wilsonian beta functions in Eq. (3.27) and (3.32), up to terms that vanish for → 0. At higher perturbative orders, even at = 0, the two schemes will differ. However, also in minimal subtraction, the β functions at all orders have the form: Fixed points. Without loss of generality we considerg > 0. While the fixed point below exists at all orders in perturbation theory, we will restrict to the first non trivial order. The beta functions at one loop: admit at = 0, for anyg =g , a fixed pointg 1 =g ,g 2 = √ 3g which is infrared attractive: The quartic operators acquire anomalous scaling ∆ a = 4∆ φ + δh a (with a = t, 1, 2). Here ∆ φ = d/4 is the dimension of the field at the fixed point. As the stability matrix ∂g a β b is triangular, δh a = ∂g a β a | hence:

Quadratic couplings
Let us consider a φ 2 perturbation of the critical theory. This comes to considering from the onset a mass parameter m 2ζ = m 2ζ c + λ φ 2 , with m 2ζ c given by (3.13). With respect to our previous discussion we now need to add bi-valet vertices in the theory, corresponding to the insertion of λ φ 2 . The back reaction of λ φ 2 on the flow of g 1 and g 2 can be neglected, as the ultraviolet behavior of the massive propagator [(p 2 + µ 2 ) ζ + λ φ 2 ] −1 is identical to that of the massless one.
Let us consider the one-particle irreducible two-point function at zero external momentum Γ R (2) (0) = G −1 µ (0). At one loop, only the tadpole with an insertion of a bi-valent vertex contributes: with D the integral of App. A. For both the Wilsonian prescription and in minimal subtraction, we get: which is valid for ≥ 0, and where 2∆ φ = d− 2 is the classical dimension of the φ 2 operator. Observe thatg φ 2 = 0 is always a fixed point of this equation. For = 0 (when the dimension of the field becomes ∆ φ = d/4), the beta function becomes: and therefore, close to the fixed pointg 2 = √ 3g the operator φ 2 acquires the anomalous scaling reproducing the anomalous dimension found in [1] by diagonalizing the four-point kernel. 5 The spin-zero bilinear operators of the type: φ(−∂ 2 ) n φ can be treated similarly. We present a detail computation of the corresponding beta functions in Appendix B. The result again reproduces the anomalous dimensions of this type of operators derived in [1] by diagonalizing the four-point kernel.

Operator Mixing
Consider a field theory described by the action: where S 0 is some free quadratic action and the perturbation is a sum over local operators O a . To simplify the discussion we assume that we do not have a wave function renormalization, which is in particular the case for our model. We denote ∆ a the dimension of the operator O a (x). As already mentioned, in minimal subtraction we replace the bare couplings by renormalized couplings plus counterterms: where the counterterms of λ a can depend on all the renormalized couplings {g b }. The renormalized action is obtained by substituting the bare couplings: and the counterterms are chosen such that the connected correlations: 4) have no poles in 1/ when expressed in terms of the renormalized couplings. The (integrated) renormalized operator O a is the derivative of the renormalized action with respect to the dimensionless renormalized coupling g a : and it is easy to check that when acting on a connected correlation the derivative with respect to g a brings down an O a operator: it has no poles in 1/ , its derivative is also finite, with the possible exception of some critical points. The bare operators are linear combinations of the renormalized ones: where M is the mixing matrix. For a free theory, that is, neglecting all the radiative corrections, the mixing matrix is M ab = δ ab µ ∆ b −d . In the interacting case, M is determined by observing that, as insertions in correlation functions, the following chain of equalities holds: where we used the fact that the renormalized operators are linearly independent. Finally, the β functions write in terms of the mixing matrix: Let us go back to the case of the action in Eq. (2.1). The four quartic operators defined in Eq. (2.6) have classical dimensions ∆ a = 4∆ φ = d − and the beta functions are β b = − a=t,1,2 λ a M ab . In particular, the bare perturbation can be written in terms of renormalized operators as: where we ignore the mixing with the φ 2 operator 6 . 6 The mixing matrix simplifies in our case. The β functions are βt = − g and βi = − gi + β i g 2 i but

First order action
Gathering the first order results in Eq. (3.40) and (3.48), the minimally-subtracted action, with bare mass m 2ζ = m 2ζ c + λ φ 2 and m 2ζ c given by (3.13), and up to quadratic order in the couplings, is: where ∆ φ = d− 4 and the bare couplings at first order are: Observe that in the second line of Eq. (4.11) we can recognize the Wick ordered double-trace interaction at leading order in N and up to a vacuum term 7 : . At this order in 1/N it is superfluous to Wick order the other quartic interactions. The renormalized operators are computed by taking the derivatives of the renormalized action: 14) where the (inverse) mixing matrix elements are: (4. 15) in order the keep the formulae slightly more general we will not substitute their explicit form. Using 0 = λ + ∂gλ βt and 0 = λi Arranging the couplings in the order t, 1, 2, we have: and the renormalized operators: The standard definition of Wick ordering [72] gives: In this section, we give a proof of the conformal symmetry of the model (2.1) at the infrared fixed point, following the discussion of Ref. [52]. The main idea is that the model (2.1) can be written with a standard short-range kinetic term by embedding it in D = d + p (where p = 2 − 2ζ) dimensional space. In this enlarged space, the action becomes: where the D-dimensional coordinates are labeled by X M = (x µ , y m ) and the original field is obtained by Φ| y→0 = φ. In this D-dimensional space, one can write down a local energy-momentum tensor: where δ M N = δ µν if both indices are in the d-dimensional space, and zero otherwise. We also introduce the orthogonal projector δ ⊥ M N = δ M N − δ M N . Now we can write the divergence and trace of the energymomentum tensor as: ab;cd + λ 2P ab;cd + λ 2P ab;cd + λ 2P ab;cd + λ 2P We note that E and E N are proportional to the equation of motion, so they vanish on-shell. Therefore, on-shell the trace of the energy momentum tensor, Eq. (5.4) is equal to a double divergence, ∂ 2 K φ 2 a , plus a term proportional to , and thus the theory is classically conformal invariant at = 0. However, due to radiative corrections, the Φ 4 | y=0 = φ 4 operators lead to infinities, hence they need to be renormalized as we discussed in the previous section. Using Eq. 4.10 we arrive at the following expression for the trace of the energy-momentum tensor: where we observe that in terms of the renormalized operators, the last terms survive at = 0, leading in general to a breaking of conformal invariance. The dilatation and special conformal transformation currents are constructed from the energy-momentum tensor as: Their divergences are found as: Inserting these operators in a renormalized n-point function and integrating over the insertion point, up to a possible boundary term we obtain: The Schwinger-Dyson equations involving E and E M are obtained from the path integral expression for the n-point function by the field redefinitions Φ → Φ(1 + δΦ) and Φ → Φ + (∂ M Φ)δΦ, respectively. The result is [56]: 14) Using these equations, together with the expression of the trace of the energy-momentum tensor (5.8), we obtain the Ward identities for the dilatation and special conformal transformation currents. Since the n-point function behaves continuously in the limit y → 0 [52], we can write down the Ward identities in the original d-dimensional space as: These Ward identities are well defined even at = 0 and the right-hand sides survives in this limit. We thus get in general a breaking of scale and special conformal invariance due to renormalization. At a fixed point, as β t (g) = β 1 (g, g 1 ) = β 2 (g, g 2 ) = 0, the invariance is restored. In fact, one should check that the integral on the right-hand side of the Ward identities does not blow up when we approach the fixed point as the inverse of the beta functions or worse. This subtlety is discussed in detail in Ref. [52] and that discussion goes through in our case. 8 The main point of repeating the derivation of the Ward identities here is to check that the structure of right-hand side (beta function times operator insertion) generalizes to our model with multiple interaction terms.
We conclude that the n-point functions of fundamental fields are conformal invariant at the fixed point. The next step would be to generalize the above result to correlators of composite operators. A possible strategy, starting from correlators of fundamental fields and using an operator-product expansion to generate the composite fields, is sketched in Ref. [52]. In the next section we will follow a different route, and compute perturbatively three-point functions of (quartic and quadratic) composite operators.

Correlation Functions
In this section we will explicitly compute large-N three-point functions among the φ 4 and φ 2 composite operators we discussed in the previous sections at the interacting IR fixed point. Since this fixed point depends parametrically on the exactly marginal coupling g, we will work at small g, at lowest perturbative order in the couplings.
Large-N and conformal limits. Before proceeding, we stress one important aspect about the interplay between large-N limit and conformal limit. Namely, it should be kept in mind that in principle we should take the conformal limit (i.e. tune the theory to the fixed point) before or at least together with the large-N limit, in order to keep non-trivial two-point functions in the limit. For example, if we were to take the large-N limit first, away from criticality, the two-point function of double-trace operators would be dominated by the diagram in Fig. 8, rather than that of Fig. 11. Away from criticality, such diagram leads to a finite contribution (after renormalization) determining the two-point function at leading order of the large-N expansion. However, at criticality its renormalized amplitude vanishes, and we would be left with a zero two-point function. The way out is to assume from the beginning that diagrams such as that in Fig. 8 are zero because we are in the massless limit, and to rescale the operators in such a way that the first non-vanishing contribution to a two-point function at criticality is always normalized to order N 0 . Making sure that such a scaling exists in the case of tensor models is non-trivial, and it is what we want to address now. Let us first review how correlators of invariants are treated in the case of matrix field theories. The analogue of our interacting action (2.1) is: In our case, we can actually understand the absence of such singularity on the right-hand side of the scale Ward identity (5.16) thanks to the presence of a small parameter, the tetrahedron coupling g. First, we can express the right-hand side of (5.16) by means of (4.6), thus reducing the scale Ward identity to a Callan-Symanzik equation. Then, for infinitesimal g, we notice that n-point functions near the fixed point are essentially polynomials in g, hence their derivative at the fixed point is finite. A singularity might instead arise at finite g, and in particular at the critical value determined by the invertibility of the relation g = λZ(λ) −2 , with Z(λ) defined by the solution of Eq. (3.12) (see Ref. [1]).
It is convenient to rescale φ → √ N φ, so that the action becomes: In such rescaled variables, the functional: and have an expansion in 1/N starting at order N 0 . Restricting to single-trace operators, the leading order graphs are connected planar ribbon graphs (or 3-colored graphs in our notation i , without explicit factors of N . They are obtained by rescaling appropriately the derivatives of W , or equivalently, by taking derivatives with respect toJ i }]. In this normalization, the connected two-point functions of single-trace operators at leading-order are of order N 0 , while their higher-point functions are suppressed in 1/N . For multi-trace operators, a subtlety rarely emphasized appears, namely, their leading-order two-point functions naively seem to be of higher order than N 0 , and higher-point functions seem to have even higher powers of N . In fact, we have: with a leading order contribution scaling as N 2+ j (t j −2) . For n = 2 and t j = 1, we get a correlator of order N 0 , as anticipated. On the other hand, if for example all the t i are greater than two, we naively seem to get an arbitrarily large power of N . However, the (naive) leading-order contributions for such multi-trace operators come from cactus diagrams, as the one of Fig. 8, which always contain a factor corresponding to the one-point function of a single-trace component of the multi-trace invariant (the end leaves of the cactus). In the conformal limit such contributions vanish, and a factorization into two-point functions of single-trace components dominates. 9 Notice that since two-point functions of single-trace operators are normalized to one, by factorization also the two-point-functions of multi-trace operators are normalized to one. Summarizing, one is typically interested in: The factorization property is often stated in a slightly different fashion in the literature [73]. Namely, considering only single trace operators, one observes that which gives two-point functions normalized to one, and higher-point functions suppressed in 1/N . We now go back to our tensor field theory. In analogy to the matrix case, we rescale the field as: so that the action becomes: where we reinstated the pillow and double trace couplings. Observe that the operators δ t , δ p and δ d , defined in Eq. (2.4), merely contract indices and do not contain explicit factors of N .
According to Appendix C, the generating functional of connected correlators admitting a large-N limit in the perturbative sense is now: with ρ b ≥ 0 chosen according to the optimal scaling defined in Ref. [44]: where F b counts the total number of cycles of alternating colors i and j with i, j ∈ {1, 2, 3}. Invariants with ρ b = 0 are called maximally single-trace (MST) [9]. Defining again correlators as derivatives with respect to non-rescaled sourcesJ b = N we have the following expansion: For MST operators we thus have an analogous result as for single-trace matrix operators, that is, their twopoint functions are of order one, and higher-point functions are suppressed. Assuming again that one-point functions are zero, for operators which are products of MST operators, which we will call maximally multitrace (MMT), we obtain a factorization property as for the matrix multi-trace operators (see Appendix C for more details). For non-MMT operators with ρ b > 0 we seem to have again two-point functions with a higher scaling than N 0 , and for those with ρ b > 3/2 we seem to have again a possible arbitrarily large power of N . However, we conjecture that, as for multi-trace matrix operators, the leading order diagrams (those with ω = 0) vanish in the conformal case, and the first non-vanishing order has ω ≥ 3 − i ( 3 2 − ρ b i ). We hope to return to this conjecture in full generality in a future publication. For the pillow invariant which is contained in φ 4 1 , and which has ρ p = 1/2, we can explicitly check that for example the two-point function at leading order is of order N 0 (see App. D). Now we proceed to the perturbative computation of the two and three-point connected correlations of the MST operators φ 4 t and φ 2 . The correlations involving the double trace φ 4 2 are obtained at leading order using the large N factorization. We do not study here the correlations involving φ 4 1 , which due to the pillow operator is neither MST nor MMT, but we will briefly discuss them in App. D. As explained above, a connected correlation function of n MST operators has an a priori scaling N 3− 3 2 n , that is N 0 for two-point correlations and N −3/2 for the three-point ones. Whenever a correlation scales less than that, we consider it suppressed in the large N limit.
Rescaled action and renormalized operators. In order to regularize the UV divergences we set > 0. We compute the correlations up to first order in the couplings in the critical theory, λ φ 2 = 0. At the relevant order the action in Eq. (4.11) becomes: where the operators are defined as before, except for the factors of N due to the rescaling (6.7) (remember the distinction between hatted and un-hatted deltas in Eq. 2.3): where, in spite of the new scaling in N we maintain the same notation for the quartic invariants. At this order we have λ = µ d−4∆ φ g, λ 1 = µ d−4∆ φ g 1 , λ 2 = µ d−4∆ φ g 2 with g, g 1 , g 2 the dimensionless couplings. Below we will use λ's for the perturbative expressions at > 0, and g's for the final expressions in dimensionless variables. At linear order, the renormalized operators are: 10 14) with Q = 2(4π) −d/2 Γ(d/2) −1 . All our correlations will be written up to terms of order g 2 which from now we omit. We will fix the pillow and double-trace couplings at their IR fixed point, which at the same order in N and g reads g 1 = g , g 2 = √ 3g .
Contact terms. Because of the melonic convolution in Eq. (3.7), we will encounter "contact" terms proportional to (∂ 2 ) n δ(x ij ). This should not come as a surprise, as such terms are for example expected in the three-point function of exactly marginal operators [61]. One should distinguish them from terms like lim →0 1 |x| d− ∼ 1 δ(x), which require regularization and lead to an anomaly [59,74]. Contact terms with finite coefficients are compatible with conformal transformations. The difference can be reformulated by the observation [60] that in d dimensions the only homogeneous distribution of dimension d is the Dirac delta, while any generalized function which coincides with 1 |x| d for |x| = 0 is an associate homogeneous distribution [75], that is, under scale transformations it transforms with an in homogeneous contact term.

Two-point functions
We denote by prime dimensionless positions, x = µx and so on. We first compute the relevant two-point functions.
The [φ 2 ][φ 2 ] correlation. Up to the first order in the coupling constant, at leading order in N only the two diagrams represented in Fig. 9 contribute. They give: where Z φ 2 is the renormalization constant of the φ 2 operator in Eq. (6.14). Computing the integral at first order in (see appendix A.1) we obtain: As expected, the 1/ pole cancels as, according to Eq. (6.14), Z φ 2 = 1 + Q g 2 . We now take to 0. The g 2 log |x − y | term combines with the constant term to give a correction to the scaling law, hence at first order in g the two-point function is: where ∆ φ = d/4 is the dimension of the field at = 0. The dimension of φ 2 is ∆ φ 2 = 2∆ φ + δh φ 2 with δh φ 2 = Q g 2 which reproduces the anomalous dimension found in Eq. (3.50).
Observe that we can normalize this two-point function to 1 by introducing the "normalized renormalized" [φ 2 ] operator: correlation. Next, we consider the two-point functions of φ 4 t operators. As the operator is complex we take the two-point function of φ 4 t with the hermitian conjugate (φ 4 t ) † = −φ 4 t . At first order in the couplings and at leading order in N only one diagram, depicted in Fig. 10, contributes, yielding: The → 0 limit is trivial. We can normalize this two-point function introducing [φ 4 an exactly marginal operator, leading to a conformal anomaly in even dimensions [63]. The appearance of an anomaly can be understood in the spirit of our earlier comment on contact terms by noticing that for d = 2n we have 1/|x| 2d− ∼ 1 (∂ 2 ) n δ(x). Since in our model we have assumed d < 4 from the beginning, the anomaly only concerns d = 2.
The [φ 4 2 ][φ 4 2 ] correlation. At the appropriate order we have two graphs (see Fig.11), where the quartic Figure 11: vertex is : φ 4 2 :, yielding: Following the same steps as before we get, after cancellation of the pole and taking the limit → 0: with δh 2 = 2 Q g 2 reproducing the perturbative result Eq. (3.46). As expected for a multitrace operator we have For a conformal field theory we expect that the two-point functions of operators with different dimensions are zero. We check this for [φ 2 ][φ 4 t ] . The leading-order contribution to this correlator is depicted in Fig. 12 (and we get the appropriate counterterm subtraction from the mixing of φ 4 t with φ 2 ) which yields: Using the convolution (3.7), one can evaluate the z-integral and obtain: which is zero in the → 0 limit. We conclude that [φ 2 ](x)[φ 4 t ](y) N −1 , hence it is suppressed at large N .

Three-point functions
We now compute the three-point functions at first order in the couplings and leading order in large N .
Up to the first order in the coupling constant, the only relevant Feynman diagram at order N −3/2 is depicted in Fig. 13 (and we get the appropriate counterterm subtraction from the mixing of φ 4 t with φ 2 ). Its Figure 13: contribution is given by: where we used again Eq. (3.7). This is of the expected form for the three-point function of exactly marginal operators [61], and up to contact terms with finite coefficients, it is consistent with the conclusion of Ref. [59], according to which a marginal operator should have vanishing three-point function with itself.
The They yield: Using Appendix A.1, the right and side of the above equation becomes: where the dots denote some finite part. As expected, the pole cancels and in the → 0 limit we get: which is the right form of the conformal three-point function with the correct anomalous dimension given in Eq. (3.50).
We have: (6.29) Using the conformal integral given in Appendix A.1, the right hand side above writes: Again the pole cancels and in the → 0 limit we get: which is of the correct conformal form.
Taking into account the operator mixing (6.14), we get: In order to obtain this equation we note that at order N −3/2 and up to linear order in the couplings the term coming from the mixing with φ 4 1 vanishes and φ 4 t φ 2 φ 2 starts at order g, hence we can ignore the renormalization of φ 2 . Only the two diagrams in Fig. 16 contribute at this order, and we find: Figure 16: where we have omitted the subtraction term in the first line, coming from the mixing of φ 4 t with φ 2 . Observe that the radiative correction (the first term in the second line) has an overall + sign which comes from taking into account the fact that the tetrahedral invariant is imaginary. For the first line, we use once more Eq. (3.7), thus obtaining contact terms proportional to δ(x − z)/|x − y| d . For the second line, we use the computations in App. A.1. The pole cancels, and in the → 0 limit we get: The contact terms have the correct conformal form at leading order and, as we expect higher-order corrections to modify them to δ(x − z)/|x − y| d+2δh φ 2 , we are not worried about the power 4∆ φ = d. However, contrary to the previous examples, the logarithms in the last line of Eq. (6.34) cannot be recombined into a correction to scaling of the form |x − y| g ∼ 1 + g ln |x − y|. In fact we miss a zeroth-order term in the coupling (the 1 from 1 + g ln |x − y|) because this term is subleading at large N , hence we are left with a logarithmic function. 11 This is a valid form of conformal three-point function if φ 4 t belongs to a logarithmic multiplet, that is, if our fixed-point theory is a logarithmic CFT [64]. Unfortunately, we have not been able to identify a logarithmic partner for the tetrahedron operator, hence we cannot unequivocally conclude that the fixed-point theory is a logarithmic CFT. On the other hand, as we see no obstruction in extending the proof of conformal invariance presented in Sec. 5 to composite operators, this must be the case.
A similar mechanism for the appearance of the logarithms in Eq. (6.34) is encountered in fishnet conformal field theory [67]: one finds that the two-point functions of certain dimension-five operators (in four spacetime dimensions) are such that the free-theory contribution is subleading in N with respect to the perturbative correction and thus some logarithmic behavior appears [68]. In that case, such phenomenon was found directly in two-point functions and the logarithmic multiplet has been identified, hence there is no doubt that the fishnet theory is a logarithmic CFT. In view of several similarities that one finds between our model and the generalization of fishnet theory to any dimension [66], such as the exactly marginal complex interaction and the structure of four-point functions of fundamental fields, we posit that a logarithmic multiplet exists also in our model. One possibility is that the logarithmic partners of the tetrahedron operator are non-invariant quartic operators, that is, operators built out of four fundamental fields, but with not all the indices contracted to form an invariant. The vector O(N ) model in the N → 0 limit is an example of logarithmic CFT in which the logarithmic multiplet is made of invariant and noninvariant operators of this sort [76]. We hope to come back to this question in the near future.

A Integrals
We compute in this appendix several integrals we encountered in this paper.
The subtracted melon integral. We first compute the subtracted melon integral: a 1 a 2 + a 1 a 3 + a 2 a 3 1 a 2 a 3  a 1 a 2 +a 1 a 3 +a 2 a 3 . (A.1) Using a Taylor expansion with integral rest we have: 1 a 2 a 3  a 1 a 2 +a 1 a 3 +a 2 a 1 a 2 +a 1 a 3 +a 2 a 3 .
The integral over t converges for d > 3ζ. Let us compute a slight generalization of this integral to q parameters α. Changing variables to β = a −1 and integrating out t yields: Introducing x = i β i and β i = s i x the integral becomes: We now use: and we finally obtain the subtracted integral: The D integral. We will repeatedly use below the integral 12 : which is convergent for 2Re(u) > Re(γ) and Re(u) > 0. In the the particular case u = ζ, γ = d/2 we get: Denoting ψ the digamma function (the logarithmic derivative of Γ) we have: 12 We have: The S 1 integral. The next integral we want to compute is: We use Mellin parameters to write: and Eq. (A.8) allows us to integrate a and b. We thus obtain: In the right half complex plane the integrand has poles at z = n, n + /2. Only the poles in 0, /2 have large residues at → 0, hence: .14) and the last line is finite in the → 0 limit. At small we get: In particular, we have The melon integral with momentum insertion. We are interested in evaluating the coefficient of p 2n in the Taylor expansion of the integral: Using the parametric representation, and observing that (q 2 1 ) n e −a 1 q 2 1 = (−∂ a 1 ) n e −a 1 q 2 1 , the momentum integrals can be computed, yielding: The coefficient of p 2n in the Taylor expansion of this integral is with: In order to compute the leading divergence of S (n) 1 , we note that: 20) and, as we encounter no singularities, we move the integration contour to z = −d/2 + i R. We thus get: and using Eq. (A.8) to integrate out the a's and b's we obtain: The only pole of the integrand in the right-half complex plane with residue of order 1/ is located at z = /2. Moving the contour across the pole we get:

A.1 The conformal integrals
We work at > 0, that is ∆ φ = (d − )/4. In the main text we encounter the following integrals involving bare propagators: and where κ is an -independent constant which we have not computed explicitly. These integrals are computed using two conformal integrals. First we have: which follows from the Fourier transform (3.5). In particular for ν 1 = ν 2 = 2∆ φ = d− 2 we get: and multiplying c(ζ) 4 and rearranging the coefficient we obtain Eq. (A.24).
In order to prove Eq .(A.25), we start from: (A. 28) In particular, we are interested in the case ν 1 = ν 2 = ∆ φ = (d − )/4, ν 3 = 2∆ φ = (d − )/2. Using the Mellin-Barnes representation [77], we rewrite the integral as: The integration contour (i.e. the constant c) is chosen to separate all poles of the first four Gamma functions from the poles of the last two Gamma functions. We close the contour to the right so that we pick up all poles of the first four Gamma functions, but none of the poles of the last two Gamma functions. The relevant poles are located at: with n 1,2 , m 1,2 = 0, 1, 2, · · · . The complete answer for the integral is given by the sum of all of these pole contributions. This is a daunting task to complete, so we look at the singular contribution in the limit → 0 with the choice ν 1 = ν 2 = ∆ φ , ν 3 = 2∆ φ . We note that for this choice one of the Gamma functions in the overall coefficient becomes: The leading behavior of the integral is O( −2 ), coming from the poles at n 2 = m 2 = 0. Namely, the poles of the third and fourth gamma functions at s = t = (3 − d)/4 lead to a O( −1 ) contribution from each of the fifth and sixth gamma functions. Overall we get: where κ is a constant. Multiplying by c(ζ) 4 , using ∆ φ = (d − )/4, and rearranging the coefficient we obtain Eq. (A.25).

B The Bilinear Operators
The spin-zero bilinear operators of the type φ(−∂ 2 ) n φ can be treated similarly to the φ 2 perturbation. We start by including a bare perturbation: and we evaluate the Taylor coefficient of p 2n in the one-particle irreducible two-point function, which we denote Γ R (n) . As the tadpole is local, only the melon with one bi-valent vertex λ (n) inserted on one of its edges contributes: where the subscript p 2n signifies that we are only interested in the coefficient of p 2n in the Taylor expansion of the integral. Using appendix A, the bare expansion becomes: Similar to the mass parameter, we obtain: where ∆ (n) = 2∆ φ + 2n is the classical dimension of the operator φ(−∂ 2 ) n φ. We note thatg (n) = 0 is always a fixed point of this equation, and that at = 0 the beta function simplifies to: At the fixed pointg the operator φ(−∂ 2 ) n φ acquires an anomalous dimension reproducing the results derived in [1] by diagonalizing the four-point kernel.

C The 1/N Expansion Revisited
In the main body of the paper we are interested in tensors of rank D = 3. However, the discussion below applies to any rank D.
A D-colored graph [7,21,22] is a graph such that: • all the vertices are D-valent • the edges have a color 1, . . . D and at any vertex we have exactly one incident edge for each color A D-colored graph is connected if any two vertices are joined by a path of (colored) edges such that two consecutive edges in the path share a vertex. For the D-colored graph h we denote V (h), E(h), C(h) and F (h) the numbers of vertices, edges, connected components and faces (i.e. bi colored cycles) of h. We also denote E c (h) the number or edges of color c and F c (h) the number of faces which contain the color c.
Invariants and Feynman graphs. The O(N ) D invariants Tr b (T ) are D-colored graphs [7] b. The vertices of b are associated to the tensors T and the edges (colored 1, . . . D) are associated to the contractions of indices: where v runs over the vertices of b and e c over its edges (c denotes the color of the edge e c ). The tetrahedral graph in D = 3 corresponds to δ t abcd T a T b T c T d . We call b the bubbles. We are interested in the partition function: where ρ b ≥ 0 are scalings chosen such that the large N limit of W exists. Observe that, contrary to [44], we allow the bubbles b to have several connected components. This is for instance the case of the double trace interaction bubble δ d ab;cd T a T b T c T d . Somewhat abusively, we some times call a bubble with several connected components a "multi-trace" interaction.
The generating function W is a sum over connect Feynman graphs G which have a new color 0 for the Wick contractions (propagators). As the propagators represent pairings of tensors, they connect vertices and G is a (D + 1)-colored graph. Denoting G0 the graph obtained from G by erasing the edges of color 0 we have: Due to the disconnected bubbles (multi trace interactions), the notion of connectivity in equation (C.3) subtle, hence the notation "i-connected" in the sum. The graph G is i-connected if any two interaction bubbles are joined by a path of edges of color 0 such that any two consecutive edge in the path are incident to the same interaction bubble. However, the graph G can be disconnected as a colored graph, C(G) > 1, because the edges in this path can be incident to different connected components in the bubbles. An example of an i-connected graph G which has C(G) > 1 is a double trace interaction decorated by two tadpole edges δ d ab;cd T a T b T c T d . It is a standard result [7,22] that the total number of faces of a (D + 1)-colored graph G is: where π runs over the D! jackets of G (that is the embedding of G corresponding to cycles over the colors) and k(π) is the non orientable genus of the jacket π. The non negative half integerω(G) is the degree of G. The degree of a disconnected graph is the sum of the degrees of its connect components. The bubbles b have only D colors therefore: The crucial property of the degree is that for any (D + 1)-colored graph G: This is a bit subtle. As G has D + 1 colors, G0 has only D colors. There is a D to 1 correspondence between the jackets π of G and the jackets π0 of G0 consisting in deleting the edges of color 0 in the jacket. As the non orientable genus can not increase by deleting edges we have k(π) ≥ k(π0) and consequently π k(π) ≥ D π0 k(π0). A (D + 1)-colored graph G has at most b∈G0 C(b) connected components. If G is i-connected, then it posses a tree of edges of color 0 connecting all the bubbles. Each edge in this tree joins two connected components on two different bubbles, hence decrease the maximal number of connect components of G by 1. Overall we get an upper bound on the number of connected components of G: Among the invariants (i.e. bubbles), an interesting subclass consists in the maximally single trace (MST) ones. They are those bubbles with only one face for each couple of colors. They are obviously connected and have exactly D(D − 1)/2 faces hence maximal possible degree: at fixed number of vertices.
D=3. Let us fix the ideas for D = 3. The bubbles b are 3 colored graphs. As such they are embedded graphs (ribbon graph, combinatorial map) with is the non orientable genus of b. Every b admits two jackets, (123) and (132), which are identical up to orientation and have non orientable genus k(b). The degree of b is its non orientable genusω(b) = k(b). The MST invariants have three faces and non orientable genus k(b) = −1 + V (b)/2. For instance the tetrahedron is MST and has non orientable genus 1. The wheel sextic interaction [78] is also MST and has non orientable genus 2. The Feynman graphs G have 4 color and 6 jackets, which are in 3 → 1 correspondence with the jackets of G0: The 1/N series.
We now chose to scale all the invariants by the "optimal scaling" introduced in [9,44]: With this optimal scaling Eq. (C.3) becomes: , (C. 12) which, due to the inequalities (C.6) (C.7), is a series in 1/N indexed by: The optimal scaling leads to a good large N limit. For some classes of interaction bubbles (like the MST or the melonic ones) the optimal scaling is the minimal scaling which still leads to a large N limit. It should be stressed however that this is not true in general: finding the minimal ρ b which still leads to a large N limit for an arbitrary interaction is a difficult open question [79].  The φ 4 1 operator is neither MST nor MMT, as the pillow is a connected invariant with ρ b = 1/2. In order to get a general idea of how it contributes to n-point functions we will consider first a simplified with only pillow operators of a single type, and then we will explicitly compute the two-point function of φ 4 1 at leading order.
Let us first consider a general correlator of n pillow operators with an arbitrary number of perturbative pillow vertices. We restrict to a single type of pillow operator, with single lines of color one: O p = φ a 1 b 1 c 1 φ a 2 b 1 c 1 φ a 1 b 2 c 2 φ a 2 b 2 c 2 . The intermediate field representation, known also as Hubbard-Stratonovich transformation, amounts to replacing it in the path integral by the integral over an auxiliary real symmetric N × N matrix field (the intermediate field), with ultralocal free covariance proportional to N (from the scaling of the pillow in the original action) and which couples to the composite matrix φ a 1 b 1 c 1 φ a 2 b 1 c 1 (see for example [15,80]). The original field appears then only quadratically in the new action, and thus it forms V q loop-vertices of valency q, for q ≥ 1, each containing two faces of the tensor model. Denoting by E the number of intermediate field propagators, and by F the number of faces that the intermediate field forms, we thus have that the connected n-point function of pillows scales as N q≥1 (2− 3 2 q)Vq+F +E−n = N q≥1 (1− q 2 )Vq+2−2g−n , (D.1) where the factor n is due to the fact that the inserted operators, unlike the perturbative vertices, carry no factor N . The amplitude would therefore be dominated by an intermediate field graph which is planar and which maximizes the number of univalent loop-vertices, that is, a usual cactus diagram. However, assuming that univalent loop-vertices (tadpoles in the original representation) have zero amplitude, we are left with dominant graphs being made of two-valent loop-vertices, joined in a planar way. Their amplitude scales like N 2−n . Comparing with Eq. (6.11), this means that such dominant graphs have ω = 1. Therefore, we conclude that also pillows at criticality have two-point functions of order N 0 , and higher-point functions suppressed in 1/N .

E Comparison to the Long-range Ising Model
We consider the [φ 4 ][φ 2 ][φ 2 ] correlation for N = 1 case, that is in the long-range Ising model. This model has a Wilson Fisher like point for > 0. We will explore why this correlation exhibits just anomalous scaling in this case, while in the large N limit of our long-range model it exhibits a logarithmic scaling. Let us first compute the correlation φ 4 φ 2 φ 2 for bare operators in the long-range Ising model [52]. This is the sum of four classes of terms, presented in Fig. 17, which yield: Figure 17: Diagrams contributing to φ 4 φ 2 φ 2 up to the first order of the coupling. φ 4 (x)φ 2 (y)φ 2 (z) = a 1 C(x − y) 2 C(x − z) 2 + a 2 g C(x − y) 2 d d u C(x − u) 2 C(z − u) 2 + (y ↔ z) where we introduced auxiliary parameters a 1 -a 3 to distinguish each diagram contribution, except for the third diagram which is order O( ). Expanding for small , we obtain: φ 4 (x)φ 2 (y)φ 2 (z) = c(ζ) 4 a 1 + Qg 4a 2 +a 3 + (2a 2 + a 3 ) log |x − y||x − z| − a 3 log |y − z| + · · · |x − y| 4∆ φ |x − z| 4∆ φ .
In order to facilitate the comparison between this case and our model, let us introduce a fictitious distinction between the classical contribution to the renormalized φ 4 operator and its quantum correction in the Ising model. That is we write: which holds in both cases. For the long-range Ising model the two φ 4 s in this equation are in fact one and the same φ 4 t = φ 4 2 = φ 4 , b = 3 and N = 1. The renormalized quartic operator in our model has the same form, with φ 4 t and φ 4 2 distinguished, N large and b = −6. In both cases, at first order in g we find φ 4 (x) φ 2 (y) φ 2 (z) = 1 + 2Q g φ 4 t (x)φ 2 (y)φ 2 (z) + bN −3/2 Q g φ 4 2 (x)φ 2 (y)φ 2 (z) , (E. 7) and adding up all the contributions (taking into account the scaling with N for our model), we obtain both for our model and for the long-range Ising model: c(ζ) 4 |x − y| 4∆ φ |x − z| 4∆ φ bN −3/2 Q g a 1 + 1 + 2Q g N −2 a 1 + Qg 4N −5/2 a 2 + N −3/2 a 3 + (2a 2 N −5/2 + N −3/2 a 3 ) log |x − y||x − z| − N −3/2 a 3 log |y − z| + · · · , (E.8) where the combinatorial factors a 1 , a 2 , a 3 and b are different in the two cases, and N = 1 for the long-range Ising model. In this form it is quite transparent why this three-point function has such drastically different behaviors in the two cases. Both for N = 1 and for N → ∞ the pole in 1/ cancels, and we get, up to overall factors: N −2 a 1 + Qg (2a 2 N −5/2 + N −3/2 a 3 ) log |x − y||x − z| − N −3/2 a 3 log |y − z| + · · · . (E.9) For N = 1 all the terms are of the same order of magnitude, and at small g the logarithms give anomalous dimensions of order g. However, for N → ∞ the constant term is suppressed an only the the logarithmic terms, scaling like N −3/2 , survive.