The conformal bootstrap at finite temperature

We initiate an approach to constraining conformal field theory (CFT) data at finite temperature using methods inspired by the conformal bootstrap for vacuum correlation functions. We focus on thermal one- and two-point functions of local operators on the plane. The KMS condition for thermal two-point functions is cast as a crossing equation. By studying the analyticity properties of thermal two-point functions, we derive a “thermal inversion formula” whose output is the set of thermal one-point functions for all operators appearing in a given OPE. This involves identifying a kinematic regime which is the analog of the Regge regime for four-point functions. We demonstrate the effectiveness of the inversion formula by recovering the spectrum and thermal one-point functions in mean field theory, and computing thermal one-point functions for all higher-spin currents in the critical O(N) model at leading order in 1/N. Furthermore, we develop a systematic perturbation theory for thermal data in the large spin, low-twist spectrum of any CFT. We explain how the inversion formula and KMS condition may be combined to algorithmically constrain CFTs at finite temperature. Throughout, we draw analogies to the bootstrap for vacuum four-point functions. Finally, we discuss future directions for the thermal conformal bootstrap program, emphasizing applications to various types of CFTs, including those with holographic duals.


JHEP10(2018)070
In [11], El-Showk and Papadodimas identified an interesting crossing equation for a two-point function on S 1 β × R d−1 (here β denotes the length of the S 1 ). Because this geometry is conformally flat, one can compute two-point functions using the operator product expansion (OPE), assuming the points are sufficiently close together. The new data entering this computation are thermal one-point functions. For example, the onepoint function of a scalar operator is The β dependence of O β is fixed by the scale symmetry, but the coefficient b O is not fixed by symmetry. The OPE gives an expression for a thermal two-point function that can be schematically written as: El-Showk and Papadodimas noted that (1.2) does not manifestly satisfy (1.3). Imposing the KMS condition therefore gives a nontrivial "thermal crossing equation". This constrains the b O 's in terms of the other data of the CFT, namely the OPE coefficients f φφO and dimensions ∆ O . Via the limit S 1 × S d−1 → S 1 × R d−1 , this equation can be understood as a consequence of the usual crossing symmetry for four-point functions where we sum over some of the "external" operators. As we explain in section 2, the one-point coefficients b O , together with the usual CFT data f ijk , ∆ i , determine all finite-temperature correlators. Thus, our focus will be on computing thermal one-point coefficients using nonperturbative methods. We should note however that many interesting finite-temperature observables, like e.g. the thermal mass (discussed in section 2.3), are difficult to extract from thermal one-point functions. Such observables are an even more challenging target for the future.
The thermal crossing equation is problematic for numerical bootstrap techniques because the expansion (1.2) has coefficients f φφO b O /c O that are not sign-definite. Signdefiniteness is crucial for the linear programming-based method of [12] and its generalizations [13][14][15][16]. In this sense, the thermal bootstrap is similar to the boundary bootstrap [17][18][19], defect bootstrap [20][21][22][23][24], and four-point bootstrap in non-unitary CFTs [25,26]. Our strategy will be to develop analytical approaches to the thermal crossing equation, with JHEP10(2018)070 the hope of eventually applying them (perhaps in conjunction with numerics) to CFTs whose spectrum and OPE coefficients are relatively well-understood, like e.g. the 3d Ising model [14,16,[27][28][29][30]. 5 We should note that most of our methods will apply with any choice of boundary conditions around the circle (perhaps with slight modifications). Although our focus will be on finite-temperature, one could also study supersymmetric compactifications, or compactifications with more general chemical potentials.
A general and powerful analytic bootstrap technique that can be applied to our problem is large-spin perturbation theory [30][31][32][33][34][35][36]. Large-spin perturbation theory was recently reformulated by Caron-Huot in terms of a Lorentzian inversion formula [37] (inspired by a classic formula of Froissart and Gribov in the context of S-matrix theory [38,39]). Caron-Huot's formula expresses OPE coefficients and dimensions in terms of an integral of a four-point function in a Lorentzian regime. Inserting the OPE expansion in the t-channel into the inversion formula, one obtains a systematic large-spin expansion for s-channel data. This process can be iterated to obtain further information about the solution to crossing symmetry [30,[40][41][42][43][44][45][46][47].
In section 3, we derive a Lorentzian inversion formula for thermal one-point functions as an integral of a thermal-two point function. The integral is over an interesting Regge-like Lorentzian regime that is more natural from the point of view of Kaluza-Klein compactification than finite-temperature physics. Our formula is very close to the Froissart-Gribov S-matrix formula, and in fact our derivation is almost identical. However, our result relies on some (well-motivated) assumptions about analyticity properties and asymptotics of thermal two-point functions that we discuss further in sections 3.4 and 3.5. Our formula shows that thermal one-point functions in conformal field theory are also analytic in spin, in the same way as OPE coefficients and operator dimensions.
In sections 4 and 5, we apply our inversion formula in some examples, including Mean Field Theory and the critical O(N ) model in d = 3 at large N . We also discuss some aspects of thermal correlators in general large-N theories, especially holographic CFTs with a large gap to single-trace higher-spin operators. For the O(N ) model, by studying the two-point function φ i φ i β we derive the thermal one-point functions b O for all single-trace operators O. This includes the singlet higher-spin currents, J ∼ φ i ∂ φ i , where is a positive even integer. The result, which to our knowledge is new for > 2, can be found in (5.19)- (5.21) and is reproduced here: 2 n n! ( − n + 1) n (2 − n + 1) n m n th Li +1−n (e −m th ). (1.4) where m th β = 2 log 1+ √ 5 2 is the thermal mass of the critical O(N ) model to leading order in 1/N . We have normalized b J by the square root of the norm of J . Interestingly, this result exhibits uniform transcendentality of weight +1, a feature that would be worth understanding more deeply. For = 2, the case of the stress tensor, the result matches that of Sachdev [48]. We also derive sums of thermal coefficients for scalar composite operators with dimension ∆ = 2, 4, 6, . . . . 5 Our methods share many similarities with the recent defect bootstrap analysis in [24].

JHEP10(2018)070
Together with the thermal mass, these higher-spin one-point functions have an interesting interpretation in the context of the holographic duality of the critical O(N ) model to Vasiliev higher-spin gravity in AdS 4 (see e.g. [49][50][51]). In particular, they determine the complete set of higher-spin charges of the putative black hole solution dual to the CFT state at finite temperature. Thus, we now have the full set of higher-spin gauge-invariant data necessary to check, or perhaps even construct, a candidate solution in the bulk.
In section 6, we use our inversion formula to develop large-spin perturbation theory for thermal one-point functions. This allows us to study the thermal data of arbitrary, stronglyinteracting CFTs. Crucially, thermal two-point functions have different OPE channels with overlapping regimes of validity. Inverting terms in one channel to the other relates thermal coefficients of operators in the theory in nontrivial ways: one-point functions determine terms in the large-spin expansion of other one-point functions. These relations can be posed to formulate an analytic bootstrap problem for the thermal data. The required calculations are similar to (but simpler than) those that arise in the context of vacuum four-point functions. For example, the one-point functions of low-twist operators at large spin are dominated by an analog of the double-lightcone limit, and one is interested in the discontinuity of the correlator (as opposed to the "double discontinuity" [37]) in this limit. In fact, we see that the large-spin perturbation theory of spectral and OPE data and of thermal data are intimately tied together.
As an example, we find a universal contribution to one-point functions of "doubletwist" operators [φφ] 0,J 6 [31,32,52], proportional to the free-energy density, Here, f = F/T d < 0, where F is the free energy density, c T is the stress-tensor two-point coefficient, c free is c T for a free boson, and γ(J) is the anomalous dimension of [φφ] 0,J . The OPE coefficients f φφ[φφ] 0 (J) and anomalous dimensions γ(J) can be computed using the lightcone bootstrap for vacuum four-point functions [30][31][32][33][34][35][36][37]. 7 Our precise result is that the two sides of (1.5) match to all orders in an expansion in 1/J. (Our inversion formula also produces nonperturbative corrections in J.) The leading term in brackets is due to the unit operator. The stress tensor contribution is the second term in brackets, and it falls off like 1/J d− 2 2 relative to the leading term. The dots represent similar contributions from other operators O that are suppressed by 1/J JHEP10(2018)070 crossing equations. As an example, we apply our large-spin technology to the 3d Ising model. We conclude with discussion and comments on future directions in section 7. In appendix A, we pursue the independent direction of studying the partition function on S 1 β × S d−1 , give a rough estimate of b T in the 3d Ising model from the β → 0 limit, and discuss further aspects of this limit in appendix B. The next appendices further elaborate on technical details in the main text.
2 CFTs at nonzero temperature 2.1 Low-point functions and the OPE Any CFT correlation function on R d can be computed using the operator product expansion (OPE). Beginning with an n-point function O 1 · · · O n , we recursively use the OPE to reduce it to a sum of 1-point functions, for example (2.1) Here the C ijk are differential operators, and we have suppressed the position dependence of the O i for brevity. Each time we apply the OPE we must find a pair of operators O i , O j and a sphere surrounding them such that all other operators lie outside this sphere. This is always possible for generic configurations of points in R d . Finally, by translation invariance and dimensional analysis, 8 one-point functions on R d are given by The same procedure works on any conformally-flat manifold M d , but with two additional complications. Firstly, non-unit operators can have nonzero one-point functions. Secondly, depending on the configuration of operator insertions, it may not always be possible to perform the OPE. More specifically, to compute O i × O j , we must find a sphere containing only O i and O j whose interior is flat (possibly after performing a Weyl transformation). However, the geometry of M d may make this impossible.
In this work, we study CFTs on the manifold M β = S 1 β × R d−1 . We use coordinates x = (τ, x) on S 1 β × R d−1 , where τ is periodic τ ∼ τ + β. One-point functions on M β are constrained by symmetries as follows. To begin, translation-invariance implies that descendant operators have vanishing one-point functions: (2.3) (The notation · · · β denotes a correlator on M β .) Thus, let us consider a primary operator O with dimension ∆ and SO(d) representation ρ. 8 Here, we assume unitarity, which implies that only the unit operator has scaling dimension 0.
is valid if the two operators lie inside a sphere. The largest possible sphere has diameter β, wrapping entirely around the S 1 such that it is tangent to itself. Here, we illustrate such a sphere (blue) in d = 2.
The geometry S 1 × R d−1 is clearly invariant under SO(d − 1). It also has a discrete symmetry under which τ ↔ −τ . In general, our CFT may not be parity-invariant, so to get a symmetry of the theory, we should accompany τ ↔ −τ with a reflection in one of the R d−1 directions. This combines with SO(d − 1) to give the symmetry group O(d − 1) ⊂ SO(d), where a reflection in O(d−1) also flips the sign of τ . 9 For O β to be nonzero, the restriction of ρ under O(d − 1) ⊂ SO(d) must contain the trivial representation This implies that ρ must be a symmetric traceless tensor (STT), with even spin J. Finally, the one-point function of a spin-J operator O is fixed by symmetry and dimensional analysis, up to a single dimensionless coefficient b O : Here, e µ is the unit vector in the τ -direction.
Here and in what follows, we are implicitly normalizing our correlators by the partition function, Z(β). We will find it convenient to study two-point functions, which encode the b O 's via the OPE. 10 Note that, unlike in R d , two-point functions of non-identical operators may be nonvanishing on M β . However, for simplicity, we focus on two-point functions of identical operators, The OPE is valid whenever both operators lie within a sphere whose interior is flat. In our case, the largest such sphere has diameter β: it wraps entirely around the S 1 and is tangent to itself (figure 1). The condition for both x and 0 to lie within such a sphere is In a parity-invariant theory, the rotational symmetry group would be O(d − 1) × Z2, and we would have the restriction that only parity-even operators could have nonzero one-point functions. 10 For previous discussions of the OPE for CFT two-point functions at finite temperature, see [3,53].

JHEP10(2018)070
Assuming |x| < β, we can use the OPE to obtain and f φφO is the three-point coefficient . (2.10) We often normalize O so that c O = 1. Note that because descendants have vanishing one-point functions, we need only the leading (non-derivative) term in the OPE for each multiplet. Plugging (2.5) into (2.8), the index contraction is given by a Gegenbauer polynomial, 11 is the Pochhammer symbol, and η = τ |x| . Thus, we obtain We can think of |x| ∆−2∆ φ C (ν) J (η) as a two-point conformal block on S 1 × R d−1 . 12

Free energy density
One of the most important thermal one-point coefficients is b T , associated to the stress tensor T µν . This is related to the free energy density of the thermal CFT as follows. From (2.5), the energy density is given by 13

JHEP10(2018)070
In particular, note that b T must be negative, by positivity of energy. By dimensional analysis, the free energy density F must take the form F = f T d , where f is a dimensionless quantity. Using the thermodynamic relations F = E − T S = E + T dF/dT , we find 14) The Ward identity fixes Consequently, the coefficient of T in the thermal block expansion of φφ β (2.12) is is stress tensor two-point coefficient for the free boson in ddimensions [56]. For a single free (real) scalar, b T = −2dζ(d)/S d , as can be checked by computing its free energy For the convenience of the reader, we now collect some known results for b T in various theories.
1. For the free scalar in three dimensions, we have b free  [48,57]. We will derive this from our inversion formula in section 5.1.

Two dimensions
In d = 2, S 1 β × R is conformal to the plane, so thermal correlators on the cylinder are determined by symmetry. All one-point functions vanish except for those of operators living in the Virasoro vacuum module: Likewise, two-point functions on S 1 β × R are determined via a conformal transformation as Thermal correlation functions are also naturally computed on S 1 β × S d−1 L , owing to the role of spherical slicing in the state-operator correspondence. Due to the presence of the S d−1 curvature radius L, these thermal correlators are less constrained by conformal invariance than their counterparts on S 1 β × R d−1 . However, they must obey the flat space limit (2.20) One-point functions are fixed by dimensional analysis and spherical symmetry to take the form is a function that is not determined by conformal symmetry; it obeys the boundary condition f O (0) = 1. On the other hand, employing radial quantization, where the sum runs over all local operators O (not just primaries) and Z(β) = O e −β∆ O is the partition function. It is useful to introduce one-point thermal conformal blocks on the sphere via . We have set L = 1, and introduced the left-and right-moving conformal weights and likewise for O . These blocks were recently computed in any d, for scalar O and scalar O , in [64]. Two-point functions may also be written using the OPE and a sum over states, although we refrain from showing the details here. Consistency of (2.23) with the flat space limit (2.5) can in principle be established by taking β → 0, which involves contributions from all high-energy states. We discuss further details of this limit and thermal blocks on S 1 × S d−1 in appendix B. The general lesson is that exact computation of O S 1 β ×R d−1 by passage from S 1 β × S d−1 L is challenging.

JHEP10(2018)070
Perhaps the simplest observable to compute using these methods is b T , and we explore this possibility in appendix A with the free boson and 3d Ising model as examples. The rest of this paper is devoted to developing new methods directly on S 1 β × R d−1 .

The KMS condition and crossing
Let us now review the derivation of the KMS condition. Consider a thermal two-point function in Euclidean time A(τ )B(0) β , and let us assume τ > 0. This is given by where H is the Hamiltonian. Note that convergence of the exponential factor e −τ H requires τ > 0 and convergence of the exponential factor e −(β−τ )H requires τ < β. Thus, the above expression defines the thermal two-point function for τ ∈ (0, β). From cyclicity of the trace, one immediately finds that This is the KMS condition. Taking Thus, we can further conclude that The fact that the scalar thermal two-point function is even in x is built into the conformal block decomposition (2.12). Another approach to understand (2.28) is to note that Euclidean thermal correlators are computed by a path integral on the geometry S 1 β × R d−1 , and then (2.28) is evident from the O(d − 1) symmetries of the geometry discussed in section 2.1.
Note that the thermal conformal block decomposition (2.12) can be constrained by the KMS condition (2.28) due to the lack of manifest periodicity for the thermal conformal block, in a similar way in which the four-point functions conformal blocks are not manifestly crossing-symmetric. This constraint is well-defined within the OPE radius of convergence, whenever both β/2+ τ , and β/2− τ ∈ [0, β] (2.7). Thus, in analogy to the crossing equation for vacuum four-point functions, we will interpret (2.28) as a constraint equation for all the thermal coefficients a O appearing in (2.12). This observation was made in [11]. The analog of expanding four-point functions around the crossing-symmetric point z = z = 1/2 is, using the reflection property (2.28), to enforce that = 0 for odd n and even m . This philosophy extends naturally to thermal n-point functions, which are expectation values of Euclidean time-ordered products 15 (2.30) The above representation of the correlator is valid if Re(τ 1 − τ n ) ≤ β. If τ n is the earliest time, then a similar manipulation to (2.25) using cyclicity of the trace implies In the second line we used that operators trivially commute inside the time-ordering symbol.) It follows that the thermal expectation value of a Euclidean time-ordered product is periodic in each of the τ i (since we can decrease τ i until it becomes the earliest time and then apply (2.31)). This is again obvious from the geometry. We may regard these periodicity conditions as crossing equations. In this work, we focus on the case n = 2.
(See [65] for recent study of KMS conditions for n-point functions.) While the KMS condition imposes constraints on the a O , there is an immediate obstacle to an efficient bootstrap: the OPE expansion (2.12) is linear in the OPE coefficients, nor must the a O be sign-definite. Thus, the resulting expression lacks manifest positivity. This is more analogous to the bootstrap in the presence of a conformal boundary or defect, rather than the vacuum four-point function bootstrap [24]. To proceed, we need to develop some complementary tools; this will be the content of section 3.

Away from the OPE regime
The OPE representation (2.12) comes from interpreting the two-point function g(τ, x) in radial quantization around a point in S 1 β × R d−1 . As discussed in section 2.1, this representation is only valid when the points satisfy |x| < β. Other ways of quantizing the theory give other representations with their own regimes of validity (possibly overlapping).
Perhaps the most familiar way to study thermal correlators is to quantize the theory on R d−1 -slices, where S 1 is interpreted as a Euclidean time direction. This quantization leads to expressions for thermal correlation functions like (2.25). It is also the most natural choice from the point of view of the limit Another way of quantizing the theory (that will prove useful in the next section) is to choose a direction in R d−1 , say x 1 , as the time direction. States then live on slices 15 Note that time ordering is the only sensible ordering when operators are at different Euclidean times (i.e. "imaginary" times, although here it corresponds to real τ ). This is because if the operators weren't ordered appropriately, the exponential factors e −(τ i −τ j )H would be divergent. By contrast, if some operators are at the same Euclidean time but different Lorentzian ("real") times, we can consider different orderings among those operators, and these orderings correspond to different analytic continuations of the Euclidean correlator. For example, in real time thermal physics (where τi = iti with ti ∈ R), one can study arbitrary orderings of the operators. with geometry S 1 × R d−2 . This quantization is natural if we imagine a Kaluza-Klein compactification of a d-dimensional QFT on a spatial S 1 . In the compactified theory, the momentum generator around the S 1 , which we call P KK , becomes a global U(1) symmetry with a discrete spectrum. The Hamiltonian H KK generates translations in x 1 . Explicitly, the generators are given by where the factor of −i in P KK and the minus sign in H KK come from Wick rotation as discussed in footnote 13. Because the charges are conserved, we can evaluate them at any value of x 1 . In H KK , we have chosen to subtract off the vacuum energy by hand so that it annihilates the vacuum on S 1 × R d−2 . 16 In our two-point function g(τ, x), we can use SO(d − 1)-invariance to set x = (x 1 , 0, . . . , 0) with x 1 > 0. Interpreting the correlator in KK quantization, we obtain where |0 is the ground-state on S 1 × R d−2 . Note that H KK is Hermitian and bounded from below, so the factor e −H KK x 1 leads to exponential suppression. We discuss the regime of validity of (2.33) in section 3.4. The behavior of the correlator at large x 1 (with fixed τ ) is determined by the mass gap of the compactified theory, i.e. the smallest nonzero eigenvalue of H KK which we call m th (the "thermal mass"). By dimensional analysis, m th is a constant times 1/β. It is a folk-theorem that dimensional reduction on a circle with thermal boundary conditions produces a massive theory, i.e. m th > 0. 17 Assuming this folk-theorem is true, the correlator approaches a factorized form exponentially quickly at large |x| Like the KMS condition, the decay (2.34) is not at all obvious from the OPE. In free theories, supersymmetric compactifications, or in the presence of nontrivial chemical potentials, we could have m gap = 0 and the behavior of the long-distance correlator would be different. Finally, let us note that the representation (2.33) does not use the full (d − 1)dimensional Poincare invariance of the compactified theory. To do so, we insert a complete set of states and classify them according to their (d − 1)-dimensional invariant mass and KK momentum. This leads to a version of the Källén-Lehmann spectral representation (2.35) 16 Note that the d − 1-dimensional vacuum energy density, equivalently the Casimir energy density of the CFT on a circle, is simply b T d 1 β d−1 = βF . In particular, it is negative since bT is negative by the discussion in section 2.1.1. 17 We thank Zohar Komargodski for discussions on this point. where n is the KK momentum and G F (x, m 2 ) is the Feynman propagator in (d − 1)dimensions. The decomposition (2.35) comes from going to momentum space in the compact direction, and then applying the usual Källén-Lehmann representation in each momentum sector, yielding a density of states ρ n (m 2 ) for each n. For real τ, x, the expression (2.35) is valid whenever |x| > 0, so it has an overlapping regime of validity with the OPE. It would be interesting to study the equality of these two representations.

A Lorentzian inversion formula
Inversion formulas provide an efficient way to study the operator content of vacuum fourpoint functions in flat space. The starting point is an expansion of the four-point function in a complete set of single-valued conformal partial waves, which are solutions to the conformal Casimir equations on R d . This basis is natural because physical four-point functions are single-valued in Euclidean space. The expansion also follows on general grounds from the Plancherel theorem for the conformal group SO(d + 1, 1) [66]. One may then invert the expansion using orthogonality and completeness, to extract the exchanged operator data from integrals of the four-point function.
In [37], Caron-Huot derived a remarkable inversion formula for four-point functions that involves an integral in Lorentzian signature. 18 In this section, we will derive a Lorentzian inversion formula for thermal two-point functions. However, our formula will not have a similarly clean group-theoretic interpretation as in the case of four-point functions. The reason is that C are not a complete set of solutions to a differential equation (for any choices of (∆, J)), simply because they are not single-valued functions on R d−1 × S 1 . 19 Still, we will make progress without a completeness relation by focusing on the OPE limit and writing a formula that picks out the data in this limit.

Euclidean inversion
In analogy with the conformal partial wave expansion for four-point functions, we complexify ∆, and write the thermal block expansion (2.12) as a spectral integral: For simplicity, we set β = 1 for the remainder of the paper. The full dependence on β can be restored by dimensional analysis. The function a(∆, J) should have simple poles at the physical operator dimensions, with residues proportional to the coefficients a O , 18 See [67] for another derivation. 19 One can restrict to a disc |x| < r and introduce boundary conditions at |x| = r to obtain a completeness relation. However, such boundary conditions are not satisfied by physical two-point functions. Alternatively, one can lift two-point functions to the universal cover of R d−1 × S 1 and study completeness relations on this space.

JHEP10(2018)070
We also require that a(∆, J) not grow exponentially in the right ∆-half-plane. When |x| < 1, we can close the ∆ contour to the right to encircle the poles clockwise (hence the minus sign in (3.2)) and recover the usual thermal conformal block decomposition. The position of the ∆ contour is arbitrary as long as the integral converges. We have chosen it to lie to the left of all physical poles, including the one from the unit operator.
It is simple to write an inversion formula that produces a(∆, J) from g(τ, x), by integrating against a Gegenbauer polynomial to pick out the contribution from spin J, and then Laplace transform in |x| to pick out poles in ∆, This choice of a(∆, J) is not unique, since we only demand that it have poles and residues consistent with (3.2). To obtain the correct poles in ∆, it suffices to integrate x over any neighborhood of the origin. We call (3.3) a "Euclidean inversion formula" because it involves an integral over Euclidean space. For simplicity, we have chosen to integrate over a circle with radius 1. 20 The factor N J is defined by where This standard normalization of the Gegenbauer polynomial is unfortunately singular when d = 2, so we will treat that as a special case.

Continuing to Lorentzian signature
The angular dependence of a two-point block on S 1 β × R d−1 is precisely the same as the angular dependence of a partial wave in a 2 → 2 scattering amplitude -both are given by Gegenbauer polynomials. In the case of amplitudes, the Froissart-Gribov formula [38,39] expresses partial wave coefficients as an integral of the amplitude over a Lorentzian regime of momenta. A standard derivation of the Froissart-Gribov formula (see e.g. [24,37,68]) carries over essentially unchanged to our case, where it gives a(∆, J) as an integral over a Lorentzian region in position space. Note that the Lorentzian region we find does not correspond to the usual real time dynamics at finite temperature (where τ gets complexified). Instead, one of the components of x gets complexified and plays the role of Lorentzian time.

Kinematics
Before giving the derivation, let us discuss the Lorentzian region that will appear in our formula. Using SO(d − 1) invariance, we can restrict x to a line and denote the coordinate along this line as x E . Let us introduce coordinates It will also be useful to introduce polar coordinates r and w = e iθ such that In Euclidean signature, w lies on the unit circle and represents the angle of the two operators relative to the τ -direction, and z,z are complex conjugates of each other. We will continue to Lorentzian signature by Wick-rotating x E = −ix L , so that z,z become independent real variables. In particular, the direction τ along the thermal circle remains Euclidean and retains the periodicity τ ∼ τ + β, and w is real. This configuration is best interpreted in terms of a Lorentzian theory, one of whose spatial directions has been compactified on S 1 . It is not the Lorentzian kinematics usually considered in thermal field theory, where one considers complex τ . Instead, x L plays the role of time. Poles in a(∆, J) corresponding to physical one-point functions will come from the regime z ∼ 0 orz ∼ 0, where one of the operators is following a lightlike trajectory around the thermal circle. These lightlike trajectories are depicted in figure 2.
The regime of small or large w will play an important role in what follows. In the limit of r fixed and w → ∞, say, we have τ → x L and x L → ∞. Given the periodicity τ ∼ τ + β, the separation between the operators approaches a lightlike-trajectory along the cylinder at asymptotically large x L . In terms of (z, z), this limit corresponds to z → ∞, z → 0 with zz fixed. This limit of large boost (w → 0 or w → ∞) is analogous to the Regge limit in flat space. Let us first present the derivation in d = 2, where it is particularly simple. As noted in section 2.1.2, thermal two-point functions in d = 2 are related by a Weyl transformation to flat-space two-point functions, so this analysis is not necessary. However, the discussion in this subsection will generalize to higher dimensions.
In two dimensions, Gegenbauer polynomials are given by 21 With this normalization, we have N J = π. Viewed as a cylinder two-point block, the first term comes from the exchange of a vacuum Virasoro descendant having weights (h, h) and spin J = h − h. The second term comes from the exchange of the conjugate state having weights (h, h) and spin −J. The Euclidean inversion formula becomes where the w-contour is along the unit circle as pictured in figure 3a. Note that J must be an integer in (3.9) in order for the integrand to be single-valued along the contour. Now, the crucial claim is that g(z = rw,z = rw −1 ) satisfies the following properties as a function of w: • It is analytic in the w plane away from the cuts (−∞, −1/r), (−r, 0) (0, r), and (1/r, ∞).
• Its growth at large w is bounded by a polynomial w J 0 for some constant J 0 . Similarly, by symmetry under w → w −1 , the growth at small w is bounded by w −J 0 .
We discuss these properties in the next section. For now, let us assume them and proceed with the derivation. By analogy with the Froissart-Gribov formula, we now deform the w contour away from the unit circle. We must do this separately for the two terms w J and w −J . The term w J dies as w → 0. Assuming J > J 0 , we can deform the contour towards zero for that term to obtain the contour 3b. Similarly, the term w −J dies as w → ∞, so we can deform the contour towards infinity for that term to obtain the contour 3c (again assuming J > J 0 ).
Let us focus on the w −J term, where we deform the contour towards infinity. By our analyticity assumption, we first encounter a branch cut at w = r −1 , or equivalently z = 1 (we comment on the contribution of the z = −1 branch-cut shortly). We thus obtain an integral of the discontinuity of the two-point function g(z,z) across this cut, The deformed contour for terms that behave as w J .
(c) The deformed contour for terms that behave as w −J . where Here, we have assumed that J > J 0 , so we can drop the arcs at infinity in figure 3c. If instead J ≤ J 0 , we must keep the contribution from these arcs. The arc contributions are the analogs of finite subtractions in the case of dispersion relations for amplitudes. Because g(z,z) = g(−z, −z), the branch cut from (−∞, −1/r) contributes the same as the cut from (1/r, ∞), up to a factor of (−1) J . Finally, because of symmetry under w → w −1 , the contribution from deforming the contour for w J towards the origin is the same as the contribution from deforming the contour for w −J towards infinity, giving an overall factor of 2.
Putting everything back in (3.9), we obtain We have explicitly indicated the presence of non-trivial contributions from the arcs when J ≤ J 0 . These are given by the large w region of (3.9). Their detailed form depends on the correlator in question. We will see some explicit examples in the next section.

The inversion formula in d > 2
To perform the same derivation in d > 2 dimensions, we must find the higher-dimensional analog of the decomposition cos(Jθ) = 1 2 (w J + w −J ). The role of w J will be played by the solution to the Gegenbauer differential equation that vanishes as w → 0 (for positive J). This is given by 22 (3.14) The Gegenbauer differential equation is symmetric under w → w −1 (because the equation depends only on cos(θ) = 1 2 (w + w −1 )), so the solution that vanishes as w → ∞ is F J (w −1 ). Because the Gegenbauer differential equation is second-order, the two functions F J (w ±1 ) span the space of solutions. In particular, a Gegenbauer polynomial can be expressed as a linear combination (3.15) The above representation is correct for w in the upper half-plane. Because F J (w) has cuts along (−∞, −1] and [1, ∞) and F J (w −1 ) has cuts along [−1, 1], the representation is different when w is in the lower half-plane (the phases in front of the two terms swap). Note that when w = e iθ is on the unit circle, the two terms are complex-conjugates of each other, so their sum is real.
Plugging this representation of the Gegenbauer polynomial into the Euclidean inversion formula (3.3), we can run the same contour argument as in d = 2. The measure contributes an extra factor of (z −z) 2ν , but otherwise the derivation is essentially unchanged. We find where It is easy to check that this agrees with (3.12) in d = 2 after accounting for the proper normalization of the d = 2 Gegenbauer polynomials. 22 The function FJ is proportional to BJ defined in [67] and QJ defined in [68]. In d = 3, it is proportional to a Legendre Q-function.

Comments on the Lorentzian formula
Like the Froissart-Gribov formula, our Lorentzian inversion formula (3.16) has the interesting property that it can be analytically continued in spin J. As explained e.g. in [37], analyticity in spin is a consequence of polynomial boundedness in the w → ∞ limitspecifically our assumption that the correlator does not grow faster than w ±J 0 . Because each partial wave with nonzero spin grows in this limit, boundedness is only possible if there is a delicate balance, due to analyticity, between each partial wave with J > J 0 . This state of affairs is precisely analogous to the Regge limit of vacuum four-point functions.
The integral (3.16) is over a Lorentzian regime of the two-point function. We will see shortly that poles in ∆ come from the regionz ∼ 0 where the factorz ∆ φ − ∆ 2 −ν is singular. The residues are then determined by a one-dimensional integral over z. Because poles come from infinitesimalz, we can replace the upper limit of the z integral with ∞. In other words, the locus that contributes to CFT one-point functions is τ ∼ x L (cf. (3.6)), which is a lightlike trajectory moving around the thermal circle while moving forward in "time" x L . This trajectory is pictured in figure 2.
The fact that physical data comes from an integral over z withz ∼ 0 is also true in Caron-Huot's Lorentzian inversion formula for four-point functions. However, in that case, the integral remains entirely within the regime of convergence of both the s and t-channel OPEs. An important difference in our case is that the z-integral extends outside the regime of convergence of any OPE. In our conventions, the s-channel OPE is an expansion around z =z = 0 and the t-channel OPE is an expansion around z =z = 1. Their regimes of convergence are: Our integral is within the regime of convergence of the t-channel OPE for 1 ≤ z < 2. But for z > 2 it exits this regime. Thus, we can only obtain partial information about a(∆, J) from the t-channel OPE expansion alone. However, as we will see in more detail in section 6, corrections coming from the region z > 2 are exponentially suppressed in J.
Another interesting similarity between our inversion formula and Caron-Huot's is the significance of a double lightcone limit. We will see in section 6 that a systematic expansion for thermal one-point functions in 1/J requires understanding the thermal two-point function in the regimez ∼ 0, z ∼ 1. This corresponds to a physical configuration where the second operator is approaching the first intersection of light rays from the first operator that wrap halfway around the thermal circle. In the context of four-point functions, the same regimez ∼ 0, z ∼ 1 corresponds to all four operators approaching the corners of the square (z, z) ∈ [0, 1], and is dubbed the "double lightcone" limit. 23 Because our limit plays a similar role in large-spin perturbation theory, we will adopt the same terminology. . For fixed r ∈ (0, 1), the s-channel OPE (expansion around z =z = 0) implies that the thermal two-point function g(z,z) is analytic in an annulus in the w plane between radii r and 1/r (shaded blue). The t-channel OPE (expansion around z =z = 1), together with symmetry under w ↔ −w, implies analyticity in the red-shaded regions, except for cuts running along (−∞, −1/r), (−r, 0), (0, r), (1/r, ∞) (indicated with zig-zags). In this section, we argue for analyticity everywhere in the upper and lower half planes.

Analyticity in the w-plane
To complete our derivation, we must justify the assumptions stated in section 3.2.3, namely analyticity of the two-point function in w outside the cuts pictured in figure 3, and polynomial boundedness in w. First note that convergence of the s-channel OPE guarantees analyticity in the annulus Convergence of the OPE around z,z = 1 and z,z = −1 additionally guarantees analyticity in more complicated regions around w = ±1 (figure 4).
To argue for analyticity in an even larger region, we will use the KK representation discussed in section 2.3, where Here, we have quantized the Euclidean theory on spatial slices with geometry The Hamiltonian H KK generates translations in the noncompact direction parameterized by x E , while P KK generates translations in τ (the periodic direction). In this way of quantizing the theory, both H KK and P KK are Hermitian. We first claim that g(τ, x E ) is bounded whenever Im(z) > 0 and Im(z) < 0. Our goal will be to relate a general configuration with Im(z) > 0 and Im(z) < 0 to a standard configuration where we know that the correlator is bounded. We begin by splitting the The Cauchy-Schwartz inequality implies This essentially allows us to assume that z,z are pure imaginary. Next, we claim that H KK ± P KK are bounded from below by a constant. To argue this, first let us forget about periodicity and consider a theory in R d . Then we have H ± P > 0 as a simple consequence of positivity of energy. (If either H ± P had a negative eigenvalue, Lorentz invariance would imply the existence of a state with negative energy.) Now consider the case where τ is periodic. The spectrum of P KK is quantized, with eigenvalues given by Kaluza-Klein (KK) momenta n ∈ Z. In general, energies of excitations with KK-momentum n are different from energies of excitations in R d with |p| = n. (For example, in the n = 0 sector, the lowest nonzero eigenvalue of H KK is the thermal mass m th , while the Hamiltonian in R d is gapless at zero momentum.) In particular, it is not obvious whether H KK ± P KK are positive operators. Positivity of H KK ± P KK does not follow immediately from positivity of energy because there is no Lorentz boost relating H KK and P KK . However, for sufficiently large |n|, periodicity of the τ direction becomes unimportant, and the spectrum of H KK approaches the flat-space spectrum. Thus, we expect H KK ± P KK are bounded from below for all n by some n-independent constant λ. 24 This is the key claim in this section, and we have not established it rigorously. However, we believe it is a physically reasonable assumption.
Thus, let us pick λ such that H KK ± P KK > λ, and let ζ = min(Im(z), − Im(z)). We have The correlator g(0, ζ) is simply a Euclidean correlator at nonzero |x| and time τ = 0. This is a nonsingular configuration, so the right-hand side is bounded.

JHEP10(2018)070
Finally, note that this derivation did not actually depend on φ(0) being primary. Thus, it applies to all correlators of descendants of φ, so all derivatives of g(τ, x) are bounded as well. It follows that g(τ, x) is analytic if ζ > 0, i.e. Im(z) > 0 and Im(z) < 0. These conditions hold for w in the upper half-plane. Symmetry under w → −w then implies that g(τ, x E ) is analytic in the lower w-half-plane as well.

Behavior at large w
Our derivation of the Lorentzian inversion formula relies on the assumption that g(z,z) grows no faster than w J 0 at large w (anywhere in the upper half plane) and fixed r, for some fixed J 0 . We have not been able to prove this claim or establish a rigorous upper bound on J 0 (analogous to the bound on chaos for thermal four-point functions [69]). In this section, we discuss the claim in more detail.
First, one can check explicitly that thermal two-point functions in d = 2, given in (2.18), are exponentially damped at large w. Thus, our inversion formula implies analyticity in spin for all J ≥ 0 in 2d. This is perhaps unsurprising given that only members of the Virasoro multiplet of the identity get thermal expectation values, and such operators lie on simple trajectories as a function of J.
For d > 2, we might hope to determine J 0 by studying perturbative examples. However, we should be aware that naïve perturbative expansions may not commute with the large-w limit. For example, in the critical O(N ) model at leading order in N (see section 5.1), we find J 0 = 0. However, at each order in 1/N the correlator may grow more quickly. 25 It would be very interesting to set up a perturbative analysis specially adapted to the large-w limit, perhaps analogous to [71].
Let us guess the behavior of the two-point function at large w by studying two interesting physical regimes. Firstly, consider w = iW with W large and real. This corresponds to a Lorentzian two-point function at finite temperature, in the limit where both operators are highly boosted. 26 In the absence of the thermal bath, such a correlator would be independent of W because a boost is a symmetry. The correlator can, roughly speaking, be interpreted as the amplitude for excitations created by the first operator to be absorbed by the second. There is no clear reason why this amplitude should be enhanced by the thermal bath, and thus we expect the correlator to grow no faster than W 0 in this regime. In fact, we might expect that the thermal bath destroys correlations between the operators, so the correlator actually decays at large W . 27 Another interesting physical regime is w = (1 + i )W with W large and real (i.e. w on top of one of the cuts in figure 4). This is the configuration that appears in the Lorentzian inversion formula. It corresponds to one of the operators moving on a nearly 25 Precisely this phenomenon happens for the Regge limit of four-point functions in large-N theories. One can show on general grounds that the four-point function is bounded in the Regge limit [69]. However, each order in 1/N contributes faster and faster growth in the Regge limit. In holographic theories, this fact is related to the necessity of regularization and renormalization in the bulk effective theory [70]. 26 We thank Juan Maldacena for suggesting we consider this regime. 27 We have checked that a boundary thermal two-point function computed in AdS4 using the geodesic approximation decays at large W . lightlike trajectory around the thermal circle, with one of the noncompact directions as increasing Lorentzian time x L (figure 2). A physical picture is that the first operator creates excitations that move around the thermal circle. They repeatedly collide at x L = β/2, β, 3β/2, . . . . Finally, some of them are absorbed by the second operator. We expect each collision to reduce the amplitude for excitations to reach the second operator. Thus, we conjecture that the correlator grows no faster than W 0 in this regime as well. If the collisions have a large inelastic component, the correlator should decay at large W .
These arguments are far from rigorous. It would be nice to understand -either from examples or a general argument -what the nonperturbative behavior of the two-point function can be in the entire upper half-plane at large w.

Applications I: mean field theory
In this and the next section, we will perform some checks of the inversion formula, derive some new results and demonstrate its mechanics. We begin here with application to mean field theory.
In MFT, the operators appearing in the φ × φ OPE for some scalar primary φ are the unit operator and double-twist operators of schematic form where is even, with dimensions ∆ n, = 2∆ φ + 2n + . Note that the free theory is the MFT with ∆ φ = ν, where the [φφ] 0, are identified with spin-currents J . The thermal two-point function can be computed by using the method of images, Using this, we will perform a brute-force expansion of the two-point functions into thermal conformal blocks and compare that with the thermal one-point coefficients generated by the inversion formula.
Expanding the thermal two-point function. We start by explicitly expanding the thermal two-point function without using the inversion formula in order to provide a nontrivial check for the entire methodology. Going back to the x and τ coordinates, we can write each term in (4.2) as 28 3) 28 Here, we use the identity where sgn(m) = m |m| and η = τ |x| . Thus, the two-point function is where (a) n = Γ(a+n) Γ(a) is the Pochhammer symbol. Plugging this into (4.4), and replacing j = 2n + , we get (4.6) This has precisely the form of the thermal conformal block decomposition given by (2.12), with support only on the unit operator and double-twist operators (4.1), whose one-point functions are given by In the free theory where ∆ φ = ν, the spin-currents J ≡ [φφ] 0, have Note that when d = 3, the coefficient a J 0 is divergent. This is because the zero mode is badly behaved under dimensional reduction to d = 2, which is related to the fact that the free boson in d = 2 with noncompact target space is pathological. We can now compare the above results to those predicted by the inversion formula, starting with the case d = 2 where the Gegenbauer polynomials take a simpler form.
Inversion in d = 2 MFT. As required in the inversion formula (3.16), we should be looking at discontinuities across the real z axis for each term in (4.2), Plugging (4.9) into the d = 2 inversion formula (3.12), we find Because we will be interested in poles coming from the region of infinitesimalz, we can replace the upper limit of the z integral with ∞. The z and z integrals in (4.10) then factorize and become respectively. As expected, (4.11) has poles at h = ∆ φ + n, corresponding to MFT operators (4.1). Computing the residues of each pole, we find (4.13) Note that we get an extra factor of two when we write the pole in h as a pole in ∆. The thermal one-point function is minus the residue of the pole in ∆. Thus, the thermal coefficient for double-twist operators of even spin can be read off as, which is in agreement with (4.7) when ν = 0.
Inversion in d > 2 MFT. While in d = 2 one can obtain the contribution of all doubletwist families through simple integral manipulations, the d > 2 case will require a more careful series of approximations to get the residues corresponding to each family's pole. Plugging in the discontinuity (4.9) into the inversion formula (3.16), we are left to compute Again, poles for double-twist operators (4.1) come from the region of integration nearz ∼ 0. Thus, we can replace the upper limit of the z integral with ∞. We are also free to rescalē z → zz and set the integration range forz back to [0, 1], since we will obtain the same pole location with the same residue. By also rescaling z → mz we find The z integral can be done explicitly, leaving thez integral We can now expand inz and get a series of poles. Let us focus on the first sets of poles in the integrand of (4.17) corresponding to the [φφ] 0, and [φφ] 1, operators, Multiplying this by the factor (1 + (−1) J )2 sin(π∆ φ )K J , left out in (4.15) for clarity, gives the full contribution of these poles to a(∆, J). The first term gives a pole at h = ∆ φ of the form where we've set ∆ = 2∆ φ + J in the last step to obtain the correct value of the residue. This agrees with (4.7) at n = 0. The orderz term in (4.18) gives . (4.20) This agrees with (4.7) at n = 1. Note that the unit operator pole was absent in the above manipulations. This is resolved by the presence of the "arc terms" in the Lorentzian inversion formula, a arcs (∆, J), which we have neglected here. In the limit |w| → ∞ the MFT correlator is simply given by, lim |w|→∞ g(r, w) = 1/r 2∆ φ . The contribution of the contours given by (3.10) precisely yields a pole corresponding to the unit operator at J = 0 , ∆ = 0 with residue equal to 1. We will witness a more intricate balance between the arc and non-arc contributions when studying the O(N ) vector model in the subsection below.

Applications II: large N CFTs
Consider a CFT with large central charge c T ∼ N 2 → ∞. In the OPE regime, we may organize two-point functions on S 1 β × R d−1 by powers of 1/N . Let us use the canonical normalization where O i are single-trace operators. Then the thermal scalar two-point function g(τ, x) receives the following types of contributions, organized by powers of 1/N appearing in the OPE coefficients: where we have again defined the double-trace composite operators [AB] n, , of schematic form 3) The first group of operators in (5.2) represents the two-point function of MFT, in which the [φφ] n, appear with the MFT OPE coefficients, which can be found in [72]; the second group represents single-trace operators; the third group represents double-trace operators, including the 1/N 2 corrections to the MFT exchanges; and so on. However, this way of organizing the contributions is not terribly useful because the one-point functions themselves scale with positive powers of N . In particular, in the normalization (5.1), one-point functions of n-trace operators exhibit the leading-order scaling This implies an infinite set of contributions to g(τ, x) at order N 0 , which poses an obvious challenge to computing g(τ, x), in contrast to the familiar 1/N counting used in vacuum four-point functions.
We now study the inversion formula in the critical O(N ) vector model, and discuss some features of its application to CFTs with weakly coupled holographic duals.

O(N ) vector model at large N
The critical O(N ) model at large N has been studied in detail before [48,73,74]. This theory has c T = N c free to leading order in 1/N . The main feature we will need is the value of the thermal mass. In the O(N ) model at large N , the thermal mass is equal to the expectation value of the IR operator σ, which appears in the action after applying a Hubbard-Stratanovich transformation to the φ 4 coupling: The critical point is obtained by taking λ → ∞ as σ 2 becomes irrelevant in the IR. In appendix C, we review the derivation of the following result [74,75],

JHEP10(2018)070
As we shall see later in this section, the above formula for the thermal mass is intimately related to correctly reproducing the O(N ) singlet spectrum from the inversion formula. Thermal properties of the O(N ) model were also studied in [76,77]. Let us enumerate the O(N ) singlets of the critical O(N ) model whose thermal expectation values we will compute (any non-singlet has vanishing thermal one-point function). The "single-trace" singlets are the scalar σ, with ∆ = 2 + O(1/N ), and the higher-spin currents J , with ∈ 2Z + and ∆ = + 1 + O(1/N ). In the φ i × φ i OPE, one generates the larger family of operators 29 where the anomalous dimensions are suppressed as γ n, ∼ O(1/N ). For n = 0, these operators are the slightly-broken higher-spin currents, The families (5.8) do not analytically continue down to = 0; instead, σ plays the role of φ i φ i in the IR. Accordingly, the most basic scalar operators are powers of σ, In what follows, we will compute thermal one-point functions of J and σ m , and exhibit the algorithm for computing the one-point functions of [φ i φ i ] n, for all (n, ). As discussed below (5.2), the thermal coefficients a O in large N CFTs receive contributions from an infinite set of operators in the φφ OPE, due to the opposite large N scaling of OPE coefficients f φφO and thermal one-point functions b O . That discussion was for singletrace operators φ, but the same scaling holds for the φ i fields in the O(N ) model, i.e. 30 where in order to derive the second set of scalings we have used the schematic operator relation ∂ 2 φ i ∼ σφ i . We emphasize that the computation of a σ m and a [φ i φ i ] n, (which we will show below) gives a window into arbitrarily high orders in 1/N perturbation theory: for instance, to derive f φ i φ i σ m would require going to (m − 1) th order in large-N , which is intractable using standard perturbative methods. [φiσφi] n,k, = φi∂ µ 1 . . . ∂ µ (∂ 2n σ)∂ 2k φi, with ∆ n,k, = 1 + 2n + 2k + 2 + . Note that such operators are degenerate in h andh for different values of n and k and may have degenerate dimensions with some n > 0 operators in (5.7); in the presence of degeneracies, the inversion method as presented here yields linear combinations of thermal one-point functions. 30 We assume canonical normalization for the operators φi and σ: φi(x)φj(0) = δij/|x| and σ(x)σ(0) = δij/|x| 4 .

JHEP10(2018)070
Thermal two-point function review. The propagator for the field φ i in fourier space is given by, At the saddle-point, the non-zero expectation of σ thus acts like a mass term which is absent when considering the MFT propagator considered in section 4. We can now use the G ij (ω n , k) to express the propagator in position-space as 31 This is similar to the MFT propagator (4.2), but with an exponentially decaying factor multiplying each term. While in the MFT study in section 4, each term in (4.2) could be expanded in Gegenbauer polynomials, to our knowledge an expansion for each term in (5.13) cannot be found in the literature. Thus, we will seek to find it using the Lorentzian inversion formula.
Inversion I: higher-spin currents. We now use the inversion formula (3.16) to recover the thermal one-point functions of the currents J , and give implicit results for the higher families [φ i φ i ] n, with n > 0. First one has to understand the discontinuities along the axis Im z = 0 with Re z > 1 for each term in (5. We now apply the inversion formula (3.16) to get the contribution of each term in (5.13). We focus on the integral, multiplying overall factors at the end. We also denote the spin as , rather than J. For terms with m ≥ 1 we find, following the same approximation scheme as in section 4 for d > 2 MFT (see around (4.16)), To derive this, we use the Poisson resummation formula to turn a sum over Matsubara frequencies ωn = 2πn into a sum over shifts in τ : Here, f is the Fourier transform of f . Thus, we can Fourier transform Gij(ω, k) (treating ω, k as continuous), which gives the Yukawa potential e −m th |x| |x| . Then we sum over integer shifts in τ to obtain (5.13).

JHEP10(2018)070
Expanding the integrand in (5.15) at smallz, Thez integral at leading order gives rise to the contribution of the first double twist family with h = 1/2, via with a pole at h = 1/2. Plugging this into (5.15), we now perform the z-integral to extract the residue at h = 1/2, which is found to be where ∆ = 1 + and K ∆− 1 2 is the modified Bessel function. The full result requires a sum over m as in (5.13); performing this sum, and appending overall factors from the inversion formula, we find that the thermal coefficient a J for the higher-spin currents J is This sum can be performed to yield the following result: We can translate this to a result for the thermal one-point function, b J , itself, using known results in the literature for the OPE coefficients f φ i φ i J , together with our φ i normalization in footnote 30. From e.g. [78], in d = 3 we have (5.20) Using the relation (2.12) between a O and b O , we find which is the ratio which is independent of the norm of J . This is an elegant result. The case = 2 corresponds to the stress tensor. In this case, the sum may be further simplified to yield

JHEP10(2018)070
As mentioned in the introduction, this result has implications for higher-spin black hole solutions of Vasiliev higher-spin gravity in AdS 4 . The translation invariance of thermal onepoint functions means that (5.21) are proportional to the higher-spin charges of the CFT at finite temperature. Together with the thermal mass m th ∼ σ β , these charges fully determine the "higher-spin hair" of the putative black hole solution dual to the CFT thermal state with vanishing higher-spin chemical potentials. This black hole has not yet been constructed, due to difficulties in interpreting and solving Vasiliev's equations. Our result provides a benchmark, both for any explicit candidate black hole, and for a physical interpretation of proposed constructions of higher-spin gauge-invariant charges (see e.g. [79][80][81][82][83]).
It is not much more difficult to derive the one-point functions of the n > 0 families appearing in (5.7). One simply has to keep higher orders inz in (5.16): a term of O(z n ) gives a pole at h = 1 2 + n, i.e. for spin-operators with ∆ = 1 + + 2n.
Inversion II: scalars. The above results are incomplete in the scalar sector, and present a small puzzle. Note that our final expression (5.18) was actually valid all the way down to = 0, which would correspond to a scalar with dimension ∆ = 1, even though such an operator is absent. The same would happen for the poles with higher n > 0, which would seem to indicate the presence of spurious scalars with odd integer dimension ∆ ∈ 2 Z + − 1. Moreover, we did not recover the ∆ = 2m scalar poles corresponding to σ m exchanges, nor the unit operator. As we now show, these issues are remedied by considering the arc contributions to the inversion formula. Following the notation in section 3.2.2 where z = rw andz = rw −1 , we are interested in computing the w → e iφ ∞ behavior for each term in the propagator (5.13). In this limit, the only surviving term is given by the m = 0 term, The contribution of the integral correction to the inversion formula is given by (3.12). When plugging in the asymptotic value of the propagator this becomes a arcs (∆, ) = 2K (1 + (−1) ) The integral over w in the limit in which |w| → ∞ is trivial and simply gives a factor of 2π when = 0 and 0 when > 0. This indeed confirms that the thermal coefficients quoted above for the currents with > 0 are correct. For = 0, we are left with an integral over r:

JHEP10(2018)070
Since Γ(x) has poles at each negative integer value of x, we can express the function a(∆, 0) around each m ∈ Z ≥0 as These poles do two things. First, they cancel all spurious scalar poles of a(∆, 0) at ∆ ∈ 2Z + − 1. Second, they give the correct poles for the actual scalar operators of the theory, which have ∆ ∈ 2Z + , as well as the unit operator pole. Let us first analyze the case m = 0. This simply returns a thermal one-point function of 1, corresponding to a correctly normalized unit operator.
Next we take m = 1. This pole at ∆ = 1 and = 0 should cancel the spurious scalar pole of the previous analysis. From (5.19), we get

Holographic CFTs
We now make some comments on large N CFTs in the context of AdS/CFT. A universal set of contributions to the OPE expansion (5.2) comes from the stress tensor, T µν , and its multi-traces, [T . . . T ], which necessarily appear in the φ × φ OPE for any φ. In a CFT with a weakly coupled gravity dual, these terms represent the purely gravitational interactions between the bulk field Φ, dual to φ, and the thermal geometry. The form of these contributions is sensitive to the gap scale to single-trace higher-spin operators (J > 2), ∆ gap . We would like to understand how.

JHEP10(2018)070
is dual to an AdS d+1 -Schwarzchild black brane geometry with inverse Hawking temperature β, 32 and the stress tensor contributions in the OPE decomposition (5.2) are dual to the exchange of arbitrary numbers of gravitons between Φ and the black brane. In a heavy probe limit 1 ∆ φ M pl L AdS , the connected two-point function may be computed as the exponential of a geodesic length, φ(x)φ(0) β ∼ e −∆ φ x . The disconnected component of the correlator, ∼ φ 2 β , is computed as an infall of each particle into the black brane horizon. 33 This disconnected contribution goes to a constant plus e −m th x corrections, and thus becomes more important at sufficiently long distances.
It is instructive to examine the classic case of strongly coupled N = 4 super-Yang-Mills (SYM), with SU(4) R symmetry. The single-trace scalar spectrum consists of the Lagrangian operator, as well as the 1/2-BPS operators O p with p = 2, 3, 4, . . ., which live in the [0, p, 0] representation of SU(4) R and have conformal dimensions ∆ = p. The R-symmetry constrains the thermal two-point functions O p O p β to take the form

+ (stress tensor terms). (5.31)
That is, in the OPE decomposition of O p O p β for any p, the stress tensor terms are the only terms besides the MFT contributions at leading order in 1/N . This follows from the absence of R-singlets in the single-trace spectrum besides the identity operator, the Lagrangian operator and the stress tensor, and the fact that the Lagrangian carries charge under a emergent U(1) Y bonus symmetry [90]. The stress tensor contribution, a T , exhibits a famous dependence on the 't Hooft coupling λ [91], In relating a T to b T , we have used that f OpOpTµν and C T are λ-independent.
In more general theories with large ∆ gap and a sparse spectrum of light operators, there can be single-trace global symmetry singlets, so contributions from operators other than the stress tensor to the scalar two-point function g(τ, x) are possible. As ∆ gap decreases, there are different possible sources of ∆ gap corrections to thermal correlation functions. First, the low-spin OPE data receive power-law corrections in ∆ gap . This includes the OPE coefficients of double-trace operators (see [92] for an N = 4 example). In addition, there are e −∆gap corrections due to new contributions of massive string states with ∆ ∼ ∆ gap . At finite ∆ gap , there are many possible behaviors.
Finally, note that if instead we examine thermal two-point functions of the stress tensor, T µν T ρσ β , the effects of large ∆ gap are more visible. For instance, T µν T ρσ T λη and T µν T ρσ O couplings scale with inverse powers of ∆ gap [84][85][86][87], thus suppressing various possible contributions to T µν T ρσ β in the OPE limit. For many reasons, it would be interesting to extend the methods discussed herein to the case of spinning external operators, 32 At infinite spatial volume, thermal AdS is thermodynamically disfavored. 33 This interaction requires a nonzero cubic coupling between two gravitons and Φ; this is forbidden at the two-derivative level, but may appear at the four-derivative level in the form of a φC 2 µνρσ coupling, where Cµνρσ is the Weyl tensor (see e.g. [88] for an application). Such couplings are, however, suppressed by the mass scale of higher-spin particles in the bulk [87] and, in more general theories of gravity, by universal bounds [87,89]. and to T µν in particular; this would allow us to study the purely gravitational physics of the thermal geometry in AdS, without the need for a probe scalar field.

Large-spin perturbation theory
So far, our discussion of thermal two-and one-point functions has been in theories where we have enough analytic control to explicitly compute the thermal two-point functions, which we can invert to obtain one-point functions. How can we analyze theories for which we don't have any direct method of computing two-or one-point functions, such as the 3d Ising CFT? Inspired by studies of CFT four-point functions in Minkowski space, we use the inversion formula to set up a bootstrap algorithm for the thermal data in any CFT. The inversion formula takes in the two-point function, and returns its decomposition in the s-channel OPE. Crucially, any presentation of the two-point function could be inserted into the inversion formula, and its inversion would yield how the given presentation is related to the s-channel data. Here, we will invert the t-channel OPE to the s-channel data. This relates the one-point functions of operators in the theory in a highly non-trivial fashion. By iterating these relations, we solve the thermal bootstrap in an all-orders asymptotic expansion in large spin, J.
We can make explicit use of our tools in the 3d Ising model. As in MFT and the O(N ) model, low-twist operators can be arranged into double-twist families with relatively small anomalous dimensions. By combining our analytic tools with previous results from the four-point function bootstrap, we will estimate the thermal one-point functions of the operators in the lowest-twist family, [σσ] 0 , as a function of b and b T .

Leading double-twist thermal coefficients
Let's study the two-point function of identical scalars g(z,z) = φ(z,z)φ(0, 0) β , (6.1) in some CFT at finite temperature. We want to understand the "contributions" to the thermal coefficients of operators in the [φφ] 0 family, a [φφ] 0,J , from other thermal data of the theory. Our starting point is the t-channel OPE (z ∼z ∼ 1), which we will systematically invert to the s-channel data a(∆, J). Expanding the Gegenbauer polynomials yields a power series in 1 − z and 1 − z: For future convenience, let's define the coefficients Let's suppose we are considering J large enough so that the contributions of the arcs in (3.16) vanish. Before inverting the t-channel OPE into a(∆, J), let's analyze a few key features of the inversion formula: • First, as discussed previously, recall that poles of a(∆, J), associated with physical operators, come from the region nearz = 0. A term likez a in the expansion around z = 0 inverts to terms of the form and gives poles at h = ∆ φ + a + m. Such poles represent infinite families of operators with unbounded spin J and scaling dimensions ∆ = 2∆ φ + 2a + 2m + J. Of course, in interacting CFTs, operators should have anomalous dimensions. We discuss the effects that shift the locations of these naïve poles to their correct values in section 6.2.3.
• Next, let's imagine a term of the formz a (1 − z) c in g(z,z) expanded around the double-lightcone limit (z = 0 and z = 1), and invert it. The z integral determines the residue of the poles in (6.7) as a function ofh (recall thath = h + J). Typical z integrals are of the form where we've set the upper limit on the z integral equal to ∞ becausez is infinitesimal. The discontinuity is Disc[(1 − z) c ] = 2 sin(−πc)(z − 1) c , so the integral gives 2 sin(−πc) Note that this is naturally a term in a large-h expansion, since (6.10) Thus, we see that terms in the double-lightcone expansion of g(z,z) correspond to power law corrections in 1/h, or equivalently in 1/J, to the thermal coefficients of families of operators (in the s-channel).

JHEP10(2018)070
• As highlighted in section 3.3, the t-channel OPE in (6.3) is valid for the range 0 ≤ z,z ≤ 2. Thus, we are only justified in integrating the t-channel OPE between 1 ≤ z ≤ 2 for the z integral. In general, we don't have an expression for g(z,z) that is valid in the region 0 ≤z ≤ 1 and z > 2. Luckily, for the integrands of interest, the z integral in the range z > 2 is exponentially suppressed in largeh, schematically as Therefore, we can work with the t-channel OPE, integrate it in the region of its validity (from z = 1 to z max = 2), and obtain an all-orders expansion in 1/h, with undetermined exponentially-suppressed corrections.
Now that we have oriented ourselves, let's calculate the contributions to a(∆, J) from a single primary operator O ∈ φ × φ in the t-channel OPE in its full glory. We take the terms in the t-channel OPE corresponding to O, and invert them to the s-channel via the inversion formula. We use a (O) (∆, J) to denote the contribution to a(∆, J) from the inversion of the contribution of O in the t-channel. Then we find The reasoning behind separating out the sums to infinite spin will become apparent in section 6.2. We have defined the function (6.14) Here B 1/zmax (h−∆−c, 1+c) is the incomplete beta function, which decays as z −h max at largeh. The factors Γ(∆ φ + m − h) in the numerator in (6.12) give poles at h = ∆ φ + n for n ∈ Z ≥0 . Naïvely, these are the poles corresponding to the [φφ] n families, at the naïve dimensions ∆ = 2∆ φ + 2n + J, without anomalous dimensions. However, the correct -36 -JHEP10(2018)070 a(∆, J) has poles in ∆ at the exact dimensions, including anomalous dimensions, so our naïve poles are shifted to their correct values, This has a subtle but important effect on the thermal coefficients. When one takes the residue, there is an extra factor dh/dJ that depends on the derivative of δ n (h), since as opposed to Note that dh/dJ = 1 when anomalous dimensions vanish. In section 6.2.3, we will provide a consistency check that the poles are indeed shifted to their correct locations. In our discussion of sums over families below, we include dh dJ for two reasons: firstly, it greatly simplifies the analysis of the asymptotics of such sums; secondly, we have in mind a situation where the anomalous dimensions δ n (h) are known through other means (e.g. the vacuum four-point function bootstrap), and we would like to use that information in the thermal bootstrap.
Finally, evaluating the residues of a ( Note thath is implicitly defined as a function of J byh = ∆ φ + n + δ(h) + J. As we have emphasized above, properties of using the OPE with the inversion formula, these contributions are naturally organized as power-law corrections in largeh. The function S c,∆ (h) behaves as at largeh, and c m (J) behaves as a constant to leading order in 1/h. So a given term in (6.18) starts at orderh −(h O −∆ φ +k+1) . Thus, we see that the contribution of an operator O in the t-channel behaves at a rate controlled by its twist. Concretely, the leading contributions in 1/h are given by the k = 0 term of the sum in (6.18), For the leading double-twist family [φφ] 0 , this further simplifies to In writing the last line, we have assumed that δ [φφ] 0 (h) grows slower thanh ash → ∞.
To help understand the examples that follow, let us introduce a diagrammatic language that helps keep track of terms in large-spin perturbation theory of thermal data. Our diagrams can be thought of as analogs of the four-point function large-spin diagrams for the thermal case. We do not have a rigorous definition of these diagrams or a complete set of rules for using them. Nevertheless, they will help organize the discussion. 34 For example, we can understand the fact that O ∈ φ × φ in the t-channel OPE inverts to give contributions to a [φφ]n proportional to a O via the diagrams in figure 5. Let's start with the t-channel diagram in figure 5a. We should read this diagram from left to right as two φ operators approach each other on one side of the thermal circle (corresponding to the t-channel), and fuse into O, which in turn gets an expectation value. The diagram illustrates that this process should be proportional to the three-point coefficient f φφO and to the one-point function b O , which is indeed the case by the definition of a O . Now, let's relate this process to the s-channel. The diagrammatic rule relating the s-and t-channels is given by taking the two external operators to the other side of the thermal circle around opposite sides. This converts the process in figure 5a to the process in figure 5b. Reading the resulting process from right to left, we interpret it as two external φ's fusing into operators in the [φφ] n families, which get expectation values proportional to a O .
Summarizing, we reiterate that the thermal coefficients of families of operators are organized into large-spin expansions, with the operators O in the OPE contributing perturbatively at order determined by their twist. Since the unit operator has the lowest twist 34 Large-spin diagrams for four-point functions can be understood as physical processes in the massive 2d effective theory defined in [52]. It would be nice to develop a similar understanding of the diagrams here. For now, our diagrams are simply mnemonic devices. in any unitary theory, it gives the leading contribution for large-spin members of doubletwist families. A second important contribution comes from the stress tensor O = T µν , which gives a universal contribution proportional to the free energy density. These two universal contributions were written in (1.5) in the introduction. Furthermore, we also see that the contribution of a given O is linear in its one-point function. This greatly simplifies perturbatively solving for the one-point functions, especially when one considers the corrections from sums of families, which we shall explore below.

Contributions of double-twist families: resumming, asymptotics, and other families
We have seen that inverting any single operator O in the φ × φ OPE gives contributions to the one-point functions of the [φφ] n families, and only these families. But many other operators appear in the φ × φ OPE, and the inversion formula must pick up their existence. How can we extract their one-point functions from φφ β ?
As seen in section 6.1, any finite number of terms in the t-channel OPE have "regular" z expansions aroundz ∼ 0, with integerz powers obtained from the Taylor series of some collection of terms of the form (1 −z) c . Inverting such terms will only give poles at h = ∆ φ + n, which thus correspond to the [φφ] n families. So, to obtain the necessary poles at other locations, we need to find terms in the t-channel expansion that behave asz c with c / ∈ Z ≥0 nearz ∼ 0. Such terms invert to poles at h = ∆ φ + c + n, and correspond to different families determined by c. We will call such terms "singular", in analogy with the Casimir-singular definition of [30,52]. Singular terms are characterized by having discontinuity aroundz ∼ 0, which means they would be picked up by the Lorentzian inversion formula that takes the two-point function and inverts it to the t-channel. The only way we can obtain singular terms is from tails of the infinite sums of families in the t-channel OPE, which, when summed up, will have differentz behavior compared to any finite number of terms. A related motivation for understanding this problem is to compute the contributions of the [φφ] n families to their own one-point functions.
We will now explain how to systematically compute the contributions of double-twist families to the one-point functions of operators that appear in the φ × φ OPE. To begin, let's focus on the contribution from summing the tail of the particular double-twist family [φφ] 0 . Once again, we start with the t-channel OPE expansion, which we now try to understand in the double-lightcone limit (z,z) ∼ (1, 0).
Let {O} be a set of operators in the φ × φ OPE with low twist. Inverting their terms in the t-channel via section 6.1, we obtain from (6.21) the leading 1/h behavior of the 1-point functions of the [φφ] 0 family, We wish to apply the inversion formula to this sum. If we naïvely invert each term in the sum, and then sum over the family, we notice that the sum over the contributions of each individual member of the family diverges. In other words, the inversion formula and the infinite sum over the family do not commute. So, we have to sum over the family first before inverting to the s-channel. This is in line with our anticipation that the poles for other families must arise from the tails of the sums over infinite families, such as [φφ] 0 . If the t-channel sum and the inversion integral commuted, we would only ever get poles for the [φφ] n families from the inversion formula.

Analytic and numerical formulae for sums over families
Let's analyze the t-channel sum over a double-twist family more carefully, in the spirit of [30]. Consider the sum over a particular term in the S c,∆ (h) expansion of one-point functions of an arbitrary double-twist family, 35 Here,h = h f + + δ(h) runs over the family with anomalous dimensions δ(h), and h(h) = h f + δ(h) are the half-twists of the operators in the family. We have switched to denoting spin by for sums over families, to avoid conflict with applying the inversion formula to these sums later on. We have left out the (−1) factor for now; we will return to it later. In general, this is a difficult sum to evaluate, and we don't yet know of an exact treatment. However, since anomalous dimensions δ(h) for the families of interest are small for largeh, we can work order by order in δ(h). Concretely, we can split the sum overh as for some large enoughh 0 = h f + 0 such that the anomalous dimensions δ(h) are sufficiently small forh >h 0 , and expand the infinite sum piece in small δ(h) log(1 − z), (6.26) 35 For he and h f , the indices e and f stand for "external" and "family". For the two-point function of two identical scalars φ, he = ∆ φ .

JHEP10(2018)070
This expansion is valid in a regime e −1/|δ(h)| < |1 − z| < e 1/|δ(h)| , which is near the doublelightcone limit for small anomalous dimensions. Now, the dependence on z -which controls the discontinuity in the inversion formula -can be factored out of theh sum.
Recall that the anomalous dimensions δ(h) are themselves analytic functions ofh, and they can be computed perturbatively in a large-h expansion [30][31][32][33][34][35][36]. For example, the anomalous dimensions of a family [φφ] 0 can be expanded at largeh, and include terms such as How can we compute the sums overh? With a more general treatment in mind, let's consider the sum with a summand f (h) which grows at most as a power law at largeh. The summands of interest for us are of the form 36 If we tried to expand in smallz and compute theh sum order by order inz, we see that for large enough powers ofz we get divergent sums inh. In fact, this is to be expected. By the existence of the inversion formula, we know that such sums must have asymptotic pieces that reproduce the singular termsz a , which could not possibly be obtained from expanding inz first (which produces only integer powers ofz). Thus, we expect the result of such a sum to be of the form with A ⊂ R\Z ≥0 some set of numbers which are not non-negative integers that we have to determine. (We can also havez a log mz terms that we will write as ∂ m ∂a mz a .) Our task is reduced to computing the coefficients c a and α k . To do this, we will separate the sum into asymptotic parts, which reproduce thez a terms, and leftover regular parts that are convergent sums inh, with which we can compute the α k coefficients. First, we determine the large-h asymptotics of f (h) in terms of the known functions Note that the set A is determined by the asymptotics of f (h), but the expansion can be written for any choice of ∆. For summands of interest like in (6.29), the asymptotic expansion is determined algorithmically from the large-h expansions of S c,∆ (h) and of δ(h) a la (6.27). Once the asymptotics in (6.31) are obtained, we can compute the singular 36 We will also be interested in sums including derivatives ∂ m c Sc,∆(h), as we will discuss later.

JHEP10(2018)070
terms from the tails of the sums over the asymptotics, by using the crucial identity of the integer-spaced sum Note that the first term is singular and the second term, proportional to 2 F 1 (· · · ,z) is regular.
We claim that the noninteger-spaced sum over S c,∆ (h) (withh determined by the anomalous dimensions δ(h)) has the same singular piece as the integer-spaced sum, This can be verified by an argument due to [30] applied to the present case. We convert the sum to a contour integral via Cauchy's residue theorem, 37 where γ is a contour along the real axis that picks up the desired poles. We can deform the contour to one that runs parallel to the imaginary axis, plus arcs at infinity. The singular terms come from the asymptotics of this integral. As long as δ(h) grows slower thanh as h → ±i∞, the asymptotic region of the integral approaches a δ(h)-independent constant exponentially quickly, since π cot(π(h −h 0 − δ(h))) → ∓1 + O(e ∓2s ) ash → ±is. (6.35) Therefore, the singular pieces are independent of δ(h). To be even more concrete, we can subtract the contour integral versions of the noninteger-and integer-spaced sums, and notice that the asymptotics vanish, or that if we expand the difference in smallz, the integrals inh are convergent term by term so the singular terms must have canceled. Thus, we see that Note that we can always choose ∆ = h e in the asymptotic expansion in (6.31) to simplify the organization of the singular terms, 37 This is also known as the Sommerfeld-Watson transform.

JHEP10(2018)070
We are left with computing the coefficients α k of the regular terms. Since we have extracted the asymptotics, we can write convergent expressions for α k by subtracting the asymptotics and expanding in smallz. This will be best done by once again converting the sum overh into a contour integral in the complexh plane, via Cauchy's theorem. We want to write a contour integral of the form where h c =h 0 + δ(h 0 ) − , for some small > 0. This contour integral will equal the sum only if f (h) decays fast enough on the arcs at infinity so that we may drop them, and if f (h) does not have any simple poles for Reh ≥ h c . The second condition is easily remedied in case f (h) does have poles, by simply removing the residues coming from those poles. The first condition is related to the more important issue that for large enough k, theh growth is divergent, and the sum over k and the contour integral do not commute. However, we can regulate the integral by subtracting the divergent asymptotics in the form of the integer-spaced sum until we get a convergent integral (for which the arcs vanish as well), and add back the known result of the integer-spaced sum. This gives the following formula Here, K should be at least k, but larger K gives a faster converging integral. In the last line, we have added back terms with r k , which is the coefficient ofz k for the integer spaced sum in (6.32), a + ∆ − n − h e k (−1) k π cot(π(a + ∆ − n −h 0 )) (−1) n n!Γ(−a)Γ(a − n + 1) .
The contour integral can be integrated numerically in Mathematica to high precision. Finally, to finish our discussion, let's consider the alternating sum, This sum is convergent order by order in thez expansion, so one does not need to subtract off asymptotics. This sum is given by where the coefficients α − k are given by the contour integral with the replacement cot → csc, Collecting our calculations, the full sum over operators with even spin is given by

Corrections to one-point functions from double-twist families
Armed with the technology to compute the sums over double-twist families, we return to understanding their contributions to one-point functions of operators. Let's recall the t-channel sum over [φφ] 0 in (6.23), and expand it in δ log(1 − z) as in (6.26),

JHEP10(2018)070
For [φφ] 0 , h f = ∆ φ , but we've kept it general here to demonstrate the general structure. Now, let's invert this piece of the two-point function to the s-channel. The integer powers z k invert to poles for [φφ] n , giving contributions to the one-point functions of these families, including [φφ] 0 itself. The contributions are controlled by which inverts to a term S in the large-h expansion of thermal coefficients. For example, including the self-corrections of [φφ] 0 to leading order in large spin yields We should remember that if we started the sum at some highh 0 , we should individually add the contributions of the low-lying members of the family that were excluded from the sum. Once we have computed the self-corrections as in (6.51), there is nothing that stops us from iterating this procedure and computing the self-corrections from the new, once-selfcorrected one-point functions. Instead of iterating indefinitely, we can solve for the fixed point of these self-corrections with a little cleverness. We have provided a method for this in appendix E.
The singular termsz a in (6.48) give poles at h = ∆ φ + a + n. We expect that these poles correspond to other families of operators with the given naïve twist -we will soon explore which families. Let's denote these families by [∆ + a] n for now. We see that the [φφ] 0 family contributes to the one-point functions of the [∆ + a] n families. The sum is over m such that S a,∆ φ appears in the asymptotic expansion of δ m m! S h O −∆ φ ,∆ φ . Of course, this is rather schematic, since for interacting CFTs, the spectrum of families of higher-twist operators is very complicated, with large anomalous dimensions and mixing among families. Regardless, these contributions are present asymptotically in large J.
In general, can we say which other families of operators appear in the asymptotics of the sum of a given family, and therefore receive contributions via (6.48)? The largespin expansion of the anomalous dimensions and OPE coefficients allows us to answer this question. Suppose O is an operator in the φ × φ OPE. Then, O corrects the anomalous dimensions of the [φφ] 0 family, via the large-spin diagram in figure 6 [30,52]. Consequently, there is a term in the asymptotic expansion of δ [φφ] 0 that goes like where δ [φφ] 0 is some coefficient. Now, imagine the contribution of the identity operator to the [φφ] 0 thermal coefficients, which goes like S −∆ φ ,∆ φ (h) to leading order. Therefore, the sum over the [φφ] 0 family to first order in δ [φφ] 0 contains the asymptotic term  figure 7 to the s-channel. As demonstrated in figure 8, the resulting s-channel process is proportional to b [OO]n , so the inversion to the s-channel must have produced poles for [OO] n . We therefore see that the intuition from the diagrams agree with concrete calculations! One can also check that the expression for b [OO] 0 obtained from the φφ β correlator agrees with the expression obtained from the OO β correlator to leading order in the large-h expansion.
In fact, the situation is much more general. For example, we can consider other terms in the asymptotic expansion of a [φφ] 0 , such as S h O −∆ φ ,∆ φ (h) coming from some other operator   O . Then, the sum over [φφ] 0 to first order in δ [φφ] 0 produces the singular termz 2h O +h O −∆ φ , naïvely corresponding to multi-twist families [OOO ]. The diagrams for the sum over this asymptotic piece and the corresponding s-channel process are given in figure 9. We could also work to higher order in the anomalous dimensions, and obtain poles for multi-twist families, and so on.
However, we are unsure what the precise rules are for which diagrams are allowed, and how to interpret them in general. We will leave deriving these diagrams from physical arguments and further generalizing them to a future project.

Corrections to pole locations
We are also in a position to address the issue of naïve versus true locations of poles of a(∆, J), raised in section 6.1. Corrections to the locations of poles essentially arise from the asymptotics of the sums over the terms S [OO] 0 (J) given in (6.55). The t-channel sum over the [OO] 0 family looks like to leading order in (1 − z). Focusing on the term a Let's assume that O = φ, so ∆ O = ∆ φ . Then, expanding to leading (constant) order in δ [OO] 0 (h), the sum becomes If O = φ is in the φ × φ OPE, we should consider the sum over δ [φφ] 0 (h)a [φφ] 0 (J), Note that we have we have used the asymptotic expansion Doing the z integral and summing over O, we obtain the contribution tog 0 (z,h) Now, let's combine this with the contribution of the unit operator tog 0 (z,h), In the last step, we reverted the asymptotic expansion in (6.64) to separate the contribution to the pole location from the residue. This looks like the first few terms in the expansion of where we have used the asymptotic expansion 38 which can be obtained from the asymptotic expansion in (6.64). Such a term inverts to poles for the multi-twist families [φφO O ] n , as we expected from the diagram. The lowest-twist family has the leading residue

JHEP10(2018)070
Now, let's reverse the diagram, which tells us to sum over [φφO O ] 0 in the t-channel to constant order in their anomalous dimensions. Performing this sum yields the singular terms which contribute the expected first-order shift to the [OO] 0 poles, sincẽ We see that we begin recovering the correct pole locations of the [OO] 0 family in the inversion of the φφ β correlator. Once again, the higher-order corrections come from diagrams with exponentiated anomalous dimensions analogous to figure 10.

Case study: [σσ] 0 β in the 3d Ising model
Our primary example for applying the above technology is the 3d Ising CFT. At this point, much is known both analytically and numerically about the spectrum and OPE data of the 3d Ising CFT. This abundance of data makes the 3d Ising CFT a natural and ideal candidate for studying thermal correlators. In [30], the low-twist spectrum of the 3d Ising CFT has been computed via the lightcone four-point function bootstrap. Especially relevant to our analysis here is the analytic computation for the anomalous dimensions and OPE coefficients of the most important double-twist family, [σσ] 0 , which has the lowest twist trajectory. Taking the spectrum and OPE data as input, we will apply the thermal bootstrap to study the thermal coefficients of the [σσ] 0 family. The most natural way to get a handle on the [σσ] 0 family is by studying the thermal correlator σσ β . Let's remind ourselves about the relevant low-twist spectrum of the 3d Ising CFT. The first few lowest-twist primary operators in the σ × σ OPE are σ × σ = 1 + + T + =4,6,...
[σσ] 0, + . . . . (6.75) Our strategy will be to determine the thermal coefficients of [σσ] 0 in terms of b and b T , which we treat as unknowns. While we do not determine the values of b and b T here, our work paves the way towards it. In a future paper [93], we will show how information about the low-twist families can be used in conjunction with the KMS condition in the Euclidean regime to "tie the knot" on the thermal bootstrap and estimate some thermal coefficients in the theory.
To numerically study the thermal coefficients in the [σσ] 0 family, we use the scaling dimensions of σ and , obtained from the numerical bootstrap study [29] ∆ σ = .5181489(10), ∆ = 1.412625 (10) . Using our result (6.21), together with these numerical values, we compute the leading contributions to the [σσ] 0 one-point functions, To emphasize the utility of this result we can write the large-spin expansion of the thermal coefficients where terms on each line come from the unit operator, the stress tensor and the operator respectively. The final ". . . " include contributions of other operators that are either suppressed in the 1/J expansion or with coefficients small enough that they can be neglected for reasonable values of the spin.
To go beyond asymptotically large spin and estimate thermal coefficients for operators with small spin, we should include higher-order corrections in 1/J. The next contributions come from the [σσ] 0 family themselves. Thus, we need to sum over the [σσ] 0 family next. We use the leading expressions in the large-spin expansion (6.27) of the anomalous dimensions of [σσ] 0 , which were computed in [30] as Upon first iteration, when considering the corrections from [σσ] 0 to itself only once, the corrected thermal coefficient is given by (6.51), We can compute the fixed point of the self corrections above using appendix E, with (6.77) as input. It turns out that the self-corrections of operators in the [σσ] 0 family is given by convergent sums over operators in the [σσ] 0 family, so one can also evaluate the sums numerically by choosing a large spin cut-off. By recursively repeating this numerical process the results converge to the fixed point determined analyticallyà la appendix E. 39 To be concrete, the table below shows a few examples for the values of the thermal coefficients a [σσ] 0, and for the thermal one-point functions b [σσ] 0, : 39 We find that for small values of the spin the contribution of the stress-energy tensor is the most affected by self-corrections, with a 20% correction for J = 4.

Conclusions and future work
Modern advances in the conformal bootstrap have focused almost entirely on constraining OPE data using CFT correlation functions in flat space. Is there potential for more? A broader perspective on the bootstrap suggests future extensions toward probing dynamical questions in CFT, which are not obviously determined by OPE data in a tractable way. As a step toward this end, we have developed an approach to bounding CFT observables at finite temperature. Treating the thermal two-point function on S 1 β × R d−1 in analogy with the flat-space four-point function, and the KMS condition as the analog of the crossing equations, one extracts constraints on the thermal one-point functions of local operators. A key intermediate tool (of independent interest) in realizing this approach was to derive a Lorentzian inversion formula (3.16) which, given a thermal two-point function, extracts thermal one-point coefficients and operator scaling dimensions. We applied this technique to the d = 3 critical O(N ) model, which yielded thermal one-point functions of higher-spin currents (5.21) and some scalar operators (5.30). More generally, we developed a large-spin perturbation theory, applicable to any CFT, in which thermal one-point functions are determined via an analytic expansion in inverse operator spin J. This included the universal contributions to thermal coefficients of double-twist operators, a [φφ] 0,J , from the presence of the unit operator and the stress tensor in the φ × φ OPE (1.5). By summing over entire families of operators and plugging back into the large spin expansion, one can solve for CFT data to increasingly high accuracy. Together with the KMS crossing condition, this suggests an iterative algorithm, discussed further below, with which to "solve" the thermal sector of an abstract CFT.
There are many future directions to explore: • In this work, we mostly consider a single thermal two-point function. However, the same one-point coefficients appear in the OPE decomposition of every two-point function in a theory (except when forbidden by symmetry). Thus, it might be very constraining to study larger systems of two-point functions simultaneously.

JHEP10(2018)070
• A more straightforward generalization of our work would be to study thermal twopoint functions of spinning operators. This is likely easier than studying spinning four-point functions on R d , due to the simplicity of the spinning thermal conformal blocks [54].
• Our Lorentzian inversion formula makes it straightforward to compute the perturbative expansion of thermal data to all orders in 1/J, using the t-channel OPE for z < 2. However, there are also nonperturbative corrections that decay exponentially in J, coming from the region z > 2 (outside the regime of validity of the t-channel OPE). How can we compute these corrections? Answering this question may require understanding the full analytic structure of thermal two-point functions better.
• It would be interesting to study more general compactifications. For example, one could study two-point functions on T n × R d−n for n ≥ 2. On the other hand, there can also be multiple one-point structures on T n for n ≥ 2, so there is more data to compute. For recent work on CFTs on spatial tori, see [94,95].
• We derived thermal one-point functions of all single-trace operators in the critical O(N ) model in d = 3. A clear target for the future is to generalize these results to other slightly broken higher-spin CFTs, such as the Chern-Simons-fundamental matter theories that are continuously connected to the O(N ) model [96,97]. The thermal mass and some current-current correlation functions at nonzero temperature have been computed in the large N limit of these theories, for arbitrary 't Hooft coupling λ [96,[98][99][100]. It would be satisfying if the thermal one-point functions b J in these Chern-Simons-matter theories take the same form as in (5.21), with the appropriate thermal mass m th (λ). More generally, we would like to understand the constraints of slightly broken higher-spin symmetry on thermal correlations, in the spirit of [34,47,101,102].
• Through the study of holographic CFTs one can get a better intuition for the applicability of the inversion formula down to small values of the spin. Such a direction would entail studying the holographic thermal two-point function in the regimes discussed in section 3.5, in which |w| → ∞. Besides offering better intuition for the applicability of the inversion formula, as discussed in section 5.2, the study of such a regime would also be illuminating for understanding the thermal properties of the stress-energy tensor as implied by black hole physics. It should also be possible to define geodesic Witten diagrams [103,104] for black hole backgrounds in AdS d≥4 , which should define an effective two-point "thermal conformal block" for d ≥ 3 CFTs with large higher-spin gap.
• In section 6, it proved useful to use diagrams to organize terms in large-spin perturbation theory for thermal correlators. It would be nice to place these diagrams on firmer footing by giving a complete specification of the rules they satisfy and what terms they correspond to. This problem is already interesting in the context of largespin perturbation theory for four-point functions [30,105], where the diagrams have an interpretation in terms of physical processes in a special conformal frame [52].

JHEP10(2018)070
• We have made predictions for thermal one-point functions in the 3d Ising CFT in terms of some unknowns, of which we expect b T (computed via Monte Carlo in [58][59][60]) and b are the most important. It would be nice to compute b . To our knowledge, it is not present in the literature, but should be straightforward to compute using e.g. using Monte Carlo simulation [106]. 40 • While in this paper we have made contact with the 3d Ising model by finding the large-spin expansion for the thermal one-point functions b [σσ] 0,J , one can imagine a more involved iterative strategy to solve the thermal bootstrap in the double-lightcone limit. This strategy can be summarized in the following diagram: Following our study in section 6.3, we start by considering the OPE presentation of the thermal two-point function σσ β , with the thermal one-point functions of a few low-twist operators as unknowns (in section 6.3, we consider b and b T as unknowns). Then, we use the inversion formula on σσ β in the double lightcone limit to determine the thermal coefficients of all remaining operators in the [σσ] 0 family as functions of the unknowns. Next, using the technology we developed in section 6.2, we sum over the [σσ] 0 family to determine the self-corrections to the thermal coefficients of the [σσ] 0 family, and also determine the thermal coefficients for the [σσ] 1 and [ ] 0 families.
In principle, this process can be iterated further by summing over more and more families, and obtaining higher terms in the large-spin expansion. Also, by studying the thermal two-point functions of other operators, we get alternative handles on the thermal coefficients of families of operators. For instance, studying β yields more direct information about the [ ] 0 family. Once the thermal coefficients of families of interest are determined to desired order, we have expressions for a large part of the low-twist spectrum, which still depend only on the unknown thermal one-point functions of the chosen low-twist operators (b and b T ). Finally, perhaps these unknowns can be determined by moving away from the double lightcone limit and applying the KMS condition, thus determining the low-twist thermal one-point functions of the 3d Ising CFT. 40 See [3] where similar quantities were computed for the O(2) model. • The eigenstate thermalization hypothesis (ETH) suggests that we can study thermal correlators as a limit of expectation values in a single eigenstate |O with sufficiently large dimension. See [107] for a recent discussion of ETH in the context of CFTs. Assuming ETH, a thermal two-point function φ(x 1 )φ(x 2 ) β is a limit of a family of four-point functions 41 It would be interesting to understand whether the ability to view thermal correlators as limits of pure correlators can bring new constraints to the thermal bootstrap. Note that certain properties of vacuum four-point functions may not survive the thermodynamic limit. For example, the analyticity structure changes, with the development of new "forbidden singularities" reflecting periodicity of the thermal circle [109].
• One big arena of physics at nonzero temperature that we have not even touched upon in this paper is transport. Quantities like the diffusivity, viscosity, electrical conductivity, and thermal conductivity are basic experimentally measurable quantities that provide a wealth of information about the low-energy excitations of a system. These transport coefficients have well-known expressions in terms of two-point functions of components of conserved currents or the stress-energy tensor [110][111][112]. The most interesting limit of the thermal two-point functions for transport phenomena is the low frequency limit, which translates to large separations in position space.
Apart from weak coupling expansions, transport has been exhaustively studied from a holographic perspective: for a recent review, see [113].
While the OPE of the thermal two-point function strongly constrains the short distance dynamics in the CFT, it does not directly constrain the long-distance behavior due to the absence of any OPE channel for |x| > β. It is easy to derive functional forms for correlators in the diffusive regime via hydrodynamics, which is the correct low-energy description [114]. Can bootstrap techniques allow us to derive this specific form of the diffusive correlator, and the value of the energy diffusion constant for the 3D Ising model? It would be very interesting to connect the OPE regime to the hydrodynamic regimes in a CFT.

JHEP10(2018)070
As discussed in section 2.1.3, estimating thermal one-point functions by taking a limit of correlation functions on S 1 β × S d−1 is challenging. In general, one needs to know the spectrum ∆ O and OPE coefficients f OO O for arbitrarily high dimension operators O (not to mention the one-point blocks for all tensor structures appearing in OO O ). In the next appendix, we give slightly more detail in d = 2.
However, in any d, the observable b T is special in that it depends only on the spectrum of the theory. 42 This is because the expectation value of the stress-tensor on S 1 β × S d−1 is proportional to a derivative of the partition function, The partition function can be expanded in characters where ∆, ρ are the dimension and SO(d) representation, respectively, of O and we sum over primaries only. In practice, even if we don't know the full spectrum of a theory, we can try to estimate b T by truncating the sum over characters at some ∆ max . More precisely, let us define We can then try to extrapolate g ∆max (β) towards β = 0. The actual value of g ∆max (0) will always be 0, because β d will dominate over the contribution of a finite number of characters at sufficiently small β. However, perhaps we can estimate b T by evaluating g ∆max (β) at a small, nonzero value of β. As a check on this idea, let us study the free boson, where we know the spectrum exactly. For concreteness, we work in d = 3. The partition function is given by Z free (q) = ∞ j=0 1 (1 − q j+1/2 ) 2j+1 , (A.5) 42 We thank Chris Beem, Scott Collier, Liam Fitzpatrick, and Slava Rychkov for discussions that inspired the calculations in this appendix.

JHEP10(2018)070
where q = e −β . 43,44 It can be decomposed into conformal characters as Z free (q) = 1 + χ free (q) + =2,4,... χ short (q) + Z long (q), χ free (q) = χ 1/2,0 (q) − χ 1/2+2,0 (q), χ short (q) = χ +1, (q) − χ +2, −1 (q), Here, χ ∆, (q) is the character of a long multiplet, and the first three terms in Z free (q) correspond to the unit operator, the boson φ itself, and a tower of higher-spin currents. The long multiplet content is Z long (q) = χ 1,0 (q) + χ 3/2,0 (q) + χ 2,0 (q) + . . . . (A.7) To determine the quantum numbers and multiplicities of long multiplets, we can include a fugacity for angular momentum and decompose the full partition function with this fugacity into conformal characters. This is a standard exercise and we do not include the details here. Using our knowledge of the spectrum, we can plot g ∆max (β) for various values of ∆ max in the free boson theory ( figure 12). The function with ∆ max = ∞ (black dotted line) decays as e −β/2 for large β (coming from the contribution of the lowest-dimension operator φ). It reaches a minimum near β = 5, and then smoothly approaches the value b T = − 3ζ (3) 2π ≈ −0.574 as β → 0. The curves with finite ∆ max move closer to the ∆ max = ∞ curve, with longer and longer plateaus near b T before eventually going to 0 at β = 0. The 3d Ising model is a nonperturbative theory where we don't know the full spectrum, but we do know a large part of it to reasonable precision from numerical bootstrap computations [14,16,[27][28][29][30]116]. In particular, the spectrum of operators appearing in the σ ×σ, σ × , and × OPEs are known up to dimension 8 [30]. Some additional low-twist families are known up to very high dimension, but these are a small portion of the highdimension spectrum. The lowest-dimension operator not appearing in the above OPEs is expected to be a Z 2 -even vector with dimension approximately 6, though its dimension is not known to high precision [117]. Thus, our knowledge of the spectrum begins to fade when ∆ max ≈ 6. Nevertheless, in figure 13, we estimate g ∆max (β) by including the known operators with dimension ∆ ≤ ∆ max .  -59 -

JHEP10(2018)070
Despite our ignorance of the high-dimension spectrum, the plot in figure 13 already shows similar structure to the free scalar case, with a plateau beginning to form near b T ≈ −0.45, close to the value b T ≈ −0.459 determined via Monte Carlo simulations [58][59][60]. It would be interesting to understand whether figure 13 can be turned into a rigorous estimate, perhaps by understanding better the analytic structure of g ∞ (β). It would also be interesting to understand whether the existence of a minimum, seen as the "dip" in these plots, is a feature shared by all CFTs. . We focus here on d = 2 for simplicity. In this case, the thermal blocks factorize, where |f (h)| 2 ≡ f (h)f (h), and (e.g. [118,119]) where q ≡ e −β . (B.1) generalizes in the obvious way to unequal left-and right-moving temperatures. Using the connection formulae for hypergeometric functions, we may rewrite the left-moving block (for generic (h O , h O )) as As β → 0, there are two branches: Combining the left-and right-moving blocks yields the full scaling behavior of the torus one-point blocks at high temperature. 45 It is remarkable that, for h O ≥ 1/2, the leading term in (B.4) exhibits the same scaling of the full one-point function on S

JHEP10(2018)070
This leads to a formal expression for the S 1 β × R d−1 one-point function b O as a sum over states, in the limit of high temperature: Note that the summand is not, in general, sign-definite, and receives contributions from all spins and arbitrarily high energies. 46 On the other hand, for low-dimension operators O, with h O < 1/2, the second term of (B.4) dominates, and even the recovery of the requisite β −2h O scaling at high temperature is sensitive to the details of the full sum. The dimensional reduction of a d-dimensional CFT on S 1 does not always give a welldefined theory in d − 1-dimensions. For example, a problem occurs if we try to compactify the free boson CFT in 3d down to 2d (with periodic boundary conditions around the S 1 ). 47 Naively, we should get the 2d free boson with noncompact target space, but this theory is pathological because correlations grow logarithmically with distance. Another way to see the problem is that if we try to compute the propagator using the method of images, the sum over images diverges. 48 We expect this issue to arise whenever we compactify a 3d CFT with a nontrivial moduli space of vacua down to 2d, as long as the boundary conditions do not destroy the moduli space. For example, supersymmetric compactifications of 3d SCFTs with Higgs or Coloumb branches should be treated with care. One way to study such theories is to introduce twisted boundary conditions that remove the zero mode from the path integral. 49 Correlation functions in the twisted setting then share many similar properties to those we discuss in this work (for example an OPE and crossing equation). It would be interesting to adapt the techniques in this work to deal with general twisted boundary conditions. In our case, thermal compactification does the job because it breaks supersymmetry and generically lifts the moduli space. (E.1) As with the main text, h f = 2h φ stands for the asymptotic h of the family, and h e = 2h φ stands for the total h of the two external operators. Self-corrections of a [φφ] 0 take the 47 Conformal invariance requires that the free boson CFT have a noncompact target space in 3d. 48 It is interesting to ask what happens if we have a physical system that has an EFT description in terms of a 3d free boson, and we place it at finite temperature. In this case, the thermal physics can depend on UV details that are not directly captured by the 3d effective CFT. For example, if the 3d boson is the Goldstone boson of a broken U(1) symmetry with symmetry breaking scale Λ, its dimensional reduction is better described as a 2d boson with compact target space, where the radius is R ∝ √ βΛ. 49 We thank Nati Seiberg for discussions on this point. 50 If desired, we could also start the sum at a higher spin, and treat the contributions of the low-spin operators exactly in their anomalous dimensions. The corresponding generalization of our algorithm is straightforward. If the anomalous dimensions are small even for low-spin operators, we can safely start the sum at low spins. In practice, one can work order-by-order in the small anomalous dimensions of the family, effectively truncating the m and n to some finite order, thus avoiding having to compute an infinite number of coefficients and to invert an infinite matrix. Finally, let's note that it is possible to generalize this method to the [φφ] n families, and even potentially to considering collections of families at once.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.