Bootstrapping the 3d Ising model at finite temperature

We estimate thermal one-point functions in the 3d Ising CFT using the operator product expansion (OPE) and the Kubo-Martin-Schwinger (KMS) condition. Several operator dimensions and OPE coefficients of the theory are known from the numerical bootstrap for flat-space four-point functions. Taking this data as input, we use a thermal Lorentzian inversion formula to compute thermal one-point coefficients of the first few Regge trajectories in terms of a small number of unknown parameters. We approximately determine the unknown parameters by imposing the KMS condition on the two-point functions 〈σσ〉 and 〈ϵϵ〉. As a result, we estimate the one-point functions of the lowest-dimension ℤ2-even scalar ϵ and the stress energy tensor Tμν. Our result for 〈σσ〉 at finite-temperature agrees with Monte-Carlo simulations within a few percent, inside the radius of convergence of the OPE.


Introduction
In [1], we initiated a study of conformal field theories at finite (i.e. nonzero) temperature in d > 2 dimensions, using techniques from the conformal bootstrap. At finite temperature, the operator product expansion (OPE) can still be used to reduce n-point correlators to sums of n−1-point correlators. However, an important new ingredient at temperature T = 1/β is that non-unit operators can have nonzero one-point functions O β . For example, the thermal one-point function of the stress tensor T µν β encodes the free-energy density.

JHEP12(2019)072
Thermal one-point functions are constrained by a type of "crossing-equation" first written down by El-Showk and Papadodimas [2]. They noted that the Kubo-Martin-Schwinger (KMS) condition for thermal two-point functions is not manifestly consistent with the OPE, and this leads to constraints on CFT data. An efficient way to study these constraints is to use the thermal Lorentzian inversion formula developed in [1], which is an analog of Caron-Huot's Lorentzian inversion formula for zero-temperature four-point functions [3][4][5].
In this work, we apply these ideas to estimate thermal one-and two-point functions in a strongly-coupled conformal field theory in d = 3 dimensions: the 3d Ising CFT. Physically, this theory describes the 2+1-dimensional quantum transverse field Ising model at nonzero temperature, and the 3-dimensional statistical Ising model with a periodic direction of length β (both at criticality). 1 Besides its physical interest, an advantage of studying the 3d Ising CFT is that we can leverage a wealth of information about its zero-temperature OPE data from the conformal bootstrap [6][7][8][9][10]. The 3d Ising CFT is a case where Monte-Carlo (MC) techniques are also very efficient for computing some finitetemperature observables [11]. However, we believe it is worthwhile to develop bootstrapbased approaches. One might hope to eventually apply these approaches to theories that are more difficult to study with MC, like fermionic theories, or non-Lagrangian CFTs.
The thermal crossing equation of El-Showk and Papadodimas is difficult to study for two reasons. Firstly, it does not enjoy the positivity conditions that are important for rigorous numerical bootstrap techniques to work [12][13][14][15][16][17]. Thus, we will not be able to compute rigorous bounds on thermal data and will have to content ourselves with estimates. Our rough strategy is to truncate the thermal crossing equation and approximate it by a finite set of linear equations for a finite set of variables. In spirit, this is similar to the "severe truncation" method initiated by Gliozzi [18,19] and applied with some success in the boundary/defect bootstrap [20][21][22][23][24][25][26][27].
However, a second difficulty is that the thermal crossing equation converges more slowly than the crossing equation for flat-space four-point functions. Thus, naïve "severe truncation" is doomed to fail, and we need a more sophisticated approach. We will use the thermal Lorentzian inversion formula and large-spin perturbation theory to estimate the behavior of a few families of operators (specifically, the first few Regge trajectories) in terms of a small number of unknown parameters. This reduces the number of unknowns in the crossing equations and allows them to be solved approximately by a least-squares fit.
In section 2, we review the conformal bootstrap at finite temperature, following [1], together with some features of the spectrum of the 3d Ising CFT [10] that play an important role in our calculation. In section 3, we outline our overall strategy and summarize the results. As a check, we perform an MC simulation of the 3d critical Ising model and find agreement with our determination of σσ β to within statistical error, inside the regime of convergence of the OPE. Section 4 presents the details of our bootstrap-based calculation. The most complicated step is the estimation of thermal one-point coefficients for subleading Regge trajectories, which we perform by adapting the "twist-Hamiltonian" procedure of [10].

The thermal bootstrap
A CFT at nonzero temperature T can equivalently be thought of as living on the space S 1 β × R d−1 , where β = 1/T is the length of the thermal circle. This space is conformally flat, so one can compute finite-temperature correlators using the OPE, just as in flat space. However, an important difference compared to flat-space is that the thermal circle introduces a scale, and as a result operators can have nonzero one-point functions. Symmetries imply that the only operators with nonzero one-point functions are primary even-spin traceless symmetric tensors O µ 1 ···µ J . For such operators, we have where ∆ is the dimension of O, e µ is a unit vector in the S 1 direction, and b O is a dynamical constant. Consider a two-point function of a real scalar primary φ at finite temperature: Here, we introduced coordinates x = (τ, x), where τ ∈ [0, β) and x ∈ R d−1 . Assuming |x| = (τ 2 + x 2 ) 1/2 < β, this two-point function can be evaluated using the OPE: Here, O runs over primary operators appearing in the φ × φ OPE, with OPE coefficients f φφO . ∆ is the scaling dimension of O, J is its spin, and ν = (d − 2)/2. We call each term in (2.3) a "thermal block." The thermal one-point coefficient b O is defined in (2.1), and we have defined the thermal coefficients a φφ O for later convenience. For simplicity, we set β = 1 in what follows. Let us use d − 1-dimensional rotational invariance to set x = (x, 0, . . . , 0) ∈ R d−1 and introduce the coordinates (2.4) Note that z, z are complex conjugates in Euclidean signature. The two-point function g(τ, x) is invariant under τ → 1−τ . In the language of thermal physics, this is the KMS condition, and it is furthermore obvious from the geometry of S 1 β × R d−1 . However, the OPE expansion (2.3) is not manifestly invariant under τ → 1 − τ . This leads to a nontrivial crossing equation that constrains thermal one-point functions b O in terms of scaling dimensions and OPE coefficients [2]. In terms of z and z, the crossing equation/KMS condition is (2.5) Here, we have also used that g is invariant under x → −x.

JHEP12(2019)072
The coefficients a φφ O can be encoded in a function a φφ (∆, J) that is meromorphic for ∆ in the right-half-plane, with residues of the form In [1], we showed that such a function can be obtained from a "thermal Lorentzian inversion formula" Here z, z are treated as independent real variables, which means that the integral is over a Lorentzian regime x → −ix L . The first term contains the discontinuity 8) and the functions K J and F J (w) are given by The second line in (2.7) represents additional contributions that are present when J < J 0 , where J 0 controls the behavior of the two-point function in a Regge-like regime. We argued in [1] that J 0 < 0 for the 3d Ising CFT. In this work, we assume this is true and ignore these contributions.

Large-spin perturbation theory
The thermal inversion formula (2.7) becomes particularly powerful in conjunction with the KMS condition (2.5). Let us call (2.3) the s-channel OPE, which in our new coordinates is an expansion around z = z = 0 and has the region of convergence s-channel OPE: |z|, |z| < 1 . (2.11) By the KMS condition, the two-point function admits another expansion around z = z = 1, which we call the t-channel: Its region of convergence is given by:

JHEP12(2019)072
We can insert the t-channel OPE into the inversion formula (2.7) to find expressions for thermal coefficients in the s-channel. In this way, we uncover non-trivial relations between the thermal coefficients of different operators in the theory. The integral in the inversion formula (2.7) is within the region of convergence of the t-channel OPE for 1 ≤ z < 2, but for z ≥ 2 it exits this region. Corrections to the residues of a(∆, J) coming from the region z ≥ 2 are exponentially suppressed in J. Thus, the t-channel OPE encodes the all-orders expansion in powers of 1/J for thermal onepoint coefficients.
Let us review how poles and residues of a(∆, J) arise from the thermal inversion formula. As an example, we study the poles and residues contributed by a single tchannel block. Individual t-channel blocks contribute poles at double-twist locations ∆ = 2∆ φ +2n+J [1]. A similar phenomenon occurs in the flat-space lightcone bootstrap, where individual t-channel blocks again contribute to OPE data of double-twist operators. To obtain poles at other locations, one must sum infinite families of t-channel blocks before plugging them into the inversion formula. (We will see several examples below.) Nevertheless, individual t-channel blocks provide an important example that will be a building block for later calculations.
Poles in ∆ come from the region z ∼ 0. Therefore, when computing residues one can simply replace the upper bound of the z integral with 1/z ∼ ∞. However, the range of the z integral must then be artificially restricted to z max = 2 when plugging in the t-channel expansion, in order for the z integral to fully be within the region of OPE convergence. This restriction is essentially an approximation that discards corrections that die exponentially in J.
The residues are determined by a one-dimensional integral over z. To see this, we first expand the function F J z/z in z in the inversion formula, where the coefficients q r (J) are and we have rewritten the inversion formula in terms of the quantum numbers The t-channel OPE can also be expanded in a power series in (1 − z) and (1 − z), The Here, the superscripts a φφ , (O) indicate that we are studying thermal coefficients for φφ , focusing on the contribution of the t-channel operator O. To go from the first equation to the second equation above we have performed the z integral and defined the function (2.20) Here     [10]. The "?" are associated to scalars whose affiliation with a certain operator family is not fully established. Errors in bold are rigorous. All other errors are non-rigorous.
Thus, we see that the contribution of the t-channel operator O dies at large J at a rate controlled by the half-twist h O = τ O /2. The unit operator has the lowest twist in any unitary theory, and thus gives the leading contribution at large J. A second important contribution comes from the stress tensor O = T µν , which gives a universal contribution proportional to the free energy density. In general, by including successively higher-twist contributions in the t-channel, we can build up a perturbative expansion for thermal coefficients in 1/J. We will review this large-spin perturbation theory of the thermal coefficients and detail how we use it for the 3d Ising CFT in section 4.

The 3d Ising CFT
In this work, we apply the thermal crossing equation and inversion formula to compute thermal one-point coefficients in the 3d Ising CFT. It will be crucial to incorporate as much information as possible about the known flat-space data (i.e. operator dimensions and OPE coefficients) of the theory. Indeed, our approach will be closely tailored to observed features of this data. We leave the question of how our approach can be generalized to arbitrary CFTs for future work. In this section, we review some features of the spectrum of the 3d Ising CFT that play an important role in what follows. The low-dimension spectrum of the 3d Ising CFT is summarized in table 1. The lowestdimension operator is a Z 2 -odd scalar σ with dimension ∆ σ ≈ 0.518. The lowest-dimension Z 2 -even scalar has dimension ∆ ≈ 1.412.
Some of the operators in table 1 are (conjecturally) identifiable as members of largespin families -i.e. families of operators whose twists τ = ∆ − accumulate at large spin.

JHEP12(2019)072
This identification works as follows. At asymptotically large spin, it is known that there exist "multi-twist" operators [O 1 · · · O k ] n, whose twists approach τ 1 + · · · + τ k + 2n as [10]. By analyticity in spin, all operators O with spin above the Regge intercept > 0 are expected to lie on curves τ i ( ) that are analytic in [3,4]. Here, i labels the Regge trajectory of the operator. If the trajectory associated to O approaches a multi-twist value τ 1 + · · · + τ k + 2n as → ∞, we say that O is in the family [O 1 · · · O k ] n . In practice, to identify a particular family in numerics, one computes Regge trajectories using the lightcone bootstrap [10,[28][29][30][31][32][33] or Lorentzian inversion formula [3][4][5] and observes which operators they pass through. 3,4 Numerical bootstrap methods reveal large OPE coefficients in the σ × σ and × OPEs for operators in the families [σσ] 0 , [ ] 0 , and [σσ] 1 . Certain other trajectories with comparable twist are not described to high precision by numerics, including for example [σσσσ] 0 . Instead, numerics indicate that these other families have relatively small OPE coefficients in the σ × σ and × OPEs. In this work, we make the approximation that we can ignore large-spin families other than [σσ] 0 , [ ] 0 , and [σσ] 1 . It is difficult to quantify the error associated with this approximation, since other families could potentially possess large thermal one-point coefficients that don't play a role in flat-space correlators, but do contribute to thermal correlators. Nevertheless, we will find a mostly-consistent picture. However, we also see some indications that other families (in particular [σσ ] 0 ) could be important for more precise calculations, see section 4.4.
Let us discuss the families [σσ] 0 , [ ] 0 , and [σσ] 1 in more detail. The lowest-twist family [σσ] 0 has twists ranging from 1 at = 2 to 2∆ σ = 1.036 as → ∞. They are increasing and concave-down as a function of , by Nachtmann's theorem [29,37,38]. The lowest-spin operator in the [σσ] 0 family is the spin-2 stress-tensor T µν . The next operator C µνρσ has spin-4 and controls the breaking of cubic symmetry when the Ising model is implemented on a cubic lattice [39]. The family [σσ] 0 is plotted up to spin 40 in figure 1. There we show both the numerical bootstrap predictions (dots) and the results of large-spin perturbation theory (curve), which agree to high precision [10,32]. The curve  [10] using the extremal functional method [7,35,36] and the numerical bootstrap. The curve shows the prediction of large-spin perturbation theory with only ∆ σ , ∆ , f σσ , c T taken from the numerical bootstrap. figure reproduced from [10].
. The OPE coefficients of [σσ] 0 in the σ × σ and × OPEs can also be approximated in large-spin perturbation theory and are given in [10].
The families [ ] 0 and [σσ] 1 are notable in that they experience large mixing with each other at small spins. For example, the operators [σσ] 1 have larger OPE coefficients than [ ] 0 in the × OPE for spins 25. This mixing can be described by supplementing large-spin perturbation theory with a procedure described in [10]. The resulting twists and OPE coefficients match well with estimates using the extremal functional method which is used in the numerical bootstrap to extract the spectrum of theories on the boundary of the allowed region [7]. We show the twists of the [σσ] 1 and [ ] 0 families in figure 2.

Summary of method
The thermal bootstrap for the 3d Ising CFT consists of two parts. In the first part, we compute the thermal coefficients of a truncated (but infinite) subset of the spectrum in terms of the thermal coefficients of a few operators -1, , and T , where a σσ and a σσ T are unknowns. Specifically, we use the thermal inversion formula to approximately determine the thermal coefficients of all operators in the [σσ] 0 , [σσ] 1 , and [ ] 0 families described in section 2.2. In the second part, we approximate σσ as a sum over the truncated spectrum with the thermal coefficients obtained in the first part. We determine the remaining unknowns by demanding that the KMS condition is satisfied in a region of the (z, z)-plane that is within the radius of converge of the s-channel OPE. The procedure is summarized graphically in figure 3. The initial steps are as follows: Operator Mixing Finally, we approximate the σσ correlator by the truncated OPE including the scalars 1 and , and the low-twist families [σσ] 0 , [σσ] 1 , and [ ] 0 . We solve for the final remaining unknown a σσ by imposing that the KMS condition is close to being satisfied for a sampling of z and z points in the interior of the square 0 ≤ z, z ≤ 1.

Results
In this section, before diving into the details of our computation, we summarize our results and compare to MC. To perform our computation, we must make some arbitrary choices and approximations. We enumerate them in section 3.2.3 and estimate the resulting errors. Overall, the results show robustness for a wide range of choices.

One-point functions
After using the thermal inversion formula together with the KMS condition we find that The values and errors quoted capture the deviations seen over several runs of our algorithm with different parameter choices. For comparison the results obtained from MC are Note that the errors for the above two observables in MC are much smaller than for the bootstrap. This is due in part to the difficulty of using the thermal crossing equation, JHEP12(2019)072  and also to the favorable behavior of finite-size effects when computing thermal correlators with MC, see appendix A. Improving the precision of thermal bootstrap results is clearly an important challenge for the future.
Our determinations for thermal coefficients in the three low-twist families, [σσ] 0 , [σσ] 1 and [ ] 0 are presented in figure 4. 5 Unfortunately, to our knowledge, there are no available MC results for the thermal one-point functions of such higher-spin operators. However, we can use the MC results for and T in (3.2) together with the thermal inversion formula to compare to the results obtained in our computation. Note that due to the strong contribution of the unit operator in the inversion formula, the standard deviations in the thermal coefficient of all higher-spin operators in all three families is much smaller than that for a σσ and a σσ T .

Two-point function of σ
In figure 5, we show the thermal two-point function σσ β computed using our algorithm and compare it to a MC simulation that we performed. The details of our simulation are described in appendix A. 5 We choose to present the thermal coefficients a  Note that we restrict the plot to the region of OPE convergence around x = 0 and τ = 0. Right: percentage difference between the two correlators, showing good agreement (within 5%) between the bootstrap and MC predictions. At small values of |x| 2 + τ 2 we expect the MC results to be inaccurate due to lattice-size effects. As |x| 2 + τ 2 → β, we exit the region of OPE convergence, and we expect inaccuracies in the bootstrap calculation.
Overall, we find good agreement between the bootstrap prediction and MC inside the regime of convergence of the OPE. In part, this is due to the fact that the unit operator gives a large contribution in this region, and its contribution is known very precisely from the four-point function bootstrap. However, the thermal OPE also correctly recovers other features of the two-point function. For example, large-spin families sum up to correctly reproduce the t-channel singularity as τ → ±1.
We also observe decay of the two-point function in the spatial direction x. Exponential decay of thermal two-point functions in x can established rigorously by expanding the correlator in states on R × S 1 , as explained in [1,44]. However, decay in x is not obvious from the OPE, where each term grows in magnitude in the x direction. The fact that we observe decay in x serves as a check on our calculation. At long distances, the correlator behaves as e −m th x , where m th is the thermal mass. It would be interesting to understand how to determine or bound m th using information in the OPE region. 6 Finally, in figure 6 we test how close we are to satisfying the KMS condition within the region of OPE convergence. As emphasized in the figure, within an 0.9β radius from the center of the OPE convergence region the deviation from satisfying KMS is < 2%.

Systematic errors
Our algorithm above involves a few choices of parameters. To check for robustness under different choices, we show the spread of results for the thermal coefficients in figure 4. Specifically, variations in our results are mainly due to the following choices: • As we explain in section 4.4, the mixing of families requires a set of z points. Figure 4 shows the results obtained when choosing different sets of z values which span a full order of magnitude. When considering our results for a σσ , the variation between JHEP12(2019)072 Figure 6. Evidence of how well the KMS condition is satisfied in the (τ, x) plane in the region of OPE convergence. We plot the difference of the two-point function and it's periodic image, , using the average thermal coefficients presented in figure 4. Note that towards the boundary of the region of OPE convergence, our estimates for the two-point function become worse and the KMS condition is further from being satisfied. For the range (τ, x) shown above the deviation from satisfying KMS is < 2%.
the set with the lowest values of z and those with the largest is at most ∼ 10%. As we will describe in section 4.4 and is already clear from figure 4, the error for higher spin thermal coefficients is significantly lower.
• In the final step of our algorithm, we choose a set of point in the (z, z)-plane, in the s-channel region of convergence, for which we require that the thermal two-point function satisfies the KMS condition approximately. When considering significantly different regions in the (z, z)-plane as exemplified in figure 13, the variation in a σσ is only ∼ 5%. Once again, the error associated to this effect for operators with higher spins is significantly lower.
Besides the choices of parameters presented above, there are several other systematic errors: • When requiring that the KMS condition is close to being satisfied at a wide variety of point in the (z, z)-plane, we truncate the OPE of the two-point function σσ to the three low-twist families σσ . For the ranges of points at which we attempt to impose the KMS condition, corrections to the two-point function are dominated by the contribution of the next Z 2 -even operator . 7 Considering that the flat-space numerical bootstrap estimates the scaling dimension of this operator to be ∆ ∼ 6.9, we can compare the contribution of the thermal conformal block for this operator to the total contribution of all other operators in [σσ] 0 , [σσ] 1 , and [ ] 0 to σσ . This helps us estimate the error associated with neglecting this operator and higher twist operators to be ∼ 4%.

JHEP12(2019)072
• The second largest systematic error which we expect comes from the fact that when using the inversion formula we truncate the range of integration to the t-channel region of convergence, z ≤ 2. As discussed in section 2.1 we expect that the correction from the region z ≥ 2 to the thermal coefficient of an operator with spin J is exponentially suppressed in J. However, since we use the inversion formula for J ≥ 4, one might worry that at small J this correction becomes large. To probe this we note that in the O(N )-model with N → ∞ the difference between the exact result and that extracted by inverting the OPE for an operator with J = 4 is only ∼ 2.8%. 8 • There are several systematic errors associated to the operator mixing procedure. The first is due to the truncation of the spectrum to operators of twist below a cut-off value. Since the contribution of operators with higher twist is visibly suppressed, such a truncation should only introduce a small error. The second is due to the fact that while multiple operator families serve as mixing inputs, we solely focus on [σσ] 0 , [σσ] 1 , and [ ] 0 as outputs. This assumes that, just like in the flat-space bootstrap, the thermal coefficients of these three double-twist families dominate over all other families with twists below the cut-off. While we have found this to be true for the thermal coefficients in the σσ correlator, there is one family -the multitwist family [σσ ] 0 -which has a contribution comparable to that of [ ] 0 in the correlator. While we will discuss the contribution of this family extensively in section 4.3, here we note that neglecting its contribution in the mixing procedure leads to an overall difference of ∼ 4% in the mixing results. Finally, we note that after mixing we assume that the [σσ] 1 family is smooth in h and we use a fit to estimate the thermal coefficients of the and T operators. We find that by varying this fit we introduce an overall error of ∼ 3% in the final results.

Details of the computation
In this section, we describe the details of the algorithm outlined in section 3.1. We will methodically iterate large-spin perturbation theory -working our way up in twist -to compute the thermal coefficients for the [σσ] 0 , [σσ] 1 , and [ ] 0 families.
In general, we will invert operators with h < 1 from the t-channel, meaning we will work to order for the thermal coefficients in the σσ correlator, dropping terms S c,∆ (h) with c > 1−2h σ , and analogously for with h σ replaced with h .

[σσ] 0
We begin by solving for the lowest-twist family of operators in the theory, [σσ] 0 . The most direct way to study this family is through the σσ β two-point function. Large-spin 8 Specifically, the exact result found in [1] predicts that as N → ∞, a perturbation theory instructs us to start by inverting the lowest-twist operators in the t-channel. The first few low-twist primary operators in the σ × σ OPE are [σσ] 0, + + . . . .

(4.2)
Note that [σσ] 0 operators are nearly killed by Disc and thus give smaller contributions than 1, T, . Thus, we will initially neglect them, but we will add them in later. We have singled out T from the rest of the [σσ] 0 family because it has the largest anomalous dimension of the family and gives the least suppressed contribution. Inverting the operators 1, , and T , we obtain a first approximation for a These contributions can be represented by the large-spin diagrams in figure 7. The next most significant contribution comes from the [σσ] 0 family itself. To compute their contributions, one needs to sum over the family in the t-channel before inverting, as discussed in [1]. The sum we need to do is 9 where h(h) = 2h σ + δ [σσ] 0 (h) and h = h(h) + . The sum is evaluated by expanding in small δ(h) log(1 − z), and then regulating the asymptotic parts of the h sum, as was explained in [1]. For the convenience of the interested reader, we review the method in appendix B. The result is as follows,

JHEP12(2019)072
where h 0 = 2h σ + 0 . Here, the set A m ⊂ R\Z ≥0 and the coefficients c a [f ] are determined by the large-h expansion of the summand f (h), via (B.5). 10 The coefficients c a [f ] do not depend on the finite part of the sum. The coefficients α k are computed via the formula (B.7), and depend on the details of the sum. We call the terms z a (and z a log m z) 'singular' terms, and the z k 'regular' terms. The singular terms have are characterized by having nonzero s-channel discontinuity (near z ∼ 0), while the regular terms have vanishing discontinuity.
The self-corrections of the [σσ] 0 family are determined by the k = 0 term on the right hand side. We are only interested in the leading large-h contribution; recalling that the power of h is controlled by the power of (1 − z), we need only consider the term with s = 0. Taking the leading thermal coefficients in (4.3) and summing over the [σσ] 0 family starting at spin 4, and inverting, we obtain the first iteration of their self-correction; Note that S In principle, we can iterate the self-correction indefinitely. The solution to this iteration is the fixed-point of the self-correction map. How to solve for this fixed-point was also explained in [1]. In practice, one needs to truncate to some order in the anomalous dimension expansion. Truncating to order δ 2 , the self-corrected thermal coefficients are

[σσ] 1 and [ ] 0
The next families that we solve for require more care. First, we compute the leading contributions to their thermal coefficients in the large-spin limit. Afterwards, we discuss subtleties that arise when considering finite spin members of the two families.

Tree level contributions
We start by computing the asymptotic contributions. Inverting the low-twist operators 1, , T , and the [σσ] 0 family in the σσ correlator gives 'tree-level' contributions to the thermal coefficients of the [σσ] 1 family. We can compute the contributions of 1, , and T via (2.21) as where the dots denote higher order terms in 1/h that we will drop. We can also sum over the rest of the [σσ] 0 family and compute it's contribution to the [σσ] 1 pole, similarly to how we computed the [σσ] 0 self-correction in (4.6). Their leading contribution is given by 0,∆σ (h).  with the first few coefficients given in (2.23) and (4.7). Of course, this is only a naive approximation of the [ ] 0 thermal coefficients which should only work for very large J.
The [ ] 0 family is more directly accessed in the correlator, where inverting any single operator gives direct contribution to this family. For example, inverting the low-twist operators 1, , and T in the correlator gives (4.14) Here, we labeled the thermal coefficients to indicate that they are the coefficients in the correlator. The relation between the thermal coefficients in the two correlators is given by the ratio of the OPE coefficients, While at large h, (4.11) and (4.16) provide good approximations for the thermal coefficients this will not be the case at small h. In this regime, the two families [σσ] 1 and

JHEP12(2019)072
[ ] 0 are very close together in twist, and have very large anomalous dimensions due to the operator mixing described in section 2.2. Naïvely, since the families are so close in twist, and strongly mix, we simply can't be sure how the residues are distributed between the families. More systematically, the presence of large anomalous dimensions means that the poles for the families are actually quite far from the naïve locations at h = 2h σ + 1 and 2h that were used to obtain (4.11) and (4.16). The effects that produce anomalous dimensions also produce corrections to the residues on a similar scale; since the anomalous dimensions are large at these intermediate h values, the contributions to the residue must also be similarly large. Finally, there are altogether other poles for multi-twist families near the twists of these families, which the residues could further mix with.
We need to develop an approach to estimate the correct, mixed thermal coefficients. In order to estimate the correct, mixed thermal coefficients, we thus need to take into account all the corrections mentioned above. Towards that end, we now turn to developing some required technology.

The half-inverted correlator
Each individual t-channel block contributes only double-twist poles in the s-channel. However, the physical correlator has poles at non-double-twist locations. Consequently, the sum over t-channel blocks cannot commute with the inversion integral when ∆ is near the physical poles. To see why, consider a contour integral around the location of a physical pole in ∆. This integral gives zero for every t-channel block, but is certainly nonzero for the full a(∆, J). By contrast, the sum over t-channel blocks does commute with the inversion integral when ∆ is imaginary. However, we would like to determine numerically what happens at real ∆.
To get a better numerical handle on how poles can shift, we will work with a more convenient object than a(∆, J). Let's imagine applying the inversion formula 'halfway', where we do the z integral to compute the residues, but leave the z integral -which produces the poles -undone. We want to define a generating function of the form (4.17) (Once again, we assume no contributions from the arcs of the inversion formula.) Now, instead of poles in h, we have powers z h . Furthermore, the anomalous dimension corrections to pole locations are of the form The idea is that (4.17) is almost the inverse Laplace transform in h of -almost due to the pesky factor of K J . The generating function we want should relate to a(h, h) along the lines ofã The inverse Laplace transform (4.20) can be performed in a region of h where the inversion integral commutes with the sum over t-channel blocks, and thus we expect it to have a convergent expansion in t-channel blocks. The idea of defining a "half-inverted" correlator was discussed in the four-point function case in [3,10]. The definitions (4.17) and (4.20) will agree if we make a few small modifications. Firstly, we should of absorb the factor of K J inside a(h, h), so the contour integral in (4.20) does not pick up unwanted poles. (At small enough twist h such that we are away from poles in K J , we can skip this step.) Secondly, we should reinterpret J inã(z, h) as an appropriate differential operator, J, as we will explain below. Thus, we definẽ which satisfies We callã(z, h) the half-inverted correlator. Inside half-inverted correlators, J should be thought of as the linear operator acting on the space of functions of the form z h log m z. Note that J appears inã(z, h) inside q r ( J), which are rational functions of J for each integer r. Therefore, we will need to invert J when acting on this space of functions. For brevity, let's denote |h, m ≡ z h log m z. Then, expressions such as 1 can be interpreted as the inverse of the appropriate linear operator acting on this space of functions. Inverting the operator z ∂ z , we have (4.29) With this interpretation, we can substitute J for J as we did in (4.22), and define the half-inverted correlator as an honest function of z and h satisfying (4.23).

Contributions to the half-inverted correlators σσ and
Returning to the Ising model, by half-inverting our low-twist operators 1, , T , and the [σσ] 0 family in the σσ and correlators, we obtain leading-order-in-large-h approximations to the respective half-inverted correlators σσ (z, h) and (z, h). The terms we compute include those that give the naïve [σσ] 1 and [ ] 0 thermal coefficients (4.11) and (4.16), but also include many other terms coming from the sum over the [σσ] 0 family.
We are not just limited to inverting the operators 1, , T , and the [σσ] 0 family. While we do not know enough about any of the other families in the theory to compute all of their contributions, there are a special set of contributions that we can compute. In particular, while the regular terms α k depend on such particulars of the family as anomalous dimensions and an exact sum over the thermal coefficients, the singular terms do not. The singular terms only depend on the asymptotic expansions. Furthermore, the leading contributions to the singular terms are to constant order in the anomalous dimensions, thus we can compute them without any knowledge of the anomalous dimensions. Therefore, we can essentially take a half-inverted correlator, and attempt to partially solve it in the large-h regime. Let's say that the sum over [σσ] 0 produced a term where h f is the asymptotic half-twist of a multitwist family f . We can safely say that p(h) is a part of the large-h asymptotics of the thermal coefficient of the family f . Now, the sum over the family f in the t-channel includes a term Note that we can determine the singular term c a [p(h)] without having to know about the small-h behavior of the thermal coefficients of the family f , or the anomalous dimensions δ f ! This is unlike the regular terms, which depend on knowing the small-h behavior of the thermal coefficients as well as the anomalous dimensions. Inverting the singular term in (4.31), we obtain a contribution to the half inverted correlator So, we take the half-inverted correlators σσ and computed from the contributions of 1, , T , and [σσ] 0 , and augment them with the singular terms (4.32) coming from all the asymptotics of thermal coefficients of other families that appear in them.
In fact, these singular terms are crucial, and augmenting by them is a natural thing to do. For example, in order to reproduce known anomalous dimensions from the thermal inversion formula -such as those of the [σσ] 0 family -one needs to sum over multitwist families in the t-channel [1]. The prototypical example of this process is illustrated in the thermal large-spin diagram in figure 9. Also, recovering the thermal coefficients of [σσ] 0 in requires summing over generically multitwist families that are generated in by the sum over [σσ] 0 , as illustrated in figure 10. We will now briefly review these relevant processes.

Generating anomalous dimensions in σσ
Let's illustrate how anomalous dimensions are generated for the half-inverted correlator σσ by an example. We saw in (4.12) that the sum over [σσ] 0 in σσ produced a pole for the [ ] 0 family. In particular, this means that the sum over [σσ] 0 contributes a term to the half-inverted correlator σσ (z, h). This implies that there is a term, given in (4.12), in the large-h expansion of a σσ [ ] 0 (h). Now, we would be wrong to say that this is a good approximation to the thermal coefficients at small h, but at large h, we know such a term is there. By crossing symmetry of figure 8, this term is responsible for generating the δ Let's consider the contributions of [ ] 0 to the thermal coefficients in σσ . To evaluate them, we need to analyze the t-channel sum over the family. This sum has the same form as the sum (4.5) over the [σσ] 0 family, where h = 2h + + δ [ ] 0 (h) and h 0 = 2h + 0 . One important difference is that since 2h − 2h σ / ∈ Z ≥0 , the terms with m = 0 have nonzero discontinuity and contribute to the inversion formula. So, we can consider the leading term m = 0 in the anomalous dimension expansion. Now, without knowledge of small-h values of a Half-inverting this term, we obtain the corresponding contribution This is a correction to the anomalous dimension (pole location) δ [σσ] 0 of the [σσ] 0 family. As expected, this is exactly the term in large-spin perturbation theory that produces the contribution of to the anomalous dimension through the crossing-symmetric process illustrated in figure 8. Other contributions arise from similar sums over other, potentially multi-twist families, as illustrated in figure 9.
One important point to highlight is that the contribution (4.36) above does not only produce the expected anomalous dimension, it also contributes to higher poles. The halfinversion of the term in (4.35) produces another term, contributing to the anomalous dimensions at the naïve location of the [σσ] 1 family, In principle, this is an important contribution when considering the [σσ] 1 family, and through mixing, the [ ] 0 family. The moral is that we should systematically generate these terms by iterating t-channel sums and subsequent half-inversions, rather than by putting the anomalous dimensions in by hand whenever they are known, as we will also generate other contributions. In summary, we put in the anomalous dimensions of [σσ] 0 and recover them, but also generate some additional terms for [σσ] n .

Generating [σσ] 0 in
Another important phenomenon is the generation of the [σσ] 0 thermal coefficients in . Using σσ , we already computed an expression for the [σσ] 0 thermal coefficients, which we believe to be accurate. One might be tempted to input them into by hand. As with the anomalous dimensions above, it's worthwhile to generate the [σσ] 0 thermal coefficients in systematically; similarly, contributions to the [σσ] 1 thermal coefficients in are also generated.
The process with which the [σσ] 0 thermal coefficients are generated in is depicted in figure 10. Our task boils down to looking at the singular terms arising from the sum over [σσ] 0 in , Using the ingredients in section 4.3 we can now explain the mixing procedure. We expect a given half-inverted correlator to have the exact form where the sum is over families f once again, with the thermal coefficients in each family given by a c f (h) and the exact half-twist given by h f (h). In the 3d Ising CFT we would like to truncate the sum of families to f ∈ F = {[σσ] 0 , [σσ] 1 , [ ] 0 }, which, due to their low twist, have the greatest contribution to the two correlators σσ and in the light-cone limit. We will denote these truncations g c F (z, h). At small z, g c (z, h) is dominated by the families f ∈ F , and therefore well approximated by g c F (z, h). We do not include multi-twist families such as [σσ ] and [σσσσ] in the sum over f for two reasons. The first is that they give a small numerical contribution to the flat-space fourpoint functions σσσσ , σσ , , so it is reasonable to guess that their contribution to thermal two-point functions is also small. The other reason is that we know much less about their anomalous dimensions and OPE coefficients, and thus wouldn't be able to write a suitable ansatz anyway. It will be important to better understand multi-twist operators to improve our techniques in the future.
The thermal coefficients appearing in different correlators are related by ratios of OPE coefficients. For each family, let's pick a thermal coefficient a u (h) from a certain correlator that we'd like to parametrize the thermal data of that family by. Given (4.43) Accordingly, we have We can now understand eq. (4.40) as an approximation to the contribution of the families correlator,g Note that at large h, due to the decrease in the anomalous dimensions for all three families in F , the terms (a c f ) naive (h) appearing in (4.40) are close to the correct thermal coefficients appearing in (4.42). However, at small values of h, as has been described in section 2.2, the anomalous dimensions of operators in the [σσ] 1 and [ ] 0 become large and thus there is a large z-power mismatch between the terms which (a c f ) naive (h) in (4.40) and those that include a c f (h) in (4.41). Thus, all the terms in the naive expansion (4.40) will mix and contribute to the accurate thermal coefficients for all three families in F . As previously mentioned, this effect is especially noticeable on families such as [σσ] 1 and [ ] 0 whose twists are close and whose naive contribution in (4.40) are difficult to distinguish at small h. For this reason, we will refer to (4.45) as the mixing equation.
In solving for the mixed coefficients a u (h) we have conveniently written (4.45) in matrix form. Thus, for each value of h that we are interested in, we can treat the mixing equation as an over-determined linear system. Concretely, we can impose that (4.45) be satisfied for several values of z from some set of values P mix . Of course, due to the truncation of the expansion (4.40), we get an overdetermined system of equations and it is impossible to satisfy the mixing equation for all values of z. However, as one can see from figure 4, when choosing,  figure 4). 12 Thus, we solve for each term proportional to each unknown a c 1, , T in a u (h) using the method of least squares 12 This remains true as long as z O(1)e −1/δ O , where δO is the average anomalous dimension at a certain value of h for the three operator families that we are considering. for each value of h. 13 To exemplify our procedure, in figure 11, we show how the coefficients multiplying a . We thus extrapolate our results for the thermal coefficients of the [σσ] 0 family down to J = 2 (see figure 12). After mixing, the thermal coefficient of T is computed in terms of the unknowns as  T ; both belong to the [σσ] 1 family, whereas the [ ] 0 family has no such operators [10]. Therefore, our mixing procedure does not work for these operators. However, it's crucial to estimate the thermal coefficients of and T for solving the KMS condition. We have found it best to extract the thermal coefficients of the low-spin members of the [σσ] 1 family by extrapolating the mixed thermal coefficients down to small h by a simple fit. This is motivated by results from the flat-space data where the OPE coefficients and anomalous dimensions of these two operators appear to lie on smooth curves with all other members of the [σσ] 1 family. The estimates for a σσ and a σσ T obtained by performing such a fit can be extrapolated using figure 11.

Solving for b O
Finally, we will input the thermal coefficients we've obtained for the three families [σσ] 0 , [σσ] 1 , and [ ] 0 into the σσ correlator, and impose the KMS condition to determine the last unknown a σσ . We do this as follows. We evaluate the correlator minus its image under crossing in various regions of the (z, z) plane, P KMS . To determine a σσ , we attempt to minimize: (4.49) By setting ∂Λ KM S (a σσ )/∂a σσ = 0 we can finally determine the results obtained in (3.1).
The thermal inversion formula guarantees that the KMS condition is satisfied in the proximity of the point (z, z) = (0, 1). Thus, if one tries to approximately impose KMS solely in that region, there would be an almost flat direction associated to the unknown a σσ JHEP12(2019)072 and, consequently, our numerical estimates would be inaccurate. However, if one imposes KMS in a region where the OPE does not converge well the results would once again be inaccurate. Thus, we try to impose that KMS is approximately satisfied in an intermediate region and check for robustness under changes of P KMS within this intermediate regime.
We find that our results are indeed robust for various choices of the (z, z) region P KMS and, as mentioned before, for the choice of z values P mix which are used to perform the mixing of the three families. To emphasize this, in figure 4, we show a spread of the thermal coefficients obtained by minimizing (4.49) for the values of P mix in (4.46) and for values of P KMS raging between the two regions showed in (13). While the value of a σσ varies by at most ∼ 10% between any two choices of P mix and P KMS , the thermal coefficients for all other operators exhibit a much lower variance. 15 For instance, the stress energy tensor thermal coefficient varies by ∼ 5%, while the thermal coefficient of the spin-4 operator [σσ] 0, =4 varies by ∼ 1%. To test how well the crossing equation is satisfied on the Euclidean thermal cylinder we plot the difference δg KMS (τ, x) = g(x, 1 + τ ) − g(x, τ ) , (4.50) in figure 6. The KMS condition is very close to being satisfied in the regime in which both the points (x, τ ) and (x, 1 + τ ) are close to the origin of the s-channel OPE, (0, 0). For instance, we find that δg KMS (−1/4, 1/4)/g(−1/4, 1/4) = 0.0037. This shows the great extent through which one could use the thermal inversion formula to systematically solve the KMS condition or, equivalently, solve the "crossing-equation" of El-Showk and Papadodimas [2].

JHEP12(2019)072
As usual, there are three main sources of error: statistical error, finite-size effects (IR), and lattice-size effects (UV). One of the nice properties of thermal correlators is that finitesize effects are much easier to control than for flat-space correlators. The reason is that we can imagine dimensionally reducing our system along the thermal circle. The result is a theory with thermal mass m th ∼ 1/β, and consequently fluctuations in the noncompact directions die off like e −x/β . Thus, on a torus with lengths β × L × L, we expect corrections from the finiteness of L to be suppressed by e −L/β ∼ 4 × 10 −6 . By contrast, to compute flat-space two-point functions, one must consider torii with size L × L × L. In that case, finite-size effects go like (L/x) −∆ O , where O is the leading operator appearing in the OPE.
Thus, we expect that finite-size effects are negligible. Our main sources of error are statistical (visible as jitteriness in figure 5) and lattice effects which cause the simulation to become inaccurate near the coincident point singularity.