Non-equilibrium scalar two point functions in AdS/CFT

In the first part of the paper, we discuss different versions of the AdS/CFT dictionary out of equilibrium. We show that the Skenderis - van Rees prescription and the"extrapolate"dictionary are equivalent at the level of"in-in"two point functions of free scalar fields in arbitrary asymptotically AdS spacetimes. In the second part of the paper, we calculate two point correlation functions in dynamical spacetimes using the"extrapolate"dictionary. These calculations are performed for conformally coupled scalar fields in examples of spacetimes undergoing gravitational collapse, the AdS$_2$-Vaidya spacetime and the AdS$_3$-Vaidya spacetime, which allow us to address the problem of thermalization following a quench in the boundary field theory. The computation of the correlators is formulated as an initial value problem in the bulk spacetime. Finally, we compare our results for AdS$_3$-Vaidya to results in the previous literature obtained using the geodesic approximation and we find qualitative agreement.

C Calculational details on the AdS 3 -Vaidya correlator 32 C. 1 Step 1: Feynman propagator across the shell 33 C. 2 Step 2: Evolve the Feynman propagator to the BTZ side 35 1 Introduction The AdS/CFT duality relates the dynamics of strongly coupled non-equilibirium quantum field theory to semiclassical gravitational dynamics [1,2]. By studying the duality in non-equilibrium settings, one can approach the long-standing problem of non-equilibrium quantum field theory. The gravitational description is particularly well suited for addressing questions of non-equilibrium physics since the bulk physics is classical at leading order in the 1/N expansion and the real time dynamics is reduced to solving differential equations, with initial data specified by the initial state in question. Using the AdS/CFT dictionary, the one point functions of boundary theory operators can be extracted using the asymptotics of the classical bulk solutions. Studying the dynamics of the one point function T µν has lead to interesting results in the dual field theory, such as the fast thermalization of classes of out of equilibrium initial conditions and early applicability of a hydrodynamic description in strongly coupled large N , N = 4 super-Yang-Mills theory [3,4].
The problem of thermalization is by now well studied in the condensed matter literature and it is known that the thermalization of one point functions differs in a qualitative way from thermalization of higher point correlation functions [5]. Even though expectation values of observables localized in a region of space might have reached their thermal values, the quantum correlations between different regions might still be non-thermal. Thus, in order to obtain a better understanding of thermalization in the field theories with gravitational duals, one can study how two point functions approach their thermal values.
One question to address before this is what kind of physical observables are sensitive to the higher point correlation functions? The one point functions tell us the expectation values of physical quantities in the quantum state in question. For example one can study the time evolution of the energy momentum tensor T µν or an order parameter of some symmetry breaking O (see e.g. [6][7][8][9][10][11][12]). Since we are dealing with quantum mechanical systems, the actually measured values of the observables fluctuate. These fluctuations are quantified by higher point correlation functions. For example, the temperature one measures in a subregion of a system has fluctuations that can be quantified by the energy-momentum tensor two point function [13]. On the other hand, (Wightman) two point functions tell us about particle creation rates, such as the photon creation rate of the quark-gluon plasma, studied in the holographic context in [14,15]. As an even simpler example, the two point function of a quantum field φ quantifies the response of a detector coupled to the field φ, as is familiar, for example, from the Unruh effect [16]. From the point of view of out of equilibrium physics in condensed matter systems, one can even directly measure real space two (and higher) point correlation functions and how they approach their thermal limits under time evolution in ultracold atomic gases [17,18].
As is well known, connected two point functions are suppressed in the 1/N expansion compared to their disconnected parts. Thus, from the bulk point of view they are quantum mechanical. The procedure for calculating them in the bulk is the following. To leading order in 1/N one has a classical bulk solution for the metric and all other bulk fields. Then, one considers fluctuations around these classical backgrounds that have to be treated quantum mechanically (see e.g. [19] for a review of the semiclassical approximation in the context of quantum field theory). In particular, one must specify a state for these fluctuations. 1 In this work, we will specialize to a set of initial states that can be prepared with a Euclidean path integral over some manifold, that can be glued to the full Lorentzian spacetime.
Next, one should ask how the bulk initial state is related to the boundary theory initial state, and what is the correct dictionary between boundary and bulk quantities. Building on earlier work [20][21][22], Skenderis and van Rees (SvR) have constructed a machinery for obtaining boundary theory correlation functions by constructing a holographic version of the Schwinger-Keldysh real time formalism [23,24]. In their prescription, the initial state of the boundary theory is prepared by a path integral on a Euclidean manifold M ∂ , while the bulk state is prepared by the path integral over the bulk Euclidean manifolds M whose boundary is M ∂ . Correlation functions are then obtained by taking functional derivatives of the on-shell action on the full glued manifold consisting of the Euclidean manifold M and a Lorentzian part for the real time evolution .
Another, a priori independent dictionary between bulk and boundary correlators is provided by the "extrapolate" version of the AdS/CFT dictionary, where one obtains boundary correlation functions as boundary limits of bulk to bulk correlation functions [25]. In Euclidean time, the "extrapolate" dictionary is equivalent to the standard dictionary where one takes functional derivatives of the on-shell action, as shown in [26], for scalar fields.
In the first part of this paper, we will show that the "extrapolate" dictionary, when generalized to real time (in an obvious way), is equivalent to the SvR dictionary, for "in-in" two point correlators of a free scalar field in an arbitrary asymptotically AdS spacetime. We will demonstrate this by directly constructing a solution to the SvR equations of motion in the glued spacetime and show that this solution is related to a boundary limit of the bulk to bulk two point function. In our opinion, the equivalence of the dictionaries can be interpreted as evidence that they provide the correct dictionary for out of equilibrium physics in AdS/CFT.
In the second part of this paper, we will consider explicit examples of out of equilibrium scalar field two point functions in spacetimes undergoing gravitational collapse. We will consider states where the boundary theory is prepared in the ground state at early times, which is dual to the bulk ground state in AdS spacetime. Then, the CFT is suddenly perturbed out of equilibrium by some operator smeared over all space. The case where the operator is a massless scalar field dual to a marginal scalar operator in the CFT was studied in [27], where it was argued that the spacetime can be well approximated by the AdS-Vaidya spacetime. The AdS-Vaidya spacetime corresponds to pressureless matter shock collapsing along an ingoing null trajectory, and forming a black hole. The massless scalar system was further studied using simulations of numerical gravity in [28,29], which gave further evidence for approximating the Vaidya spacetime to model the CFT process. As a further example, the AdS 4 -Vaidya spacetime is known to arise as an exact dual spacetime of a CFT process where one turns on a time dependent spatially constant electric field in one of the field theory directions [30].
We will consider the AdS vacuum state as the initial state for the scalar field in the bulk. This initial state can in principle be obtained by performing a Euclidean path integral over Euclidean AdS. However, this is not essential for us since the ground state wavefunctional of a scalar field in AdS is explicitly known. As we will review later, knowing the wavefunctional of a free field is equivalent to knowing the equal time two point function of the field. To study how the correlation functions evolve in time, we use the fact that the two point function satisfies the Klein-Gordon equation. This allows us to evolve the two point function in time by solving the Klein-Gordon equation. As the two point function depends on two bulk points x 1 and x 2 , we must solve the equations of motion with respect to both of these coordinates.
This paper is structured as follows. In Section 2, we discuss the AdS/CFT dictionary in out of equilibrium situations. In Section 3, we show how the calculation of boundary two point functions can be written as a bulk initial value problem using the "extrapolate" dictionary. In Section 4, we review the Schwinger-Keldysh formalism and the SvR prescription to calculate out of equilibrium correlation functions and show that it is equivalent to the "extrapolate" dictionary. In Section 5, we consider the example of a massless scalar in the AdS 2 -Vaidya spacetime. In Section 6, we study a conformally coupled scalar in AdS 3 -Vaidya spacetime. In Section 7, we compare the results of our computations in AdS 3 -Vaidya to earlier results obtained using the geodesic approximation.
2 On the holographic dictionary out of equilibrium In our practical calculations, we will use the "extrapolate" version of the AdS/CFT dictionary [25,31]. It directly relates correlation functions in the bulk to correlation functions of the dual boundary operators. For example, for a scalar field of mass m, that is dual to a scalar operator O of dimension ∆ = (d − 1)/2 + (d − 1) 2 /4 + m 2 , the dictionary in Euclidean time is where C n is a constant related to the normalization of the CFT operators. The more conventional version of the AdS/CFT dictionary identifies the boundary generating functional of connected correlators with the Euclidean on-shell action of the bulk theory [32,33]. Then, correlation functions are obtained by differentiating the generating functional of connected correlators where againC n is a constant coefficient related to the normalization of the CFT operators and φ 0 is the coefficient of the non-normalizable solution of the bulk equations of motion, i.e. φ 0 (x) = lim z→0 z −∆− φ(x, z), ∆ − = d−1−∆. For Euclidean asymptotically AdS spacetimes, the two version of the dictionary have been shown to be equivalent [26,31].
For the case of real time out of equilibrium correlation functions, where analytic continuation to a Euclidean spacetime is not available, the situation is less clear. In particular, it is not obvious how to generalize the "differentiate" dictionary (2.3) to such situations to obtain all the correlation functions of interest. For retarded two point correlation functions, the out of equilibrium version of (2.3) can be straightforwardly developed [34] building on previous work [35,36]. Nevertheless, the generalization to other correlation functions, such as multi point Wightman correlation functions, is missing. For a class of states whose wavefunctionals can be obtained as Euclidean path integrals a generalization of the "differentiate" dictionary was constructed by Skenderis and van Rees (SvR) [23,24] (see [20][21][22], for earlier work in this direction). SvR constructed a holographic version of the Schwinger-Keldysh generating functional, from which all the desired correlation functions can be obtained.
On the other hand, the "extrapolate" dictionary (2.1) has an obvious generalization to out of equilibrium correlation functions. First, one should choose the CFT state of interest. We will assume this state can be prepared with a path integral over some Euclidean manifold. We perform this Euclidean path integral in the bulk in a semiclassical approximation to prepare the state of the bulk quantum fields. Then, one should choose the CFT correlation function of interest and make the replacement This way, the calculation of the boundary two point function becomes a standard problem of quantum field theory in curved spacetime.
This procedure can be shown to give correct results for correlation functions in the vacuum and in thermal equilibrium, as is shown for example in Appendix A of [37] for the BTZ black hole. Issues such as choosing ingoing versus outgoing boundary conditions and the problem of boundary terms from the horizon [35] simply do not exist in this formalism. Thus, we view the "extrapolate" dictionary as a natural starting point for out of equilibrium AdS/CFT.

Two point functions as initial value problems
Our method for calculating the bulk two point function is to view the time evolution of the correlation function as an initial value problem. In this work, we will only consider free scalar fields in the bulk. When the one point function of the scalar vanishes, this is sufficient for the purpose of calculating two point functions to leading order in the 1/N expansion. 2 To fix our conventions, we use the following action for a free scalar field In the Heisenberg picture, the bulk scalar quantum field satisfies its Heisenberg equation of motion Moreover, the state |ψ of the system does not vary in time and can be fixed by an initial condition. A straightforward exercise shows that the bulk Feynman two point correlator Thus, if we know G F (x 2 , x 1 ) and its first time derivative on some initial time slice, we can use the equations of motion to evolve it forwards in time. For free bulk fields, this implies that the specification of the initial state is identical to specifying the two point correlation function and its first time derivatives with respect to t 1 and t 2 as initial data.
How do we solve the equations of motion (3.4) given some initial data? A standard solution to this problem is to use the method of Green's functions [34,38,39] where D t = √ −gg tµ ∂ µ and the integral is over all space on a constant time slice. In the following, we will use bold letters for spatial vectors. Thus, the integral measure dx 2 = d d−1 x 2 denotes an integral over a constant time slice. Here, we have chosen the times in the order t 1 < t 2 < t 3 . This formula means that if we know the two point function G F (x 2 , x 1 ) we can obtain it at some later time t 3 by using (3.5). The formula is visualized in Figure 1 in a diagrammatic form.
A derivation of (3.5) goes as follows where we used We will choose the integration region in (3.6) to run over all space from time t = t 2 to some final time t f > t 3 , where t 1 < t 2 < t 3 . Then, using the identity and (3.4), we obtain which gives (3.5) after noting that G R (x 3 , x f ) ∝ θ(t 3 − t f ) = 0 while the boundary term from the AdS boundary vanishes as the two point functions approach zero with the rate z ∆ .
A convenient fact is that the retarded correlator can be obtained from the Feynman correlator using the identity This identity simply follows from the operator definition of G F . Some relations like this among real time correlation functions are summarized in Appendix A. Another convenient fact is that for a free field G R does not depend on the state of the system. This can be understood as follows.
Since the operators φ(x) commute at equal time it follows that G R = 0 on an equal time slice. The first time derivatives of G R are proportional to equal time canonical commutators, e.g.
. Thus, also the first time derivatives of G R are independent of the state on an equal time slice. This means that the initial data for G R is independent of the state at an initial time slice and since G R satisfies (3.7), a second order differential equation, it is independent of the state at all times. 3 The same procedure can be used to propagate the point x 1 in G F (x 3 , x 1 ) forwards in time. For technical reasons that we will explain in a moment, we have found it necessary to rather use the joining formula (3.5) three times in the explicit calculations (3.11) This is visualized in Figure 2 in a diagrammatic form for the case where the bulk spacetime is the AdS-Vaidya spacetime. This is the basis of our numerical calculations. Also we will choose to use an ingoing null coordinate v as the time coordinate, since the Vaidya spacetime is particularly simple in terms of this coordinate. It is convenient to choose the integrals dx 1 and dx 2 to be located at the value of v, when the collapse in Vaidya happens, while the integral dx 0 is on some (arbitrary) earlier value of lightcone time v. Now, we can see why using only two joining integrals fails in the case of Vaidya spacetime. We could have instead performed two joining integrals dx 1 and dx 2 on the lightcone time slice coincident with the shock wave, where the initial data would then be given by the Feynman two point function G F (x 2 , x 1 ) at the time of the infalling shock wave. The problem is that the Feynman two point function has divergences when the points x 1 and x 2 are lightlike separated. In the case when v 1 = v 2 , we only encounter divergences with power −1/2 and one divergence with power −1 at a single point where two lightcone divergences merge. In contrast, for v 1 = v 2 there is a line of divergences with power −1. Thus, we use a third joining integral at some earlier time v 0 to circumvent these lightcone divergences with power −1 in the initial data.
4 Equivalence of two out of equilibrium dictionaries

Review of the Schwinger-Keldysh formalism and the Skenderis-van Rees prescription
Since we are interested in calculating expectation values of operators in a given state |ψ , a standard Feynman path integral is not sufficient as it computes transition amplitudes from initial to final states. This is not suitable for the non-equilibrium setting as we are interested in studying the time evolution of observables without presupposing a final state. To find the correct generalization to non-equilibrium situations, we study the following two point function as an example where U (t 2 , t 1 ) is the time evolution operator from t 1 to t 2 . For concreteness, let us suppose that t 2 > t 1 . The operator O is some local operator built out of the fields in the theory ϕ. A standard time slicing argument (see e.g. [40]) shows that the time evolution operators may be written in a path integral form as where |ϕ is a field operator eigenstate. Replacing the time evolution operators in (4.1) with (4.2), we end up with the path integral where we introduced the wavefunctional Ψ[ϕ] = ϕ|ψ . The ϕ + path integral corresponds to time evolution forwards in time until the time of the latest operator, 4 while the path integral over ϕ − corresponds to time evolution backwards in time towards the initial time. Thus, one can imagine that the path integral is over fields living on a closed time contour, which runs from the initial time t i to some time t m later than the last operator insertion and then back to t i . This time contour is visualized in Figure 3.
Correlation functions can also be computed by functional differentiation of the following Schwinger-Keldysh generating functional (4.4) This generating functional calculates correlation functions that are ordered along the time contour in Figure 3. Continuing earlier work [20][21][22], SvR suggested that the bulk dual of the field theory Schwinger-Keldysh generating functional is given by the bulk path integral [23,24] where the boundary sources are, as usual, identified as the asymptotic values of bulk fields  where Φ collectively denotes all the bulk fields. Using the semiclassical approximation in the bulk, the path integral (4.5) becomes where the bulk fields Φ ± satisfy the saddle point condition Furthermore, SvR considered the case where the bulk initial state Ψ is obtained from a path integral over Euclidean manifolds, that can be glued in to the Lorentzian part of the bulk spacetime. The semiclassical approximation reduces this path integral to a single Euclidean manifold. Then, (4.8) gives rise to the equations of motion for all the bulk fields and conditions for how the fields are glued to the Euclidean manifold, at time t i . In the following, we will assume that the metric and possible other fields are consistently glued. The details of the gluing procedure can be found in [24].
One of the bulk fields Φ is the scalar field φ, with the action (3.1). We are interested in calculating correlation functions of the CFT operator dual to φ. Thus, our task is to solve (4.8) for the scalar field. The wavefunctional of the scalar field in the semiclassical approximation is given by the on-shell Euclidean action (see e.g. [41][42][43] for discussions on wavefunctionals in quantum field theory) where we used the fact that the on-shell action can be written as a boundary term at the gluing surface and that it is a quadratic functional of φ i after using the Euclidean equations of motion. 5 The last equality in (4.9) defines K in terms of the Euclidean manifold and the action functional. This wavefunctional is defined on an initial Cauchy slice, i.e. on a spatial slice t = t i , which we assume to exist. Above, N is a normalization factor chosen so that the norm of the wavefunction is one, i.e.
where the integral over field configurations integrates over spatial configurations at fixed time t = t i . The kernel K(x 1 , x 2 ) is related to the initial equal time two point function of the bulk field φ through The path integral is Gaussian and can be easily performed to give The reader unfamiliar with this fact should consult Problem 3, in Section F of [44].
Thus, K is simply the inverse of the initial equal time two point function and is implicitly determined in terms of G by (4.14) In the following, we will not need the explicit form of the kernel K(x 1 , x 2 ), except that we will assume that it is a real function which satisfies (4.14).
For the scalar field, the saddle point condition (4.8) becomes where The variation of (4.15) with respect to φ ± away from the endpoints of the contour leads to the bulk equation of motion The variation of φ + at the initial time t i gives an initial condition Similarly, the variation of φ − at the initial time gives an initial condition for φ − , which differs from the boundary condition for φ + by a single minus sign. Variation at the "turning point" of the contour, t = t m , gives The delta functional in the original path integral gives one more condition at t = t m , which is the continuity of the fields

Proof of equivalence of the two dictionaries
Instead of finding a general solution to the set of equations, that is obtained in the SvR prescription, we start from an ansatz motivated by the "extrapolate" dictionary and show that it solves all the equations (4.18), (4.19), (4.20), (4.21) and (4.22) As in Euclidean time AdS/CFT, the most general solution to the bulk equations of motion, satisfying the boundary conditions can be written in terms of bulk to boundary propagators where x B is a point at the boundary. In the SvR prescription, the boundary theory correlation functions are determined by taking functional derivatives of the generating functional log Z, i.e. of the on-shell action. This way, the boundary two point functions of interest are given by The other correlation functions can be obtained from these by taking complex conjugates. We will first consider the case where J − = 0 and J + = 0, in which case we need to construct the bulk to boundary propagators K ++ and K −+ . The opposite case J − = 0 and J + = 0 can be obtained from the former one, by performing a time reversal transformation on the complex time contour. Then, due to the linearity of the problem, the general case is a superposition of these two.
Claim: K ++ and K −+ are given by boundary limits of the following bulk correlators which are the bulk Feynman and Wightman functions, computed in the state |ψ . Thus, they satisfy the equations of motion with respect to both of their arguments. Also they satisfy the initial conditions Proof: We want to show that K ++ and K −+ , constructed using the "extrapolate" dictionary, (4.29) and (4.30), are equivalent to the bulk to boundary propagators obtained using the SvR prescription. The equations determining the bulk fields in the SvR prescription, namely the equations of motion (4.18) for φ + and φ − with the boundary conditions (4.19), (4.20), (4.21) and (4.22), translate into conditions for the bulk to boundary propagators using (4.24) and (4.25). We show that (4.29) and (4.30) satisfy these conditions and the correct boundary conditions (4.28).
First, consider the initial condition (4.19). We have to show that Since G F (x, x B ) satisfies the Klein-Gordon equation with a delta function source, we can use a "joining formula" where we choose t 1 = t i + and the integral is performed on a constant time slice. We need to separate t i and t 1 a bit to make sense of the distributions that appear in the following manipulations. The final result is still unambiguous. To check (4.34), we need to take a time derivative To proceed, we need the derivatives These can be calculated by using path integrals and the initial wavefunctional since time derivatives of the field are given by the conjugate momentum Π = −D t φ and the conjugate momentum operator in the path integral becomes a functional derivative with respect to the conjugate field. These path integrals are performed in Appendix B and the results are The above results can be also obtained directly from the Klein-Gordon equation, but we find the path integral derivation more convenient. Using these we obtain where all the terms involving have disappeared, showing that the result is independent of the regulator. This is the left hand side of (4.34). The right hand side of (4.34) is then given by where we again used a joining integral to express G F as an integral over an initial time slice at t 3 = t 1 + δ. Using (4.37) and , which follows from the initial condition for G F (4.33) and from (4.14), it is easy to see that (4.40) is identical to (4.39). This proves that (4.29) satisfies the correct initial condition (4.19).
The same calculation can be repeated for G + (x, x B ). According to (4.20), we have to show that According to the results of Appendix B, the early time derivatives are given by 43) and the joining formula is given by where the integral is performed over the constant time slice at t 1 = t i . Plugging the joining formula into (4.20) and using (4.42) and (4.43), it is easy to see that (4.41) is satisfied.
Thus, we have shown that (4.29) and (4.30) both satisfy the correct initial conditions. It is clear that (4.30) satisfies the right equation of motion, while (4.29) satisfies the right equation of motion everywhere except at the point near the boundary x = x B , where there is a delta function. We will study this point closer in a moment.
Both of the bulk to bulk correlators in (4.29) and (4.30) satisfy normalizable boundary conditions at the AdS boundary, except possibly at x = x B , where the equation of motion for the Feynman propagator has a singularity. Thus, we are lead to conclude that (4.29) and (4.30) satisfy the correct boundary conditions (4.28) at the AdS boundary except possibly at x = x B , where the bulk to boundary correlator has a delta function non-normalizable mode. Above, we assumed the wellknown fact that, given a solution with a vanishing non-normalizable part, time evolution according to the unsourced equation of motion will not generate one. Our next task is to show that (4.29) also has a delta function non-normalizable mode at x = x B . This follows directly from the delta function on the right hand side of the equation of motion (4.32) for G F (see e.g. [39] for a similar proof on a bulk string worldsheet). To see this, we consider the equation of motion in the region of small z, where the geometry is well approximated by AdS At small z, the derivative terms in the boundary directions are subleading and can be neglected.
Rewriting the z derivatives in a convenient way and multiplying the equation with z ∆−d in both sides, we obtain Integrating both sides over z, from zero to z > z B , and using the fact that for z < z B , the Feynman correlator satisfies the normalizable boundary conditions, leads to While G F (x, x B ) decays as z ∆ for z z 0 , it can have a part which decays as z ∆− , for z > z B . Denoting this part as and plugging into (4.47), gives Thus, as long as we first take z B → 0, before z → 0, we obtain which is the correct boundary condition for a bulk to boundary propagator (4.28).
Finally, we should show that (4.29) and (4.30) satisfy the right boundary conditions at the "turning point" of the Schwinger-Keldysh contour. Using the definition of the Feynman correlator we learn that when t > t B we have the identity Thus, (4.29) and (4.30) satisfy the continuity conditions (4.21) and (4.22) at t = t m , for any choice of t m > t B .
Finally, by using the SvR dictionary (4.26) and (4.27), we obtain the boundary correlators which proves the equivalence of the two dictionaries.

Two point functions in AdS 2 -Vaidya
In the rest of this paper, we study examples of correlation functions in non-equilibrium settings in AdS/CFT. We consider a simple class of states, where the system is prepared in the vacuum state and is then suddenly perturbed out of equilibrium. The perturbation backreacts on the geometry and generically leads to black hole formation in the bulk. A simple model of this process is provided by the Vaidya spacetime. In the following, our main interest is the AdS 3 -Vaidya spacetime, which is dual to a 1+1 dimensional CFT. In this Section, we first work out the simpler example of the AdS 2 -Vaidya collapse as a warm-up problem. This serves as an illustration of the calculational method we use in the AdS 3 -Vaidya case, with the advantage that the results in the AdS 2 -Vaidya case are analytic. The results of this Section have been obtained earlier in [38] and another closely related calculation can be found is [45].

AdS 2 -Vaidya spacetime
The AdS 2 metric is given by For future convenience, we introduced an ingoing null coordinate v. Thus, v = v 0 = constant is an ingoing null ray that starts from the boundary at time t = v 0 . Gravity in 1+1 dimensions is in certain sense trivial. For example, a well-known fact is that there are no gravitational waves in dimensions smaller than four. However, there are black hole solutions even in 1+1 dimensions (see e.g. [46] for a more thorough discussion). The one we will be considering has the metric where we use an ingoing null coordinate labelled by the same letter v. This solution has a horizon at z = 1 and a Hawking temperature T H = 1/2π. Due to triviality of gravitational dynamics, this spacetime is in fact locally identical to AdS 2 . Indeed, we can perform the following change of variablesz which leads to It is important to note that the black hole coordinates only cover a part of the AdS 2 spacetime [46], which one might call the Rindler wedge of AdS 2 . Thus, the thermal nature of quantum fields in the black hole spacetime can be understood as a manifestation of the Unruh effect [16] in AdS 2 . In particular, correlation functions of quantum fields in the AdS 2 vacuum state correspond to thermal correlation functions in the black hole coordinates, with temperature given by the Hawking temperature T H = 1/2π.
When 1+1 dimensional gravity is coupled to matter, one can find solutions that describe the collapse of matter into a black hole. Perhaps the simplest case is the AdS 2 -Vaidya spacetime This metric is identical to AdS 2 when v < 0 and to the black hole metric when v > 0. The metric functions have a discontinuity at v = 0, where the collapsing matter is located. Thus, the matter is falling along a null geodesic.
Above, we have only given the 1+1 dimensional metrics of interest, while we have not even specified what the gravitational theory is. For explicit theories leading to the above spacetimes, we refer the reader to the recent work [47] and the references therein.

Calculation of correlation functions
Next, we turn to study correlation functions of quantum fields in the AdS 2 -Vaidya spacetime. For simplicity, we will consider the correlation functions of a massless scalar field with the action (3.1), with m = 0 and d = 2. This action is invariant under Weyl transformations g µν → Ω(x)g µν , which simplifies the task of calculating the two point functions.
The Feynman correlator for a massless scalar in AdS 2 in the above coordinates is given by Its form can be easily understood as follows. Since AdS 2 can be transformed to the upper half plane of Minkowski space (imagining that time runs horizontally) by a Weyl transformation, the two point function in AdS 2 can be calculated in Minkowski space with a Dirichlet boundary condition at z = 0. The correlator of a massless scalar in 1+1 dimensions is given by − log(x 2 + i )/4π. To satisfy Dirichlet boundary conditions at z = 0, one must add a mirror contribution to the two point function from the mirror point at (z m , v m ) = (−z, v + 2z). 6 These two together, give the two point function (5.6). For another derivation of this result see e.g. [46].
Given the coordinate transformation (5.3) from the black hole metric to AdS 2 , the Feynman correlator on the black hole side of the shell is As we explained, this is the two point function in a thermal state in the black hole spacetime, with temperature T H = 1/2π. In the following, we will use the barred coordinates in the v > 0 part of the Vaidya spacetime. This is merely to simplify the algebra. The same result can be obtained by working in the original (z, v) coordinates in the whole Vaidya spacetime. At v = 0, the AdS 2 coordinates are matched to the barred coordinates as follows As initial data, we know the two point functions in the v < 0 part of the AdS-Vaidya spacetime. As discussed in Section 3, the two point function from the AdS to the BTZ side of the AdS-Vaidya spacetime can be obtained by using the joining formula (3.5), which in this case reads where the point (z 0 , v 0 ) is in the AdS part, v 0 < 0, and the point (z 1 , v 1 ) is in the black hole part of the AdS-Vaidya spacetime, v 1 > 0. The superscript a in (5.9) refers to the fact that this correlator is between points across the shell. Integrating (5.9) by parts we obtain a more convenient form where the boundary term vanishes, since The retarded correlator can be determined from the Feynman correlator using the following identity and therefore We can evaluate the above integral to obtain In the next step, we use the fact that the correlator G a F also satisfies the Klein-Gordon equation of motion (with a delta function source) with respect to the variable (z 0 , v 0 ). Thus, we can propagate that point past the shell to the black hole spacetime by using the joining formula (3.11), visualized in Figure 2, We will again only obtain δ-function contributions by first observing and therefore Evaluating the z 0 -integral and using (5.14), we find We can therefore conclude that the bulk Feynman correlator becomes instantly thermal when both points of the correlation function are outside of the infalling shock wave, i.e. for v 1 > 0 and v 2 > 0. This result was indeed obtained in [38]. We can also define boundary correlation functions using the "extrapolate" dictionary leading to One interesting thing to note about the result is that it agrees exactly with the results obtained in [34] using the (complex) geodesic approximation and after setting ∆ = 1. While the calculation in [34] was obtained in the context of the AdS 3 -Vaidya spacetime, the same geodesics can be used in AdS 2 -Vaidya leading to the same result.
6 Two point functions in AdS 3 -Vaidya spacetime

AdS 3 -Vaidya spacetime
Gravity in 2+1 dimensions is trivial in a similar sense as 1+1 dimensional gravity. There are no local graviton excitations and the dynamics of spacetime is directly related to the dynamics of the matter in the theory. 2+1 dimensional gravity still has black hole solutions, which makes it an interesting system to study. The black hole we will be studying is the BTZ black hole [48] (or more precisely black string) with the metric The BTZ black hole has an event horizon at z = 1 and a Hawking temperature T H = 1/2π. As in the case of 1+1 gravity, there is a coordinate transformation that maps the BTZ metric to the with A crucial feature of the coordinate transformation (6.2) is that the black hole coordinates again cover only a part of the AdS 3 spacetime. This part is often called, the AdS-Rindler wedge (see e.g. [49]). As in the AdS 2 case, quantum fields in the vacuum state in AdS 3 will be in a thermal state in the BTZ spacetime, due to the Unruh effect.
The BTZ black hole can be formed by collapse of matter. In the following, we will be studying the AdS 3 -Vaidya spacetime that corresponds to the collapse of null shock wave starting from the boundary at time v = 0. The metric of the AdS 3 -Vaidya spacetime is given by It is easy to see that the metric (6.4) is identical to the AdS 3 metric with v = t − z for v < 0. For v > 0, the metric (6.4) is identical to the BTZ metric (6.1) with v = t − 1 2 log((1 − z)/(1 + z)).

Correlation functions in AdS 3 and BTZ backgrounds
Again, we will consider a free scalar field with the action (3.1) with m 2 = −3/4 and d = 3. In the boundary theory φ is dual to a scalar operator O with scaling dimension ∆ = 3/2. We have chosen the value of the mass m 2 = −3/4 since, with this choice, the action (3.1) can be written in a form invariant under Weyl transformations g µν → Ω 2 (x)g µν and φ → Ω −1/2 (x)φ (see e.g. [50]). This symmetry simplifies the form of the bulk correlation functions.
The Feynman correlator for a scalar field with scaling dimension ∆ = 3/2 in AdS 3 is The result can be understood as follows. By a Weyl transformation with Ω = z, one can transform AdS 3 into the upper half-plane of Minkowski spacetime. The transformation takes the equation of motion of φ into the equation of motion of a massless scalar. The two point function of a massless scalar in Minkowski space is 1/(4π √ x 2 + i ). To satisfy Dirichlet boundary conditions at z = 0, one must include a contribution from a mirror point at (z m , v m ) = (−z, v + 2z). Combining these two leads to the correlation function (6.5), where the overall factor of √ z 1 z 2 appears from the Weyl transformation of the operator φ when one transforms back to AdS 3 . Another derivation of (6.5) can be found in [51].
Translational invariance along the x-direction allows us to Fourier transform the two point function, where K 0 is a modified Bessel function.
As we argued earlier, the ground state in AdS 3 gets mapped in the coordinate transformation (6.2) to a thermal state with the temperature given by the Hawking temperature T H = 1/2π. Thus, the two point function in a thermal state in the BTZ spacetime can be obtained by a coordinate transformation of the AdS 3 two point function 7 (6.5) To compute the Fourier transform of the BTZ-correlator, we use the following integral identity As a check, one can easily verify that (6.7) is periodic in imaginary time, with a period 2π.
We can use this to get the Fourier transform of (6.7): Finally, the retarded correlator in the BTZ spacetime is given by where we replaced θ(t 2 − t 1 ) by θ(v 2 − v 1 ) which can be justified by Im G F vanishing for spacelike separated points.

Calculation of the AdS 3 -Vaidya correlator
As we know the ground state correlator in AdS 3 and the retarded correlator (6.12) in the BTZ spacetime, we are ready to use the joining formula (3.11) to obtain the two point function in the full AdS 3 -Vaidya spacetime. For the purpose of numerics, (3.11) is somewhat problematic. First of all, there are six numerical integrations one should perform, three over the z-coordinates of the joining points and three over the corresponding x-coordinates. Thus, we have a six dimensional integral to perform. Another problem is that the two point functions that one should integrate over have branch point singularities at the lightcone. These singularities are of the type 1/ √ s 2 and 1/(s 2 ) 3/2 , where s 2 is the square of a Lorentzian distance. Thus, a direct numerical integration over the lightcone seems difficult. However, we can use x translational invariance to Fourier transform the joining formula (3.11), which in momentum space becomes where the two point functions appearing in the integrand are given by (6.22) and (6.10). This form (6.13) has two advantages as compared to the position space integral (3.11). Firstly, it has only three integrations instead of six. Secondly, the lightcone divergences in the momentum space two point functions are only logarithmic ∝ log s. Integrals over logarithms can be straightforwardly integrated numerically. We have found it useful to use Mathematica's NIntegrate, together with the command Exclusions, which allows us to exclude the singular points from the numerical integrals. Then, the singular points have to be dealt with analytically. We describe the details of the procedure in Appendix C. The numerical accuracy of our methods is studied in Appendix D. Next, we move on to discuss the results of the numerical integration of (6.13).

Thermalization of the Feynman correlator
Using (6.13), we can now compute the Feynman correlator on the BTZ side of the Vaidya geometry. Firstly, we see that the imaginary part of the propagator thermalizes instantly after the quench, as shown on the right hand side of Figure 4 for the bulk to bulk correlator. This is expected as the imaginary part of the Feynman propagator is related to the retarded correlator, which does not depend on the quantum state for a free field. Therefore, the imaginary part is not affected by the excitation due to the collapse and immediately settles for the BTZ thermal state value. The real part of the Feynman correlator approaches its thermal value at different rates depending on the value of the momentum. On the left hand side of Figure 4, one can see that the real part still equals the AdS value right after the collapse, which is in great contrast to the imaginary part. As time evolves the imaginary part stays constant at the thermal value. Moreover, the real part approaches the thermal value as can be seen from Figure 5. There, we plot the real part of the Vaidya correlator for different values of the lightcone time v 1 after the shock wave while keeping δv = v 2 − v 1 constant. We can see that the real part of the Vaidya correlator, being close to the AdS correlator right behind the shock wave, approaches the thermal correlator at later times.
In Figure 6, we plot the difference of the real part of the boundary Feynman correlator from its thermal value. We consider correlators between points very close to the boundary and rescale the correlators according to the "extrapolate" dictionary (2.5) to obtain the boundary field theory correlator even though we cannot take the strict boundary limit. Thus, we work with the finite cutoff version where the overall factor of 2π is chosen so that the vacuum two point function has a canonical normalization 1/|x 1 − x 2 | 2∆ . The different curves correspond to different momenta, while T is the average boundary time We see that the two point function approaches the thermal one at a different rate for different values of the momentum. The slowest approach is found for the smallest value of the momentum k = 0.01. Defining a strict thermalization time for the correlation function seems difficult, as the correlator is found to oscillate around the thermal value. The higher momenta fall of more rapidly at early times and then undershoot the thermal value and oscillate around it.
One can also ask how the approach to the thermal limit depends on the time difference δt = t 2 − t 1 in the boundary correlation function. Again, the imaginary part is thermal as long as both of the points are at positive times. The time evolution of the Feynman correlator is shown in Figure 7 as a function of the average time T , while the different curves correspond to different values of δt. Figure 7 shows that the time evolution of the correlator is fairly insensitive to the value of δt at sufficiently late times, as all the different curves essentially collapse to one curve. In contrast, the early time dependence, when δt is close to T , is sensitive to the value of δt.

Why does the Feynman correlator thermalize?
We have seen in the previous Section 6. 4  a time dependent spacetime, it might at first seem surprising that the correlator approaches the thermal value at all as free field theories generically do not thermalize. 8 Thermalization of free fields in a collapsing black hole spacetime at very late times was explained by Hawking in [52]. In this Section, we will see how thermalization appears using a formalism more similar to the one in the current paper. Closely related discussions have appeared e.g. in [39].
The retarded correlator propagates the initial data on the shell to the black hole part of spacetime. The retarded BTZ correlator between a point (v 1 = 0, z 1 ) on the shell and a point outside the black hole (v 2 > 0, z 2 ) falls off away from the lightcone as (6.16) As we are interested in the correlators at the late time limit v 2 → ∞, this contribution can be neglected. On the other hand, there are logarithmic divergences when (0, z 1 ) and (v 2 , z 2 ) are separated by radial null geodesics, around which the retarded correlator is order one (as a function of v 2 ). Thus, in order to find the region of space at v = 0 where the initial data is most relevant for late times, we should understand the relevant radial null geodesics. The ingoing null geodesics are simply lines of constant v. On the other hand, the outgoing geodesics are given by where z 1 is the z-coordinate on the initial slice v = 0. A congruence of outgoing null geodesics starting at v = 0 is shown in Figure 8, where the initial radial position z 1 approaches the horizon exponentially for the plotted geodesics. Once the geodesics reach the boundary, they are reflected and fall into the horizon. Therefore, (v 2 , z 2 ) is reached by a direct and a reflected geodesic, which started out at the shell at z 1 = z ± (v 2 , z 2 ), where z ± are defined in (C.4). These two points on the shell, which are lightlike separated from (v 2 , z 2 ), approach the horizon exponentially for growing v 2 as Thus, the part of the initial data that is relevant for the late time boundary two point function is located exponentially close to the horizon. As we will soon review, the near horizon correlator splits into ingoing and outgoing parts. The part relevant for the late time dynamics is the outgoing part only, as the ingoing part rapidly falls through the horizon into the black hole. Consequently, if we show that the outgoing part of the near-horizon Feynman correlator in Vaidya spacetime is thermal after the collapse, then we expect the full correlator in Vaidya spacetime to thermalize at late times.
Our next task is to calculate the outgoing part of the near horizon correlator. To do this, we first note that for v > 0 the Klein-Gordon equation reduces to the equation of motion of a massless 1+1 dimensional scalar near the horizon. This can be seen by using the coordinate 19) in terms of which the Klein-Gordon equation satisfied by the Feynman propagator becomes where we use the Schwarzschild time coordinate t.
, the near-horizon propagator can be split into two parts, one of which only depends on z i * − t i = v i and is therefore ingoing and the other one only depends on z i * + t i = v i + 2z i * and is therefore outgoing. This applies to both i = 1, 2 and we can therefore split the near-horizon Feynman two point function according to where the indices i and o indicate ingoing and outgoing parts with respect to the two coordinates. In order to calculate the outgoing part G F,oo right after the shock wave, we need the initial data from the AdS vacuum correlator. The Fourier-transformed AdS correlator is, as we have seen in (6.22), given by where The horizon in the AdS-part of Vaidya spacetime is given by the outermost outgoing null geodesics that reach z = 1 at v = 0. The outgoing geodesics are parametrized by v + 2z = const. Therefore, the horizon is located at v + 2z = 2, i.e. for fixed v 1 < 0 and v 2 < 0, the near-horizon limit is where the two limits are approached at the same rate. In the near-horizon limit a 1 → 0 while a 2 is finite. Therefore, only the first of the Bessel functions contributes to the near-horizon correlator as it diverges logarithmically when a 1 → 0 . Therefore, the dominant term is This is sufficient to obtain the outgoing part of the initial data. Now, by imposing continuity of the correlator at v = 0 we obtain 9 The above equation implies that the outgoing part is fixed up to a constant to be Near the horizon, we can approximate z ≈ 1 − 2e −2z * which gives A straightforward manipulation of hypergeometric functions shows that the outgoing part of the thermal BTZ propagator (6.10) indeed agrees with (6.28). Another way of seeing thermality in (6.28) is to note that the correlator is periodic in imaginary time v → v + 2πi, recalling that in our units β = 2π.
To summarize, we have shown that the outgoing part of the near horizon correlator in the Vaidya spacetime is identical to the outgoing part of the near horizon thermal correlator in the BTZ background. Thus, as the value of the correlator at very late times is only sensitive to the near horizon correlator, we can conclude that the correlator thermalizes at late times.
In Figure 9, we show that indeed our numerical calculation of the full Vaidya correlator supports the above conclusion. To isolate the outgoing part, we consider the derivative ∂ 2 G F /∂z 1 ∂z 2 , where the derivatives get rid of the ingoing and mixed terms in the near horizon correlator (6.21). We compare this quantity to the corresponding thermal one as a function of z. Indeed, what we find is that just outside the shock wave the outgoing Vaidya correlator agrees with the thermal correlator with good accuracy near the horizon, while the difference between the two increases further away from the horizon.

Comparison with geodesic approximation
The bulk to bulk correlator of a free scalar field is given by an inverse of the operator − m 2 . This inverse operator ( − m 2 ) −1 has a well-known path integral representation as a path integral of a relativistic particle (see e.g. [53,54]) When m is large, we can use a saddle point approximation in the path integral. This way, we obtain an approximate expression for the bulk to bulk correlator as where x cl (τ ) minimizes the length functional L = τ2 τ1 dτ |ẋ cl (τ )|, i.e. it is a geodesic. This procedure is called the geodesic approximation.  Figure 9: The quantity is plotted as a function of z with parameters v 1 = 0.01, v 2 = 0.26 and k = 1. Finite differences have been employed to approximate the derivatives. We can see that the relative difference between the derivatives of the Vaidya and the thermal correlator vanishes close to the horizon signalling that the outgoing part of the near-horizon correlator thermalizes quickly behind the shock wave, while it is not thermal closer to the boundary yet.
In Euclidean time, the inverse in (7.1) is unique and the resulting length functional is real and the approximation above can be justified. The case of real time is more subtle, as discussed in [54]. For real time, the inverse of − m 2 is no longer unique as one can add to it any solution of the equation of motion ( − m 2 )G = 0. This reflects the fact that in real time, one must specify the state of the quantum field φ, while in the semiclassical Euclidean setting, the state is fixed by regularity. Even when ignoring this issue, the saddle point approximation of the resulting path integral is in general somewhat problematic. For spacelike geodesics, the action in the exponent is real so a saddle point approximation would seem viable. Recalling that in a Lorentzian spacetime a spacelike geodesic is not a curve of minimal length, it becomes clear that in order to apply a stationary phase approximation, one must rotate the integration contour in the path integral to complex values of the coordinates to reach the stationary phase contour. Thus, one has to assume that the spacetime metric is analytic in order to justify such a rotation. The metric of the Vaidya spacetime is not analytic as a function of v. Thus, it is not clear whether the geodesic approximation can be applied in this case. Despite these uncertainties, several works have used the geodesic approximation to calculate correlation functions in the BTZ-Vaidya spacetime, finding physically reasonable results for the boundary correlators [55,56]. In this Section, we will present a comparison of how our results for the correlation functions compare to results obtained using the geodesic approximation.
To apply the geodesic approximation, one has to find the spacelike geodesics connecting the corresponding boundary points. For the BTZ-Vaidya spacetime, this was done in [55,56], where we refer the reader to for details of the calculations. The relevant results we need for the following are summarized in Appendix E. As we want to compare the thermalization of the geodesic correlators to our results, we need to calculate the Fourier transform where G th (x, t 2 ; 0, t 1 ) is the thermal correlation function, which is the same for the geodesic approximation as for the exact result (6.10). 10 To obtain some intuition about how the quantity δG is expected to behave, we will first cook up a toy model for the geodesic correlator. For |x| < t 1 + t 2 the geodesic correlator is known to be given by the thermal correlator. On the other hand for |x| t 1 + t 2 , it is given by [56] G(x, t 2 ; 0, t 1 ) ≈ 1 |x| 2∆ (cosh(t 1 /2) cosh(t 2 /2)) 2∆ . (7.4) Thus, our toy model for the geodesic correlator is The Fourier transformed correlator is then simply Some basic features of the integral can be understood as follows. As t 1 and t 2 increase, the integration region over x shrinks as the lower integration limit is given by t 1 + t 2 . As the integrand decays as 1/(t 1 + t 2 ) 2∆ , the whole integral is decaying as a power of (t 1 + t 2 ). Also since the integrand oscillates, the whole integral can be expected to oscillate as cos((t 1 + t 2 )k) up to a phase shift. These expectations can be made more precise by evaluating the integral at large t 1 + t 2 , where the result was obtained by integration by parts. This toy model correlator has the property that short distance correlations thermalize before the long distance correlations. In momentum space this statement is a bit less clear than in position space. From (7.7) we see that at very large times the correlations for arbitrary non-vanishing k decay with the same rate in time e −∆(t1+t2) (t 1 + t 2 ) −2∆ . This is quite different from the position space result. In position space the correlator reaches the thermal value exactly at a time t 1 + t 2 = |x|. This pattern shows up in the rate of oscillation of the Fourier transform, as non-analyticity in real space translates into oscillation in momentum space.
Next, we compare the result of the full geodesic approximation given by first numerically evaluating the correlator in x space using (E.7) and then Fourier transforming the corresponding correlator to obtain δG. Again the integration region starts from |x| = t 1 + t 2 so that we avoid the problem of having to integrate across the short distance singularity as long as t 1 −t 2 t 1 +t 2 . We will consider values for t 1 and t 2 here for which this is true. Results from the geodesic approximation, together with our earlier results, are shown in Figure 10. As the figure shows, the geodesic approximation agrees well with our calculation as far as the qualitative features are concerned. Quantitatively, the results generically differ by an overall order one factor.  Figure 6. The dashed curves are obtained with our numerical method using (6.13), while the solid curves are obtained using the geodesic approximation.

Conclusions
In the first part of this work, we studied different versions of the AdS/CFT dictionary for computing out of equilibrium two point correlation functions. In the first version, developed by Skenderis and van Rees in [23,24], one constructs a holographic version of the Schwinger-Keldysh generating functional. This procedure amounts to calculating the on-shell action for solutions of the bulk equations of motion in a spacetime that is obtained by gluing together Euclidean and Lorentzian spacetimes to construct the Schwinger-Keldysh contour. In this formalism, correlation functions are obtained by taking functional derivatives of the on-shell action. The second version of the AdS/CFT dictionary we discussed was a non-equilibrium version of the "extrapolate" dictionary [25]. In this dictionary, boundary correlation functions are obtained from the bulk correlation functions with the operator replacement O(x) = lim z→0 z −∆ φ(x, z). In Section 4, we explicitly showed that the two dictionaries are equivalent, by showing that the bulk to boundary propagators following from the "extrapolate" dictionary satisfy all the equations of motion and boundary conditions following from the SvR prescription. Thus, the two dictionaries give the same "in-in" two point correlation functions in the boundary CFT. We believe that the equivalence might hold beyond two point functions, but we have not yet constructed a proof of this statement. One approach might be to generalize the path integral approach of [26] to include real time and external wavefunctions. We will leave this problem for future work.
In the second part of the paper, we studied examples of two point functions in dynamical spacetimes corresponding to sudden quenches of the dual CFT. In this kind of quench, one starts from the ground state of the CFT and suddenly perturbs it out of equilibrium by turning on time dependent sources for some operators. This procedure can be modelled with the Vaidya spacetime.
In this work, we studied two point correlation functions of scalar fields in the AdS 2 -Vaidya and AdS 3 -Vaidya spacetimes. As far as we know, this provides the first study of correlation functions in a collapsing spacetime that starts from the AdS d vacuum (with d > 2) and that does not use the geodesic approximation (for correlators in collapse starting from an initial black hole see [57,58]). We provided a straightforward numerical procedure for calculating the correlation functions in the Vaidya spacetime. To simplify the calculational steps, we considered conformally coupled scalar fields. We do not believe that there is any in principle obstruction for generalizing our calculations to other values of scalar field masses and to other bulk fields. In particular, it would be interesting to study how the correlation functions depend on the mass of the bulk scalar m and see whether at large m, the results approach the geodesic approximation. This is a non-trivial question, as the geodesic approximation suffers from several problems discussed in [54] and in our Section 7, so that it is not clear whether the approximation is well justified.
For AdS 2 -Vaidya, the two point correlators for a massless scalar were computed in [38]. Here, we used this example as a warm-up problem, to show how our calculational procedure works. In this case, one finds that the two point correlator thermalizes immediately as both of the points are evolved past the collapsing shell.
For AdS 3 -Vaidya, we found that the Feynman two point function is approaching the thermal two point function. This approach was found to depend on the momentum k. In particular, we found that for small k, the two point function approaches the thermal value slower than for larger k. This is consistent with the general pattern found from many different systems, that longer distance correlations thermalize slower. This has been observed in the geodesic approximation, in the AdS 3 -Vaidya spacetime [55,56] as well as in CFT calculations in certain type of quenches [5] and even in cold atom experiments [18].
In more detail, we compared our results for the two point function to the Fourier transformed geodesic two point functions. The basic qualitative features were found to agree, as can be seen in Figure 10.

A Identities between two point functions
The Feynman correlator is given by where we defined the Wightman function and we used The retarded correlator is given by which can be compactly written as Using (A.1), we also get the relation Another convenient relation is On the other hand the Wightman function G + can be obtained from the Feynman correlator as

B Initial state path integrals
We need the following expectation values where Ψ is the wavefunctional given in (4.9). Taking the functional derivatives gives On the other hand Using these two results we also get where we used D t φ = −Π. Another useful expectation value is Using this result, we obtain The order terms in (B.4), (B.6) and (B.7), can be systematically computed as follows. For a general operator A, one can write Then, the commutators can be determined from the Heisenberg equations of motion. We will not perform this calculation here as we do not need the O( ) terms explicitly.

C Calculational details on the AdS 3 -Vaidya correlator
Following an approach from [34] we compute the Feynman propagator on the BTZ-side of the Vaidya spacetime by time evolution of the propagator on the AdS-side. This is done in two steps: 1. Join two propagators on the AdS-and BTZ-side of the Vaidya spacetime to get the Feynman propagator across the shell.
2. Evolve the Feynman propagator across the shell to get the propagator on the BTZ side of the spacetime.
These two steps can be understood as performing the integrals in (6.13) in two steps. The integrals over z 1 and z 2 yield propagators across the shell. The final integration over z 0 gives the Feynman propagator fully on the BTZ side of the geometry.
C. 1 Step 1: Feynman propagator across the shell As shown in [34], the Feynman propagator across the shell in Vaidya spacetime, i.e. for v 0 < 0 and v 1 > 0, can be computed as where we have used integration by parts with the boundary term vanishing. We can restrict the integration region to the interval [0, 1] as the retarded correlator G BTZ R (v 1 , z 1 ; 0, z; k) vanishes for z > 1, i.e. for (v = 0, z) behind the horizon, due to causality. As our initial state is the pure AdS vacuum, we know the Feynman correlator in (C.1). The same is true for the retarded correlator on the BTZ-side, since the retarded correlator for a free field does not depend on the quantum state after the collapse (whereas the Feynman correlator does such that we cannot infer it from pure BTZ). Since we have computed the pure AdS and BTZ correlators in (6.22) and (6.10), we can compute the integral numerically. However, we need to extract δ-function contributions coming from the derivative of the retarded correlator on the BTZ side. The retarded correlator is has logarithmic divergences at b 1 (v 1 , z 1 ; 0, z) = −1 and b 2 (v 1 , z 1 ; 0, z) = −1, where (v 1 , z 1 ) and (v 0 , z 0 ) are lightlike separated (see (6.11) for the definition of b 1 and b 2 ). At these points, the imaginary part behaves like a Heaviside step function. Consequently, the derivative contributes a δ-function to the integrand. Since this δ-function is not taken into account by numerical integration, we compute it analytically: where z + (v 1 , z 1 ) and z − (v 1 , z 1 ) are the solutions of b 1 (v 1 , z 1 ; 0, z) = −1 and b 2 (v 1 , z 1 ; 0, z) = −1 respectively, C. 2 Step 2: Evolve the Feynman propagator to the BTZ side We can finally compute the Feynman propagator on the BTZ side of the Vaidya spacetime given by for some v 1 , v 2 > 0 and v 0 < 0. Using the results from (C.5), (C.8) and (C.9), we can express the above integral with all δ-functions extracted, While the evaluation of (C.14) is straightforward as all the functions are known analytically, we want to comment on the numerical method we use to compute (C.13). The second and third terms are single integrals specified in (C.6). We evaluate them in Mathematica using the NIntegrate command and Exclusions for the points where the logarithmic divergences occur. The first term is more involved since we have to integrate the functions G a F and G a R which are themselves not known analytically. We evaluate the integrals contributing to G a F and G a R numerically for N different values of z 0 with equal spacing and upper bound z + (v 1 , z 1 )−v 0 /2 as higher values do not contribute to the integral due to causality. Next, we use the command Interpolation to fit a function to the N points. Then we can use the resulting interpolating functions to integrate over z 0 . This integral is again performed with the command NIntegrate including Exclusions to treat the logarithmic divergences. The numerical error and its dependence on N is studied in Appendix D.

D Estimating the numerical error
In this Appendix, we estimate how accurate our numerical procedure for solving the AdS 3 -Vaidya propagator is. As a first check, we can apply the calculational method to evolve the AdS 3 propagator forwards in time. In this case, we know the exact result analytically, which allows us to test the numerical accuracy of our method. In our numerical procedure, we have two sources of error. The first one appears from the numerical accuracy of Mathematica routine NIntegrate where we use Mathematica's Machineprecision, which is of the order 10 −15 . Another source of error appears as we discretize the z space to perform the final integral over z 0 in (6.13). The number of points we use in the discretization is denoted by N . The first measure of error we use is The above measure is plotted in Figure 11 as a function of the number of steps N for a typical values of the parameters specified in the caption. The error is seen to decay approximately exponentially until it reaches a plateau at order 10 −10 . From Figure 11, we see that N controls the numerical accuracy until at some point we reach a plateau, where increasing N does not help anymore. We believe that this plateau appears when the error from NIntegrate becomes larger than the discretization error. When defining the error measure E 1 , we are dividing by the value of the exact correlation function which is by itself of the order 10 −2 , thus amplifying the error from Machineprecision.
As the two point function in AdS 3 -Vaidya is not known analytically, we do not have an obvious quantity to compare it to. The imaginary part of G F is expected to thermalize immediately after both points of the correlation function are in the region v > 0. Thus, as a second estimate of the numerical error we use where now G numerical F is the correlation function in the Vaidya spacetime. This error estimate is shown in Figure 11 as a function of the number of steps N . We see the same effect here as in the estimate E 1 . The error first decays until it reaches a plateau at order 10 −9 . Again, we believe that the plateau appears for the same reason, namely due to the error in Mathematica's NIntegrate. One should note that here the value of Im(G F,th ) is of the order 10 −6 which explains why the plateau in E 2 is larger than that of E 1 by a few orders of magnitude.
As a third measure of the numerical error, we consider where we displayed the dependence of G F on the number of steps explicitly. This estimate tells us how fast the discretized result is converging to a continuum result. We see here that there is no clear plateau appearing in this measure. Thus, the discretized integral is approaching a continuum limit, which as E 2 implies, is not quite the right continuum result, but differs from it by an order of magnitude given by the error in NIntegrate.

E Summary of geodesic approximation results
Let us begin from the results for the equal time correlators in [55]. The regulated geodesic length is given by with c = √ 1 − s 2 . In our numerical approach, we are forced to work with unequal time correlators and thus, we should see how the above result changes when the boundary times t 1 and t 2 are slightly different. The geodesics for this case were calculated in full generality (at spacelike separation) in [56]. Here, we will perform a simpler calculation by assuming that δt is small, which is sufficient for our purposes. We will denote t 1 = t and t 2 = t + δt. Changing the endpoint of the geodesic by a small amount δt will induce a change in the geodesic (x 0 (z), v 0 (z)) → (x 0 (z) + δx(z), v 0 (z) + δv(z)), (E. 4) where (x 0 (z), v 0 (z)) denotes the equal time geodesic parametrized as a function of z. To first order, the length of the geodesic changes by a boundary term δL = ∂L ∂v δt, (E. 5) where L is defined through L = dz L. Using the results in [55], this can be written as 11 δL = c 2ρ δt.