Analytic and Numerical Bootstrap of CFTs with $O(m)\times O(n)$ Global Symmetry in 3D

Motivated by applications to critical phenomena and open theoretical questions, we study conformal field theories with $O(m)\times O(n)$ global symmetry in $d=3$ spacetime dimensions. We use both analytic and numerical bootstrap techniques. Using the analytic bootstrap, we calculate anomalous dimensions and OPE coefficients as power series in $\varepsilon=4-d$ and in $1/n$, with a method that generalizes to arbitrary global symmetry. Whenever comparison is possible, our results agree with earlier results obtained with diagrammatic methods in the literature. Using the numerical bootstrap, we obtain a wide variety of operator dimension bounds, and we find several islands (isolated allowed regions) in parameter space for $O(2)\times O(n)$ theories for various values of $n$. Some of these islands can be attributed to fixed points predicted by perturbative methods like the $\varepsilon$ and large-$n$ expansions, while others appear to arise due to fixed points that have been claimed to exist in resummations of perturbative beta functions.


Introduction
The Landau theory of phase transitions [1] has provided a powerful framework for the search of emergent critical behavior for decades. It has served as a solid foundation for the development of renormalization group (RG) methods like the ε expansion [2,3], Monte Carlo [4] and functional RG [5,6]. These methods have been very successful in predicting critical behavior in a wide variety of situations, but there is still a surprising number of discrepancies and disagreements in the literature, pointing to potentially deep underlying physical principles. More recently, following pioneering work of [7][8][9] and especially [10], old conformal bootstrap ideas have morphed into an an entirely new and computationally rigorous approach for the study of conformal field theories (CFTs)-for a review see [11] and for an introduction [12].
In this work we use the conformal bootstrap method, both analytically and numerically, to  [13,14], and indeed such transitions have been claimed to be observed in various experiments. Our work is also of pure theoretical interest, since it attempts to shed light on the possible existence of fixed points that arise due to resummations of perturbative beta functions. The existence of such fixed points has been questioned by some functional RG studies [15][16][17]. Thus, theoretically the state of affairs regarding these fixed points has remained murky despite decades of effort, with conflicting results obtained by different RG methods.
The most unequivocal results for O(m) × O(n) CFTs have been obtained by taking m small and n large. The existence of a well-defined large-n expansion was established in [18], and the strongest results to date have appeared in [19][20][21]. The ε expansion has been widely used as well, see e.g. [18,22,13,14,23,24], with the highest-loop study (six loops) performed recently in [25]. Here we will use the analytic bootstrap method of large spin perturbation theory, introduced in [26,27] and developed further in [28][29][30], to confirm existing results in the literature and obtain some new large-n ones. The analytic bootstrap gives the same type of results as diagrammatic methods, but simplifies the computation of certain quantities, such as scaling dimensions of spinning operators, OPE coefficients and central charges.

Analytic bootstrap
The analytic bootstrap relies on the fact that conformal four-point correlators can be computed from their double-discontinuities, up to potential contributions from operators with spin zero or one.
The double-discontinuity measures the singularities that arise in the lightcone limit, corresponding to operators approaching pairwise null separation in Lorentzian signature, and is sensitive to operators of large spin in the operator product expansion (OPE). Spinning conformal primary operators group into twist families, where the scaling dimensions and OPE coefficients, collectively the CFT-data, are given by functions analytic in spin extracted from the double-discontinuity. All operators in a twist family have approximately equal value of the twist, defined as the difference between scaling dimension and spin: τ = ∆ − . In [31,32] it was shown that such twist families must exist in any CFT in dimension d > 2, and that the CFT-data is sourced by operators appearing in the OPE decomposition of the crossed channel. Specifically, the identity operator 1 in the φ four-point function gives rise to the leading order OPE coefficients for double-twist operators φ∂ φ with τ → 2∆ φ . Other crossed-channel operators induce corrections to the CFT-data in the twist families.
Large spin perturbation theory [26,27] combined with the Lorentzian inversion formula [33] constitutes a systematic framework for analytic bootstrap for theories with a small expansion parameter. It applies to both week coupling and strong coupling expansions, as well as to expansions in inverse number of degrees of freedom. At each order in the expansion, the whole double-discontinuity can be generated from an ansatz of contributions from a small set of crossedchannel operators, and the undetermined constants of this ansatz can later be fixed by consistency conditions, for instance conservation of symmetry currents. The method applies to a wide range of theories, and in particular it has been used to study the ε expansion for the Wilson-Fisher fixed point [28] and in the large-N expansion for the O(N ) model [30]. 1 In this paper we show how to generalize these implementations to critical φ 4 theories with any global symmetry group.
Consider a field φ transforming in some (vector) representation V of the global symmetry.
For the case of O(m) × O(n) we will take φ = φ ar in the bifundamental representation, where a = 1, . . . , m and r = 1, . . . , n. Operators and twist families in the OPE decomposition of the φ four-point function will transform in all irreducible representations R in the tensor product V ⊗ V . Looking first at the ε expansion, the identity 1 and bilinear operators φ 2 R will source the complete CFT-data of double-twist families in all representations to order ε 3 . Moreover, these bilinear operators are the spin zero operators of the same twist families. Despite the fact that spin zero is beyond the formal validity of the inversion formula, it was shown in [28] that it is possible to analytically continue to the formula for the scaling dimensions in the twist families to include the scalar operators. Encouraged by this observation, we conjecture that the same is true for φ 4 theories with any global symmetry, which leads to a set of quadratic equations for ∆ φ 2 R . By solving these equations we can determine all perturbative fixed points in the ε expansion consistent with the given global symmetry. 2 For the φ 4 theories with O(N ) symmetry, it is well-known that the large-N limit admits a 1 The critical O(N ) model has also been studied from the bootstrap perspective by the methods of multiplet recombination [34,35] and Mellin space bootstrap [36]. We briefly revisit the Mellin space bootstrap in Appendix D. 2 Apart from O(m) × O(n) symmetry considered here, we have checked this reproduces all known fixed points also in the case of MN [37], hypercubic and hypertetrahedral symmetry [38]. description in terms of a Hubbard-Stratonovich transform. In this description, the bilinear singlet operator φ 2 S gets replaced by an auxiliary field σ of approximate dimension 2. At criticality, the large-N expansion and the ε expansion are compatible-using the perturbative results one can for instance confirm that ∆ φ 2 S → 2 + O(1/N ). For certain global symmetry groups, the Hubbard-Stratonovich transformation can be generalized, where N corresponds to a specific group parameter. 3 More precisely, if certain bilinear operators φ 2 R have dimensions approaching 2, they should be promoted to auxiliary fields R. In large spin perturbation theory, these auxiliary fields will, together with 1, source the generalized 1/N expansion, and play the role same role as σ in the treatment of the O(N ) model in [30].
For the symmetry group of this paper, O(m) × O(n), we will expand in 1/n for fixed m.
From the results in the ε expansion, we note that two representations furnish scalars that can be promoted to auxiliary fields: S in the singlet representation, and W ab in the irrep that is a traceless symmetric tensor in O(m) and singlet in O(n). Like in the ε expansion, we can derive quadratic equations, and the solutions will determine all CFT-data at order 1/n. Apart from the free theory, we find that these equations generate all three non-trivial fixed points from the

Numerical bootstrap
The numerical implementation of our work involves the study of both single and mixed correlator bootstrap systems. In our study we exclude values of scaling dimensions of various operators that are not compatible with the combined requirements of unitarity and crossing symmetry; the latter is also known as associativity of the OPE. First, we probe the correlator φ ar φ bs φ ct φ du for self consistency, this is our single correlator system. Previous studies of O(2) × O(n) single correlator systems have appeared in [40] and [41]. We extend their results and match with analytic predictions from the large-n expansion.
For m = 2 and sufficiently large n (e.g. n 10), we find excellent agreement between the numerical bootstrap predictions and the analytic ones. This can be clearly seen in Figs. 2 and 3. operators. This reinforces the empirical notion that kinks in bootstrap bounds correspond to the position in parameter space of actual CFTs. We find that the antichiral fixed points appear as kinks in our W sector, which is a representation furnished by operators that transform in the two-index traceless symmetric representation of O(2) and the singlet representation of O(n). The chiral fixed points coincide with kinks in our X sector; operators in this representation transform as singlets of O(2) and two-index traceless symmetric tensors of O(n). As expected, for smaller values of n the agreement between large-n and numerical predictions becomes progressively worse.
For n = 2 we find a pronounced kink that appears to correspond to a known fixed point of the ε expansion as discussed in [37]. In the n = 3 case, there exist mild kinks that we study extensively with a mixed correlator system.
Bootstrap for the mixed correlator system for O(2) × O(n) CFTs is studied for the first time in this work. It consists of probing self consistency for four-point functions involving both φ and S, where S denotes the smallest dimension scalar operator (above 1) in the singlet representation.
Our goal is to obtain closed isolated regions (islands) in parameter space, which may correspond to physical CFTs. This method has so far produced extremely accurate calculations of critical exponents in the Ising and O(2) critical theories [42,43]. Islands have also been discovered in other theories with relevance to three dimensional statistical field theory [44][45][46][47]-for a more comprehensive list of references we refer the reader to [11].
We find two sets of islands, which we identify with two qualitatively different types of fixed points. The first set corresponds to the theories predicted by the large-n and ε expansions; these are found by saturating bounds in the W and X sectors as discussed in the previous paragraph.
We have found these islands for n as low as 6, which appears to be in agreement with the predictions of [25]-below n = 6 these fixed points are expected to be nonunitary. The second set of islands appears to correspond to fixed points that have been claimed to arise after resummations of perturbative beta functions [48][49][50]. These are not the same fixed points that are found in the standard large n and ε expansions, and their existence is not widely accepted [15][16][17]. In our bootstrap studies these islands are found by saturating bounds in the W and Z sectors, where operators in the Z sector transform in the antisymmetric representation of both the O(2) and the O(n). We find such islands for n = 3, which is an experimentally relevant value of n. chiral and collinear fixed points, we believe that further research is required to conclusively settle outstanding issues related to criticisms of some authors in the functional RG community [15][16][17].
Another issue that remains unresolved is related to the assertion of some authors that the O(2)×O(3) chiral and collinear fixed points are of the focus type and thus nonunitary [51,52,50,53].
We note that in the O(2) × O(3) case, the two islands we find are consistent with a second scalar singlet that has a scaling dimension above three, i.e. the corresponding fixed points are both stable in the context of the RG. This is something that cannot hold for fixed points found in the ε expansion [54,24,55]. We also note that all previously mentioned islands are obtained by making assumptions on the second B sector operator, which contains odd-spin operators among which the first spin-one operator is the conserved vector of the O(n) contained in O(2) × O(n).
The fact that the allowed region presents a sensitivity to assumptions specifically on the B sector was observed empirically. The dependence of bootstrap bounds on assumptions in sectors that contain conserved operators have been studied in other cases in e.g. [56][57][58][59].
The structure of this paper is as follows. In section 2 we review known perturbative results obtained in the ε and large n expansions. In section 3 we lay out the general formalism of the analytical bootstrap, applicable to any φ 4 theory. In section 4 we apply the formalism of section The mn scalar fields are arranged into a matrix, φ ar , with row indices running from 1 to m and column indices from 1 to n. The standard summation convention for repeated indices is used in where φ a are m vectors of size n each. The couplings u, v of (2.2) are related to the couplings λ, g of (2.1) by Fixed points with v < 0 are called collinear or sinusoidal and fixed points with v > 0 chiral, helical or noncollinear. With m n stability of the scalar potential requires u > 0 if v 0 and For every m there is a value of n, indicated by n + (m) above, for which the chiral and antichiral fixed points collide in the real u-v plane and subsequently move to the complex u-v plane. For n > n + (m) the chiral fixed point is stable, but for n < n + (m) there is no stable fixed point.
However, for some n − (m) < n + (m) two fixed points reappear in the u-v plane-this time they are called sinusoidal and antisinusoidal because they have v < 0, and the sinusoidal fixed point is stable. Furthermore, there is a value n H (m) < n − (m) below which the O(mn) fixed point is stable, for one of the fully interacting fixed points of the n H (m) < n < n − (m) regime crosses the v = 0 line and acquires v > 0 (chiral), while the other remains with v < 0 (sinusoidal). The situation is summarized in Fig. 1.
Numerical estimates for other values of m can be found in [25, Table 8]. In this work we will attempt to use our bootstrap bounds to independently estimate these quantities, particularly n + (2), and compare with (2.5). For m = 3 a similar study was performed in [40], while the existence of O(2) × O(n) theories in d = 5 was examined with perturbative methods in [60].
As mentioned in the introduction, it has been suggested that fixed points beyond the ones we just reviewed exist in O(2) × O(2) and O(2) × O(3) theories. Confusingly, the terminology "chiral", "collinear", etc., is still used for those fixed points, depending on their sign of the coupling v.
For scalar theories in the ε expansion below d = 4, it is a theorem that a stable fixed point, if it exists, is unique [54,24,55]. For the numerical studies in this work, we will fix m to a small value, specifically m = 2, and obtain bounds for increasing n. Thus, we expect that if kinks appear at large n, they will be due to fixed points of regime (I). In that case, we expect from the ε expansion that since the chiral fixed point is stable, the antichiral is unstable. This prediction is also expected to hold in the large-n limit in d = 3, to which we now turn.

Large n
As mentioned in the introduction, it was realized a long time ago that O m,n theories admit a large-n expansion [18]. For the chiral fixed point in d = 3, large-n computations give [19][20][21] These computations are done in the more general setting of arbitrary d at large n, but here we present the d = 3 results only. In [20], the operator C corresponds to our Z, and ηc computed there is given by ηc = ∆Z − 1. In [21], the operator T corresponds to our W , and χT computed there is given by χT = 3 − 2∆ φ − ∆W .
while for the antichiral fixed point in d = 3 at large n the results are [19,21] Here we denote singlet operators in the φ × φ OPE with the letter S, operators that transform as two-index traceless symmetric tensors under O(m) and singlets under O(n) by the letter W , and operators that transform as two-index antisymmetric tensors under both O(m) and O(n) with the letter Z. 6 As usual, primes denote the order in scaling dimension of these operators, i.e. S is the leading singlet, S the next-to-leading singlet and so on. We have not found large-n results for ∆ S − in the literature, but it is widely believed that ∆ S − < 3, i.e. the antichiral fixed point is unstable. By explicitly constructing singlet operators we find the ones in Table 1, where we tabulate the spectrum of the lowest dimension scalar singlet operators in the three fixed points at large n.
The results in (2.6) and (2.7) were obtained with the use of a Hubbard-Stratonovich transformation, extending a procedure used first in the O(N ) models by [61,62]. Two Hubbard-Stratonovich auxiliary fields, S and W ab , are introduced in this case. S is a singlet, while W ab transforms as a traceless symmetric tensor under O(m) and a singlet under O(n). The Lagrangian is [19,21] where repeated indices are summed over and w = u + (1 − 1/m)v with u, v the couplings in (2.2).
The equations of motion for S and W ab can be used to go back to (2.2). Below we will reproduce (2.6) and (2.7) and obtain more results at order 1/n following the analytic bootstrap logic of [30].

Analytic bootstrap for any global symmetry
In this section we outline an implementation of the analytic bootstrap that can be applied to φ 4 theories with any global symmetry. We begin with a brief review of large spin perturbation theory. This is followed by a summary of relevant results from the literature, in terms of a toolbox containing the explicit solution to the inversion problems we will encounter. We then give a general recipe for applying these tools, first to the ε expansion and then to the large-N expansion, where the latter is applicable to φ 4 theories which admit a Hubbard-Stratonovich description.
For the analytic bootstrap, we will consider the four-point function of φ i ∈ V , written in the In this expression, T ijkl R are tensor structures for the irreducible representations R ∈ V ⊗ V , 7 and u, v are the usual cross-ratios defined by The crossing equation follows from exchanging operators at x 2 and x 4 , and can we written as where the explicit form of the matrix M can be worked out from the tensor structures for the symmetry group under consideration. Here we choose normalizations in agreement with [30], so 7 More precisely, the T ijkl R are the projectors that can be used to decompose the four-point function into invariant subspaces labeled by R.
that the matrix in the O(N ) case takes the form where the representations are the singlet S, rank-two traceless symmetric T and antisymmetric A representations of O(N ). For any symmetry group, we reserve the letter S for the singlet representation.
Each of the functions G R (u, v) admits a decomposition in conformal blocks, where the sum runs over conformal primary operators O appearing with OPE coefficient c φφO in the OPE φ i × φ j | R , and the conformal blocks G ∆, (u, v) are functions which sum up the contributions of descendants to that primary.
The OPE expansion (3.5) is regular in the limit z,z → 0. However, we will expand in the lightcone limit, defined by z 1 −z 1. Taking z → 0, the conformal blocks as functions of z,z simplify as Lagrangian theories these operators have the schematic form O 1 ∂ µ 1 · · · ∂ µ n O 2 , up to contributions from descendants. In the theories we consider, φ is near the unitarity bound, ∆ φ = d−2 2 + γ φ , which means that the leading double-twist operators are weakly broken currents Our main objective is to determine the CFT-data of these operators, which consist of their scaling dimensions ∆ R, and their OPE coefficients a R, = c 2 φφJ R, . Of particular interest are the conserved currents: the stress-energy tensor T µν = J S,2 , and in the case of continuous global symmetry, Noether currents J µ R = J R,1 in one or several representations R. They have conserved dimensions, ∆ S,2 = d, ∆ R,1 = d − 1 and their OPE coefficients are related to central charges C T and C J R by the relations 8 following from conformal Ward identities [63,7].
Our main tool is the Lorentzian inversion formula, derived in [33]: where the double-discontinuity dDisc is defined as the difference between the correlator and its In particular, a single conformal block has vanishing double-discontinuity. The sign in (3.8) depends on the symmetry of G R (z,z) under exchanging x 1 and x 2 , and the normalization constants are given by The poles in ∆ of the function C R ( , ∆) are located at the scaling dimensions of the operators O R , with OPE coefficients given by the residue; more precisely The inversion formula is valid for > 1. If a power z τ /2 appears in dDisc[G(z,z)], it signals the existence of a family of operators of twist near that value. This follows from the scaling , which induces poles in C R ( , ∆) from the z → 0 limit of the z integral: . (3.10) We now focus on the leading twist family in each representation, and assume that τ R, = 2∆ φ + γ R, for small anomalous dimensions γ R, . In that case, following the manipulations of [33,Sec. 4], the z integral can be computed and (3.8) reduces to the one-dimensional inversion problem given in [64]: the CFT-data a R, and γ R, are given bŷ These expressions were derived in [64], by assuming that γ R, ∼ g 1 and expanding all quantities in g. This expansion generates terms proportional to z ∆ φ log p z in the double-discontinuity, which under the z integral are converted to higher order poles in ∆ in (3.10). These poles are responsible for the derivatives ∂h appearing in (3.11), following from changing variables from (∆, ) to (τ,h) in (3.9). An alternative heuristic derivation of (3.12) starting from the collinear conformal blocks is given in [28].
The success of large spin perturbation theory stems from the fact that dDisc[G R (z,z)] can be computed through the crossing equation (3.3). At each order in the expansion parameter, the whole double-discontinuity is computed from the conformal blocks of a small set of crossed-channel operators. In particular, the double-twist operators [φ, φ] R,n, themselves do not contribute at leading order. This is because the power (1 −z) −∆ φ from crossing is cancelled by (1 −z) τ n, 2 from the conformal blocks, as seen from the expansion (3.6) evaluated in the crossed channel, z → 1 −z.
Expanding τ R,n, = 2∆ φ + n + γ R,n, , we note that the first non-zero double-discontinuity arises from the term In the two cases relevant in this paper, the leading double-discontinuities will be generated by scalar operators, and the weakly broken currents J R, = [φ, φ] R,0, will only contribute at subleading order.

Inversion toolbox
In this section we collect all inversion formulas from the literature that are required to find the leading order CFT-data in the ε expansion and in the large-N expansion. The entries take the a statement referred to in the literature as reciprocity and proved in the context of CFT in [65].
Inversion 1. The identity operator 1 contributes with only. This holds for generic ∆ φ , and applies in both the ε expansion and in the large-N expansion. This result can be directly computed from the integral (3.12) using an integral representation for kh(z) (see e.g. eq. (4.7) of [33]).
This was derived in eq. (2.34) of [28] using the explicit form of the scalar conformal block as an infinite sum [7] and taking the small z limit. Here S 1 (h − 1) denotes the harmonic numbers.
Inversion 4. The leading contribution from an infinite sum over ∈ I ± of broken currents J with anomalous dimensions γ = κ J 2 and OPE coefficients αa GFF 0, | µ−1 (where a GFF n, | ∆ φ are the generalized free field OPE coefficients given in (B.2)) is Here E ± are lengthy expressions given explicitly in (B.3). This formula was derived in [30] by summing over blocks on the unitarity bound and subsequently inverting the sum.

General solution in the ε = 4 − d expansion
Consider first the contribution from the identity operator, appearing in the singlet (S) representation. This will give rise to the leading contribution to U (0) R,h in all representations. Since this is the only operator contributing until order ε 2 , we get, using Inversion 1, . From this expression, the leading order OPE coefficients can be extracted: Here takes even (odd) values for R being an even (odd) representation. The scalar bilinears φ 2 R in the even representations have OPE coefficients These scalars are the next operators to contribute to the double-discontinuity. Assume that they . Then, using Inversion 2 we get the order ε 2 corrections Using (3.11) we can thus write down the leading correction to the anomalous dimension, Next, as observed in [28], we assume that it is possible to analytically continue the result γ R (h) to spin zero, by making the change of variablesh →h f = ∆+ 2 , i.e. we replace the bare with the full conformal spin (eigenvalue of the quadratic Casimir). For spin zero we should evaluate . This leads to a system of equations R even , (3.27) at order ε, where now one power of ε in the γ R (h) cancels against the factorh f − 1 = (g R − 1)ε/2 in the denominator. This simplifies to which is a system of k quadratic equations for the k constants g R , where k is the number of even representations, or equivalently the number of scalar bilinears. Solving (3.28) gives all possible fixed points in the ε expansion with the given symmetry.
As an example, consider the O(N ) model with crossing matrix (3.4). The even representations are S and T , and the bilinear scalars are φ 2 S = φ i φ i and φ 2 T = φ {i φ j} . There are two solutions to (3.28), g S = g T = 0, which is the Gaussian theory, and which is exactly the critical O(N ) model [29]. With these values we have Γ The singlet spin-two current in any global symmetry group is the stress-energy tensor with . Using this we write down the full dimension of the broken currents to order ε 2 : The OPE coefficients are extracted using (3.11), From the spin two singlet OPE coefficient we can extract the central charge correction given by (3.7) Conservation of the stress-energy tensor allows one to compute the leading (order ε 2 ) anomalous dimension of φ.
In one or several of the odd representations R, the current at spin = 1 may be conserved, being a generator of global symmetry. This therefore gives further constraints ∆ R, = d − 1, which must be explicitly checked. By (3.7) the corresponding OPE coefficient is related to the C J of that symmetry current: Let us discuss the extension to higher orders in the ε expansion. To order ε 3 , the operators contributing with a nonzero double-discontinuity are the same as at the previous order, namely the bilinear scalars φ 2 R . At higher orders, infinite families of operators contribute. In the O(N ) model, the only such families at order ε 4 are operators of approximate twist 2 and 4, and subsequently the problem was solved there. We expect that this generalises to any global symmetry. However, to compute the contribution from approximate twist 4 requires detailed knowledge of the operator content of the theory in question. This was done in the case of the O(N ) model in [29]. To order ε 5 the same operators will contribute but now with subleading corrections. To work this out, even in the N = 1 (Ising) case, is still an open problem.
We have seen that at order ε 2 all constants that enter the problem can be fixed using continuation to spin zero and conservation of the stress-energy tensor. This is no longer true at higher orders. At order ε 3 a total of 2k + 1 new constants appear: γ R ε) + . . . , and the corrections α R to the OPE coefficients defined by Based on experience from the O(N ) model [29], the order ε 2 continuation to spin zero requires order ε 4 results for the currents, so the only new equations at order ε 3 are the conservation of the symmetry currents (including the stress-energy tensor). In general this will not provide enough equations to fix all constants, but in many cases we can still make progress. Firstly, we will use that the correction to the OPE coefficient takes the form α R = −g R . This holds for any global symmetry, and follows from analytic bootstrap in Mellin space, as we show in Appendix D.
Secondly, the second order corrections g (2) R to the bilinear scalar dimensions are in many cases known from the literature, and one can proceed using these as input.
Using Inversion 2, it is straightforward to derive expressions for U (p) R,h at order ε 3 . The anomalous dimensions extracted from these expressions take the form From the corresponding expression for the OPE coefficients using Inversion 2, we can extract the central charge correction: Here we used that the stress-energy tensor conservation eliminates the dependence on g R in favour of γ (3) φ . For two-coupling theories as considered in [23] we may find where λ * , g * are the coefficients of the order-ε values of the two couplings at the fixed point, i.e. Similarly, for the current central charges we derive the expression

General solution in the large N expansion
Let us now describe the computation of CFT-data in the large-N expansion for a generic symmetry group, parametrised by some number N . Compared to the ε expansion the situation is a bit more complicated, since the parameter N enters in the crossing matrix M R R itself. In a given even representation R, we assume that there are two possibilities for the smallest dimension scalar. It is either a scalar bilinear φ 2 R with dimension 2∆ φ + O(N −1 ), or a Hubbard-Stratonovich field R with dimension 2 + O(N −1 ). If one has access to results in the ε expansion, one can assess the situation by taking the large N limit of the order ε scalar dimensions. For instance, in the O(N ) model we get, using (3.29), so we see that the singlet representation admits a Hubbard-Stratonovich field S (in the literature denoted by σ), but not the traceless symmetric representation (whence we keep the notation φ 2 T ). We assume that the Hubbard-Stratonovich fields R have dimension ∆ R = 2 + O(N −1 ) and OPE . In order to provide some structure of the subsequent computations we define the following subsets of the representations in V ⊗ V = I ∪ II: • Group I: Representations whose only corrections at order 1/N come from crossed channel Hubbard-Stratonovich fields.
• Group II: Representations where the corrections at order 1/N come from Hubbard-Stratonovich fields as well as from broken currents in Group I representations in the crossed channel.
• Group III: Representations that admit a Hubbard-Stratonovich field. Typically III ⊂ II.
As an example, in the O(N ) model we have S ∈ II ∩ III and T, A ∈ I. Our strategy will then be the following. First, as in the ε expansion, the identity operator creates the leading contribution to U Hubbard-Stratonovich fields will give the order 1/N anomalous dimensions in these representations.
Using Inversion 3 we see that these corrections will be proportional to 1/J 2 . Finally we turn to representations in the Group II. Here we get contributions from both the Hubbard-Stratonovich fields, using Inversion 3, and from the currents in Group I. Due to the particular form of the anomalous dimensions of these currents, we can use Inversion 4 to find the complete order 1/N CFT-data.
The expressions will depend on |III| + 1 free parameters: the OPE coefficients a R = c 2 φφR for R ∈ III, and the leading order anomalous dimension of φ. The only consistency conditions available to fix these constants are the conservation of the symmetry currents (including the stress-energy tensor), and depending on the number of conserved currents this may or may not be enough.
. For the representations in Group I we get the contributions from Hubbard-Stratonovich fields in Group III. Using Inversion 3 we get and a corresponding expression for U R,h . From this we extract the order N −1 anomalous dimensions of currents in Group I representations: where the scaling dimensions are given by ∆ R, = 2∆ φ + + γ R, In step 3 we consider the second group of operators, II. They get contributions both from R for R ∈ III and from J R, for R ∈ I.
We get where the + (−) sign is used if the operators in the R representations have even (odd) spin. This means that the anomalous dimensions of the group II double-twist operators are In the above expressions J 2 = (µ − 1 + )(µ − 2 + ) and As an example, let us explicitly evaluate K R and K R in the O(N ) model. We get . (3.50) The extension to subleading orders in 1/N is a complicated task, which was achieved in [30] for the T and A representations. This involved computing the contributions from operators [σ, σ] n, ,which was found in the form of a Mellin space amplitude, using OPE coefficients derived from the mixed correlator φφσσ . We do not attempt to generalize it for generic global symmetry group.

Analytic bootstrap of O(m) × O(n) CFTs
In this section we will apply the methods

Results in the ε expansion
In the ε expansion, the operators φ 2 R in the five even representations R introduce corrections to weakly broken currents in all nine representations. Solving equations (3.28) for the constants g R we find four sets of solutions, corresponding to the free theory and to the O(mn) (Heisenberg), chiral and antichiral fixed points. For the latter two fixed points, of interest to this paper, the explicit expressions for the g R are rather complicated, containing square roots. We give complete results in an ancillary data file, which we describe in Appendix A. For presentation purposes the expressions in the ε expansion in this section are expanded for at large n, but at each order in ε presented here the complete function of m and n has been determined.
The constants g R correspond to the scaling dimensions of the scalar operators, which take the for the chiral fixed point, and for the antichiral fixed point. These results agree with [23, Eqs. (5.92) and (5.93)]. 9 The operators φ 2 S , φ 2 W and φ 2 Z correspond to S, W and Z, respectively, in (2.6) and (2.7). Having identified the fixed points we move on to a determination of the CFT-data to order ε 3 .
As described in the previous section, we need to take as input the second order corrections g (2) R to the anomalous dimensions γ φ 2 R of bilinear scalars, given in [23].
agrees with the literature values [19], whereas the results of the spinning operators are new, as far as we are aware. These include the dimensions of the non-conserved spin-one operators and the central charges More results can be found in the ancillary data file, as described in Appendix A.

Results at large n
As mentioned in the introduction, we can give a description at large n by introducing Hubbard-Stratonovich operators in the S and W representations. This is in agreement with the results (4.3) and (4.4) in the ε expansion above. As a starting point for the analytic bootstrap at large n we will therefore assume that these two representations contain scalar operators S and W, with dimensions and OPE coefficients given by We will take these representations to consitute group III in our implementation of the recipe of section 3.3. The next task is to determine what operators contribute at order 1/n to the broken currents in all of the nine reprensentations in (4.1). This is done by expanding the crossing matrix Having identified the groups I, II and III, we follow the implementation of section 3.3 and generate CFT-data at order 1/n in all representations. In particular, the scaling dimensions are given by ∆ R, = 2(µ − 1) + + 2γ for γ R, given by (3.45) and (3.47) and µ = d/2. These expressions depend on three undetermined constants, a S , a W and γ (1) φ . Fortunately, the stress-energy tensor and the two global symmetry currents J µ A and J µ B provide three consistency equations for these unknowns, namely There are four solutions to these equations, which we can identify with the four fixed points of the ε expansion, In Table 2 we summarize the twist families of the O(m) × O(n) symmetric theory at large n in the chiral and antichiral fixed points. We give the leading twist family in each representation, and we also display a couple of subleading families in the singlet representation. The existence of each of these subleading families follows from the initial analytic bootstrap considerations of [31,32], since they are the double-twist operators in a suitable four-point function. Importantly, these families contain more than one operator at each spin and therefore participate in mixing. In the cases where there is an operator at spin zero, we match it with the scalar singlets presented in Table 1.
In Table 2 we also explain how the scaling dimension of each twist family relates to the corresponding scalar. In similarity with the O(N ) model, we assume that the expressions (4.9) can be analytically continued to spin zero, giving the dimension or the shadow dimension of the corresponding scalar. Including also ∆ φ , this gives for the chiral fixed point and for the antichiral fixed point  The values for φ, S, W and Z agree with those quoted in section 2.2, whereas we are not aware of any previous results for the remaining operators. We also give results for the non-conserved spin one operators 14) The computation of the OPE coefficients provides results for the central charges, by (3.7). For the chiral fixed point we get (4.15) and for the antichiral fixed point where the precise form of the constants c i is given in (B.5) in Appendix B.

Single correlator
In the single-correlator bootstrap, for which the crossing equations are discussed in Appendix C, we have obtained bounds on the dimensions of the leading scalar operators in the representations S, W, X, Y, Z as functions of the dimension of φ. The most interesting results are obtained in the W and X plots. More specifically, using large-n results we can see that the (mild) kinks that appear in the X-bounds (see Fig. 2) are due to the chiral fixed points, while the kinks that appear in the W bounds (see Fig. 3) are due to the antichiral fixed points. While this is clear at large n only, we will assign the same meaning to the kinks at smaller n, but only down to n = 6 where, as we will see below, the situation becomes more subtle. As seen in Fig. 3, obtained with m = 2 for different values of n, there are very sharp kinks at large n that get smoothed out as n decreases.
In Fig. 2 the kinks are much milder. They certainly exist at large n, but they are not so clear at low n-see Fig. 4-even when at the same n, e.g. n = 4, there is a clear kink in Fig. 3.
At n = 6 we can see from The persistence of the kinks in the ∆ W bounds even at small n (n = 4, 5 and n = 3 although the kink is much milder for n = 3), combined with their absence in the ∆ X bounds, is rather puzzling. After all, intuition from the ε expansion dictates that the antichiral and chiral fixed points, to which we have attributed the kinks in the ∆ W and ∆ X bounds, respectively, annihilate and become complex fixed points at some n (in this case n + (2), whose estimated value in the ε expansion is 5.96(19) [25]). Since the bootstrap excludes nonunitary theories, both kinks are expected to disappear at some n around 6, and indeed this is borne out to some extent for the kinks in the ∆ X bounds, as we discussed in the previous paragraph.
One explanation for the persistence of the ∆ W kinks at small n is that our numerical bounds are insensitive to the putative nonunitarity of the antichiral fixed point for small n. This scenario could be further examined by estimating the size of this nonunitarity in perturbation theory, in a properly quantified sense that we do not discuss here, and comparing it with that of the chiral fixed point. We do not pursue this direction here, but it is worth investigating in the future.
Another possibility is that the kinks in the ∆ W bounds at small n are due to another fixed point, which is not the naive continuation of the antichiral fixed point to small n. Evidence for the existence of such a fixed point, belonging to the chiral universality class, exists in the literature; see [14,Sec. 11.5.3] and references therein. A further universality class, typically called collinear, is supposed to exist for small n. 10 These chiral and collinear fixed points arise after resummations of the perturbative beta functions. However, this approach has been criticized in [15], and the functional RG predicts that they do not exist [16,17]. In section 5.3 below we will see that, consistently with the conclusions of [41],

Mixed correlators
In the mixed correlator system, to find an island around the chiral fixed point we use the following (A-O 2,10 -5) dimension of next-to-leading bifundamental operator, φ , above 1.6, i.e. ∆ φ 1.6.
Let us note here that even with ∆ S 3 we obtain islands around antichiral fixed points, so long as we keep the gap on ∆ φ small. This is inconsistent with the fact that the antichiral fixed point is unstable, but with our numerical power (see Appendix C) we cannot see the inconsistency.
However, when we increase the gap on ∆ φ sufficiently, we do see that the antichiral island disappears with ∆ S 3. This is presumably due to the crossing equations that arise from the φφSS four-point function, which in the 12 → 34 channel are sensitive to both φ and S, while in the 14 → 32 channel they are sensitive to φ but not S. As we have already mentioned, the ε expansion predicts that n + (2) = 5.96 (19), meaning that a unitary chiral fixed point may exist for n = 6. The state-of-the-art analysis of the O 2,6 theory with the ε expansion was performed recently in [25]. It turns out that we can also find an island with our nonperturbative numerical bootstrap methods, and so we can compare our results with large-n and ε expansion results. This is done in Fig. 9.  Table 9]. This island is obtained with the assumptions (C-1)-(C-5) using saturation of the O 2,6 X bound in Fig. 5.

Mixed correlators in the O(2) × O(3) case
The case m = 2, n = 3 is of particular interest since it is believed to appear as a symmetry in the continuum limit of frustrated spin models at criticality. We remind that the fixed points for m = 2, n = 3 arise after resummations of the perturbative beta functions, i.e. they are not found in the standard perturbative ε expansion. There are two fixed points, called chiral and collinear, with the chiral being of relevance for the experimentally observed phase transitions in stacked triangular antiferromagnets [14,Sec. 11.5], and the collinear potentially relevant for the normal-to-planar superfluid transition in 3 He [52]. Monte Carlo simulations support the existence of these fixed points, but functional RG methods come to a different conclusion, namely that these fixed points do not exist and that experiments in frustrated magnets are actually seeing weakly first-order phase transitions [15][16][17].
In what follows we present results regarding putative theories that live on the W and Z sector single correlator bounds. The W (resp. Z) sector bound has a mild kink that appears to correspond to the chiral (resp. collinear) fixed point. This observation was first made in [41]. Here we take the next logical step and perform a mixed correlator bootstrap around these kinks. We also compare to theoretical predictions and experimental data where applicable. We note that the general trend of sensitivity to assumptions in the B sector persists here as well.
Let us start by noticing that, however mild, there seems to be a kink on the W sector single correlator plot of Fig. 6 around ∆ φ = 0.539. The kink can be seen more clearly in Fig. 10. As mentioned above, evidence for the existence of the O 2,3 chiral fixed point has appeared before in the literature, based on resummations of perturbative beta functions. Such computations have been performed at six loops using the massive zero momentum (MZM) scheme [51] and at five loops using the modified minimal subtraction (MS) scheme [49] (see also [48]). Monte Carlo simulations have also shown signs of a critical theory, with the most recent analysis performed in [53]. According to results of [50, Eq. (2.9)] and [48, Table III], the theory at the chiral fixed point has ∆ φ = 0.545 (20) and ∆ W = 1.79 (9) in the MS scheme. 11 In the MZM scheme, [51, Table   III] and [48, Table III] give ∆ φ = 0.55(5) and ∆ W = 1.91 (5). The MS scheme result is more consistent with the location of the kink in Fig. 10. (C-O 2,3 -4) dimension of the next-to-leading scalar singlet, S , above 2, i.e. ∆ S 2.
(C-O 2,3 -5) dimension of next-to-leading bifundamental operator, φ , above 1.5, i.e. ∆ φ 1.5, then we obtain Fig. 11. Let us note that the island in Fig. 11 remains even if we make the more constraining assumption ∆ S 3, which is compatible with the RG stability of the O 2,3 chiral fixed point. If instead of (C-O 2,3 -3) we assume that J µ B can appear with dimension below 2.4, then both the island and the peninsula are part of a bigger, continuous peninsula that includes  Even with the large error bars in (5.1), we see that agreement is best with the results of [50], mainly due to ∆ S . 12 In conjunction with the ∆ W results mentioned earlier, it is clear that the MS results of [50] and [48] agree best with our bootstrap results for the chiral fixed point.
Experimental results for transitions described by the O 2,3 chiral fixed point can be found in [14, Table 37]. The agreement of our results for the critical exponent ν = 1/(3 − ∆ S ) is very good, although the same cannot be said for the critical exponent β = ∆ φ /(3 − ∆ S ).
Refs. [51,50,53] indicate that their fixed points are of the focus type, which means that they are nonunitary, while our study pertains to unitary theories. It is unclear to us how sizable nonunitarities of the type discussed in [51,50,53] could have been missed by our bootstrap results.
We note that the numerical bootstrap has previously found islands for theories that are believed to be nonunitary, namely the five-dimensional O(N ) models [66] and the Ising model in d = 4 − ε [67].
In the former case, increasing the constraining power of the numerics led to the disappearance of the allowed region. It is possible that also our island in Fig. 11 will disappear with stronger numerics.
The O 2,3 collinear fixed point corresponds to a kink in the bound of the first scalar operators in the Z irrep; see Fig. 12. According to results of [52] and [48,Table III], the theory at the collinear fixed point has ∆ φ = 0.543 (12) and ∆ Z = 1.8(1) in the MS scheme. 13 In the MZM scheme, [52] and [48, Table III] give ∆ φ = 0.5395 (35) and ∆ Z = 1.75 (10). The consistency of the MZM scheme result with the location of the kink in Fig. 12   (Coll-3) dimension of next-to-leading vector operator in the B sector, J µ B , above 2.5, i.e. ∆ J µ B 2.5, (Coll-4) dimension of next-to-leading scalar singlet, S , above 3, i.e. ∆ S 3.
(Coll-5) dimension of next-to-leading bifundamental operator, φ , slightly above ∆ φ , i.e. ∆ φ ∆ φ + 0.01, we find a rather large island that appears to end slightly to the left of the position of the kink in island and peninsula, but rather a continuous peninsula that gets narrow in the region between the island and peninsula of Fig. 13. According to [52], in the MS scheme the O 2,3 collinear fixed point has ∆ φ = 0.543 (12) and ∆ S = 1.41 (20), while in the MZM scheme it has ∆ φ = 0.5395 (35) and ∆ S = 1.31 (11). Here the MZM scheme result for ∆ φ and the MS scheme result for ∆ S appear to agree better with our island in Fig. 13

Single correlator in the O(2) × O(2) case
In the O 2,2 case, the collinear fixed point is equivalent to the O(2) fixed point [48], and so we will only be examining the proposed chiral fixed point [51,48]. As discussed in Appendix C, the single-correlator crossing equations of the O 2,2 = O(2) 2 /Z 2 theory are equivalent to those of the MN 2,2 = O(2) 2 S 2 theory discussed in [37]. A strong kink was obtained in the X sector of [37] (see Fig. 1 there), which corresponds to the Z sector in (C.6) below. This bound is shown in Fig. 14. According to [51, (5), ∆ S = 1.25 (9), ∆ WX = 1.75(4), ∆ Y = 0.97 (7) and ∆ Z = 0.46 (12), while, in the MS scheme, [50,Eq. (2.8)] and [48, Table III] give ∆ φ = 0.545 (20), ∆ S = 1.46 (14), ∆ WX = 1.66 (15), ∆ Y = 1.00 (15) and ∆ Z = 0.63 (15). 14 Therefore, we conclude that the kink in Fig. 14 does not correspond to the O 2,2 chiral fixed point. In [37] it was suggested that this kink may correspond to the fully-interacting theory of the ε expansion analyzed in [68][69][70][71]. To continue our search for the O 2,2 chiral fixed point, we obtain a bound on the dimension of the first scalar operator in the WX representation. The bound is shown in Fig. 15. An extremely mild change in the slope of the bound is observed at ∆ φ = 0.547 (2), at which point ∆ WX = 1.507 (10). These numbers are in relatively good agreement with the MS scheme results mentioned in the previous paragraph. The mildness of the kink in Fig. 15, however, indicates that we need stronger numerics in order to reach solid conclusions. We leave this for future work. also obtained an island (see Fig. 9). In all these cases we used a mixed correlator bootstrap. In the m = 2, n = 6 case we were able to compare our nonperturbative bootstrap results with the state-of-the-art resummed ε expansion results of [25]. The agreement is reasonably good, but our island is rather large. We need stronger numerics and more refined bootstrap techniques in order to make our island smaller and obtain more accurate determinations of critical exponents.  [48][49][50] and Monte Carlo computations [53] indicate the presence of fixed points beyond the ones found with perturbative methods, while the functional RG concludes that such fixed points do not exist [15][16][17].

Summary and conclusion
Correspondingly, there are two conflicting suggestions, namely that experiments are seeing critical (second-order) or weakly first-order behavior. It is important to note here that the beta function resummation and Monte Carlo fixed points are suggested to be nonunitary (focus type), with next-to-leading scalar singlets of complex scaling dimensions. We attempted to address these issues in this work. Our results provided support for the existence of these fixed points, but we saw no signs of nonunitarity. Overall, we were unable to provide conclusive answers, but we believe that more dedicated bootstrap work with stronger numerics will be able to reach definitive conclusions in the near future. In the ancillary data file of the arXiv submission, we give full results for a number of quantities, presented on the form quantity fx-pt expansion , where fx-pt ranges over the fixed points ON, Chiral and Anti; and expansion over the expansions Eps and N. In Table 3 we list the options for quantity .

Acknowledgments
The expressions make use of the constants implemented as sqRmnval and eta1val respectively (see (3.50)). The values that we present for deltaScalars fx-pt Eps contain results from the literature at order ε 2 , extracted from [23].  Data for ope f x-pt Eps include long expressions that we do not include in the ancillary data file but can be made available upon request.

Appendix B. Explicit formulas used in the main text
Here we give a few explicit formulas used in the main text. In (3.32) we have which is an expansion of the general formula for the OPE coefficients of generalized free fields [72], where (a) n = Γ(a+n) Γ(a) is the Pochhammer symbol. In Inversion 4, E ± is given by and in expressions (4.15) and (4.16) we have where the subscript N is introduced to capture the relation δ ij δ ij = N . The four-point function φ ar φ bs φ ct φ du can be decomposed as where the sum over R runs over the representations S, W, X, Y, Z, A, B, C, D, For the crossing equation we define and we find 16 In our numerical bootstrap considerations we define conformal blocks using the conventions of [67]. 16 In (C.5) we suppress the labeling of the F ∆, 's and c 2 O 's with the appropriate index I. The appropriate labeling, however, is obvious from the overall sum in each term.
The signs that appear as superscripts in the various irrep symbols indicate the spins of the operators we sum over in the corresponding term: even when positive and odd when negative.
When m = n there is a reduction in the number of crossing equations. In this case, instead of separate projectors P W arbsctdu and P X arbsctdu , we only have the projector P WX arbsctdu = P W arbsctdu +P X arbsctdu , and similarly P AB arbsctdu = P A arbsctdu +P B arbsctdu and P CD arbsctdu = P C arbsctdu +P D arbsctdu , always with m = n. The crossing equation in the m = n case is antiferromagnets, helimagnets and structural phase transitions has been discussed extensively in [37] and references therein.
In this paper we also consider a mixed correlator bootstrap involving the four-point functions φφφφ , φφSS and SSSS , where S is the first scalar singlet in the φ × φ OPE. The crossing equations for this system are straightforward to derive due to the fact that S transforms in the singlet representation, and so in the OPE φ × S we find only operators that transform in φ's representation under the global symmetry.
As usual, our numerical treatment involves two steps, namely generation of an xml file that encodes the problem and subsequently its solution with a numerical algorithm. For the first step we use PyCFTBoot [67], and for the second SDPB [73].

Appendix D. Mellin bootstrap for any global symmetry
In this appendix we will revisit some equations from the analytic boostrap in Mellin space and apply them to the ε expansion for φ 4 theories with arbitrary global symmetry. The framework is described in detail in [74] and was generalized to global symmetry in [36], focusing on O(N ) and hypercubic symmetry. We will show how the Mellin space bootstrap reproduces the system of equations (3.28), which we used to find the perturbative fixed points for a given global symmetry.
In addition, we will show that valid for all t, where we have generalized the notation of [74] by adding global symmetry indices.
From (D.2), one derives the constraint equations, by projecting to a given spin in the s-channel of the R representation, using the fact that the M ∆, (s, t) can be expressed in terms of continuous Hahn polynomials. The first term becomes simply ∆ c R ∆, , but for the other two terms, operators of any spin contribute. We arrive at constraint equations of the form ∆ c R ∆, q R(2,s) ∆, where, again, the exact form of the involved quantities can be found in [74]. We have collected the t and u channel contributions under the label t, and the appearance of the matrix M R R follows from using the projectors of (3.  [36]. We will evaluate them for = 0, assuming that ∆ φ 2 and c 2 φφφ 2 R = a R,0 (1 + α R ε) + O(ε 2 ), and expanding to order ε. Only = 0 contributes to the sum, and following [74] we substitute In [36], the CFT-data was computed to order ε 3 for O(N ) and hypercubic symmetry, and we believe that this will generalize to arbitrary global symmetry. From such implementation, one can derive the order ε 3 correction to the OPE coefficient c 2 as input), an observable that is inaccessible from large spin perturbation theory in its present formulation.