Open quantum systems and Schwinger-Keldysh holograms

We initiate the study of open quantum field theories using holographic methods. Specifically, we consider a quantum field theory (the system) coupled to a holographic field theory at finite temperature (the environment). We investigate the effects of integrating out the holographic environment with an aim of obtaining an effective dynamics for the resulting open quantum field theory. The influence functionals which enter this open effective action are determined by the real-time (Schwinger-Keldysh) correlation functions of the holographic thermal environment. To evaluate the latter, we exploit recent developments, wherein the semiclassical gravitational Schwinger-Keldysh saddle geometries were identified as complexified black hole spacetimes. We compute real-time correlation functions using holographic methods in these geometries, and argue that they lead to a sensible open effective quantum dynamics for the system in question, a question that hitherto had been left unanswered. In addition to shedding light on open quantum systems coupled to strongly correlated thermal environments, our results also provide a principled computation of Schwinger-Keldysh observables in gravity and holography. In particular, these influence functionals we compute capture both the dissipative physics of black hole quasi- normal modes, as well as that of the fluctuations encoded in outgoing Hawking quanta, and interactions between them. We obtain results for these observables at leading order in a low frequency and momentum expansion in general dimensions, in addition to determining explicit results for two dimensional holographic CFT environments.


Introduction
The study of quantum fields in curved spacetimes, especially in geometries with horizons, such as black holes or cosmological spacetimes, has been an immensely valuable window into the semiclassical aspects of gravitational dynamics. Such investigations have been instrumental not only for understanding the effects of vacuum polarization, particle production, etc., but also have played an important role in the holographic AdS/CFT correspondence. In the latter context the semiclassical gravitational computations give us JHEP07(2020)242 insight into observables e.g., correlation functions, von Neumann entropy, etc., of the dual strongly coupled QFT. One natural set of observables in a quantum system are thermal correlation functions, capturing real-time response to perturbations and the attendant thermal fluctuations. These are especially interesting in strongly correlated systems where analytic techniques to compute response functions are limited. Here, holography provides a valuable avenue: one can extract dynamical response of strongly correlated thermal plasma by performing classical calculations in the dual black hole geometry. This approach has paid rich dividends over the past two decades: ranging from understanding thermalization [1], linear response and hydrodynamics [2], to real-time transport computations [3]. Much of this success owes to the fact that black holes in asymptotically AdS spacetimes are dual to thermal QFTs.
We wish to argue here that black holes in fact provide a playground for the exploration of a much richer set of dynamics, viz., that of an open quantum field theory, where the degrees of freedom of the quantum system are non-trivially entangled with some external environment (or bath) degrees of freedom. The set-up we have in mind is the following: consider a QFT say with one bosonic degree of freedom Ψ(t, x) (for simplicity) which is our system of interest. We will model the environment by another field theory, now with many degrees of freedom X i (t, x). The unitary microscopic theory is of the form: (1.1) The combined system and environment is prepared in some initial state, which we may even take to be factorized between the respective degrees of freedom. Integrating out the environment degrees of freedom X i we end up with a non-unitary evolution of our system. The basic paradigm for such an effective theory was described by Feynman and Vernon [4] who noticed that the natural way to describe the system is in terms of a doubled set of degrees of freedom for the system, together with a non-trivial interaction between them, which they dubbed influence functionals. (1.2) where S IF [Ψ R , Ψ L ] is the aforementioned influence functional, induced onto the system owing to the coupling with the environment. This paradigm is well understood and tested for Gaussian dynamics in quantum mechanics, as exemplified by the Caldeira-Leggett description of quantum Brownian motion [5]. For an overview of developments in the study of open quantum systems see [6][7][8].
One major question in this scenario is to find a set of sufficient conditions for a local effective field theory to emerge (for our system degree of freedom Ψ). This question is quite difficult to address within weakly coupled systems for the following reason: often a local description emerges at time scales longer than the 'environmental memory' time scale τ m . Here, τ m should be thought of as the time that the environment fields take to forget the information of the initial state; it is inversely proportional to the interactions within the environment. Consequently, one has to often wait for a non-perturbatively long JHEP07(2020)242 time for a local description to be valid. This necessarily means that derivation of open quantum field theories is inevitably a non-perturbative question. This explains to some extent why to date there are no simple microscopic models from which a local non-unitary open quantum effective field theory has been systematically derived. In particular, as far as the authors are aware, currently there are no microscopic QFTs from which a local open EFT with interactions can be derived. 1 Our goal in this work is to address this lacuna by using holography.
The set-up we have in mind is semi-holographic, cf., [15]. Say we wish to understand the dynamics of a single bosonic degree of freedom which we continue to call Ψ(x) in d spacetime dimensions. We imagine coupling this to a strongly coupled thermal environment comprising of some intrinsic microscopic degrees of freedom. For concreteness, one can imagine the environment to be the thermal large N , strongly coupled, N = 4 Super Yang-Mills (SYM) theory (gauge group SU(N )) in d = 4, or a large c thermal CFT in d = 2. Following common practice, we will often refer to the system and the environment as the probe and the bath, respectively.
The coupling of the probe/system degree of freedom Ψ to the thermal bath is via a local coupling d d x Ψ(x) O(x)where O(x) is a simple operator in the bath/environment theory. In the aforementioned examples O could be a low lying single trace conformal primary operator. Crucially, the thermal bath theory is assumed to have a holographic dual. This will enable us model the environment by a dual black hole geometry. The influence phases of our system are then encoded in the real-time or Schwinger-Keldysh (SK) correlation functions of the environmental degrees of freedom. The computation of the latter, thanks to the holographic map, is one which can be carried out using classical fields propagating in a classical black hole background. From the holographic standpoint, the problem therefore boils down to developing a formalism for computing real-time observables in black hole backgrounds (or more generally in spacetimes with horizons).
Before commenting on the real-time computation let us first recall a well-known, but remarkable, fact of the Euclidean gravitational path integral. Thermal boundary conditions require that the Euclidean time (t E ) direction be compact with period given by the inverse temperature asymptotically. For the gravitational path integral these boundary conditions pick out the Wick rotated black hole solution as the Gibbons-Hawking saddle point solution [16]. This is unlike any non-gravitational system, where the Euclidean thermal circle is more of a computational aid, the background geometry being non-dynamical. Said differently, gravitational dynamics interplays non-trivially with thermal boundary conditions. Given this solution, one way to pass to a real-time description is to slice open the Euclidean solution at some instant of real-time, say at t = 0, which exposes quite naturally two copies of the asymptotic region at (t E = 0 and t E = β 2 , respectively). This initial data can then be evolved in real Lorentzian time to give the (future half) of the eternal black hole solution. This is the gravitational preparation of the thermofield double state in the JHEP07(2020)242 doubled field theory Hilbert space, which in energy eigenbasis is expressed as where the R(L) refers to the asymptotic boundaries at t E = 0 (t E = β 2 ). Insofar as thermal equilibrium properties are concerned, the thermofield double construction proves ample. One gets to ask questions about correlation functions with operators inserted in either copy of the doubled system. One can therefore view the gravitational path integral as computing the following generating function for asymptotic observers (in asymptotically AdS spacetimes) The fact that we slice open the functional integral midway is what is responsible for the fractionation of the thermal density matrix ρ β .
To obtain real-time response one then has to analytically continue these results to the real-time domain. In the absence of any sources deforming the state away from equilibrium there is in principle no obstacle to carrying out this analytic continuation. However, once one moves to the physics of systems in local equilibrium, or more generally outof-equilibrium, which more pertinently applies to our discussion of computing influence functionals, the thermofield double state proves less useful. In this context the Schwinger-Keldysh formalism provides a cleaner and more natural way for computing real-time observables by keeping manifest causality and unitarity, without relying on the aforesaid analytic continuation. The Schwinger-Keldysh generating function does not fractionate the density matrix, but rather computes the generating function: (

1.5)
This fact is well-known in non-equilibrium QFTs where real-time observables, say linear response captured by viscosity, conductivity, etc., are always computed (using Kubo formulae) from the Schwinger-Keldysh formalism. The pictorial representation of the path integral contour shown in figure 1 provides a quick way to see the difference between the two constructions (cf., [17] for further discussion). Given this fact, it is natural to ask how the Schwinger-Keldysh construction can be adopted to gravitational theories. We will primarily focus on situations where the system is prepared in a thermal state. We wish to know how the gravitational dynamics fills in an asymptotic Schwinger-Keldysh contour shown in figure 1(b). This question has been considered by many authors in the AdS/CFT context over the years. In [18] the first proposal for computing real-time correlation functions was given. These authors posited that one should consider the future half of the domain of outer communication of a Lorentzian black hole spacetime, and impose ingoing boundary conditions on the future horizon to extract causal observables (retarded Green's functions). In addition they also argued for the absence of any boundary contribution from the horizon. This prescription was justified JHEP07(2020)242 shortly thereafter in [19] using the maximal Kruskal extension of the black hole geometry (the logic was to exploit the analytic structure taking inspiration from the Euclidean thermofield double construction). One limitation of this approach was that it was well adapted to the computation of two-point functions, but it left implicit how to obtain higher point functions. Nevertheless, over the years, various authors have attempted to use this prescription for various applications [20][21][22][23][24][25].
In order to explain the subtlety let us remind the reader how Euclidean n-point correlation functions in AdS/CFT are computed via Witten diagrams [26]. One first computes the appropriate bulk-boundary propagators for various external operator insertions. These enable us to 'evolve' the corresponding fields into the bulk with sources set by the boundary conditions. One convolves the bulk fields thus obtained using the bulk interaction vertices which are dictated by gravitational dynamics. Since the fields interact locally, one has to integrate the position of the vertex over the entire Euclidean bulk manifold. Modulo the choice of temporal boundary conditions one expects something similar for the computation of real-time retarded observables. However, it is was a-priori unclear what domain of the bulk geometry one ought to integrate the bulk interaction vertex over, even assuming that the ingoing boundary conditions serve to pick out the appropriate bulk-boundary propagator.
This issue was addressed in [27,28] who gave a more detailed prescription for real-time computations, arguing that one should fill in asymptotic Schwinger-Keldysh contours with piecewise smooth geometries: real-time evolution sections of the boundary contour get filled in with Lorentzian geometries, and imaginary-time segments with Euclidean geometries. These geometries are glued together continuously along codimension-1 spacelike slices. The authors developed a robust holographic renormalization scheme for asymptotically AdS geometries [28]. Furthermore, the ingoing boundary condition for retarded correlation functions was derived quite cleanly using this prescription in [29]. This prescription was employed in the derivation of covariant holographic entanglement entropy proposal [30].
Per se, it is then clear that in order to carry out the computation of real-time correlation function of probe operators in a fixed state, one could simply use the prescription of [27,28].

JHEP07(2020)242
However, if one were to ask questions about dynamically evolving geometries, one realizes that the piecewise smooth geometries pose a potential issue in the presence of horizons. Physically, the question essentially becomes one of coming up with a prescription ensuring that effects of the outgoing Hawking quanta are correctly accounted for (see [22,23] for attempts in this direction). While a complete answer to this question is still unclear, an interesting prescription was recently given by [31] to address this lacunae in the probe limit (see also [32]).
The prescription of [31] postulates that the gravitational dual of the asymptotic Schwinger-Keldysh contour is given by a complex two-sheeted spacetime. This geometry is made of two copies of the black hole exterior smoothly glued together across the future horizon, with a particular monodromy condition. One may view this as the statement that the asymptotic Schwinger-Keldysh contour gets filled in by a complex geometry. This can be motivated by the slicing open the Euclidean Gibbons-Hawking saddle [16] whilst incorporating the Schwinger-Keldysh boundary conditions, as we have attempted to illustrate in figure 2. The authors of [31] already demonstrated the efficacy of their prescription by obtaining the quadratic influence functional in a low frequency and momentum expansion for a bulk scalar field in the probe limit. 2 Subsequently, [33] explored how this prescription may be employed to study nonlinear Langevin dynamics of a single particle degree of freedom in quantum mechanics. Their analysis subjected the prescription to the test of dealing with self-interacting fields in the bulk and demonstrated that it continues to give reasonable answers. Specifically, the authors of [33] modeled the Brownian particle using the holographic construction described in [21,34]. Using the complex holographic geometry of [31] to compute influence functionals, they demonstrated that general expectations of non-linear Langevin dynamics explored in earlier [35,36] was borne out. We will demonstrate in the course of our analysis that the prescription correctly captures real-time observables and can be used in the semi-holographic setting to model open effective field theories.
We have motivated the discussion in terms of invoking the holographic duality to learn about influence functionals in an open quantum system. It is instructive to also keep in mind that these influence functionals for holographic thermal baths in fact encode interesting semiclassical gravitational information about black holes. In the semi-holographic set-up we motivated above, the system degree of freedom Ψ gets imprinted upon by the characteristics of the environment. Since our thermal environment is provided by a black hole, we would effectively be encoding not only the dissipative behaviour of the horizon, which we know to be characterized linearly by quasinormal modes, but also the fluctuations of the horizon. The latter are nothing but the outgoing Hawking radiation. Since the quasinormal modes refer to the physical response of infalling matter we would effectively be capturing, in the non-Gaussian influence functionals, the interaction between infalling matter and the Hawking radiation. 3 2 The authors of [31,32] also addressed the problem of Maxwell gauge field in the probe limit. However, as pointed out in [31] many aspects of how the prescription works for gauge theories are still unclear. 3 The non-Gaussian correlations we compute holographically using the semiclassical gravity approximation are suppressed in the planar expansion by powers of 1/N , as indeed are all higher point functions in a large N (or large central charge) environment.

JHEP07(2020)242
The coupling of AdS black holes to external systems has recently been of active interest in the context of the black hole information paradox [37,38] (cf., [39]). In these examples the external system is treated as a passive reservoir wherein one captures the Hawking radiation. Our discussion applies in this context as well; the external system's observables will faithfully be able to diagnose the interaction of Hawking radiation with infalling matter. The structure we get to probe however is only the leading semiclassical pieces of the interaction. We are considering a single gravitational saddle point configuration, and not including contributions from non-trivial replica saddles which have been important in understanding the purification of the Hawking radiation and the reproduction of the Page curve for AdS black holes [40,41].
In this paper, we will be considering the coupling of our probe/system (modeled by a single bosonic field) to a scalar operator in the holographic thermal system. Most of the technical computations we report will involve computing Schwinger-Keldysh correlation functions of a scalar field in an asymptotically AdS black hole background using the holographic Schwinger-Keldysh geometry of [31]. This work had already considered the computation of the two point function in a low-energy gradient expansion, i.e., perturbatively at low frequencies and momenta. We will extend this to higher point functions, but also show how to get results outside the gradient expansion in two dimensional CFTs.
The outline of the paper is as follows. We will review the gravitational prescription for filling in the Schwinger-Keldysh contour in section 2, illustrating in the process the generalizations to an arbitrary spacetime with a Killing horizon. In section 3 we outline the specifics of the open quantum system we wish to study and describe the holographic thermal environment we are coupling it to. In section 4 we then turn to the task of solving the scalar wave equation (which we generically do at long wavelengths), and describe various propagators of interest that appear in the computation of the Schwinger-Keldysh Witten diagrams. We compute the influence functionals using holographic methods in section 5 and furthermore demonstrate in section 6 that we can use our results to provide a stochastic description of the effective open quantum field theory. We end with a brief discussion in section 7.
Several technical steps are outlined in various appendices. In section A we explain various aspects of the long-wavelength gradient expansion we work in, while section B gives specifics of the Green's function in Schwarzschild-AdS d+1 geometries obtained in this approximation. We review in section C why the standard Witten diagram technique continues to work for computing influence functionals. In section D we describe how non-Gaussian influence functions in 2d CFTs can be computed. Finally, in section E contains details of the divergence structure of bulk Witten diagrams and a desciption of counterterms that enter into the influence functionals.

grSK: the gravitational Schwinger-Keldysh saddle
The Schwinger-Keldysh contour for a thermal state is a complex time path running from t = 0 to t = T and thence to t = 0 − i β as depicted in figure 1(b). Thus we have a contour in a complex time plane for the temporal part of the action at the boundary for a JHEP07(2020)242 Re(r) Figure 3. The complex r plane with the locations of the two boundaries and the horizon marked. The grSK contour is a codimension-1 surface in this plane (drawn at fixed v). The direction of the contour is as indicated counter-clockwise encircling the branch point at the horizon.
holographic field theory. The proposal of [31] is to extend this contour to a codimension-1 hypersurface in the complexified bulk spacetime in gravity. To be specific, let us first introduce this geometry for stationary configurations with a timelike Killing field, such as the planar Schwarzschild-AdS d+1 black hole. We start with the metric written in ingoing Eddington-Finkelstein coordinates, which are regular at the future horizon, viz., The coordinate v is identified with the time coordinate t on the boundary of the spacetime r → ∞, which as we have argued, is to be interpreted as a curve in the complex plane. The idea is to also upgrade the radial coordinate to the complex domain and pick a codimension-1 slice through the resulting complex spacetime. Operationally, we in fact upgrade radial tortoise coordinate to a complex variable, which we refer to as the mock tortoise coordinate, ζ. We define this coordinate by the differential relation: where β = 4π d r h is the inverse temperature of the black hole. The rationale for introducing this coordinate is that ζ picks up a logarithmic branch cut from the integral about the zero of the emblackening function f (r) of the black hole. The choice of normalization is such that the monodromy around this cut is set to unity. The coordinate ζ can be viewed as parameterizing a two-sheeted surface, each of which can be thought of as the bulk extension of the Schwinger-Keldysh contour. On each sheet ζ has an imaginary part running from 0 at the AdS boundary to ∞ at the horizon. In addition it has a real part which differentiates the two sheets and is given by the monodromy around the horizon. By convention we will choose one of the sheets to have vanishing real part, and the other to have unit real part (this is based on our choice of normalization). We will also cut-off the AdS geometry at a radial cut-off r = r c for computational ease. With this choice, we have the two branches on which the mock tortoise coordinate asymptotes to A section of geometry in the mock tortoise complex plane is illustrated in figure 3.

JHEP07(2020)242
The metric then takes the form where we treat r(ζ) using (2.2). We will refer to this geometry as the gravitational Schwinger-Keldysh (grSK) saddle or geometry. One can integrate (2.2) explicitly to find that the mock tortoise coordinate for Schwarzschild-AdS d+1 geometries in terms of a hypergeometric function, viz., where ζ c is chosen to make ζ = 0 at r = r c + i ε. The branch-cut of the hypergeometric function is taken to run from r = r h to ∞. 4 One can motivate the grSK geometry by recalling the Schwinger-Keldysh generating function (1.5) which computes casual response of the thermal state (as well as the fluctuations thereabout). In particular, it requires that we do not fractionate the thermal density matrix. Since the excursion into the complex time domain between the two segments of the contour on the boundary in figure 1(b) is infinitesimal, the gravity saddle point should respect this separation. One way to achieve this is to prepare the thermal state using the Euclidean path integral. The real-time part of the contour can thence be obtained by slicing open the Euclidean black hole solution around t = i ε and t = −i β + i ε and continue the geometry into the Lorentzian section. The evolution of the slice at t = i ε will give a section of the domain of outer communication of the Lorentzian black hole geometry. The slice at t = −i β + i ε will lead to something similar with a reversed temporal direction. These are the two slices of the geometry living at ζ = 0, 1, respectively.
This basic picture was first espoused in [29], but the prescription of [31] has the added advantage of smoothly connecting the two slices across a 'horizon-cap'. The resulting geometry is smooth and is coordinatized by (2.4). Pictorially, the construction is depicted in figure 2. Now that we have identified the grSK geometry we can explain how to study dynamics thereupon. One should think of all the fields as residing on a complex ζ contour and upgrade the classical bulk action to a contour integral over the mock tortoise coordinate. We write: where x µ are the boundary coordinates. We will use this form of the action for computations of influence functionals. Before we get to those however, we describe how to generalize our construction covariantly to spacetimes with a Killing horizon, and also explain some useful properties of the grSK geometries. 4 It is actually convenient in explicit computations to work with a redefined radial coordinate = r d which is the argument of the hypergeometric function.

JHEP07(2020)242
Covariant grSK spacetime: while we described the construction in a coordinate dependent manner, one can give a more covariant presentation of the same. We can describe this for any spacetime with a smooth future horizon. Pick some intrinsic coordinates on the spatial sections of the horizon, call them x. One can let the temporal evolution be determined by the affine parameter, v, along the horizon generators. A natural radial coordinate r can be chosen by demanding that it be generated by the null normal to the horizon, normalized with respect to the horizon generators, i.e., ∂ v · ∂ r = 1. For a nondegenerate horizon one has ∂ v being timelike outside the codimension-1 null hypersurface (the future horizon). In a local neighbourhood of the horizon one would then end up with a metric of the form (2.1) with f (r) having a simple zero. The details of the rest of the geometry will depend on the asymptotics etc., but the part of the construction that matters for us is indeed the neighbourhood of the future horizon. One can now convert this classical geometry to a two-sheeted geometry with a gluing condition across the horizon as described above by replacing r → ζ.
To illustrate the construction more generally, it is sufficient to consider the nearhorizon region of the spacetime. For non-degenerate black holes this is given by the Rindler geometry. In this case we have the familiar form of the geometry as well as the ingoing coordinatization to be given as: (2.7) Here we have defined ρ ≡ 1 2 r 2 . It is then easy to explicitly identify the mock tortoise coordinate, ζ = 1 2πi log ρ on the primary branch (we work with the Rindler temperature normalized to be 2π). The grSK Rindler geometry would then take the form ds 2 = 2e 2πiζ dv(2πidζ − dv).
(2.8) grSK time reversal: before proceeding further, it is useful to note one useful feature of the grSK geometries. These geometries are not time reversal invariant as is indeed appropriate for the Schwinger-Keldysh dual. However, the Schwinger-Keldysh construction has a Z 2 involution that can be thought of as time-reversal (cf., [17] for a discussion). As described in [33] the geometry does indeed have an involution which can be used to map ingoing solutions to outgoing ones. In the coordinates used in (2.4) the transformation takes the form: where ω is the frequency conjugate to v. More generally, on tensor valued fields the map acts via an idempotent (1, 1) tensor: Time reversal on tensors by contracting indices appropriately, which can be inferred from the action on one-forms and vectors, respectively. These are given to be: (2.11)

JHEP07(2020)242
It is useful to write these equations in terms of ingoing and outgoing modes explicitly, which we denote with superscripts '±'. So one has leading to the general tensor transformation: We will exploit this symmetry to construct solutions of the scalar wave equations in an asymptotically AdS grSK geometry to obtain the corresponding boundary-bulk Green's functions.
3 Open scalar field theory and holographic baths As described in section 1, our goal is to construct the open effective field theory of a single scalar degree of freedom coupled to a holographic thermal field theory. We will start by describing the general set-up and then specialize to the case of two dimensional theories. We will first begin our description by focusing on the 'bath/environment' theory with a holographic dual and how its real-time correlators can be computed via AdS/CFT. We will then describe how this computation amounts to deriving the open effective theory for the probe.

General set-up
Let us consider a scalar probe Ψ(x) coupled to a d-dimensional field theory with fields denoted collectively by X. The latter is taken to be in a thermal state and we let O ≡ O[X(x)] be a local gauge invariant operator in this theory. The action for our system is then We wish to integrate out the thermal degrees of freedom characterized by X and derive an effective action for Ψ. The couplings in S eff [Ψ] are determined by standard sum rules in terms of the Schwinger-Keldysh thermal correlators of the environment variables X [4,5] (see also [11]). So in what follows we will focus on computing thermal Schwinger-Keldysh observables for the environment, with the understanding that these feed into the effective action of our open quantum system. The implications for the open effective theory will be described below in section 3.2.
The thermal field theory we consider will be taken to be holographic. For example we can consider S[X] to refer to a strongly coupled, planar gauge theory in d > 3 (like N = 4 SYM) or a large c 2d CFT. Concretely, in the familiar duality between SU(N ) N = 4 Super Yang-Mills (SYM) and string theory on AdS 5 ×S 5 , the map between parameters is where s is the string length scale, and P the five-dimensional Planck scale. We will work in the regime N 1 and g 2 Y M N 1 whence the holographic system can be described by classical gravitational dynamics on AdS 5 ×S 5 .
The gravitational dual description is in terms of a planar Schwarzschild-AdS d+1 black hole, whose metric in ingoing coordinates was presented in (2.1). The scalar operator, O, is characterized by its conformal dimension ∆ and maps, under the AdS/CFT dictionary, to a scalar field Φ propagating on this black hole background with mass m 2 2 AdS = ∆(∆ − d). In particular, the marginal case ∆ = d corresponds to m 2 = 0 in AdS. The correlation functions of the operator O can be computed by studying the dynamics of Φ in the gravitational theory. For the purposes of our discussion we will model the scalar dynamics by a minimally coupled scalar with a contact self-interaction. For the most part the self-interacting scalar action we work with takes the form: While this will be sufficient to illustrate the general features, it should be borne in mind that one can obtain in actual (top-down) holographic models, such effective action for the bulk fields using dimensional reduction from 10 or 11 dimensional supergravity. The above discussion describes the basic set-up for any AdS/CFT computation. However, for the purposes of our real-time computation we would need to upgrade this action to reside on the gravitational Schwinger-Keldysh black hole geometry (2.4). This is readily done, for we simply change coordinates and rewrite the scalar action as a contour integral over the grSK geometry. To wit, We will study this problem in general dimensions, constructing first the boundary-bulk propagators on the grSK geometry. This involves only the quadratic part of the action; we essentially need to invert the kinetic terms on the grSK geometry. The boundary to bulk propagators will be specified by suitable boundary conditions around the horizon-cap in the geometry (2.4) and non-normalizable boundary conditions characterizing sources on the AdS boundary. We will find that there are two types of propagators: retarded (ingoing) and advanced (outgoing). They will be related in a simple manner by the time-reversal involution identified at the end of section 2. These Green's functions can be obtained by solving the linear scalar wave equation on the grSK geometry. Working in momentum space variables allows mode decoupling as usual. For a general field f on the grSK geometry we adopt the following notational contrivance for Fourier transforms: In general d dimensions we will not be able to solve for the propagator in an explicit analytic manner, as the wave equation on the Schwarzschild-AdS d+1 black hole is not JHEP07(2020)242 known to admit closed form solutions. Hence we will resort to a gradient expansion in the frequencies and momenta, aiming for a low-frequency, long-wavelength expansion. In d = 2 however we will be able to obtain closed form expressions for the propagators. Once we have obtained the boundary to bulk propagators we can compute higher point functions using standard Witten diagram technology in the grSK geometry. Consider the computation of the 4-point function of the operator O in the boundary theory on the Schwinger-Keldysh contour. Of the 4! Wightman functions with various time-orderings, 8 are computed by the Schwinger-Keldysh time-ordering. This is clear from the generating functional (1.5); on the Schwinger-Keldysh contour operators can be inserted either in the forward (R) or backward (L) segments, effectively doubling the number of correlators. These correlation functions are simply related via the Keldysh rules to sequences of nested commutators and anti-commutators of the operator O with suitable time-ordering stepfunctions, see [17,42]. Furthermore, the thermal KMS relations group the correlation functions into orbits of 4 elements (coming from cyclic symmetry of the thermal trace) [43].
It is convenient to introduce a couple of different basis of operators and sources which are convenient in the Schwinger-Keldysh formalism for various computations. First, we introduce the average-difference or Keldysh basis Within the Schwinger-Keldysh literature, the average/difference fields are also sometimes referred to as classical/quantum components respectively. The Keldysh basis naturally separates out an averaged mean field vs the corrections due to quantum/statistical fluctuations. Another useful basis which manifests the KMS relations is the retarded-advanced (RA) basis 5J where n ω is the Bose-Einstein statistical factor: As we shall see below the distinction between the P and F combinations naturally shows up on the holographic side as the distinction between the ingoing modes and the outgoing modes.
The structure of the influence functionals we wish to extract can be now succinctly summarized as follows. Having solved for the field Φ in terms of sources J R and J L on the JHEP07(2020)242 asymptotic boundaries of the AdS grSK geometry the influence functionals are obtained from the generating functional, assuming an n-point bulk contact interaction: In case there are lower-point contact interactions in the bulk, then we should also include tree level Witten diagrams where we have bulk-bulk propagators between the vertices of such lower order contact terms. While self-evident we nevertheless provide a brief argument for this prescription in appendix C. We have given the influence functionals both in the average-difference basis, denoted by I as well as in the retarded-advanced basis, denoted I. We will later write down expressions for them directly in terms of scalar Green's functions. However, there are some general statements that one can make regarding their structural properties prior to any explicit computation, which follow directly from the generating function (1.5). In the average-difference basis the fact that Z SK [J, J] = Tr(ρ β ) implies that This is reflection of the Schwinger-Keldysh collapse rule: the difference operators (equivalently the average source) cannot be futuremost.
In the retarded-advanced basis the fact that we have folded in the statistical factors makes not only the Schwinger-Keldysh collapse rule manifest, but also incorporates the KMS condition. We have We shall see the holographic SK geometry naturally incorporates these relations through the smoothness of the solution across the future horizon cap region that glues together the two sheets of grSK geometry.

Deriving an open EFT from holography
In section 3.1 we have described how one could get a generating functional for Schwinger-Keldysh correlators from holography. This generating functional is evaluated in the presence of a source J for an operator O in the theory with holographic dual. This selfsame JHEP07(2020)242 generating functional has another physical interpretation as emphasized by Feynman and Vernon [4] in their seminal work on open quantum systems. The form of the coupling between the probe and the bath systems in (3.1) suggests that we can interpret the source J of the operator O as the probe field Ψ itself. Further assuming that the dynamics of our probe system is slow, we can then integrate out the 'fast' degrees of freedom which make up the environment. This will obtain for us an effective SK action for this probe, which was termed as the influence functional of the probe in [4]. The Schwinger-Keldysh generating functional of the CFT Z SK [J R , J L ] can hence be given an alternate interpretation as the influence functional of the probe Various structural features we described above for the generating functional can then be re-interpreted as necessary features for the probe influence functional. As described in section 1 if the environment is sufficiently forgetful, one would expect that a local influence functional could be written down for the probe. We will indeed see that the holographic influence functionals derived this way do satisfy this property. Given the difficulty of constructing local influence functionals from perturbative methods (see section 1), this is indeed a fortunate circumstance. We will show that a variety of open EFTs can be derived this way by using grSK geometry and holography.
The KMS conditions and Schwinger-Keldysh collapse rule then have a definite counterpart from the influence functional perspective. One can first of all construct a stochastic field theory that is dual to the influence functional of the open system. In this stochastic field theory, the structural features described above get re-interpreted as non-linear generalizations of fluctuation-dissipation theorems (FDTs). Our holographic construction naturally leads to these non-linear FDTs via the physics of the Hawking radiation and its interaction with the ingoing modes. As far as the authors are aware, this is the first instance where these non-linear FDTs are derived within a field theoretic setup by integrating out bath degrees of freedom. 6 The stochastic description for the dynamics of the open system degree of freedom which follows from the influence functionals S IF [Ψ R , Ψ L ] can be obtained in the following manner. We work in the average/difference basis Ψ a and Ψ d defined analogously to (3.6) for the system variables. The idea is to write the dynamics of the average field Ψ a as a Langevin equation after eliminating the difference or fluctuation field Ψ d . We recall that the average field is the 'classical' variable, where the difference field is the 'stochastic/fluctuation' variable. One starts with the following ansatz for the Langevin dynamics (3.12) The parameters {K, D, γ, θ k ,θ k } are coupling constants, to be determined in terms of the influence phase parameters. The variable η here is the thermal/stochastic noise, with JHEP07(2020)242 strength f ; it is drawn from a non-Gaussian probability distribution To relate the Langevin ansatz (3.12) to an effective action arising from the influence functionals (see eq. (6.1) for an explicit action) we follow the Martin-Siggia-Rose (MSR) trick [46], wherein one converts the stochastic Langevin equation into an effective action using a Lagrange multiplier field Ψ d and the functional integral identity: (3.14) Integrating out the noise field η one obtains the SK effective action for Ψ a and Ψ d . To capture the leading order influence phase we only need to shift η → η + i Ψ d and take the limit η → 0. We then arrive at the Schwinger-Keldysh effective action for our probe system Recall that the standard (linear) FDT relates the friction term γ and stochastic noise f via We find that one obtains non-linear FDTs which relate the non-Gaussian couplings θ k and θ k amongst each other in the form.
This is the advertised FDT for which we will give a derivation once we have derived the influence functionals from holography in section 6.

Scalar propagation in grSK geometries
We now turn to the solutions of the wave equations in diverse dimensions. We will first describe the general set-up, identifying various Green's functions of interest. Our goal is to determine the full solution for the scalar field with prescribed sources J R and J L on the two boundary segments of the grSK geometry. It will turn out that a useful way to proceed is to first identify the ingoing propagator G + , which solves the wave equation with infalling boundary conditions, and thence use time reversal to obtain the outgoing propagator G − (or a linear combination, the Hawking propagator, G H , that we introduce below). Once we have the general framework we will exploit the relative simplicity in d = 2 where the BTZ black hole is a quotient of AdS 3 to find an explicit expression for our propagators. For d > 2 the Schwarzschild-AdS d+1 geometries do not admit closed form JHEP07(2020)242 solutions to the scalar wave equation. However, one can solve for the propagators explicitly in a gradient expansion, i.e., order by order perturbatively in frequencies and momenta. We further elaborate on this gradient expansion scheme in section A exposing some useful technical tricks to organize the solution and extract the propagators.
In what follows we impose Dirichlet (standard) boundary conditions for our scalar field asymptotically (for simplicity), so the conformal dimension of the CFT operator and the mass of the field propagating in the grSK geometry are related by ∆ = d 2 + d 2 4 + m 2 2 AdS .

Scalar boundary to bulk propagators in grSK geometry
The classical equation of motion for a minimally coupled scalar field in the grSK geometry (2.4). Written out explicitly in We can now proceed to solve this problem, but it is useful to understand some structural aspects first. We wish to identify the boundary to bulk propagators for the wave equation (4.2), which allows us to evolve a non-normalizable source at the boundary of the spacetime into a field value at a bulk locale. We will start by first identifying the retarded or ingoing bulk to boundary Green's function G + and subsequently use information about the time reversal symmetry to extract the outgoing or advanced Green's function.
Let us first introduce a new pair of radial derivatives [33] which conjugate to each other in the form: Notice that D ± ζ are related to each other by time reversal. The rationale for introducing them is that these derivatives allow us to absorb the odd powers of ω in the wave equation into themselves. Indeed, in terms of these derivations the scalar equation of motion takes the form: Ingoing boundary to bulk propagator: the ingoing Green's function G in ≡ G + is a solution to (4.5) satisfying a regularity condition at the horizon and normalized to unity at the cut-off boundary of the spacetime.

JHEP07(2020)242
The choice of boundary conditions is such that we are looking at infalling modes across the future horizon, which isolates for us the quasinormal modes (general solution being a superposition of these modes in a linear theory). As such the ingoing Green's function will the one that was obtained in [18] who argue for the ingoing boundary conditions to compute the retarded propagator. We will shortly demonstrate how to obtain G + perturbatively in ω and |k|, i.e., in a gradient expansion in general d and also obtain an explicit analytic form in d = 2.
Outgoing boundary to bulk propagator: once one knows ingoing Green's function G + the advanced or outgoing Green function should be obtained by suitably time reversing it. As argued at the end of section 2, while our coodinatization of the geometry is not time reversal invariant, there indeed an involution realized by the diffeomorphism (2.9). Let us see how this acts on the equation of motion (4.5). First, we note that after reversing the frequency dependence we obtain G − (ω, k) ≡ G + (−ω, k). Using the conjugation relation (4.4) we can then infer that the function G − (ω, k)e −β ω ζ satisfies the wave equation provided Note that (4.7) differs from (4.5) only in the signs of the temporal derivatives i.e., through ω → −ω. It therefore follows that the outgoing Green's function can be obtained in Fourier domain as Full solution and boundary conditions: now that we have the formal expressions for the ingoing and outgoing Green's functions, we can take a suitable superposition to write down the general solution for the linear wave equation. The explicit form of the full solution takes the form We can now impose boundary conditions at the conformal boundary r = r c ± i ε of the grSK geometry. We demand: Using the fact that G + and G − are normalized to unity at these boundaries, we find which results in the solution Hence the general solution to the scalar wave equation (4.9) takes the form

JHEP07(2020)242
where we have used the Bose-Einstein identity: 1 + n ω = e βω n ω . (4.14) It is helpful, before proceeding further, to rewrite the result in terms of linear combinations of the L/R Schwinger-Keldysh sources. For instance, in the retarded-advanced basis (3.7) we find More interesting to us is the solution in the average-difference basis:  Explicitly, it is given by:

Propagators in d = 2
As a warm up we start with the BTZ geometry where d = 2. The metric we recall is where we have written the metric in ingoing coordinates both in the standard AdS radial coordinate as well as in the BTZ adapted global coordinate r = r h cosh ρ. We have either by direct integration or by simplifying (2.5) the following expression for the mock tortoise coordinate where accounted for the cut-off surface where we are imposing our boundary conditions.

JHEP07(2020)242
We can solve the massive, minimally coupled, scalar wave equation (4.2) in terms of hypergeometric functions (4.21) The first of these is the solution that satisfies ingoing boundary conditions and is regular at the future horizon (near ρ ∼ 0, we see that the linearly independent solutions are φ(ρ) = c 1 + c 2 ρ i βω of which the constant behaviour is the correct ingoing mode). To keep expressions compact, we have introduced the lightcone like dimensionless combination of frequencies and momenta, including contributions from the dimension which will appear in the solutions below: The ingoing Green's function of interest, normalized to unit at the AdS boundary ζ = 0 or ρ c + i can then be immediately inferred to be: .

(4.23)
We have written the answer in the mock tortoise coordinate which makes clear that the solution is continuous across the horizon-cap in the grSK geometry. The knowledge of the retarded Green's function is sufficient to obtain the full solution for the field Φ on the Schwinger-Keldysh contour using (4.16).

Propagators in d > 2: gradient expansion
In dimensions d > 2 the scalar wave equation (4.2) does not admit a simple closed form solution in the Schwarzschild-AdS d+1 backgrounds. One can however make progress by solving the equations order by order in a low-energy, long-wavelength limit, i.e., we can expand our Green's function in the limit where βω, β|k| 1. We consider where we continue to work with the dimensionless frequency and momenta. We have isolated the leading order term in the gradient expansion with some hindsight to simplify the computations.

JHEP07(2020)242
To determine the solution for the scalar field, we first solve the equation (4.2) with ingoing boundary conditions. This amounts to imposing the following boundary conditions on the Green's function G + : For the series coefficients in the gradient expansion this translates to the requirement Some of the coefficients above are trivial, spatial reflection symmetry sets G + n,2m+1 = 0, so we do not encounter any terms which are odd in momenta.
It transpires that the knowledge of G + n,m suffices to obtain the solution to the scalar wave equation in the average-difference basis in the long-wavelength gradient expansion. The final expression can be compactly summarized as (4.27) An explicit derivation of the above as is outlined in section A, but the basic strategy is easy to describe. We essentially introduce the bulk analog of the retarded-advanced sources and carry out the gradient expansion both for the Green's function and for the statistical factors that enter in the construction (3.7). We then use the time-reflection symmetry (2.9) to determine the outgoing Green's function, and employ (4.16) to assemble the pieces to give the solution to the wave equation. Closed form solutions for the leading order terms in the gradient expansion can be obtained (see section B). We present below of the basic data that we will use in the rest of the discussion.

Influence functionals
We now have all the pieces necessary to compute the influence functionals of interest. We first outline the computation of the quadratic effective action in section 5.1. As is usual JHEP07(2020)242 in AdS/CFT this is obtained as a boundary term in the on-shell action computation. For the higher order influence functionals we need to employ the standard Witten diagram technology on the grSK geometry (justified in section C). Armed with these results we compute the n-point contact influence functions in section 5.2. Along the way we will argue that the non-Gaussian contributions are well defined after a renormalization of the sources, and determine the appropriate counterterm action necessary to obtain the physical influence functionals. This will turn out to be crucial when our system couples to a marginal operator of the environment theory. We will give some explicit results for cubic and quartic self-interactions in section 5.2. In sections D and E we compile various technical details underlying the results we present in this section.

Quadratic effective action
Let us begin with the evaluation of the quadratic part of the influence functional. Since we have solved for the field Φ on the grSK contour, it follows that result should be given by a boundary term. This is indeed the case, for starting with (3.4) with λ = 0 we have upon passing to momentum space In the second line we substituted the equation of motion (4.5) and performed the ζ contour integral in the final line, expressing the result as a pure boundary term on the grSK contour. The general structure of the answer from evaluating the boundary term takes the following form in the average-difference basis: We see that I aa = 0 a consequence of Schwinger-Keldysh unitarity of the microscopic theory. This is because the coefficient of the average source is the ingoing Green's function which is manifestly regular on the grSK contour.
On-shell action in d = 2: given our solution in the BTZ geometry it is straightforward to evaluate the boundary term. Using the normalized wavefunction built from (4.21), the boundary term contribution evaluates to a, b, c, ξ). The subscript at the end instructs us to extract the coefficient of r 4−2∆ c , which is the end result of carrying out a counterterm subtraction using standard holographic renormalization methods. 7 To understand this recall that we have normalized G + (ζ c , ω, k) = 1 which means that the two point function is obtained from the term that scales like 1 rc 2∆−4 leading to the aforementioned prescription. A short calculation results in:

JHEP07(2020)242
This result is indeed the correct expression for the retarded Green's function G R (ω, k) for a 2d CFT on the infinite line. One can obtain it directly by starting from the conformal 2-point function on the plane, conformally mapping it to the cylinder to obtain the (Euclidean) thermal correlator, and thence take the discontinuity across the lightcone branch cut while analytically continuing it to the timelike Lorentzian domain (using appropriate i prescription to do so). The result has been obtained in various places in the literature before. In [47] the computation of the Fourier transform was described and the imaginary part of the Green's function obtained (in the context of computing the absorption crosssection of D-branes). The first, principled, holographic derivation of the Green's function was given in [18].
The retarded Green's function G R (ω, k) has poles in the lower half complex ω plane, at These poles of course correspond to the BTZ quasinormal modes and set the characteristic scale for the decay of the response function in the time domain [1,48]. Finally, we also note that the expression can be written in a form that is symmetric between the operator O of dimension ∆ and its shadow O s of dimension ∆ = d − ∆ (with d = 2 here). To make the operator dimension's contribution to p ± defined in (4.22) explicit we redefine the light-cone momenta: We can then write the 2-point influence functional in terms of the function 7) 7 We assume, for simplicity, that the field Φ satisfies standard (Dirichlet) boundary condition at infinity.
While this restricts us to ∆ ≥ d 2 , the result is unchanged for ∆ ∈ ( d 2 − 1, d 2 ) once we include additional boundary terms to account for the alternate (Neumann) boundary conditions.

JHEP07(2020)242
as We will find this notation useful in simplifying the analysis of the 3-point influence functional.
Having understood the computation of the influence functional I ad we next can compute I dd . The computation proceeds along similar fashion lines and we obtain Im (I ad (ω, k)) .

(5.9)
We have expressed the final result making manifest the fluctuation/dissipation relation. We recall that Im(I ad ) gives us the spectral function at finite temperature, see [42,43] for further details.
On-shell action in gradient expansion: the analysis of the influence functionals in a gradient expansion is straightforward given our explicit expressions from before. It is however convenient in actuality to assemble the pieces somewhat differently and work in a basis that is better adapted to the bulk field Φ. We describe such a We have been able to evaluate the expressions analytically to leading order in the gradient expansion. We recall that we have normalized G + 0,0 (0) = 1 and use the solution for G + 1,0 given in (4.28), to learn that

JHEP07(2020)242
Expanding out the Legendre polynomial we pick up the a coefficient of the desired power of r c , and find the influence functionals to be While we are retaining terms only to leading order in the gradient expansion one can nevertheless see that the linearized version of the fluctuation dissipation relation (5.9) continues to hold quite generally (as is in fact readily inferred from (5.11)).
The subleading terms in the influence functional I ad are computable. While they do not seem to be amendable to a closed form analytic expression, we can easily extract the dependence on the physical parameters. Using the results for the higher order terms in the gradient expansion given in section B we can show that the quadratic order terms in I ad are given by the following: (5.14) Here g 0,2 and g 2,0 are purely numerical coefficients and computed by the integral expressions involving Legendre functions; see (B.20) and (B.21), respectively. Specifically, they are We have used (B.13) and expressed the answer as integrals over Legendre functions which can be evaluated numerically. For ∆ ∈ d 2 , d 2 + 1 the integrals are absolutely convergent and thus may be determined without need for a detailed counterterm analysis. Representative data for these quantities as a function of ∆ in various dimensions is plotted in figure 4 within this window of conformal dimensions. We will exploit this structural form when writing down the effective action for the open quantum degree of freedom in section 6.

Interactions: contact self-interaction
Using the standard Witten diagrams on the grSK contour, illustrated in figure 5 (we motivate this briefly in section C), the influence functionals can be straightforwardly written down. A contact n-point self-interaction vertex in the bulk leads to the following contribution to the influence functional:  We can in general simplify expressions such as the above by using the fact that the contour integral over the mock tortoise coordinate can be done by basically integrating the discontinuity across the branch cut extending from the horizon over the radial coordinate. To wit, for any function on the grSK geometry L(ζ) where we used √ −g dζ dr = r d−1 . In the second line above, we have assumed that the integrand doesn't have a simple pole at r = r h and hence the horizon itself gives no contribution to the integral. Consequently, entire contribution arises from the discontinuity across the branch cut that extends from the horizon to the conformal boundary.

Influence functionals in the advanced-retarded basis
Using the explicit solution on the grSK contour (4.15) we identify from (3.9) the influence functionals in the retarded-advanced basis to be JHEP07(2020)242 One can evaluate the integrals directly given the solution to the Green's function on the grSK contour. Note that the causal structure of the influence function is completely manifest. The ingoing Green's function G + (ζ, ω, k) is manifestly regular in the upper-half ω plane, while its conjugate G − (ζ, ω, k) is regular on the lower-half plane by frequency reversal. The time-reversal also correctly incorporates the statistical factor which enters the final expression.
In the holographic setting, we expect the only singularities in G + (ζ, ω, k) to arise from quasinormal poles, while those of G − (ζ, ω, k) are attributable to the time-reversed antiquasinormal poles. So the influence functional I F ···F P ···P (k 1 , · · · , k n ) will have quasinormal poles for all the ω i with i ∈ F and anti-quasinormal poles for ω k with k ∈ P . We shall see an explicit illustration of this below.
The radial integral in the expression above, in general needs to be regulated. However, for a certain range of operator dimensions, specifically, for ∆ < (1 − 1 n ) d the integrals are absolutely convergent and need no regulating. 8 In the discussion below we will focus on the case where the operator coupling to our system is sufficiently relevant to avoid introducing any regulators.

JHEP07(2020)242
As an illustration consider the result for the 3-point influence function in a 2d CFT which can be obtained to be (using energy-momentum conservation to eliminate the ω 3 and k 3 ) where the function J δ F F P is given as an infinite sum involving generalized hypergeometric functions: This result is valid for 1 < ∆ < 4 3 , chosen so that the radial integral was absolutely convergent. We have employed a notational contrivance to keep the expression simple -the above is actually a sum over eight terms with similar structure, where each of the three operators enters either as itself, or its shadow. This is indicated by using δ to sample the operator and shadow operator dimension as indicated in the summation. The lightcone momentum of the outgoing mode is conjugated because of time-reversal. We can use energy-momentum conservation to rewrite K * 3± = K 12± as a function of ω 1 + ω 2 and k 1 + k 2 . Finally, we have chosen to write the final expression as a infinite sum, but one can equivalently have presented the result as set of contour integrals (which is where the sum is obtained from). We give a detailed derivation of this result in section D.
Let us focus on some of the features that the readily visible from the above parameterization of the influence functional. Firstly, the result is factorized into right and left movers. This is not manifest in when we consider the AdS radial integral representation, but is expected on grounds that the Euclidean thermal 3-point correlator factorizes into holomorphic and anti-holomorphic parts. Secondly, the influence functional, correctly captures the causality requirements noted above. The only singularities which can arise in the complex frequency plane are from the Gamma functions for ω 1 and ω 2 . One can also check that the generalized hypergeometric functions do not contribute any singularities (this is easier to see directly from the infinite sum or contour integral representation (D.15)).
Of the eight possible choices of δ i we note that there can only be poles when δ i = ∆. When δ i = 2 − ∆ the Gamma function in the denominator also a pole which cancels out the behaviour of the numerator leaving a finite answer. We conclude that the correlator is analytic in the upper half of the complex ω 1 and ω 2 planes and encounters the usual JHEP07(2020)242 quasinormal type poles in the lower half-planes. So the influence functional only has simple poles at are at the following locations determined by the quasinormal modes.
The last set translates upon using energy-momentum conservation to the anti-quasinormal modes of the advanced operator with frequency ω 3 .
We note in closing that one can give a closed form expression for the thermal real-time 3-point functions in 2d CFT in terms of Meijer G-functions. This has been derived by Fourier transforming the cylinder correlator of a 2d CFT to momentum space with an appropriate i prescription in [50]. 9 We have not attempted to simplify our expression (5.19) to this form, though we note that the analytic structure does match between the two expressions (as it must).

Influence functionals in the average-difference basis
If we wish to evaluate the influence functionals in the average-difference basis we can simply implement the basis transformation to go from (3.7) to (3.6). While working in the longwavelength gradient expansion however, one must also account for the statistical factor, which has been folded into (4.27). We will now use this expression to evaluate the influence functionals in the average-difference basis.
As noted above, depending on the nature of the operator O we may need to incorporate counterterms to evaluate the influence functional which enter into the effective field theory of open system degrees of freedom. For the influence functionals to respect the microscopic unitarity, these counterterms must be both state independent and suitably factorized across the SK contour. We see that this is indeed possible, provided that we suitably renormalize the sources for the holographic operator O obtain finite correlation functions. Operationally, from our discussion in section 3.2 this means that we must renormalize our system degrees of freedom while constructing the open effective field theory. 10 We will first describe the structure of the divergences encountered in the computation of the influence functionals themselves. These will appear from the bulk calculation in UV divergent terms from the radial integral, scaling with our radial cut-off r c . We confine our attention to the leading order terms in the gradient expansion, so the results will be accurate to O(βω), for convenience. A more thorough analysis at arbitrary orders in the gradient expansion can be carried out along similar lines, and the details will appear elsewhere.
To begin with let us make the following two assertions about the non-linear influence functional.

JHEP07(2020)242
1. First, the anharmonic influence phase to linear order in the coupling λ n takes the form: Here F b n,k denotes the integral where G + m,n to the desired order are given by the expressions in (4.28). We have anticipated the need for regulating the sources and the integrals and used the superscript b to denote that the above expressions are the bare results, prior to any renormalization prescription.

Second, the integrals F b
n,k appearing in the computation could be divergent. This depends on the conformal dimension ∆ of the operator O. We find: • For ∆ < (1 − 1 n )d the integrals involved in evaluating F b n,k are convergent as noted in section 5.2.1.
• For relevant operators with ∆ ∈ [(1 − 1 n )d, d) we need to regulate some of the integrals. In particular, the integrals F b n,2k+1 with an odd argument are divergent. We can estimate them to behave as follows: where F r n,k denote the renomalized integrals which are completely finite and well-defined as r c → ∞.
• Finally, for a marginal operator ∆ = d all the integrals are divergent. The functions F b n,2k+1 are power-law divergent, but now for even argument F b n,2k one encounters a logarithmic divergence.
The divergences Λ d and Λ l are given by We prove the aforementioned assertions regarding the structure of the bare influence functional and the divergences of the integrals appearing therein in section E. Armed with these statements we can now assert the following result: define the renomalized probes J r a and J r d via

JHEP07(2020)242
While the integrals are divergent for some range of relevant operators, we will argue below that the divergences can be canceled by a standard counterterm. However, we see that there is a need for a non-trivial renormalization of the sources in the case of a marginal operator, arising primarily from the logarithmic divergence encountered in F b n,2k . Note that as r c → ∞ we have lim rc→∞ Λ l Λ d = 0, i.e., the bare and the renormalized sources agree when the cutoff is removed. At finite cutoff however we need to renormalize the sources slightly; this is achieved by mixing them with difference probes in a temperature dependent manner. If the influence functional was completely finite, this deformation would disappear as we take r c → ∞ limit. But given the UV divergences in the case under consideration, this small renormalization, when ignored, can lead to new divergences which look like temperature-dependent divergences which cannot be canceled by unitary, stateindependent counterterms factorized across the SK contour.
However, with respect to the renormalized sources there is no ambiguity and the result is consistent with the requirements delineated above. To see this, define the following counterterm action in terms of the renormalized probes that is unitary, state-independent, and factorized across the SK contour (5.28) With this choice we obtain a manifestly finite influence functional to linear order in couplings λ n and to linear order in derivatives, given by the expression

(5.29)
Once we have performed the renormalization described, we find that we have to evaluate the renormalized integrals F r n,k . We have been able to obtain closed form expressions for a marginal operator with ∆ = d (see section E.1). For the quartic influence functional, passing to momentum space and letting δ(k) = (2π) d δ d 4 i=1 k i , we find: Im(F 4,4 ) Figure 6. The numerical values of the quantities F 3,k and F 4,k as a function of the conformal dimension ∆ (normalized by dimension). We have confined attention to the case of relevant operators ∆ ∈ ( d 2 , n−1 n d) when the integrals are convergent without need of additional counterterms. We note that F n,2k are purely imaginary and F n,2k+1 are real -we have thus plotted just the imaginary or real parts in the corresponding case. sets the thermal scale) and obtain a result up to a numerical factor which is a function of the scaling dimension. These coefficients are given in terms of integrals of Legendre functions. More specifically, we find that the influence functionals are given in terms of the renormalized integrals To derive this expression we used the fact that one can solve for the combination ζ +

Stochastic description of the open effective field theory
Given that we have computed the influence functionals we can return to the problem of constructing the open effective field theory for our system degree of freedom Ψ(x). As remarked in section 3.2 this is easy, since all we need to do at the linear order is replace the sources J a and J d by the average and difference field Ψ a (x) and Ψ d (x), respectively. The renormalized effective action for when the operator Ψ couples to has a n-point contact interaction in the holographic set-up is then immediate to write down, using the results

JHEP07(2020)242
of section 5.1 and section 5.2. Collecting the results from (5.11), and (5.29) the effective action in position space reads: The quantities g 0,2 and g 2,0 are finite numbers (they are functions of ∆ and d) obtained from integral expressions in (5.14), as are F n,k (∆) which is given in (5.32). Sample values of these are also plotted in figure 4 and figure 6, respectively. We can compare the effective action we have derived with the general expectation based on stochastic dynamics of the open system (3.12). The coupling constants of appearing in the Langevin effective action obtained earlier in (3.15) can be read off from (6.1). At the Gaussian order we have the relations while the non-Gaussian terms lead to for k ∈ {1, · · · , n} , for k ∈ {1, · · · , n − 1} . (6.4) The first set of these relations (6.3) capture the standard fluctuation dissipation relations expected from a Gaussian noise source. The fact that θ k andθ k also owe their origins to the same set of influence functionals (and thus to the same underlying microscopic dynamics of the bath fields we integrated out) leads to a set of generalized FDTs quoted earlier in (3.17), which we reproduce here for convenience:

Discussion
We have initiated the study of open quantum systems, where the environment/bath is modeled by a strongly correlated thermal medium with short relaxation times. We modeled the latter using holography, and used standard arguments to relate the real-time thermal observables of the bath to the influence functionals of the open effective field theory. For JHEP07(2020)242 simplicity, our focus was on systems comprising of a single scalar degree of freedom, coupled to a gauge invariant operator of the bath theory of scaling dimension ∆. As we have endeavored to explain, it hitherto has been an unanswered question, at least within the remit of perturbative techniques, whether the dynamics of the system is described by a local influence functional, even assuming that the environment is sufficiently scrambling. As one might suspect, holographic systems, dual to black holes are maximally efficient in scrambling, and thus ought to be able to obtain local influence functions. This is indeed borne out by our explicit analysis, whereby the holographic influence functionals are manifestly local, and provide a proof-of-existence of such local open EFTs. We have furthermore shown that this open EFT can be given a stochastic interpretation, satisfying the required non-linear fluctuation-dissipation relations.
From a gravitational standpoint, the discussion also sheds further light on the semiclassical geometries that compute real-time observables in black hole backgrounds, providing further evidence to the proposal of [31]. In this context, we capture in the influence functionals, not only the known dissipative physics contained in response functions, gravitationally encoded in the quasinormal modes, but also at the same time, learn about their interactions with the outgoing Hawking quanta. While earlier works, [33], demonstrated this in the context of a single Brownian particle degree of freedom, the present discussion generalizes this to scalar probes of the black hole environment. We should note that earlier attempts to understand fluctuations of the outgoing Hawking quanta, as for example in [51], where the authors studied an interacting scalar model in an ersatz two-dimensional black hole were plagued by infra-red divergences. These are due to the fact that probe field has a large effect on the background in the solution chosen, as well as, the fact that the horizon boundary conditions are not imposed appropriately. 11 The calculation we have described in the preceding sections suffers no such ambiguities and indeed reproduces results that are consistent with expectation of the dual boundary field theory.
From the viewpoint of holography, the grSK geometries we have employed provide a satisfactory answer to an old conundrum. How does one efficiently compute higher-point correlation functions in black hole backgrounds, while manifestly respecting causality properties of response functions and fluctuations thereof? One question of particular interest was whether one was required to integrate the interaction vertex of the Witten diagrams throughout the maximal Kruskal extension of the black hole geometry, including regions near the singularity. The grSK geometries give a satisfactory answer: observables are computed in smooth two-sheeted geometries, each extending only up to the future horizon and joined together across a smooth horizon-cap, without any information about the black hole interior, per se. To be clear however, our analysis probes a stationary thermal medium with external sources, which are not themselves backreacting on the medium. We have utilized this to allow the real-time part of the Schwinger-Keldysh contour to extend all the way to t → ∞, despite the fact that evolution future of the latest operator insertion in Lorentzian time is redundant (by the collapse rules). We are as yet unaware of a geometric construction that makes this redundancy manifest, an issue that will be important when we discuss certain generalizations below. 12

JHEP07(2020)242
There are several straightforward generalizations of our analysis such as discussing open systems with spin degrees of freedom coupled to holographic matter [52]. More interesting are situation where the dynamics involves coupling the open system to conserved currents of the holographic theory. For instance, the open system could be coupled to a conserved R-current of the holographic field theory, or to the energy-momentum tensor. In either case, one has to examine the dynamics of gauge fields in the grSK geometries and compute the Schwinger-Keldysh correlators of these conserved currents. As mentioned in section 1, there already exists an analysis in [31] for the dynamics of a Maxwell field in the grSK Schwarzschild-AdS geometry. However, there are several peculiarities of the gauge dynamics which deserve clarification (as already noted in the reference cited). Understanding probe gauge field dynamics is but a precursor to the more interesting question: how does one account for gravitational backreaction and consider holographic environments with non-trivial temporal evolution. A non-trivial example which would be worth analyzing is the couping of our system to a holographic plasma which is spatially inhomogeneous and temporally evolving towards equilibrium. This would require us to understand the gravitational Schwinger-Keldysh construction for the fluid/gravity spacetimes [53].
While we have focused on understanding the influence functionals that are computable by the Schwinger-Keldysh temporal ordering, some of these thermal observables are in turn related by the general KMS relations to certain-out-of-time-order (OTO) observables [43]. For instance, all OTO three-point correlation functions can be generated from two Schwinger-Keldysh correlators, while many of the OTO four-point functions are related to ones with Schwinger-Keldysh ordering (though not the oft studied chaos correlator, which lies in a separate KMS orbit). A natural question is to encode these into the influence functionals, and use the system to probe OTO correlation functions of the holographic bath. This requires a suitable generalization of the grSK geometry to an grOTO saddle, which would be interesting to analyze. Note that the imprints of OTO observables on a Brownian particle probe has been analyzed hitherto in [35].
Finally, in the holographic setting, we have argued that one needs to carry out a suitable renormalization of the open system's degrees of freedom, in order to derive a local effective field theory, whilst respecting microscopic unitarity, state-independence, and factorization across the two legs of the Schwinger-Keldysh contour. This was starkly visible at leading order in the gradient expansion for the coupling of a bosonic degree of freedom to a marginal operator of the holographic thermal field theory. While we have analyzed the renormalization effects at leading order, developing a system holographic renormalization procedure, and understanding the requisite counter-terms, along the lines of [28], would be very helpful. For one these would be invaluable in attempting to construct a Wilsonian effective field theory for open systems (analogous to [54,55]). Of particular interest in this context is to understand how strongly correlated thermal environments evade the issue of non-trivial (infra-red) divergences that appear to arise in perturbative open field theory computations [11,12,14,56].

JHEP07(2020)242
Akhil Sivakumar, and Spenta Wadia for helpful discussions. We would also like to thank Yogesh Dandekar for initial collaboration on the project. CJ and RL would like to thank the hospitality of Indian Institute of Science Education and Research (Bhopal) during National Strings Meeting 2019 (NSM2019) where this work was presented. RL would like to thank the organizers of the "Workshop on Holography, Entanglement and Complexity" at Ashoka University and the workshop on "Recent advances in QFT and gravity, Amplitudes and CFT correlators" at Saha Institute of Nuclear Physics. CJ and RL would also like to acknowledge our debt to the people of India for their steady and generous support to research in the basic sciences. MR would like to thank KITP, UCSB for hospitality during the workshop "Gravitational Holography", where the research was supported in part by the National Science Foundation under Grant No. NSF PHY1748958 to the KITP. MR was supported by U.S. Department of Energy grant DE-SC0009999 and by funds from the University of California.

A Gradient expansion on the grSK contour
In the main text we defined the retarded-advanced basis for the sources on the boundary. It is useful to extend this to the bulk and define two related combinations: The solution for Φ can then be written as with this simple change of basis as Note that G + even , G + odd are obtained by separating the ingoing Green function G + into even and odd parts under frequency reversal once we recall that G − (ζ, ω, k) = G + (ζ, −ω, k).
The advantage of the combinations J even and J odd lies in the following fact: given the ingoing Green function G + in a boundary gradient expansion, one can obtain the solution on the full grSK contour by a simple substitution. We simply multiply the even frequency part of the Green's function by J even and the odd frequency part by J odd , i.e., Φ(ζ, ω, k) =

JHEP07(2020)242
where we resort to the truncated notation: There is another useful property obeyed by the even and odd combinations which is worth highlighting: This can be established by a direct differentiation of the definitions above. It follows that J even and J odd are solutions of the time-reversal invariant differential equation This equation involves no spatial boundary derivatives and hence should be thought of as a differential equation on each ingoing Eddington-Finkelstein tube in the grSK contour (the terminology comes from the fluid/gravity correspondence, cf., [53]). We can then define J even as the solution that interpolates between J L at ζ = 0 to J R at ζ = 1. The odd combination J odd is obtained by differentiation D + ζ J even = w J odd as given from (A.6). More explicitly, we have For computations involving the full solution written down in a derivative expansion we can use the following gradient expansions of the even and odd sources: (A.9) Here we have introduced a new bulk variable y defined as to simplify the expressions somewhat. Using the gradient expansion of J even given above, in (A.8) we end up the result (4.27) quoted in the main text. For certain computations it is helpful to have the solution in the gradient expansion directly in position space. One can verify that (4.27) can be simplified to the following form, accurate to linear order in the derivative expansion (A.11)

JHEP07(2020)242
It is also useful to record here a simple expression for the D ζ derivative of the field which enters into the computation of the quadratic influence functional. It is with G + −1,m ≡ 0. Note that, by construction, J odd occurs in the full solution with at least one βω factor multiplying it. It follows that if G + has a derivative expansion, so does the full solution on the grSK contour written in average-difference/Keldysh basis. Note also that when the even part of the ingoing solution is lifted to the holographic SK contour, the order of the derivative expansion is maintained. In contrast, when the odd part of the ingoing solution is lifted to the holographic SK contour, the source J d occurs with one less time derivative.

B Gradient expansion of the Green's functions
For a massive scalar field m 2 = ∆(∆ − d) = 0 on the grSK geometry we now analyze the wave equation (4.1) in a gradient expansion (4.24). It is helpful to introduce a new radial coordinate and re-express (4.1) as in terms of the following operators: We have also defined a rescaled conformal dimension which will be the only parameter entering our expressions. Upon plugging in the expansion (4.24) we immediately find the recursion relation: with the understanding that G + m,n = 0 for either m < 0 or n < 0. We will now solve this order by order sequentially.
We first note that the leading order term in the expansion, the zero-mode or the DC part, G + 0,0 satisfies the homogeneous equation with no sources:

JHEP07(2020)242
We can exploit this fact to simplify our recursion relation, remove the contribution arising from the mass, and rewrite the differential part of the expression as a total derivative, amenable to integration by quadratures. Define the ratio which one can check satisfies the following equation where we have exploited the fact that G + 0,0 is in the kernel of the differential operator. We impose regular boundary condition at the horizon and normalization at the conformal boundary at each order in perturbation theory as described in (4.26). The general solution can then be immediately written down: where the inner integral has its constant of integration chosen to ensure that the pole at = 1 is canceled.
Zeroth order solution: let us examine some of the leading order terms in the gradient expansion explicitly. The differential operator appearing on the l.h.s. of (B.2) is the Legendre differential operator. Consequently, G + 0,0 satisfies the homogeneous Legendre differential equation (B.6) which has the general solution 13 It is easy to check that the Legendre function of second kind Q −ν (2 − 1) diverges at the horizon, = 1. Hence the normalized solution obeying our boundary conditions is simply This is the result quoted in (4.28). 13 We define the associated Legendre function of the second kind as in [57, eq. 14.3.7] which disambiguates the various definition for this function employed in the literature. In particular, we define it in terms of the regularized hypergeometric function as This expression is well defined for x ∈ (1, ∞) which is the domain of interest to us.

JHEP07(2020)242
First order solution: at the next order it is immediate to see that G + 0,1 (or for that matter any G + m,2n+1 ) vanishes. We then have to solve for G + 1,0 which can be ascertained to satisfy the inhomogeneous equation which gives a solution by quadratures where constants of integration have been chosen to ensure regularity at the horizon = 1 and vanishing of the function at the cut-off surface. We have in addition also made use (B.11) to write the result as an integral over Legendre functions. It is useful to massage this to the form: (B.14) which will enter in the computation of the influence functionals in the gradient expansion. This integral can be explicitly evaluated in terms of Legendre functions of the second kind. In terms of x = 2 − 1 we find after a small algebraic manipulation, the following simple relation where we have identified the Wronskian of the Legendre functions, see [57, eq. 14.2.10]. It then follows that where we fixed the constant by matching the asymptotics of the integral and the Legendre functions. Note that the integral is convergent for ν > 1 2 , though in the main text we will often restrict attention to ν ∈ ( 1 2 , 1] . The bound on ν arises from our focus on relevant and marginal operators (a detailed analysis of convergence properties is given in section E).
Second order solution: at the next order in the gradient expansion we first solve for G + 0,2 which satisfies the inhomogeneous equation (B.17) .
(B.19) For the computation of the quadratic influence functional, we do not need the precise form of these functions. Having their derivatives at the cut-off surface suffices to determine the quantitiesġ 0,2 andġ 2,0 entering I ad (ω, k) in (5.11). These can be computed straightforwardly as we have: where we have integrated by parts and used (B.13) to simplify the integral. These are the expressions that are compiled in (5.14) and (5.15).
Explicit solutions for massless fields: the expressions are in fact simplest for d|∆, i.e., ν ∈ Z + , since then the Legendre functions are simple polynomials (note that P −ν (x) = P ν−1 (x) from the Legendre differential equation). For instance, carrying out the exercise JHEP07(2020)242 we find for a massless, minimally coupled field m 2 = 0 or ∆ = d:

(B.22)
Solutions in the BTZ geometry: it would not be a surprise to reader to know that the above expressions can be integrated in d = 2. We can alternately use the knowledge of the hypergeometric series to aid in expanding out (4.23).
to aid the expansion, but note that we have resum the terms to get the perturbative contributions. The resummation can be done in terms of polylogarithms. For example: Using this we can expand G + (ω, k) to the desired order, for G + (ζ, ω, k) = 1 + i 1 π w log 1 + e 2πi(ζ+ζc) 1 + e 2πi ζc − 1 2 π 2 w 2 log 1 + e 2πi(ζ+ζc) 1 + e 2πi ζc 2 + · · · × 1 + 1 4 π 2 q 2 − w 2 Li 2 sec 2 π(ζ + ζ c ) − Li 2 sec 2 πζ c + · · · (B.25) where we have kept terms to O w 2 and O q 2 , respectively. In other words, the analog of (B.22) now reads to quadratic order: In this appendix we give a short argument in favour of using Witten diagrams on the grSK contour. Per se the argument is no different from the one used in the standard AdS/CFT context on a single sheeted geometry, but for completeness we give a brief account. Consider the action (3.4) whose equation of motion is: We can solve this equation perturbatively with λ as a perturbation parameter, viz., Here Φ 0 solves (4.2) while Φ 1 satisfies the following inhomogeneous equation 3) It will be helpful for our analysis to fix the boundary conditions in a manner that is well adapted to this perturbation theory, so that a minimal number of terms contribute in the evaluation of the action. We pick: Assuming that we have the solution to a desired order we can then evaluate the on-shell action. For instance, if we compute the action up to O(λ n ), as appropriate for a contact bulk interaction, we will simply find the contribution from the zeroth order solution, for: In passing from the first to the second line we plugged in the equation of motion. We then see that as expected the first term in the second line is a total derivative, and thus a boundary term. Since our boundary conditions (4.26) set Φ 1 = 0 at the conformal boundary, so its contribution to the boundary term vanishes. The other term in the expression is the Φ 0 equation of motion which vanishes on-shell. The action with the surviving terms is the expression which only cares about Φ 0 and indeed the contribution to the interaction can be obtained by writing out Φ 0 (ζ, x) in terms of the boundary values using boundary-bulk propagators. Given this, we have chosen drop the subscript '0' from Φ 0 in the main text.

JHEP07(2020)242 D Cubic influence functionals in 2d CFTs
The computation of the influence functional in the RA basis in a two-dimensional CFT is quite general, since conformal invariance on the plane fixes the 3-point function in Euclidean space. One can obtain the thermal correlation function by conformally mapping the complex plane onto the cylinder. A further analytic continuation with suitable i prescription gives the result in Lorentz signature. For a certain choice of operator ordering the result is given in momentum space in [50]. We will now derive the general influence functions in momentum space from the bulk BTZ geometry. It is worth emphasizing that while holography provides a simple way to get the answer, the result is not particularly holographic, since it is constrained entirely by the underlying conformal symmetry.
For definiteness we will compute the result for three scalar primary operators O i with conformal dimensions ∆ i with i = 1, 2, 3. We will assume that the OPE coefficient is C 123 which sets the strength of the bulk vertex (the result presented in the main text corresponds to the special case ∆ 1 = ∆ 2 = ∆ 3 = ∆). We take the bulk action to be given by The cubic influence functional I F F P is then given by where the BTZ Green's function in radial coordinates is given by where we have included the normalization factors which set the source term to unit on the boundary of the spacetime. 14 In the second line we have written out the discontinuity 14 This normalization is different from what is used in the main part of the discussion where we have left the cut-off explicit. Here we choose to normalize the propagator by demanding that limr→∞ r d−∆ G + ∆ = 1. If one uses this normalization then the computation of the two point function proceeds as described in section 5.1 with some minor changes. Firstly, we should pick up the constant term in the asymptotic expansion to compute S (2) . Secondly, there is a renormalization of the operator -the 2-point functions computed in this fashion are rescaled by a factor of ∆ 2(∆−1) relative to the result in (5.4). For the computation of the 3-point function, we should supply a factor r ∆ 1 +∆ 2 +∆ 3 −6 h at the end to ensure that we have the correct scaling of the correlator.

JHEP07(2020)242
across the branch cut which includes a factor of e −2βω 3 ζ in terms of the standard BTZ radial coordinate.
The non-normalizable mode in G ± ∆ leads to a r ∆−2 fall-off. This can be seen from the fact that (with = r h r ) This implies that the radial integral is absolutely convergent 15 if ∆ 1 +∆ 2 +∆ 3 < 4. We will restrict to this range of conformal dimensions to avoid introducing UV regulators in the computation for the present (a discussion of the regulators can be found in section E).
On way to proceed is to use the contour integral representation of the hypergeometric function. A variant of the Barnes integral representation [57, 15.6.7] gives .
This is valid for | arg(1−z)| < π and the choice of the contour C is as follows. One separates the poles of the Gamma functions in the integrand into two sets. The first set includes the poles of Γ(s) and Γ(c − a − b + s) i.e., s ∈ P left = {−n , −n + a + b − c | n = 0, 1, 2, · · · }. The second set comprises the poles of Γ(a − s) and Γ(b − s) i.e., s ∈ P right = {a − n , b − n | n = 0, 1, 2, · · · }. The contour C is chosen as a separatrix between P left and P right . Depending on the range of z and the parameters a, b, c one can either take it to be a vertical contour betweenP left and P right , or a contour that encircles the poles of one of the sets. For the BTZ Green's function we can write the ingoing Green's function using (5.6) and (5.7) as We also have shortened the notation for functions of both K + and K − and written the dependence into simply K ± to allow for more compact formulae below. We will do likewise for G(K ± , ∆). The contour we want to pick is one which encircles the poles in the left set P left = {−n , −n + ∆ − 1 | n = 0, 1, 2, · · · } , (D. 7) and is anchored at −∞. As expected the leading divergence comes from the rightmost pole located at s = ∆ − 1 which gives us the asymptotic behaviour in (D.4). Taking then JHEP07(2020)242 s i ≤ ∆ i − 1 for i = 1, 2, 3 we see that the influence functional receives contribution from the following radial integral We have explicitly used the fact that the integral is absolutely convergent when as argued earlier and have employed momentum conservation to eliminate w 2 . The integral as defined appears to have poles when w 3 = i π, but note that the prefactor (1 − e βω 3 ) also vanishes at that point rendering this harmless. Putting all the pieces together we find: (D.10) In writing the above expression we have used the fact that the sign reversal of ω 3 owing to the contribution coming from outgoing Green's function G − (r, ω 3 , k 3 ) can be expressed by complex conjugation of the lightcone momenta, for K * ± (ω, k) = K ∓ (−ω, k). The main point to note is that the choice of parameters we have made is such that the contribution from the hypergeometric function arising from the radial integral is completely regular in s i . This means that we can close the contours C i to the left picking up the contributions from the poles in P i + for i = 1, 2, 3. This implies that we can write the final result as a triple sum: × G (K 1± , δ 1 + 2n 1 ) G (K 2± , δ 2 + 2n 2 ) G K * 3± , δ 3 + 2n 3 . (D.11) We have written this result in a shorthand notation: accounting for the poles from s i = −n i and from s i = −n i + ∆ i − 1 is tantamount to picking contributions involving either the JHEP07(2020)242 operator dimension or that of its shadow for each i. This is encoded in binary choice of our parameter δ i , viz., The result is expressed in terms of the function G introduced during the analysis of the 2-point function, eq. (5.7), for brevity. To be explicit (D.13) Note that the first term in the parenthesis in the second line of (D.11) contains the residues from the poles and includes a ratio arising from our having expressed the result in terms of the functions G defined above.
To understand the analytic structure it is helpful to reinstate ω 2 and use energymomentum conservation to eliminate ω 3 instead. Noting that we can write the three-point influence functional by replacing G K * 3+ , K * 3− , δ 3 + 2n 3 with the simpler expression G (K 12+ , K 12− , δ 3 + 2n 3 ).
(D. 15) The analytic structure of the influence function can be read off without further effort. The terms on the second line of I F F P (k 1 , k 2 ) involving the sum of n i are regular as a function of ω 1 and ω 2 . The only singularities are from the Gamma functions containing p k± = K k± + ∆ 2 , i.e., the pieces that occur already in the 2-point function. Of the eight possible choices of δ i we note that there can only be poles when δ i = ∆ i . When δ i = 2 − ∆ i the Gamma functions in the denominator factor also have a pole which cancels against that of the numerator leaving a finite answer. We conclude that the correlator is analytic in the upper half of the complex ω 1 and ω 2 planes and encounters the usual quasinormal type poles in the lower half-planes.

E Counterterm analysis for influence functionals
We prove the statements made in section 5.2 regarding the renormalization of the nonlinear influence functional. We will demonstrate the divergence of the integrals noted in the main text and then argue that a suitable temperature dependent mixing of the difference source into the average source serves to give a counterterm action that is consistent with microscopic unitarity. The staring point of our analysis is the time-domain solution of the free massless scalar equation on the gravitational SK contour (A.11) which we reproduce here with the bare sources explicitly marked.
where we are using (B.7). We work to linear order in the gradients, leaving a more detailed analysis of the gradient expansion for the future.

JHEP07(2020)242
To compute the influence phase to linear order in perturbation theory, this unperturbed solution is sufficient. Taking the n th power of this solution, we obtain Φ n n! = n k=0 1 (n − k)! G + 0,0 Note that in the above expression we have expanded the solution accurately to linear order in time-derivatives and we have combined all the contributions coming from i 2 G + 1,0 β∂ t into a total derivative using Leibniz rule.
We now compute the integral over the bulk gravitational SK contour. Only the terms with branch cuts can contribute to the radial integral: this implies we can drop the k = 0 term in the sum above since it is analytic. As expected we get no contribution involving only the average sources. We can also drop the total time derivative in the last line since it gives a boundary contribution to the influence functional (this is partially the reason for working in coordinate space). Thus, we have

(E.3)
Defining the integrals we get the following contribution to the influence functional (E.5) To proceed we need to estimate the integrals. Since we have explicit expressions for the massless scalar field we will first describe the computation in that case, before outlining the general story.

E.1 Divergence structure for a marginal operator
For the massless scalar, for which the radial functions appearing in the gradient expansion are given in (B.22) we can simplify the integrals F n,k defined in (E.4) since G + 0,0 = 1. Dropping the subscript n, since it is unnecessary, we focus on the integrals F defined in (5.23), viz., and immediately obtain ζ + G + 1,0 ∼ 1 r d . Furthermore, since G + 1,0 is continuous across the grSK contour the integrals F b k are simply: To understand the divergence structure it is sufficient to use the asymptotic expansion, though we will use the explicit form.
First consider the situation with the argument k being an odd integer. Then,