The Cubic Fixed Point at Large $N$

By considering the renormalization group flow between $N$ coupled Ising models in the UV and the cubic fixed point in the IR, we study the large $N$ behavior of the cubic fixed points in three dimensions. We derive a diagrammatic expansion for the $1/N$ corrections to correlation functions. Leading large $N$ corrections to conformal dimensions at the cubic fixed point are then evaluated using numeric conformal bootstrap data for the 3d Ising model.


Introduction
Theories often simplify when there are a large number of fields. For theories with vector-like large N limits, the infinite N behavior can often be solved exactly, and finite corrections computed in a systematic 1/N expansion. The most famous such class of theories are the O(N ) vector models, as reviewed for instance in [1], which have been studied in great detail using both traditional Feynman diagram methods and analytic bootstrap techniques [2][3][4][5].
Other theories which can be solved at large N include scalar O(M )×O(N ) and S M O(N ) M models [6,7], the Gross-Neveu-Yukawa model [8] and its supersymmetric generalizations [9], conformal QED [10], and the Abelian Higgs model [11]. In each case, a renormalization group flow exists between a free field theory in the UV and an interacting fixed point in the IR.
The flow becomes tractable at infinite N , and 1/N corrections to correlators in the IR can be computed systematically. Similar flows exist in holographic theories when the boundary is perturbed by a "double-trace" deformation, which changes the boundary conditions of the bulk fields [12,13].
Techniques used to solve vector-like large N limits have rarely been applied to cases where the UV theory is interacting. At strictly infinite N it has long been known that the IR theory is almost identical to the UV theory, up to a shift in the conformal dimension of the perturbing operator [14][15][16]. Computing finite N corrections to this limit, however, requires a quantitative understanding of the initial UV theory, and beyond d = 2 this has proven challenging to acquire. Over the last decade, the numeric conformal bootstrap has emerged as a practical tool for non-perturbative computations in interacting conformal field theories [17,18]. For the Ising model in d = 3, the numeric bootstrap provides the most precise determinations of both conformal dimensions and OPE coefficients [19][20][21][22]: ∆ σ = 0.5181489 (10) , ∆ = 1.412625(10) , f σσ = 1.0518537 (41) , f = 1.532435 (19) , (1.1) where the error bars are rigorous bounds. Further, though less rigorous, results are available for many other conformal dimensions, OPE coefficients, and even some four-point functions [23,24]. Computations competitive with the best Monte Carlo results for critical exponents have been recently achieved for the O(2) and O(3) models [22,[25][26][27], while also yielding OPE coefficients which are otherwise difficult to determine using traditional methods.
Our aim in this paper is to use the 3d Ising model bootstrap data to study the cubic fixed points at large N . We will do so by coupling N Ising models via the interactions where i (x) is lowest dimension Z 2 -even operator in i th copy of the 3d Ising model. These are the only relevant interactions preserving the S N Z N 2 hypercubic symmetry of the N decoupled Ising models. As we flow into the IR while suitably fine-tuning µ, the theory flows to a conformal field theory known as the cubic fixed point [16,28]. Our strategy is similar to that in [29,30], in which conformal perturbation theory was combined with bootstrap data in order to study the random-bond Ising model (which is the N → 0 limit of the cubic fixed point) and the long-ranged Ising model respectively. Unlike conformal perturbation theory, in which a controlled approximation is possible only given an exactly marginal operator, the large N limit naturally gives us a controlled approximation for the cubic fixed points without this onerous requirement.
Recall that the cubic model is described within the 4 − expansion by the Lagrangian where we use · 0 to denote expectation values in the decoupled theory. The Hubbard-Stratonovich trick allows us to rewriting the cubic coupling using an auxiliary field: Note that to properly make sense of this equation we must regulate the φ(x) kinetic term, a subtlety we discuss in Appendix A. Applying (2.2) to our generating functional and then redefining φ(x) → φ(x) + µ, we find that Here we have used the fact that each Ising model is decoupled when λ = µ = 0 to factor the correlator. Each expectation value in (2.3) is then the generating functional for correlators of (x) and O(x) in a single Ising model: It will prove more convenient, however, to work with the free-energy which is the generator of connected correlators. We therefore write The operators i (x) form a reducible representation of the symmetric group S N , splitting into the two representations We will find it useful to rewrite the sources J i (x) to match this decomposition, defining so that the generating function becomes Performing the change of variables φ → φ − 1 N J , we find (2.10) The quadratic term in J (x) does not contribute to correlation functions at separated points, and so we can ignore it. We then see that correlators of φ(x) are related to those of E(x) by the relationship and so for simplicity we shall redefine our path integral so that J (x) is the source term for φ(x) rather than E(x), defining: We can now take λ → ∞ while keepingμ = µ λ fixed, and so arrive at the expression Nowμ is the only dimensionful parameter left. When it is suitably fine-tuned to criticality we can compute correlators of E(x), i (x), and O i (x) at the cubic fixed point.
Nothing we have done so far has required N to be large. In order to actually compute (2.12) however, we will need to consider the large N limit. For simplicity, let us first set J i = K i = 0 and so focus only on correlators of φ(x). To compute these we expand the Ising model free-energy as a sum of connected correlation functions Because the Ising model is conformally invariant, the first few connected correlators take a simple form: while the four-point function can ne computed numerically using conformal bootstrap results.
We can thus expand (2. 16) At large N the path-integral integral is dominated by the quadratic term, so that as N → ∞ the correlators of φ(x) are given by those of a mean-field theory. For this theory to be conformally invariant we must set the dimensionful parameterμ to zero. Computing the φ(x) two-point function is then straightforward (see for instance [34]): As expected for a double-trace deformation, φ(x) has conformal dimension d − ∆ [14][15][16].
This should be compared to E(x) in the UV theory, which has dimension ∆ .
Large N corrections are generated by the higher degree terms in (2.16). We can compute them using a Feynman diagram expansion. For each k ≥ 3 we have a k-point vertex, whose value is given by the k-point connected correlator in the Ising model: Each of these vertices corresponds to one term of the infinite sum in the second line of (2.16); we can think of each as a non-local φ k interaction. In a Feynman diagram each leg of a vertex is either connected to another leg, or it is connected to an external operator. Internal legs are connected to other legs via the propagator and we then integrate over both z 1 and z 2 . Legs connected to external operators take a similar form: 20) but this time we only integrate over z, but not x. Finally, we can always connected two external operators to each other: (2.21) To compute an n-point correlator φ(x 1 ) . . . φ(x n ) , we sum over all Feynman diagrams which have an external leg φ(x i ) for each i = 1 , . . . , n. As in any Feynman diagram expansion, we must divide by the symmetry factor; this purely combinatorial factor is identical to that found in a local theory with φ k interactions. Counting the overall order of a diagram is straightforward: each vertex contributes a power of N , while each propagator gives a power of N −1 .
As an example, the leading corrections to the two-point function are given by: where we have dropped any disconnected and tadpole diagram which may occur at O(N −1 ).
To illustrate our rules let us write the first diagram in (2.22) explicitly in terms of Ising model correlators. Rewriting the diagram to make the internal spacetime points explicit, we see that: We will evaluate this integral in the next section, along with the other diagram (2.22).
Together, they allow us to compute the leading correction to the anomalous dimension of φ at large N .
Let us next turn to correlators of O i (x). For simplicity we begin with the 1pt function, which we can compute by differentiating the generating functional: We can then expand the derivative of the free energy as a sum of connected correlators and so find that We can then write the first few terms in (2.26) at large N as (2.28) In practice we can always fix O i (x) = 0 through a field redefinition, and so can automatically subtract of these tadpole diagrams.
Next consider the higher point correlator If all of the i k 's are distinct, each derivative acts on a different copy of the Ising model and (2.30) As an example, the first few contributions to the two-point function If multiple O i (x) operators have the same index i, additional terms will appear in (2.29) where a single Ising model free energy is differentiated multiple times. For instance, We have already seen how to compute the first term. To compute the second, we note that To represent these, we introduce vertices with two external lines, along with k internal φ(z) lines, whose value is the connected correlator O( Because of the Kronecker delta δ ij , this vertex only appears if i = j, so that both Os belong to the same copy of the Ising model. With this addition, we can write the first few corrections : . . . Consider now more general correlators . is a different kind of local operator in the 3d Ising model (again excluding the (x) operator, which we will discuss shortly). So for instance, we might want to consider correlators including both the σ i (x) operators and the stress-tensor T µν i (x) operators. Generalizing our diagrams to this case is straightforward: we simply introduce vertices with more general external legs . . .
We compute correlators by summing over all diagrams with the appropriate external legs.
Finally, let us turn to correlators of i (x), which we can compute by differentiating with respect to J i . Even though these satisfy the constraint i J i = 0, we can treat the sources so that, up to contact terms, the operator i i (x) automatically vanishes. We can therefore treat the i (x) akin to other local operators, introducing external vertices such as . . . In this paper we will compute leading corrections to anomalous dimension, and so we will only need to compute one-loop diagrams. Because of this, we can subtract off divergences in an ad hoc fashion. Higher-loop calculations however require a more careful treatment. A common strategy used to study the large N limit of the O(N ) models is to introduce a longranged kinetic term for the auxiliary field [35][36][37], a technique we expect should generalize without difficulty to the case considered here.

Anomalous Dimensions
Having derive expressions for cubic correlators at large N expansion, we will now evaluate the leading corrections to conformal dimensions of operators at the cubic fixed point.

Single-trace operators
The σ i (x) operators form a single irreducible representation of the S N Z N 2 hypercubic symmetry group. When suitably normalized, their two-point function is where ∆ * σ (N ) is the conformal dimension of σ i (x) at the cubic fixed point. Using (2.35) and noting that σ i (x)σ j (x) automatically vanishes due to hypercubic symmetry, we find that and so are left with a single diagram to compute: The integral diverges as both z 1 and z 2 approach either x or 0. To regulate it, we cut out little spheres of size Λ around both x and the origin when performing the z 1 integral. We then note that under conformal inversion, Performing the change of variables we thus find that where we integrate over the region C for which z 1 satisfies As we take Λ → 0, however, we can replace the latter condition with the simpler one which is equivalent to the original condition up to O(Λ) corrections. We then perform a second change of variables Now observe that the integral is rotationally invariant, so that we can set z 1 = rê 1 without loss of generality, whereê 1 is the unit vector (1, 0, . . . , 0). Performing a final change of variables z 2 → rw and then using scale invariance, we can write is the surface area of the unit sphere in d dimensions. Defining and subtracting off the logarithmic divergence which appears as Λ → 0, we arrive at the Comparing this correction to the exact result (3.1) for the two-point function, we conclude All that remains for us to compute is ζ σ . Integrals of this form were considered in [29,30], where they arise in calculations of β-functions and anomalous dimensions in conformal perturbation theory. Our integral (3.11) suffers from power-law divergences when w →ê, which we must subtract. As shown in [29], the integral (3.11) can be rewritten as where R = {w : |w| < 1 , |w| < |w −ê|}. Only the second term in this expression has a power-law divergence, which appears as w → 0 and must be subtracted. We now expand σσ as a sum of conformal blocks: where z = w 1 + i w · w − w 2 1 , and where g ∆ O , O (z,z) are the conformal blocks computed as in [21]. The low-lying conformal dimensions and OPE coefficients contributing to these correlators have been computed using the numeric conformal bootstrap [23], and so we can sum the operators with ∆ < ∆ * = 12 in order to approximate these four-point functions.
The truncation error is approximately bound by [30,38]: (3.16) We then apply the method described in Appendix C of [29] to numerically evaluate (3.14) block by block. We find that ζ σ ≈ −6.40 ± 0.02 (3.17) for the 3d Ising model, where the error is dominated by the uncertainty in the OPE coefficients computed in [23], rather than from the truncation bound (3.16) or from the numerical error in the integration. However, the error estimate should be taken with a grain of salt, as it assumes both that the OPE uncertainties are independent and that no operators are missing from [23] with ∆ < 12. Plugging this into (3.13), we conclude that The leading large N correction to ∆ * σ is much smaller than the typical error estimate for ∆ * σ when computed using either the 4 − or d = 3 perturbative expansions, and so we cannot usefully compare our results to these calculations. 1 The only value of ∆ * σ we have a values, computed using various resummations of the 4 − and d = 3 expansions, can be found in [39]. For N = ∞, the estimate with the least uncertainty is ∆ * σ,∞ = 0.516(1), which both considerably smaller and much more uncertain then the bootstrap result ∆ * σ,∞ = ∆ σ = 0.5181489 (10). Finite N estimates suffer from similar issues.
precise handle on is when N = 3. For this special case we expect the conformal dimension to be close to that of the O(3) model, with estimated from the d = 3 perturbative expansion [40]. Combining this with the current best result for ∆ σ, O(3) = 0.518920 (25), computed using Monte Carlo [41], we conclude that ∆ * σ, N =3 cubic = 0.51887 (6) .

E(x) and i (x)
Our next task is to compute the anomalous dimension of i (x). To this end, we expand The first diagram is an vertex with only two external 1 (x) legs, which is simply the two-point function of (x) in the Ising model, The second diagram is given by where we have introduced δ as a regulator, evaluated the integrals by using . (3.25) twice, and then taken δ → 0. The result is finite, and so does not contribute to the anomalous dimension of i .
The final diagram for us to consider is: This diagram is almost identical to (3.3) except that the external legs are now operators rather than σ operators. We evaluate it in an identical fashion, finding that (3.28) Just like ζ σ , we compute ζ by expanding as a sum of conformal blocks using [23], and then integrate the blocks numerically using [29], finding that Combining everything together, and so we conclude that Now let us turn to E(x), which at the cubic fixed point becomes the operator φ(x).
We already wrote the leading corrections to φ(x)φ(0) conn in (2.22), and now must simply compute the two diagrams. We begin with the first diagram in (2.22): (z 24 z 35 |x − z 1 ||z 6 |) 2d−2∆ (z 12 z 13 z 23 z 45 z 46 z 56 ) ∆ . (3.32) To evaluate the integrals, we make repeated use of the star-triangle relation which holds when a 1 + a 2 + a 3 = d. We first integrate z 1 and z 6 using (3.33), then next integrate z 2 and z 5 , before finally integrating z 3 . Once the dust settles, we find that We can perform our final integral over z 4 using (3.25), however, when we do so we find it diverges. To extract the finite contribution we use dimensional regularization: Subtracting the divergent term, we conclude where we define (3.37) Now we turn to the second diagram in (2.22): To evaluate this, we note that this diagram is just the shadow transform of (3.26), and so We can now compute the integrals over z 1 and z 2 using (3.25), and so find that (3.40) Having evaluated both diagrams in (2.22), we can combine them to find that: We thus conclude that at large N , the dimension of (3.42) In Table 1, we compare our large N expression for ∆ * E to calculations performed in both the 4 − expansion [42] and the d = 3 fixed dimension expansion [43,44]. These expansions give series which are only asymptotically convergent and so require resummation. As we can  This value is a little larger than the estimates from the perturbative expansions, and a little smaller than the estimates from the large N expansion.

Stress Tensor
We now turn to the stress tensor. In the Ising model, or indeed, any local CFT, the stress tensor satisfies the Ward identity [47] while the two-point function is given by (3.46) With this choice of normalization, the quantity c T , sometimes called the "central charge" of the theory, is defined such that c T = 1 for a free massless scalar. Numeric bootstrap has been used to compute c T ≈ 0.946534 (11) for the 3d Ising model [20].
Each of the N copies of the Ising models contributes an operator T µν i (x), but once we couple the models together only the sum of these stress tensors remains conserved. The other N − 1 linearly-independent operators will acquire an anomalous dimension, so that ∂T µν i (x) ∝ V ν eats some vector multiplet V ν . As we shall see, this multiplet recombination allows us to compute the anomalous dimension of T µν i without computing the TT two-point function directly. Our calculation is similar to one performed in [30], in which the anomalous dimension of the stress tensor was computed in the presence of a long-ranged interaction.
We begin by using conformal invariance to fix 49) and so find that (3.50) But we know that because ∂ µ T µν i (x) vanishes at infinite N , in order for it to be non-zero at finite N it must eat some independent conformal primary V ν (x) of dimension 4. The only such operator available is and so we conclude that for some function b(N ). To fix b(N ), we first compute that (3.53) But at large N we can alternatively calculate ∂ µ T µν i (x)V λ (0) diagrammatically. We begin by computing 54) where in passing from the second to the third line we have made use of the Ward identity (3.45). We now use (3.54) to compute: (3.55) Comparing this identity to (3.53), we deduce that We now use the multiplet recombination to compute (3.57) and by comparing to the general result (3.50), conclude that Specializing to the case d = 3, we find that

Multi-trace Operators
Having computed the leading corrections to the scaling dimensions of all single trace operators in the cubic fixed point with dimension ∆ ≤ 3, we now will move to consider multi-trace operators.
Let us begin by considering the double-trace operators σ i σ j for i = j, which together form an irreducible representation of the S N Z N 2 hypercubic symmetry. To compute their anomalous dimension, we expand where we have introduced δ as a regulator. Using (3.33) to evaluate the z 1 integral and then (3.25) to evaluate the z 2 integral, we find that (3.61) Combining this with the expression for σ 1 (x)σ 1 (0) conn derived in Section 3.13, we find that and so conclude that

(3.63)
We should note that the splitting between single and double-trace operators under a doubletrace deformation was previously computed in [48], by first considering the leading 1/N corrections to four-point functions and then extracting double-trace anomalous dimension by decomposing these four-point functions into conformal blocks. Their results applied to the cubic fixed point imply that which matches with our own calculation. Now let us consider the more general tensor σ i 1 . . . σ i k (x) for k < N , where each i m = i n for any m = n. Such tensors again from an irreducible representation of S N Z N 2 . It is not hard to see that at leading order in the 1/N expansion, the only diagrams contributing to σ i 1 . . . σ i k (x)σ i 1 . . . σ i k (0) conn are products of (3.3) or (3.60) with disconnected σσ two-point functions, so that (3.65) We therefore conclude that

(3.66)
Next we turn to the tensor i j (x), calculating .

(3.67)
This leading correction is almost identical to the first diagram in (3.60), differing only in the OPE coefficient and power of |x| out front: We thus conclude that

(3.69)
Once again, we can easily extend the calculation to the more general tensor i 1 . . . in , finding

Discussion
In this work we studied the large N limit of the cubic fixed point. We derived a diagrammatic expansion for the 1/N corrections to correlators at the fixed point. Combining this with results from the 3d Ising model numeric bootstrap, we were able to compute leading corrections to the anomalous dimensions of a number of operators.
Our results for ∆ * E could be compared to perturbative calculations and found to be in good agreement when N ≥ 4. Unfortunately, the other conformal dimensions considered in this paper have to our knowledge never been accurately computed before, 3 whether that be through perturbation theory, Monte Carlo simulation, or numeric conformal bootstrap.
Recent numeric conformal bootstrap studies have considered theories with hypercubic symmetry [50][51][52][53], but these have so far been unable to isolate the cubic fixed point. It would be very interesting to try to use the large N results derived here to determine reasonable assumptions that could be fed into the conformal bootstrap, hopefully allowing the cubic fixed point to be bootstrapped and hence our large N results to be tested.
While we focused in this paper on the cubic fixed points, our methods are general and can apply to any theory with a vector-like large N limit. The most straightforward generalization is to the M N models [54], which have S N O(M ) N symmetry. These should be related to 3 As noted in Section 3.1.1, ∆ σ has been studied in both the 4 − and d = 3 perturbative expansions, but the uncertainties are large relative to the size of the 1/N correction to ∆ * σ derived in this paper. A number of other conformal dimensions have been computed at to 2-loop in the 4 − expansion [49], but this does not suffice to reliably extrapolate to = 1. It is straightforward to generalize our results to theories with double-trace deformations of the form where O I (x) are operators in a conformal field theory which are charged under some global symmetry G. For instance, the 3d O(M ) models for M ≥ 2 contain a symmetric tensor t ij (x) whose conformal dimension ∆ t < 3 2 . We can use this operator to construct a number of relevant double-trace deformations which break the S N O(M ) N symmetry of N decoupled O(M ) models down to discrete subgroups, and so can construct a large class of conformal field theories with tractable large N limits.

A Regulating the Hubbard-Stratonovich Trick
In this appendix we will discuss proper meaning of (2.1). Naively, we would expect that however, this last expression is problematic as i (x) 2 is not a well-defined operators. To circumvent this, we modify the path-integral of φ(x) so that We then write Using the OPE expansion, we can then evaluate We can ignore descendants, as these becomes total derivatives upon integration over x. In the limit → 0, terms for which ∆ A > 2∆ will go to zero, whereas if ∆ A < 2∆ , the coefficient diverge and need to be subtracted off. This means that we only have to worry about operators satisfying ∆ A = 2∆ , but of course no such operator exists in the 3d Ising model. We thus conclude that (2.1) is true so long as we regulate and subtract off divergences.