Bootstrapping 3D fermions with global symmetries

We study the conformal bootstrap for 4-point functions of fermions 〈ψiψjψkψℓ〉 in parity-preserving 3d CFTs, where ψi transforms as a vector under an O(N ) global symmetry. We compute bounds on scaling dimensions and central charges, finding features in our bounds that appear to coincide with the O(N ) symmetric Gross-Neveu-Yukawa fixed points. Our computations are in perfect agreement with the 1/N expansion at large N and allow us to make nontrivial predictions at small N . For values of N for which the Gross-Neveu-Yukawa universality classes are relevant to condensed-matter systems, we compare our results to previous analytic and numerical results.

In [13] we initiated a study of the conformal bootstrap applied to the 4-point function ψψψψ of a single Majorana fermion ψ of scaling dimension ∆ ψ that belongs to a 3D CFT with parity symmetry. We found that by varying the gap between the scaling dimensions ∆ σ and ∆ σ of the first (σ) and the second (σ ) parity-odd scalar appearing in the ψ × ψ OPE, the resulting allowed region in the {∆ ψ , ∆ σ } plane showed features (kinks) that tracked the GNY fixed points at large N . However, the precise value of N corresponding to each of these kinks could be inferred only approximately from the value of ∆ σ . From tracking the kinks as we varied ∆ σ one can thus extract functional relations ∆ σ (∆ σ ) and ∆ ψ (∆ σ ), which we expect should hold, conjecturally, even for the GNY theories away from the large N limit. At some value of ∆ σ , the family of kinks intersects the supersymmetry line ∆ ψ = ∆ σ + 1/2, and we suspect that it is this value of ∆ σ corresponds to the N = 1 theory, which is expected to have enhanced N = 1 supersymmetry [27,28].
In order to study the GNY theories more precisely, we should make use of the information that these theories have O(N ) global symmetry. In the present work, we will therefore generalize the study of [13] to 4-point functions of Majorana fermions ψ i ψ j ψ k ψ , where ψ i transforms in the vector representation of an O(N ) global symmetry. In particular, we will use semidefinite programming methods to compute universal bounds on the leading O(N ) singlets and symmetric tensors appearing in the ψ i × ψ j OPE. As we will see, we observe a sequence of kinks in the space of allowed theories which match precisely with the GNY fixed points at large values of N . At small integer values of N , we use the locations of these kinks to make predictions for the critical exponents of these models. We additionally compute universal bounds for the stress-energy tensor and current central charges as a function of N , again finding features that coincide with the GNY models.
When making no assumptions about the spectrum, the bounds on the scaling dimensions of O(N ) singlet operators do not make contact with the GNY models. (They do show, however, some evidence for the existence of new "dead-end" CFTs, as we will discuss.) On the other hand, when we do make some plausible assumptions about the gap above the first parity odd singlet, we do find a set of kinks that correspond to expected scaling dimensions in the GNY models. Furthermore, when imposing such a gap for small values of N , a second set of kinks is visible. This is reminiscent of structure found in -expansion studies, e.g. recently discussed in [28], where there is an additional fixed point besides the GNY model which for small values of N has a possibility of becoming a unitary CFT. They are also similar to a second kink observed in [13], appearing close to the feature conjectured to coincide with the N = 1 supersymmetric Ising model. Our bootstrap results may support the existence of this second set of CFTs (which we refer to as GNY * ), though the story is currently unclear. The possible existence of GNY * theories is worth further study. This paper is organized as follows. In Section 2 we review the Gross-Neveu-Yukawa Model. In Section 3 we discuss the crossing relations applicable to fermionic correlators with O(N ) symmetry and present the theoretical set-up for bootstrap applications. In Section 4 we specifically focus on the Gross-Neveu-Yukawa Model by either studying universal bounds on scaling dimensions or by placing bounds after imposing further gaps in the spectrum. In the latter case, we comment on a second set of kinks which solely appears for small values of N . Next, in Section 5 we study an expanded set of universal bounds on scaling dimensions of operators and in Section 6 we study universal bounds on central charges. Finally, in Section 7 we discuss future directions.

Review of Gross-Neveu-Yukawa Model
While most of the bootstrap bounds presented in this paper will have universal implications for the space of all CFTs, as mentioned in the Introduction, our main focus will be to study the CFT data of the Gross-Neveu model at criticality [12]. The Gross-Neveu-Yukawa description contains N Majorana fermions ψ i and a parity-odd scalar field φ, with Lagrangian Here g and λ are coupling constants. The theory has an O(N ) global symmetry with ψ i transforming in the fundamental representation (i = 1, . . . , N ) in addition to a parity symmetry. For even values of N , this theory can be studied perturbatively in d = 4 − dimensions. It has a critical point that can be achieved by appropriately tuning the scalar mass m 2 , while parity symmetry forbids a fermionic mass term. This critical point, thought to survive down to d = 3, has been extensively studied using perturbative, analytic and numerical methods. Specifically, the model has been studied in the 1/N expansion [14][15][16][17][18]28], and in the -expansion (see e.g. [28][29][30][31][32][33] and references therein). The models have also been recently studied from the perspective of weakly-broken higher spin symmetry in [34]. The critical exponents at this fixed point have also been estimated through nonperturbative RG methods [35][36][37] and numerically, using Monte Carlo simulations [38][39][40][41][42].
Beyond serving as one of the most basic models for scalar-fermion interactions, the GNY-model and its variations frequently appear as universality classes for quantum phase transitions in condensed-matter systems with emergent Lorentz symmetry. It has been employed as a model for describing phase transitions in graphene [20,21,33], the Hubbard model on the honeycomb and π-flux lattice [19], models of time-reversal symmetry-breaking Perturbative estimates for lowest-lying scalars in the GNY model Perturbative estimates for central charges in the GNY model  [14][15][16], while the dimension of ψ (i ψ j) was computed in Appendix B of [13]. We list no corrections for dimension of the lowest operator in the O(N ) anti-symmetric representation, as we could not find any discussion of it in the literature. The most recent -expansion estimates for the scaling dimensions are available in [33], where O( 3 ) corrections are available for the dimensions associated to ψ i and φ. Prior work on the -expansion for the GNY models can be found in [28][29][30][31][32][33]. Bottom: Representations and one-loop values for the central charges for the O(N ) conserved current and for the stress-energy tensor at large N as determined in [45].
in d-wave superconductors [22,23], models of 3-dimensional gapless semiconductors [43,44] and models that exhibit emergent supersymmetry on the boundary of topological superconductors [27]. Thus, we are optimistic that our results for this model can find future applications in a number of different experimentally-interesting systems. Specifically, for those interested in the application of bootstrap results to such systems, in Table 3 we will list the critical exponents for the universality classes relevant to the metal-insulator transition for spinless fermions on the honeycomb lattice (N = 4) and for the semi-metallic to insulator transition in graphene (N = 8). For convenience, Table 3 also lists the values for the critical exponents obtained through several other methods.
In the context of this work, we consider a four-point function of fermionic operators ψ i . We can distinguish the operators appearing in the ψ i ×ψ j OPE by their O(N ) representation: they can be in either the singlet representation (S), the two-index symmetric traceless representation (T ), or the two-index anti-symmetric representation (A). The dimensions of these operators can be estimated at large N or in the 4 − expansion-see Table 1 for such estimates for a few of the lowest dimension operators. In Table 1 we also give the large-N estimates for the central charges of the stress-energy tensor and of the conserved O(N ) current, whose normalization we define in Section 6.
As briefly mentioned in the Introduction, in the 4 − expansion, besides the critical GNY model one finds an additional fixed point (GNY * ) [28]. For large values of N , the additional fixed point has a negative λ coupling in the Lagrangian (2.1) and thus it is expected that GNY * is non-unitary at large N . In fact, when computing scaling dimensions at such a critical point, one finds that for large N , there are scaling dimensions in the theory that become negative, thus violating the unitarity bound. However, as N is decreased, the scaling dimensions of such operators grow and eventually the unitarity bound may be satisfied. Thus, it is possible that at small values of N , the GNY * theories become unitary and could be detected using the conformal bootstrap. Generically, one expects that the scaling dimension of the scalar φ in the theory is lower in GNY * than in GNY, such that the operator φ 4 is relevant in GNY * , but irrelevant in GNY, and one could thus flow from GNY * in the UV to the GNY model in the IR. As we discuss below, the GNY * models may appear in bootstrap bounds on the dimension of low-lying parity-odd operators.

Crossing and Bootstrap with O(N ) Symmetry
In this section we set up the conformal bootstrap constraints for 4-point functions of Majorana fermions, ψ i , transforming in the fundamental vector of a global O(N ) symmetry of a parity-preserving 3d CFT. To keep our formulas compact, we follow the notation of [13] and contract the fermions with auxiliary commuting polarization variables s α . By explicitly imposing the O(N ) global symmetry, we will thus generalize the analysis of [13], where we focused on the case of four identical Majorana fermions.
Three 3-point functions between two fermions and a spin-operator O , have four possible tensor structures r a with independent coefficients λ a O , of which two (a = 1, 2) have even parity and two (a = 3, 4) have odd parity. 1 In addition, we can work in a basis such that structures a = 1, 2, 3 are anti-symmetric under the exchange 1 ↔ 2 if is even and symmetric if is odd, while the structure a = 4 is symmetric under 1 ↔ 2 if is even and anti-symmetric if is odd.
Here S, T and A denote O(N ) singlets, two-index symmetric traceless tensors and twoindex antisymmetric tensors. The superscript ± denotes whether the operators appearing are parity even or parity odd. So, for example, the first sum in the second line runs over {S + , even}, which are parity even O(N ) singlet operators of even spin. Note that we work in a basis such that the structures t I are symmetric under exchange 1 ↔ 3 for I = 1, 2, 3, 4 and antisymmetric for I = 5, 6, 7, 8. Additionally, when external dimensions are equal, the t 3,4,8 contributions vanish. For detailed definitions of the tensor structures r a and t I , see [13].
The crossing equation under the exchange 1 ↔ 3 can then be written as: and Similar definitions apply to T I A/B and A I A/B , which sum the contributions of O(N ) symmetric tensors and antisymmetric tensors, following the pattern of quantum numbers appearing in (3.2). The functions F ±I ab,∆, are defined as We can now exclude assumptions about the spectrum by applying a set of functionals where the sum is over all O(N ) representations R = S, T, or A. Here the sets of possible spins allowed in each representation, R and their complements R , are given by S = T = A = {all even spins} and where S I A/B,ij,∆,l , T I A/B,ij,∆,l , and A I A/B,ij,∆,l correspond to the contribution of the conformal block F I ij,∆,l to S I A/B , T I A/B , and A I A/B , respectively. Concretely, we look for a vector of functionals α I that satisfies the inequalities: where we apply the functionals to blocks corresponding to all representations R = S, T, or A. We now search for a functional α such that: For instance, if we want to bound the OPE coefficient of the stress-energy tensor, as we do in Section 6, we will have R O = S, and since the operator is parity-even, we have In this case, the canonically normalized OPE coefficients will be determined by the Ward identity for the stress-energy tensor.
Eqs. (3.9) and (3.12) then imply the inequality: where the expression on the RHS corresponds to the identity operator contribution. Finding a functional α I obeying (3.12) places an upper bound on (λ i O /λ i Ocan ) 2 . To make the bound as strong as possible, we search for an α I satisfying the relations (3.12) that minimizes − α I · V I 11,0,0 . We search for functionals satisfying either the (3.11) or (3.12) constraints by approximating the search as a semidefinite program and implementing it in the solver SDPB [7], following the steps described in [13]. A more detailed description of our SDPB implementation is presented in Appendix A, while the resulting constraints on the space of CFTs are presented below.

Bootstrapping GNY Models from Scaling Dimension Bounds
In this section we start off by presenting several selected bounds on scaling dimensions of scalars appearing in the ψ i × ψ j OPE. Let us first focus on the most interesting bounds in the context of GNY models and leave a more systematic presentation of general bounds for Section 5. Due to the intensive computation required for each individual run, in this work  We label the lowest dimension O(N )-invariant scalars in ψ i × ψ j OPE by σ (parity odd) and (parity even). As before, we use transcripts T and A for operators transforming as O(N ) symmetric traceless tensors and O(N ) antisymmetric tensors, respectively. We have seen in Section 3 that σ, σ T , σ A , , T can all appear in ψ i × ψ j OPE, but A cannot. Higher dimension scalars in a given representation are labelled by increasing number of primes (e.g., σ , σ , . . . ).
In Figure 1 we show the most general upper bounds on the dimension of the symmetric tensor σ T as a function of ∆ ψ for different values of N . For a range of ∆ ψ near 1, the bounds overlap and then depart from the shared curve at a critical value of ∆ ψ . In this sense, all of the bounds seem to have the feature of a "kink," reminiscent of the Ising model kink observed in [4]. Using the large-N results for the scaling dimension ∆ ψ and ∆ σ T , we can identify the bottom two kinks in Figure 1 as GNY models with N = 10, 20. For lower values of N , it is plausible to conjecture that the other kinks in Figure 1 correspond to GNY models as well. Using this as a conjecture, Figure 1 then gives a non-perturbative estimate of ∆ ψ and ∆ σ T in the GNY models at all N .  Figure 2: Bound on the scaling dimension of the lowest-lying parity-odd O(N ) singlet, ∆ σ , for a unitary CFT containing fermions with scaling dimension ∆ ψ , when imposing a gap above this operator up to ∆ σ > 3. Once again, we focus on the cases N = 2, 3, 4, and 10. We notice that when imposing such a gap, we observe features that are close to the values of ∆ ψ found for the GNY kinks in Figure 1. Estimates from the large-N expansion (black markers) for the scaling dimension of ∆ σ also agree well with the position of the kink for N = 10, but is inaccurate at smaller values of N . The red markers are the three-loop -expansion results for the dimensions of the GNY models after performing a Pade [2,1] approximation (see Eqs. (11)- (13) and Table II in [33]). They are reasonably close to the upper kinks in the bounds for small N . The lower kinks appear close to the three-loop -expansion estimates for the GNY * models (green hollow markers), obtained after performing a Pade [1,2] approximation, following the methods in [28] and [33]. While these second kinks are close to the -expansion estimates for small N , for N = 10 the second kink does not exist at all.
While it is satisfying to see that the GNY-model may saturate the universal bounds in the space of allowed scaling dimensions for this sector, we would want to determine the CFT data for operators in all the sectors of the low-lying spectrum. In order to learn about more of the spectrum we will however need to restrict the space of CFTs that we study. There are two possibilities in the present context: we can either assume gaps for the scaling dimension of operators in certain O(N ) representations and obtain tighter bounds, or, using the extremal functional method [47,48], we can reconstruct the spectra of theories close to the kinks seen in Figure 1. In the present work we will mostly focus on the former strategy, though we will mention some preliminary results from studying extremal functionals below.
As an example, we impose gaps in the O(N ) singlet parity-odd scalar sector of the  Table 2: Dimensions of some of the low-dimensional operators in GNY models obtained in this work. Dimensions of ψ, σ and σ T can be obtained from Figures 1 and 2 directly. Other dimensions were obtained from zeros of the extremal functionals and are thus given to less precision.
theory. We will thus assume the existence of an operator in this sector with dimension ∆ σ and assume a gap above it until the operator σ with the dimension ∆ σ . We aim here to make a plot analogous to Figure 3 of Ref. [13]. There we did not impose O(N ) symmetry and used the gap to the second parity odd scalar (which corresponds to σ T ) as a proxy for N . For the each value of the gap, we observed a kink coinciding with the corresponding GNY values for ∆ ψ and ∆ σ . In the present work, with the O(N ) symmetry imposed, we have the possibility to find a stable GNY kink for a range of gap assumptions. As long as the assumed gap on ∆ σ is not larger than the correct value in the GNY model the kink can exist, and it should disappear only once we choose ∆ σ too large. Figure 2 explicitly shows the consequence of imposing a gap on ∆ σ : we plot the allowed region in the space of scaling dimensions (∆ ψ , ∆ σ ) when imposing the gap ∆ σ > 3. This assumption certainly holds for large-N GNY models, where, as shown in Table 1, ∆ σ = 3+64/(π 2 N )+. . ., but is not a priori justified for small values of N . We take it as a working hypothesis, as it has an appealing interpretation that there is only one relevant scalar in this particular sector. By imposing the gap we carve out the allowed region below the free theory revealing new smoothed out "kinks" on the boundary of the allowed region for the scaling dimension for each value of N . At small values of N , for each value, we observe the existence of two distinct kinks (besides that corresponding to the free theory, in the upper left corner). We will first comment on the association of the top set of kinks with the GNY models, and then discuss the connection between the bottom set of kinks and the GNY * models.
For the set of kinks with higher values of ∆ σ , the position of this feature once again sits near the value predicted by the large-N expansion for the GNY theory: e.g. for N = 10, indicated in Figure 2 through a black square. We conjecture that for lower values of N the kinks can be used to read off the scaling dimensions ∆ ψ and ∆ σ for the GNY-model at strong coupling. This conjecture is further supported by results from the -expansion (shown by the red shapes) whose estimates for the scaling dimensions ∆ ψ and ∆ σ are somewhat close to the top set of kinks even for small values of N .
Interestingly, the -expansion also provides a possible explanation for the bottom set of  Table 3: Anomalous dimensions, η σ = 2 ∆ σ − 1 2 and η ψ = 2(∆ ψ − 1), and correlation length exponent, 1/ν = 3 − ∆ , for N = 1, 2, 4, 8. Conformal bootstrap results for N = 1 are taken from our previous work in [13] by intersecting the curve of kink found by imposing a gap on the dimension ∆ σ with the SUSY line ∆ ψ = ∆ σ + 1/2. The bootstrap results for 1/ν are obtained from estimating ∆ using the extremal functional method. The results shown from [33] for the anomalous dimensions and for the correlation length exponent are the Padé [2,1] approximations obtained from the -expansion at three loops (see Eqs. (11)- (13) and Table  II in [33]). The results shown from [28] for the same exponents are the two-sided Padé [4,2] , Padé [4,2] and, respectively, the Padé [1,5] obtained from the two-loop 2 + and 4 − expansion (see Table 1 in [28]). The results from [32] come from the -expansion at four loops with no Padé resummation performed (see Table 2 in [32]). kinks in Figure 2: they may be alternative fixed-points, dubbed GNY * models. Specifically, the values for the scaling dimensions ∆ ψ and ∆ σ estimated by the -expansion expansion, shown by the green shapes, are close to the lower set of kinks. 2 While this is encouraging, our assumption that ∆ σ > 3 is disfavored by the -expansion results for the GNY * theories which suggest, for instance for N = 4, the estimate ∆ σ ≈ 1.9. 3 It will be important to perform a bootstrap study better tailored to the GNY * models in the future.
Thus, by imposing a minimal set of gaps, we find kinks corresponding to the GNY models that give information about some operator dimensions in those theories. However, we can potentially study all sectors of the GNY model using the extremal functional method [5,9,47,48]. To be effective, the extremal functional method requires a point on the boundary of the allowed region of CFT data that is close to the theory of interest. (For example, because the bounds on the Ising model are so strong, the extremal functional method gives good estimates for several operators in that theory [5,9].) We have performed a preliminary study of the extremal functionals in the case of GNY models. The resulting spectra are somewhat noisy, but they allow an estimate of low dimension scalar operators, which we have included in Table 2.
In Table 3, we give a comparison of critical exponents obtained from our bootstrap study as well as results from the -expansion, the functional RG method, and Monte Carlo simulations. (The relation between the anomalous dimensions η σ and η ψ , the correlation length exponent ν −1 , and the scaling dimensions computed in this paper are given by where one should set = 1 in order to extrapolate to three dimensions.) We note that our results for the anomalous dimensions are closest to some of the -expansion Padé approximations performed in [28]. However, both our bootstrap results and the perturbative results differ significantly from the Monte Carlo simulations.

Universal Bounds on Scaling Dimensions
In this section we present additional general bounds on low-lying scalars in the ψ i ×ψ j OPE. We consider scalars of different parities in various representations of O(N ). We explore the bounds at somewhat higher values of the fermion dimension ∆ ψ than in Section 4. The results exhibit some unusual jumps reminiscent of the features previously seen in the fermion bootstrap of [13] where no global symmetries were imposed.

The Lowest Dimension Parity-Odd Scalar Singlet
In Figure 3, we plot universal upper bounds on the scaling dimension of the lowest-lying parity-odd scalar singlet ∆ σ as a function of ∆ ψ , for any 3D parity-invariant CFT with a global O(N ) symmetry. The bound starts at the point (∆ ψ , ∆ σ ) = (1, 2), corresponding to the free theory with N fermions. For each value of N , the bound then increases monotonically up to the point where it intersects the horizontal line ∆ σ = 3 (corresponding to a marginal operator). At these intersection points, a sharp vertical discontinuity occurs, after which the bounds plateau at much higher values of ∆ σ . For instance, at N = 2 we find a plateau around ∆ σ ≈ 9.5, while for N = 10 we find ∆ σ ≈ 16. The intersection points occur at lower values of ∆ ψ for larger N , while the jump in ∆ σ increases with N .
From these bounds, we can at least determine the values of ∆ ψ for which the CFTs must have a relevant parity-odd singlet in the ψ i × ψ j OPE. For instance, for N = 2 we conclude that 3D parity invariant CFTs with an O(2) global symmetry that have ∆ ψ < 1.175 must Reminiscent of the jump noticed in our previous work [13] when bounding the scaling dimension of the lowest-lying parity-odd operator in a theory with one fermion, we observe a different jump for each value of N that we study. Once again, all jumps occur once the bound intersects the horizontal line on which σ is precisely marginal, ∆ σ = 3.
also have at least one relevant parity-odd singlet scalar. Conversely, all such theories that have no relevant parity-odd singlet scalar must have ∆ ψ > 1.175.
The jumps shown in Figure 3 are reminiscent of features previously encountered in studies of fermions without a continuous global symmetry [13]. In Figure 1 from [13], we observed that when ∆ ψ ≈ 1.27, the bound for the dimension of the lowest-lying parity-odd scalar jumps from ∆ σ ≈ 3 to ∆ σ ≈ 7.7. When bounding dimensions in the parity-even sector for such theories, one finds a kink at the same value of ∆ ψ with the scaling of the lowestlying parity-even scalar ∆ ≈ 5.1. As noted in [13], the observed jump is similar to ones corresponding to the 3D Ising model, when studying scalar dimension bounds using mixedcorrelators [6]. By analogy with the existence of the 3D Ising model at the location of the jumps, we therefore conjectured that there exists a "dead-end" parity-invariant CFT without any relevant scalar operators, with ∆ ψ ≈ 1.27 which would have very large anomalous dimensions in both the parity-even and odd sectors, ∆ ≈ 5.1 and 3 < ∆ σ < 7.7. It is therefore tempting to extend the line of reasoning from [6] and conjecture the existence of a family of "dead-end" theories with an O(N ) global symmetry. For each value of N , the dimension of the scalar singlet in the parity-odd scalar in such a theory would satisfy monotonically with N . In the future it will be interesting, for instance, to study the extremal spectra of these jumps and try to find further evidence for the existence of new O(N )symmetric "dead-end" CFTs.

The Lowest Dimension Parity-Even Scalar Singlet
We now focus on the parity-even singlet sector. Figure 4 shows an upper bound on the scaling dimension of the lowest-lying operator in this sector as a function of ∆ ψ . The bound starts at (∆ ψ , ∆ ) = (1, 3), 4 and increases monotonically as ∆ ψ is increased, with an inflection point that gets closer to ∆ ψ → 1 as N is increased. In [13], the jump in the parity odd sector coincided with a kink at the same value of ∆ ψ in the parity even sector. However, this is not the case for the bounds presented in Figure 4. Thus, if the conjectured family of "dead-end" CFTs from Section 5.1 were to exist, they should lie within the allowed region, for instance implying that the theories with O(2) symmetry should have ∆ < 5.7.
Interestingly, the values of ∆ ψ at which we find the inflection point for the boundary curve from Figure 4 coincidentally seems to match the value predicted by the large-N expansion for the GNY-model. However, as can be seen in top of Table 1 symmetric-traceless tensor, ∆ σ T , for a unitary CFT containing fermions with scaling dimension ∆ ψ . Once again, we focus on the cases N = 2, 3, 4, 10 and 20. As ∆ ψ → 1, the bound approaches the free theory value of ∆ σ T = 3. There are two sets of kinks. The first set of kinks correspond to the GNY-models and are located in the square on the lower left. We zoom onto this lower left square in Figure 1 and discuss the meaning of the kinks in more detail in Section 4.
for the dimension of the parity-even scalar ∆ , which is identified with the scaling dimension of φ 2 , is well within the bootstrap-determined universal bounds.
Thus, the connection between these features and the GNY models is uncertain. However, as discussed in the next section, when placing universal bounds on the scaling dimension of the parity-odd symmetric traceless scalar, the presence of the GNY model on the boundary becomes apparent.

The Lowest Dimension Parity-Odd Symmetric Traceless Scalar
In Figure 5 we show an upper bound for the scaling dimension of the lowest-lying parity-odd O(N ) symmetric traceless scalar, ∆ σ T , as a function of ∆ ψ . Once again, the bound starts at the point corresponding to the free theory with (∆ ψ , ∆ σ T ) = (1,2). For all values of N , each bound increases monotonically with ∆ ψ , encountering a first kink at points with values of ∆ ψ and ∆ σ T ranging from (∆ ψ , ∆ σ T ) ≈ (1.069, 2.45) for N = 2, to (∆ ψ , ∆ σ T ) ≈ (1.008, 2.06) for N = 20. The location of these kinks can be seen more clearly in Figure 1, where we zoom in on the region of interest from Figure 5. For large values of N we find a strong agreement between the position of the kinks and large-N expansion estimates for the values of the scaling dimensions, ∆ ψ and ∆ σ T , in GNY-theories.
As we move towards larger values of ∆ ψ we find yet another discontinuity: for all values of N we once again find a jump in the bound on ∆ σ T that occurs precisely where the bound intersects the horizontal line where σ T is precisely marginal. The value of ∆ ψ at which the jump occurs increases monotonically with N , ranging from ∆ ψ ≈ 1.21 for N = 2 to ∆ ψ ≈ 1.28 for N = 20. The range of the jump in the value of ∆ σ T decreases with N , as for N = 2 one finds ∆ σ T ≤ 8.6, while for the highest value of N (N = 20) that we have used numerics for we find ∆ σ T ≤ 7.2.
Once again following the reasoning presented in Section 5.1 as well as in our previous work [13], we are led to consider the possibility that the set of jumps follows a second family of theories, namely a second extension of the "dead-end" CFT seen when placing universal bounds on fermionic theories with no O(N ) global symmetry.
While we cannot yet check the existence of these conjectured CFTs, Figure 5 allows us to conservatively claim for what values of ∆ ψ there must be at least one relevant operator in this sector: for instance, for N = 2 we find that if ∆ ψ < 1.21 then the theory needs at least one relevant operator in the symmetric traceless representation, and conversely, if we want to study CFTs with no relevant operators in this sector, we need ∆ ψ > 1.21.

Universal Bounds on Central Charges
In this section we place lower bounds on the central charge of the O(N ) conserved current C J , and on the stress-energy tensor central charge C T . These quantities are defined as coefficients in the the two-point function of the corresponding currents, Here J µ can and T µν can denote the canonically normalized O(N ) conserved current and stress-energy tensor, defined such that a theory of free N = 1 Majorana fermions has C free J = 2 and C free T = 1. The conserved current J µ can is the spin-1, parity-even, antisymmetric operator with dimension ∆ J = 2, while the operator T µν can is the spin-2, parity-even, O(N ) singlet operator with scaling dimension ∆ T = 3.

Bounds on the O(N ) Current Central Charge
Let us first determine how C J appears in the conformal block decomposition. We start by using the Ward identity for O(N )-transformations, where J µ can,ij denotes the canonically normalized conserved current (6.1), antisymmetric in the indices i and j. From this Ward identity, we find a condition on the OPE coefficients λ a J,can of the current J µν can,ij appearing in the ψ i × ψ j OPE: Here we introduced a new basis for 3-point functions and corresponding OPE coefficients µ a : In our bootstrap setup, the normalization of operators appearing in the ψ i × ψ j OPE depends only on ∆ and and is otherwise independent of the details of the CFT we study. Thus, we need to relate the canonically normalized current to the operators we use in our setup, which have normalization fixed by The comparison with (6.1) gives: We can now rearrange the terms in the crossing equation as follows: where we suppressed the u, v-dependence and used the hat on V to indicate that we switched to a 3-point function basis corresponding to µ a 's. We now search for functionals α i satisfying conditions (3.12). Notice, however, that only µ 1 J,can is fixed by Ward identity, while µ 2 value of µ 2 J , and to avoid the costly scanning over all possible values of µ 2 J , we look for the functionals which satisfy the additional condition: We can now proceed in exactly the same way as described in Section 3 to put a lower bound on C J . Figure 6 shows universal lower bounds for the current central charge C J in the units of free field value C free J = 2, as a function of the fermion scaling dimension ∆ ψ . At the unitarity bound ∆ ψ = 1 the bound goes to C J = C free J . In other words, the bound is saturated by the theory of N free Majorana fermions for all values of N . For large values of N the bound has a sharp local minimum at scaling dimensions ∆ ψ corresponding to GNY models. The values of C J at the minimum agree with the large-N estimates of C J for N = 10, 20. This is analogous to the phenomena observed in [4,5] where the Ising model was found lying in the local minimum of the C T bound, and similarly in [10] where the same was found for O(N ) vector models under certain assumptions.
For smaller values of N we do not observe any local minima or even sharp changes in the slope of the bound. It is possible that increasing the derivative cutoff Λ in our computations (see Appendix A) would improve the bounds significantly and result in local minima even for small values of N . For large N at least, we can see that the GNY models saturate the bootstrap bound and furthermore lie at a special point (minimum) of the bound, similarly to the bounds on the scaling dimension ∆ σ T in Section 4. Note that we could try to use the functional obtained in minimization of C J to extract the spectrum of the GNY model in question. In our preliminary investigations we have found that functionals obtained from scaling dimension bounds seem to give more precise results.

Bounds on the Stress-Energy Tensor Central Charge
Using the Ward identity for translations, Once again we must relate the canonically normalized stress-energy tensor to the one we use in our bootstrap program, written in Eq. (6.6), which is equivalent to imposing that the 2-point function of the stress-energy tensor is normalized as (I µρ (x 12 )I νσ (x 12 ) + I µσ (x 12 )I νρ (x 12 )) − 1 3 η µν η ρσ . (6.12) From (6.2) and the equation above, one finds We can now put a lower bound on C T by bounding the OPE operators λ 1,2 T . We do this by isolating the contribution of the parity-even spin-2 operator with ∆ = 3 from the singlet sector in the crossing Eqs. (3.3)-(3.5), where the summation in the second term on the right-hand side now excludes the stress energy tensor and the identity operator, whose contributions we wrote separately. We now  Figure 7 shows universal lower bounds on C T as a function of ∆ ψ . For all values of N , the bound starts at (∆ ψ , C T ) = (1, N ), which is saturated by the free theory with N Majorana fermions. The bound then decreases as a function of ∆ ψ until ∆ ψ ≈ 1.465, at which point all values of C T consistent with unitarity become allowed. While the GNY-model does not saturate the bounds in Figure 7, it is important to check that the recently obtained large-N estimate for the central charge C T is in the allowed region.

Discussion
In this work we computed universal bounds on scaling dimensions and central charges of 3d parity-invariant CFTs containing fermions charged under an O(N ) symmetry. We observed a sequence of "kinks" that match the O(N ) Gross-Neveu fixed points at large N and can potentially be used to learn about the theories at small N . At larger values of the fermion dimension, we also observed a sequence of discontinuous jumps in the bounds occurring when some scalar operator dimension passes through marginality. It will be important in future work to clarify two questions: 1) Why do these jumps always coincide with an operator becoming marginal and can this be understood analytically? 2) Do these jumps coincide with new physical "dead-end" CFTs? One avenue for making progress on these questions is to more carefully study extremal spectra as one passes through these discontinuities.
We would also like to understand how to isolate the GNY (and possibly GNY*) fixed points as closed islands in the allowed space of scaling dimensions, similar to how the O(N ) vector models were isolated in [8,11]. This will likely require extending the bootstrap to mixed correlators containing both fermions and scalars, for which the needed conformal blocks were worked out in [50]. An important question is whether the GNY* fixed points are fully unitary in 3d or whether they are only approximately unitary. If the GNY* kinks could be turned into closed islands then a way to probe this question would be to push the bootstrap to higher derivative order and check if the islands eventually disappear. Along these lines it will also be interesting to isolate and learn more about the 3d N = 1 supersymmetric extension of the Ising model, and to find ways to probe important variants of the GNY models with multiple scalar order parameters (e.g., the variant with 3 scalars is connected to the Hubbard model on a honeycomb lattice [51]) or additional supersymmetry.
Finally, it could be interesting to relax the assumption of parity and probe 3d fermionic CFTs with parity violation. It should also be possible to perform a similar study of the fermion bootstrap in 4d, making contact with BSM ideas related to partial compositeness [52][53][54]. We hope that our study helps to illustrate that the spinning bootstrap, even with additional global symmetry structure, is currently viable and that many more nontrivial results may be attained by studying other external spinning operators such as fermions in other global symmetry representations, global symmetry currents, stress-energy tensors, and more.  for a, b ∈ {1, 2} or (a, b) = (3, 3), (a, b) = (4, 4). In our applications, we take the operator O 0 used for the normalization of the functionals to be either the identity operator or the stressenergy tensor. Note that because of the multiplication of crossing equation by (z −z) 5 , some of the constraints in (A.4) are identically zero, or their linear combinations are identically zero, i.e. the set of constraints is not linearly independent. This can cause instabilities in SDPB, making it run indefinitely. We want to remove such "flat directions" and give only linearly independent constraints to SDPB. This can be done numerically. We can view the set of constraints (A.4) as a matrix with rows labeling the constraints and columns labeling the components of a functional, a I ± mn . We then only need to find the linearly independent rows of the matrix. That can be done for example in Mathematica using the built-in RowReduce function. Notice that this step needs to be done only once for a given Λ.
The full description of implementing the polynomial matrix program required to find a I ± mn can be found in the SDPB manual [7]. We have used a Mathematica script to manipulate the fermionic conformal blocks to obtain the matrix input for SDPB.
In order to obtain numerically accurate results we have used the parameters presented in Table 4 in our SDPB implementation. For Λ = 19 generating the input file required by SDPB takes about 30 minutes (on a single core), while solving each semi-definite program takes about 1 hour (allowed points) or 10 hours (disallowed points) on an 16 core machine.