Holographic open quantum systems: Toy models and analytic properties of thermal correlators

We present a unified picture of open quantum systems, the theory of a system probing a noisy thermal environment, distilling lessons learnt from previous holographic analyses. Our treatment is applicable both when the system is coupled to short-lived (Markovian), and long-lived (non-Markovian) environmental degrees of freedom. The thermal environment is modeled using an asymptotically AdS black hole, and the systems of interest are simple probe field theories. The effective stochastic dynamics of the system is governed by real-time thermal correlators, which we compute using the gravitational Schwinger-Keldysh (grSK) geometry. We describe the structure of arbitrary tree-level contact and exchange Witten diagrams in the grSK geometry. In particular, we argue, that all such diagrams reduce to integrals supported on a single copy of the exterior of the black hole. The integrand is obtained as a multiple discontinuity of a function comprising ingoing boundary-bulk propagators, monodromy functions which appear as radial Boltzmann weights, and vertex factors. These results allow us to deduce the analytic structure of real-time thermal n-point functions in holographic CFTs. We illustrate the general statements by a two-dimensional toy model, dual to fields in the BTZ background, which we argue captures many of the essential features of generic open holographic QFTs.


Introduction
The basic principles underlying the effective dynamics of open quantum systems are wellunderstood, cf., [1]. However, actual progress in the field theoretic context has been limited, in part, for technical reasons. In weakly coupled systems there is no obvious scale separation (environmental dynamics in such situations is slow) leading to long-time effects in the effective description. Motivated by this, [2] argued for using holographic field theories as thermal environments to model open quantum dynamics. We shall further explore this framework and provide an illustration of the general lessons learned from the holographic analysis in the context of an analytically tractable low-dimensional toy-model.
Our set-up is the following: a finite temperature holographic field theory provides a strongly correlated environment, which we probe by a simple quantum system, a free scalar field coupled to a gauge invariant single-trace operator. The effective field theory of the scalar probe depends on the real-time thermal correlation functions of the holographic theory. These are computed on the complex gravitational Schwinger-Keldysh geometry introduced in [3] (important earlier work on holographic computations of Schwinger-Keldysh correlators includes [4][5][6][7]).
The authors of [2] focused on the case where the single trace operator was a conformal primary. These operators typically relax within a thermal timescale. This fast relaxation was the motivating factor behind the choice of holographic environments as an ideal test-bed for analyzing open quantum dynamics. However, holographic environments also have slowly relaxing modes, generically corresponding to components of conserved current operators. They lead to hydrodynamic behaviour of the holographic plasma. Our aim is to give a unified treatment which subsumes both long-lived and short-lived operators, which following [8] we refer to as non-Markovian and Markovian, respectively.
For example, the energy-momentum tensor of a thermal system has both such modes. Tensor polarizations are Markovian (short-lived), while transverse vectors describing momentum diffusion, and longitudinal scalar polarizations capturing the phonon mode are non-Markovian (long-lived). 1 In [8,9] the authors constructed the Gaussian effective action for these modes. 2 The hydrodynamic behaviour of course has been well-known for a long time at the level of dispersion relations [12,13].
At this stage we should clarify our terminology a bit. Conventionally, in the open quantum system literature, non-Markovian dynamics refers to situations where there is long-time temporal correlation leading to memory effects. Typically, this arises because the environment that is being integrated out has such long-time dynamics. In our examples, we are explicitly coupling our probe systems to operators of the thermal environment whose correlation functions exhibit long-distance and large-time effects. Furthermore, these arise because of the presence of low-lying Goldstone type collective modes. We will therefore adapt the definition non-Markovian systems to characterize modes of the environment itself. In particular, we refer to modes of the environment that have large infra-red effects as non-Markovian. As noted above, a single (gauge invariant) operator of the environment might itself have contained Markovian and non-Markovian components, which one would have to disambiguate carefully to obtain a local open effective action.
The holographic analysis was aided by repackaging different polarizations of currents in terms of diffeomorphism invariant combination of gravitational perturbations [14]. These combinations end up being non-minimally coupled scalar fields in the background with their kinetic terms modulated by an auxiliary dilaton. For massless fields like the graviton, the Markovian versus non-Markovian nature can be deduced from the asymptotic behaviour of 1 The polarization decomposition is with respect to the little group of rotations in the space orthogonal to a fixed vector (the direction of motion). 2 An analogous analysis for finite charge density environments was carried out in [10,11]. the dilaton. 3 It was useful to characterize the distinction using a single parameter, the Markovianity index M, as in the aforementioned works. As explained in [9], for massless bulk fields, M ≥ 1 corresponds to Markovian operators, M ≤ −1 to non-Markovian operators, and M ∈ (−1, 1) could be either depending on the boundary conditions. 4 For energy-momentum dynamics in a d-dimensional conformal plasma, gravitational dynamics breaks up into a set of d(d−3) 2 Markovian modes with M = d − 1 (d − 2) non-Markovian modes with M = 1 − d, and a single scalar mode with M = 3 − d. These modes can be characterized by their transformation under the SO(d − 2) little group, corresponding to rotations about a fixed spatial vector (taken to be the direction of spatial momentum). The Markovian modes correspond to the transverse traceless spin-2 polarizations of the energy-momentum tensor (roughly T ij ), the d − 2 non-Markovian mode are the transverse vector polarizations (roughly T 0i ), and the single longitudinal mode is the energy density (roughly T 00 ). The quadratic action for these modes has been derived from Einstein-Hilbert in [8,9] (and for the Einstein-Maxwell action in [10,11]). Higher order graviton vertices are necessary to extract the non-Gaussian influence functionals. Rather than undertake this exercise (which is straightforward but somewhat technical) we identify a simple toy model where we can exhibit all the essential features with the added benefit of analytic tractability.
Our discussion broadly comprises two parts. We first present an abstract model which distills the analysis of gravitational fluctuations. Here we write down an interacting theory of designer scalar fields, analyze the quantization conditions, and give the general rules for computing boundary observables, generalizing [2]. A new element in our presentation is an evaluation of bulk exchange diagrams for real-time Schwinger-Keldysh correlators. In fact, we will argue that a general tree-level Witten diagram can be computed given the knowledge of the ingoing boundary-bulk propagator (similar observations have been made earlier in [15]). We furthermore argue that the general diagram can be evaluated on a single copy of the exterior of the black hole, with the integrand given as a multiple discontinuity [16]. We illustrate this explicitly for low-point correlators, giving expressions for four-point Schwinger-Keldysh correlators in terms of a single master bulk integral.
In the second part, we specify this model to the particular case of the BTZ geometry. This has the advantage that the boundary-bulk propagator may be obtained analytically in closed form in terms of hypergeometric functions. 5 We will furthermore highlight the fact that the designer scalars can have long-lived diffusive quasinormal modes even in this low dimension. Recall that in a two-dimensional thermal CFT energy-momentum dynamics is fixed by the Virasoro symmetry; the system only has left and right moving modes, and thus no hydrodynamic behaviour. While the designer scalars are not components of the conserved current (bulk graviton dynamics in AdS 3 is trivial), the model captures the essential features of the higher dimensional physics in a tractable setting. 6 With benefit of hindsight from our model, we also describe the analytic structure of higher-point Schwinger-Keldysh thermal correlators. We argue that a generic n-point thermal correlator will be meromorphic, with quasinormal or anti-quasinormal poles depending on whether the operator insertion is retarded or advanced. The poles appear in the frequencies corresponding to the external operator insertions, and also in the internal operator exchanges. The analysis is facilitated by use of thermally adapted advanced retarded basis introduced in [18], which judiciously factors out the thermal statistical factors associated with the KMS constraints on real-time correlators. In addition to these quasinormal poles, we also find that exchanged operator frequencies can have Mastubara poles outside their domain of analyticity. These occur at discrete multiples of 2πT , with retarded (advanced) operators potentially supporting such in the lower (upper) half-plane.
The outline of the paper is as follows: In §2 we introduce the general class of models that are of interest and work out the appropriate generalization of the holographic GKPW dictionary. Specifically, in §2.2 we outline the computation of bulk exchange diagrams in the grSK geometry. We also take the opportunity to explain the analytic structure of realtime holographic thermal correlators in § 2.3. In particular, we argue that even higherpoint functions only have quasinormal poles (and potentially Matsubara poles in exchanged momenta).
§3 is devoted to our three-dimensional toy model. We first study designer fields in the BTZ background in §3.1, obtaining analytic expressions for the bulk Green's functions and the boundary two-point correlators §3.2. The analysis of the correlators reveals interesting linear (mode) instability domains in our model ( §3.3). In §3.4 we compute three and fourpoint functions. While the explicit expressions are somewhat complicated, we express them optimally to extract physical lessons for open quantum systems, which are summarized §4. We conclude with a broader discussion and potential generalizations in §5. Some technical details for computing exchange diagrams in the grsK geometry can be found in Appendix A.

Designer fields in grSK spacetime
The problem we want to consider is that of a simple scalar field Ψ coupled to a holographic thermal field theory. Let O a be a set of single trace operators of the holographic theory. We will be interested in both short-lived (Markovian) and long-lived (non-Markovian) operators. For the present, we are going to elide over the tensor indices of the operator. The reader might find it helpful to view O a as specific polarizations of a single tensor operator. As discussed in §1 this is the case for the energy-momentum tensor, with the different components having Markovian, or non-Markovian characteristics.
The combined dynamics is specified by the action of the schematic form We have chosen to represent the coupling between the holographic theory and our system field Ψ using a functional α a [Ψ] to keep track of the dependence on tensor components, polarizations, etc. 7 Since the coupling is directly to the gauge invariant operators, it follows that the functionals α a [Ψ] are simply sources for O a . Therefore, insofar as the effective field theory of Ψ is concerned, we first need the data of the real-time correlation functions of O a , which as noted above are computed holographically in the grSK geometry. This will allow us to write down the leading order (in amplitudes) the effective stochastic theory for Ψ, cf., [2]. The grSK spacetime is a two-sheeted geometry characterized by a complex line element. In the conventions of [2], a stationary asymptotically AdS black hole geometry in ingoing Eddington-Finkelstein coordinates is extended to the complex domain. Consider parameterized by the mock tortoise coordinate ζ. Assuming the emblackening function f (r) to have a simple zero at r = r + , the location of the horizon, ζ is defined by , (2.3) subject to the following boundary conditions at the cut-off surface r = r c This geometry has two asymptotic boundaries which are the backward and forward segments of the boundary Schwinger-Keldysh time contour. The bulk spacetime can be simply characterized by a contour definition for ζ which encircles the horizon at r = r + , cf., Fig. 1.
Re(r) Figure 1: The contour picked by the grSK geometry in the complexified r plane with the ingoing time coordinate v kept fixed. The contour encircles the branch point at the horizon r = r+ counter-clockwise starting from left SK boundary (denoted L) at rc + iε and ending up at the right boundary (denoted R) at rc − iε.
The operators O a are dual to fields ϕ a in the bulk geometry. These fields are not minimally coupled but have a modulation of their gravitational interaction by an auxiliary dilaton χ a (r). We consider a number of these operators having different bulk modulations interacting via a cubic coupling, a bulk three-point vertex for the fields ϕ a . Thus, the action for these fields takes the form The asymptotic behaviour of the dilaton has implications for the boundary conditions, and in certain cases also for the nature of the operator O a . We characterize the asymptotic fall-off of this auxiliary dilaton by an index M a lim r→∞ e χa = r + r We have chosen a minimally coupled scalar field to have index M = d−1. In addition, we have allowed ourselves the freedom for the vertex function to also have a radial modulation through another auxiliary dilaton χ λ (r), which is uncorrelated to the modulation in the kinetic terms. The boundary terms in S bdy include the appropriate boundary terms to ensure the stationarity of the action (discussed below), while the counterterms S ct provide the regularization of the bulk computations.
The motivation for this model comes from the structure of the gravitational fluctuations around a neutral AdS black hole decomposed in gauge invariant variables. As explained in [8,9] (following the original derivation of [14]) the energy-momentum tensor can be decomposed into polarizations based on the little group orthogonal to a chosen spatial momentum k in SO(d − 2) irreps. Transverse traceless tensor polarization are Markovian with M = d − 1, transverse vectors and the longitudinal scalar are non-Markovian with indices M = 1 − d and M = 3 − d, respectively. If we decompose the Einstein-Hilbert action in terms of the gauge invariant combinations of gravitational fluctuations the resulting action truncated to quadratic order will be of the form given in (2.5). The quadratic action for the fields can be found in the aforementioned references, but the cubic vertices have not yet been evaluated directly. We have simply distilled the essential features into the simple model above. 8 In general the dilatonic couplings in the kinetic terms χ a (r) are functions of the radial coordinate -either simple power laws, or some functions constructed from the background metric data. An exception is energy density operator, the scalar polarization of the gravitons, where the modulation also depends on the spatial momentum, see [9,11].
For the present, we will focus on this model and explain the general features for the computation of the real-time correlation functions of O a , which can then be translated into an open stochastic effective action for the field Ψ. Much of this has already been explained in [2] and [8], so we will be brief, and only highlight the salient features. The main new ingredient not present in these works is a discussion of the bulk-bulk propagator, the computation of exchange diagrams in the bulk, and a discussion of the analytic structure of the correlators.

Scalar propagation in grSK geometries
The essential data we need for setting up the computation of real-time thermal correlators from the grSK geometry is the boundary-bulk propagator with ingoing boundary conditions at the horizon. This will suffice to construct the solution of the homogeneous wave equation on the grSK geometry with suitable boundary conditions for the field, and will also determine the bulk-bulk propagator.
A scalar field ϕ with dilaton χ(r) satisfies the free wave equation 9 Using the asymptotic behaviour of the dilaton (2.6) we determine the linearly independent solutions near the AdS boundary where we defined 'conformal dimensions' ∆ to satisfy 10 The fall-offs reduce to the familiar expressions for minimally coupled fields whence M = d−1.
Note that with our choice for the dilatonic modulation (2.6) we effectively can think of the scalars as propagating on an AdS spacetime with effective dimension d eff = 1 + M (similar features were observed earlier in a different context in [20]). We now have to decide how to quantize the fields given these fall-offs. As such, we require the masses to satisfy m 2 ≥ − (M+1) 2 The possible boundary conditions now depend on choice of M and ∆. We do not want to enforce that the ∆ satisfy the unitarity bound relevant for regular conformal dimensions. Since the operator O is to be viewed as a component of a gauge invariant operator satisfying unitarity, it by itself does not need to satisfy the usual restrictions. So apart from requiring ∆ ∈ R we won't prejudice the model with further restrictions. Later, in our two-dimensional example we will see a constraint from mode stability which cuts-off a part of the (M, ∆) space.
To proceed, we also note that conjugate momentum π behaves as (2.10) We can therefore postulate the following boundary conditions for computing the generating functional of correlators: • For m 2 = 0 we have the modes falling off as a constant and r −1−M , respectively. This singles out M = −1 as the separatrix between two possibilities. For M > −1 and M < −1 we impose Dirichlet and Neumann boundary conditions, respectively. In the window M ∈ [−1, 1] we are free to impose either (or even mixed boundary conditions) based on the choice of boundary terms. This was the situation encountered for the designer scalars in higher dimensions [8,9] (and perhaps most physical).
• When m 2 = 0 the fall-offs are separated across the locus ∆ = 1+M 2 , irrespective of the magnitude and sign of M a . We should impose Dirichlet boundary conditions when ∆ ≥ 1+M 2 , identifying the mode falling off as r ∆−1−Ma as the non-normalizable mode. Following the usual rules of quantization in AdS/CFT we will therefore define the source of the operator O in this case to be (2.11) we should instead impose Neumann boundary conditions (2.12) • We need to include suitable boundary terms in order to ensure that the variational principle is upheld. For the Dirichlet boundary condition S bdy,D = 0 since the variation of the bulk action gives a boundary term proportional to π δϕ. On the other hand, to impose the Neumann boundary condition we need to include a boundary term S bdy,N = −´d 2 x π ϕ.
The Markovian or non-Markovian nature of the operator O is not a-priori dictated from the boundary conditions alone. We define an operator to be Markovian if the boundary retarded Green's function is analytic in Fourier domain at low momenta and frequencies.
Otherwise, it will be designated to be non-Markovian. In particular, non-Markovian operators will have thermal correlators which have quasinormal poles with dispersion relation ω(k) → 0 as k → 0.
Markovian operators encountered in the literature are dual to massless fields (m 2 = 0) with M > −1 quantized with Dirichlet boundary conditions, while non-Markovian operators correspond to fields with M < 1 (and m 2 = 0) quantized with Neumann boundary conditions [8,9]. However, is worth emphasizing that the choice of boundary conditions is independent of whether the operator has long-lived or short-lived characteristics. Indeed, we shall see examples of this below, for once we have a mass term, there is a possibility of encountering non-Markovian behaviour with either Dirichlet or Neumann boundary conditions.

Boundary-bulk propagator:
With the boundary conditions specified, we can now give the prescription for computing correlation functions. The only piece of data we need is the ingoing boundary-bulk propagator, G in (ζ, v, x). It will be convenient to work in Fourier domain where 11 We have given the definition for Dirichlet boundary conditions, which can straightforwardly be generalized to the Neumann case. 12 It suffices to determine the ingoing Green's function G in on one of the sheets of the grSK geometry. It is regular as ζ jumps across the sheets for it satisfies (2.14) Given the ingoing boundary-bulk Green's function we can obtain the outgoing Green's function by conjugation [2] G out (ζ, k) = e −βω ζ G in (ζ,k) , (2.15) wherek µ is the d-momentum with frequency reversed,k µ ≡ (−ω, k). The exponential factor e βωζ , which arises as the monodromy around the horizon along the grSK contour, acts as a radial Boltzmann weight. We will adapt this terminology in the sequel. It is also useful to record the time-reversed propagator (2.16) 11 Fourier conventions: We define the d-momentum k µ = (ω, k), but refrain from indicating the Lorentz index. Spatial momentum magnitude will be denoted by |k|. The frequency reversed d-momentum is singled out by a bar decoration:k µ = (−ω, k). Finally, the Fourier transforms is defined as We also abbreviate the momentum integrals by´k as indicated above. 12 For the Neumann boundary condition we would require the source, determined by the conjugate momentum, to be suitably normalized at the boundary and regular at the horizon.
The general solution to the homogeneous wave equation (2.7) with sources J R and J L prescribed on the two boundaries (R and L, respectively) of the grSK geometry is then Here n B is the Bose-Einstein statistical factor The combination of sources appearing above can be used to define a variant of the retardedadvanced basis introduced in [18] As noted in §1 an advantage of this basis is that it factorizes out (for the external operators) the statistical factors n B which arise from the KMS conditions on thermal correlators. It then follows that the solution for the scalar field on the grSK geometry is given as (2.20) We will use this form of the solution in what follows since the Schwinger-Keldysh and KMS conditions imply that any correlator with all operators solely of either the P or F-type vanishes. The only non-vanishing correlators are of the mixed type [18]. From the boundary-bulk ingoing propagator we can obtain the boundary retarded Green's function K(k), using the fact that it is given by the asymptotic value of the field and the conjugate momentum, viz., K(k) = lim r→∞ G in (ζ, k) π(ζ, k). This was justified in the grSK geometry in [2].

Bulk-bulk propagator:
The other piece of data we require is the bulk-bulk propagator G bb (X, X ), where X = {ζ, v, x} for brevity. This satisfies the wave equation with a delta function source One can obtain this Green's function using the variation of parameters trick, which exploits the solutions of the homogeneous wave equation (2.7). To write this efficiently, let us introduce a linear combination of boundary-bulk Green's functions that are normalizable at one or the other boundary of the grSK geometry. Denoting these as G R (ζ, k) and G L (ζ, k), respectively, we have The function G R has a source on the right boundary (ζ = 1), while G L has a source on the left boundary (ζ = 0), and they are respectively normalizable at the other end viz., The bulk-bulk propagator can be seen to be given by (2.24) Here Θ(ζ − ζ ) is a contour ordered step function along the contour depicted in Fig. 1. The prefactor can be obtained from the Wronskian of the two linearly independent solutions to the homogeneous equation which we have taken to be the left and right normalizable boundarybulk propagators. We have separated this out into a normalization factor N(k) and a piece that depends on the radial position of the delta function source. Of these, the e βω ζ factor will be crucial -it ensures that the observables satisfy the Schwinger-Keldysh and KMS constraints. We can furthermore show that the normalization factor N(k) can be determined in terms of the boundary retarded Green's function K(k) as . (2.25) To derive these statements, first, recall that the overall normalization of the bulk-bulk propagator is given in terms of the Wronskian between G R and G L . This in turn is related to the Wronskian between G in and G out through Wr(G R , G L ) = (n B + 1) Wr(G in , G out ) . (2.26) The latter can be computed as follows: To obtain this we used instead of the radial derivative, the derivative operators D ± ζ = ∂ ζ ± βω 2 , introduced in [2]. They importantly satisfy the conjugation e βωζ D + ζ e −βωζ = D − ζ , and also define the conjugate momentum in the non-orthonormal ingoing coordinate basis. The final step is realizing that the product of the ingoing propagator and the conjugate momentum π is the boundary Green's function K(k). This justifies our claims for the normalization of the bulk-bulk Green's function.
The bulk-bulk Green's function (2.24) is normalizable on both ends of the grSK geometry. For fixed ζ , taking ζ → 0, picks up the G R (ζ) G L (ζ ) term, which vanishes by (2.23). In the opposite limit ζ → 1 we pick up G L (ζ) G R (ζ ) which is again normalizable. Using the shorthand notation we can express the bulk-bulk propagator in a compact form as This completes the basic data necessary for computing correlation functions. Note that we only need to obtain the boundary-bulk propagator -all the remaining Green's functions can be expressed in terms of it. This is indeed expected; for instance [5] argued for an analogous structure in the thermofield double state. Their argument ought to apply to the Schwinger-Keldysh geometry. Likewise, we should highlight the fact that the bulk-bulk Schwinger-Keldysh propagator for Schwarzschild-AdS black holes was written down in [15] (see also [21]). Their result satisfies the inhomogeneous equation, but appears to be missing a radial Boltzmann factor, e βω ζ . This factor, as we argue below, is crucial for ensuring that the boundary thermal correlators obey the Schwinger-Keldysh and KMS constraints. 13

Witten diagrams on the grSK geometry
We now have all the ingredients necessary for computing a general n-point function with boundary Schwinger-Keldysh ordering. The recipe for computing contact diagrams has been described earlier in [2]. We will supplement it here with the rules for computing the bulk exchange diagram. Additionally, we will argue that in neither case is there any contribution localized on the horizon (modulo an assumption about vertex factors which we explain below).
The general rule for computing diagrams is to start with the contour integral over the mock tortoise coordinate ζ and convert it to an integral over a single sheet of the spacetime outside the horizon. In doing so, we have to account for the radial Boltzmann factors, e βωζ , arising from the outgoing and bulk-bulk propagator (including the piece originating from the Wronskian). Modulo these pieces, the rest of the integrand can be written in terms of G in (ζ), which is periodic (2.14).
Since the contour in Fig. 1 encircles the branch cut emanating from the horizon we typically are computing the integral of the discontinuity of a function evaluated on the two sheets. For instance, given a integrand L(ζ), which is regular on the horizon, it is easy to show that (nb: dζ Let us suppose for simplicity that we have an n-point self-interaction of a single field, say ϕ n . The bulk action contains a term of the form Expanding out the fields ϕ using the solution (2.20) we can immediately write down an integral expression for influence functional I F···FP···P , the term with say p F sources, J F and (n−p) P sources, J P in the retarded-advanced basis. The integrand is a particular combination of the ingoing propagators and radial Boltzmann factors of the mock tortoise coordinate.
Assembling the pieces we see that contour integral simply picks up a discontinuity, which owing to the periodicity of the ingoing Green's function (2.14) leads to an expression for the influence functional as a single-sheeted integral of the form [2] 14 (2.32) Note that the result depends on the analytic structure of the interaction vertex factor e χn(r) . In writing (2.32) we assumed that this factor is regular the horizon. This is, for instance, the case for the minimally coupled fields, with non-derivative polynomial interactions. Vertex factors which have singularities at the horizon will lead to additional localized contributions. This, in fact, does occur. For transverse fluctuations of a Nambu-Goto string probing the black hole such a vertex exists [19]. One can argue that this extends to derivative interactions, which is the case for conserved currents. In addition, in those cases, it is also possible for there to be special kinematic loci where we find singular vertex functions, potentially along the integration contour. For simplicity, since we will assume the vertex functions to be regular in our analysis, an assumption, which is reasonable to make for non-derivative polynomial interactions. The reader can find further commentary on this issue in §5 where we indicate some generalizations.
Let us next turn to exchange diagrams. The integrand in this case has not only the boundary-bulk propagators, but also the bulk-bulk propagator (2.24). If we expand out the latter, we can assemble the integrand as a product of ingoing boundary-bulk Green's functions (with some reversed frequencies) and radial Boltzmann factors of e βωζ . We additionally have a product of contour-ordered step functions with various radial orderings, i.e., terms of the form Θ(ζ − ζ ) Θ(ζ − ζ ). We decompose the calculation of the integral into a sum of terms, each with a single string of contour-order step function product.
The number of these summands depends on the number of bulk-bulk propagators. A diagram with bulk-bulk propagators has + 1 distinguished radial interaction vertices which need to be ordered. Each bulk-bulk propagator has a binary ordering choice, so altogether we have 2 terms each with a product of ( − 1) contour-ordered step functions arising when we expand out the integrand.
For instance, a single bulk exchange diagram, which has one bulk-bulk propagator has one contour-ordered step function, splits into a sum of two integrals, viz., where the functions F 1 and F 2 are products of the propagators, radial Boltzmann factors, and vertex functions. A two-exchange diagram will analogously result in four integrals as explained in Appendix A.2. Each of these contour integrals needs to be evaluated by respecting the ordering specified. We use the contour-ordered step functions to select the relative ordering and convert the contour integral to a single copy integral. For a given contour-ordered step function there are a-priori three sets of contributions: both the vertices on the L segment, both in the R segment, and one where there lie on opposite segments. In the latter case the two radial integrals run from the boundary to the horizon without constraint. However, when both are on the same sheet, one of the radial integrals is constrained by the other. To disentangle this, we adopt the following convention: • Θ(ζ − ζ ) is the contour ordered step function, with the flow dictated from the boundary SK contour, so out from Re(ζ) = 0 towards Re(ζ) = 1 encircling the cut.
• θ(r − r ) is the ordinary step function defined on a single sheet. We express the radial integrals starting at the boundary and running to the horizon, to maintain consistency with the contour direction. Hence, θ(ζ − ζ ) = 1 when ζ c < ζ < ζ and analogously for θ(ζ − ζ).
With this convention, we can convert the contour-ordering into radial ordering in (2.33). We are furthermore free to use the step function identities and convert all the constrained integrals into unconstrained ones. For instance, the reader can easily verify that the two integrals in (2.33) may be decomposed as The important fact to note is that the integrand is picking up appropriate discontinuities coming from the two sets of contour integrals folded down to a single-sheeted integral. We have organized the latter by the relative radial ordering of the vertices. This process can be iterated to compute any exchange diagram (see Appendix A for further details). Additional features of exchange diagrams are also explained in [16]. Schwinger-Keldysh and KMS conditions require that the influence functionals with purely retarded or advanced operator (i.e., purely P or F-type operators) should vanish. As noted in [2], for contact diagrams this trivially follows from the periodicity of the ingoing boundarybulk propagator (2.14) (upon using momentum conservation for the case of all P sources). Therefore, the corresponding contribution from exchange diagrams also should vanish identically. This should follow again from the periodicity property and not the details of the boundary-bulk propagator. This can only happen if the combinations multiplying the radial step function combinations in the single-sheeted integrals are identically zero. We have checked this to be the case for the four and five point functions with one and two exchanges, respectively, see Appendix A. This, in fact, suffices for any integrand with one or two exchanges since the additional pieces are simply the boundary-bulk propagators (using momentum conservation to remove the exponential factor from the P sources), which are themselves periodic. It should be possible to generalize this argument to prove that the result holds for any number of exchanges, and also for bulk loops, which we leave for the future.
One can also give a unified presentation of the non-vanishing Schwinger-Keldysh correlators. For instance, the single-exchange four-point functions involving four different external The integrand is expressed solely in terms of the ingoing boundary-bulk propagator using (2.15), (2.16), and (2.24). The two factors of G e in originate from the bulk-bulk propagator and are kept distinct by the use of the different frequency labels. The superscript species label indicates the dependence on the index M and the dimension ∆ of the external and internal operators, while the explicit e χ(ζ) factors allow for potential asymmetry of the vertices. The radial Boltzmann weight factors e ω 7 ζ and e ω 8 ζ originate from either outgoing propagators or the Wronskian factor of the bulk-bulk propagator. In physical examples, k 5 , k 6 , ω 7 , ω 8 are functions of the external momenta, but it is helpful to distinguish them for clarity.
We give here the final expressions for the four-point functions in terms of the master integral defined in (2.36). To keep the expressions compact we will write them in terms of the exchanged momenta k = k 3 + k 4 = −(k 1 + k 2 ) whose frequency is ω.
Assuming the external operators to be all distinct, there are a-priori three correlators we should consider three orderings: FFFP, FFPP, and FPFP. All others can be obtained by exchanges F ↔ P . First we note that the FFFP influence function takes the form (2.37) The FFPP correlator is similarly simple, and can be shown to be Finally, the FPFP correlator, which happens to be the most involved, turns out to be (2. 39) In writing the expressions we have repeatedly used momentum conservation for simplification. By assuming that the external operators are all distinct, and allowing the vertices to be non-identical, we have avoided having to sum over different channels. In the case of identical external operators we should add contributions form different channels as appropriate.
The arguments ω 7 and ω 8 of the master integral (2.36) which enter into the expressions for the four-point correlators (2.37)-(2.39) are not generic, but are constrained to satisfy This fact will be helpful when we attempt an explicit evaluation of the integral for the twodimensional model we introduce.

Analytic structure of the correlators
We are now in a position to understand the general features of thermal real-time correlators computed using holography. We delineate the analytic structure of the higher-point functions computed using contact and exchange diagrams. We spell out the assumptions we are making to deduce these results as we go along.
Consider the ingoing boundary-bulk propagator G in , which obeys (2.13). In particular, the source has been normalized to unity. Black hole quasinormal modes, on the other hand, are defined to be normalizable modes that are ingoing at the horizon. The ingoing propagator therefore ought to have poles at the quasinormal frequencies. Based on this intuition it will be useful to write down a factorized form for the boundary-bulk propagator (2.41) The function K(k) is the retarded thermal two-point function on the boundary. K(k) is meromorphic, its poles are the quasinormal modes, which for a sensible thermal system, reside in the lower half of the complex frequency plane. The function G in (ζ, k) is generically regular. It has zeros in its non-normalizable part to compensate for the K we factored out. This is because we have chosen to normalize the source to unity on the boundary. Its normalizable part is clearly analytic, as it actually corresponds to the boundary correlator.
When the boundary-bulk propagator is used to compute the boundary two-point function, one ends up evaluating the difference of the product of the field and its conjugate momentum, ϕπ, at the two asymptotic boundaries of the grSK geometry. Using the asymptotic behaviour (2.13) and (2.10) one then derives the retarded boundary correlator K(k). This typically takes the form of a rational function. While this is guaranteed by virtue of meromorphicity, in examples one finds that the structure to be more specific. The numerator is a function G(k) which depends on the conformal dimension, and the denominator ends up being the same function, but now evaluated as a function of the shadow dimension. The Schwinger-Keldysh structure is completely captured by the fact that only the J F J d combination of the source term is present in the product, cf., (3.14). For minimally coupled scalar fields, expressions for K(k) have been known in the AdS 3 /CFT 2 context for a long while. They were originally discussed in [22] by analytically continuing CFT results, but also later obtained using realtime methods in [4]. 15 Likewise, the result for K(k) is now also known in four dimensions [23] using connection between the minimally coupled scalar wave equation in Schwarzschild-AdS 5 and differential equations satisfied by supersymmetric instanton partition functions.
This data is sufficient to deduce the analytic structure of the higher-point correlation functions. As discussed in §2.2, we will assume here that the vertex functions e χ λ are nonsingular along the ray Re(r) ∈ (r + , ∞) running from the horizon to the boundary.
Contact diagrams: Consider first bulk contact interactions; we can use (2.32) to deduce the analytic structure of the influence functionals. The integrand is a product of ingoing and outgoing propagators. Using the decomposition (2.41), ingoing propagators have quasinormal poles and outgoing propagators have anti-quasinormal poles (due to frequency reversal). This implies we can rewrite (2.32) as (2.42) We have factored out all the non-analytic pieces in the first line. The integrand of I in the second line is a product of analytic functions G in (k). As long as the radial integral converges in the limit r c → ∞, and at the lower limit r + (the latter is at times more constraining), we expect no additional non-trivial singularities. By this we mean that we should not encounter singularities that depend on the dynamical data, viz., momenta and dimensions of the external operators.
It is however possible that owing to the presence of the radial Boltzmann weights we end up with a function that is cognizant of the statistical factor n B . First, we note that Given a particular sequence of F or P operators, by causality the correlator for the corresponding ω i should be analytic in the upper or lower half-plane, respectively. In (2.42) the Boltzmann factors involve a sum of advanced frequencies, and thus any contribution proportional to them should be analytic in the lower-half ω j planes, with j = p + 1, · · · , n plane. Should the radial integral produce a non-analytic term, it has to, for consistency respect this. We indeed find that the radial integrals generically produces a factor of Γ 1 + iβ 2π n j=p+1 ω j . This has the correct analytic structure -its poles are at Matsubara frequencies of ω j in the upper half-plane.
However, the presence of such Matsubara poles in the external operator kinematics is forbidden by the fact that the FP-basis factors out the KMS constraints. Indeed, this was the primary motivation for the introduction of the FP basis in [24]. They in particular, argued that the Schwinger-Keldysh correlators can be obtained from a suitable n-point spectral function, dressed up with Boltzmann factors to obtain the desired correlator specified by F and P labels.
In contact diagram, these Gamma factors with Matsubara poles are rendered harmless by the overall statistical factor out front in (2.42), whose zeros nullify the poles. This, in fact, serves as a nice cross-check for the vanishing of the all F or all P correlator. In summary, the contact diagrams indicate that the n-point correlator only has quasinormal poles corresponding to the retarded (F) operators, and anti-quasinormal poles corresponding to the advanced (P) operators.
Exchange diagrams: These ought to behave similarly, but here we need to account for the analog of unitarity cuts arising from the poles of the bulk-bulk propagator. Let us first use the factorized form to rewrite the bulk-bulk propagator: Then note that the normalization factor N(k), which is given by (2.25) has a factor of K(k)−K(k) in the denominator. It also has a set of Matsubara zeros owing to the Boltzmann factor n B + 1. This implies that N(k) has both quasinormal and anti-quasinormal poles -this follows since G in and G out furnish a basis of solutions. 16 For purposes of ascertaining the analytic structure, we can replace (K(k) − K(k)) −1 by (K(k) K(k)) −1 . With this understanding, it then follows that one of the step-function terms in G bb takes the form (2.45) Therefore, the only terms from the bulk-bulk propagator which have singularities are when both factors are ingoing or outgoing. The former has quasinormal poles, while the latter has anti-quasinormal poles.
It is easy to see what this property of the bulk-bulk propagator implies for an arbitrary tree level exchange. To do so, let us generalize the master integral by introducing a primitive integrand for cubic bulk vertices 17 (2.46) The function B is now analytic, with the poles factored out. A general tree level exchange requires the evaluation of a recursive integral of the form 1 , k 2 , k This integral, should it converge, is analytic, even with some external or internal frequencies reversed. All the factors of K(k) which account for the quasinormal poles, neatly factor out. We then learn that a tree-level exchange has 16 Equivalently, we could have used the fact that G L and G R have both quasinormal and anti-quasinormal poles, which is necessary since neither set by itself is a complete basis of mode solutions [25]. 17 The generalization to higher order bulk interactions is straightforward; we simply upgrade the primitive to include as many factors of the boundary-bulk propagator as the valency of the vertex.
• Quasinormal (or anti-quasinormal) poles from the factors K(k) multiplying the external boundary-bulk propagators.
• A set of quasinormal and anti-quasinormal poles for each of the internal bulk-bulk propagators, arising from two of the terms in (2.45). This occurs when successive factors of B in the integral (2.46) have their third momentum argument coincide. The poles are in the lower half-plane when both the propagators from the bulk-bulk Green's function are ingoing, and in the upper half-plane when they are outgoing.
• Potentially additional Matsubara poles in the exchanged frequencies. These are not forbidden by the spectral analysis of [18], but have to respect the causality properties of the correlator. If the exchanged frequency arises from all F (or all P) operators, then we could have Matsubara poles in the lower (or upper) half-planes. If both F and P type operators are involved, then the poles can occur in the entire frequency plane.
We shall verify these statements explicitly below when we compute the correlators in our two-dimensional example. The discussion above does not rely on specifics of the model, and only rests on the structure of the Green's functions and bulk interaction vertices.

A low dimensional toy model
The discussion above has been quite general and serves to highlight the fact that the computation of real-time thermal correlation functions from holography is in principle on as firm a footing as computing vacuum correlators. It would however be useful to have an example where we can compute the correlation functions explicitly. In dimensions d > 2, there are no analytic expressions available for the essential ingredient, the boundary-bulk propagator (for asymptotically AdS d+1 black holes). However, the BTZ black hole in d = 2 provides a setting which is analytically tractable. This was already exploited in [2] to compute the two and three point influence functionals for a minimally coupled massive scalar field. Generalization to the designer fields is straightforward, and as we shall see below, allows us to even mimic the behaviour of non-Markovian dynamics observed for conserved currents in the higher dimensional examples. In the rest of the paper we will use this simple toy model to illustrate the general principles.

The set-up
We aim to study a simple model of designer scalar fields in the BTZ geometry. We will work with the action (2.5) for the fields ϕ a , characterized by indices M a and dimension ∆ a , respectively. As described in §2.1 we treat the ∆ a as a proxy for the masses of the bulk fields. While they do not necessarily indicate the conformal dimension of a boundary operator in the dual 2d CFT, we will find it convenient to refer to them as such.
The line element for the BTZ black hole in ingoing coordinates is given by (3.1) In the second line we have introduced the inverse radial coordinate z ≡ r + r . The solution has a Hawking temperature T = r + 2π . The mock tortoise coordinate can be determined by integrating (2.3), and we find The dynamics of a designer field ϕ with index M and dimension ∆ is determined by the equation of motion The auxiliary dilaton is taken to be a power law e χ = z 1−M . We have written the equation in terms of time-reversal covariant derivative D + , cf., [8], and introduced the dimensionless frequency and momenta w and q, respectively. These are defined to be The solution to the wave equation (3.3) can be readily obtained in terms of hypergeometric functions. Imposing regularity at the horizon and a unit normalized source on the boundary for a field ϕ a obeying Dirichlet boundary conditions, we find the ingoing boundarybulk Green's function to be 18 This solution can be uplifted to the full grSK geometry with the aid of (3.2).
In presenting the solution, we introduced the combination p ± , which are defined as For a minimally coupled scalar field, which has M = 1 (in d = 2) these combinations reduce to the light-cone momenta. For non-minimally coupled designer scalars, however, we notice that the dependence on the spatial momentum is modified to a surdic form. This observation will be of importance in modeling non-Markovian dynamics with our model. For later reference let us also record the symbol for the frequency reversed version of p ± , which we denote with a bar decoration (as for the 3-momentum k) The final piece of data we need is the bulk-bulk propagator which we can express in terms of the ingoing boundary-bulk propagator following (2.24). We find with the normalization determined to be . (3.9) We will confirm later that this agrees with the form quoted in (2.25).
The results above are derived for the case where the field ϕ M is quantized with Dirichlet boundary condition. If we were to instead quantize it using Neumann boundary condition (depending on the relative values of M and ∆ as described in §2.1), the only change is that we replace the conformal dimension ∆ with the appropriate generalization of the shadow dimension∆ = 1 + M − ∆ . In particular, the ingoing boundary-bulk Green's function with Neumann boundary condition is obtained from (3.5) by swapping ∆ and∆. We therefore for the most part will focus on the Dirichlet choice, and use the replacement rule ∆ ↔ 1 + M − ∆ to infer the answer for the Neumann case. The answers obtained here for the dilaton modulated scalar in the BTZ geometry, curiously enough, have also been realized earlier. In [17] a three-dimensional theory with broken spatial translations was studied. The authors studied an Einstein-scalar system, with two scalars having profiles that preserve the spatial symmetries in the bulk, while having nonvanishing, symmetry breaking, sources along the boundary. The dual is planar black hole with non-degenerate horizon. Unlike the Schwarzschild-AdS 4 geometry, the presence of nonvanishing scalar profiles, causes the metric to fall-off more slowly, with the emblackening function being of the BTZ form. 19 If we study probe fields in this geometry that are insensitive to the broken translations (i.e., those that do not couple to the scalar fields sourcing the geometry) then the probe equations will reduce to (3.3) for some M. [17] analyzed probe Maxwell dynamics in this spacetime, which indeed behaves as we have seen above -the wavefunctions are hypergeometric, and the momentum dependence is of the surdic form in p ± . The paper also examined certain components of the stress tensor correlators, but here we are mildly confused about the harmonic plane wave decomposition in a system with broken translations. This model later used to illustrate the convergence of the hydrodynamic expansion of U (1) current correlators in [27].

Gaussian influence functionals
Armed with the boundary-bulk propagators we can evaluate the on-shell action at quadratic order to obtain the real-time two-point functions of the dual operator O.
First, consider the case where we choose Dirichlet boundary conditions for the field ϕ. We find that the on-shell action evaluates to a boundary contribution (3.11) The conjugate momentum π has an asymptotic behavior (3.12) The kernel K M,∆ (k), which is the boundary retarded correlator, is given as the ratio of Gamma 19 The bulk metric has f (r) = 1 − r 2 + r 2 in d = 3. Therefore, technically, the geometry is not asymptotically locally AdS4 in the standard sense, since the fall-offs are too slow. However, they can be viewed as such, since accounting for scalar counterterms the boundary observables are finite. (3.13) Indeed, the generating function for the boundary correlators reads: This can be expressed in a more familiar form by switching to the average difference basis of operators The quadratic effective action then takes the form making it transparent that the kernel indeed gives the retarded two-point function (it is the coefficient of the source term J a J d ). Now that we have the boundary retarded correlator we can confirm that the normalization of the bulk-bulk propagator quoted in (3.9) has the correct analytic structure. We find from which (2.25) follows.
The thermal two-point function has an interesting structure as the ratio of functions that are characterized by the dimensional ∆ and the shadow dimension∆. While this was known for the minimally coupled fields in the BTZ background (see [2]), it also appears to hold for the thermal correlators in higher dimensions [23].
A similar analysis can be carried out if we quantize the field with Neumann boundary conditions when ∆ < 1+M 2 . In this case we fix the conjugate momentum as the source (2.12). As noted at the end of §3 ingoing boundary-bulk Green's function for this can be obtained from the result for the Dirichlet boundary condition (3.5) by the replacement ∆ ↔ 1+M−∆. This implies that the retarded Green's functions for Dirichlet and Neumann boundary conditions are inverses of each other viz., 21 .
(3.18) 20 The overall normalization factor in K M,∆ should be computed with care. Naively, one expects the conjugate momentum to have a factor of ∆ from the radial derivative. However, this gets converted to the quoted factor when a proper counterterm renormalization is carried out. 21 We are normalizing the sources for the Dirichlet and Neumann boundary conditions with a slight difference to enable this simple inversion of the retarded propagator.

Analytic structure of the retarded correlator and stability
We now analyze the analytic structure of the boundary retarded correlator uncovering a surprising feature. There exists a region of the (M, ∆) parameter space in which the BTZ black hole has a linear instability! This is somewhat curious for one would have thought that being a quotient of AdS 3 the solution is linearly stable towards perturbation by designer scalar probes. The change in the boundary conditions however appears to have a strong effect. A similar phenomenon for minimally coupled scalars with Robin boundary conditions was discovered in [28].

Dirichlet quantization:
The retarded two-point functions have a series of simple poles inherited from the Γ functions in G M (k, ∆). Owing to the rational form of the expression we find a set of at frequencies w D conditioned as follows 22 (3.19) That is, we pick up the only those poles from the Gamma functions appearing in the numerator of (3.13), which are not simultaneously zeros of the denominator. Recalling that our solution is valid for generic values of M and ∆ satisfying ∆ −∆ = Z, we find a set of discrete poles of the two-point retarded correlator for generic q at Taking w and q to be complex-valued, the above defines the quasinormal spectral curve. At special discrete points on this curve, however, the correlator is analytic. For the present case, this happens at the loci (w * D , q * ) characterized by m, n ∈ Z ≥0 w * D = −i (1 + n + m) , q * = ±i M + (∆ + n − m)(M + 1 + ∆ + n − m) . (3.21) We will refer to this set of points on the quasinormal spectral curve as apparent quasinormal modes. 23 Note that this phenomenon is already present for thermal correlators of quasiprimaries in a two-dimensional CFT as discussed in [27,29]. Usually, in discussions of quasinormal modes, one takes the momenta to be real, as one is interested in late-time decay of linearized perturbations (which are superpositions of the quasinormal mode functions). The analytic properties of the retarded two-point function, however, are naturally discussed with both w, q ∈ C. 22 We have placed a subscript D to remind us that the field was quantized with Dirichlet boundary conditions. 23 In the literature this phenomenon is often referred to as the 'pole-skipping' behaviour of the thermal Green's function. We find the terminology confusing when applied to K M,∆ (k) since there is no pole to be skipped -the correlator is analytic at these loci. We will return to this point later in §5.
For our model to be physical, the poles of the retarded Green's functions should lie in the lower half of the complex frequency plane. Since ∆ > M+1 2 with this choice of boundary condition, stability is guaranteed at high momenta q where the quantity inside the square root is positive. For low momenta, we expand out the answer in powers of q 2 , obtaining There are no apparent quasinormal modes for low momenta, as all modes at these discrete loci have w ∼ O(1), i.e., the frequencies are of order the thermal scale. The poles of the retarded correlator lie on the lower half-plane provided Note, in particular, that the minimally coupled case M = 1 is always stable for scalar primary operators satisfying unitary bound ∆ > 0. Interestingly, for M ∈ (−∞, −1) ∪ (3, ∞) there is a choice of ∆ for which the n = 0 mode has diffusive dispersion relation. This is the characteristic feature of non-Markovian behaviour. In this situation the Green's function does not have an analytic low energy expansion, owing to the presence of this ungapped mode. The gray regions correspond to domains where the lowest quasinormal mode is unstable at low momenta. At its boundary indicated by the blue lines, we encounter the existence of diffusive modes, analogous to the hydrodynamic modes encountered for higher dimensional black holes. We have also indicated the massless case, whence ∆ = M+1, which we note is always stable in the Dirichlet regime and has hydrodynamic behaviour in the Neumann regime.
It is curious to see the appearance of non-Markovian operators with Dirichlet boundary conditions. All the examples encountered in higher dimensions in [8][9][10][11] the non-Markovian operators were massless and quantized with Neumann boundary conditions (for computing correlators). However, once we allow for massive fields, we see that the relative independence of the conformal dimension ∆ and the dialtonic index M allows for a more intricate interplay. We illustrate the general features of the boundary conditions, stability and presence of gapless modes in the (M, ∆) parameter space in Fig. 2.
Since we have an exact expression we can continue to work with the correlators even in the presence of an ungapped mode. However, should we be interested in having low energy effective description, then this would not suffice. In that case, as advocated in [8], rather than deriving the generating functional of correlation functions parameterized by the sources {J R , J L }, one computes a Legendre transformed object, the Wilsonian influence functional, parameterized by the operator expectation values {Φ R ,Φ L }.
For the operator O, whose dual field ϕ satisfies Dirichlet boundary conditions, the field valuesΦ are given in terms of the conjugate momentum The Wilsonian influence functional can therefore be computed by quantizing the field with Neumann boundary conditions instead. All the Legendre transform does is invert the Green's function, so the Wilsonian effective action for non-Markovian fields reads

24)
We will for the most part focus on computing the generating function of correlators in what follows. However, once we have the result, we will explain the salient features which we expect for the Wilsonian influence functionals. The analytic control of our model makes it quite straightforward to translate between the two pictures. 24 Neumann quantization: The above discussion can be generalized immediately to the case where we quantize the fields with Neumann boundary conditions. The quadratic generating function of correlators can be immediately written down While the result appears to be similar to the Wilsonian influence functional for operators quantized with Dirichlet boundary condition owing to (3.18), we emphasize that the two expressions (3.24) and (3.25) are qualitatively different. The quasinormal poles for the retarded Green's function (subscript N) are correspondingly located for generic momenta at The apparent quasinormal modes are present at the discrete points (w * N , q * ) Therefore, the model is stable to linear perturbations only when The stability domain can pictorially read-off from Fig. 2.
Once again it is possible to choose ∆ ∈ (−∞, −1)∪(3, ∞) such that the n = 0 quasinormal mode is long-lived with diffusive dispersion. In particular, massless fields with M < −1 always have a non-Markovian mode. This is exactly the class of designer scalars that have been encountered in the higher dimensional black hole context. While the index in those cases is integral, we will refrain from making that choice to respect the genericity condition ∆ −∆ / ∈ Z. 25 In the presence of a gapless quasinormal mode, we can construct the Wilsonian influence functional by Legendre transforming the generating function. For fields with ∆ < 1+M 2 , which quantized with Neumann boundary conditions, this amounts to instead quantizing with renormalized Dirichlet boundary conditions. For the Gaussian effective action the effect is to invert the kernel in (3.25). This has been extensively discussed in the higher dimensional examples, where the only analytic expressions available are in a low energy gradient expansion.

Non-Gaussian influence functionals
We have thus far computed the ingredients that enter the open effective action at the quadratic order. Our next step is to consider the cubic interaction term in (2.5) and compute the non-Gaussian corrections. Note that we have posited that the bulk cubic vertex is itself modulated by a dilaton χ λ . For simplicity, we will take this to be a power law as well, letting the cubic vertex be The cubic vertex leads at leading order to two types of contributions to the effective action: cubic terms which arise from a bulk contact diagrams, and quartic contributions, arising from a bulk exchange diagram. We are going to assume for simplicity that the three fields interacting at the vertex have been quantized with Dirichlet boundary conditions. As noted above the switch to Neumann boundary conditions can be easily achieved by replacing the operator dimension by the shadow dimension.

Cubic influence functional
The contact diagram for a bulk cubic interaction has already been evaluated in [2] for the case of a minimally coupled scalar field in the BTZ geometry (whence M 1,2,3 = 1 and M λ = 0). Allowing for non-trivial dilatons does not change the qualitative nature of the bulk integrals to evaluate, so we will be brief, mainly quoting the result in what follows.
By the Schwinger-Keldysh and KMS conditions the FFF and PPP correlators vanish in the retarded-advanced basis. Furthermore, FFP and PPF are complex conjugates of each other. Thus, we only have one diagram to compute. Using (2.32) we find that we need to the following single-sheet integral The ingoing boundary-bulk Green's function in (3.5) can be conveniently expressed as an integral using the Barnes' representation of the hypergeometric function as The contour C runs parallel to the imaginary axis and is chosen to as to separate the poles at s = −n and s = M+1 2 −∆ − n from those at s = p ± + n. We write each of the three boundary-bulk propagators by this integral representation. The radial z integral turns out to converge provided (3.32) We assume this to hold and complete the radial integral, which gives us a ratio of Gamma functions. The last step is to carry out the three contour integrals from the Barnes' representation. This can be done by closing the contours to the left picking up the poles at s i = −n i and s i = ∆ a i − 1+M i 2 − n i , respectively. Most of the sums can be carried out and the final answer as a single sum over a product of generalized hypergeometric functions as in [2]. Factoring out the boundary Green's functions, we can express the result as follows

33)
We define here a function J δ , which controls the residues at the poles. The notation is as follows: the parameter δ can either be the dimension or the shadow dimension for each of the three external operators. The result is a sum over eight choices indicated by the summation in the second line. The function J δ is itself given as J δ 3c (n) = Accounting for the Schwinger-Keldysh and KMS conditions, which imply the vanishing of FFFF and PPPP correlators, we are left with computing the FFPP, FPFP, and FPPP correlators. We are distinguishing the FFPP and FPFP correlators for the present since we are assuming that the operators are distinct. From the discussion in §2.2 it suffices to quote the master integral (2.36) for our model, as all 4-point orderings can be recovered from it. For completeness, let us record one of the integrands that enters the double bulk-integral (2.35), say for the FFPP ordering . (3.36) In writing this expression we have introduced the 2-momentum k = k 3 + k 4 (and hence w = w 3 + w 4 , and q = q 3 + q 4 , respectively). We can still use momentum conservation to eliminate w (since k 1 + k 2 + k 3 + k 4 = 0) and write the expression as a function of k 1 , k 2 , k 3 alone, but will refrain from doing so. The reader can directly verify that one reproduces the result quoted in (2.38). Thus, we are left with evaluating the master integral (2.36). We use again the Barnes' representation ingoing boundary-bulk propagator (3.31), and write the master integral as the following nested single-sheet integral: where we defined the exponents 26 This general expression turns out to be hard to evaluate. The inner integral can be completed in terms of Appell functions, but that leaves us with a complicated outer integral. However, for the computation of the influence functionals we do not need the expressions in full generality. As noted in §2.2 it suffices to analyze the integral subject to the constraints (2.40). With this assumption, (3.37) simplifies considerably. Focusing on the radial integrals we need to evaluateˆ1 (3.39) 26 We have introduced Ma 5 = Ma 6 = Me and ∆a 5 = ∆a 6 = ∆6 to write the expressions compactly.
These nested integrals can be done in terms of generalized hypergeometric functions. Note that the convergence of the integral demands a constraint on α 1,2 as well as one on the frequencies: 27 We will assume this to hold for the present. With its aid, we can show that the master integral reduces to the following set of contour integrals Note that the parameters α 1 and α 2 depend on the s i which are integrated over. The convergence condition (3.40) is analogous to the one for contact diagrams (3.32) and constrains the external operator dimensions. The constraints on the frequencies (3.41) simply is a statement of causality; one can check that it requires the correlator to follow the analyticity properties dictated by the F and P labels. The contour integrals can be done as before using residue calculus, and the expression written as a sixfold sum. We are able to complete two of the sums in terms of hypergeometric functions and record here the final answer as a fourfold sum over n i with i = 1, · · · 4.
J δ 4ex (n 1 , n 2 , n 3 , n 4 ) . (3.43) We have once again factored out the essential pieces involving the boundary retarded Green's functions. Now we not only have the factors corresponding to the external operators, but also additional contributions which depends on the internal bulk exchange, the k 5 and k 6 terms. We recall that these momenta take values k ork. In the actual influence functionals there is an additional factor of N Me,∆e (k), the normalization factor of the bulk-bulk propagator (3.9), which cancels some poles from these additional factors of K (ensuring that there are no double poles). There are additionally Matsubara poles from the factor Γ(1 − iw 7 ), which we can check always corresponds to the exchanged frequency when it remains uncanceled by the statistical factors.
The function that captures the residue information J δ 4ex is in its turn given by J δ 4ex (n 1 , n 2 , n 3 , n 4 ) (3.44) We have defined In carrying out the sums over n 5 and n 6 (the poles of s 5 and s 6 contour integrals) we have made some choices for pairing them with the remaining summation variables. This completes the derivation of the master integral in terms of which the various 4-point functions are given in (2.37)-(2.39). The reader can check from these expressions that the analytic properties of the correlators delineated at the end of §2.3 are confirmed by these expressions.

Physical lessons for open quantum systems
We have all the necessary ingredients to extract some general lessons for open quantum systems, thanks to the general arguments in §2.3 and the explicit results in our two-dimensional toy model §3. We now examine some specific features for both Markovian and non-Markovian modes. To keep the discussion organized, we first explain features when only one of these types of fields is present, and then turn to the case where they interact.

Markovian self-interactions
The simplest case in question is the self-interaction of a set of fields, all of whose modes are short-lived. This is the case for minimally coupled massive scalars. This was already explored in [2] correlators computed using contact diagrams. The exchange diagrams do not substantially alter the picture.
Consider for the sake of simplicity, a single field φ, whose dual operator O has no longlived modes. We will also consider only the scalar correlators, and take the field to have a cubic vertex, with coupling λ and no vertex function χ λ (r) = 0. The bulk Lagrangian is then To write the generating functions in a compact form, we introduce the following notation: With the aid of (4.2) the generating function for the correlators takes the following form: We have singled out the quadratic part of the generating function to make the dependence on the retarded Green's function of the operator O manifest. In the process, we defined The cubic and quartic influence functionals for our model can be read-off from the previous section. For the quartic case, the sum over channels for the bulk exchange has been performed in writing the above, so the functions I FFFP and I FFPP are suitable combinations of the master integrals (3.43). The exchange of F and P labels, can be achieved by frequency reversal on the corresponding 2-momentum, viz., k →k (for parity even systems). The reader can confirm the analytic structure of the correlators are as predicted in §2.3. In particular, (nb: k = k 3 + k 4 ) The functions R i control the residues -they have dependence on the momentum labels, which are not indicating. What is now manifest is that this is the structure expected from the field theoretic Schwinger-Keldysh and KMS conditions. Real-time diagrammatics, involves only a FP 2-point function and a FFP and PPF 3-point functions. Therefore, the singularities of the 4-point function in the composite momentum k signal the particular channel in which we are decomposing the correlator.
A scalar primary operator of a two-dimensional CFT is a particular example of the type of operator we are considering. In terms of the parameters of the model (4.1), we set M = 1; in this case ∆ is indeed the conformal dimension. The two-point function (3.13) is a wellknown result dating back to [30] and was derived holographically in [4]. The three-point function was first analyzed in [31]; they computed the Fourier transform of the Euclidean (cylinder) correlator and expressed it in terms of Mathieu functions. In [2] it was computed holographically and expressed in the form quoted in (3.33). The four-point function, as far as we are aware, is new and has not been obtained in the literature before. We describe some further applications of our analysis in this particular context in §5.
This data can now be used to construct the effective action of the open quantum field theory. Consider for simplicity, the quantum system, which we use to probe the holographic environment, to be a free scalar field Ψ. We start with the system-environment action (in 2d Minkowski spacetime) (4.6) The combined system is initialized in the product state (|0 0|) Ψ ⊗ (ρ β ) CFT , and evolved with the joint Hamiltonian deduced from the action above. Integrating out the holographic environment, we end up with the effective action for the open Ψ system, which takes the form The factorized kinetic term encodes the bare part of the system action, but the influence functions are obtained by replacing the sources for O by the corresponding retarded and advanced combinations of the field Ψ. One can read off from this action the effective couplings and deduce a classical stochastic model for the open system along the lines described in [2].

Non-Markovian self-interactions
Let us now turn to the situation where we have a single field ψ, whose dual operator P has a long-lived mode. Within the context of our models of §3, such long-lived modes have diffusive dynamics. In higher dimensions, not only do we have diffusive dynamics, but also attenuated phonon modes. For the latter, the dilatonic modulation is more complicated. For simplicity, we will focus on a particularly simple example of a massless field, with Markovianity index M < −1 quantized with Neumann boundary conditions, since this is the situation that arises naturally in higher dimensional examples. The bulk dynamics is characterized by a single cubic coupling, We have indicated the explicit Neumann boundary term necessary to compute the generating function of correlators. Now the asymptotic fall-offs are r 0 and r |M|−1 , with the latter defining the non-normalizable mode corresponding to the source J for P. The operator P has dimension∆ = 0 from the faster fall-off mode. With this data we can write down the generating functional S[ J F , J P ]. This takes the same structural form as for the Markovian case (4.3), which schematically we will write as S[ J F , J P ]. We define this with a different sign from the Dirichlet case to account for the Neumann boundary term. So S gen [ J F , J P ] = S (2) [ J F , J P ] + S (3) [ J F , J P ] + S (4) [ J F , J P ] , (4.9) where now (4.11) The non-Gaussian terms are given as before, cf., (4.3). While this generating functional has all the information one needs, owing to K P (k) having gapless mode, the generating functional is non-local. These poles are explicit in the Gaussian term (4.9), but are also present in the non-Gaussian correlators, as can be discerned from our general discussion §2.3 or directly read-off from (4.5). The origin of this non-local behaviour is easy to discern: in deriving S[ J F , J P ] we have integrated out the low-lying non-Markovian quasinormal mode. The fix, as described in [8], is obvious. We follow the Wilsonian logic, and retain the gapless mode in the low-energy description. One way to implement this is to Legendre transform the generating action to a Wilsonian influence functional parameterized by the expectation values of the operator P. To wit, letting we define 28 The Legendre transform is straightforward to carry-out. At leading order, we recover the expected relation between the sources and fields 29 (4.14) In particular, note that the fields and sources are correctly related by the retarded and advanced sources (F and P, respectively). This relation will get corrected perturbatively (in bulk coupling parameter λ) by the non-Gaussian terms. To get results to quartic order, it suffices to work out the correction from the cubic terms, viz., by solving the system 28 The coupling between sources and operators in the retarded-advanced (FP) basis follows directly from the couplings J R P R − J L P L on the boundary Schwinger-Keldysh contour. 29 In writing the expression, we have used the fact that overall 2-momentum reversal effectively is a frequency (or time) reversal in a system that is parity invariant.
for the sources.
The reader can deduce that the effect of Legendre transform at quadratic and cubic orders simply substitutes the classical relation between sources and fields obtained the Gaussian action (4.14). At quartic order, however, we have in addition a contribution from convolution of 3-point contributions arising from the correction captured by (4.15). We write the result somewhat schematically, but in a suggestive form using a replacement rule, as The last term corrects the quartic influence functional. Its effect is to ensure that the Wilsonian influence functional has a well-behaved low energy expansion.
One can see this as follows: while the cubic influence functions' analytic structure was governed solely by those of the boundary Green's function, the quartic influence functional had additional singularities from the intermediate factorization channels (4.5). First, we note that all the singularities in the external operator positions are removed by the leading part of the solution between the sources and fields (4.14). This is the rationale for writing the result using the replacement rule. To address the second set of singularities, realize that the contribution in δS (4) WIF [Φ F ,Φ P ] is schematically proportional to either I FFP K P I PPF , or I FFP K P I PPF . The specific structure is dictated by the factorization channel of the term under consideration. This term has the same set of singularities at the external operator insertions, but also has singularities in the intermediate channels. These intermediate factorization singularities cancel between the two terms. As with the construction of 1PI actions in QFTs, this is exactly what the Legendre transform is supposed to achieve.
The overall structure can be discerned from how one would organize the Schwinger-Keldysh perturbation theory in the boundary. The retarded Green's function K P (k) acts as the 'kinetic term' for the field variableΦ and one has cubic vertices set by I FFP and I PPF . The Gaussian part has been computed for a conserved U (1) current, and for the energy-momentum tensor in arbitrary dimensions (in both neutral and charged plasmas). The expressions for the non-Gaussian terms can be written down in our toy model (though we have chosen not to explicitly do so).
Let us turn to the open quantum system of a scalar Ψ coupled to a non-Markovian operator P. The dynamics is specified as in (4.6), with the replacement O ← P. This time we write down the effective field theory as the field Ψ coupled to a gapless fieldΦ. The effective action takes the form The coupling in the second line inverts back the Legendre transform. If we carry it out, we end up imprinting the long-lived dynamics of P into non-local terms in the open effective field theory of Ψ alone. However, going by Wilsonian intuition, it is more natural to leave the result in the form given in (4.17), which is manifestly local, and admits a sensible low energy expansion.

Interaction of Markovian and non-Markovian modes
Let us finally turn to the case where we have two fields, one with a Markovian mode (φ) and another with a non-Markovian mode (ψ). We can have two types of cubic interactions between the φ and ψ, so the bulk dynamics can be modeled as (4.18) The fields φ and ψ are also taken to the of the form introduced in §4.1 and §4.2, respectively. We, however, have switched off the self-interactions of the fields for simplicity, to focus on the physics of the interaction between short and long-lived modes. It is once again straightforward to obtain down the generating function for the correlators of the operators O and P, dual to φ and ψ, respectively, in terms of the sources J and J, viz., S gen [J F , J P , J F , J P ], as delineated above. We want to Legendre transform this data to S WIF [J F , J P ,Φ F ,Φ P ] and obtain the influence functional parameterized by the expectation values of the non-Markovian field. We now describe the salient features for the two types of cubic couplings in turn.
Consider first the case where λ 1 = 0, λ 2 = 0. We have a non- Since the dependence on the non-Markovian sources is at most linear at cubic order, we can directly solve for them in terms of the non-Markovian field expectation values and Markovian sources. We obtain In the Wilsonian influence functional, the quadratic and cubic terms are again obtained by substituting (4.14). This also holds for the quartic couplings which mix the two fields, i.e., terms like I PPFF . However, the purely Markovian influence functionals acquire a correction from the cubic pieces in (4.19). Let us again write the Wilsonian influence functional reads S WIF [J F , J P ,Φ F ,Φ P ] = S gen J F , J P , J F → K PΦF , J P → K PΦP + δS (4) WIF . (4.20) Then, for instance, the Markovian four point function I FFFP gets corrected by a term of the form: . (4.21) This is again easy to intuit from the presence of only an FP propagator for the fieldΦ. The corresponding change in the 4-point Markovian correlator I FFPP arises similarly from the FP channel.
A similar exercise can be carried out for the case where the coupling λ 2 = 0. One just has to account for the corrections to the quartic terms arising from the Legendre transform. The cubic couplings correct the I FFPP type correlators. With this understanding it is again easy to write down the effective open quantum description for a field Ψ coupled to the holographic environment. In this case we model the system-environment action as This kind of coupling of a single field to both the Markovian and non-Markovian operators can arise if we couple our system to a conserved current operator of a holographic CFT. For example, the coupling ∂ µ Ψ J µ to a conserved U (1) current J µ is of this form, with the transverse photons being Markovian, and the longitudinal modes being diffusive. Similar statements hold for coupling to the energy-momentum tensor. We can write down following the preceding discussion the effective action for the open system. We integrate out the Markovian field O, treating Ψ as its source, whilst retainingΦ the field parameterizing the expectation value of P. One ends up with

Discussion
We have described a broad class of models for studying open quantum field theories, both with and without long-lived gapless modes. Our construction was broadly inspired by earlier holographic analysis, in particular, the ability to model long-lived modes using designer scalars. While realistic examples have addition features, such as a more complicated radial mass or potential term, the essential point is the overall simplicity afforded by the holographic constructions.
The open effective field theory is governed by real-time thermal correlators of the holographic environment. To compute these, we exploit proposal of [3]. In the process we strengthen and amplify a point made in [2], viz., that the grSK geometry provides the natural background for computing higher point functions using Witten diagrams. In the present work we have explained how to compute exchange diagrams. These, in principle, could be exploited to also compute bulk loop effects. A useful corollary of our analysis is that thermal n-point functions are computed as radial integrals in a single copy of the exterior region of the stationary black hole spacetime. Moreover, the only data necessary is that of the ingoing boundary-bulk propagator. The integrands for the n-point functions are obtained as multiple discontinuity of a function built from them, and a radial extension of the Boltzmann weight, the ubiquitous factors e βωζ , which capture the monodromy picked up as we cross the horizon. These features can also be interpreted in terms of a 'bulk open effective field theory' as will be discussed in [16].
As noted in the main text some of these aspects have been touched upon in the literature earlier. For instance, [15] discussed the computation of exchange diagrams in the black hole background. Their analysis relied on using advanced/retarded propagators in the bulk, but by working directly with our bulk-bulk propagator we have established that the grSK geometry respects all the Schwinger-Keldysh and KMS conditions.
Our analysis here also touches upon several themes that have been discussed in the literature in related contexts. We will outline some lessons and open questions in these contexts, organized thematically, below.
The grSK contour & horizon localization: As we have shown, the computation of bulk contact and exchange diagrams in the grSK geometry, reduces to the computation of a radial integral in the domain r ∈ [r + , ∞). The integrand for the contact diagrams is a single discontinuity of a combination of ingoing boundary-bulk propagators, radial Boltzmann weights, and vertex factors. A similar structure pertains to the exchanges, though now we have to take a multiple discontinuity. The ingoing propagators and e βωζ factors are regular at the horizon. Therefore, as noted below (2.32) and in §2.3, as long as the vertex factors are analytic at the horizon, there is no localized contribution in the grSK geometry.
In general for non-derivative polynomial interaction of fields, it would be non-covariant to have bulk vertices with explicit factors of 1/f . Hence, horizon localized contributions are precluded in such situations. However, as noted in footnote 8 these can occur when we have derivative couplings. This is for instance the case for fluctuations of a probe string studied in the context of Brownian motion in [19]. There localized contributions arose owing to the derivative interactions from the Nambu-Goto action (leading to a double pole at the horizon).
Our interest in such localized contributions stems from [32], who were analyzing 4-point out-of-time-order correlators and encountered such in the eikonal limit. While their analysis is a-priori not directly related to ours, it is interesting to inquire whether the there are situations in the grSK geometry with horizon localized contributions, and what they imply for boundary observables.

Hydrodynamic correlators:
In §3 we pointed out the similarity between the two-dimensional toy model we analyzed, and a three-dimensional model with broken translational symmetry. As noted there [17] analyzed the behaviour of a probe Maxwell field and computed the twopoint function of the boundary global current. Since the bulk Maxwell theory is Gaussian, this is the only observable. However, a variant of the model, an abelian-Higgs system, with a charged scalar field (also a probe), offers interesting opportunities to examine the interaction between fields that different infra-red behaviour. The Maxwell field has a non-Markovian piece in the charge density mode, which interacts now with the charged scalar field, along the lines sketched in §4.3.
Another context where there is a natural interaction between a Markovian and non-Markovian mode, is when a spinless primary field interacts with the energy density operator. In the gravitational dual, this maps to the interaction between the scalar field dual to the spin-0 primary interacts with the scalar polarizations of the gravitons. There are some interesting aspects to this, owing in particular to the fact that the energy density mode has a momentum dependent modulation. We hope to report on this in the near future as it appears to have a useful lesson for horizon localized contributions.

Relation to thermal bootstrap:
As noted at the end of §4.1, an interesting corollary of our analysis is the computation of thermal four-point functions for scalar primaries of a 2d CFT. While their analytic structure is clear from the general discussion of §2.3, it would be useful to decompose the result in terms of thermal blocks on the cylinder and obtain an inversion formula along the lines of [33] for the thermal OPE data. The conformal bootstrap at finite temperature was broadly analyzed in [34]. In the specific case of two-dimensional Euclidean CFTs [35] examined torus conformal blocks. The aim here would be to set up the corresponding problem in real-time, using the momentum space data obtained for the influence functionals. In this context, it is also worth noting that thermal CFTs in momentum space were analyzed in general dimensions in [36], who also obtained the spectral density from the (cylinder) conformal blocks in two-dimensions. 30 Comments on apparent quasinormal modes: Our analysis gives a clear picture of the analytic structure for the boundary correlators. In the C 2 parameterized by ω, |k| we encounter a codimension-1 loci of quasinormal spectral curves (labeled by a non-negative integer n). However, at discrete codimension-2 points, viz., at special kinematics, the correlators are analytic. These we decided to refer to as apparent quasinormal poles. By virtue of the arguments in §2.3, the same behaviour holds for the bulk Green's functions G in , G out , and G bb .
This observation helped us disambiguate some statements and terminology which we have found confusing in the literature, as noted in footnote 23. In attempting to find a relation between the physics of scrambling and many-body chaos, which is encoded in out-of-time-ordered four-point functions, and thermal energy density two-point functions, [38] argued for the phenomenon of 'pole-skipping'. The phrase refers to (potential) poles of the retarded two-point functions, which are skipped because of an accidental zero. The particular locus nevertheless is claimed to control scrambling features of the higher-point correlator.
As one might expect, there is a story from the bulk analysis, which is revelatory. Initially, [39] observed in the wave equation governing energy density correlators, special codimension-2 loci in C 2 , additional modes that are normalizable at infinity and analytic at the horizon. This was point was fully appreciated by [40] who did a careful analysis of the wave equation, but as pointed out in [27,29] the phenomenon is quite generic.
Modes that are normalizable at the boundary do not however register in the boundarybulk propagators. Moreover, no poles are technically skipped in the boundary Green's function (3.20). The function K M,∆ (k) is meromorphic with simple poles at these loci (nb: meromorphic functions are rational functions). So if one was to examine, say the Mittag-Leffler form of the thermal two-point function, one would find no indication of a special values of frequency and momenta. One would simply read off the physical quasinormal modes where the correlator is singular. Consequently, there is no clear notion of pole-skipping if one examines just the retarded two-point functions alone. However, given the original motivation for the connection, this deserves further investigation.
Bulk loops and analytic structure: Our analysis of non-Gaussian influence functionals was restricted to tree level Witten diagrams in the bulk. This captures the leading planar contribution to the correlation functions. Bulk loops, which can in principle be computed within our formalism, will give subleading corrections. It is interesting to ask whether the analytic structure found for the correlations, viz., that they are meromorphic with quasinormal (and anti-quasinormal) poles, is robust to such non-planar corrections. 31 We speculate that in planar perturbation theory the analytic structure is not modified.
To motivate this, consider the vacuum 2-point function, where the bulk 1-loop diagram gives the leading non-planar correction to the anomalous dimension. Focusing for simplicity on primaries of dimension ∆ in a 2d holographic CFT, the thermal 2-point function can be obtained from analytically continuing this result. In Fourier space, the resulting 1-loop answer takes the form quoted in (3.13), albeit with a corrected conformal dimension, ∆ → ∆+O c −1 . Expanding out this in the anomalous dimension, we find the 1-loop answer can be expressed as a product of the tree level result and polygamma functions, which are meromorphic. While suggestive, it remains to confirm whether this expectation is borne out for d > 2. One also ought to be able to address this without recourse to the analytic continuation directly within the grSK formalism.
The vanishing condition is strongly contingent on the presence of the exponential factor of e βω ζ in the bulk-bulk propagator. It is easy to check that without it, the above integrands are non-vanishing. As we noted in the text, this factor is deduced by examining the Wronskian between the properly normalized basis of solutions. We have independently tried to constrain the bulk-bulk propagator directly by demanding the vanishing condition for all F or all P correlators (cf., footnote 13). At the level of four-point functions, we find the constraints do not suffice fix to fix the functional dependence on the source point (ζ above) completely. It however is possible that a similar exercise with a few higher point correlators could suffice. We will not undertake this exercise here, but argue below for consistency checks at the level of two bulk exchange processes.

A.2 Two bulk exchange diagrams
We now turn to the situation where we have two bulk exchanges on the grSK geometry. This is mostly to exemplify the general structure and provide further evidence for the consistency of our identification of the bulk-bulk propagator. For simplicity, we focus on five-point functions computed in a theory with cubic bulk vertices.
With this assumption, topologically, the diagram is concatenation of two bulk-bulk propagators with five boundary-bulk propagators. Generally we have integrands of the form etc., where we have binary choices of G ϕ in and G ϕ out to attach to the sources J F and J P , respectively.

(A.3)
For the vertex positions ζ, ζ , ζ we have 3! orderings. But there are only 2 2 = 4 combinations of contour ordered theta functions arising from multiplying out the bulk-boundary propagators. The coefficient functions F 1 , · · · F 4 denoted above can be viewed as follows. The F 1 and F 4 terms are fully contour ordered (the former will be referred to as ordered, and the latter as anti-ordered), but the integrands F 2 and F 3 are only partially ordered. We use the contour step function identity (2.34) to complete such partial orderings. This gives us the six ordering expected for the 3 vertex integration positions. We then pick out as the following integrands: Our aim is to combine these using the contour ordering and reduce the result to an integral on a single-sheet, in analogy with (2.35).
To do so, we realize that we also have the freedom of placement of the operators in either leg of the bulk SK contour. Since the absolute ordering is fixed, we get only 4 choices for each of the six cases (altogether 24 possibilities). Operationally, imagine fixing a permutation of the three vertices, say the fully ordered ζ > ζ > ζ , and then side the positions consistent with the ordering on each of the sheets of the grSK contour cyclically (i.e., from L to R through the horizon cap). Finally, we decompose these, after placement of the vertices, correctly ordered, into single copy integrands, with the standard step functions. At the end of the day we find for each of the 6 single copy orderings eight possibilities, leading to altogether 48 terms. The count of 8 per single-copy collapsed order is easy to see as a triple discontinuity; each discontinuity/contour integral gives a pair of terms with relative sign.