A Numerical Approach to Virasoro Blocks and the Information Paradox

We chart the breakdown of semiclassical gravity by analyzing the Virasoro conformal blocks to high numerical precision, focusing on the heavy-light limit corresponding to a light probe propagating in a BTZ black hole background. In the Lorentzian regime, we find empirically that the initial exponential time-dependence of the blocks transitions to a universal $t^{-\frac{3}{2}}$ power-law decay. For the vacuum block the transition occurs at $t \approx \frac{\pi c}{6 h_L}$, confirming analytic predictions. In the Euclidean regime, due to Stokes phenomena the naive semiclassical approximation fails completely in a finite region enclosing the `forbidden singularities'. We emphasize that limitations on the reconstruction of a local bulk should ultimately stem from distinctions between semiclassical and exact correlators.


Introduction and Summary
Many of the most challenging conceptual problems in theoretical physics were only resolved after physicists discovered how to 'shut up and calculate' a large variety of observables to high precision. For example, our modern understanding of quantum field theory was only developed after the physics community had decades of experience with perturbative calculations. And it is hard to imagine how decoherence could have been understood without the temporary crutch provided by the Copenhagen interpretation and its instrumental approach to the Born rule.
Though we have struggled with the black hole information paradox for decades, major progress has been possible through the development of AdS/CFT. Resolving the information paradox in AdS/CFT will require a precise understanding of bulk reconstruction and its limitations. Although reconstruction presents thorny conceptual problems, the limitations on reconstruction should ultimately stem from discrepancies between the predictions of gravitational effective field theory in AdS and conformal field theory. This means that to make progress, it will be crucial to be able to directly compare the approximate correlation functions of bulk EFT and the exact correlators of the CFT.
AdS 3 /CFT 2 may provide the best opportunity for such comparisons. Many features of quantum gravity in AdS 3 can be understood 'kinematically' as a consequence of the structure of the Virasoro algebra. To be specific, the Virasoro conformal blocks have a semiclassical large central charge limit that precisely accords with expectations from AdS 3 gravity, reproducing the physics of light objects probing BTZ black holes. In the semiclassical approximation, the Virasoro blocks exhibit information loss in the form of 'forbidden singularities' and exponential decay at late times [1][2][3][4][5][6]. Moreover, these problems can be partially addressed by performing explicit analytic calculations [7]. The blocks can also be computed directly from AdS 3 [8][9][10][11][12][13][14][15].
In this work we will investigate the discrepancies between semiclassical gravity and the exact CFT by computing the Virasoro blocks numerically to very high precision. This is possible via a slightly non-trivial implementation of the Zamolodchikov recursion relations [16][17][18]. We discuss the blocks and the algorithm in detail in section 2 and appendix A. For the remainder of the introduction we will explain the physics questions to be addressed and summarize the results.

When is the Semiclassical Approximation Valid?
The Virasoro conformal blocks have a semiclassical limit. CFT 2 correlators can be written in a Virasoro block decomposition as It is natural to ask about the range of validity of this approximation -how does it depend on the kinematic variable z and the ratios h i /c and h/c?
One reason to ask is simultaneously speculative and pragmatic -one might like to know if it is possible to explore AdS 3 quantum gravity in the lab by engineering an appropriate CFT 2 (for a concrete idea see [19]). But gravity will only be a good description if the semiclassical limit provides a reasonable approximation at accessible values of c. Unfortunately, even in the semiclassical limit the Virasoro blocks are not known in closed form for general parameters. But we can partially test the validity of this limit by computing c 2 c 1 log V(c 1 ) log V(c 2 ) for c 2 ≈ c 1 , as this ratio will be 1 when the semiclassical limit holds. We plot this ratio in figure 7, which shows that the blocks adhere to the semiclassical form of equation (1.2) remarkably well (up to an important caveat to be discussed later).
We would also like to understand if the semiclassical limit breaks down in specific kinematic regimes associated with quantum gravitational effects in AdS 3 . This is what we will explore next.

Information Loss and OPE Convergence in a New Regime
Information loss can be probed using the correlators of light operators in a black hole background [20], as illustrated in 1. When computed in a BTZ or AdS-Schwarzschild geometry, these correlators decay exponentially as we increase the time separation t between the light operators, even as we take t → ∞. This behavior represents a violation of unitarity for a theory with a finite number of local degrees of freedom on a compact space. Thus it is interesting to see how it is resolved by the exact CFT description.
As a first step, one would like to understand how unitary CFTs are able to mimic bulk gravity, including the appearance of information loss. This arises from the heavylight large central charge approximation [5][6][7] of the Virasoro blocks. Thus it seems to be very universal, as it is largely independent of CFT data. The next step is to understand how finite c physics corrects this approximation, and what CFT data are involved in resolving information loss.
It is useful to think about the Fourier representation of both the correlator and the individual Virasoro blocks. For the full correlator this is where λ(E) is the OPE coefficient density and we have taken z = 1 − e −it to study Lorentzian time separations between the light probe operators. At large t we probe the fine structure of λ(E), which means that the least analytic features of λ(E) dominate the late time limit. Practically speaking, this means that the very late time limit probes the discrete nature of the spectrum, and we become sensitive to the fact that λ 2 (E) is a sum of delta functions. At early or intermediate times we only discern the coarse features of λ(E). There are at least five different timescales associated with black holes in AdS/CFT.
The inverse temperature β = 2π |α H | where α H ≡ 1 − 24h H c sets the shortest relevant scale, where h H is the holomorphic dimension of O H . The scale β log c estimates the time it takes for infalling matter to be scrambled [21,22]. At times of order the entropy S = π 2 c 3β , heavy-light correlators cease their exponential decay; this is also the evaporation timescale for black holes in flat spacetime. We expect that the typical energy splittings among neighboring eigenstates to be of order e −S , which means that at times of order e S we will be sensitive to the discreteness of the spectrum. Finally, on timescales of order e e S the phases of the eigenestates can come back into approximate alignment, leading to recurrences.
As discussed in detail in section 3, what we find is that the Virasoro blocks with h L < c 24 < h H behave very differently at early and late times, as was presaged by analytic results [7]: • The blocks with intermediate operator dimension h c 24 are well-described by their semiclassical limit [1,2,5] for When h > c 24 the blocks are also well-described by the semiclassical limit at early times, but we do not have a precise formula quantifying 'early'.
• Heavy-light blocks with h h H initially grow, as was found from a semiclassical analysis [5]. We find that they reach a maximum and then subsequently decay. The factor 5 2 comes from empirical fits; the function A t ( c h H ) is always order one and is approximately linear in c h H . Other sub-leading behavior is discussed in section 3.
• Numerical evidence indicates that all heavy-light Virasoro blocks decay as at late times, independent of h, as long as h H > c 24 > h L . We present evidence that this decay persists beyond the exponentially long timescale ∼ e S , so we believe that it represents the true asymptotic behavior of the heavy-light blocks.
From the point of view of the 1 c ∝ G N expansion, the universal late-time power-law decay comes from non-perturbative effects. If this behavior persists to all times, as our empirical evidence indicates, then the late time behavior of CFT 2 correlators must come from an infinite sum over Virasoro blocks in the heavy-light channel. 1 From a pure CFT perspective, the late Lorentzian time behavior represents a new limit in which the bootstrap may be analytically tractable. Most analytic bootstrap results, including the Cardy formula [23], OPE convergence [24], and the lightcone OPE limit [25,26] arise in a similar way. In fact, because the expansion of CFT 2 correlators in the uniformizing q-variable converges everywhere, including in deeply Lorentzian regimes, it affords the opportunity to explore many new 'analytic bootstrap' limits.

Forbidden Singularities and Bulk Reconstruction
It is interesting to understand when exact CFT correlators differ markedly from predictions obtained from a semiclassical AdS description. The late time regime we discussed Figure 2. This figure shows the Penrose diagram for an energy eigenstate black hole in AdS, suggesting the role of ingoing and outgoing modes behind the horizon and their relationship with local CFT operators. Analytic continuation provides a painfully naive but instrumentally effective method for studying correlators behind the horizon.
above provides one example of this phenomenon. As we discuss here and in section 4, there are also Euclidean regimes where the semiclassical approximation to the Virasoro blocks fails completely. Correlation functions in CFT 2 must be non-singular away from the OPE limits where local operators collide [7,27]. The Virasoro conformal blocks must have this same property [5]. But in the semiclassical approximation, the blocks develop additional 'forbidden singularities' [7] that represent a violation of unitarity. These singularities are a signature of semiclassical black hole physics in AdS 3 . They arise because thermal correlators exhibit a Euclidean-time periodicity under t → t + iβ, and so the OPE singularities have an infinite sequence of periodic images. The exact Virasoro blocks are not periodic, but in the semiclassical approximation they develop a periodicity at the inverse Hawking temperature β = 1 T H associated with a BTZ black hole in AdS 3 . By studying the Virasoro vacuum block in the vicinity of potential forbidden singularities, one can show that at finite c the singularities are resolved in a universal way [7] via an analytic computation. This method predicts the kinematic regimes where non-perturbative effects should become important; it can be extracted from equation 4.2 and the results are displayed in figure 17. Thus it is interesting to investigate the divergence between the exact blocks and their semiclassical approximation more generally. We study this question numerically in section 4.
Discrepancies between exact and semiclassical CFT correlators near the forbidden singularities could have implications for the reconstruction of AdS dynamics. Bulk reconstruction in black hole backgrounds is rather subtle [28][29][30][31], and perhaps requires some understanding of the analytic continuation of CFT correlators. But there is also a very simple and physical reason to expect that the analyic properties of correlators could have something to do with black hole interiors [28].
As emphasized by Raju and Papadodimas [32,33], a field operator behind the horizon consists of both ingoing and outgoing modes, but only the ingoing modes can be immediately associated with local CFT operators. This issue is portrayed in figure 2. The analytic continuation of local operators by t → t + iβ 2 provides a naive, instrumental source for the outgoing modes. 2 Thus it is natural to ask whether the exact and semiclassical correlators differ significantly at t + iβ 2 , which is 'halfway' to the first forbidden singularity.
We will observe in section 4 that the exact and semiclassical correlators behave very similarly at these points, though they seem to differ markedly both very near (within 1 √ c ) and beyond the first forbidden singularity. The results can be seen in figure 15. As previously discussed [7], we expect that Stokes and anti-Stokes lines emanate from the locations of the forbidden singlarities, so that different semiclassical saddle points dominate in different regions of the q-unit disk. It appears that different saddles dominate as we cross the locations of the forbidden singularities, so that the naive semiclassical blocks (the saddles that dominate near q = 0) are not a good approximation beyond the first singularity. In fact the semiclassical approximation appears to break down in a finite kinematic region, as shown in figure 17. Furthermore, the existence of such regions seems to depend in an essential way on the presence of a black hole, ie a state with energy above the Planck scale (h H > c 24 ), as semiclassical/exact agreement is excellent when h H < c 24 , as we see in figure 18. Perhaps future investigations will uncover bulk observables that are sensitive to Stokes phenomena in the large c expansion of the Virasoro blocks. We hope that the black hole information paradox can be understood with more precision and detail through such calculations. This work takes steps in that direction by identifying new kinematic regimes where the semiclassical limit breaks down badly and by providing results for the correct non-perturbative Virasoro blocks.

Kinematics, Convergence, and the Semiclassical Limit
A great deal of information about the behavior of CFT 2 correlation functions is encoded in the structure of the Virasoro conformal blocks. We are interested in 4-pt correlators of primary operators, which can be written as where the P h,h are products of OPE coefficients. The V h i ,h,c (z) are the holomorphic Virasoro blocks, which will be the main object of study in this work. The blocks are uniquely fixed by Virasoro symmetry and depend only on the external dimensions h i , the exchanged primary operator dimension h, and the central charge c. Often it will convenient to write z = 4ρ (1+ρ) 2 so that the full z-plane lies inside the ρ unit circle [24]. The Virasoro blocks are not known in closed form, but they can be computed orderby-order in a series expansion using recursion relations. We provide a brief summary here, leaving the details to appendix A.
There are two versions of the Zamolodchikov recursion relations (for a nice review see [35]). The first [16] is based on writing V h i ,h,c as a sum over poles in the central charge c, plus a remainder term that survives when c → ∞ with operator dimensions fixed. The second [17], which is more powerful, arises from expanding the blocks as a sum of poles in the intermediate dimension h plus a remainder term that survives as h → ∞. The remainder term can be computed from the large h limit of the Virasoro blocks [17,18]. This large h limit of the blocks takes a simple form when written in terms of the uniformizing variable where K(z) is the elliptic function The q-coordinate can be derived from the accessory parameter/monodromy method in the semiclassical limit [36] or from a quantization of the theory on the pillow metric [27]. It has the remarkable feature that q(z) covers the full multisheeted z-plane (the sphere with punctures at 0, 1, ∞), as depicted in figure 3. The Virasoro blocks can then be written in the form where H (c, h, h i , q) can be obtained from the recursion relation: We note that this recursion relation naturally produces a series expansion in the variable q. For more details along with the definitions of the quantities appearing in these equations see appendix A.
In this work, we will be using the recursion relations to obtain the q-expansion of the Virasoro blocks to very high orders. It appears that most prior implementations of the Zamolodchikov recursion relations could not reach the N ∼ 1000 that we will study. 3 Our improvements are fairly elementary, and are based on computing and storing the specific coefficients of powers of q in H(c, h m,n + mn, h i , q), as we describe in more detail in appendix A. The computational time complexity of our algorithm is roughly O(N 3 (log N ) 2 ), while it seems that some earlier implementations scaled exponentially with N . The maximum N is limited by memory consumption, with memory usage scaling roughly as O (N 3 log N ). We have verified our code by comparing to a number of previous results, including prior implementations, the semiclassical blocks, blocks computed by brute force from the Virasoro algebra, and the special case of degenerate external operators.

Kinematics and Convergence of the q-Expansion
Both the correlator and the Virasoro blocks in equation (2.1) can have singularities in the OPE limits, which occur when z → 0, 1, ∞. Generically we expect branch cuts in the z-plane running between these three singularities. So for our purposes, the most remarkable feature of the variable q(z) is that the region |q| < 1 covers not only the complex z-plane, but also every sheet of its cover. The relationship of the z plane and its branch cuts to the region |q| < 1 [27] is depicted in figure 3. The Zamolodchikov recursion relations provide an expansion for the Virasoro blocks that converges for all |q| < 1, which means that they can provide a good approximation to the 4-pt correlator in any kinematic configuration. In particular, we can use the q-expansion to study the Lorentzian regime with arbitrary time-orderings for the operators.
The existence of the q-variable implies that in CFT 2 , there are an infinite number of distinct regimes where the bootstrap equation may be analytically tractable. In the case of d ≥ 3, one can study the OPE limit z → 1 using conformal blocks expanded in  Figure 3. The q(z) map takes the universal cover of the z-plane (the sphere with punctures at 0, 1, ∞) to |q| < 1. This figure suggests the relationship between the z plane, the unit ρ disk, and the unit q disk, with branch cuts indicated with colored lines [27]. The relations between these variables are q = e and z = 4ρ (1+ρ) 2 , and the inverse transformations are z = θ 2 (q) The Virasoro blocks converge throughout |q| < 1, with OPE limits occurring on the q unit circle.
the OPE limit of small z, and this implies various exact results about the properties of large spin operators [25,26]. However, because the Euclidean OPE in d ≥ 3 does not converge deep in the Lorentzian region, one cannot study other OPE channels in the same way. This obstruction disappears in d = 2, where one must be able to reproduce all of the distinct OPE limits |q| → 1 pictured in figure 3 using the small q expansion. The large Lorentzian time limit that we will discuss in section 3 provides a physically motivated example of this idea.
We will be studying numerical approximations to the Virasoro blocks based on a large-order expansion in the q variable. Thus to understand the convergence properties of our expansion, it may be useful to map out the regions of constant |q|. For this purpose we can use the coordinate ρ(z) defined via z = 4ρ (1+ρ) 2 [24], because the entire z-plane can be easily visualized as the region |ρ| < 1. The operators at z = 1 and ∞ are mapped to ρ = 1 and −1, respectively. In figure 4 we have plotted contours of constant |q| in the ρ-coordinate system. In figure 5 we present results on the convergence region of the q-expansion of the Virasoro blocks for various values of the dimensions and central charge.
A kinematic configuration that will be of particular interest represents z = 1−re −it (andz = 1−re −it as well) and is depicted in the AdS/CFT context in figure 1. With this Since this is only the first sheet of the z-plane, it corresponds to the region in the q-disk enclosed by the two blue lines connecting ±1 in figure  3. The correlator can have singularities in the OPE limits ρ → −1, 1 and these correspond to q → −1, 1 as well. Away from these limits |q| < |ρ| and the q-expansion converges much more rapidly than the ρ expansion.
setup we can study the correlator of light operators O L (z)O L (0) at timelike separation in the background created by a heavy operators O H . At large times t, this correlator can be used as a probe of information loss in pure state black hole backgrounds, as we will discuss in section 3. On the z plane, the late time behavior is obtained by analytically continuing the conformal block around the branch-cut starting at z = 1 multiple times. Explicitly, the Lorentzian value of the q variable is obtained with the following analytic continuation of the elliptic integral: where the elliptic functions on the right-hand side are evaluated on the principle sheet with the branch-cut chosen to be z ∈ [1, ∞). At large t we have with the elliptic function K(z) are taken on their principle sheet, so that g(r, t) is periodic in t. This means that |q| 2 ≈ 1− π 3 t 2 (g+g * )+· · · and the real part Re[g(r, t)] > 0, so that |q| < 1 for all t, as expected. In the limit that r 1, we have g(r, t) ≈ 1 π log 16 r , which leads to the estimate in the limit of r 1 and t → ∞. Thus we can translate between convergence in |q| and t; very roughly, we expect that working to order q N will allow us to probe t ∝ √ N at large N . We can visualize the trajectory of q(r, t) for various r in figure 6.

Review of Blocks and Adherence to the Semiclassical Form
Much is known about the Virasoro blocks in various limits. In the limit c → ∞ with all dimensions held fixed, the Virasoro blocks simply reduce to global conformal blocks, which are hypergeometric functions. Corrections to this result up to order 1/c 3 are known explicitly [39]. In the heavy-light limit, where we take c → ∞ with two 'heavy' where T H is the Hawking temperature of a corresponding BTZ black hole. In the case of the vacuum block, which is h = 0, the 1/c corrections to this limit are also known explicitly [4] for any h H /c. Finally, in the semiclassical large c limit, where all dimensions h i , h ∝ c, there is overwhelming evidence that the blocks take the form as though they are derived from a semiclassical path integral (and in fact they have an sl(2) Chern-Simons path integral representation [14]). The semiclassical saddle points have been classified [5], and in some kinematic limits we can determine the behavior of f analytically. In particular, the large Lorentzian time behavior of f with the kinematics of figure 1 and h L < c 24 < h H has been determined [5]. The result is that the leading semiclassical contribution always decay exponentially at sufficiently large times 4 at the rate and α H = 2πiT H with T H the corresponding Hawking temperature. As we will review in section 3, this demonstrates that information loss due to black hole physics [20] occurs as a consequence of the behavior of the individual Virasoro blocks [5,7]. Finally, some exact information about the behavior of the Virasoro blocks can be obtained by studying degenerate states [7].
Most of these approximations hold in the large central charge limit when the kinematic configuration is held fixed. But the deviations between the exact and semiclassical Virasoro blocks may depend importantly on the kinematics. As we will discuss in detail below, the semiclassical blocks have 'forbidden singularities' that are absent from the exact blocks [7]. We also find that as expected [5,7], the exact and semiclassical blocks have very different behavior at large Lorentzian times. More generally, we would like to map out the kinematic regimes where non-perturbative corrections to the semiclassical Virasoro blocks become large. But at a more basic level, it is interesting to ask how large c must be before the semiclassical limit of the Virasoro blocks provides a reasonable approximation to their behavior. This has immediate implications for the possibility of constructing a 2d CFT and probing quantum gravity in an experimental lab. A natural way to probe the existence of the semiclassical limit is by studying the ratio of logarithms of blocks at somewhat different central charges c 1 and c 2 . If the semiclassical limit of equation (2.11) is a good approximation, then this quantity will be 1, but otherwise we expect it to deviate from 1 by effects of order 1 c . In figure 7 we explore this ratio and find that the semiclassical form V ≈ e cf provides a remarkably good approximation for very small values of c.
There is an important caveat that we will return to in section 4. An infinite number of distinct semiclassical saddle points can contribute to the Virasoro blocks in the large c limit [5]. Thus it is possible that V ≈ e cf for some f , but that due to Stokes phenomena, the dominant saddle f changes as we move in the q unit disk. So although the semiclassical limit may appear to describe the blocks well for all q, as indicated by  11 10 c in order to test the semiclassical limit. As we increase c, the semiclassical limit becomes a better approximation and R → 1, but even for c = 2.1 the blocks are remarkably well approximated by the semiclassical form. For the larger choices of c the functions have similar shapes up to an overall rescaling; this suggests that the first 1/c correction is dominating the discrepancy R − 1. In the OPE limit q → 0 the semiclassical limit always applies. We find similar results for non-vacuum blocks. figure 7, in fact the saddle that is leading near q ≈ 0 may be sub-leading at general q. Thus the naive semiclassical blocks may differ greatly from the exact blocks; in fact we will find this to be the case in section 4.
Nevertheless, figure 7 suggests that we should be optimistic about probing semiclassical CFT 2 correlators in the lab! It would be very interesting to engineer a CFT 2 with c > 1 and no conserved currents aside from the stress tensor [19].

Late Time Behavior and Information Loss
One sharp signature of information loss in AdS/CFT [20] is the exponential decay of correlation functions at large time separations in a black hole background. This can be studied using heavy-light 4-point functions in CFT [1]. As portrayed in figure 1, we can interpret this correlator as the creation and subsequent measurement of a small perturbation to an initial high-energy state. In a unitary theory on a compact space with a finite number of local degrees of freedom, this initial perturbation cannot completely disappear. But a computation in the black hole background displays eternal exponential decay, capturing the physical effect of the signal falling into the black hole. At a more technical level, the exponential decay rate can be obtained from the quasinormal mode spectrum of fields propagating in the black hole geometry.
The simplest way to see that heavy-light correlators cannot decay forever is to where λ 2 (E) is a product of OPE coefficients. Because the sum on the right-hand side is discrete, the correlator must have a finite average absolute value at late times. When h H c 24 h L , we expect the states contributing in (3.9) to be a chaotic collection of e S blackhole microstates with energy near that of O H , and with S = π 2 3 T H c. The amplitude will initially decay due to cancellations between the essentially random phases, but these cancellations cannot become arbitrarily precise. Roughly speaking, the decay should stop when the correlator reaches ∼ e −S and begins to oscillate chaotically. At a more detailed level, the time dependence can change qualitatively on timescales of order S and e S as different features of λ 2 (E) come into play [40][41][42].
In this work, we will not study the which is related to the first channel by the bootstrap equation (or by modular invariance in the case of the partition function [42]). In this channel we are sensitive to the exchange of states between the heavy and light operators. For example, pure 'graviton' states in AdS 3 correspond to the exchange of the Virasoro descendants of the vacuum, which are encapsulated by the Virasoro vacuum block. Other heavy-light Virasoro blocks include a specific primary state along with its Virasoro descendants, which one can think of as gravitational dressing. We are interested in this channel because heavy-light Virasoro blocks encode many of the most interesting features of semiclassical gravity. We would like to understand to what extent the exact Virasoro blocks know about the resolution of information loss.
It is convenient to think of the time dependence of the Virasoro blocks as coming from a potentially continuous λ h (E) associated with each block, via where h labels the dimension of an intermediate Virasoro primary OPEs. Roughly speaking, the late time dependence of V h (t) will come from the least analytic features of λ 2 h (E). For example, in the leading semiclassical limit, heavy-light Virasoro blocks decay exponentially at late times at a universal rate given in equation (2.12). This semiclassical behavior comes from a function λ 2 h (E) that is smooth on the real axis, but has poles in the complex E-plane. In AdS 3 these poles can be interpreted as the quasinormal modes of a BTZ black hole background (at least for small h). A straightforward contour deformation of equation (3.2) connects these poles to the exponential decay.
At sufficiently late times, the physics of the quasinormal modes will be subdominant to less analytic features in λ 2 h (E). For example, if λ 2 h (E) exhibits thresholds of the form (E − E * ) p−1 with E * real, then V(t) will inherit a power-law behavior t −p at late times. And if λ 2 h (E) receives delta function type contributions, then V(t) will have a finite average absolute value at late times. If such features are present in V h (t), then it is natural to investigate the timescale where V h (t) transitions from exponential decay to some other late-time behavior.
The full CFT 2 correlator should not become much smaller than ∼ e −S . Since Virasoro blocks associated with light operators initially decay exponentially, one might naively expect that V h (t) should change qualitatively after a time of order S. More specifically, for heavy-light correlators dominated by the vacuum block, we would expect a departure from exponential decay by a time up to an unknown order one factor. This argument is rather weak, since the full correlator might not behave like the light-operator Virasoro blocks. However, the same prediction for t D was derived from an analysis of non-perturbative effects [7] in the vacuum block. We discuss the equation that led to that prediction in section 4.3.
We will see empirically that Virasoro blocks with small h do undergo a transition at a timescale remarkably close to t D . Furthermore, at late times the behavior of the heavy-light Virasoro blocks appears to be a universal power-law: where we require h H ≥ 1 24 , so that at least one external operator is heavy enough to create a blackhole. When the intermediate dimension h h H the late time powerlaw behavior remains the same, although the transition time then also depends on h (and we do not have an analytic prediction to compare to). This universal behavior suggests a threshold √ E − E * in λ 2 h (E), which seems to correspond with random matrix behavior [41,43]. Our results indicate that the t − 3 2 power-law persists to timescales ∼ e S , so individual heavy-light Virasoro blocks are not sensitive to the discreteness of the spectrum.
These results show that the time-dependence of the heavy-light Virasoro blocks has some qualitative similarities with that of the Virasoro vacuum character after an S transformation and the analytic continuation β → β + it [42]. Both the heavy-light blocks with small h and the vacuum character have an initial exponential-type decay, though the precise time-dependence is rather different. The heavy-light blocks and the vacuum character have the same power-law decay at late times, though non-vacuum characters decay with a different late-time power-law [42].
In what follows we will study the heavy-light blocks V h (t) empirically to establish the robust features of their time-dependence. We also translate the late-time t −3/2 behavior into a statement about the coefficients of q N in V h (q) at large orders in the q-expansion, as one might hope to derive this asymptotic behavior for the coefficients using the Zamolodchikov recursion relations. One might also compute λ 2 h (E) directly using the crossing relation [44,45]. Finally we discuss the implication of our results for the late time behavior of the correlator.

Vacuum Virasoro Blocks
Using the methods discussed in section 2, we compute the vacuum Virasoro blocks at late times. Figure 8 shows the result along with a comparison to the semiclassical blocks computed using semi-analytic methods [5]. For numerical convenience we avoid certain rational values of c to prevent singularities in intermediate steps of the computation.
Using the numerical result of the full Virasoro blocks, we can measure the departure time t d when the semiclassical block drops below the exact block. We compare this measured value to the prediction of (3.3) in figure 9. The logic leading up to (3.3) is only valid parametrically, so it is remarkable that it agrees with the measured t d up to a small constant shift. Note that we parameterize the time dependence via z = 1 − re −it , and this constant shift depends on r. We have also checked that t d is primarily controlled by the ratio h L c and has a very weak dependence on h H and c. Around the time t D , all vacuum blocks show an obvious change of behavior from an initial exponential decay to a much slower power law decay. To very good accuracy, the power of this decay seems to be t − 3 2 universally in all of the parameter space we were able to explore with an external operator with dimension h H > 1 24 .
A few examples are provided in figure 8, but we tested this behavior with hundreds of different parameter choices.

General Virasoro Blocks
The non-vacuum blocks also exhibit universal t − 3 2 late-time decay. The difference from the vacuum case is that we no longer have a simple estimate for the time scale of the transition. In particular, we find that generically the non-vacuum blocks grow at early times, reach a maximum at time t max , and then start to decay, finally settling down to the t − 3 2 power law behavior. These features are illustrated by examples in figure 10. From the data plotted in figure 11, we see that beyond the blackhole threshold h > c 24 , the timescale t max has a simple dependence on parameters. We can fit it to the ansatz with α h = 1 − 24h c and obtain A t and b time empirically. The parameter A t is almost a linear function of c h H , as can be seen in figure 12 figure  22 in the appendix, which displays the variation in A t .
On the left of figure 11 we plot |Ṽ max |, which is the maximum of the absolute value of the block after extracting a universal prefactor via |V max | = 16 h− c−1 24 |Ṽ max |. We see that |Ṽ max | also has a simple dependence on h c . We can perform a similar fit for |Ṽ max |,  and t max as linear functions of log h c and |α h | respectively (we've used points with h c = n 3 for n = 1, 2, · · · , 30). We find that both plots are robustly c-independent for c 5, as expected in the semiclassical limit. We see explicitly that there is little dependence on h L ; in the a height plot the variation with h L is almost invisible. and we find that Empirically we obtain a height ≈ −2.5 from the fit in figure 12. The b time and b height parameters do not fit a simple pattern; we provide some data on these parameters in figure 21 in the appendix. These fits led to the result summarized by equation (1.5) in the introduction, which neglects the small offsets from the b-parameters. We expect that |V max | and t max are controlled by semiclassical physics (for example, see figure  20), so it would be interesting to try to prove these empirical relations using analytic results [5] on the semiclassical time-dependence. In principle these results could also be obtained from an AdS calculation involving black holes and deficit angles.

Probing Exponentially Large Timescales
Formally, we are interested in high-energy pure states corresponding to BTZ black holes, which have a large entropy S = π 2 3 cT H in the large central charge limit where AdS gravity provides a reliable description. This suggests that timescales of order e S will be unreachably large. Nevertheless, by considering either small c or small T H , we can probe order one S, and thus reach t ∼ e S within the range of convergence of our numerics.
In fact, the plot on the bottom-right of figure 8 is already in this regime. Due to its low temperature of T H ≈ 0.03 in AdS units, we have S ≈ 3.3 so that times of order e S ≈ 27 are within the range of convergence. Thus this plot already suggests that the t − 3 2 power-law decay persists to exponentially large timescales. In figure 13 we have displayed four other choices of parameters where timescales of order e S , and even e e S , are visible within the range of convergence. Two examples have order one T H and small c, while two others have very small T H and relatively large c. In all cases we see that the t − 3 2 late-time decay persists on these exponentially large timescales. This provides good evidence that the heavy-light Virasoro blocks really do decay in this way at very late times. This means that these blocks are not sensitive to the discreteness of the spectrum in other channels.

Power Law Behavior of q-Expansion Coefficients
We have observed an apparently universal late-time power-law behavior in the heavylight Virasoro blocks V h (t). One might try to derive this behavior by studying its implications for the q-expansion. In fact, for a large region of parameter space, the t − 3 2 decay translates to a power law growth of the coefficients in the q expansion.
To see this, we note that at late times q approaches 1 with a rate given by (2.7). This implies that θ 3 (q) ∼ √ t, which means that the prefactor in (2.4) behaves like t A power law in the late time behavior of the H can be directly related to the large order behavior of the q-expansion coefficients c n . We find that c n ∼ n s with where s is the dominant power of the coefficient growth, and we are assuming that H(t) does grow at large t, which roughly requires h H > c 16 . Examples of this behavior are shown in figure 14. If H(t) decays at late times, then there must be cancellations in the sum over q n , and we cannot predict such a simple power-law. So in addition to directly computing the late time values of the Virasoro blocks, we can test whether the blocks follow the t − 3 2 decay simply by comparing the coefficients of the q-expansion of H to the prediction (3.8). This is actually a more efficient method that allows us to access certain regimes, such as larger c and h H of the parameter space where the direct Virasoro block calculation converges poorly.
However, the prediction (3.8) is less universal than the t − 3 2 behavior. For example, outside the regime where H(t) grows, the coefficients c n can have alternating signs, so that there are large cancellations between different terms in the q-expansion. Then the magnitude of the coefficients will no-longer follow the simple pattern depicted in figure 3.8. Empirically, another example is when h L c is small. In this case the coefficients are pretty small and show complicated irregular behaviors. Examples can be see in figure  23 in the appendix. Yet in all cases the overall late time behavior of the heavy-light Virasoro blocks is still the t − 3 2 power law. One would hope to derive the power-law behavior c n ∼ n s using the Zamolodchikov recursion relations. Unfortunately, it appears that this behavior arises from a large number of cancellations between much larger terms. Thus we leave this problem to future work.

Implications for Information Loss and the Bootstrap
In the semiclassical limit, heavy-light Virasoro blocks decay exponentially at late times. We do not expect that perturbative corrections in G N = 3 2c will alter this conclusion, and to first order this has been demonstrated explicitly [4]. Thus the late time powerlaw behavior of the exact blocks represents a non-perturbative correction that ameliorates information loss (insofar as information loss is tantamount to late-time decay). However, since the Virasoro blocks continue to decay, albeit much more slowly, this effect does not solve the information loss problem. For this we need an infinite sum over Virasoro blocks in the O L O L OPE channel. 5 Let us examine the correlator as a sum over blocks from the point of view of the bootstrap equation [46][47][48]. This equation dictates that 6 Here we have equated a sum over energies in the O H O L OPE channel with a sum over heavy-light Virasoro blocks in the O L O L OPE channel. In d > 2 dimensions this equation would be meaningless at large t, because we would be well outside the regime of convergence of the OPE expansion on the right-hand side. Remarkably, as discussed in section 2.1, the Virasoro block decomposition converges for all values of t, so it is possible to try to 'solve' for the coefficients P h,h by equating the large t behavior of both sides. More generally, one could take the limit |q| → 1 with various phases for q and derive new, potentially analytic regimes for the bootstrap (this is non-trivial because it could enable a partial analytic treatment without requiring a complete solution to the bootstrap equation). The only obvious obstruction to this procedure is that we do not have simple analytic formulas for the Virasoro blocks in such limits. As we have already noted, equation (3.9) can only be satisfied at late times if we have an infinite number of Virasoro blocks contributing on the right-hand side. Such infinite sums are compulsary in order to reproduce conventional OPE limits [23][24][25][26]. But it is easy to see that the Cardy formula and the asymptotic expectations on P h,h from Euclidean crossing or the light-cone OPE limit are insufficient to account for the late-time behavior. The reason is that conventional arguments require the large h,h terms in equation (3.9) to reproduce either the identity (vacuum) or perhaps the contribution of low dimension or low twist operators in the crossed channel. These would correspond to the very small E region of λ LH (E). But the late time behavior arises from the collective contributions of ∼ e S states with large E ∼ h H +h H , not from the small E states. 7 In this regard there is an amusing connection with Maldacena's original discussion [20] of the large time behavior. He suggested that in a black hole background, contributions from the vacuum, corresponding to the E = 0 term in equation (3.9), might resolve the information loss problem. But the vacuum in the O H O L OPE channel just corresponds with the Cardy-type growth (or more precisely OPE convergence [24] type growth) of P h,h . So this simple OPE convergence growth fails to account for the late time behavior for the same reason that Maldacena's suggestion did not resolve the information loss problem.
In summary, the late-time bootstrap equation (3.9) cannot be solved without providing a more refined asymptotic formula for P h,h at large h,h. However, it does not appear that a discrete spectrum in the O L O L channel is required to obtain the correct late-time behavior. We will not pursue this in detail since we only have some rough empirical information about the behavior of V h (t), but it might be interesting to study this bootstrap equation for the case of the partition function [42] where the Virasoro characters are known in closed form.

Some Philosophy
Eventually, we hope to learn about bulk reconstruction -and its limitations -by comparing exact CFT correlators to their semiclassical approximations. It is not clear whether this is possible, even in principle, due to ambiguities in the reconstruction process associated with bulk gauge redundancies (see e.g. [50] for a recent discussion). For now we will take a very instrumental approach, or in other words, we will try to 'shut up and calculate' some potentially interesting observables.
The information paradox pits local bulk effective field theory in the vicinity of a horizon against quantum mechanical unitarity. But in the strict semiclassical limit, information is lost and the (approximate) CFT correlators agree precisely with perturba-tive AdS field theory or string theory. Thus one would expect that bulk reconstruction should be possible in this approximation, since we have allowed the local bulk theory to 'win' the fight, at the expense of unitarity. 8 But even in the semiclassical limit, bulk reconstruction has been controversial [28][29][30][31]. On an intuitive level, this is because correlators at infinity must have exponential sensitivity to 'observe' physics near or behind a black hole horizon. At an instrumental level, this means that there may be obstructions to the existence of smearing functions mapping boundary to bulk observables. These issues can be avoided by going to momentum space [32,33], or perhaps via an appropriate analytic continuation [28] or cutoff procedure [31].
Another elementary issue with semiclassical bulk reconstruction is pictured in figure  2. The problem is that only the ingoing modes behind the horizon can be reconstructed in an obvious way from the degrees of freedom of a single CFT [28]. This can be understood by considering the extended AdS-Schwarzschild spacetime, or simply by studying Rindler space. Field theory degrees of freedom behind the horizon appear as a linear combination of modes from the left and right 'wedges', but in a single-sided black hole, only one asymptotic region is present.
If the goal is simply to compute correlators behind the horizon of a single-sided black hole, then there is a naive, instrumental way to obtain outgoing modes. One can obtain correlators that behave like those of the other asymptotic region by analytically continuing [28] CFT operators O(t, x) in Euclidean time toÕ(t, x) ≡ O t + iβ 2 , x . This procedure has an important flaw -operators on opposite sides of the black hole should commute, but O andÕ may not. Nevertheless, we can force O andÕ to commute (by definition) if we choose an appropriate but ad hoc analytic continuation procedure for correlators involving O andÕ. Conceptually, this does not seem to be an improvement on state-dependent mirror operators [32,33], which represent a modification of quantum mechanics. In fact, our procedure implements its own form of state-dependence, since the analytic continuations will depend on all of the other local operators inserted into the correlator. However, the prescription does have the simple advantage of being relatively precise and unambiguous.
In any case, we are led to a very simple question -do the correlators of operators like O t + iβ 2 , x receive large non-perturbative corrections? Do the semiclassical Virasoro blocks provide a good approximation to the exact blocks with these kinematics?

Forbidden Singularities and Thermofield Doubles
The questions raised in the previous section can be explored using the methods of this paper. They are also closely related to observations about information loss [7]. Finite-temperature correlation functions must satisfy the KMS condition, which for identical operators just means that O(t, x)O(0) β must be periodic in Euclidean time with period β. It has been shown that in the large central charge limit with h H > c 24 , heavy-light Virasoro blocks appear thermal. 9 Since the 4-point correlator has an OPE singularity as z → 1, in the heavy-light semiclassical limit, it will also have singularities at z n = 1 − e nβ for all integers n. While such singularities are permissible for correlators in the canonical ensemble, they are forbidden [7,27] from 4-point correlators of local operators in unitary CFTs. They are also forbidden from individual Virasoro blocks at finite central charge [5,7]. Thus exact Virasoro blocks completely disagree with their semiclassical counterparts at z n = 1 − e nβ , the locations of the singularities. So to summarize, we know that the exact and semiclassical blocks match at z = 0, and completely disagree at z = 1 − e nβ for n = 0. Thus it is natural to wonder whether the semiclassical blocks are a good approximation at z = 1 − e − β 2 −it , which corresponds to the location of O(t + iβ 2 ). More generally we would like to understand the kinematical regimes where the (leading) semiclassical approximation breaks down.
We observe from figure 15 that as expected, the exact Virasoro blocks do not have forbidden singularities. Nevertheless one might have expected to see bumps or local maxima at z n = 1 − e nβ , whereas the exact correlator simply grows as a function of z ∈ [0, 1). In fact local maxima are prohibited because the exact blocks are analytic functions of q and z away from the true OPE singularities. 10 Thus the semiclassical approximation breaks down badly beyond the first forbidden singularity.  Figure 15. In this plot, we compare the exact and semiclassical blocks. One can see that at the positions of the semiclassical forbidden singularities, the exact blocks are smooth. Fixing h L and h H c as we increase c, the exact blocks approach the semiclassical block in the region between the origin and the first forbidden singularity. However, beyond the first forbidden singularity the exact blocks deviate greatly as we increase c. This indicates that we have passed a Stokes line (emanating from the forbidden singularity) and some other semiclassical saddle dominates the exact blocks in the large c limit. The gray line is the position of t = iβ 2 .
We compare the exact and semiclassical blocks at finite time in figure 16. We see that the semiclassical blocks remain a good approximation to correlators of O(t + iβ 2 ) as long as we avoid the long-time region of t ∝ S that was discussed in section 3. In particular, there is not a significant difference between the quality of the semiclassical approximation to correlators of O(t + iβ 2 ) and O(t). The most naive interpretation of this fact is that non-perturbative quantum gravitational effects do not obstruct local physics across the horizon of pure, energy-eigenstate black holes. A qualitatively similar conclusion was reached for late-time deviations [52] from the semiclassical limit. This result was also anticipated by the analytic analysis of [7], which only suggested large non-perturbative corrections within 1 √ c of the forbidden singularites. In the next section we will discuss that analysis and compare it with our numerical results.

Fate of the Semiclassical Approximation from Analytics and Numerics
We do not have to rely entirely on numerics to explore the regime of validity of the semiclassical limit. It has been shown that the vacuum block's forbidden singularities have a universal resolution due to non-perturbative effects in central charge. Specifically, the heavy-light vacuum block (with h L and h H c held fixed at large c) should obey  figure 6. Apparently the semiclassical approximation works well at t + iβ 2 . an approximate differential equation [7] where τ = − log(1 − z) is a Euclidean time variable, and this equation neglects terms of order 1/c 2 and higher as well as effects that are less singular near the forbidden singularities. We provide the functions g H and Σ H in appendix B.1. This differential equation also predicts [7] that the semiclassical vacuum block will receive large nonperturbative corrections after a Lorentzian time of order S BH h L T H ∝ c h L . That prediction was corroborated in section 3.
Neglecting the term proportional to 1 c on the right-hand side, equation (4.2) is solved by the semiclassical heavy-light vacuum block. But when the right-hand side of this equation becomes large, non-perturbative effects come into play, resolving the forbidden singularities. We plot contours of the function |Σ H V V | for h L = 1 in figure  17. We see that this function becomes large and makes important contributions in the immediate vicinity of the forbidden singularities, though at sufficiently large c the right-hand side of equation (4.2) will remain small a finite distance away from these singularities. At a more detailed level, the function |Σ H V V | can be compared directly to the deviation of the numerical and semiclassical vacuum block. We plot contours of the ratio of the exact and semiclassical blocks in the ρ unit disk, corresponding to the entire Euclidean z-plane in figure 17 (recall that we compared various kinematic variables in figures 3 and 4).
Our numerical results demonstrate that the semiclassical approximation breaks down in a finite region enclosing the forbidden singularities. We believe this phenomenon occurs because Stokes and anti-Stokes lines (for review see e.g. [53]) emanate from the forbidden singularities, as has been demonstrated for the correlators of degenerate operators [7]. As we cross Stokes lines, the coefficients of semiclassical saddles change by discrete jumps. Across anti-Stokes lines saddles exchange dominance.
Near the OPE configuration z ∝ ρ ∝ q ≈ 0 where the light operators collide, a special 'original' semiclassical saddle dominates the large c limit [7] of the Virasoro blocks. But in a finite region near the forbidden singularities, different semiclassical saddles [5] can come to dominate, and the original saddle may become sub-leading. In other words, analytic continuation in the kinematic variables does not commute with the large c limit. Non-perturbative effects can dramatically alter the behavior of CFT 2 correlation functions with these kinematics, supplanting the naive semiclassical limit and the perturbation expansion around it.
It would be fascinating if the black hole interior depends in some way on the behavior of CFT correlation functions in these regimes. Note that when h H < c 24 , so that the heavy background state does not correspond to a black hole, the original semiclassical approximation remains good throughout the Euclidean region. We demonstrate this explicitly in figure 18. So the breakdown of the semiclassical limit exhibited in figure  17 really does depend on the presence of a black hole, and is not a general feature of all Virasoro blocks at large central charge.

Discussion
We would eventually like to resolve the black hole information paradox by doing the right calculation. In the context of AdS/CFT, this means discerning under what circumstances, if any, bulk reconstruction is possible near and behind black hole horizons.
If firewalls [51] are completely generic, or if bulk reconstruction is sufficiently ambiguous, then this could be a fools errand. But even in this case, one can still hope for a more constructive argument rather than various reductio ad absurdums [34]. For example, one would like to reconstruct the 'experience' of a collapsing spherical shell, and explicitly compute the timescale beyond which subsequent infallers will not see a smooth (or well-defined) geometry.
But let us imagine that the strict semiclassical limit is not misleading and black holes often have smooth interiors. In this case, violations of bulk locality should arise from the difference between computations in the semiclassical limit and the exact CFT observables (or perhaps meta-observables). This sort of approach has already been successfully pursued in the context of local bulk scattering [27]. We have identified gross differences between exact and semiclassical CFT correlators in both the late Lorentzian time and the Euclidean regime. These do not seem to affect a certain naive bulk reconstruction algorithm, but perhaps they do afflict more sophisticated methods yet to be developed. Hopefully we have done some of the right calculations but do not yet know how to give them the right interpretation. In the case of quantum mechanics and QFT, we were in that sort of boat for decades.

A Details of Recursion Relations and Our Algorithm
In this appendix we will present more details about Zamolodchikov's recursion relations and the algorithm we used to compute with them.

A.1 Zamolodchikov's Recursion Relations
There are actually two Zamolodchikov recursion relations, based on viewing the Virasoro blocks as either a sum over poles in the central charge c or the intermediate state dimension h. The latter is more powerful and will be our focus. The and the inverse transformations is If we parametrize the central charge c, the external operator dimensions h i and the degenerate operator dimensions h mn as follows then the function H (b, h i , h, q) can be calculated using the following recursion realtion where R m,n is given by and the ranges of p, q, k, and l are: p = −m + 1, −m + 3, · · · , m − 3, m − 1, q = −n + 1, −n + 3, · · · , n − 3, n − 1, k = −m + 1, −m + 2, · · · , m, l = −n + 1, −n + 2, · · · , n.
The prime on the product in the denominator means that (k, l) = (0, 0) and (k, l) = (m, n) are excluded. Note that our definition of λ p,q differs by a factor of − i 2 from the original paper.
In each iteration of the recursion relation A.6, the only thing that changes is the value of the intermediate state dimension h → h m,n + mn, which only depends on the values of m and n. For simplicity we'll omit the arguments and denote H(b, h i , h, q) as H and H(b, h m,n + mn, h i , q) as H m,n in the following discussion.
This recursion relation was derived by viewing the Virasoro block V h as a function of the intermediate dimension h, so it can be written as a remainder term that survives when h → ∞ plus a sum over poles at h = h m,n , where h m,n are the dimensions of the degenerate operators. The prefactor in front of H in A.1 is the h → ∞ limit of V h , as can be derived from [16][17][18]. The reason that V h has poles at h = h m,n is because of the existence of the null-operator (whose norm is zero) at level mn of the descendants of O hm,n , which usually will make V h diverge when h → h m,n . 11 The residue of the pole at h m,n will be proportional to the block V hm,n+mn with intermediate operator being the null-operator with dimension h m,n + mn. Thus, these residues will have high powers of q, which accounts for the q mn factor in front of H m,n and naturally makes the Virasoro block V h a series expansion in q.
The numerator of the factor R m,n is constructed such that it vanishes when O 1 (or O 3 ) belongs to the set of operators allowed by the fusion rule of O 2 O hm,n (or O 4 O hm,n ). The denominator of R m,n comes from the norm of the null-state when h → h m,n (factoring out h − h m,n ); as far as we know, although it has passed numerous checks, it's never been derived from first principles.

A.2 Algorithm
In this paper, we only consider the case that h 1 = h 2 = h L and h 3 = h 4 = h H . Under this circumstance, R m,n becomes directly proportional to λ 2 p,q , so R m,n = 0 whenever (m, n) are both odd, because (p, q) can then be (0, 0). This means that every H m,n with odd mn is also zero, as every term contributing to it contains at least one R m l ,n l with odd m l n l . As a consequence of this, only even powers of q ever appear, and there's no need to compute anything with odd mn. This provides some simplification for the calculation, but it's easy to generalize the following discussion to the case that all h i s are different. Now we turn to the algorithm we used to compute the recursion relation. The main idea is to sort every contribution to the functions H and H m,n by its order in q. By doing this from the beginning of the computation, we are able to use lower-level terms as partial sums for the higher-level terms, saving a great deal of computation.
Denote the coefficient of q k in any function f as f (k) . Then the recursion relation for the coefficients of q k in the function H is H (12) There are several other tricks that one can do to even speed up the calculation. For example, one can precompute all of the residue prefactors Rp,q hm,n+mn−hp,q ≡ C m,n,p,q in A.9. There are only O(N (log N ) 2 ) of these, so we can save time by computing them in advance and reusing them. Although precomputation dramatically improves performance, it also doubles memory consumption; but since we store the H (i) m,n in diagonals, this can be ameliorated by deleting them after they're used, as shown in 19.
Precomputing C m,n,p,q can only improve overall speed if each of its terms can be computed in constant time. This is potentially problematic, since R p,q contains two products of O(pq) complexity, but it can be solved by filling R p,q recursively -R p,q can be computed in O(p) time from R p,q−2 , and there are only O(N log N ) of them, so the computational complexity of filling all R p,q is just O (N 2 log N ). These can be further sped up by pairing up terms to rewrite all of the defining equations in terms of b 2 and λ 2 m,n instead of b and λ m,n . In addition to the reduced number of multiplications, this also allows the entire computation to be done using real numbers when c > 25, which is generally an order of magnitude faster. When c < 25, b 2 becomes complex, and even though the final coefficients must be real by unitarity, this only occurs at the very last step in the form of a b 2 ↔ 1 b 2 = (b 2 ) * symmetry. We have implemented this algorithm in both Mathematica and C++(with Mathematica integration). The Mathematica notebook is included as a companion to this paper, while the C++ implementation is maintained at https://github.com/chussong/ virasoro. The C++ implementation is about one order of magnitude faster, and the coefficients used in this paper were obtained using it. The C++ implementation has used the GMP [54], MPFR [55], MPC [56], and MPFR C++ [57] numerical libraries. On standard personal computers we were able to compute the H (k) to k = 1000 in around two minutes or k = 2000 in about 22 minutes (for c > 25 so that b is real); the main barrier to going higher is memory consumption, which grows roughly as N 3 log N : we need to remember O(N 2 log N ) numbers and they need to be kept at O(N ) bits of precision due to the increasingly large cancellations between different H m,n , which often reach into the thousands of binary orders of magnitude.
Using a cluster with 128 GB of RAM, we estimate that we could reach order of 6000 in a few hours. We also find that the coefficients of q i approach a power law in i well before the limits of our desktop computation, and expect that a numerical fit for this power law would be good enough to get higher order coefficients.
At the end of this section, we want to mention an issue about the recursion relation if b 2 is a rational number. Notice that the denominator in A.9 and the denominator of R m,n in A.7 can be zero: Both of these will eventually occur for any rational choice of b 2 . This would appear to preclude numerical computation entirely (since for numerical calculation, b provided to the computer will always be rational), but actually for almost all rational numbers they will not appear until very high orders in the computation, so they can be ignored as long as the numerator or denominator of b 2 (as a irreducible fraction) is very large. In this paper, we've chose √ c to be irrational (and set b to be a very high-precision number) to avoid this problem.

B.1 A Non-Perturbative Differential Equation for the Vacuum Block
Here we describe the functions appearing in the differential equation (4.2). Note that although the equation itself is perturbative, its solution includes non-perturbative corrections to the heavy-light vacuum Virasoro block. The equation was derived [7] by studying the general differential equations satisfied by degenerate operators and then analytically continuing these equations in the integer index r labeling the degenerate operators. We should also note that although equation (4.2) only includes some of the first 1/c corrections, if one zooms in on the vicinity of the forbidden singularities by holding √ c(z − z n ) fixed at large c, then the equation incorporates all of the leading effects at large c. As discussed in [7], there are both general arguments and consistency checks on the validity of equation (4.2).
We identify the parameter r = 2πiT H = 1 − 24h H c , so that T H is the Hawking temperature associated with the heavy operator. We also are using a Euclidean time variable τ = − log(1−z). Then the functions included in equation (4.2) are g H ≡ g 2πiT H with g r (τ ) = coth τ 2 − r coth For derivations and more complete descriptions see [7].  figure 10, but includes a match to the semiclassical blocks obtained using the methods of [5], which allow for h, h L ∝ c. The poorly fitted dashed line is the approximation of equation (2.10), which assumes h, h L c, and clearly provides a much less reliable fit for these parameter values.

B.2 Some Extra Plots
In this section we have included some extra plots for readers who might like to some more details and examples. These include the semiclassical fit to our numerical results for h, h L ∝ c using [5]  We also show some plots of the more complicated coefficient behavior which was alluded to in section 3.2, with the sign of the coefficients corresponding to the color of plotted points. Figure 23 illustrates a very common scenario where the coefficients are chaotic at low c, but as c increases they coalesce into distinct positive and negative lines. A spike-shaped feature then appears at low order and moves upward, turning the coefficients that it passes positive, until all (visible) coefficients have become positive. The two lines then gradually merge into a single power law similar to those shown in figure 14.