A thermal product formula

We show that holographic thermal two-sided two-point correlators take the form of a product over quasi-normal modes (QNMs). Due to this fact, the two-point function admits a natural dispersive representation with a positive discontinuity at the location of QNMs. We explore the general constraints on the structure of QNMs that follow from the operator product expansion, the presence of the singularity inside the black hole, and the hydrodynamic expansion of the correlator. We illustrate these constraints through concrete examples. We suggest that the product formula for thermal correlators may hold for more general large N chaotic systems, and we check this hypothesis in several models.


Introduction
A thermal two-point function describes the response of a system at finite temperature to a small perturbation [1]. In theories with simple holographic duals, it is captured by the wave equation in a black hole background [2][3][4][5][6].
In this paper we mostly (but not only!) concern ourselves with holographic CFTs, for which the relevant geometry is a black hole in AdS. For this case, the problem of computing the thermal response is reduced to a quantum-mechanical scattering problem in a potential V (z) which encodes the black hole geometry. Very generally for non-extremal black holes, this potential has the following features, 1 V (z) ∼ 1 z 2 z → 0 (AdS boundary), ∞ n=1 a n e − 4πn β z z → ∞ (BH horizon). (1.1) The exponential decay close to the horizon is consistent with the periodicity of the potential V (z) = V (z + i β 2 ) under complex shifts of z. Let us denote the two-sided thermal function that can be obtained by solving the wave equation as G 12 (ω), where we keep the angular/spatial momentum dependence implicit. There are two basic facts about G 12 (ω) that will be important for us in the present paper: • It is a meromorphic function of ω [7][8][9]. Its singularities, known as quasi-normal modes (QNMs), encode the characteristic decay of perturbations in time [10,11].
• It has no zeros, or in other words, 1 G 12 (ω) is an entire function of ω. This property is intimately related to the behavior of the potential V (z) close to the black hole horizon, see [8] and Appendix C.
Even though these properties have been extracted by studying classical wave propagation on a black hole background, at the level of the thermal two-point function they can be considered more generally in situations where there is no well-defined notion of a classical geometry or black hole horizon. In fact, extrapolating beyond holography, we will find that both properties hold at finite coupling in several low-dimensional examples of large N systems.
The two properties above imply, modulo simple technical details deferred to the main text, that we can write a product formula representation for the two-sided correlator, 2 , G 12 (0) > 0. (1.2) Another characteristic feature of the black hole geometry in AdS d+1>3 /CF T d>2 is the presence of the curvature singularity. In [9] it was argued that the singularity leads to exponential decay of the two-sided two-point function at large imaginary ω. An illuminating way to encode this property is to say that the two-sided correlator obeys the following dispersive sum rules BH singularity : (1.4) see Figure 1. Again, even though (1.4) was extracted by studying the classical wave equation, it is interesting to contemplate the possibility that it holds more generally and consider its violation as some sort of "resolution" of the singularity. The purpose of this paper is to analyze the general properties of the product formula (2.9), as well as to explore its workings in various examples. In particular, we study the constraints imposed on the structure of QNMs by the OPE; we discuss how the singularity sum rules (1.4) are satisfied and violated; and finally, we explore the relation between the structure of QNMs and the low-energy (hydrodynamic) expansion of the correlator. In the latter case, we find that the hydrodynamic expansion encodes certain moments of the QNM density.
After examining the general constraints, we specialize to several analytically solvable examples, including BTZ and Rindler space. We then proceed to discuss several examples in pure GR and beyond for which no analytic expression is available. To compute the QNMs for these cases, we use the package QNMSpectral [12].
In the final section of the paper, we show that the two-point functions in several SYK-type models satisfy the holographic properties mentioned above. We conclude by presenting a hypothesis regarding the two-point function in planar chaotic theories, and with a discussion of various open directions.

Analytic structure of thermal two-point functions
In this paper we are mainly interested in the two-sided thermal two-point function, defined as follows, where β is the inverse temperature and the spatial coordinate ⃗ x labels a point on the spatial manifold on which the CFT lives. We will mostly consider the case where ⃗ x ∈ R d−1 in the present paper.
We will be considering G 12 (t, ⃗ x) in Fourier space where k ≡ | ⃗ k|. We will be mostly concerned with the dependence of the two-sided correlator on ω, while keeping k fixed. We will therefore sometimes write G 12 (ω) with the spatial dependence kept implicit. We have chosen to work with G 12 (ω, k) because of its nice properties, but all other two-point functions can be easily obtained once G 12 (ω, k) is known, see Appendix A.

Definitions and general properties
Let us review some basic properties of the two-sided thermal two-point function G 12 (ω) as a function of ω: • Unitarity: on the real axis the two-sided correlator is real and non-negative This property simply follows from the insertion of a complete set of states between the two operators in the definition of the two-sided correlator and unitarity of the underlying theory. 3 Assuming that G 12 (ω) is real-analytic for real ω, we get in the complex plane G * 12 (ω) = G 12 (ω * ). (2.4) When G 12 (ω) = 0 on the real axis something special happens, namely the local operator with a given energy ω annihilates the thermal state. We expect that this is only possible in free (or integrable) theories, and that in interacting theories G 12 (ω) > 0 on the real axis.
(2.5) 3 In the most general setting, G 12 (ω) is a distribution, and therefore (2.3) should be understood in the distributional sense. In the present paper we will treat G 12 (ω) as a function. In particular, we will study its analytic properties in the complex ω plane. We believe that this is justified for interacting CFTs on S 1 × R d−1 and large N interacting CFTs on S 1 × S d−1 above the Hawking-Page phase transition [13]. It is true by inspection for holographic CFTs.
• OPE limit: the large frequency limit of the correlator is universal and is controlled by the unit operator, lim ω→+∞ G 12 (ω) ∼ e − βω 2 ω 2∆−d , (2.6) where ∆ is the conformal dimension of O and the proportionality coefficient is known exactly in terms of d and ∆. As we will see later, corrections to this formula at large ω can be systematically computed in the 1/ω expansion, and are determined by the OPE coefficients and thermal one-point functions of operators that appear in the O × O OPE expansion.
The properties above hold in any CFT and do not rely on holography.

Properties of holographic thermal correlators
Next we list the properties of G 12 (ω) which are specific for theories with a classical gravity dual: • Meromorphy: the only singularities of G 12 (ω) are isolated simple poles in the complex plane [7][8][9]. 4 These singularities correspond to quasi-normal modes (QNMs), and they characterize the real-time decay of a perturbation to the thermal plasma [10,11]. By virtue of (2.4) and (2.5), QNMs come in families (ω n , −ω n , ω * n , −ω * n ). Similarly, the residues of G 12 (ω) for various frequencies within one family are all related to each other, Res ωn G 12 (ω) = −Res −ωn G 12 (ω) = (Res ω * n G 12 (ω)) * = −(Res −ω * n G 12 (ω)) * . (2.7) • No zeros: the two-sided correlator does not have zeros in the complex ω-plane [8]. This property follows from the fact that the two-sided correlator is obtained by solving the wave equation in the black hole background and from the universal properties of the scattering potential close to the black hole horizon, see Appendix C. Put differently, 1/G 12 (ω) is an entire function.
The three properties above imply the following product formula for holographic thermal correlators 5 , . (2.9) To see this, we use Hadamard's factorization theorem, see e.g. [15], which states that an entire function f (z) of order m and zeroes a n can be written as ∞ n=1 E ⌊m⌋ (z/a n ), (2.10) where P (z) is a polynomial of degree q ≤ m and ⌊m⌋ k=1 e z k /k (2.11) are elementary factors. By assumption, 1/G 12 is of order m = 1. Imposing evenness, ℓ = 0, and noting that E 1 (z)E 1 (−z) = E 0 (z)E 0 (−z) = 1 − z 2 , the product formula follows.
In writing (2.9) we assumed that ω * n ̸ = −ω n , which is not true for purely imaginary QNMs. It is trivial to extend the ansatz by including the product of purely imaginary simple modes ω n = iω n as follows, ∞ n=1 1 + ω 2 ω 2 n −1 .

Black hole singularity sum rules
In [8,9] it was shown that for the Schwarzschild black hole, G 12 (ω) at large imaginary ω is controlled by a real geodesic that bounces off the black hole singularity [16,17]. As a result, the correlator decays exponentially as ω → ±i∞ [8,9], (2.12) The two-point function is also exponentially decaying in the region ω → ±∞, where it is controlled by the OPE behavior (2.6). Therefore in the presence of the curvature singularity there exists a positive number κ such that (avoiding the poles) This is a stronger condition than (2.8), which allows for a vanishing decay rate λ(θ) = 0 along some complex ray.

Figure 1:
The contour integral along C∞ in (2.14) vanishes due to the asymptotic decay of the correlator. By deforming the integration contour inside and picking up the contribution of QNMs we get the sum rules (2.15).
We can use the condition (2.13) to write down a simple set of sum rules, 6 BH singularity : (2.14) see Figure 1. The asymptotic behavior (2.13) guarantees that the contour integrals (2.14) vanish. We can then deform the contour and rewrite these sum rules in terms of the residues of G 12 (ω), λ n , at the QNMs ω = ω n . For m even the sum vanishes identically, but for m odd we find a nontrivial constraint: Re n ω m n λ n = 0, for odd m > 0, (2.15) where we used the relations (2.7) between residues of poles in the same family. Note that the sum rules (2.15) immediately imply that the residues λ n decay exponentially fast in n. As we will see in Section 5.3, the singularity sum rules follow directly from the product formula (2.9) with some mild assumptions on the asymptotic structure of QNMs. In Section 5.3 we will use the exponential decay at ω → i∞ along with the OPE to uniquely fix the leading high-energy asymptotics of QNMs.
Let us also notice that we can consider finite energy singularity sum rules by placing the integration contour C Λ in (2.14) at some finite radius |ω| = Λ. 7 In this case the sum rules are satisfied up to exponentially small corrections e −cΛ . The advantage of considering sum rules at finite ω is that we expect them to hold also at finite 't Hooft coupling (or string length) as long as Λ ≲ 1 ls .

QNMs and black hole geometry
In the next section we will analyze various constraints placed on the QNMs by the product formula (2.9). To set the stage for these constraints, let us first give a brief overview of QNMs and their connection to the geometry of the black hole background, closely following [8,21]. For a detailed review, see [11]. We consider a spherically symmetric asymptotically AdS d+1 black hole, with metric In order to write the scalar field wave equation (□ − ∆(∆ − d))ϕ = 0 in a convenient form, let us Fourier expand where Y J ⃗ m are spherical harmonics on S d−1 . In terms of ψ, the wave equation takes the form of the Schrodinger equation for a particle in a potential, where dz = − dr f (r) is the tortoise coordinate. The horizon is located at z = ∞ and the boundary at z = 0. The potential is given by where ν = ∆ − d/2. The potential for a black brane takes the same form, with J(J + d − 2) replaced with k 2 . Now let us discuss the behavior of the potential (3.4) near the horizon and infinity. For a general AdS black hole, the asymptotic behavior is universally given by The behavior at the boundary follows from the asymptotics f (r) ∼ r 2 at large r. The structure of the potential near the horizon can be justified as follows. Since the redshift factor f (r) has a simple zero at the horizon r S , we can expand The tortoise coordinate is then Solving for r(z), we find (3.9) where the higher corrections are of the form exp (−4πn(z − z 0 )/β). Since the potential is an analytic function of r, we therefore find that the potential takes the form shown in (3.5).
Because perturbations can fall into the horizon, we cannot define normal modes of the scalar field. Instead, we consider quasi-normal modes, which are ingoing at the horizon and normalizable at the boundary, Such resonances can only exist at a discrete set of complex energies, which satisfy Im ω < 0 for real J [10]. Note that Im ω < 0 implies that perturbations decay in time, whereas poles in the upper half plane would lead to a growing mode, signifying an instability of the black hole. There are also constraints on QNMs from causality [22][23][24], but we do not explore them in this paper. In order to obtain a qualitative understanding of the structure of QNMs, it is useful to consider the WKB limit ∆ → ∞ with ω/∆ and J/∆ fixed. The potential (3.4) then takes the form For z 0 extremizing the potential, there exists a line of QNM poles emanating from V (z 0 ). However, not all extrema contribute. To find the relevant values of z 0 , one (a) (b) Figure 2: (a) A potential with a metastable minimum outside the horizon, leading to a line of weakly damped QNMs. The potential goes to +∞ at the boundary, to −∞ at the singularity, and is zero at the horizon r = r S . (b) A potential with a global stable minimum, corresponding to a line of virtual bound states. The minimum is in between the inner horizon r − and the outer horizon r + .
must first solve the equation V (r) = u 2 for the turning point r(u). The physical turning point is defined to be the solution r(u) which behaves as r(u) ∼ u/∆ for large real u.
Then the extrema of the potential corresponding to QNMs are those at which the physical turning point merges with another turning point at some value of u in the lower half plane. We refer the reader to [21] for further details. The first few modes of this line of poles are given by the expression The branches of the square roots should be chosen so that ω n is in the lower half plane. Let us now discuss several important cases of the Bohr-Sommerfeld approximation to QNMs.
• Weakly damped modes: The first example we consider is where V (z) has a metastable minimum outside the horizon, as depicted in Figure 2a. Then in the approximation (3.13), the associated QNMs are purely real. The imaginary part is exponentially small since it is related to tunneling over the potential barrier, and it is not captured by the leading WKB approximation. Since these modes decay very slowly, they represent the leading contribution to the two-point function at late times. One case where such a metastable minimum arises is at large J. The corresponding modes are related to stable orbits around the black hole [8,21,[25][26][27].
• Virtual bound states: In the previous example, the minimum of the potential was outside the horizon. We can also consider the case where the minimum is behind the horizon in a region where V (z 0 ) < 0, see Figure 2b. We see from (3.13) that the modes are then purely imaginary. Charged black holes provide an example of this phenomenon, as we will see later. In this case, there is a stable minimum of the potential between the inner and outer horizon.
• Complex extrema: Generically, extrema of the potential are at some arbitrary point in the complex plane off the real r axis. In this case there is no clean physical interpretation of the QNMs in terms of the real section of the geometry. Note that the extrema come in complex conjugate pairs, since V * (r) = V (r * ). Therefore a complex extremum leads to two lines of poles related by ω n → −ω * n , at an angle determined by (3.13). For example, consider the QNMs of the black brane metric f (r) = r 2 − 1/r d−2 at zero spatial momentum. There are d complex extrema at r d = 1 − d/2, but the only two that contribute according to the prescription described above have phase e ±iπ/d . This leads to two lines of poles in the lower half plane at angles e −πi/d and e −(d−1)πi/d [21].
Let us now comment on the QNMs beyond the lowest lying modes (3.13). In the case of a metastable minimum, the low-lying QNMs have exponentially small imaginary part, but at some n max the modes go off into the complex plane at a finite angle. Therefore there are only finitely many weakly damped modes. In contrast, the virtual bound states and modes associated to complex extrema go on forever. As we will see later, there must be infinitely many QNMs in order to reproduce the operator product expansion.
In this section we have discussed the QNMs in the large ∆ limit. The highly damped QNMs can also be computed for order one ∆ by solving the wave equation near the singularity and near the boundary and then matching the two solutions. See [28,29] for details.

Asymptotic Minkowski OPE
Let us consider a CFT two-point function of scalar primary operators on S 1 β × R d−1 . We can write the following OPE expansion, (cos θ) are the usual Gegenbauer polynomials and τ is the Euclidean time β > τ ≥ 0. The sum goes over the primary operators that appear in the OPE of O ×O, and the expansion coefficients a ∆,J are given by the product of the corresponding threepoint function and thermal expectation value ⟨O ∆,J ⟩ β , see [30] for details.
In this paper we are interested in the Lorentzian correlator. More precisely, we consider the two-sided correlator, which is related to the Euclidean correlator above as follows, We would like to understand how the OPE expansion (4.1) constrains the form of G 12 (ω, k). We will be interested in the limit ωβ ≫ 1 with kept fixed. It has been suggested in [31] that in this limit the two-sided correlator admits the following asymptotic Minkowski OPE expansion where the conformal blocks G ∆,J (ζ) take the following form (see Appendix D for the derivation), The support of G 12 (ω, k) in (4.4), namely the presence of θ(1−ζ), was discussed in [32]. This expansion is supported by the known holographic and perturbative examples and we will assume it to be true in this paper. Note that G 12 (ω, k) is nonzero for k > ω, but it is expected to be exponentially suppressed in this region [6,33,34]. The basic idea behind (4.4) is the following. The coordinate space OPE leads to the large frequency expansion of the Euclidean correlator. The behavior of the Euclidean correlator at large Matsubara frequencies controls the behavior of the retarded correlator via the relation G E (ω n ) = G R (iω n ). Consider next the subtracted dispersion relation for the retarded correlator. The large frequency behavior at imaginary frequencies controls the large behavior of the spectral density or, equivalently, the two-sided correlator via complex Tauberian theorems for the Stieltjes transform [19].
Assuming that the G 12 (ω, k) itself admits a power-like expansion leads to the formulas above. More generally, rigorous formulas can be derived as in [19]. This is particularly relevant for finite c T CFTs on a sphere, where G 12 (ω) is given by a sum of δ-functions.
We next discuss the universal contribution to the OPE as well as its structure in holographic theories.

The unit operator and the stress tensor
The leading contribution at large frequencies comes from the unit operator ∆ = J = 0. The first subleading correction comes from the stress tensor. The leading terms in the expansion of the two-point function thus take the following form where the coefficient of the stress tensor is fixed by the Ward identities as follows, Here is the free boson two-point function of the stress tensor. Finally, f is related to the free energy density as follows For holographic theories a convenient formula was found in [35] where we used that f = − s d , and s is the entropy density. Plugging this into (4.7) we get (4.10) We have tested this formula in the case of an AdS black brane by numerically solving the wave equation at large ω. Let us also notice that for ∆ O = d, d − 1, ..., the contribution of the stress tensor drops out.

Multi-stress tensor operators
Next we discuss the contribution of multi-stress tensor operators T n to the OPE. These have scaling dimension ∆ T n = nd and spin 0 ≤ J ≤ 2n. In a generic CFT it is not known how to compute the contribution of these operators to the OPE. For theories with gravity duals it can be done, see [36]. Indeed, their contribution is controlled by the physics close to the AdS boundary. To keep the discussion simple, let us present the result for the contribution of T 2 operators in d = 4. Combining the results above with the ones in [36] we get 8 Notice that the effect that we observed at the level of a single stress tensor, namely vanishing of its contribution for ∆ O = 3 and ∆ O = 4, continues here as well. We see that for ∆ O = 3, 4, 5 the expression above vanishes. We believe the same effect continues for T n , and therefore we can write simple closed form expressions (up to nonperturbative corrections) for the OPE expansion. It follows that (4.6) represents the full answer for ∆ O = 3, 4, 5 in d = 4.

Leading nonperturbative correction
We now consider the leading nonperturbative correction to the large frequency expansion above. In coordinate space this correction has a regular expansion at small separations and therefore we expect that it could be related to the double trace operators. These are operators schematically of the form O□ n ∂ µ 1 ∂ µ 2 . . . ∂ µ J O. In the limit we are working their scaling dimension is 2∆ O + J + 2n. Therefore, we see from (4.1) that they will contribute to the OPE with analytic terms.
For theories with gravity duals these corrections can be computed by analyzing the wave equation in the large ω regime, see [8] and references therein. The solution to the wave equation close to both the singularity and the boundary can be approximated by Bessel functions in the tortoise coordinate z. For large ω the regions of validity of these approximations overlap and it is possible to obtain a solution valid for all z. Finally, one can use such a solution to extract the leading order contribution to (4.4) with its nonperturbative corrections (but no power law corrections). The leading term with the first nonperturbative correction reads (4.12) In Appendix E we show that this correction can be reproduced from the exact expression for G 12 (ω, 0) found in [37].

QNMs and the OPE
We have argued that the quasi-normal modes determine the holographic Wightman function up to a constant. By providing additional input for G 12 , we can place constraints on the spectrum of QNMs.
In this section, we will input the constraint from the operator product expansion at large real frequency. As ω → +∞, we should recover the zero temperature answer up to an overall factor of e −βω/2 , see [32,38,39], where the proportionality constant and subleading in 1 ω corrections can be found in the previous section.
The goal is to use this constraint to derive conditions on the asymptotic behavior of QNMs. The first step is to convert the product over QNMs into a convenient sum. From (2.9), we have On the other hand, from (5.1) we have Let us now compare these two expressions. It is immediately clear that (5.3) implies that the number of QNMs has to be infinite and that they must extend all the way to infinity. Indeed, imagine that there are a finite number of QNMs. Then ∂ ω log G 12 (ω) ∼ 1 ω at large ω, which is not consistent with (5.3). Next we consider several simple ways in which QNMs can approach infinity and see how this approach is constrained by the OPE.

A single line of asymptotic QNMs
We start with the simplest case, where the QNMs are asymptotically organized into a single line at angle θ, see Figure 3. We consider the following ansatz 9 Computing the large ω behavior of the sum we get where the leading asymptotic does not depend on the lower limit in the sum.
For a line of simple poles on the imaginary axis, one should divide the right hand side of (5.7) by two, since there is only one image. Let us now consider the first subleading term at large n. For this purpose, we take the improved ansatz ω n = re iθ n + se iϕ + . . . .
(5.8) 9 Strictly speaking these are the images of QNMs under KMS symmetry, not QNMs themselves.
Using the Euler-Maclaurin formula, the sum (5.2) to order 1/ω is Comparing with (5.3), we find To shed more light on this result, let us imagine that we add an extra QNM at ω 0 . A single QNM contributes at large ω as ∂ ω log G 12 (ω) ∼ − 4 ω . For the computation above adding a single QNM corresponds to effectively relabeling n → n − 1, which shifts the subleading constant as se iϕ → se iϕ − re iθ . Plugging this expression into (5.9) indeed leads to the expected shift − 4 ω . In this sense the subleading term in (5.8) is sensitive to the global structure of QNMs, but it is not sensitive to the details of their distributions at low energies.
To summarize, we have learned that for a single line of QNMs, the spectrum must be asymptotically linearly spaced. Moreover, the spacing and angle must be related in terms of the temperature, and the first subleading term in the expansion is related to the conformal dimension of the operator under consideration. These conditions apply for any black hole solution in which there is only one line of asymptotic QNMs. For example, the QNMs of massless fields with ∆ = d in large AdS d+1 /Schwarzschild black holes asymptotically approach the line 10 [28,29] which is easily seen to satisfy (5.7) and (5.10). Notice that the OPE expansion effectively maps to the 1 n expansion of the QNMs. Moreover, we see that at each order in the expansion we have two real parameters to fix, and we have only one equation.
To proceed to subleading orders we can write If γ is an integer, then it is straightforward to check that subleading terms in ∂ ω log G 12 (ω) receive contributions from all n in the sum, so we cannot fix t purely in terms of the asymptotics of the QNMs. Let us now consider the case where γ is fractional. At large ω, the sum (5.2) is dominated by n ∼ ω, so we can expand the summand in n −γ . Expanding ∂ ω log G 12 at large ω and plugging in (5.10) and (5.7) gives (5.13) Since fractional powers of ω do not appear in the OPE expansion of ∂ ω log G 12 (ω) for holographic correlators, it must be the case that In Section 6, we will check this relation for AdS black branes. We have seen that fractional powers can appear in the expansion of the QNMs at subleading orders. In fact, if there are two asymptotic lines of QNMs, fractional powers can appear in the leading behavior of one of the lines. We turn to this case next.

Adding another asymptotic line
We now consider the case of several asymptotic lines of poles. If all the lines are asymptotically linearly spaced, then the previous calculation generalizes immediately, and we find As mentioned above, a line of simple imaginary poles needs to be treated separately, and contributes 2π/r i to this sum. Similarly, we have the subleading constraint The more interesting situation is when not all lines are asymptotically linearly spaced. Let us consider the case of two asymptotic lines. We take one of the lines along the imaginary axis for simplicity, and choose the ansatz Proceeding as above, we find We now match powers of ω to the OPE expansion (5.3). The first term in (5.19) reproduces the result (5.7) from before. The other terms must cancel, giving the conditions γ = 1/α and . (5.20) We will confirm this prediction later by analyzing the QNMs in the presence of a higher derivative correction W 2 ϕ 2 , where W is the Weyl tensor. This higher derivative term introduces a new line of QNMs on the imaginary axis with asymptotic n 3 scaling, and we will see numerically that (5.20) is satisfied to high accuracy in this case.

Behavior in the complex ω plane
So far, we have analyzed the behavior of the correlator at large real ω. Let us now discuss what happens when ω is taken to infinity along some complex ray. In the case of a single line of QNMs, the lines of poles divide the complex plane into four asymptotic regions. For example, as ω → i∞ with the ansatz (5.4) with α = 1, we find It follows that the correlator exponentially decays at ω = i∞ with rate 2π cos θ/r. As long as θ ̸ = π/2, this decay rate is nonzero. Similar remarks apply for the other three regions of the complex plane, so the singularity sum rules (2.15) are satisfied when θ ̸ = π/2. We conclude that meromorphy and no zeroes imply the singularity sum rules, given a line of asymptotic QNMs that is not parallel to the imaginary axis.
In fact, the same conclusion holds if there are several lines of asymptotic QNMs, as long as at least one line is at angle θ ̸ = π/2. The lines divide the complex ω plane into different asymptotic regions, and in each region the correlator exponentially decays. Therefore the singularity sum rules will be satisfied.
It is instructive to compare the large frequency expansions along the real and imaginary axis. To this end we consider the ratio where we wrote only the leading contribution to the ratio and we definedβ through the leading asymptotic of log G 12 (iω) ≃ −β ω 2 . Bothβ and C 0 are sensitive only to the large n tails of the large QNM expansion, and therefore can be safely computed using the same methods that we used above.
Using ω n = re iθ n + se iϕ we get The holographic computation in the AdS-Schwarzschild background gives [8,9] β = β cot π d , which leads to the following expression for the asymptotic behavior of the QNMs, For ∆ = 4 this agrees with (5.11).

OPE sum rules
In the previous section we explored the basic relationship between the leading terms in the OPE expansion of the correlator and the structure of the QNMs. In this section we show how inputting more data about the asymptotic expansion for the quasi-normal modes leads to exact sum rules for QNMs. These sum rules encode the subleading terms in the OPE expansion, going beyond the previous section.

OPE sum rules for QNMs
Let us imagine that we have worked out the asymptotic expansion of the QNMs as above to a certain order in 1 n , such that ω n = ω asy n + δω n , We can then use the OPE to write down a set of sum rules for δω n . More precisely, we write where G asy 12 (ω) is obtained by replacing ω n with ω asy n in the product formula (2.9) without the prefactor G 12 (0). Given the OPE and ω asy n we can compute the large ω expansion of the left hand side. Expanding the right hand side at large ω under the sum we get Such an expansion only makes sense if the sums over n converge. The most dangerous contribution comes from δω n (ω asy n ) 2k−1 ∼ δω n n 2k−1 . For this to converge we require For example, to write down the sum rule corresponding to the contribution of the stress tensor in d = 4 we need k = 2, which requires k * = 4. Let us write more explicitly the first sum rule. We set where c is computed in Appendix F. Note that the phase of the last term in this expression agrees with the prediction (5.14). Using (6.6), we can write down the first nontrivial sum rule for k = 1. The left hand side of (6.4) takes the form where δc 3 is given explicitly in (F.5). Therefore we see that by combining the asymptotic large n expansion of the QNMs with the asymptotic OPE expansion we get exact nonperturbative QNM sum rules, which are sensitive to all QNMs.
Let us consider the k = 1 sum rule explicitly for ∆ = 4 and d = 4. Setting β = 1, we get the following equation . For m → ∞, the partial sums s m approach 0. We can compute the first few QNMs numerically using QNMSpectral and check how well the sum rules are satisfied. The reader can find these modes in the accompanying file modes.txt. We plot s m in Figure  5.
We also checked numerically that the leading large n asymptotic of δω n in this case takes the form ≃ 0.340 e −i π 48 n 7/3 and this tail correctly accounts for the ∼ 23% of the sum rule left for the k > 60 modes in Figure 5.

Overall constant
There is one last piece of data that comes from the OPE that we have not used. It is the overall normalization of the correlator at large ω, as given by (4.6). Similarly, on the product formula side (2.9), G 12 (0) has not entered in any of the equations above because it does not contribute to ∂ ω log G 12 (ω). Likewise, it canceled in the ratio G 12 (ω)/G 12 (iω).
As above, we can imagine that the asymptotic large n expansion of QNMs is known. We then can write . (6.9) We can now take the large ω limit to get the following equation where we set β = 1. We would now like to understand how this result is reproduced from the sum rule (6.10). The asymptotic form of the QNMs takes the form see (6.6) for d = 4. Computing the large ω limit and using (4.6) we get Therefore we get the following sum rule , so that lim m→∞ s m = 0. We compute ω n numerically using QNMSpectral, and we plot s m in Figure 6.
In principle we can repeat the same computation as above by including the next term in the large n expansion, δω asy n = 6π With this correction the sum rule becomes ∞ n=1 log |ω asy n + δω asy n | 2 |ω n | 2 ≃ −0.03873.

Positive moments for hydrodynamics
It is interesting to see how the location of QNMs translates into the low-energy (hydrodynamic) expansion of the correlator. In this section, we find that by restricting the location of QNMs to a subregion of the complex plane, we get nontrivial bounds on the hydrodynamic expansion coefficients.
The basic observation is that the no-zero property of the thermal correlator can be thought as positivity of the discontinuity of log G 12 (ω). The product formula (2.9) corresponds to the dispersive representation of log G 12 (ω), with log G 12 (0) being the subtraction constant. We then find that the coefficients of the low-energy expansion of the two-sided correlator are related to certain moments of the non-negative density of QNMs.
The mathematical structure that emerges is similar to bounds that appear in the context of scattering amplitudes [41]. The role of Wilson coefficients is played by the hydrodynamic expansion, whereas the role of dispersion relations is played by the product formula (2.9).
To follow this idea, it is convenient to introduce a positive-definite measure where ω n are QNMs which satisfy Re ω n > 0, Im ω n > 0. As before it is also convenient to separately consider purely imaginary QNMs, for which ω = iω n for realω n . In terms of these densities we can write Let us consider next the low energy expansion of the two-sided correlator. It takes the following form, where We immediately see that the hydrodynamic expansion is naturally organized in powers of the smallest QNM. To this end, let us introduce ω min andω min for the QNMs closest to the origin. It is convenient to define the moments in units of |ω min |. For this purpose we can write where the moments have been defined as follows, We therefore see that the low-energy expansion is related to moments of the QNM density. In fact, the moment problem forμ k is nothing but the simplest Hausdorff moment problem, which can be seen by switching toỹ = 1 ω 2 ∈ [0, 1]. The moment problem for µ k is more nontrivial since it involves the geometry of the 2d plane. Below we restrict our considerations to the cases where the QNMs are supported in the region Re ω n ≤ Im ω n . For this configuration µ 1 ≥ 0 and it is natural to consider bounds on the ratios ( µ 2 µ 1 , µ 3 µ 1 ). We have derived the bounds in two steps. First, we derived the bounds numerically by discretizing the region of support of QNMs and then solving the coupling maximization problem using LinearProgramming in Mathematica. We then have identified the boundary of the allowed regions analytically, finding that for the cases considered here they are very simple.

Purely imaginary QNMs
We first consider the simplest case when QNMs are purely imaginary so that µ i = 0. Restricting our attention to (μ 1 ,μ 2 ,μ 3 ), we can write down a set of optimal constraints [42]:μ The allowed region is depicted in Figure 7. An example of this type is provided by the BTZ correlator at zero momentum, for which we havẽ Doing a shift ∆ → ∆ − d−2 2 one recovers the Rindler space correlator, and by replacing ∆ → 2 q we get the result in the SYK model.

QNMs on a line
For the next model we consider the case where all QNMs are located on a single line in the complex plane. In other words, we consider ω n = e iϕ + a n e iθ , (7.9) where a n ≥ 0 and real. This case again effectively reduces to a one-dimensional moment problem. This time, however, it is not easily treated analytically. It is nevertheless straightforward to work out the desired bounds numerically, by reformulating it as a linear programming problem.
To be concrete let us consider the particular case ϕ = π 4 and θ = π 2 . We first notice that in this case µ 1 ≥ 0, where µ 1 = 0 corresponds to a degenerate case where we have a single QNM at a = 0, where a is defined in (7.9). We can then derive the bounds for ( µ 2 µ 1 , µ 3 µ 1 ) shown in Figure 8 . ) of the first three low-frequency expansion coefficients for the case when all QNMs lie on the line (7.9) with ϕ = π 4 and θ = π 2 . The boundary of the allowed region can be realized with two QNMs and one QNM. The red curve corresponds to the k = ∆ momentum BTZ correlator (7.8) as ∆ is varied. For small ∆ → 0 the moments approach µ 2 µ 1 → −∞ and Let us compare the derived bounds with the explicit results for the BTZ black hole for momentum k = ∆ (which corresponds to ϕ = π 4 ) We see that we get a two-sided bound for µ 3 µ 1 , but only a one-sided bound for µ 2 µ 1 . Moreover, the simplest example of the BTZ black hole correlator realizes arbitrarily negative µ 2 µ 1 . Let us explain the structure of the boundary of the allowed region in Figure 8 touches the blue boundary which is described by

Sectorial QNMs
Finally, we consider the case where the allowed region of support of QNMs is given by a sectorial region: Im ω n ≥ 1, Re ω n ≤ Im ω n . (7.11) It includes the two cases considered above but also allows for two-dimensional configurations of QNMs (not necessarily located on a line). We plot the results in Figure  9.
Let us comment on the structure of the boundary in Figure 9. The upper green line corresponds to a pair of QNMs As α → 0 the point on the curve moves to the left, while α → 1 corresponds to the cusp at (1, 1), where we have only one QNM at i. The blue segment corresponds to a single QNM at ω 0 + i, where 0 ≤ ω 0 ≤ √ 7 − 2. Finally, the lower green line is described by a pair of QNMs ρ(ω) = αδ(ω − ( An obvious question is: what happens to the bounds as we relax the condition (7.11)? For example, a natural minimal condition to consider is Im ω n ≥ 1 which simply corresponds to the statement that there is a minimal characteristic timescale for perturbations to decay. We have not explored thoroughly what can be said about the structure of the hydrodynamic expansion in this case. However, we observed that for ( µ 2 µ 1 , µ 3 µ 1 ) we do not get any nontrivial bounds in this case.

BTZ and Rindler space
Next we would like to discuss how the general structure discussed in the previous sections is realized in various examples. We start by considering the simplest examples in which the correlator is known analytically.
The two-sided correlator in the BTZ background is given by [6] (here β = 2π) ). The poles are located at ω n = k + i(2n − 2 + ∆) for n = 1, 2, . . . and the three images (−ω n , ω * n , −ω * n ), with residues The pole spacing is r = 2 at angle θ = π/2, and the subleading constant correction is se iϕ = k + (∆ − 2)i, so the OPE constraints (5.7) and (5.10) indeed hold. Note that the singularity sum rules (2.15) are not satisfied, since there is no curvature singularity in the BTZ black hole. For instance, truncating the first sum rule at n = n max gives This diverges as n max → ∞, reflecting the fact that the contour integral along C ∞ does not vanish. It is straightforward to verify that the BTZ two-sided correlator (8.1) can equivalently be written using the product formula (2.9) with ω n = k + i(2n + ∆) (note that G 12 (0, k) ̸ = 1 in (8.1), so that the product formula will differ from (8.1) by a normalization). For k = 0 the simple poles merge into double poles on the imaginary axis.
An example of known higher-dimensional thermal correlators comes from CFTs on S 1 β × H d−1 with β = 2π which can be conformally mapped to a Rindler wedge in R d , see e.g. [43]. The scalar two-point function on S 1 β × H d−1 was studied in [44] and is given by (β = 2π) which reduces to the scalar correlator in BTZ (8.1) for d = 2. Here k labels the principal series representation of SO(1, d − 1), see [44] for details. This has no zeroes, and has poles at ω = k + i − d 2 + ∆ + 2n − 1 for n = 1, 2, . . . and the corresponding images. Note also that the OPE expansion of (8.4) is different from what we derived in Section 4 because the spatial slice H d−1 has non-zero curvature. Given the poles, we can reconstruct the two-sided correlator using the product formula.

Examples in pure GR
In this section we consider examples that originate from matter minimally coupled to general relativity in AdS. The utility of the product formula (2.9) comes from the fact that it fixes the two-sided correlator in terms of the QNMs. The latter can be rather effectively found numerically [12] (whenever they are not known analytically), thereby paving the way to straightforwardly computing the two-sided correlator G 12 (ω).

General strategy
In this section we study QNMs numerically for various wave equations and use the (truncated) product formula to reconstruct the correlator. More precisely, we can truncate the product formula at a fixed number n max of QNMs, and then systematically improve it by increasing n max and/or attaching to it a tail of the asymptotic QNMs. The error made in this way can be estimated as follows. As discussed in Section 5, one way to reproduce the OPE (which is realized in many examples that we will consider) is if the poles are asymptotically linear. If asymptotically we have a single line of QNMs ω n ≈ re iθ n + . . . for 0 ≤ θ ≤ π 2 , then for ω ≪ ω nmax we have 11 For the special case θ = π/4 (which holds for AdS 5 /Schwarzschild), the error is order 1/n 2 max . If there are multiple lines of QNMs one has to sum over all of them. Now let us introduce an improved truncation scheme which improves the error estimate (9.2). Let us assume as above that we have ω n ≈ ω nmax + re iθ (n − n max ) for n max large and n ≥ n max . Assuming that we choose n max big enough, we can then improve the truncated solution bŷ 3) can be generalized to purely imaginary modes with asymptotically linear scaling with the mode number n. The improved correlator (9.3) correctly captures the large n behavior of QNMs, so that for ω n = re iθ (n − n max ) + ω nmax + f n γ e iρ + . . . we instead find which is improved compared to (9.2) by a factor of n −(1+γ) max for generic values of θ and ρ.
We begin this section with a toy model, the two-point function of R-currents in N = 4 SYM with zero spatial momentum. The correlator and the QNMs are known analytically and can be compared with the truncated counterpart. We then move on to scalar QNMs in a black brane background, which can be calculated using a variety of methods (see [10,11] and references therein). In practice we use the QNMSpectral package [12]. The truncated correlator is then obtained using (9.1). We further show the applicability to stress tensor correlators by studying metric fluctuations, and also consider the scalar correlator in a charged brane background. In these examples, the QNMs grow asymptotically linearly with the mode number n, so the truncated product formula can be improved using (9.3). For ω ∈ R we further numerically solve the bulk wave equation to extract the correlator and compare against the product formula with perfect agreement.

R-currents
Our first example is the two-point function of R-currents with zero spatial momentum 12 in N = 4 SYM [46], which is the simplest known example of a two-point function in a higher dimensional black hole. The retarded Green's function is where ψ is the digamma function and β = 2π. Using (A.7), we find the two-sided Wightman function The spacing between QNMs is r = √ 2 at angle θ = π/4, and the subleading constant correction vanishes, so the OPE predictions (5.7) and (5.10) are indeed satisfied with ∆ = 3.
-33 - Let us notice that the sum converges exponentially fast due to the 1 sinh(πn) factor. We now assume that the QNMs ω n = (1 + i)n are given as input. Using the product formula (2.9), we reproduce the expected result On the other hand, truncating to include the first n max QNMs, we plot in Figure 10 log(G nmax 12 (ω, 0)/G 12 (ω, 0)) at ω = 1 as a function of n max . In this case, the modes are given by ω n = √ 2e i π 4 n, with no subleading terms in 1/n. One then finds along the lines of (9.2) that For ω = 1 this decays as 1 12 n −3 max , in agreement with the observed values in Fig. 10.

Scalar QNMs in black branes
We now consider a more nontrivial example for which we resort to numerics: the scalar two-point function with ∆ = 4 on S 1 β × R 3 dual to a massless scalar field in an AdS black brane background 13 . This is equivalent to the scalar perturbation of the metric in pure GR and the correlator will therefore be the same as for ⟨T xy T xy ⟩ 14 . The metric is given by 15 where f (u) = 1 − u 2 and u = r 2 0 r 2 , r 0 being the horizon radius. In (9.12) without loss of generality we set the AdS radius and r 0 to 1, which fixes the temperature β = π. The wave equation for a scalar field ϕ(u, t, z) = ϕ(u)e −iωt+ikz dual to an operator with ∆ = 4 is then given by Imposing infalling boundary conditions ϕ(u) ∼ (1 − u) −iω/4 at the horizon u → 1, we can extract the retarded correlator and hence the two-sided correlator using (A.7). For further details see e.g. [6,47,48].
In Figure 11a we display the first 7 QNMs obtained from (H.2) for a scalar field with k = 0 and k = 2. The truncated product formula with n max = 27 modes and improved with r = 2 √ 2 and θ = π 4 is shown in Figure 11b, where we also compared against the result obtained from numerically solving the ODE (9.13) in Mathematica. 13 It is straightforward to do the same computations for generic values of ∆. 14 Up to overall normalization. 15 We follow the conventions in [47].
The picture discussed here does not change qualitatively as we change ∆ or d. As we vary ∆ away from 4 we get more contributions in the large ω OPE expansion. We considered ∆s for which both the stress tensor and the double-trace stress tensor operators appear, and checked that the product formula correctly reproduces the numerical solution. One difference as we change d is that the asymptotic angle at which the QNMs approach infinity changes to e −i π d , see (6.6). As a noteworthy example, it has been recently shown [49] that QNMs of the gravity dual to the BFSS matrix model [50] can be computed using the black brane geometry in d = 14 5 .

Metric fluctuations
In this section, we will observe numerically that the product formula (2.9) correctly reproduces the correlator for metric fluctuations. The metric perturbations h µν are classified according to their symmetry properties with respect to rotations along the transverse directions of the propagation. This leads to three gauge invariant combinations Z i = Z i (h µν ), corresponding to three independent scalar functions determining the stress tensor two-point functions at the boundary. However, note that the analysis of the zeroes of the Wightman function in Appendix C is only applicable for the scalar wave equation, since the potential for metric fluctuations contains singularities at positions depending on the spatial momentum and frequency (see [51] for a recent discussion). It would be interesting to generalize the analysis to metric fluctuations as well. We leave this problem for future work. , where the lines were obtained from the improved truncated product formula with nmax = 13, r = 2 √ 2 and θ = π 4 . The points were obtained from NDSolve. For small k the peak is located at ω = k 2 4 , due to the hydrodynamic shear mode [48]. As above, we have chosen the normalization for the product formula to agree with the result from NDSolve at the lowest value of ω used in the latter, as well as chosen the absolute normalization arbitrarily.
The scalar perturbation is equivalent to a ∆ = 4 scalar, which was studied in Section 9.3. For the shear channel Z 1 and sound channel Z 2 , the wave equations are given by [47] where Here f (u) = 1 − u 2 and we have set β = π. Given the gauge-invariant observables Z i , the stress tensor two-point functions for various polarisations can be extracted, see [47,48] for details.
In Figure 12a we have plotted the QNMs in the shear channel for k = 6 10 and k = 12 10 . In Figure 12b we plotted the imaginary part of the retarded correlator G tx,tx for ω ∈ R using the improved truncated product formula and compared against the numerical solution. Similarly, in Figure 13a we plotted the QNMs in the sound channel due to the hydrodynamic sound mode [48]. Here we have chosen the normalization for the product formula to agree with the result from NDSolve at the lowest value of ω used in the latter, as well as chosen the absolute normalization arbitrarily.
for k = 0.6 and k = 3, and the comparison between the improved truncated product formula and the numerical solution of the bulk wave equation in Figure 13b. The asymptotic behavior of the QNMs for metric fluctuations was obtained in [28,29] in general dimensions. For both the shear and the sound channel these are given by which is the same as a scalar with ∆ = d. Note that in both cases the hydrodynamic modes are not included. It is then seen that the OPE of T tx T tx is correctly reproduced, since (9.23) where we used (5.2), (5.7) and (5.10) together with (9.22), ∆ = d and the fact that there is a purely imaginary hydrodynamic mode and its image contributing as − 2 ω . Likewise, we find that . . , (9.24) in agreement with expectations for the T tt T tt OPE. The term − 4 ω is again due to the hydrodynamic mode and its images.

Charged black brane
In this section we consider a scalar correlator in a charged state, dual to a charged black brane. The QNMs in a charged AdS black hole were analyzed in [28,52,53]. The main new feature in this case compared to the uncharged case is the presence of an infinite line with purely imaginary modes. By modifying the truncated product ansatz (9.1) to include such imaginary modes we will again verify good agreement with the numerical solution on the real line. The metric for a charged black brane is given by, see e.g. [54], where f (r) = 1 − µr −4 + Q 2 r −6 . We fix the temperature T = 1 π (1 − Q 2 2 ) by setting the horizon at r = 1, which further implies that µ = 1 + Q 2 .
In Figure 14a we plot the QNMs with Q = 1 and k = 0, 2, which contain a new infinite line of imaginary modes that were absent in the uncharged case. In Figure  14b we compare the numerical solution to the improved truncated product formula, modified to include pure imaginary modes.
To compare with (5.15) we fit the two lines of QNMs. For the k = 0 case, using about 50 complex modes and 80 imaginary modes we find

Higher derivative corrections
The examples considered so far have all been at infinite λ, but the product formula (2.9) is equally applicable when higher derivative corrections are taken into account.
In this section we analyze several instructive examples of higher derivative terms in the Lagrangian, and confirm the product formula and the predictions from the OPE.

Gauss-Bonnet black holes
Let us first calculate the QNMs in Gauss-Bonnet gravity, from which we compute the two-point function and compare with the numerical solution to the wave equation. The Gauss-Bonnet QNMs were previously discussed in [55,56]. Consider a black brane in Gauss-Bonnet gravity, with the metric [57,58] , λ GB is the Gauss-Bonnet coupling and The AdS radius is given byL 2 = L 2 /f ∞ . The Hawking temperature is From now on we set L = 1 and further fix the temperature in terms of λ GB by setting r + = 1.
The QNMs can be calculated as above by passing to Eddington-Finkelstein coordinates 16 , and are shown for a scalar operator with ∆ = 4 in Figure 15a. Likewise, the numerical solution of the correlator is shown in Figure 15b and compared against the improved truncated product formula.
A numerical fit gives for the two lines of QNMs (again we restrict to the k = 0 case) ω 1,n ≃ 7.77in , ω 2,n ≃ (2.53 + 2.54i)n . This gives 4π sin θ 2 r 2 + 2π r 1 ≃ 3.30, (10.5) while the inverse temperature reads consistently with (5.15). To check the subleading sum rule more QNMs are needed, and we leave this problem to future work. It would be also curious to generalize this analysis to metric perturbations [59,60].

ϕ 2 W 2 coupling
We now introduce a higher derivative term that displays novel behavior for QNMs. In particular, we will find a line of imaginary modes which asymptotically behave as ω n ∝ in 3 . Let us consider the Lagrangian where W is the Weyl tensor. A similar interaction W 2 ϕ was analyzed in [61,62] as a model for higher derivative corrections to the one-point function. We consider again the black brane metric (9.12) with β = π. The squared Weyl tensor in these coordinates is

The equation of motion reads
The coupling α is measured in units of the AdS radius L = 1. We plot the QNMs for The OPE expansion predicts that the line of complex QNMs takes the form ω 2,n ≃ re iθ n + se iϕ n 1/3 + c . (10.11) A numerical fit involving about 70 complex modes confirms the predictions and sets (10.13) consistently. Finally, we also confirm the prediction made in (5.20), that is Note that b, c will contribute at higher orders in the ω −1 expansion. Here we only considered QNMs for α ≲ 0. More generally, the structure of the QNMs as a function of α displays a rich structure. To gain a qualitative understanding, let us consider the WKB limit ω, ν, α → ∞ with α/ν =ᾱ and ω/ν fixed. The WKB potential reads where r = 1/u 2 is the usual radial coordinate. We start by consideringᾱ < 0. In this case, since the α term dominates close to the singularity, the potential goes to +∞ as r → 0. We can distinguish 3 cases: The WKB potential is positive definite and has a global minimum at the horizon r = 1 (z → ∞). Here Therefore we can still apply (3.13), and we find a new line of poles at ω n = −(2 + 4n)i for n ≥ 0.
•ᾱ < α c . The global minimum moves outside the horizon to a finite value of z where V (z) is negative. This minimum leads to linearly spaced (∂ 2 z V (z min ) ̸ = 0) purely imaginary modes in the upper half plane. These modes correspond to bound states, since for ω n = i|ω n | the ingoing solution behaves as e −|ωn|z as z → ∞. As mentioned in Section 3, they are associated to an instability of the black hole.
• α c <ᾱ < 0. The minimum moves inside the horizon, and the potential hosts virtual bound states corresponding to linearly spaced QNMs on the lower imaginary axis.
On the other hand whenᾱ > α * > 0, with α * ≃ 0.12, the potential develops a metastable minimum outside the horizon, corresponding to weakly damped QNMs. If α * >ᾱ > 0 the metastable minimum becomes complex and the QNM spectrum is qualitatively unchanged compared to the α = 0 case. Note that the WKB analysis holds for ω ∼ ν, while to compare with the OPE predictions we need to analyze the asymptotic behavior of QNMs in the limit ω ≫ 1, ∆, k. The qualitative features of the QNMs are the same as discussed above. The main difference is in the spacing between modes in the case α c <ᾱ < 0. As shown in Figure 16b, in the large ν regime the separation between these modes indeed starts linearly, but increases gradually as ω becomes larger than ν, until we find ω n ∝ in 3 as noted previously.

Beyond the strong coupling limit
In the bulk of this paper we have discussed the analytic properties of the two-point function in the holographic regime of infinite N and infinite λ. In order to understand the range of applicability of our results, it is important to analyze which of these properties survives at finite coupling and finite N .
Let us first discuss the effect of 1/N corrections. At one loop, the two-point function develops a branch cut, corresponding to late-time hydrodynamic tails [39]. As a result, the Wightman function is not an analytic function of frequency, so its poles and zeroes are no longer enough to determine the full function. It follows that the results presented here have limited applicability at finite N (at least without some modification).
We now turn to the case where the coupling is finite but N is infinite. In the rest of this section, we will present several instructive examples of lower-dimensional models where the holographic properties (meromorphy and no zeroes) extend to arbitrary coupling. We first consider the SYK model [63][64][65][66], a chaotic 0+1 dimensional theory with a built-in disorder average. We then discuss a 1+1 dimensional generalization of the SYK model known as the SYK chain [67,68]. In both cases, G 12 is meromorphic and has no zeroes, which shows that the analytic properties of a holographic CFT have a chance of extending beyond the holographic regime. We conclude with several examples of non-chaotic theories, in which both properties are violated.

SYK
The SYK model is a quantum mechanical system of N fermions, with a random coupling involving an even number q fermions at a time. The Hamiltonian is given by [63][64][65][66] We consider the limit N → ∞ with q fixed, followed by the limit q → ∞. In this regime, the two-sided Wightman function is given by [63,69] (11.2) Fourier transforming, we find coupling. In Figure 17 we plot G 12 (ω) for several values of the parameters. The poles of this function are all close to, but not necessarily on, the imaginary axis. Moreover, there are no zeroes, so this model provides a nontrivial example where the properties of the holographic Green's function extend beyond the holographic regime. We have checked numerically that (11.5) has no zeroes for various parameter values. It would be interesting to prove this analytically.

Non-chaotic theories
In theories which are not chaotic, we do not expect the no-zeroes or meromorphy properties to hold. For example, consider N = 4 SYM at zero 't Hooft coupling, as analyzed in [7]. Setting β = 2π, the Wightman function for O = Tr F 2 at zero spatial momentum is given by . (11.9) This is an analytic function of ω. However, we see that there is a zero at ω = 0, so the no zeroes property no longer holds. When a nonzero spatial momentum is introduced, the correlator develops branch cuts in the complex ω plane [7], (11.10) It follows that the analyticity property is broken at finite momentum. Let us briefly discuss the properties of this function as we increase k. First of all we see that it has a zero at ω = k (the same that we observed for k = 0). At this point we have ω − k = 0 and therefore our local operator carries a null momentum, and thus can be interpreted as a light-ray operator that annihilates the vacuum [70]. Resolving the identity as in [71], we can write the thermal expectation value as a sum of double commutators which have zeros at certain integer-spaced values of scaling dimensions. Due to the integrability of the free theory the whole spectrum consists of integer-valued scaling dimension operators and we get zero. For the same reason we do not expect G 12 (ω, k) to have zeros for real ω in interacting theories. For k > ω the correlator decays exponentially fast with the rate e −πk as expected on general grounds [6,34]. It is important to understand whether analyticity and no-zeroes are restored at small but finite 't Hooft coupling. A resummation is likely necessary to address this problem, since perturbation theory in λ breaks down at late times [72]. An alternative possibility is that there is a transition between the holographic and free field theory behaviors at some intermediate coupling (see also [56,73,74] for related discussions).
Another simple example of a nonchaotic system is generalized free field theory. The generalized free field result takes a very simple form for the case S 1 × R d−1 , where the two-sided correlator is given by the vacuum block in (4.4), see [32]. 17 Ignoring the θ(ω 2 − k 2 ), we see that the two-sided correlator exhibits branch cuts at ω = ±k and a power-law behavior along the imaginary axis ω → ±i∞ in contrast to the result in the black hole phase. Let us also mention that below the Hawking-Page transition the correlator to leading order is given by the generalized free field result on S 1 × S d−1 .
It is an interesting question what happens for interacting QFTs in AdS at finite temperature. Perturbatively, we do not expect to see any change [75]. However it could be that perturbation theory breaks down at large (but not too large) times, which happens in other related cases [72].
Finally, vector models in d = 3 at large N are characterized by weakly broken higher spin symmetry [76], and in this sense are not chaotic. This manifests itself for example through the fact that higher spin currents have anomalous dimension O(1/N ). Another way to see that these theories are not chaotic is that the Lyapunov exponent is O(1/N ) in this case [77]. Thermal correlators in the O(N ) models have been studied for example in [78][79][80]. The basic result is that they do not exhibit quasi-normal modes as expected in chaotic theories, and they have branch cuts in the ω-plane. Similarly, we expect Chern-Simons matter theories at finite temperature [81] not to exhibit meromorphy and no-zero structure, but we have not explicitly checked this.

Conclusions and further discussion
The emergence of semi-classical gravity in holography is associated with strongly coupled quantum dynamics. Thermal correlators are particularly interesting observables in this regard since the dual geometry contains a black hole. However, understanding thermal correlators in strongly coupled systems is a challenging task. In particular, developing efficient bootstrap methods for probing them is an open problem.
In this paper we have considered the holographic thermal two-point function and we noted that it exhibits an intriguing and non-obvious property: thermal two-sided correlators do not have zeros in the complex energy plane. From the dual geometry point of view, this property is very closely associated with the presence of the black hole horizon. Together with meromorphy, the no-zero property leads to the product representation of the thermal correlator (2.9). We have derived this representation for holographic systems. However, we believe it may hold more generally.

Thermal product hypothesis (TPH):
The two-sided thermal two-point correlator in a chaotic large N system is a meromorphic function with no zeros in the ω-plane. As such it is given by a product over its poles.
From holography we expect that large N chaotic theories at finite temperature are described by some kind of stringy black holes. Then the content of TPH is that at the level of the thermal two-point function, stringy black holes behave like ordinary black holes. We checked this property analytically in the SYK model and numerically in the SYK chain in Section 11. In particular, in the latter case it looks quite nontrivial; it would be instructive to derive the no-zero property of correlators in the SYK chain using an analytic argument. We also checked in Section 10 that TPH holds if we include certain higher derivative corrections to GR.
A few comments are in order. In the context of CFTs, we expect that TPH can only be valid for theories where the anomalous dimensions of higher spin currents are O(1), as opposed to O(1/N ). In particular, TPH does not apply to theories with slightly broken higher spin symmetry [76]. The requirement that the theory is chaotic is necessary since there are branch cuts for the correlator in free field theory or in the planar O(N ) model [80].
We expect that TPH may hold for CFTs on S 1 × R d−1 , or S 1 × S d−1 above the Hawking-Page phase transition. On S 1 × S d−1 , the assumption of large N is required since otherwise the two-sided correlator is given by a sum of δ-functions. On S 1 × R d−1 , large N is required since otherwise we expect to have branch cuts due to long-time tails, see e.g. [80]. We require the theory to be above the Hawking-Page transition so that it is described by the black hole geometry. Finally, cuts appear in a simple kinetic theory description of thermal correlators in weakly coupled gauge theories, see e.g. [74]. However, to the best of our knowledge their existence has not been rigorously established, see the discussion in [7].
We have also reviewed the observation of [9], related to the previous work [16,17], that holographic two-sided correlators decay exponentially fast as ω → ±i∞ in d > 2. This behavior is associated with a light-like geodesic that bounces off the black hole singularity. We noticed that when combined with the prediction from the OPE there is a dispersive way to express this behavior through the black hole singularity sum rules (1.4). Assuming the product representation of the thermal correlator, we see that the sum rules are trivially satisfied by having a family of QNMs that go to infinity in the complex plane at some angle ω ∼ e iθ with θ ̸ = π 2 . At finite string coupling λ we expect the sum rules to still hold in some range of the ω-plane for which the gravity approximation applies. In this case we can consider finite energy singularity sum rules: where Λ is a finite energy scale. We can imagine that a similar version of the sum rules might hold at large but finite N .
Another open problem is to develop computational techniques to probe the region of large imaginary frequencies in string theory (or gauge theories, see e.g. [72]). In particular, tidal effects should become important near the singularity. These effects are responsible for resolving singularities in the two-point function on the bulk light-cone in one-sided correlators [82,83], and it seems plausible that tidal effects change the behavior of the two-sided correlator at large imaginary frequencies as well. However, the black hole geometry receives large corrections at a string length away from the singularity, so new insights are likely needed to address this problem.
In Section 5 and Section 6, we explored how the structure of the QNMs is constrained by the OPE. In Section 7 we studied how the known structure of the QNMs puts constraints on the low-energy or hydrodynamic expansion of the correlator. The basic observation here is that the product representation of the correlator and its nozero property lead to a convenient dispersive representation of the two-point function, see (5.2). Moreover, the no-zero property translates into the fact that the discontinuity of the correlator is non-negative. The hydrodynamic expansion then turns into a moment problem for QNMs, and depending on the structure of QNMs bounds on the hydrodynamic expansion can be derived. It would be very interesting to understand such constraints more systematically. Conversely, the known data on the hydrodynamic coefficients can be used to put bounds on the geometry of QNMs (again assuming the no-zero property).
Another challenge that could provide important insights into the finite coupling structure of the correlators and the dual geometry is to find models in which the structure of the QNMs can be computed in a controllable setting and the emergence of gravitational features can be understood, see e.g. [84]. It would also be interesting to see if the product representation of thermal correlators generalizes beyond the two-point function [85,86].
Finally, the results of [87,88] suggest that instanton corrections to weakly coupled thermal correlation functions can be interpreted in terms of AdS black holes. This is reminiscent of the fact that instanton corrections to vacuum correlators at weak coupling can be written as an integral over AdS, where the AdS space is simply the moduli space of instantons [89][90][91][92]. It would be interesting to explicitly compute the first instanton correction to the thermal two-point function at weak coupling in N = 4 SYM, and to analyze its properties.
The two-sided Wightman correlator is In frequency space, the two-point functions are related as where in (A.5) ω approaches the real axis from above. Note that for real ω, G 12 (ω) = ρ(ω)/(2 sinh(βω/2)), where the spectral density ρ is defined by ρ(ω) = 2Im G R (ω). In general, (A.5) defines a distribution and not a function, but for holographic correlators G R (ω) is a meromorphic function with no poles on the real axis, and we can rewrite (A.5) as follows This relation can now be analytically continued to ω ∈ C and defines the meromorphic G 12 (ω) studied in this paper. The "no-zero" property discussed in the main text translates to the statement that the equation has solutions only for Matsubara frequencies The Euclidean correlator defined at Matsubara frequencies ω n = 2πi β n is related to the retarded correlator as follows (A.10)

B Probing the ω-plane in real time
From the point of view of an experimentalist, the complex ω-plane discussed in this paper is not easily accessible. On the other hand, QNMs are known to control the late time behavior of the retarded two-point function which can be accessed more easily. In fact, for the correlators we consider QNMs control correlators at all times. By this we mean the following. Let us consider the retarded correlator for some positive non-zero time t > 0, We can close the contour into the lower half-plane to get where the condition t > 0 is necessary to drop the arc at infinity. In particular, there could be distributional terms localized at t = 0 which are not captured by QNMs. These correspond to subtractions in dispersion relations for G R (ω). Let us next express the residues in terms of QNMs. For this purpose we write In writing the formula above we used meromorphy of G R (ω), see Appendix A and (A.7). Considering next ω in the lower half-plane we get Restricting thus to QNMs with Re ω k > 0 we get the following expression for the retarded two-point function in terms of the QNMs, The sum converges for any t > 0, with more and more QNMs becoming important at earlier times.
As an example we consider the thermal two-point function of R-currents, see (9.6). In Figure 18 we plot G R (t) and its approximation by a few QNMs. We see that as we decrease t from infinity, higher QNMs become important in an ordered fashion. Note also that by measuring the residue of the first two QNMs we get access to some information about the high-energy tail through (B.5). Indeed, in the ratio of the residues G 12 (0) cancels and we get some prediction for the infinite products.

C Analytic properties of AdS wave equations
Here we review some analytic properties of holographic correlators based on the bulk equations of motion [8,9]. In particular, we will recall the argument that Wightman correlators have no zeroes, which is a consequence of scattering theory in quantum mechanics [94].
Consider the scalar wave equation (3.3) in an AdS black hole background. The normalizable mode g and non-normalizable modeg are specified by the boundary conditions at z → 0 18 g(ω, z) ∼ z Various properties under conjugation and ω → −ω can readily be found from these boundary conditions. These solutions can further be expressed in terms of each other, in particular where f (ω) is the so-called Jost function given by the Wronskian Since the Wronskian is independent of z, we can evaluate it at z = 0. This gives Therefore the analytic properties of the Jost function are the same as those of the ingoing mode h R (ω, z). We further consider the physical solution proportional to the normalizable mode where C(ω) is fixed by the normalization and δ is the phase shift. It was shown in [9] that and that the Wightman correlator is given by (C.9) The analytic properties of the Wightman function can therefore be read off from those of the Jost function f (ω). In particular, for a regular potential with the asymptotic behavior (3.5), the Jost function is a meromorphic function with simple poles at the Matsubara frequencies ω = −i 2πn β with n = 1, 2, . . .. It therefore follows that the Wightman function has no zeroes. Moreover, the zeroes of the Jost function correspond to the QNMs.
Let us review the argument that the only poles of the Jost function appear at the Matsubara frequencies. For simplicity, we begin by assuming that ν 2 = 1 4 so that the potential behaves like V ∼ z ϵ−2 as z → 0 with ϵ > 0. In this case we can borrow the techniques of scattering theory with zero angular momentum l = 0. It is straightforward to generalize the results to generic real values of ν 2 > 0, as we will see later. The ingoing solution can be written as a Volterra equation as follows, where we have multiplied V (z) by a parameter γ, which we will eventually set to one. The solution to (C.10) can be shown to define an absolutely convergent power series in γ at finite z, see e.g. Sec 12.1 in [94], if Assuming the potential is regular for finite z, only the large z region in (C.11) could give rise to a divergence. In particular, for any decaying potential α is finite for Im ω > 0. Because of the presence of the horizon, we restrict to potentials which are exponentially decaying as V = n a n e − 4πn β z as z → ∞. For the purposes of studying whether α is finite or not, we can then choose z > z 0 for some large enough z 0 and replace the potential with the sum of exponentials. Consider therefore The integral converges for Im ω > − 2πn β , corresponding to the nth Matsubara frequency. The region of analyticity can however be extended beyond this point with simple poles at ω = −i 2πnm β for m = 1, 2, . . .. This can be seen explicitly by solving (C.10) order-by-order in γ, with a new simple pole arising at each order, see Section 12.1.1 in [94] for details 19 . It then follows from (C.5) that the Jost function f (ω) is meromorphic with simple poles at the Matsubara frequencies.
We note that at isolated points in parameter space, the residue of the pole of the Jost function at one of the Matsubara frequencies might vanish. For example, for the black brane potential we have It follows that a 1 = 0 when This coincides with the first pole-skipping point [95][96][97][98][99]. For the special value (C.14) of k, G 12 (ω) has a pole at the first Matsubara frequency. Consider now the case ν 2 ̸ = 1 4 . The solution for h R (ω , z) can be found from the where G is the Green's function One can show that G is an analytic function of ω and that the solution to (C.15) has the same analytic structure as h R (ω, z). The solution h R,ν then has the same poles as h R .
Consider the Jost function f (ω) = W (h R , g) at spectral points ω = ω n ̸ = 0, where f (ω n ) = 0 by definition. One can show that f ′ (ω n ) ̸ = 0 for ω n in the lower half plane and on the positive imaginary axis, corresponding to resonances and bound states respectively, and the Jost function therefore has simple zeroes at those spectral points [100]. This follows from the use of the wave equation to relate f ′ (ω 0 ) to an integral of the square of h R (ω 0 , z), which converges and is non-zero for resonances and bound states. For virtual states with ω n = −i|ω n | this no longer is true and it is in principle possible that the zeroes can be of higher order. This arose in the study of higherderivative corrections [73], which found that QNMs on the imaginary axis for certain values of the coupling collided and branched into complex QNMs.

D OPE in momentum space
Let us consider the OPE for a pair of identical scalar operators at finite temperature, where we set β = 1. We would like to perform the Fourier transform of this OPE expansion to derive the asymptotic expansion of the two-sided correlator G 12 (ω, k). We start with the spatial The expression for G 12 for a theory dual to a black brane geometry in the bulk reads are the parameters of the wave equation of a scalar in the black brane background in the variables The parameters u, a are related to each other through F N S (a i , a, t) as follows: F N S is given as a convergent power series in t. Since for the black brane t = 1/2, (E.1) cannot be analytically expanded perturbatively for t ∼ 0. However since 1/2 is a small number and the series is convergent, one can still truncate the series to a given order and compute (E.1) numerically. In Figure 19 we compare Im G R obtained from the product formula with the numerical evaluation of (E.1). Equation (E.1) dramatically simplifies as |ω| → ∞. The key observation is that for large ω (E.3) reduces to an hypergeometric equation in the variablex = (x − t)/(1 − t), In order for this to be consistent with (E.3) and (E.4), at leading order in ω we must have Figure 19: (a) A plot of Im G R (βω) in the complex βω plane for k = 0 and ∆ = 5 2 from the product formula. Here the first 3 modes have been used. (b) Here Im G R (βω) is computed from a numerical evaluation of the exact result found in [37]. The F N S series is truncated at order t 5 . The slight mismatch in the position of the second pole is expected to be resolved by truncating the series at higher orders. The white spots along the imaginary axis are due to the presence of unphysical poles in the relation (E.4).
Note that the previous expression gives F N S at a specific value of a = a t , so it cannot be used to evaluate ∂ a F N S in (E.1). By solving (E.4) order by order in t and then expanding for large ω we find ∂ a F N S = icβω + O(ω 0 ) , c ∈ R . (E.7) Substituting (E.6) and (E.7) in (E.1) and expanding for large positive ω we get the leading order with its nonperturbative corrections. On the other hand, to compute power law corrections we would need corrections to (E.6).
Up to the first nonperturbative correction we get (E.8) Comparing with (4.12) with d = 4 and β =β we find the following prediction for c, The prediction is confirmed numerically, see Figure 20.

F Subleading corrections to black brane QNMs
In this appendix we compute the large n asymptotics of QNMs in the AdS black brane to the next subleading order beyond the constant term. We start with the black brane potential at k = 0, Following [28,29,116], we need to compute the series expansion of the potential near the singularity. We work in the conventions of [116], so that the tortoise coordinate is defined by dz = dr/f (r), with z = 0 corresponding to r = 0. The expansion of the tortoise coordinate around r = 0 is Plugging into (F.1), we find The asymptotic expansion of the QNMs can now be computed to first order in perturbation theory, by matching the solution at the singularity to the normalizable solution at the boundary and imposing ingoing boundary conditions. This computation was done in [116], and we can read off the answer from Equation 100 in that paper.

(G.5)
Let us compare this with the Wightman function of the normalized scalar primaries ⟨O(x)O(0)⟩ = c O x 2d of scaling dimension ∆ = d, . (G.6) From this we find that This is the normalization for the scalars that we need to use in our computations for ∆ = d in order to compare to the stress-energy computation. The shear viscosity η is defined as [118] η = lim In a theory with a gravity dual it is given by where σ abs (0) is the black hole cross-section at zero frequencies, which is given by the area of the horizon (per unit CFT volume) σ abs (0) = Area. (G.10) Recall that the black hole entropy density is in the same way s = Area 4G N . This produces the famous relation η/s = 1 4π [119]. Therefore the prediction is where G R (ω, 0) is the scalar two-point function with the normalization used in this paper, namely ⟨O(x)O(0)⟩ = 1 x 2d . Recall that we have [35] where we set R AdS = 1.
In this way we get the following prediction for the holographic correlator for ∆ = d

H Computation of QNMs
In the main body of this work we have numerically calculated the QNMs in various examples using the publicly available QNMSpectral package for Mathematica [12]. The method used is described in detail in [12], which we briefly review here. In holographic theories, the calculation of QNMs boils down to solving wave equations on top of a black hole background in AdS with specified boundary conditions, namely ingoing boundary conditions at the horizon and normalizable boundary conditions at the AdS boundary. As discussed in detail in Appendix C, the QNMs are solutions for discrete values of ω ∈ C such that the solution with ingoing boundary conditions at the horizon is proportional to the normalizable mode and the Wronskian therefore vanishes. The method of [12] finds a solution to the wave equation numerically by discretizing the radial direction on a grid of n + 1 points, reducing the wave equation to a generalized eigenvalue problem of (n + 1) × (n + 1)-matrices. The boundary conditions are imposed by expanding in a set of functions which manifestly obey the specified boundary conditions. The output is n + 1 eigenvalues, some of which correspond to physical QNMs while some are unphysical. To filter out unphysical solutions, following [12] we compute the QNMs twice with different grid sizes 20 and keep only those that appear in both cases.
The computational time and the number of physical QNMs that are obtained depend on the grid size n and the precision used when solving the generalized eigenvalue equation. These aspects were studied for AdS-Schwarzschild in Appendix A of [12]. It was found that when the precision is set to n/2, the computational time grows roughly like t ∼ t 0 n 3.3 for large enough n, and the number of physical QNMs grows linearly with n.
Let us review in a bit more detail how this works for the ∆ = 4 scalar ϕ(U, t, x) = ϕ(U )e −iωt+ikz in a black brane background. A convenient way to impose the correct the boundary conditions for QNMs is to pass to Eddington-Finkelstein coordinates where F (U ) = 1 U 2 − U 2 and G(U ) = − 1 U 2 . We further consider the wave equation in terms of the field ψ(U ) = U −3 ϕ(U ), which is given by The main points are the following: 1) after going to EF coordinates, at the horizon the ingoing mode approaches a constant ψ ∝ 1 + . . ., while the outgoing mode oscillates rapidly ψ ∝ (1 − U ) iω 2 + . . . and 2) after rescaling by U −3 , at the boundary the nonnormalizable mode diverges as ψ ∝ U −3 + . . ., while the normalizable mode goes to zero linearly, ψ ∝ U + . . ..
In order to discretize and solve (H.2) numerically, [12] uses a pseudo-spectral method where a function f (x) is approximated on a grid x i with i = 0, 1, 2 . . . n as where C j (x i ) = δ ij . In particular, the so-called cardinal function C j (x) is given by Inserting the choice of grid points (H.5) into (H.4), the functions C i can be written as a linear combination of Chebyshev polynomials T i (x), see [12] for further details.
In particular, with this choice the interpolation (H.3) can never approximate a function that diverges close to the boundary and oscillates rapidly as we approach the horizon. The boundary conditions corresponding to QNMs are therefore automatically implemented. By separating the ω 0 and the ω 1 terms in (H.2), the wave equation evaluated at the grid points can be put into the form of an (n + 1) × (n + 1) generalized eigenvalue equation. The n + 1 eigenvalues then corresponds to candidate QNMs ω n . However, for a fixed number of grid points n, only a subset of the (n + 1) eigenvalues correspond to actual QNMs. To select the physical solutions, one can do the computation with different grid sizes n 1 and n 2 , and keep only the solutions which agree between the two [12]. 22 21 These live in [−1, 1] but can shifted and rescaled to [0, 1] to apply to (H.2). 22 This is not guaranteed to always correctly select the QNMs but in practice seems to work well.