Precision Bootstrap for the $\mathcal{N}=1$ Super-Ising Model

In this note we report an improved determination of the scaling dimensions and OPE coefficients of the minimal supersymmetric extension of the 3d Ising model using the conformal bootstrap. We also show how this data can be used as input to the Lorentzian inversion formula, finding good agreement between analytic calculations and numerical extremal spectra once mixing effects are resolved.


Introduction
The modern conformal bootstrap has yielded the world's most precise critical exponents in the 3d Ising model [1]. Another theory that appears to be solvable by bootstrap methods is the 3d N = 1 supersymmetric extension of the Ising model, or the N = 1 super-Ising model [2,3]. The critical exponents of this theory can also be determined to high precision using bootstrap methods. In this paper, we continue the study of this model both numerically and analytically.
This model is also of relevance to condensed matter physics. In particular, it was argued in [4] that the corresponding superconformal fixed point can be realized as a quantum critical point at the boundary of a topological superconductor. The N = 1 super-Ising model contains only one relevant operator that is invariant under time-reversal symmetry. Physically, such a property means that non-supersymmetric renormalization groups flows can reach the fixed point by just tuning the coupling of the corresponding operator to zero. This property is called "emergent supersymmetry" and is critical for experimental realization.
In the first part of the paper, we push the numerical calculation of the critical exponents of the N = 1 Ising model to higher precision. Determining the critical exponents of the N = 1 super-Ising model is desirable for many reasons. This model is the first entry of a family of models called Gross-Neveu-Yukawa models, whose Lagrangian is These models and their variations have many interesting applications in condensed matter physics. In particular, a model with N = 8 fermions describes the quantum phase transition from the semi-metal phase to the charge density instability phase in graphene [5]. There is a large literature on theoretical studies of these models using perturbative methods (see for example [6][7][8][9][10][11]), quantum Monte Carlo simulations [12][13][14][15][16][17], and the (non-supersymmetric) numerical bootstrap [18][19][20]. Reproducing the precise value of the critical exponents from the superconformal numerical bootstrap will be an important consistency check on other methods. Knowing the precise critical exponents of the N =1 model also allows one to perform a two-sided Padé approximation of the large-N perturbative series [8][9][10] and improves the predictions of critical exponents for models with higher N . See e.g. [2].
The development of the numerical bootstrap program has been complemented by recent progress in the analytic bootstrap [21][22][23][24][25][26]. In particular, the Lorentzian inversion formula [26,27] allows one to make precise predictions for the conformal data of operators that belong to the low twist Regge trajectories, using the scaling dimensions and OPE coefficients of low-dimension operators as input [28][29][30]. In the second part of this paper, we focus on developing the analytic bootstrap for the N = 1 super-Ising model. At large spin, the spectra of the N = 1 super-Ising model approaches the spectra of generalized free fields. Due to the existence of the Majorana fermion, there is operator degeneracy. For example, the operators σ∂ µ 1 . . . ∂ µ l σ and ψγ (µ 1 ∂ µ 2 . . . ∂ µ l ) ψ have the same twist and spin. In interacting theories, these two operators mix and one needs to resolve this mixing effect to obtain accurate conformal data. We show that this can be done with the help of superconformal relations, which then allow us to use the non-supersymmetric Lorentzian inversion formula to make predictions for the conformal data of the super-multiplets that belong to the leading Regge trajectories. We then compare the resulting predictions with numerical data from the extremal functional method (EFM) and find good agreement.
The paper is organized as follows. In section 2, we push the numerical calculation to very high precision and make new precision determinations of critical exponents for this model. In section 3, we develop the analytic bootstrap, first working out the generalized free theory solution for our setup, then computing the asymptotic behavior of the leading anomalous dimensions, and finally analyzing the non-perturbative predictions from the Lorentzian inversion formula. In section 4, we use the extremal function method (EFM) to extract data about the low-twist spectrum and compare with the analytic predictions.

Numerical Bootstrap
Our setup for the conformal bootstrap builds on the previous works [31,2,3] and is identical to the one described in [2], which is an application of the "long multiplet bootstrap" idea initiated in [32]. In particular, we use the 4 crossing relations arising from the correlators σσσσ and σσ , where σ and are the parity-odd and parity-even scalars contained in the leading B (0) − multiplet Σ. The second parity-odd scalar σ ∈ Σ is assumed to be isolated. All other B Under these assumptions we have used a Delaunay triangulation search [33] and the convex optimization solver sdpb [34,35] to compute islands in the {∆ σ , ∆ σ } plane at derivative orders Λ = 35, 43, 51, 59, improving on the Λ = 27 computation performed in [2]. The parameters used for these computations are described in Appendix A. At Λ = 59 we have computed a total of 92 primal points, which include a fine scan to improve the resolution of the upper-right tip.
We will now briefly comment on numerology. We first note that our result is still barely compatible with the possible relation λ /λ σσ = tan(1) [3], which corresponds to the very upper-right tip of the allowed region. Our value of ∆ σ is also compatible with the possible analytic formula ∆ σ = (Γ(5/24) − 4)/(Γ(1/3) − 2) ≈ 0.58444186, which sits in the middle of the allowed region. A sharp goal for future work is to see if we can rule out either of these possibilities. The Inverse Symbolic Calculator [36] is an interesting tool to check whether numerical expressions are compatible with analytic formulas. With higher precision determinations for the scaling dimensions, it will be very interesting to see whether one can identify an analytic expression with a physical meaning behind it.
Finally, we have computed the minimum and maximum values of the OPE coefficient λ σσ at derivative order Λ = 51, by sampling a set of 50 primal Λ = 59 points across the island. The result is that the OPE coefficient lives in the range λ σσ = 1.072125 (9) .

Supersymmetric Generalized Free Fields
Next we will work out the predictions from the analytic bootstrap for N = 1 SCFT. To begin we need to discuss the correlators of N = 1 supersymmetric generalized free fields.
As a reminder, the 4-point function of a non-supersymmetric generalized free field σσσσ can be decomposed into conformal blocks, where exchanged operators of even spin exist at dimension ∆ n, = 2∆ σ + 2n + with coefficients [37] 1 Now we will re-interpret this coefficient in a generalized free theory with N = 1 superconformal symmetry, as describing multiple degenerate contributions.
Assuming it lives with σ in a supermultiplet, = Q 2 σ, and hence couples to the same operators, we expect a decomposition of the form Here the ratios of OPE coefficients in the second line are fixed by superconformal symmetry and given in Appendix B.
We can also consider the mixed correlator σσ , which in the (12)(34) channel only contains the identity operator, thus the total non-identity contribution must vanish: Finally we can consider the ordering σ σ , which in the (12)(34) channel yields the GFF coefficient .
Taking to be even this can then be decomposed as These equations have the solution Note that the coefficients of the [B

Inversion Formula
As demonstrated in [28][29][30], one can use input from the numerical bootstrap along with the Lorentzian inversion formula [26,27] to make analytic predictions for the spectrum of the theory.
In the non-supersymmetric case, one defines a generating function whose small z behavior captures the anomalous dimension of the leading-twist trajectory. Assuming this is given by the double-twist tower [φ 1 φ 1 ] 0, , it will give a leading contribution to the generating function associated to the correlator φ 1 φ 1 φ 2 φ 2 : where (h,h) ≡ ∆− 2 , ∆+ 2 = τ 2 , τ 2 + and for convenience we will freely interchange these variables. On the other hand, the inversion formula predicts where we follow the notation of [29] for the prefactor κ 2h and function k r,s h (z). Similarly, φ 1 φ 2 φ 2 φ 1 has a generating function: Thez integrals can be straightforwardly done by expanding the conformal blocks contributing to the integrands in terms of 2d or SL 2 blocks, each of which give a contribution expressible in terms of 4 F 3 hypergeometric functions [28,29,[39][40][41][42][43][44][45][46].
Applying this formalism to an N = 1 SCFT, for σσσσ and , one should write while for σσ one should write and for σ σ we should write To a first approximation, we can expand the exponents to linear order in the anomalous dimensions and match z ∆σ and z ∆σ log(z) terms on both sides of the equations. After using the SUSY relations between the OPE coefficients and expanding the exponents to linear order, the σσσσ correlator gives the conditions while (making use of the analytic expressions for OPE coefficient ratios in Appendix B) the σσ correlator gives the conditions and similarly the correlator gives The double bracket notation on the left-hand side of these formulas denotes the averaged values of the analogous non-supersymmetric quantity and can be computed by extracting the coefficients of the z ∆σ or z ∆σ log z terms from the appropriate inversion integrals. Note that one must take great care with the asymptotic behavior of the combinations δh O ×h, which we will see shortly are not suppressed at largeh.
As discussed in [25,[28][29][30], a more refined analysis can be performed by evaluating the equations and their derivatives at finite values of z. Concretely, a strategy that works well is to solve (31) and (33)

Anomalous Dimensions
Next let us use our linear approximations to solve for the anomalous dimensions in terms of the OPE coefficients, which will help us extract their asymptotic behavior. We can use (37) and (38) to obtain the solution To a first approximation the OPE coefficients can be replaced by the GFF coefficients (22) and (25), and in general are expected to have the same asymptotic behavior at largeh, given by This behavior is required from consistent inversion of the identity operator in all correlators.
Then the leading behavior of the anomalous dimensions at largeh will be given by The first term gets a dominant contribution from σ-exchange and leads to a behavior which falls off like ∼ 1/h ∆σ , while the second term gets contributions from stress-tensor and exchange (the lowest-twist operators in the σ × σ OPE) which contribute subleading effects ∼ 1/h and ∼ 1/h ∆ .
More precisely, we have where we used Eq. (3.20) in [29] to express the integral in terms of Ω and A as defined in that paper. We conclude that we have asymptotic behavior with

Extension to Higher-Twist Trajectories
We can similarly solve for the asymptotic behavior of the higher-twist trajectories. Focusing on the twist 2∆ σ +1 trajectories, the largeh limit of the GFF coefficients satisfy the relations Using the condition then gives Similarly, gives The inversion integral then gives the dominant contribution where we have factored out the asymptotic behavior of the OPE coefficient λ 2 σ [σ ] 0, in order to make theh −∆σ suppression manifest. After evaluating (53) and (56) and simplifying we obtain the asymptotic behavior with

Double-twist Improvement
By inverting isolated blocks we can approximate the double discontinuities entering Eqs. ] 0 towers inside the generating function for σσσσ . We will start the sum over spins at an even intermediate value 0 , anticipating that we will later include operators with spins < 0 as isolated contributions. Then we have h,∆σ+p,∆σ , similar to Eq. (3.30) in [29]. We can now compute the asymptotic contribution from the sum using Eqs. (44) and (48): where we wrote the coefficient in terms of used Eq. (6.23) of [25], and dropped all terms of order z ∆σ+1 or higher. Here the resulting coefficients are given by where ψ (n) (x) is the polygamma function and H n = ψ (0) (n + 1) + γ is a harmonic number.
This yields corrections to the left-hand side of (36) and (37) of the form For , a similar computation gives Using Eqs. (44) and (48) we find a contribution of the form The coefficient inside the sum is proportional to S 0,0 ∆σ−1 (h ) so this will contribute terms of order z ∆ +∆σ−1 = z 2∆σ , as well as terms of order z ∆ and z ∆ log(z) to the generating function, giving no contribution to the twist 2∆ σ towers but requiring the appearance of operators with twist approaching 4∆ σ and 2∆ .
Finally, let us consider the generating function of the mixed correlator σσ , which receives contributions from the operators of asymptotic twist 2∆ σ + 1: In the last line we extracted the leading log(z) and regular terms using Eq. (4.47) of [25]. Note that in the first line we assumed 0 was even for concreteness, but the final formula is valid for both even and odd 0 .
Thus, the leading term on the left-hand side of (39) is given by while the correction to the left-hand side of (38) is given by

Extremal Spectrum
Using the extremal functional method [47,48,25] we can obtain numerical estimates of the higher spectrum of the theory and compare these to predictions from the analytic bootstrap. We have applied this method using the script spectrum.py [49,50] to extract extremal spectra correponding to the minima and maxima of the OPE coefficient λ σσ described in the previous section. Overall, the resulting spectra are quite unstable compared to, e.g., the Ising [25] and O(2) models [29]. We attribute this to the significant amount of operator mixing in the theory which is difficult to numerically disentangle. However, we have been able to extract a few seemingly robust features from these spectra, described below.

Leading Scalars and Low Spin
In Fig. 3 we show the locations of low-spin operators in the extremal spectra, where each bubble denotes a collection of operators appearing near that scaling dimension across the extremal spectra and the size of the bubble is proportional to the number of extremal spectra in which the operator is found.
In the B ( ) + sectors we typically find a scalar operator at dimension 4.38(1) and a spin-2 operator at dimension 3.28 (1). In the B ( ) − sectors we consistently find the next scalar operator after Σ around dimension 6.0(1) and the leading spin-2 operator at dimension 4.58 (1). In the F ( −1/2) + sectors the spectra show significant fluctuations, with the leading = 1 multiplet typically containing a scalar component with dimension ∼ 5 − 9. The extremal spectra are also consistent with a F (3/2) − stress-tensor multiplet at ∆ = 5/2, as expected for a local CFT.   Figure 3: Scaling dimensions of low-spin operators appearing in the extremal spectra, where the size of each bubble is proportional to the number of extremal spectra in which the operator is found.

Leading-Twist Trajectories
Many of the extremal spectra also show clear trajectories in the F ( ) − and B ( ) + sectors, where the component appearing in the σ × σ OPE has twist near 2∆ σ . These trajectories can be compared with the predictions from the Lorentzian inversion formula discussed in the previous section. While some of the extremal spectra only have partial trajectories and behave somewhat erratically, others are nearly complete and show a smooth behavior.
We show an example of such a trajectory in Fig. 4, 2 where one can see an excellent agreement with the Lorentzian inversion formula after resolving the mixing effects as described in the previous section. In this plot we show the result from including the exchanged operators {1, σ, } using the leading log(z) expansion, along with the results from matching at the finite value z = .05 with various choices of exchanged operators. All computations are performed using dimensional reduction and truncating the sum over p at order p max = 4.
The notation [σσ] 0 =100 2−98 and [σ ] 0 =100 0−99 means we have included the double-twist improvements described in the previous section with the specified value of 0 , along with isolated contributions in the trajectories of asymptotic twist 2∆ σ and 2∆ σ + 1 for spins below 100, obtained from the analytic solution from the previous set {1, σ, , [σσ] 0 =2 , [σ ] 0 =2 }.   − ] 1 where we used the same value as was used to compute the extremal spectrum, ∆ σ = 2.888214659. For their OPE coefficients we used the approximations These approximations can be improved in future analyses, but already one can see excellent agreement. We also see the importance of resumming the leading twist trajectories in order to obtain good predictions at lowh. In the most precise spectrum where we have included isolated contributions for spins up to 0 = 100, we see excellent agreement with the existence of a twist τ = 1 stress tensor when the Q + [F ] 0 trajectory is extrapolated down toh = 5/2, where the inversion formula gives τ ∼ 0.9996.
The corresponding predictions for the OPE coefficients (normalized to the SUSY generalized free values) are shown in Fig. 5. In order to relate the coefficients λ 2 (h) as a function ofh to the coefficients at physical spins , we must numerically solve the equation =h− τ (h) 2 for integer and include the Jacobian factor λ 2 ( ) = 1 − τ (h) 2 −1 λ 2 (h). We can see that the OPE coefficients are more sensitive than the spectrum to the precise operators included, particularly at lowh. Converting the OPE coefficient of the stress tensor into the central charge, our most precise spectrum gives the estimate C T /C free T ≈ 1.737, close to the previous estimates C T /C free T ≈ 1.684 from the numerical bootstrap [2] and in excellent agreement with C T /C free T ≈ 1.73 from a 2-sided Padé [1,1] approximation applied to the -expansion [16]. We show the corresponding CFT data in table 1, which represents our current best analytic computation of the spectrum.
Unfortunately, our extremal spectrum OPE coefficient data at Λ = 51 seemed to fluctuate by O(1) amounts and did not appear reliable. We suspect this is due to a numerical difficulty in resolving the mixing effects along with the sharing effect described in [29]. Further study will be needed in order for us to understand the most reliable way to minimize these effects at high derivative order.
We have also investigated the CFT data using the navigator method [51]. In our navigator computations, we imposed the following gaps in the spin-2 sectors: ∆ B (2) , in addition to the scalar gaps described in section 2; gaps in other channels are the corresponding unitarity bound shifted by 10 −10 . To better understand the data associated to the stress tensor, we computed a rigorous bound on its OPE coefficient across the island at Λ = 19, i.e. we maximized/minimized λ σσQ inside the island. 3 The result is that the OPE coefficient must live in the range C T /C free T = 1.68414 (14) .
Comparing with the data in table 1, we can see that the analytic calculation of the spin-2 coefficient has an error at the ∼ 1% level.
We have also extracted the extremal data from minimizing the navigator function over the island, so far computed up to Λ = 27. 4 As shown in [51], the CFT data at the minimum navigator point can give a better estimation, compared to the data at an arbitrary point inside the island. In particular, we extracted the OPE data at the minimum navigator point using the script spectrum.py [49,50]. The OPE data is included in Fig. 5, 5 and for the most part it sits between the two most precise analytic curves, giving confidence that the analytic calculations are converging towards physically sensible results.
The specific locations for final points in the navigator computations is shown in Fig. 6.
In future work it could be helpful to use the numerical bootstrap to compute and 3 We used the Σ-navigator [51], i.e. the normalization vector is given by summation of a few crossing vectors at discrete {∆ i , i } in the channels B  4 We used the same gaps and Σ-navigator as before. The specific location of the minimum point of the navigator function can vary slightly if we choose a different normalization vector for the Σ-navigator. However in practice the difference is very small. 5 Here we have dropped a few data points at large spins > 18 which started randomly fluctuating. In cases where an operator appears at the imposed gap we plot the averaged OPE coefficients in order to reduce the sharing effect. Let us also comment that the scaling dimension data overlaps very closely with the points shown in Fig. 4 so we have not shown them explicitly.  − ] 1 multiplets 6 in order to further improve our understanding of the mixing effects at low spin. It will also be important to perform a careful study of the z dependence of the analytic trajectories and explore the extrapolation down to spin 0 and the leading Regge intercepts. We can also incorporate the trajectories with asymptotic twist ∆ σ + 2 as exchanged operators in order to further improve our results. Nevertheless, it should be clear that the N = 1 super-Ising model is an excellent laboratory for further development of analytic bootstrap methods.

Conclusion
In this work, we have pushed the numerical bootstrap calculation of the critical exponents of the N = 1 supersymmetric extension of Ising model up to Λ = 59, improving on the Λ = 27 computation performed in [2]. As can be seen from Figure 1, the rate at which the size of the islands shrink as we increase Λ has became quite slow at Λ = 59. To further improve the numerical precision of the critical exponents, naïvely increasing Λ will not be the best strategy, which is quite cost intensive. One obvious thing to attempt is to study mixed correlators containing the superfields Σ and Σ . Experience from studying mixed correlators of the non-supersymmetric Ising and O(N ) vector models tells us that this strategy may lead to considerable improvements in the precision of the numerical calculation. We leave this study for future work.
We also worked out the analytic bootstrap for the N = 1 super-Ising model and obtained precise formulas for the conformal data of super-multiplets that belong to the leading Regge trajectories. In addition to giving a precise picture of the spectrum, our formulas may be useful in future attempts at combining numerical and analytic methods.
It will be interesting to generalize this work to theories with higher supersymmetries. One obvious target is the N = 2 super-Ising model studied in [52,53]. The operator mixing we have encountered here is not unique to N = 1 supersymmetric models. This phenomena is in fact quite general in superconformal field theories. For example, superconformal blocks appearing in chiral 4-point functions in N = 2 SCFTs take the form The constants a i are fixed by superconformal symmetry and are given explicitly in [52,53]. In case of the N = 2 super-Ising model, this superconformal block corresponds to long multiplets that appear in the OPE of a chiral superfield Φ and an anti-chiral superfield Φ † . The four point function of the chiral primaries can be written as  [54]. By studying the extremal spectra of the numerical bootstrap kink observed in [52,53], one might then be able to compare the numerical result and the analytic predictions. Moreover, a SUSY inversion formula was derived in [54]. (See also [55] for an earlier formula of four dimensional N = 4 superconformal field theories.) These SUSY inversion formulas allow the superconformal OPE coefficients f 2 φφ † O , instead of averaged quantities like (36), to be be obtained directly by integrating over the double discontinuity of the crossed channels' correlation functions. It will be interesting to also derive similar SUSY inversion formulas for N = 1 theories, making the solution of the operator mixing problem much easier. We leave such explorations for future work.

A Parameters
We used the following choices for the set of spins to compute the islands at each value of Λ: The sdpb parameters used in our computations of the bootstrap islands are given in     used a finite difference of 10 −20 in each argument of the navigator function. The gradient computation can be done using the approx objective program in SDPB version 2.5.