The Gross-Neveu-Yukawa Archipelago

We perform a bootstrap analysis of a mixed system of four-point functions of bosonic and fermionic operators in parity-preserving 3d CFTs with O(N) global symmetry. Our results provide rigorous bounds on the scaling dimensions of the O(N)-symmetric Gross-Neveu-Yukawa (GNY) fixed points, constraining these theories to live in isolated islands in the space of CFT data. We focus on the cases N = 1, 2, 4, 8, which have applications to phase transitions in condensed matter systems, and compare our bounds to previous analytical and numerical results.


Introduction
The conformal bootstrap [1][2][3] has emerged as a powerful tool to rigorously 1 constrain CFT data of strongly-coupled fixed points. Based solely on unitarity, symmetry, and assumptions about gaps in the spectrum of scaling dimensions, the bootstrap has produced stringent bounds on critical exponents of several universality classes describing real-world statistical and quantum phase transitions. These include the 3d critical Ising model [4][5][6][7][8][9] which describes liquid-vapor transitions and uniaxial magnets, and the O(N ) models [8,10,11] -in particular the O(2) model [12] which describes the superfluid transition in 4 He, and the O(3) model [13] which describes classical Heisenberg ferromagnets. The Ising and O(N ) models are perhaps the simplest 3d universality classes that can be reached via a renormalization group (RG) flow from a scalar theory. In this work, we focus on perhaps the simplest 3d universality class involving fermions: the O(N )-symmetric Gross-Neveu-Yukawa (GNY) model. The GNY model contains N Majorana fermions ψ i transforming in the vector representation of O(N ), interacting with an O(N )-singlet pseudoscalar ϕ [14]. The Lagrangian is 2 This theory has a critical value of m 2 below which ϕ spontaneously gets a nonzero vacuum expectation value (VEV), giving a mass to ψ i and, consequently, breaking parity. Above the critical value of m 2 , the VEV of ϕ vanishes, parity is preserved, and the fermions are massless. At the critical value of m 2 , this theory is expected to flow to a CFT with a single relevant parity-even O(N )-singlet scalar operator ϵ ∼ ϕ 2 .
Beyond serving as one of the simplest models of scalar-fermion interactions in quantum field theory, the GNY universality classes have been proposed to describe a variety of quantum phase transitions in condensed matter systems with emergent Lorentz symmetry. For example, this model and some of its variations have been proposed to describe phase transitions in graphene (using N = 8) [16][17][18][19], time-reversal symmetry breaking in d-wave superconductors (also for N = 8) [20,21], and time-reversal-symmetry breaking transition of edge-modes in topological superconductors (for several low values of N ) [22]. For the special case N = 1, the GNY critical point is expected to exhibit emergent supersymmetry at the transition [23].
The GNY models have been studied previously with the conformal bootstrap in [15,24]. Those works performed a bootstrap analysis of a single four-point correlator of fermionic 1 Throughout this paper, when stating that our bootstrap results are rigorous we mean that they do not rely on any unstated assumptions about QFTs. However, the bounds are not completely rigorous in the mathematical sense since they rely on some technical assumptions about our search algorithms (for example, over OPE space). 2 Here we are following the conventions in appendix A of [15] and contracting the indices of the two components of the Majorana fermions by ψiψi = Ω αβ ψi,αψ i,β . Note that ψiψi is parity-odd -hence the Yukawa term ϕψiψi preserves parity, since ϕ is a pesudoscalar.
operators ⟨ψ i ψ j ψ k ψ l ⟩. 3 The resulting bounds exhibited a sequence of kinks on the boundary of the space of allowed CFT data, which showed good agreement with perturbative estimates of the scaling dimensions of the GNY critical points (such as the ϵ [19,[25][26][27][28][29][30][31] and large-N expansions [15,29,[32][33][34][35][36][37][38]). However, bootstrapping the four-fermion correlator was not enough to constrain the GNY theories to lie within isolated islands. In this paper we expand this study to include a mixed system of scalar and fermionic operators, characterized by their representations under parity and the global O(N ) symmetry. We focus on the critical points with N = 1, 2, 4, and 8, which capture the experimentallyrelevant transitions described above and are outside the perturbative control of the large-N expansion. The four-point functions we analyze include all combinations allowed by the symmetries of ψ i (the lowest dimension fermion in the vector representation of O(N )), σ ∼ ϕ (the lowest dimension, parity-odd, O(N ) singlet scalar), and ϵ ∼ ϕ 2 (the lowest dimension, parity-even, O(N ) singlet scalar). We refer to σ, ψ i , and ϵ as "external operators" (as opposed to the "internal operators" that appear in their OPEs). By simultaneously imposing crossing symmetry and unitarity for all four-point functions of external operators, we show that each GNY critical point lies within an isolated island that severely constrains its low-lying scaling dimensions.   [31] (represented by the orange boxes), and the location of the N = 1 island from the N = 1 supersymmetric Ising bootstrap [42][43][44] (represented by the x).
Our bounds on the scaling dimensions of σ, ψ i , and ϵ for the theories with N = 2, 4, and 8 are illustrated in figures 1 and 2, and we give a numerical summary in table 1. Results for the N = 1 theory using similar methods are shown later in figure 8. For all values of N , our results are close to perturbative estimates from resummations of the ϵ-expansion. For N = 8, the large-N estimates are also close to the island that we find. Monte Carlo estimates for some of the scaling dimensions are also available for N = 2, 4, 8 [39,41,45]; these results, while close to our islands, have error bars that, in most cases, are disallowed by the rigorous bounds obtained from the bootstrap (see section 4). In all cases, our work presents a significant jump in the precision of scaling dimension determinations.
An important subtlety is that there are in fact two different GNY models that are often confused in the literature, having the same number of fermions but different global symmetry   [42][43][44] (represented by the x). We also superpose the general σ-ϵ bootstrap bounds with the assumption ∆ σ ′ > 3 from [43].
groups. In addition to the O(N ) GNY models discussed above, there are also "chiral" GNY models that possess an O(N/2) 2 ⋊ Z 2 global symmetry. These models, sometimes referred to as being in the "chiral Ising" universality class, are nearly degenerate with the O(N ) GNY models for the most common low-lying operators, being only distinguishable at high perturbative order. The Monte Carlo estimates mentioned above at N = 4, 8 are believed to fall in this class. However, due to the expected near-degeneracy, we posit that our bootstrap results for the leading operators also provide good (albeit non-rigorous) estimates of the scaling dimensions in the O(N/2) 2 ⋊ Z 2 GNY models. In this work we review some of the existing perturbative estimates for scaling dimensions in both models and provide a number of new ones that will be useful in our bootstrap study. We will also do our best to differentiate which of the two critical models is known to describe various phase transitions discussed in the condensed matter literature. Another interesting case is the N = 1 GNY model which, as previously mentioned, is believed to have emergent supersymmetry at criticality. The island that we find with mild assumptions about gaps in various sectors of this model, shown later in figure 8, is fully consistent with this picture. Without making any assumptions about supersymmetry in the bootstrap setup, we find: (i) The island lies right on a line along which the low-lying scaling dimensions are related due to supersymmetry, ∆ ϵ = ∆ ψ + 1 2 = ∆ σ + 1. (ii) The lowest dimension operator with spin-3/2 is very close to the unitarity bound across the entire island, suggesting the existence of (at least) an approximate supercurrent for any theory that lies within it. (iii) The N = 1 super-Ising CFT, whose scaling dimensions have been tightly constrained in previous bootstrap studies by a priori assuming supersymmetry, can also be seen to live right at a tip of the island. These computations provide a nice consistency check of our implementation and suggest that there are no non-supersymmetric fixed points at N = 1 (consistent with our gap assumptions).
Our work is a culmination of a series of developments in the numerical bootstrap that extend practical limits to allow studies of a wider set of 3d CFTs. To set up the crossing equations, we used the formalisms developed in [46,47] and the conformal block computation algorithms developed and implemented in blocks 3d [47][48][49][50]. We also implemented several Haskell libraries (described in appendix A) to efficiently and robustly set up systems of mixed correlators and compute the resulting bounds. We explored the space of external scaling dimensions and external OPE coefficients using the Delaunay search and "cutting surface" OPE search algorithms developed in [12]. Finally, we solved large-scale SDPs using the solver SDPB [7,51]. We hope that these technologies can be used to place comparable constraints on more complicated 3d CFTs with fermionic degrees of freedom, including extensions of the GNY models with different global symmetries or Chern-Simons matter theories whose monopole operators carry half-integer spin.
The structure of this paper is as follows. In section 2, we give theoretical background, including details about the perturbative expansions of the GNY models. We also discuss the gap assumptions that we impose as well as the differences between the O(N ) GNY models and the O(N/2) 2 ⋊ Z 2 GNY models. In section 3, we discuss the numerical setup of our bootstrap problem; the reader interested solely in results for the GNY model can continue straight to the next section. There, in section 4 we discuss the main results of the paper and show a series of bounds on scaling dimensions of the low-lying operators in the theory. We discuss possible future directions in section 5.

Large-N and ϵ-expansion
Since we need a baseline expectation for the scaling dimensions of the external operators σ, ψ i , ϵ, and we must also impose gaps for some low-lying internal operators, we collect the leading estimates for scaling dimensions of the O(N ) GNY model obtained from the large-N Operator Parity O(N ) ∆ at large N . ∆ in ϵ-exp. Table 2: Large-N [15,29,[32][33][34][35][36][37][38] and ϵ-expansion estimates [19,[25][26][27][28][29][30][31] for the scaling dimensions in the GNY model. The results for ∆ ψ ′ and ∆ χ are new and derived in appendix B. The numerators denoted by # indicate known expressions that have been omitted here for concision. ∆ ψ has a positive correction at O(1/N 3 ) that can be found in [34] and ∆ ϕ k has a negative correction (for k > 1) at O(1/N 2 ) that can be found in [37].
expansion and the ϵ-expansion in table 2. In addition to the known perturbative results, we give new calculations of the leading correction at large-N for ∆ ψ ′ i and ∆ χ i in appendix B. Some of these scaling dimensions have been computed to higher order, primarily in the context of the closely related O(N/2) 2 ⋊ Z 2 GNY models (discussed further below). However, the results for the leading scaling dimensions ∆ ψ , ∆ σ , ∆ ϵ are degenerate between the O(N ) theory and the O(N/2) 2 ⋊Z 2 theory up to 3-loop order. For large-N estimates see [15,29,[32][33][34][35][36][37][38], for ϵ-expansion estimates see [19,[25][26][27][28][29][30][31], or see [29] for a two-sided Pade expansion in the 2 + ϵ and 4 − ϵ expansion. We will primarily compare our results to the 4-loop ϵ-expansion resummations performed in [31], done using a computation scheme that is believed to be applicable to the O(N ) GNY models.

The N = 1 theory and emergent supersymmetry
We also get additional information for our gap assumptions from the N = 1 case which is of interest by itself. This fixed point is expected to exhibit emergent supersymmetry in the IR and, while this has not been rigorously checked, there have been numerous perturbative tests of this proposal. The basic argument relies on the N = 1 critical point having a single relevant singlet scalar in the IR. That is, we expect that only a single coupling -the mass m -needs to be tuned in (1.1) to m = m * (g, λ) in order to reach the O(1) GNY IR fixed point, and the same fixed point is reached regardless of the values of g and λ. On the other hand, for N = 1 and λ = g 2 /2, the interaction in (1.1) is explicitly N = 1 supersymmetric. After tuning the scalar mass to the critical point the RG flow will preserve N = 1 supersymmetry. 4 Since, as all other critical RG flows of (1.1), it terminates at the O(1) GNY CFT, this implies that this CFT has N = 1 supersymmetry. This argument has been more concretely probed by observing the expected supermultiplet relations ∆ ϵ = ∆ ψ + 1/2 = ∆ σ + 1 between the ϵ-expansion results for ∆ σ , ∆ ψ , and ∆ ϵ obtained at four loop order [30]. These relations were also probed using a 2-sided Padé approximation in the 2 + ϵ and 4 − ϵ expansions [29].
The emergence of supersymmetry was also probed non-perturbatively with the conformal bootstrap in [24], where, without explicitly imposing supersymmety, a kink was seen close to the line relating ∆ ψ = ∆ σ + 1/2 when imposing the appropriate bound for ∆ σ ′ . Nevertheless, since the CFT associated to the O(1) GNY critical point has no theoretical guarantee to live at the kink and since the location of the kink in [24] changed heavily with the imposed gap on ∆ σ ′ , it is interesting to try to constrain the fixed point to lie within an island through which the line ∆ ϵ = ∆ ψ + 1/2 = ∆ σ + 1 would pass. In this paper we will indeed find such a bound by making some mild assumptions about the gaps in the various sectors of the theory, providing additional evidence for supersymmetry in N = 1 case, see section 4 for details.
By using assuming N = 1 supersymmetry for several scaling dimensions and OPE coefficients and using a mixed system of scalar operators [42][43][44] found highly accurate estimates for the CFT data in the N = 1 super-Ising critical theory, the fixed point to which it was suggested that the O(1) GNY model flows in the IR. For instance, the scaling dimension of ∆ σ was found to be [44], ∆ σ = 0.5844435(83), from which the scaling dimensions of ∆ ϵ and ∆ ψ can also be found. By imposing relations between scaling dimensions within the same supermultiplet, [44] also found accurate estimates for the scaling dimensions of ∆ ψ ′ = 3.3869(25), ∆ χ = 4.88, ∆ ϵ ′ = 3.8869 (25), and ∆ σ ′ = 2.8869 (25). We will use these results to inform our gap assumptions for other values of N as we describe below. We now discuss what gaps we can reasonably assume in the spectrum when trying to find islands for the GNY fixed points. Given the abundance of evidence in favor of supersymmetry for N = 1 it will be convenient to assume it when studying the theories with other values of N . Specifically, we can combine the highly accurate superconformal bootstrap results for the scaling dimensions of N = 1 theory with the results from the large-N expansion in a two-sided Padé approximation as shown in figures 3 and 4. This allows us to better estimate the scaling dimensions of the various low-lying operators in the theory and to better assess what gap assumptions we can make in the various sectors of the theory for the other values of N that we study. We want to stress that while these approximations are by no means rigorous, we only use them to motivate our gap assumptions. Given the gap assumptions, which we state explicitly below, our bounds are rigorous.
The plots in figures 3 and 4, together with the ϵ-expansion and other large-N results, lead us to make the following gap assumptions: Large-N   • Since ϵ, σ and ψ i will serve as the external operators in our bootstrap search, we will   not make any assumption about their scaling dimensions.
• We assume that ϵ is the only relevant neutral scalar operator in the theory. Consequently, we will assume ∆ ϵ ′ > 3. This assumption is substantiated by the Padé approximation obtained in the bottom-right plot in figure 3 as well as by ϵ-expansion estimates.
• Motivated by the large-N equation of motion which removes the operator ψ i ψ i from the spectrum of primaries, it will also be useful to impose bounds on ∆ σ ′ . The Padé approximation for this scaling dimension is shown in the bottom-left plot of figure 3. It suggests that this operator is always irrelevant for N ≥ 2. For most of our plots we will make the assumption of irrelevance ∆ σ ′ > 3, but for N = 2 we will also study the more conservative assumption ∆ σ ′ > 2.5.
• For operators that are not part of the spectrum for N = 1, such as σ T , we can no longer rely on the two-sided Padé approximations and will instead use the estimates from the large-N expansion as well as the past bootstrap results [24], which showed a kink in the (∆ ψ , ∆ σ T ) at the expected location for the GNY critical point. There it was found that ∆ σ T > 2 for all N which, emboldened by the large-N estimates, we will take to be our gap value in this sector.
• We will also assume gaps for the low-lying fermionic operators. In particular, due to the equation of motion / ∂ψ i = −gϕψ i , 5 χ i has a larger scaling dimension than could be naïvely expected. We will conservatively assume that ∆ χ > 3.5. This assumption is well within the expectations from the two-sided Padé interpolation that is shown in the lower-right plot of figure 4, which in fact estimates that ∆ χ > 4 for all N ≥ 1.
• We will also assume ∆ ψ ′ > 2, based on the two-sided Padé approximation shown in the lower-left plot of figure 4.
• We also assume in all cases a small twist gap of 10 −6 to improve numerical stability of the semidefinite programming algorithm. That is, we assume that all operators except the identity, the stress tensor, and the conserved O(N ) current have twist at least 10 −6 higher than allowed by the unitarity bound. This a safe assumption because, based on the existing estimates of the scaling dimensions of σ and ψ, we expect the smallest twists to be on the order of 10 −2 or more above the unitarity bound [52,53].
Let us briefly comment on another interpretation of these gaps and why some of them may be helpful for isolating the GNY model. One can consider nonlocal deformations of the GNY fixed point, obtained by coupling its low-dimension operators to generalized free fields. Similar nonlocal deformations of the GNY models were discussed in [54] and are generalizations of the description of the long-range Ising model developed in [55][56][57], e.g. one can couple ψ to a generalized free field χ GFF of dimension ∼ 3 − ∆ ψ ≈ 2 and potentially flow to a nearby nonlocal fixed point. Our gap ∆ χ > 3.5 then excludes any such nearby solution to the bootstrap equations. Similarly, one could couple σ to a generalized free field σ GFF of dimension ∼ 3 − ∆ σ ≈ 2.2 − 2.4. Our gaps ∆ σ ′ > 2.5 or 3 similarly exclude these possible solutions. One could also consider nonlocal deformations involving the ϵ, σ T , or / ∂ψ operators, which would also be excluded by our gap assumptions. It would be interesting to give a more systematic study of the nonlocal fixed points (and their dualities) that could be reached by deforming the GNY models.

A distinction between two GNY fixed points
Before diving into the detailed analysis of the condensed matter applications of the O(N ) GNY model in (1.1), we would like to compare this model to another theory, referred to in the literature as the "chiral" GNY model, which possesses an O(N/2) 2 ⋊ Z 2 symmetry. The distinction between the two models and physical examples of each universality class is not always clearly described in the literature. Below we shall show that while the two models flow to different fixed points, many scaling dimensions and OPE coefficients at the two fixed points agree up to a high order in a perturbative expansion. To make the distinction clear, we will henceforth refer to the two theories by their global symmetry groups. 6 In the O(N/2) 2 ⋊ Z 2 GNY model there are two species of two-component Majorana fermions, ψ L i and ψ R i , such that each species has N/2 flavor components. Both species have a Yukawa coupling to a pseudoscalar ϕ, but with opposite signs: (not what we study) . (2.1) Here i = 1 . . . N 2 and A = L, R. Each species of fermion has its own O(N/2) symmetry. Additionally, there is a discrete Z 2 "chiral" symmetry of ψ L i ↔ ψ R i , ϕ → −ϕ which exchanges the fermion species. Note that this symmetry is not really chiral in (2+1)d since there is no notion of "left" or "right" fermions. 7 In total, the flavor symmetry is O(N/2) 2 ⋊ Z 2 . When the fermions spontaneously generate mass due to ϕ getting a VEV, they preserve a parity and a time-reversal symmetry but break the Z 2 symmetry. The Z 2 symmetry breaking of this theory is characteristic of the so-called chiral Ising universality classes, and it has been studied extensively [19, 25, 26, 29-34, 37, 40, 41, 45, 60-72].
At the critical value of m 2 perturbative calculations show that the model (2.1) should also be described by a CFT whose scaling dimensions precisely agree at low perturbative orders with that of the critical model (1.1). Due to this seeming coincidence of the perturbative estimates of the CFT data, the distinction between the O(N ) GNY universality class and the 6 In discussing the related gauged QED3-GN(Y) theories, there is some literature which describes a similar distinction with different nomenclature. The [58,59]. As far as we are aware, the analogous notation has not been used regularly in the gauge-free GN(Y) theories. 7 It is worth noting that the chiral symmetry is indeed related to the spacetime chiral symmetry of a (3+1)d fermionic theory with 4-component fermions. This connection is laid out explicitly in appendix C.
O(N/2) 2 ⋊ Z 2 GNY a.k.a. chiral Ising universality class has been unresolved. 8 We will now clarify this ambiguity and show that the two models are distinct when computing observables to higher perturbative orders. Let us thus compare the two models in the large-N expansion. To compute correlators in the two models, we consider the Feynman diagrams with the leading order propagators of ϕ and ψ i (or ψ L,R i ) denoted by In both: (2. 2) The only relevant interactions in the two models are of the form ϕψψ and ϕ(ψ L ψ L − ψ R ψ R ), respectively, and corrections to the two-point function of ϕ k are given by fermion loops whose vertices involve such interactions. Since in a fermionic loop, both the fermions ψ L and ψ R can always propagate, the only distinction between the large-N Feynman diagrams of the two models can come from fermionic loops with an odd number of vertices. In the O(N/2) 2 ⋊ Z 2 GNY fixed point, the contribution of such loops is always vanishing since for each loop in which L fermions propagate, there is a loop in which R fermions propagate that has the opposite sign, This cancellation is due to the presence of the chiral Z 2 symmetry under which the field ϕ is charged. In the O(N ) model, however, such loops are not guaranteed to vanish when the odd number k of fermions propagating through the loop is k ≥ 5. 9 This leads to a five-point point function for the pseudoscalar ϕ In fact it has been stated that the two models can be related through a field redefinition of ψ L and ψ R [17]. However, while that redefinition makes the Yukawa interaction terms in (2.1) to be the same as in (1.1), the fermionic kinetic terms then differ. 9 For k = 3 such loops can explicitly be shown to vanish due to Tr(/ x 12 / x 23 / x 31 ) = 0. Alternatively one can use the fact that the three-point function of pseudo-scalars in 3d is always vanishing to arrive at the same conclusion.
that can be explicitly checked to be non-zero at order 1/N 3/2 . 10 Such higher k-point functions of pseudoscalars are in principle non-zero in 3d CFTs due to the existence of parity-odd kpoint structures for k ≥ 5 [46]. On the other hand, due to (2.3) the five-point point function of ϕ in the O(N/2) 2 ⋊ Z 2 GNY models vanishes. While this already proves that the fixed points in the two models are distinguishable, we would like to see how these differences are manifest in more commonly discussed observables such as the scaling dimensions of ψ i , ϕ, or ϕ 2 . For this we simply have to find the leading diagram that includes such loops with an odd number of fermion vertices. For instance, we find where we show an example diagram (rather than all the diagrams) contributing to the leading order at which the distinction between the two theories is present. Consequently, even in these low-lying scaling dimensions the two models are different at a fairly high order in 1/N . We can similarly determine the large-N expansion for all the other various operators in the O(N/2) 2 ⋊Z 2 GNY theory. Since these calculations require a lengthy discussion of the irreducible representation of O(N/2) 2 ⋊ Z 2 -which is not the symmetry for the theory we rigorously constrain in this paper using the bootstrap method -we discuss these calculations in appendix C.
Similar logic shows that the models differ only at high order in the ϵ-expansion. For instance, it was noted in [30] for the O(N/2) 2 ⋊ Z 2 GNY model at 4-loop order that the ϵ-expansion, continued to N = 1, is inconsistent with the expected emergence of supersymmetry; the authors of [30] found that manually adding a 5-fermion loop diagram contribution restored the superscaling relation. As in the large-N expansion, in the O(N/2) 2 ⋊ Z 2 GNY model, loops with an odd number of propagating fermions have vanishing contributions owing to the sign difference in the Yukawa coupling between fermion species. Conversely, these odd-fermion loops should be generically nonzero for the O(N ) theory starting at loops with 5-fermions (such as (2.4)). Such a prescription of adding back diagrams with an odd-number of propagating fermions can in principle be used when computing the scaling dimension of any operator in the theory and was named DREG 3 by the authors of [30]. Since the scaling dimensions are only affected at high loop order, the differences between the estimates of {∆ ψ , ∆ σ , ∆ ϵ } for the two models are very small (≲ 3 × 10 −6 ) for all values of N that we consider. 11 Nevertheless, when comparing our results to those from the ϵ-expansion we will use the estimates from [31], which rely on the DREG 3 prescription.
Materially for this paper, we see it fit to compare this work's results with previous results for the chiral Ising universality class, since the scaling dimensions are expected to be close to those of the O(N ) GNY models even at low values of N . Nevertheless, our bounds will not be able to rigorously constrain the chiral Ising GNY models. Rigorous bounds for these models, in which the O(N/2) 2 ⋊ Z 2 global symmetry would be explicitly implemented in the bootstrap equations, represent a separate target for the bootstrap to be studied in the future.

Condensed matter applications
As mentioned in the introduction, the universality class of the GNY model is used to describe a variety of quantum phase transitions in condensed matter systems. In this subsection we will list some of the proposals in the literature. Since there can be confusion regarding the two universality classes associated to the O(N ) GNY critical point and the O(N/2) 2 ⋊ Z 2 GNY critical point, we will revisit this point for each quantum phase transition.
D-wave superconductors [20,21]: In [20], the possible quantum critical points in d-wave superconductors were classified according to their order parameter in an effort to describe anomalous behavior in cuprate superconductors. The two transitions relevant to the discussion in this paper are the transition to d x 2 −y 2 + is pairing which is described by the universality class with symmetry O(4) 2 ⋊ Z 2 (the chiral GNY model for N = 8) and the transition to d x 2 −y 2 + id xy pairing which can be seen to correspond to the universality class with symmetry O(8) (the non-chiral N = 8 GNY model). As shown in [20], these are the only two transitions that have a nodal quasiparticle momentum distribution curve with a width proportional to k B T . The additional requirement that the superconductor exhibits negligible scattering along the (1, 0) and (0, 1) directions uniquely isolates the transition d x 2 −y 2 + id xy which thus underpins the importance of the O(8) GNY model for phase transitions in such superconductors.
Chern insulators, and topological superconductors [23]: There are more systems like the d-wave superconductors which have come into focus in the condensed matter community in recent years, which all possess the common factor of time-reversal symmetry breaking (TRSB). TRSB has been known to the condensed matter community for decades, and systems with a broken time-reversal symmetry are known to have the integer quantum anomalous Hall effect. Prototypical examples of these systems include the Haldane model (a Chern insulator in Cartan symmetry class A) [73]. The O(N )-symmetric critical point represents a phase transition between the time-reversal-preserving class DIII and the timereversal-breaking class D, which in the case N = 1 is discussed in [23]. More specifically, the universality class of the O(1) GNY model is expected to describe a quantum phase transition, with emergent supersymmetry, at the boundary of a topological superconductor where the superconducting Majorana edge modes begin to gap out. Phases associated to this transition are expected to be found for a thin film of superfluid He 3 -B.
In general, it is understood that the cases N > 1 also should have the same phase transition between symmetry classes DIII and D, though we are not aware of any proposed experimental design for empirical observation. 12 Graphene [16][17][18][19]: The distinction between the O(N ) and O(N/2) 2 ⋊ Z 2 theories is also apparent in the case of graphene lattices. In [17] the authors exhaustively enumerated the various order parameters that are expected to exist in graphene lattice theories. It's understood that the chiral Ising universality class, with symmetry O(N/2) 2 ⋊Z 2 , describes the semimetal to charge density wave insulator transition on a honeycomb lattice [16]; specifically for N = 8 it describes a theory of spinful fermions and for N = 4 it describes spinless fermions. In these two cases, time-reversal symmetry is preserved, but the fermions are gapped out, leading to an insulating phase. The corresponding time-reversal-breaking O(N ) case is also expected to exist, but with a different order parameter which does break timereversal symmetry.

SDP formulation of general crossing equations
The system of crossing equations that we study in this paper is fairly complicated: it involves many correlation functions, several of which involve operators with non-trivial spin and/or flavor charges. On top of it, we are working with fermions and need to keep track of lots of minus signs associated with permutations. As a result, rewriting our system of equations in an SDP form suitable for numerical analysis is a rather non-trivial task and is especially prone to human error.
Motivated by this, and also with a view towards future applications, we developed a computer code bootstrap-bounds (see appendix A) which handles most of the bookkeeping associated with passing from physically-transparent crossing equations to the SDP form. In this section we describe the basic algorithm that it uses, in the context of a general conformal bootstrap problem.
Let us consider a general conformal bootstrap problem for a set of external primary operators O 1 , · · · , O n . We allow these operators to have arbitrary spins and flavor symmetry representations, and package their dependence on space-time coordinates and the various polarization indices into an abstract argument p, i.e. we write O i (p i ). In particular we assume that the action of all known symmetries on O i is expressed in terms of the action on p i . For simplicity, we assume that all operators can be chosen to be Hermitian in Lorentzian signature (as is the case in our setup), (3.1) 12 The notion of an interaction spontaneously modifying topological order is a matter of great interest, as it is believed that interactions break the class D free fermion classification of the topological invariant from Z to Z16 [74]. In other words, the topological invariants of the O(N )-symmetric, time-reversal-breaking theories should have topological order defined by ν = N mod 16, since there is an adiabatic way to get from a nontrivial invariant with ν = 16 to a trivial invariant with ν = 0. However, this procedure requires breaking the O(N ) symmetry while deforming the theory from the two topological phases.
Furthermore, we introduce the notation O ∆,ρ for the primary operators that can appear in the OPE of O i but are not among the O i , with ∆ denoting the scaling dimension, and ρ all the other quantum numbers (such as spin, flavor symmetry representation, space parity, etc.). First, we consider the three-point functions. We will need the following two types, Let us focus on the former. It can be written as where Q a,ijk (p 1 , p 2 , p 3 ) is some basis of three-point tensor structures which we are free to choose. Intuitively, this statement is clear, but it turns out we need to formalize this a little in order to have a well-defined algorithm. Formally, choosing a basis of Q a,ijk amounts to the following: • For any choice of the ordered triple i, j, k and the index a, specify a function Q a,ijk in three variables p 1 , p 2 , p 3 , which is invariant under all the available symmetries when transformed with the quantum numbers of O i at p 1 , of O j at p 2 , and of O k at p 3 .
Note that the order of the arguments p 1 , p 2 , p 3 is fixed: for a function, we know what is the first, what is the second, and what is the third argument. This means that for fixed i, j, k the functions Q a,ijk , Q a,jik , Q a,kji , · · · all require separate choices. But of course, the corresponding physical correlation functions are all related to each other in the obvious way (taking into the account fermionic permutation signs), and this induces a relation between the corresponding OPE coefficients λ a ijk , λ a jik , λ a kji , · · · . We then demand, as is always possible to do, that the functions Q a,ijk are chosen so that the OPE coefficients λ a ijk are functions of the unordered triple (i, j, k). In other words, so that In terms of Q a,ijk this means where ± account for fermionic permutation signs. Since these signs are included here, (3.5) is true even when some O i are fermions. This choice of Q a,ijk is provided to the algorithm by the user; in this way, the algorithm does not have to know about the permutation properties of the operators, and can reason in terms of the simple coefficients λ a ijk . Furthermore, we require Q a,ijk to be chosen so that λ a ijk ∈ R, which is again always possible to ensure.
Next we consider the three-point functions where the convention is the same. Concretely, we write Now we only have to worry about the ordering in the pair i, j: we agree to always keep the generic operator O ∆,ρ at p 3 . 13 So in this case we need to provide only two structures for each pair i, j: Q a,ij;∆,ρ and Q a,ji;∆,ρ . We again choose them in a way such that The above convention allows us to interface the general conformal block code blocks 3d [50] from our algorithm. In order to produce a conformal block for the four-point function 14 this code requires the user to specify the three-point structures for Our algorithm can simply look up the functions Q a,ij;∆,ρ and Q b,lk;∆,ρ and pass them to blocks 3d. As a result of our conventions, the expansion of the four-point function becomes λ a ij;∆,ρ λ b kl;∆,ρ G ab,ijkl,∆,ρ (p 1 , · · · p 4 ), (3.11) where G ab,ijkl,∆,ρ is the block returned by blocks 3d, and ′ ∆,ρ = ∆,ρ + n denotes the sum over O ∆,ρ appearing in O i × O j OPE plus the sum over the O n appearing in the same OPE.
The above describes the conventions for the three-point functions, but we also need to specify the crossing equations. Any crossing equation involves only one four-point function, expanded in different channels. For a given four-point function there are at most 3 distinct channels, depending on which operator out of O j , O k , O l we take the OPE of O i with. In any case, an equality of two channels in a four-point function takes the form

12)
13 Recall that we have agreed not to use the label O∆,ρ for any of the Oi. 14 The code blocks 3d can only compute the 3d conformal blocks; the flavor structure is then added by an additional layer of code.
where we used | to separate the groups of operators between which we take the OPE. The ± sign is chosen based on the statistics of the operators. Note that expresses the same equality. We will use the convention in which we order operators in such a way that the two terms in the crossing equation differ by swapping the operators at p 1 and p 3 , like in (3.12). We choose some complete and independent set of crossing equations, written in this convention. We can expand each four-point function in a basis of four-point tensor structures, 14) The conformal block expansions take the form where G denotes the blocks computed by blocks 3d. Again, we have a choice to make for the functions T I,ijkl . We don't constrain it in any particular way. However, since the functions above differ only by the order of operators, we must have for some matrix M J I (which depends on the choice of i, j, k, l). The crossing equation then takes the form The appearance of 1 − z is due to our convention for the crossing equation and the choice of cross-ratios. We can furthermore expand the crossing equations (3.19) in a power series around z = z = 1 2 . In general, not all the Taylor coefficients are going to be linearly-independent [46,75,76]. We therefore choose some independent set, with a cut-off n max on the order. We then introduce a combined label I which runs over all crossing equations, all four-point tensor structures I and the independent Taylor coefficients up to the cutoff n max . We write g I ijkl for the contribution of the four point function OPE channel ⟨O i O j |O k O l ⟩, in this specific ordering of indices, to the scalar crossing equation labeled by I. 15 For example, if I corresponds to I = 0 and the (z − 1 2 ) m (z − 1 2 ) n term in the equation (3.19), then we have with all other g I ··· vanishing. With this notation, the full set of crossing equations takes the form where we sum over all choices of ijkl. Ultimately, the coefficients g I ijkl are provided to the algorithm by the user. In our implementation the user provides the equations (3.19) and the list of Taylor coefficients that need to be considered, from which the code reads off the g I ijkl . We can extend this notation in the obvious way to G, so that We now consider three types of contributions to (3.24) Generic contributions We start with the contributions of generic operators O ∆,ρ . After imposing a cut-off on the spin of O ∆,ρ there are finitely many distinct ρ appearing in (3.24). We call such ρ "operator channels." To each operator channel we associate the set I ρ = {(a, (i, j))|Q a,ij;∆,ρ ̸ = 0}, (3.25) where (i, j) denotes an unordered pair. In other words, the elements of I ρ label the OPE coefficients λ a ij;∆,ρ which are allowed by symmetries. We now introduce the PSD matrix (P ∆,ρ ) αβ with α, β ∈ I ρ (P ∆,ρ ) αβ = degeneracies λ a ij;∆,ρ λ b kl;∆,ρ ⪰ 0, (3.26) 15 Note that for some orderings of ijkl it might be that g I ijkl vanishes for all I due to the choice of the crossing equations. Furthermore, some four-point functions vanish by symmetry and the corresponding g I ijkl are also 0.
With this notation and using (3.5), the contribution of generic operators to the crossing equations is (3.28) External contributions We now focus on the contribution to the crossing equations of O i themselves. We define the index set where (i, j, k) denotes unordered triples. That is, E labels the OPE coefficients λ a ijk which are not forced to be 0 by symmetries. We now define the rank-1 PSD matrix, assuming no degeneracies among the quantum numbers of O i , where α = (a, (i, j, n)) and β = (b, (k, l, m)).
We define the symmetric matrix G I ext by the requirement that the contribution of the external operators to the crossing equation (3.24) is (recall (3.5)) 0 = Tr P ext G I ext + · · · . (3.31) The matrix G I ext has a straightforward expression in terms of G I ab,ijkl,∆,ρ which is however awkward to describe.
Special operators There are sometimes special exchanged operators such as the stresstensor T µν or the identity operator 1. These operators contribute to the crossing equations as where G I is determined by the special block for the given operator. It may or may not depend on parameters such at OPE coefficients or some scaling dimensions. For example, for the identity operator the contribution is completely fixed by the normalization of two-point functions of O i .

Final crossing equation
The final crossing equation then takes the form 34) and potentially additional contributions from T or other special operators. Here, P ext , P ∆,ρ ⪰ 0. So, for example, in a feasibility study we look for functionals F I such that Crucially, these conditions can be formed automatically once the user provides the three-point structures Q and the crossing equations in the form (3.19).

Three-point functions
Returning to the bootstrap problem for the GNY model, each local primary operator is characterized by the scaling dimension and three other quantum numbers, namely, spin j, space parity P , and an O(N ) irreducible representation µ. 16 We consider the mixed system of the lowest dimension operators ψ, ϵ and σ, where their quantum numbers (j, P, µ) are When N = 1, the global symmetry becomes O(1) = Z 2 . Since the non-trivial element sends ψ → −ψ and ϕ → ϕ, it coincides with (−1) F and should not be considered as a separate global symmetry. Hence, in this case the operators are labeled by (j, P ) only. 17 Consequently, there are no flavor structures to consider. There is a corresponding reduction in the number of equations and OPE channels.
The crossing equations under study involve the following set of four point functions {⟨ψψψψ⟩, ⟨ϵϵϵϵ⟩, ⟨σσσσ⟩, ⟨ψψϵϵ⟩, ⟨ψψσσ⟩, ⟨σϵψψ⟩, ⟨σσϵϵ⟩}. The conformal block expansions of these correlation functions involve OPEs between all combinations of ψ, ϵ, σ. 16 Note that the parity makes sense both for integer and half-integer j (we use the transformation defined in [46]), although for the latter the notions of "even" and "odd" parities can be exchanged by redefinition of the parity transformation by (−1) F . This allows us to choose the parity for one fermionic operator at will, so we choose ψ to be parity-even. 17 In our discussion, which is valid for generic N , one can take • → • and → • while all other representations are omitted.
We list the tensor structures appearing in relevant three-point functions in table 3 and  table 4. We choose to represent the conformal structures in the SO(3) r basis following the convention in [47,50], 18 and we define the O(N ) flavor structures as follows, Some conformal structures might disappear for low spin of the exchanged operator O. The selection rule is that |j 12 , j 123 ⟩ is present in ⟨O 1 O 2 O⟩ if j 12 ∈ j 1 ⊗ j 2 and j 123 ∈ j 12 ⊗ l, where j 1 and j 2 are the spins of O 1 and O 2 . In practice, to obtain the conformal structures, we did calculations in the q-basis [46] and then converted them to the SO(3) r basis as discussed in [50]. Tensor structures for three external operators are obtained from those in tables 3 and 4 by restricting O to the relevant special case. The only exception is the structure for ⟨ψψσ⟩ which we take to be −|1, 1⟩ (differs by a factor of −1 from table 4).
Our choice of the tensor structure basis ensures the equality of the following OPE coefficients, as required by the algorithm in section 3.1: and similarly λ ψσO = λ σψO etc. A slight difference from 3.1 is that the OPE coefficients for our structures are purely imaginary if they involve fermionic operators. This difference is accounted for in the software. The stress-tensor T and the conserved O(N ) current J transform with (l, P, µ) equal to (2, even, •) and (1, even, ). Their OPE coefficients are constrained by Ward identities as OPE O ∈ (l, P, µ)  where T = C (3.52) In practice, the Ward identity for J does not affect the numerics because it only gives an interpretation of an OPE coefficient in terms of C J . The Ward identity for T doesn't affect the numerics unless a gap above T is assumed. As explained in section 2.3, we do assume a gap of 10 −6 above T . However, since this gap is small, it is likely that imposing the T Ward identity doesn't have a significant effect on our results.

Four-point functions and crossing equations
We now study the four-point functions of our mixed system. Following the procedures outlined by (3.12), (3.14), and (3.15), we construct the four-point tensor structures T I,ijkl and find the crossing equations in the form of (3.19).
In general, one single four-point structure is the product of a flavor and a conformal structure. The full structure should be invariant under kinematic permutations, which is the group that preserves all conformal cross-ratios formed by the coordinates of the four operators. We refer to [46] for the details on this group.
When the four operators are identical, ρ 1 = ρ 2 = ρ 3 = ρ 4 and the kinematic permutation group is Z 2 ×Z 2 = {e, (12) To ensure that the full structure is invariant under the kinematic permutations, the flavor and conformal structures should transform in the same irrep of the kinematic permutation group.
There are two additional requirements for the four-point tensor structures. Firstly, the structures must have space parity consistent with the parity of the operators in the four-point function. And secondly, they must have definite parity under the transformation z ↔ z, which we refer to as the t-parity transformation. 19 This is needed to simplify the corresponding symmetry of the coefficient functions g(z, z), see [46].

⟨ψψψψ⟩
We first consider the four point function of four fermions ⟨ψ i ψ j ψ k ψ l ⟩ = I,a t I T ijkl a g I,a ψψψψ (z, z), (3.53) where T a stands for the flavor structures and t I for the conformal structures. It participates in the crossing equations of the form understood in the sense of equation (3.12). There are three possible flavor structures for ⟨ψ i ψ j ψ k ψ l ⟩, namely δ ij δ kl , δ ik δ jl , and δ il δ jk . They all are invariant under the kinematic permutation group Z 2 × Z 2 , i.e. are all in the (++) representation. For future convenience, we define the following linear combinations of these structures, (3.55) 19 Not to be confused with time reversal.

(3.69)
When Taylor-expanding these equations, one should keep in mind the t-parity of the conformal structures which determines the parity of the above functions under z ↔ z. Finally, the crossing equations for g ⟨↑↑↑↑⟩ + ,Ta (z, z) at z = z are redundant with other crossing equations and should not be imposed in order to avoid numerical instabilities (see appendix A of [46]). In practice this is done by requiring that n > 0 derivatives are taken in the direction orthogonal to z = z for this structure.

⟨ψψσσ⟩ and ⟨ψψϵϵ⟩
Mixed four-point functions containing both σ and ψ give rise to two crossing channels in the sense of (3.12), We are slightly abusing the notation since the conformal structures are different for each of the three orderings (since different operators are inserted at different points x i ). Only one flavor structure T ij a = δ ij exists for each of these four-point functions, and it is kinematic-symmetric. It is also mapped to itself under the crossing permutation (13). Hence, we will ignore the flavor structure in the following discussion. Furthermore, since the products σσ and ϵϵ have the same parity, the analysis below is the same for the correlation functions ⟨ψψσσ⟩ and ⟨ψψϵϵ⟩ . We will focus on the correlator ⟨ψψσσ⟩ for concreteness.
The two parity-even conformal structures for the ordering ⟨σσψψ⟩ are [0, 0, 1 2 , 1 2 ] and [0, 0, -1 2 , -1 2 ], and similar expressions apply for the other two orderings. These structures are automatically kinematically symmetric. We will simplify the notation and write them as [↑↑] and [↓↓] for each of the orderings (i.e. omitting the 0 charges), and define the structures with definite t-parity as (3.73) Applying (3.60) and factors from fermion permutations, we find that conformal structures with the same label are mapped into each other under the (13) permutation. Hence, we can form crossing-symmetric and anti-symmetric combinations as The convention for the crossing equations is the same as in the previous subsection. This time, there are no redundancies between the crossing equations.

⟨σϵψψ⟩
The independent crossing equations in this case are Similar to the case above, we have a trivial flavor structure that can be ignored. However, the overall parity ⟨σϵψψ⟩ is now odd and so a separate analysis of the conformal structures is required. For all of the orderings of the operators, the parity-odd conformal structures are [↑↓] and [↓↑], using the same notation as in the previous subsection. There is no kinematic symmetry to consider. To form structures with definite t-parity, we write (3.78) Applying (3.60) and taking into account the factors from fermion permutation, we can form the crossing-symmetric and crossing-antisymmetric linear combinations as before,

Scalar four-point functions
Since both flavor and conformal structures are trivial, it is a straightforward exercise to form crossing-symmetric and crossing-antisymmetric functions, symmetric: g σσσσ , g ϵϵϵϵ , g σϵσϵ , g σσϵϵ + g ϵσσϵ , (3.81) The t-parity of all these functions is +1. In total, we end up with 38 crossing equations for N ≥ 2 and 28 equations for N = 1.

Numerical computations
After setting up the crossing equations and OPE channels, our workflow is automated by softwares described in appendix A to search for the functionals F I to satisfy (3.35-3.37). The space of CFT data that we searched over is 6-dimensional, but is better understood as two separate searches done together. We have three scaling dimensions and three OPE coefficient ratios that we search over; the latter is done by a cutting-surface search algorithm, and the former is done by a Delaunay mesh search algorithm. We will now briefly go over these search algorithms.

OPE scan
As noted in [8], by including assumptions of OPE coefficient ratios in our bootstrap computations, the allowed region in the space of scaling dimensions can be improved at the cost of also having to search over those OPE coefficient ratios. To that end, we employ the algorithm described in [12] to include constraints imposed by the OPE coefficients involving the only the external operators, λ ext . The non-vanishing OPE coefficients are (3.83) With (3.36), we require that the contribution of external scalar OPE coefficients to the crossing equation has a definite sign after applying the functional, independent of the values of those coefficients.
However, (3.36) is strong enough that the conclusion above stands for any matrix with the decomposition M ext = A † A. We, on the other hand, are only interested in the case M ext = P ext = ⃗ λ ext ⃗ λ T ext , which is a rank-1 matrix. Therefore, we instead look for functionals F I 's that satisfy for each [λ ext ] ∈ RP 3 (independent of the magnitude or sign of the vector), along with (3.35) and (3.37). A point in the dimensional space (∆ ψ , ∆ σ , ∆ ϵ ) is ruled out if such F I 's exist for all [λ ext ], and is allowed otherwise. Hence, by scanning over the OPE space, we obtain a union of allowed regions in dimension space, each permitted by some OPE direction: Note that this union of allowed regions is contained in the allowed region obtained by imposing the stronger (3.36), as this is the special case when imposing the condition that M ext is rank-1. We should note that for more than one OPE coefficient ratio, the cutting surface algorithm is non-rigorous [12]. For each point in the dimension that the cutting surface algorithm doesn't exclude, it outputs a direction in the OPE space that could not be excluded by our constraints. Combined over all the allowed points in the dimension space, this produces a list of OPE coefficient ratios with a relatively small variance. We report the full range of values of the OPE coefficient ratios that we found in our searches in section 4 as estimates of the real ratios. It should be noted however that these estimates are not rigorous. In particular, we cannot exclude a systematic error that could arise from the way the cutting surface algorithm samples the OPE ratios. Similarly, the error bars are not rigorous.

Delaunay mesh search
The islands that we show were computed by the Delaunay mesh search algorithm; we will refer for details on the algorithm to the original paper [12]. However, in order to properly interpret how we have chosen to represent our results, we will provide a brief qualitative review of this method.
The principle of the algorithm is to divide the search space into a suitable simplicial complex, known as a Delaunay mesh, where each vertex is a point that has already been computed as being allowed or disallowed. A simplex that has entirely allowed or disallowed vertices is assumed to be completely on the interior or exterior of our island. If a simplex has both an allowed and disallowed vertex, then the simplex is deemed to be on the boundary of the island; we can think of the volume of the simplex as the region of uncertainty between allowed and disallowed. This implied boundary can be further refined by computing the feasibility of the point at the simplex's centroid. 22 Thus, with each iteration we get an increasingly refined mesh.
After a search has finished, we can determine the island's boundary in one of two ways. The first is to take the centroids of the boundary simplicies and compute their convex hull. This method produces fairly smooth islands, but whose bounds are not strictly rigorous. The second, more conservative, approach is to take the convex hull of all boundary simplices, which is equivalent to taking the convex hull of all vertices that neighbor allowed vertices. This includes the interiors of all boundary simplices into our island so we can be sure that the area outside this hull is strictly disallowed by our bootstrap constraints, up to assumptions of convexity. Thus, the latter method produces islands and bounds that are rigorous. However, it should be noted that the latter method's islands tend to be more jagged than islands computed with the centroid method given the same set of points. This distinction is heightened when 22 The Delaunay search algorithm selects the centroid, which is also the mean of the vertices, as the next point to be computed in the numerical bootstrap [12]. Interpreting the simplex as the uncertain region of the boundary, the centroid is the mean of that uncertain region.  Figure 5: The allowed region for the N = 2 critical GNY model when imposing that ∆ σ ′ > 3 computed at n max = 18, projected to the (∆ σ , ∆ ϵ ) and (∆ ψ , ∆ σ ) planes. This should be compared to the Borel re-summed result obtained from the ϵ-expansion [31] (shown in orange) and with the Monte Carlo results from [39] (shown in purple).
working with relatively small numbers of points computed. In this work, we have opted to take the more conservative approach.

Results
In this section, we present the results of our numerical bootstrap computations for various values of N . While the space of CFT data that we searched over is 6-dimensional (as discussed in section 3), in most cases, we have projected the results into the planes (∆ ψ , ∆ σ ) and (∆ σ , ∆ ϵ ) for ease of visualization. We first discuss the bounds obtained for N = 2, 4, and 8 and compare them with the existing studies from ϵ-expansions (after Borel resummation) and from Monte Carlo simulations [31,[39][40][41]. We then focus on the N = 1 theory and discuss how our results, without assuming supersymmetry a priori, are strong evidence for the emergence of supersymmetry in the IR in the N = 1 critical GNY model. For N = 2 at n max = 18 we can report rigorous estimates of our external scaling dimensions of ∆ ψ = 1.06861 (12), ∆ σ = 0.6500 (12), ∆ ϵ = 1.725(7), given the assumptions as discussed in section 2.3 of ∆ ϵ ′ > 3, ∆ σ ′ > 3, ∆ σ T > 2, ∆ ψ ′ > 2, ∆ χ > 3.5. Figure 5 shows the allowed  Figure 6: The allowed region for the N = 4 critical GNY model when imposing that ∆ σ ′ > 3 computed at n max = 18, projected to the (∆ ψ , ∆ σ ) and (∆ σ , ∆ ϵ ) planes. This should be compared to the Borel re-summed result obtained from the ϵ-expansion [31] obtained using the DREG 3 regularization scheme (shown in orange) and with the Monte Carlo results from [40] (shown in purple). Results for ∆ ψ are not available from the Monte Carlo study cited. regions 23 after projections to the (∆ σ , ∆ ϵ ) and (∆ ψ , ∆ σ ) planes. As shown in the (∆ σ , ∆ ϵ )plane, the bootstrap results exclude the reported error bars from earlier Monte Carlo studies from [39]; as can be seen in both plots, the bootstrap results also marginally exclude the reported error bars from the ϵ-expansion results after Borel resummation [31]. For both studies, the bootstrap results improve the precision of some of these estimates by orders of magnitude.
We also report (nonrigorous) estimates at N = 2 at n max = 14 24 of the OPE coefficient ratios of λ ψψσ λ σσψ = 0.5087(10), λ ψψϵ λ σσψ = 0.2392 (6), and λϵϵϵ λ σσψ = 1.629 (13). Note that these estimates are nonrigorous as discussed in section 3.4.1; for the reader's convenience they are also reported in table 5. 23 Readers familiar with bootstrap results may be concerned as to the jagged nature of the islands reported.
The nmax = 18 computations are very computationally expensive, so we only have a relatively small number of allowed points in each island. However, because we have a lot of information as to the disallowed points from lower nmax, and because the precision of these islands are still improved, we have elected to report the superficially more jagged islands as they represent our best results. We have taken pains to ensure that these results are rigorous, as outlined in section 3.4.2. 24 We report our nmax = 14 estimates because we ran out of computing resources while computing nmax = 18.
While we have enough data to feel confident in our rigorous scaling dimension estimates, we deferred to our lower nmax results for the OPE coefficient ratio estimates.  Figure 7: The allowed region for the N = 8 critical GNY model when imposing that ∆ σ ′ > 3 computed at n max = 18, projected to the (∆ σ , ∆ ϵ ) and (∆ ψ , ∆ σ ) planes. This should be compared to the Borel re-summed result obtained from the ϵ-expansion [31], once again obtained using the DREG 3 regularization scheme (and shown in orange), and with the Monte Carlo results from [41] (shown in purple).
As noted in section 2.3, the assumption of ∆ σ ′ > 3 for N = 2 is perhaps at risk of being too strong because of the N = 1 value which violates this assumption. We note that despite this assumption we are still able to find feasible points. Moreover, the two-sided Padé interpolations shown in figure 3 support this assumption, so we think this gap is justified. In future work we plan to study the CFT data in the σ ′ sector more rigorously using the navigator method [77]. For now, we report our estimates using both ∆ σ ′ > 3 as well as the more conservative assumption of ∆ σ ′ > 2.5 in tables 1 and 5.

N = 1 and emergent supersymmetry
We also computed islands at n max = 6, 10 for the N = 1 case, shown in figure 8. Shown also in the plot is the expected relation between scaling dimensions for an N = 1 SCFT; our three external operators are expected to all be in a supermultiplet with each other. We can see that the very tip of the island indeed is consistent with the assumption of supersymmetry. There is a long tail, however, which seems to get cut away as n max is increased. The intersection of the tip of the island with the supersymmetric line is very narrow, and at n max = 10 we can report a rigorous estimate of ∆ σ of 0.58444 (8), which is both completely consistent though roughly an order of magnitude less precise than the n max = 30 superconformal bootstrap results reported in [44]. 26 One advantage of our mixed fermion bootstrap setup is that we now have access to half-integer spin exchange channels. In this case, the parity-odd ℓ = 3/2 channel would be expected to have a supercurrent at the unitarity bound of ∆ = 2.5 if the solution to crossing corresponds to a supersymmetric CFT. We can therefore investigate whether a given solution to crossing in our N = 1 island must have supersymmetry by imposing gaps in the supercurrent channel and seeing what the upper bound of that gap is. If the feasibility of a solution to crossing is sensitive to this gap assumption, we can say that in that solution there must be a low-lying operator in that spectrum.
Specifically, we can perform a binary search in the supercurrent gap to find precisely what is the upper bound of the scaling dimension of the leading operator in that channel. Scanning along an axis of our n max = 10 island, which goes from the tip through the center (shown in figure 9), we find that for the entire scan ∆ SC < 2.54. In particular, at the tip of the island (corresponding to the N = 1 super-Ising model), the upper bound drops to ∆ SC < 2.5003219. This implies that any CFT with these parameters must be, to a high degree of precision, supersymmetric. 27

Discussion
In this work we have obtained the first rigorous and precise islands for the conformal data of the 3d O(N ) Gross-Neveu-Yukawa fixed points from the conformal bootstrap. Much like the 3d Ising and O(N ) vector models, these theories appear to be readily amenable to bootstrap methods. In particular, we have shown that one can obtain small islands in the parameter space of CFT data after studying the crossing relations for the operators {ψ, σ, ϵ} and imposing gaps in the spectrum which isolate the leading scalar and spin-1/2 operators. The gaps we have chosen are motivated by perturbative calculations in the large-N and ϵ-expansions. They can also be viewed as being necessary in order to exclude potential nearby nonlocal fixed points obtained by coupling the leading operators in the GNY theories to generalized free fields. In the χ and σ ′ sectors the gaps can be viewed as imposing a consequence of the equations of motion for the fundamental fields. We saw a particular sensitivity of our results to the gap in ∆ σ ′ , and establishing its irrelevance down to N = 2 appears to be an important open problem. In future work we plan to study the allowed values of this gap more rigorously using navigator methods [77] and use the results to further improve our islands.
In the present study, we focused on the O(N )-invariant GNY theories with N = 1, 2, 4, 8, but it is clear that the results can be extended to any N . In the case of N = 1, the fixed 26 While the nmax values reported here are quite different, we should note that the system of crossing equations studied in [44] has only 4 crossing equations while our N = 1 system has 28. 27 We also did a preliminary exploration of how the bound on lowest spin-3/2 operator changes as a function of the gap in ∆ σ ′ . As this gap is increased towards the value that it takes in the N = 1 super-critical Ising model, we saw that the upper-bound on the dimension of the spin-3/2 operator becomes even stronger, as the right part of the plot (the tail of the island) shrinks away.  Figure 8: The allowed region for the N = 1 critical GNY model when imposing that ∆ σ ′ > 2.5, projected to the (∆ ψ , ∆ σ ) and (∆ σ , ∆ ϵ ) planes. The solid lines capture the expected relation between scaling dimensions in an N = 1 SCFT. The tip intersects with the supersymmetric constraint on scaling dimensions, though note that due to the projection and visualization method (described in section 3.4.2) that the extent of that intersection is exaggerated, as the line and island are actually in 3-dimensional space. We separately computed a binary search along the supersymmetric line in all three external dimensions and found a rigorous estimate of ∆ σ ′ = 0.58444 (8). This agrees exactly with those found for the N = 1 super-Ising model in [44].
point is believed to have emergent 3d N = 1 supersymmetry. If one assumes supersymmetry then one can perform a precise bootstrap of this model using only external scalar fields as was demonstrated in [44]. If one does not assume supersymmetry, then we have seen in the present work that the numerical bootstrap with external fermions still forces the solution to be approximately supersymmetric, requiring a spin-3/2 "supercurrent" operator in the spectrum that is near the unitarity bound.
In this paper we also highlighted the distinction between the O(N ) GNY fixed points and the O(N/2) 2 ⋊ Z 2 GNY fixed points. For the leading operators {ψ, σ, ϵ}, differences in their scaling dimensions only show up at 4-loop order in the ϵ and large-N expansions and they are expected to be extremely close to each other. On the other hand, the models have more significant differences in other parts of the spectrum, e.g. in the number of conserved currents and spectrum of fermion bilinear operators. In the future it will be interesting to study the O(N/2) 2 ⋊ Z 2 GNY fixed points using numerical bootstrap methods and see if we can clearly resolve the difference between these models. It will additionally be interesting to  Figure 9: On the right we show the upper-bound at n max = 10 on the scaling dimension of the lowest spin-3/2 operator as a function of ∆ ϵ obtained along an axis of the N = 1 island when imposing ∆ σ ′ > 2.5. The axis was selected such that it interpolates between the location of the N = 1 super-Ising model as determined in previous literature [42][43][44] and the center of the n max = 10 island. We show this axis projected into (∆ σ , ∆ ϵ ) space in orange in the figure on the left. We see that all points along this axis are forced to have a spin-3/2 operator that is close to the unitarity bound, where such operator becomes a supercurrent. The red circle in the plot on the right gives the location of the critical N = 1 super-ising model as determined from the N = 1 superconformal bootstrap [44]. study bounds on OPE coefficients and central charges in both models. These fixed points can also readily be generalized to models with multiple scalar fields, e.g. the chiral-XY and chiral-Heisenberg GNY fixed points, which will also be interesting targets for the bootstrap. Our work was the first significant application of the new software tool blocks 3d [50], which enables systematic and efficient calculations of 3d conformal blocks with arbitary spinning operators. Along with it, we have developed an extensible and modular software stack, as described in appendix A. The success of our study makes it clear that this approach can be used in other studies of interest, e.g. mixed correlators containing various combinations of scalars, fermions, currents, and stress tensors. We hope that future systematic studies of the bootstrap constraints for such correlators will lead to the discovery of exciting new islands in the vast ocean of possible CFTs.
B ∆ χ and ∆ ψ ′ at large N To better isolate the O(N ) GNY-model using the conformal bootstrap it is useful to get a better estimate for the scaling dimension of various operators in the theory, especially for those that are involved in the equation of motion or in the new fermionic channels that are involved in a fermionic-scalar OPE. In particular, we would like to compute the 1/N correction to the dimension of the fermionic operators ϕ 2 ψ i (the lowest dimensional operator above ψ i ) and ϕ 3 ψ i (the lowest dimensional primary in its parity sector and O(N ) representation). We will determine these scaling dimension by extracting the logarithmic divergence coefficient in the propagator of the more generic operator O i = ϕ k ψ i . The bare dimension of this operator is ∆ ϕ k ψ i = k + 1 + O(1/N ). We will follow a similar logic to that used to determine the dimension of ϕ k in appendix B of [15].
Up to order 1/N , this propagator (in the concrete, but generalizable, example with k = 2) is given by the following diagrams: The first diagram gives the leading bare propagator. The next two subleading terms (i.e. the one proportional to k and the one proportional to k(k−1)/2) capture the anomalous dimension of ϕ k , which was found to be δ ϕ k = 1 N 16k(3k−5) 3π 2 [15]. The next term yields the 1/N correction to the fermionic propagator and gives the contribution of the anomalous dimension of the fermionic field, δ ψ i = 1 N 4 3π 2 . Finally, the last term is a new diagram which yields a correction kη to the anomalous dimension. Thus, the scaling dimension of the operator is given by, To determine η we note that in the special case when k = 1, following from the equations of motion, the operator O i is a descendant and therefore, ∆ ϕψ i = 1 + ∆ ψ i . This implies that the logarithmic divergences for the following diagrams should cancel each other: from which we conclude that η = −δ ϕ = 32 3π 2 N . Consequently, we find: In particular, we find C More on the differences between the O(N ) and O(N/2) 2 ⋊Z 2 GNY models

C.1 2-vs 4-component
We will note that in the literature sometimes authors will describe the theories discussed in this paper using the language of 4-component fermions as opposed to the natively 3dimensional 2-component language which we use in this paper. This is especially true in the case of 4 − ϵ calculations. A Lagrangian may look like where Ψ i are 4-component Dirac spinors 29 and i = 1 . . . N 4 , N 4 ≡ N/4. Notably, the Lagrangian has an additional discrete "chiral" symmetry of Ψ i → γ 5 Ψ i , ϕ → −ϕ, inhereted from the four-dimensional theory. 30 Models of this form have been studied extensively to understand the chiral Ising universality class, which for N = 4, 8 (i.e. N 4 = 1, 2) describe spinless/spinful critical points for the semimetal-to-CDW transition in graphene. The critical points describe breaking of the symmetry while preserving time-reversal symmetry [30,31].
To make the distinction between this theory more clear from the theory that is the focus of our paper, it is helpful to consider decomposing our 4-component Dirac spinors into 2component Majorana spinors. In a suitably convenient 4 × 4 basis of the gamma matrices, such as the following defined in [59,78] in terms of the 2 × 2 γ µ we have used throughout this paper (defined in [15]), we have: Here we define Ψi ≡ Ψ † i γ 0 , where γ 0 is defined in C.2. The Clifford algebra satisfies { γ µ , γ ν } = 2η µν , just like the 2 × 2 gamma matrices defined in [15] that we use in this paper. 30 We denote the "chiral" symmetry in quotes as there is no inherent spacetime notion of chirality in three (or indeed, any odd number of) dimensions. This symmetry amounts to an internal flavor symmetry of the fermions.
We should note that the γ 3 will not play a role in the 3d theory; thus this gamma matrix basis reduces to a block-diagonal form in 3d. Therefore, the 4-component spinors can be broken down into 2-component spinors as We get (This can be further decomposed into the requisite Majorana spinors without significant difference in form, save for the index i going from 1 to 2N 4 .) In this form, we see that we have two distinct fermion species with opposite signs on their Yukawa couplings; the chiral symmetry becomes We can see that, if there are N Majorana fermions total, the symmetry of this Lagrangian is O(N/2) 2 ⋊ Z 2 , where Z 2 is the chiral symmetry. For completeness, we will note that the 4-component spinor notation for the fermion bilinear which appears in (1.1) is iΨ i γ 3 γ 5 Ψ i .

C.2 Irreps of O(N/2) 2 ⋊ Z 2
In the next subsection we will classify the various low-lying primary operators of the O(N/2) 2 ⋊ Z 2 GNY theory. To attempt such a classification, we should first discuss the irreps of O(N/2) 2 ⋊ Z 2 . We will start with a more general discussion. Given a compact simple Lie group G, let H = (G × G) ⋊ Z 2 . The group H is generated by (g L , g R ) ∈ G × G, together with an element s such that s 2 = 1 and s(g L , g R ) = (g R , g L )s. (C.5) For each irrep ρ of H, we can consider its restriction to G × G. The irreps of G × G have the form ρ 1 ⊠ ρ 2 where ρ 1 , ρ 2 are irreps of G. The symbol ⊠ means we take a tensor product as vector spaces, but not as G representations. The first G acts on the left-tensor factor and the second G acts on the right tensor factor.
• Suppose that ρ 1 = ρ 2 . In this case, we claim that ρ ± = (ϕ ± ϕ ′ )(ρ 1 ⊠ ρ 1 ) are Hinvariant subspaces of ρ. Clearly they are G × G-invariant subspaces. To prove they are also H-invariant, note that By irreducibility, one of the ρ ± must be all of ρ, and the other must vanish. When ρ = ρ + or ρ = ρ − , we denote the H-representation by ⟨ρ 1 , ρ 1 ⟩ ± . A basis of ⟨ρ 1 , ρ 1 ⟩ ± is given by u ⊗ v with u, v ∈ ρ 1 , with the s-action s(u ⊗ v) = ±v ⊗ u. (C.12) Below, when listing the perturbative estimates for the scaling dimensions for some of the low-lying operators in the theory we will be interested in ρ i ∈ {•, , , } which corresponds to the singlet, vector, symmetric traceless tensor and antisymmetric representations.
where we have again listed examples of diagrams (rather than all the diagrams) contributing to the leading non-vanishing order. As explained in section 2, these diagrams cancel in the O(N/2) 2 ⋊ Z 2 GNY models. More generally, since such loops only contribute at higher order we see that at leading order where O can represent an operator with any spin and parity representation, ρ can be one of the { , , . If that is not the case, a more elaborate analysis is needed. For example, consider the operator σ ⟨•,•⟩ + ∼ ψ L i ψ L i + ψ R i ψ R i . 31 The diagrams contributing to the two-point function of this operator are precisely the same as those contributing to the two-point function of σ ⟨ ,•⟩ . The fact that the O(N/2) indices are contracted differently does not matter at low enough order when computing the scaling dimension of the operator and only affects the overall normalization of the two-point function. Therefore, we find that (C. 15) This approximate relation between the scaling dimensions in the singlet sector and those in the ⟨ , •⟩ irrep can be extended to the parity even sector. For instance, consider the 1 + ⟨ , ⟩ − 2 + 16 3π 2 N + . . . − Table 6: (C.16) Next, we discuss some of the operators in the O(N/2) 2 ⋊ Z 2 model whose scaling dimensions are not found among operators in the O(N ) model even at leading order in the large-N expansion. One example of such an operator is σ ⟨ , ⟩ + ∼ ψ L i ψ R j + ψ R i ψ L j . The two-point function of such an operator up to order O(1/N 2 ) is given by To extract the anomalous dimension of the operator, we need to read-off the coefficient of the logarithmic divergence in (C.17) for the last three diagrams δ σ ⟨ , ⟩ + = 2δ ψ − η vertex which are responsible for the 1/N correction. The first two diagrams contribute the same amount since δ ψ = δ ψ L = δ ψ R , while the last diagram was computed in appendix B of [15] in order to compute the anomalous dimension δ ϵ in the O(N ) model. The calculation is almost identical to there, the only difference is that the last diagram in (C.17) has the opposite sign due to the difference in sign between the Yukawa coupling of ψ L and ψ R . Therefore, Similarly, we see that the diagrams needed to compute the scaling dimension of ϵ ⟨ , ⟩ − ∼ ϕ(ψ L i ψ R j + ψ R i ψ L j ) and ϵ ⟨ , in the antisymmetric representation of each O(N/2) subgroup. However, to further distinguish the two models and to ensure that the O(N/2) 2 ⋊ Z 2 does not have a symmetry enhancement in the IR to a Lie group with a greater number of generators, we can compute the anomalous dimension of the bifundamental current j µ ⟨ , ⟩ − ∼ ψ L i γ µ ψ R j − ψ R i γ µ ψ L j . The computation for the anomalous dimension of this spin-1 operator follows from the computation which shows that the conserved current j µ