Correlation functions of the Bjorken flow in the holographic Schwinger-Keldysh approach

One of the outstanding problems in the holographic approach to many-body physics is the explicit computation of correlation functions in nonequilibrium states. We provide a new and simple proof that the horizon cap prescription of Crossley-Glorioso-Liu for implementing the thermal Schwinger-Keldysh contour in the bulk is consistent with the Kubo-Martin-Schwinger periodicity and the ingoing boundary condition for the retarded propagator at any arbitrary frequency and momentum. The generalization to the hydrodynamic Bjorken flow is achieved by a Weyl rescaling in which the dual black hole's event horizon attains a constant surface gravity and area at late time although the directions longitudinal and transverse to the flow expands and contract respectively. The dual state's temperature and entropy density thus become constants (instead of the perfect fluid expansion) although no time-translation symmetry emerges at late time. Undoing the Weyl rescaling, the correlation functions can be computed systematically in a large proper time expansion in inverse powers of the average of the two reparametrized proper time arguments. The horizon cap has to be pinned to the nonequilibrium event horizon so that regularity and consistency conditions are satisfied. Consequently, in the limit of perfect fluid expansion, the Schwinger-Keldysh correlation functions with space-time reparametrized arguments are simply thermal at an appropriate temperature. A generalized bilocal thermal structure holds to all orders. We argue that the Stokes data (which are functions rather than constants) for the hydrodynamic correlation functions can decode the quantum fluctuations behind the horizon cap pinned to the evolving event horizon, and thus the initial data.


Motivation and aims
One of the outstanding issues in the holographic duality, that maps strongly interacting quantum systems to semi-classical gravity in one higher dimension, is to understand the dictionary in real time. The applications of the holographic approach to many-body physics are especially limited without explicit and implementable prescriptions for computing out-of-equilibrium correlation functions. Even in the weak-coupling limit, these correlation functions are fundamental tools for studying decoherence and thermalization, e.g. to understand how the commutator and the anticommutator evolve to satisfy the fluctuation-dissipation relation leading to the emergence of the Kubo-Martin-Schwinger periodicity at an appropriate temperature, and how the occupation numbers of quasi-particles equilibrate or evolve to new fixed points [1][2][3][4]. These correlation functions are actually indispensable for understanding dynamics far away from equilibrium in the strong coupling limit where the system cannot be described by quasi-particles. Although one can compute the one-point functions such as the energy-momentum tensor in real time using the correspondence between time-dependent geometries with regular horizons and states in the dual theory, with the remarkable fluid-gravity correspondence [5][6][7][8] providing a primary example, and numerical relativity [9] providing a powerful tool, an explicit computation of a generating functional for hydrodynamic Schwinger-Keldysh correlation functions, as for instance, has not been achieved yet. 1 The object of interest is the generating functional in the dual field theory, whereρ 0 denotes the initial density matrix, the closed time Schwinger-Keldysh contour composed of the forward and backward arms -= ∞ −∞ + −∞ ∞ , where the source J(t, x) is specified such that it is J 1 (t, x) and J 2 (t, x) on the forward and backward arms of the contour, respectively, and T c denotes contour ordering. Formally, we can rewrite the densitymatrix asρ 0 = [Dφ 1 ][Dφ 2 ] ρ 0 (φ 1 , φ 2 ) |φ 1 φ 2 | (1. 2) in terms of a basis |φ of field configurations, and construct the kernel Whenρ 0 is the thermal density matrix, this computation simplifies drastically. One needs to just add an appendage of length T −1 at the end of the closed real time contour parallel to the negative imaginary axis, and impose periodic boundary conditions on the full contour implementing Kubo-Martin-Schwinger (KMS) periodicity [12]. One can also expect a similar simplification for the Bjorken flow which provides the simplest example for evolution of an expanding system on the forward light cone (see Sec 3.1 for details). The state is assumed to have boost invariance, and also translational and rotational invariance along the transverse plane so that the energy-momentum tensor can be expressed only in terms of the energy density T ττ = (τ) via Ward identities, where τ = √ t 2 − z 2 is the proper time of an observer co-moving with the flow and z being the longitudinal coordinate along which the expansion happens. At late time, (τ) is described by hydrodynamics and thus it reaches a perfect fluid expansion, so that (1.6) The full hydrodynamic series for (τ) in powers of τ − d−2 d−1 (essentially a derivative expansion) is given in terms of µ (which is determined by initial conditions) and the transport coefficients which are determined by the fundamental microscopic theory. In this case, it is natural to ask whether there can be a simpler computation of the Schwinger-Keldysh partition function of W[J 1 , J 2 ] in the hydrodynamic limit, since just like the thermal case, the state can be essentially captured by a single parameter.
More generally, we would expect that general methods for computing W[J 1 , J 2 ] would exist in the hydrodynamic regime where the energy-momentum tensor and conserved currents are described only by the hydrodynamic variables, namely the four velocity u µ (t, x), the energy density (t, x) (or equivalently the temperature T (t, x)), etc., and we would not require the knowledge of the detailed (off-diagonal) matrix elements ρ(φ 1 , φ 2 ) of the state or the kernel K(φ 1 , φ 2 ; J 1 , J 2 ) explicitly. In fact, W[J 1 , J 2 ] is related to the generalization of the thermodynamic free energy to hydrodynamics via Legendre transform, and the latter especially in the context of macroscopic spacetime configurations of conserved currents is also known as the large deviation functional [13] which can be computed in many models studied in classical non-equilibrium statistical mechanics. 2 The primary aim of this work is to show how the explicit computation of W[J 1 , J 2 ] can be achieved by holographic methods in the hydrodynamic limit of the Bjorken flow. We will also present concrete steps for understanding how to go beyond the hydrodynamic limit and recover the initial state.

A brief historical review and summary of results
The first major advance in understanding thermal real-time correlation functions in holography was the Son-Starinets prescription for computing the retarded correlation function, according to which the ingoing boundary condition at the horizon implements the causal linear response in the classical gravity (large N and infinite strong coupling) approximation [15]. Using the Chesler and Yaffe method for causal time evolution in the bulk [9], this approach was suitably generalized to compute the out-of-equilibirum retarded correlation function in holography [10]. The first concrete implementation of the Schwinger-Keldysh contour in holography is due to Son and Herzog utilizing the eternal black hole geometry [16]. The two boundaries of the eternal black hole were shown to provide the forward and backward arms of the Schwinger-Keldysh closed time contour with the backward arm displaced by −iβ/2 (note β = T −1 ) along the imaginary axis.
The most concrete prescription for real time gauge-gravity duality for general initial states is due to Skenderis and van Rees [17,18]. This however requires detailed understanding of the state in terms of semiclassical field configurations of dual gravity, and is best defined for states which can be constructed using Euclidean path integrals. In this prescription, one explicitly constructs the bulk geometry corresponding to the boundary Schwinger-Keldysh contour with specified sources, and extends data on the the field theory contour into the bulk in an appropriate manner. It is however, not easy to apply this approach to realistic computations for generic initial states. Furthermore, as mentioned above, we would expect a simpler approach in the hydrodynamic regime. We will compare this approach with ours in Sec 6.
Our method is based on generalizing the recently proposed horizon cap prescription due to Crossley, Glorioso and Liu (CGL) [19] for the static black brane dual to the thermal state. Here the Schwinger-Keldysh contour is realized by a horizon cap, in which the ingoing Eddington-Finkelstein radial coordinate goes around the horizon in the complex plane in a little circle of radius (not to be confused with the energy density) before going back to the real axis and reaching the second boundary. Thus the two arms of the Schwinger-Keldysh contour at the two boundaries are connected continuously in the bulk through the bulk radial contour. The horizon cap implements the appropriate analytic continuation of bulk fields from one arm of the bulk geometry to the other with the sources J 1 and J 2 are specified independently at the two boundaries.
We provide a novel, elegant and simple proof that the CGL horizon cap prescription reproduces the KMS periodicity which relates the different Schwinger-Keldysh correlation functions, and furthermore also reproduces the Son-Starinets prescription for the retarded correlation function at arbitrary frequency and momenta. These were earlier verified only to the quadratic order in frequency [19]. Our proof relies on a simple and novel matrix factorization of thermal correlation functions. This matrix factorization is crucial for the extension of the CGL method to the out-ofequilibrium hydrodynamic Bjorken flow and also for performing crucial consistency checks of the prescription as mentioned below.
The primary tool for extending the CGL prescription into the hydrodynamic regime is the existence of a Weyl rescaling of the Bjorken flow such that the temperature along with the energy and entropy densities becomes a constant at late time. Thus the late time perfect fluid expansion is mapped to a constant temperature. However, there is no time-reversal symmetry because the longitudinal direction expands and the transverse directions contract with proper time evolution in a way such that the total spatial volume density is a constant. This Weyl transformation can be lifted to a bulk diffeomorphism. At late proper time, the dual black hole's event horizon (and also the apparent horizon which coincides with it at late time) reaches a fixed location. The area and surface gravity of the event horizon (and therefore the entropy density and the the temperature of the dual state) remain constant at late time, but the event horizon shrinks in the directions transverse to the flow while expanding in the longitudinal direction. The latter necessitates the viscous and higher order corrections in hydrodynamics (in this Weyl rescaled version) so that the horizon is smooth (these corrections are rather obvious for the original Bjorken flow). The computation of viscous and higher order corrections is just a special case of the fluid/gravity correspondence [6] and has been worked out up to very high orders for the Bjorken flow [20]. The correlation functions for the Bjorken flow can be obtained simply by undoing the Weyl transformation. The Weyl anomaly could contribute to contact terms, which are diagonal in the Schwinger-Keldysh basis and state independent.
Crucially at late proper time, the Schwinger-Keldysh correlation functions can be mapped to the thermal correlation function after spacetime reparametrizations. This follows from the result that at leading order, the Klein-Gordon equation for the massive scalar field dual to the operator can be mapped to that in the static black brane, thus generalizing the result of Janik and Peschanski [21] for the homogeneous case after adequate spacetime reparametrizations manifested via the Weyl rescaling mentioned above. Furthermore the corrections due to viscous and higher order derivative effects in the background can be incorporated systematically, however this requires a further non-trivial log correction to the time-reparametrization.
Let us present the crucial result for the limit of the perfect fluid expansion in some details. The d-dimensional background metric for the Bjorken flow is written in the Milne coordinates: τ, the rapidity ζ = arctanh(z/t) and the transverse coordinates x ⊥ , in which the Minkowski metric takes the form (1.8) Then the Milne metric above is Weyl equivalent to Note that the spatial volume factor is unity although there is expansion and contraction along longitudinal and transverse directions respectively. The above metric should be set as the boundary metric of the dual AdS d+1 geometry (d > 2) which describes the Weyl rescaled Bjorken flow in the field theory, to implement our horizon cap prescription. Although there is no time-like Killing vector even at late time, the dual black hole's event horizon eventually attains a constant surface gravity and area, implying that the temperature, entropy and energy densities in the dual field theory become a constant. In order to find the scalar Schwinger-Keldysh correlation function at late time, it is useful to define (1.10) and take the limit σ → ∞ with other variables fixed. In this limit, the correlation function is essentially governed by the Weyl-transformed perfect fluid expansion in which the temperature and entropy density become constant although time-translation symmetry is not recovered as mentioned above. (Note due to boost invariance, and the transverse translational and rotational symmetries of the Bjorken flow, the correlation function should be functions of σ 1 , σ 2 , ζ r and x ⊥r only. ) One of our main results is that the Schwinger-Keldysh correlation functions in this limit for holographic conformal field theories reduces to thermal correlation functions. 3 Undoing the Weyl transformation, we obtain the following correlation functions of the perfect fluid expansion in terms of those (denoted as G β ) of a static black brane dual to a thermal state: (1.11) (Note G is 2 × 2 matrix in the Schwinger-Keldysh basis and ∆ O is the scaling dimension of the scalar operator). Thus with the manifest spacetime reparametrizations, the correlation functions are thermal at a fixed temperature given by the constant surface gravity of the (Weyl-transformed) dual d + 1-dimensional black hole horizon at late time when it reaches its final radial location ε 1/d 0 . Explicitly, this temperature (which is equal to the static black brane temperature appearing in G β ) is given by with 0 defined in (1.5). The d + 1-dimensional gravitational constant is related to the rank of the gauge group of the dual theory, as for instance in N = 4 super Yang-Mills theory with S U(N) gauge group in 4 dimensions, G −1 N = 2N 2 /π. Remarkably, we obtain a thermal fluctuation-dissipation relation after the necessary spacetime reparametrizations when the state approaches perfect fluid expansion. 3 The factors (d − 1)/(d − 2), σ/τ 0 and (σ/τ 0 ) −1/(d−2) which appear with σ r ,ζ 1 −ζ 2 and x ⊥1 − x ⊥2 below respecitively can be obtained readily from the Weyl rescaled metric (1.9) metric. These are √ −g σσ , √ gˆζˆζ and √ g ii respectively with g ii denoting the components of the diagonal metric in the transverse spatial directions.
Since τ 0 has been chosen arbitrarily, let's consider a scaling τ 0 → ξτ 0 to see if there is any ambiguity in the above result. Note that G β in (1.11) depends only on T 0 σ r and T 0 ζ 2 r + x ⊥ 2 r as it is a thermal correlation matrix of a CFT. Under this scaling, it is evident from (1.8) that σ → ξ 1 d−1 σ, and therefore σ r → ξ (1.5) and (1.12) give Together these imply that G β is invariant under τ 0 → ξτ 0 since T 0 σ r and T 0 ζ r and T 0 x ⊥r are invariant this scaling. However, the Weyl factor in (1.11) scales as ξ 2∆ O d−1 implying that the dimensionless correlation function σ −2∆ 0 G is invariant under the scaling of τ 0 . This holds to all orders.
We can systematically include viscous and higher order corrections to the correlation functions. The correlation functions for the hydrodynamic Bjorken flow can be systematically obtained in an expansion of the form where G 0 coincides with the thermal correlation function G β given by (1.11). Remarkably, the full series (for the hydrodynamic background) can be written as a bi-local generalization of the thermal form in the language of the matrix factorization mentioned above for the static black brane. The key points are 1. At each order in the late proper time expansion in σ −1 , the behavior of the bulk scalar field at the horizon cap must be same as at the zeroth order with the latter in turn mapped to that in case of thermal equilibrium. This is required for consistency with field-theoretic identities.
2. The above is possible only when the horizon cap is pinned to the non-equilibrium event horizon. 4 Our proposal passes many consistency tests. Firstly, using the general identity we can show that even out of equilibrium the retarded response is obtained only from the ingoing mode which is analytic at the horizon cap in consistency with the general prescription of [10]. Secondly, the ansatz naturally reproduces the homogeneous transients (a non-trivial extension of the result of Janik and Peschanski [21] to higher orders). Also our prescription ensures that the advanced response is always given by the outgoing mode which is regular (it has the same kind of branch cut at all orders as in thermal equilibrium, after Weyl rescaling of the Bjorken flow). We will show that this is needed for consistency in the gravitational setup also. The quadratic on-shell action and thus W[J 1 , J 2 ] originates only from the crossterms between the in-going and out-going modes to all orders. This ensures a bi-local generalization of the thermal structure of the correlation functions in terms of the matrix factorization mentioned above, to all orders.
Finally, our result that the horizon cap for the hydrodynamic Bjorken flow is pinned to the non-equilibrium event horizon captures the causal nature of the Schwinger-Dyson equations for the real time (out-of-equilibrium) correlation functions in the dual field theory. 5 The series (1.13) is not expected to be convergent, and therefore requires a trans-series completion with appropriate Stokes data which should be actually functions of σ r , ζ r and x ⊥r . We discuss their physical role in deciphering the information of the initial state which is lost in hydrodynamization, and also how they can be used to decode the interior of the event horizon.

Plan
The paper is organized as follows. In Section 2, we introduce the Crossley-Glorioiso-Liu (CGL) horizon cap prescription for the thermal Schwinger-Keldysh correlation functions in holography. We prove that the prescription indeed reproduces the KMS periodicity so that they are given just in terms of the retarded correlation function, and that the latter is exactly what we obtain from the Son-Starinets prescription. As mentioned, we use a new matrix factorization of thermal correlation functions. In Section 3, we review the Bjorken flow and its holographic dual. We also introduce the Weyl rescaling of the Bjorken flow along with the dual bulk diffeomorphism such that the final state has a fixed temperature and entropy density. As mentioned, although the event horizon has a constant surface gravity and area at late time, it stretches and expands in the directions longitudinal and transverse to the flow respectively. Additionally, we discuss the proper residual gauge transformation corresponding to radial reparametrization.
In Section 4, we study the probe bulk scalar field in the gravitational background dual to the hydrodynamic Bjorken flow and show how we can preserve the analytic structure of the horizon cap to all orders in the proper time expansion. Crucially, we find that it requires the horizon cap to be pinned to the non-equilibrium event horizon. In Section 5, we use these results to extract the real-time correlation functions of the hydrodynamic Bjorken flow. After presenting the result for the perfect fluid limit in terms of a thermal propagator with spacetime reparametrizations, we show how we systematically obtain the corrections in a proper time expansion. We also discuss many non-trivial consistency checks of our results. In Section 6, we present a discussion on how a trans-series completion of this expansion can lead to seeing the quantum fluctuations behind the non-equilibrium event horizon, and matching with initial data lost during hydrodynamization.
Finally, we conclude in Section 7 with an outlook.

The CGL horizon cap of the thermal black brane
The Crossley-Glorioso-Liu (CGL) horizon cap prescription is a simple proposal for the holographic realization of the Schwinger-Keldysh contour at thermal equilibrium [19]. The thermal nature of the correlation functions obtained from this prescription, including their consistency with the Kubo-Martin-Schwinger (KMS) periodicity, has been explicitly verified up to quadratic order in the small frequency expansion. Furthermore, it has been verified that up to this order, the retarded correlation function is implied by the ingoing boundary condition, as demanded by the Son-Starinets prescription [15]. These were sufficient to obtain a rudimentary effective theory of diffusion and dissipative hydrodynamics [22][23][24][25][26] from holography [19,[27][28][29]. Here, we will present a novel and elegant proof that the CGL horizon cap indeed gives thermal correlation functions satisfying KMS periodicity, and that it also implies the Son-Starinets prescription for the retarded correlation function, at any arbitrary frequency and momenta. The generating functional for the thermal Schwinger-Keldysh correlation functions in a quantum field theory is 6 with ρ β denoting the thermal density matrix, 1 and 2 denoting the forward and backward arms of the contour, and T c denoting (time) contour ordering. The contour ordering implies that (with · ≡ Tr(ρ β ·)) Succinctly we can write the above as with (i, j) = (1, 2). The above holds even out of equilibrium with ρ β in (2.1) replaced by an arbitrary initial state ρ 0 . It can readily be shown that the KMS periodicity (arising from the Tr in (2.1) after extending the contour along the negative imaginary axis by β as shown in Fig. 1) implies that the thermal correlation functions in Fourier space (with i and j standing for the 1 (forward) or 2 (backwards) arms of the contour) defined as assume the form is the retarded propagator, is the advanced propagator, and n(ω) = 1/(e βω − 1) is the Bose-Einstein distribution function. It is easy to see from these definitions that The crucial element of the proof of why the CGL prescription works is a simple and general factorization property of thermal correlation functions in field theory (irrespective of whether the theory is holographic or not). The Schwinger-Keldysh thermal correlation functions (2.5) obtained by differentiating the real-time partition function at a temperature T = β −1 can be factorized as shown below is the third Pauli matrix, and Clearly A → λA, B →λB, a → λa and b →λb gives the same thermal matrix, so the factorization is unique up to the multiplicative complex functions λ(ω, k) andλ(ω, k). The CGL horizon cap glues two copies of the black brane geometry, whose boundaries represent the forward and backward arms of the Schwinger-Keldysh time contour respectively, at the horizon as shown in Fig. 2. For reasons to become clear later, this prescription is easily implemented in the in-going Eddington-Finkelstein (EF) coordinates. The in-going EF radial coordinates of the two geometries, representing the forward (1) and backward (2) arms of the time contour respectively, are displaced along the imaginary axis by ∓ (i.e. r 1 → r 1 − i and r 2 → r 2 + i ). The smooth gluing is achieved by the encircling of the complexified radial coordinate around the horizon r = r h clockwise along a circle of radius O( ) as it is analytically continued from the (first) copy dual to the forward contour to the (second) copy dual to the backward contour. The direction of time in the second copy has to be reversed so that full complexified bulk geometry has a single orientation. Therefore, the analytic continuation of the radial coordinate automatically necessitates the closed Schwinger-Keldysh time contour.
Explicitly, the AdS d+1 static black brane geometry in the ingoing Eddington-Finkelstein coordinates is where r is the bulk radial coordinate, t is the Eddington-Finkelstein time and the horizon is at The on-shell action for bulk fields in this geometry is identified with the generating functional of connected real-time correlation functions of the dual operators at the boundary. A bulk scalar field configuration can be written in the form (2.12) On-shell, Φ(r, ω, k) is a sum of two linearly independent solutions φ in (r, ω, k) and φ out (r, ω, k) which are in-going and out-going at the horizon respectively. Therefore, Φ(r, ω, k) = φ in (r, ω, k)p(ω, k) + φ out (r, ω, k)q(ω, k) (2.13) generally with p(ω, k) and q(ω, k) representing the arbitrary Fourier coefficients of the solutions which are in-going and out-going at the horizon respectively. The latter thus provide a basis of solutions for given ω and k, and can be uniquely defined via the following conditions where is the inverse Hawking temperature of the black brane, and r h = ε −1/d 0 is the radial location of the horizon. Indeed, near the horizon (r ≈ r h ), as should follow from the universal validity of the geometrical optics approximation at the horizon. The CGL horizon cap prescription for the analytic continuation of the radial coordinate from one copy of the bulk spacetime to another then implies that the Fourier coefficients of the on-shell solutions in the two copies are related by with 1 and 2 denoting the copies ending on the forward and backward arms of the time contour respectively at the their boundaries. The on-shell solution in the full geometry can therefore be written in the following matrix form: providing a basis of solutions for the entire complexified spacetime comprising of the two copies smoothly glued at the horizon. The sources J 1 (ω, k) and J 2 (ω, k) specified at the two boundaries (see below) implement the Dirichlet boundary conditions that determine p(ω, k) and q(ω, k) uniquely for real frequencies and momenta, and thus yielding a unique bulk field configuration in the full complexified spacetime. According to the holographic dictionary, the generating functional for the connected correlation functions is identified with the on-shell action for the scalar field dual to the operatorÔ, on the full complexified spacetime, i.e.
W[J 1 , J 2 ] = iS on−shell . (2.20) Assuming minimal coupling to gravity, the on-shell quadratic action for the bulk scalar field Φ dual to a scalar operator takes the form The first piece S in−in is quadratic in the ingoing mode. Since the ingoing mode is analytic at the horizon, the contributions from the two arms cancel each other out (as the solutions are the same on the two arms) while the circle around the horizon does not contribute as well. Therefore, S in−in = 0. Note if we keep the in-going mode alone, then J 1 = J 2 . In the field theory, W[J 1 = J 2 ] = 0 because the partition function e W with the same unitary evolution forward and backward in time, equals unity. Therefore, S in−in = 0 is consistent with field theory. The second piece S in−out , which is the sum of cross-terms between the in and out-going modes, has a branch point at the horizon. Integrating over the two arms amounts to integrating around a branch cut, and results in the two boundary terms on-shell, i.e. S in−out = S bdy1 + S bdy2 . (2.22) The third piece S out−out has possibility of a pole at the horizon, i.e. (r h − r) −1 terms in the Lagrangian density arising from the radial derivative acting on the non-analytic piece (r h − r) iβω 2π , which we denote collectively as S . Essentially S gets contributions from the following two terms: Remarkably, the poles originating from these two terms cancel each other out (note β is given by (2.15)) resulting in S = 0. The remaining terms quadratic in the outgoing mode are analytic, so that S out−out is also the sum of two boundary terms. These two boundary terms cancel each other out as in the in-going case. This is less obvious though. On the gravitational side, the easy way to see S out−out = 0 is to note that the in and out-going solutions are gauge-invariant up to overall multiplicative terms. In the outgoing Eddington-Finkelstein gauge, the outgoing solution is analytic at the horizon, and therefore S out−out = 0 exactly for the same reason that S in−in = 0 in the in-going Eddington-Finkelstein coordinates. Therefore, we should have S out−out = 0 in the in-going Eddington-Finkelstein gauge too. 7 If we keep the outgoing modes only, then J 2 (ω, k) = J 1 (ω, k)e βω . In any field theory [26] W where W T stands for the time-reversed process in which we specify with the same density matrix in the future instead of the initial time. 8 For J 1 (t, x) = J 2 (t, x) = J(ω, k), this amounts to The LHS of the above equation vanishes because once again the forward and backward evolution with the same source are inverses of each other (there is no operator insertion in the past now although the state is specified in the future). Therefore, the RHS of the above equation should vanish too, implying that Therefore, S out−out = 0 is consistent with field theory.
The upshot is that we obtain only two boundary contributions from the cross-term between the in-going and out-going modes, so that where S bdy1 and S bdy2 are the contributions from the two boundaries after taking into account counter-terms necessary for holographic renormalization [30]. This implies that This argument is subtle and actually assumes that the KMS periodicity has been implemented. It mirrors the fieldtheory argument presented below. Note that if we go from one gauge to another, then the backward arm of the Schwinger-Keldysh contour shifts along the negative imaginary axis as noted later. In the out-going Eddington-Finkelstein gauge, the backward arm is actually pulled by −iβ along the negative imaginary axis. However, this does not affect this argument. 8 Succinctly, The matrices R and S are defined as follows. Let the asymptotic (r ≈ 0) expansions of the in-going and out-going modes be 9 (2.30) The · · · above stands for (state-independent) contact terms which we ignore. Denoting Therefore, the identification (2.20) together with (2.2) implies that, 12 From the matrix factorization of thermal correlation functions given by (2.9), we readily find from (2.29), (2.30) and (2.31) that the correlation functions obtained by differentiating the on-shell gravitational action are thermal, i.e. assume the form (2.5) provided 13 Remarkably, the above are exactly the Son-Starinets prescriptions [15] for the retarded and advanced propagators according to which they are obtained from the in-going and out-going boundary conditions at the horizon respectively. Furthermore, since the out-going mode is time reverse of the in-going mode (which is not manifest in the Eddington-Finkelstein gauge but can be evident from transforming to Schwarzchild-like coordinates), 14 we should have We therefore conclude that the CGL horizon cap prescription reproduces the Son-Starinets prescription for the retarded propagator together with KMS periodicity and the thermal structure of the correlation functions at any frequency and momentum. A similar approach was adopted earlier by Son and Herzog by identifying the two sides of the eternal black hole with the forward and backward arms of the Schwinger-Keldysh contour [31]. However, in this case, the backward part of the time contour needs to be shifted by β/2 along the negative imaginary axis. The main advantage of the CGL prescription is that we do not need an eternal black hole geometry for its implementation suggesting that its non-equilibrium generalization would be generically more feasible. Furthermore, it is also not clear if out-of-equilibrium correlation functions can be analytically continued in their time arguments as required by the Son and Herzog implementation of the Schwinger-Keldysh contour. Also it should be possible to define integration over bulk vertices and bulk quantum loops in the CGL prescription as well via the analytic structure of the complexified spacetime with the horizon cap. However, this is outside the scope of the present work, and therefore we do not further discuss about this issue. Finally, we note that the arguments presented here are simpler compared to [19] since we do not employ an expansion about ω = 0 which obscures the analytic continuation at the horizon cap by producing (log ω) n terms.
z=t z=-t Figure 3. The schematic diagram of the Bjorken flow illustrates the evolution of an expanding system on the forward light cone. Here z is the longitudinal direction along which the system expands. The transverse directions have been suppressed. Initial data is specified on a constant τ hyperboloid.
proper time τ and the rapidity ζ which are related to the Minkowski (lab frame) time t and the longitudinal coordinate z (along which the system is expanding) as The transverse coordinates x ⊥ are the same in both Milne and Minkowski coordinate systems. In the Milne coordinates, the Minkowski metric takes the form where ds ⊥ is the line element in the transverse plane. The symmetries of the Bjorken flow imply that the expectation value of any operator depends only on the proper time τ. Thus an ansatz for the expectation value of the energy-momentum tensor in a d-dimensional theory, which is also consistent with the transverse translational and rotational symmetries of the Bjorken flow, can take the form in the Milne coordinates for d > 2. Clearly, , p L and p T denote the energy density, longitudinal and transverse pressures respectively. The local conservation of energy and momentum ∇ µ T µν = 0 implies that and the conformal Ward identity T µ µ = 0 imposes Consequently, the evolution of the energy-momentum tensor is determined by (τ) in a conformal field theory. At large proper time τ, the Bjorken flow admits a hydrodynamic description [33]. To explicitly map the energy-momentum tensor (3.2) to that of a fluid, we need to set the flow velocity as i.e. dτ is co-moving with the flow in the Milne coordinates. In a conformal system, the large proper time expansion of (τ) is given by a single parameter, namely which is determined by the initial conditions -0 is a constant energy density, and τ 0 can be chosen to be the value of τ where we intialize. The large proper time expansion of (τ) takes the form where λ n are (state-independent) constants that are determined by the transport coefficients of the microscopic theory. As for instance, λ 1 is related to the shear viscosity η as which should indeed be a constant in a conformal theory. The leading term of the expansion gives an exact solution of the Euler equations, and thus represents the expansion of a conformal perfect fluid.
In what follows, we will need a Weyl transformation of the Bjorken flow. In a conformal theory, the hydrodynamic equations are Weyl covariant [8]. We are ignoring the Weyl anomaly for the moment, but we will explicitly mention it later. Under a Weyl transformation which transforms the metric and the energy-momentum tensor as the new solutions of the hydrodynamic equations are given by in any conformal theory. Consider the combined operation of the time reparametrization and the Weyl scaling with under which the Milne metric (3.1) transforms to (withζ = ζτ 0 ) 12) and the energy-momentum tensor given by (3

.2), (3.3) and (3.4) transforms to
with T ii denoting the diagonal transverse components and denoting the derivative w.r.t. the argument of the corresponding function. It follows that in the hydrodynamic limit, the Bjorken expansion (3.6) takes the resultant form (3.14) The Weyl scaled metric (3.12) has the property that is a constant, and the spatial volume factor is unity, same as in the Minkowski coordinates. However, the longitudinal volume expands, while the transverse volume contracts with the evolution. Also note that for the Weyl scaled Bjorken flow (3.14), we have Therefore, instead of a perfect fluid expansion, the flow attains a constant temperature, energy and entropy densities at late time although no time-like Killing vector exists in the background metric. The latter feature leads to viscous and higher order corrections. The large (reparametrised) proper time expansion is determined by 0 , the final thermal value of the energy density, while τ 0 appears in the Weyl scaling factor Ω as explicit in Eq. (3.11).

Gravitational setup
The gravitational dual of the Bjorken flow [20,[34][35][36][37] has been extensively studied in the literature with the late time evolution providing a primary example of the fluid/gravity correspondence [5][6][7][8] where large order resummation of the hydrodynamic series [20] has been explicitly carried out revealing the hydrodynamization [38] of a far-from-equilibrium state. When a state hydrodynamizes, the energy-momentum tensor can be described an optimally truncated (divergent and asympotic) hydrodynamic series even when it is far from equilibrium [20,38]. In the context of the Bjorken flow, the evolution of the energy density approaches a hydrodynamic attractor [20,38]. This is a generic property of a many-body relativistic system irrespective of whether its degrees of freedom interact weakly or strongly (see [39] for a recent review).
Here we will review the gravitational dual of the Bjorken flow in the hydrodynamic limit and then describe its Weyl transformation in detail. This Weyl transformation is what has been described in the previous subsection. In the bulk it is implemented by an appropriate diffeomorphism. As a result of this transformation, the state reaches a constant temperature and entropy density at late proper time instead of attaining perfect fluid expansion. The dual black hole also attains a horizon with constant surface gravity and area. However, even at late proper time there is time-like Killing vector -the directions longitudinal and transverse to the flow expand and contract respectively such that the horizon area remains constant at late proper time. Along with the Weyl transformation of the metric and the energy-momentum tensor described in the previous subsection, the holographic dual also produces the Weyl anomaly. The Weyl transformation will be an important tool in implementing the horizon cap prescription out of equilibrium.
Additionally, we will focus on the residual gauge freedom which allows us to fix the nonequilibrium event or apparent horizon at a fixed radial location. We will see that it is crucial to pin the horizon cap at the non-equilibrium event horizon for regularity, and therefore this gauge freedom will play an important role. We will explicitly show that this gauge freedom does not affect the dual metric or the dual energy-momentum tensor (and is thus a proper gauge transformation).
The holographic dual of the Bjorken flow in a d-dimensional conformal theory is a (d + 1)dimensional geometry which satisfies the Einstein's equations with a negative cosmological constant: In what follows, we will set L = 1 for convenience. In addition to the field theory coordinates, we need an extra radial coordinate to describe the dual geometry. The state of the conformal theory dual to a specific solution of (3.17), lives at the boundary (ρ = 0) in the boundary metric, which is defined as where a and b stand for the field theory indices. Since we are considering the Bjorken flow in the Milne metric (3.1), g b µν should coincide with it. Similarly, if we consider the Weyl scaled version of the Bjorken flow, the boundary metric should coincide with (3.12).
Before considering the Bjorken flow, it is useful to first understand the vacuum solution, which is pure (maximally symmetric) AdS d+1 spacetime with the desired boundary metric. In the ingoing Eddington-Finkelstein gauge, the vacuum state in the Milne metric (3.1) is thus dual to where r is the radial coordinate. Similarly, the vacuum in Weyl scaled metric (3.12) is dual to where v is the radial coordinate. These bulk metrics (3.19) and (3.20) are related by the diffeomorphism (3.21) For both cases, (3.19) and (3.20), we obtain the boundary metrics (3.1) and (3.12) from (3.18), after replacing ρ with r and v, respectively. Any Weyl transformation at the boundary is dual to a bulk diffeomorphism. Since the boundary metrics (3.1) and (3.12) are related by a Weyl transformation, (3.21) is simply a specific instance of this general feature of holographic duality. Note that τ and σ are related exactly by the time reparametrization (3.10) at the boundary. 15 Holographic renormalization [30,[40][41][42] provides the framework for extracting the T µν corresponding to the state in the field theory dual to a specific asymptotically AdS d+1 bulk geometry. The procedure essentially amounts to covariantly regularizing the Brown-York tensor on a cut-off hypersurface with local counterterms built out of the induced metric, and then taking this surface to the boundary. For the bulk geometry (3.19), T µν = 0 in the dual vacuum state living in the flat Milne metric (3.1) at the boundary. On the other hand, for the vacuum state living on the Weyl transformed Milne metric (3.12) which is dual to the bulk geometry (3.20), T µν = 0 only if d is odd. For even d, holographic renormalization reproduces the Weyl anomaly of the dual field theory. In the case of d = 4, we obtain (using minimal subtraction scheme) whereg denotes the Weyl rescaled background metric (3.12), R µν is the Ricci tensor built out of it, etc. It is easy to verify that i.e. energy and momentum is conserved in the Weyl rescaled background metric (3.12) (with ∇ being the covariant derivative built out of it), and (3.25) we can readily find that (3.25) reproduces the Weyl anomaly of S U(N) N = 4 super-symmetric Yang-Mills theory [40,41]. An asymptotically AdS d+1 metric dual to a Bjorken flow on the flat Milne metric (3.1) at the boundary, is a solution to the vacuum Einstein's equations (3.17) which takes the form with the following Dirichlet asymptotic boundary conditions These boundary conditions ensure that the boundary metric (3.18) (with r being the radial coordinate in place of the generic ρ) coincides with the Milne metric (3.1). The Einstein equations (3.17) can be readily solved in the late-time expansion, as functions of the scaling variable and with the expansion parameter being 0 is a constant which will be related to the single parameter µ of the Bjorken flow defined in (3.5) (or equivalently to 0 ) below. Explicitly, The functions a (i) , l (i) and k (i) satisfy ordinary differential equations with source terms at each order. We require that these functions do not blow up at the perturbative horizon which is at Together with the Dirichlet boundary conditions (3.28), these finiteness conditions ensure that we obtain solutions which are free of naked singularities in the perturbative expansion and which are unique up to terms which are determined by a single coefficient [20]. This coefficient captures the residual gauge freedom of the ingoing Eddington-Finkelstein coordinates which is the reparametrization of the radial coordinate r (without spoiling the manifest translational and rotational symmetries along the transverse directions). Usually, this residual gauge freedom is fixed by setting the radial location of the apparent or event horizon at (3.32) to all orders in the late proper time expansion [9]. However, in what follows, we will show that the residual gauge freedom should actually be fixed by the regularity of the horizon cap. Note that although the requirement that the horizon cap should be pinned to the evolving event horizon for regularity is a gauge-invariant statement, the residual gauge freedom will be crucial to implement the prescription in the Eddington-Finekelstein coordinates. We therefore keep this gauge freedom unfixed here and later show how it is fixed by regularity of the horizon cap such that the latter is pinned to the evolving event horizon (which lies at a fixed radial location after the residual gauge fixing).
It is also crucial to emphasize that the residual gauge freedom involving the reparametrization of the radial coordinate is a proper diffeomorphism, i.e. it leaves both the boundary metric (which is the flat Milne background (3.1)) and also the T µν of the dual Bjorken flow extracted from holographic renormalization is invariant. It's useful to see this explicitly. For illustration, let's consider the AdS 5 case. The asymptotic expansions take the form: Above, the function a 1 (τ) is related to the residual gauge freedom, and can be chosen arbitrarily. Furthermore, the constraints of Einstein's equations (3.17) impose Using these constraints, one can find via holographic renormalization that Firstly, the above result is exact to all orders in the late proper time expansion. Secondly, we readily find that T µν (τ) is determined by a 4 (τ) alone and is independent of the arbitrary function a 1 (τ) capturing the residual gauge freedom in the asymptotic expansion after utilizing the constraints where a d (τ) is the coefficient of r d in the asymptotic expansions of A(r, τ). Furthermore, extracting a d (τ) from (3.31) we obtain that at the leading order in the late proper time expansion where we have used (3.32), and also (2.15) to define an instantaneous Hawking temperature T (τ) given by Once again comparing with the general (hydrodynamic) late proper time expansion (3.6), we find that For the case of AdS 5 , the identification (3.26) implies that at late proper time For any d, (τ) is given by a perfect fluid expansion given by (3.38) at late time with Thus, at late proper time, the energy density (τ) and the pressures p L (τ) and p T (τ) are thus given by the thermal equation of state obtained from a static black brane geometry, but with a timedependent temperature (3.39) which satisfy the Euler equations.
In order to construct a regular horizon cap, it is useful to change coordinates from r and τ to v and σ following (3.21) as in the case of the spacetime dual to the vacuum state. Note, it follows from (3.29) that v and s are the same. In these new coordinates, the metric (3.27) takes the form: 43) and is dual to Bjorken flow on the Weyl rescaled background metric (3.12) at the boundary. Holographic renormalization and the constraints of the Einstein's equations (3.17) imply that the energymomentum tensor of the dual Bjorken flow takes the form where T W µν takes the general Bjorken form with non-vanishing components given by (3.13) in which and A µν is the Weyl anomaly appearing for even d. Comparing with (3.37), we indeed verify that the bulk diffeomorphism (3.21) implements the Weyl transformation and time reparametrization in the dual theory (with the Weyl factor given by (3.11)) and also reproduces its Weyl anomaly. Particularly, for d = 4, we recall that A µν is simply given by (3.23). The anomalous term is stateindependent (and is always the same as in the Weyl transformed vacuum state). We again note that T µν is invariant under the residual gauge symmetry since it is independent of a 1 (τ) after we implement the gravitational constraints following the previous discussion. Obviously, the late proper time expansion (3.31) takes the form Explicitly, for d = 4, where Above α 1 is the dimensionless parameter associated with the residual gauge freedom. At any order in the late proper time expansion, the terms multiplying α 1 in (3.47) remain the same, however we should replace α 1 by α n at the n-th order. It is also easy to see that g(x) is finite at x = 1 implying that the metric is regular (with no naked singularity) at the perturbative horizon The above metric reproduces the late proper time expansion of (σ) which takes the form (3.14) with 0 given by (3.37) and λ n taking specific values for a given d. Particularly, for any d, we obtain It is easy to verify from (3.7) and the equation of state (see (3.38)) for any d > 2. More details of the perturbative expansion are in Appendix A.

The bulk scalar field and the horizon cap of the Bjorken flow
The key to obtaining the real time correlation functions is solving the dynamics of the scalar field in the gravitational background dual to the Bjorken flow. The starting point, however, is to construct the analogue of the bulk Schwinger-Keldysh contour with the horizon cap for the gravitational background itself. This is straightforward. The metric dual to the the Weyl scaled Bjorken flow (given by Eqs. (3.43) and (3.46)) reaches a constant horizon temperature at late time although the boundary metric has time-dependent spatial components. Since we would be working perturbatively in the late proper time expansion, we will fix the horizon cap at the constant late-time value v = ε −1/d 0 to all orders in the perturbative late proper time expansion while keeping the residual gauge freedom of radial reparametrization unfixed as mentioned above. The metric is analytic to all orders at the horizon cap. Therefore, there is no modification to the metric on the other arm of the bulk spacetime as it is reached via the complexified v contour encircling the horizon as shown in Fig. 2 (as emphasized earlier, there is no analytic continuation in σ and other coordinates). Exactly the same gravitational background is valid on both arms of the complex v contour. It is also easy to see that the on-shell Einstein-Hilbert action on the two arms cancel each other out implying that the dual (non-)equilibrium partition function in absence of additional sources (and with same boundary metric on the two arms) is exactly zero, as should be the case. 16 There is a crucial subtlety to this rather simple construction. We should worry about the residual gauge symmetries at first and higher orders in the late proper time expansion. In what follows, we will show that the analytic behavior of the sourced scalar field retains its equilibrium nature to all orders in the late proper time expansion at the horizon cap, provided the residual gauge symmetry is fixed in a unique way at each order. This addresses an issue which would have arisen if we had fixed the residual gauge symmetries to keep the apparent or the event horizon at v = ε −1/d 0 to all orders in the late proper time expansion. However the location of these two horizons differ at second and higher orders in the proper time expansion. We find that the gauge fixing which implements the regularity of the horizon cap is exactly the same which fixes the event (but not the apparent) horizon at v h = ε 1/d 0 up to third order in the proper time expansion in the case of AdS 4 16 In the field theory, this is equivalent to the statement that W[J 1 = J 2 = 0] = 0. and AdS 5 . Although we do not have an analytic proof that this feature will continue to hold at higher orders, we expect it to be the case as we explicitly find that the gauge fixing is independent of the mass of the bulk scalar field, and it should also hold for fermion, vector and higher rank tensor fields. To see the main advantage of working in the v and σ coordinates, note that the explicit form of the leading order metric (as evident from Eqs. (3.43) and (3.46)) is: A naive ansatz for the scalar field is The above is readily motivated by the Weyl rescaled background metric of the dual state (the boundary metric of the bulk gravitational geometry) given by (3.12) which provides the additional factors multiplying ω, k L and k T (such as (d − 1)/(d − 2) in the front of ω) from √ −g σσ , gζζ and √g ii respectively. However, it turns out that it is more useful to define new momentum variables which depend on σ such that it takes care of the expansion of the longitudinal expansion and transverse contraction. We replace the anstaz (4.2) with a different and a more conventional one shown below: where remarkably there is no explicit σ dependence in f (only implicitly through κ L and κ T ) while in the phase we have the usual fixed conjugate momenta k L and k T . With this new ansatz, we obtain a systematic expansion of the equations of motion in powers of σ −1 . Particularly, the Klein Gordon equation (2 − m 2 )Φ = 0 for the bulk scalar field at the leading order σ 0 is then 17 Remarkably, the left hand side of (4.5) is just the Klein Gordon equation for the massive bulk scalar in the AdS d+1 static black brane geometry (2.11) at temprature β −1 given by (3.51), with v substituting for r, and ω, κ L and κ T identified with the canonical frequency and momenta. The implication is that f will have exactly the same solutions at the leading order as in the static black brane geometry and therefore the same analytic structure at the horizon cap at the leading order. For the homogeneous (k L = k T = 0) and massless case, this result reduces to the observation made by Janik and Peschanski [21] in the context of the transients of the Bjorken flow. The ansatz (4.4) for the scalar field also works after incorporating the viscous and all higher order corrections to the gravitational Bjorken flow background with simple additional modifications that ensure the regularity of the horizon cap after fixing the residual gauge freedom. The ansatz for the bulk scalar field in the full complexified spacetime with 1 and 2 labelling the sheets of the bulk spacetime ending at the forward and backward arms of the Schwinger-Keldysh contour respectively at their boundaries is where Above, in (4.7), a non-analytic term proportional to σ iγ 0 has been introduced with γ 0 being a new dimensionless function of only ω/ε 1/d 0 . This non-analytic factor is crucial to have a regular horizon cap as shown below. However, it does not affect the zeroth order equation of motion which takes the form of the Laplacian on a static black brane as discussed above. For any ω, (4.7) has the appearance of a trans-series in σ −1 with ω characterizing the continuous instanton exponent. Note p and q are functions of ω, k L and k T (and not the redefined momenta κ L and κ T ), since we are integrating over ω, k L and k T , and utilizing the superposition principle to obtain the general solution with right behaviour at the horizon cap. However M n depend on κ L and κ T , and therefore they have non-trivial derivatives w.r.t. σ. The coefficients p and q are determined by imposing Dirichlet boundary conditions at the two boundaries.
The structure of M n in (4.8) is determined as follows. Note that at the zeroth order, i.e. at n = 0, we obtain exactly the solutions of the static black brane as noted above. We can choose the basis of solutions which are in-going and out-going at the horizon, and satisfying the normalization conditions given by (2.14). Explicitly, near the perturbative horizon v = ε 1/d 0 where the horizon cap is located, they take the forms 18 respectively. Note we have used (3.51) to set βω/(2π) = 2ω/(dε 1/d 0 ). For n ≥ 1, φ n,in and φ n,out represent corrections to these zeroth order solutions as discussed below. In (4.8), we have assumed that φ n,in and φ n,out have the same behavior at the horizon cap v = v h = ε −1/d 0 for n ≥ 1as in the case of the equilibrium (or equivalently at the zeroth order). This indeed turns out to be the case with appropriate residual gauge fixing as mentioned before and explicitly shown below. 19 At higher orders, the equations determining φ n,in and φ n,out can be obtained from expanding (2 − m 2 )Φ = 0 in the late proper time expansion after substituting k L and k T in (4.7) by the redefined momenta (4.3), and isolating the σ −n term. We obtain that φ n,in and φ n,out satisfies the inhomogeneous ordinary differential equations Dφ n,in = S n,in , Dφ n,in = S n,out (4.10) with D being the same operator (4.6), which is simply that corresponding to the Klein-Gordon equation for the static black brane at temperature β −1 given by (3.51), at all orders. The sources S n,in and S n,out are functions of v, ω, κ L and κ T , and depend on φ m,in and φ m,out respectively for m < n. Both of these sources also depend on the functions a (i) , k (i) and l (i) , which appear in the late proper time expansion (3.46) of the background metric, with i ≤ n. Note that the source is linear in the bulk field Φ, and therefore splits into S n,in and S n,out at each order in the σ −1 expansion for n ≥ 1. Since φ n,in for n ≥ 1 corrects φ 0,in , we include the particular solution with only φ m,in (and m < n) appearing in the source term S n,in in it. The particular solutions sourced by φ m,out with 0 ≤ m < n which appear in S n,out are added by definition to φ n,out similarly. The regularity of the horizon cap implies that at first and higher orders in the proper time expansion, we should have with γ n,out s and γ n,in s as new functions of ω, κ L and κ T for n ≥ 1. This will ensure that the analytical dependence of the field on v at the horizon cap v = v h = ε −1/d 0 is the same to all orders in the proper time expansion. Note that at the zeroth order, the analogous conditions for φ 0,in and φ 0,out are simply given by 1 on the RHS in both cases by choice (see (4.9)). We will show that the γ n,out s and γ n,in s can be determined uniquely for n ≥ 1 via horizon cap regularity and field theory identities.
Firstly, note that (4.11) implies that φ n,in and φ n,out can be written in the form φ n,in = φ n(p),in − iγ n,in φ 0,in , φ n,out = φ n(p),out − iγ n,out φ 0,out (4.12) for n ≥ 1 such that φ n(p),in and φ n(p),out are the particular solutions of the inhomogeneous ordinary differential equations (4.10). Both φ n(p),in and φ n(p),out are determined by the sources S n,in and S n,out , and are proportional to the coefficients p and q, respectively (i.e. they vanish when p = q = 0). Near the horizon cap v ≈ ε 1/d 0 , they behave as 19 In the following section, we show that preserving the near-horizon behavior to all orders is required for satisfying many consistency conditions, such as ensuring that the retarded correlation function is given by causal response. respectively. Equivalently, where the coefficients should be determined by the equations of motion. For the outgoing solution, we find that this behavior is possible only when the residual gauge parameter α n appearing in the background metric and γ n,out are chosen appropriately to cancel double and single poles appearing in the equation of motion (4.10) at the horizon cap for each n ≥ 1. Furthermore, α n s are simply numerical constants (as they are defined to be), and γ n,out s are linear functions of ω/ε 1/d 0 only. Thus the outgoing solutions appearing in M n in our ansatz (4.7) are determined uniquely. As discussed before, and will be explicitly shown again in the next section this completely determines the advanced propagator of the Bjorken flow. The non-equilibrium retarded propagator then is also determined uniquely, since even out of equilibrium, the advanced and retarded propagators are related by the exchange of the spatial and temporal arguments. Utilizing this, we can determine γ n,in s uniquely as well for n ≥ 1 as will be shown in the next section. In this section, we will focus mainly on the outgoing mode.
It is easy to see that (4.13) (equivalently (4.14)) implies that γ n,in and γ n,out are simply the coefficients of the homogeneous solutions of the equations of motion (4.10) for n ≥ 1. Therefore γ n,in and γ n,out appear in S m,in and S m,out respectively for m > n. We will illustrate by example how requiring the regularity condition (4.13) (equivalently (4.14)) at the n-th order determines γ n−1,out along with the gauge parameter α n (which is a constant) recursively.
Unlike the case of the outgoing mode, the ingoing mode is always analytic at the horizon cap. Therefore, we need to use consistency conditions for the Schwinger-Keldysh correlation functions to determine γ n,in s for n ≥ 1. However, γ 0 in the ansatz (4.7) appears in both the in-going and out-going modes. Our construction passes a significant consistency test that the same function γ 0 (ω/ε 1/d 0 ) determines the homogeneous transients (sourceless solutions which are ingoing at the horizon) with the argument ω/ε 1/d 0 taking values corresponding to the appropriately rescaled quasinormal mode frequencies of the static black hole, as will be discussed in Section 5.5. This is remarkable as we determine the function γ 0 analytically by imposing the regularity condition (4.14) on the outgoing mode at the horizon cap.
In what follows, we illustrate how we determine the outgoing mode and the residual gauge fixing uniquely in the case of AdS 5 . Let us first see how we determine γ 0 and α 1 at the first order in the proper time expansion. This requires the first order correction to the background metric given by (3.47) and (3.48). We find that the equation of motion (4.10) for φ 1,out explicitly takes the following form near the horizon cap v ≈ ε 1/4 0 up to overall proportionality factors: with · · · denoting terms which are regular at v = ε 1/4 0 if φ 1,out is of the form (4.14). See Appendix C for more details. We readily see that in order to have a solution of the desired form (4.13) we must impose so that the double and single pole terms of the equation of motion appearing at the horizon cap vanish. The double pole term determines α 1 and the single pole term determines γ 0 . As claimed before, we find that γ 0 is indeed a simple linear function of ω/ε 1/4 0 (rather just proportional to it) while the gauge parameter α 1 is a numerical constant as it should be.
At the second order in the proper time expansion, similarly α 2 and γ 1 are determined by the vanishing of the double and single pole terms in the equation of motion for φ 2,out respectively. Here we have to utilize the explicit second order correction to the background metric given in Appendix A. Explicitly, We find once again α 2 is just a numerical constant as it should be and γ 1 is a linear function of ω/ε 1/4 0 . It is indeed crucial that the γ n,out s for n ≥ 0 are functions of ω only and are independent of κ L and κ T which depend on σ. Otherwise, the central assumption (4.11) (and thus (4.14)) are not valid for the outgoing mode at the horizon cap, and should be corrected by log terms. The latter would have implied that the behavior near the horizon cap at first and higher orders is different from the zeroth order which is the same as in thermal equilibrium. Note that the ansatz for the ingoing mode in (4.14) remains valid to all orders even if γ n,in s depend on κ L and κ T . The coefficients p n,k in (4.14) also involve derivatives of γ m,in s w.r.t. κ L and κ T with m < n. See Appendix C for more details.
Finally, we have explicitly verified that the values of the residual gauge parameters α 1 and α 2 are such that the event horizon is pinned to the horizon cap v = v h = ε −1/d 0 at the first and second orders respectively for both AdS 4 and AdS 5 . Note that the apparent horizon differs from the event horizon from second order onwards, so the evolving apparent horizon is behind the horizon cap. See Appendix B for details. We expect that this feature persists to all orders so that although the interior of the event horizon is excised, the full double sheeted geometry with the horizon cap still covers the entire bulk regions which can send signals to the boundary. 20 This feature mirrors the causal nature [1] of the Schwinger-Dyson equations for the correlation functions in the field theory. We will discuss more about the consistency of this result in Sec 5.5.
The quadratic on-shell action for the bulk scalar field is the sum of three pieces, namely S in−in and S out−out which are quadratic in the in-going and out-going modes respectively, and the crossterm S in−out . As in the thermal case discussed in Section 2, S in−in = 0. The in-going mode is analytic at the horizon and the contributions from the forward and backward arms of the radial contour cancel out. Once again this is required for consistency, as if we keep only the in-going mode by setting q = 0 in (4.7), then J 1 = J 2 and W[J 1 = J 2 ] = 0 for an arbitrary initial (nonthermal) state. (Recall W is identified with iS on−shell .) The cross term S in−out has a branch point at the horizon cap and the integration over the radial contour results in the two boundary terms like in the thermal case discussed in Section 2. S out−out potentially has a single pole (v h − v) −1 divergence which we denote as S . Explicitly, at the n-th order in the late proper time expansion. We have verified that S = 0 for the solution with the regular behavior at the horizon cap given by (4.11), obtained for the appropriate choices of α n and γ n−1,out as discussed above. One can also check that terms like φ * k,out φ l,out with k + l ≤ 2n also do not contribute to S . However, our arguments for S out−out = 0 for the thermal case in Section 2 do not go through here since they rely on the KMS boundary condition. Nevertheless, since the pole at the horizon vanishes, S out−out is the sum of two boundary terms as well. As shown in the following section, we need to set S out−out = 0 for consistency, particularly to ensure that the (non-equilibrium) retarded correlation function is given by linear causal response. This can be done by simply constraining the initial conditions for the bulk scalar field on the two arms of the bulk geometry such that S out−out = 0 after integrating over the full contour time. We will discuss the physical implications in Section 6.
The upshot is that just as in the thermal case the on-shell action is simply the sum of two boundary terms (including the counter-terms for holographic renormalization) obtained from S in−out . We can then readily differentiate this on-shell action to obtain the Schwinger-Keldysh correlation functions of the Bjorken flow.

Some useful relations and their consequences
Some crucial identities are valid for the Schwinger-Keldysh correlation functions even in out-ofequilibrium states. The first such identity of interest is Similarly, we can arrive at the identity These identities give the retarded and advanced correlation functions in the usual Schwinger-Keldysh basis. The definitions of the correlation functions also imply that in any arbitrary state Finally, it is obvious also that in any arbitrary (non-equilibrium) state We will show that these identities imply the following for the horizon cap of the Bjorken flow. Firstly, will use (5.1) to show that the retarded correlation function is always obtained from the ingoing mode, and (5.2) to show that the advanced correlation function is always obtained from the outgoing mode. Secondly, due to (5.4), we can uniquely determine γ n,in for n ≥ 1. These coefficients cannot be determined by the regularity condition (4.11) of the horizon cap unlike γ n,out as dicussed previously, but together with γ n,out they provide unique solution of the bulk scalar field for specified sources at the two boundaries via our ansatz (4.7). We will also see that by satisfying (5.4) we can also ensure the validity of (5.3). Another important issue is the Weyl transformation. Consider a background metric g µν and its Weyl rescaled version Ω 2 (x)g µν such that Ω(x) is a non-vanishing function. Disregarding Weyl anomaly, the correlations functions of a scalar primary operatorÔ of conformal dimension ∆ O in a conformal field theory living in these two background metrics would be related by Tr(ρÔ(x 1 ) · · ·Ô(x n )) Ω 2 g µν = Ω −∆ O (x 1 ) · · · Ω −∆ O (x n )Tr(ρÔ(x 1 ) · · ·Ô(x n )) g µν (5.5) whereρ = U † ρU with U denoting the unitary operator implementing the Weyl transformation. In the context of the holographic Bjorken flow, the states ρ andρ would be described by the holographic geometries whose boundary metrics are the Milne metric and its Weyl rescaled version respectively, and the respective energy-momentum tensors are also appropriately Weyl transformed including the correct holographic Weyl anomaly. This has been described in detail already in Section 3.2.
We have seen that it is easier to implement the out-of-equilibrium horizon cap prescription in the Weyl transformed geometry (state) in which the temperature and entropy density become a constant at late proper time, and the Klein-Gordon equation in the bulk assumes the form of that in a static black brane at the leading order. It is convenient to compute the correlation functions in this background first, and then go back to the usual Bjorken flow background by Weyl transformation.
In this case, we should use with G(x 1 , x 2 ) being the correlation functions computed after Weyl rescaling to obtain the correlation functions in the usual Bjorken flow background with Ω given by (3.10) and (3.11). We disregard Weyl anomalies since they are not state-dependent (same as in vacuum). Note that the full correlation functions to all orders in the late proper time expansion should depend onζ 1 andζ 2 only through the relative separationζ 1 −ζ 2 in rapidity due to the boost invariance of the background, and similarly only on | x ⊥1 − x ⊥2 | due to translation and rotation symmetries of the background in the transverse spatial plane.

General structure of the hydrodynamic correlation functions
It is useful to define a new variable withγ 0 = γ 0 ε 1/d 0 /ω being a numerical constant (independent of ω). We have already seen that the smoothness of the horizon cap requires that in the case of AdS 5 ,γ 0 = 1/4 (see (4.16)). 21 In terms of this variable, we can readily see from (4.7) that the non-normalizable mode of the scalar field in the gravitational background dual to the Bjorken flow takes the form: and similarly the expectation value of the dual operator (up to state-independent contact terms) is where S n and R n can be defined from the asymptotic expansion of M n : 9) and the labels 1 and 2 stand for the sheets of the bulk spacetime ending on the forward and backward arms of the Schwinger-Keldysh contour respectively. It is also useful to define Clearly (4.8) leads to gives us a(σ, ω, κ L , κ T ), A(σ, ω, κ L , κ T ), b(σ, ω, κ L , κ T ) and B(σ, ω, κ L , κ T ) from their following near boundary expansions It should also be obvious that if we define a n , A n , b n and B n via The correlation function can be extracted simply from the on-shell action, which as shown in the previous section, is the sum of the two boundary terms obtained from S in−out . The computations are similar to the case of the thermal equilibrium discussed before. The general structure of the correlation function with {a, b} = {1, 2} standing for the forward and backward arms of the Schwinger-Keldysh contour (using (3.15)) is (see Appendix D for more details) where (compare with the thermal case (2.31) -recall σ 3 is the Pauli matrix) .
As mentioned, we will show in Section 5.4 that the above can be satisfied by appropriate choices of γ n,in s at each order in the late proper time expansion where we go to large values of the average reparametrized proper time σ = (σ 1 + σ 2 )/2 with fixed difference σ r = σ 1 − σ 2 . Thus we see that the correlation functions of the hydrodynamic Bjorken flow has a hidden and simple bi-local thermal structure to all orders in the late proper time expansion. We will show this satisfies crucial consistency tests. Using (5.19) and (5.20), and the identities (5.1) and (5.2) which are valid out-of-equilibrium, we readily see that the actual retarded and advanced propagators of the Weyl transformed Bjorken flow is Therefore, it follows from (5.20) that indeed the retarded propagator is obtained purely from the out-of-equilibrium ingoing mode and the advanced propagator is also obtained from the out-ofequilibrium outgoing mode as claimed before. Note that at the zeroth order these results are automatic as it reduces to the thermal case described earlier. 22 Finally, to obtain the Schwinger-Keldysh correlation function of the Bjorken flow, we should implement the Weyl transformation. Using (3.11), the Weyl transformation finally yields the correlation function of the Bjorken flow:

In the limit of the perfect fluid expansion
At very late proper time, the Bjorken flow is simply a perfect fluid expansion. For the Weyl transformed Bjorken flow which reaches a specific final temperature, the scalar field in the dual gravitational geometry can be mapped to that in the static thermal black brane spacetime at the leading order in the late proper time expansion. This naturally implies that the Schwinger-Keldysh correlation functions at late proper time should be related to the corresponding thermal correlaton functions via appropriate spacetime-reparametrizations. After undoing the Weyl transformation, we should obtain the Schwinger-Keldysh correlation functions of the perfect fluid expansion.

(5.26)
As shown in the previous section, both φ in and φ out assume the form in the static black brane geometry, and therefore given by the corresponding static black brane results and therefore all Schwinger-Keldysh correlation functions after the momentum integrals shown in (5.16) should assume the thermal form, i.e.
with G β denoting the thermal correlation functions. In order to obtain the above form, we change variables in the momentum integrals from k L and k T to κ L and κ T . Note the Jacobian of this transformation is identity. We also need to change the integration variable ω to ω = ω(d −1)/(d −2) which yields a Jacobian (d − 2)/(d − 1). Furthermore, we have used The factor (d − 2) 2 /(d − 1) 2 in (5.16) is cancelled by one factor of (d − 1)/(d − 2) each from the Jacobian and s . Finally, the Weyl transformation in the late proper time limit gives the Schwinger-Keldysh correlation functions of the perfect fluid expansion. Since the static black brane has full rotational invariance in terms of the reparametrized spacetime coordinates, we may also write for the correlation functions in the limit of the late time perfect fluid expansion.
It is obvious that we are actually resumming over the associated σ factors in the spatial factorŝ ζ 1 −ζ 2 and x ⊥1 − x ⊥2 . For brevity we will denote the variables as σ r = σ 1 − σ 2 (as defined before) and The late proper time expansion of the Schwinger-Keldysh correlation function then amounts to the following series with · · · denoting trans-series type completion about which we will discuss more below.
Since the choice of τ 0 determines the effective final temperature T = β −1 via (3.51), let's consider a scaling τ 0 → ξτ 0 to see if there is any ambiguity in the above result. Note that G β in (5.31) depends only on T σ r and T ζ 2 r + x ⊥ 2 r as it is a thermal correlation matrix of a CFT. Under this scaling, it is evident from (1.8) that σ → ξ 1 d−1 σ, and therefore σ r → ξ T as discussed earlier in Section 3.1. Together these imply that G β and also G n in (5.33) are invariant under τ 0 → ξτ 0 since T σ r and T ζ r and T x ⊥r are invariant this scaling. However, the Weyl factor in (1.11) scales as ξ implying that the dimensionless correlation function σ −2∆ 0 G is invariant under the scaling of τ 0 .

First and higher orders in the late proper time expansion
From the form of the Schwinger-Keldysh correlation functions given by Eqs. (5.16)-(5.20) which are valid to all orders, we can readily go beyond the perfect fluid limit systematically and construct the late proper time expansion (5.33). The contribution to the first order correction comes from the following terms: 1. The phase factor, 3. Late-time expansion of the matrices R and S given by (5.11) in σ −1 . This also includes the κ L and κ T dependence. For instance, consider any generic function f (κ L1 ) and f ( κ T 1 ), and similarly for f (κ L2 ) and f ( κ T 2 ).
From (5.17), we obtain An easier way to read off the first order correction is to use (5.19) and (5.20), and consider similar σ −1 expansion of G R and G A . Collecting first order terms from (5.20), as for instance, we obtain and This allows us to finally determine γ 1,in and similarly γ n,in for n > 1 as follows. After integrating over the frequency and momentum integrals (i.e. doing integrals over ω, k L and k T ) as in (5.23), we obtain the first order correction in the series expansion of the type (5.33): 40) and similarly for G A (x 1 , x 2 ). For the identities (5.1) and (5.2) to be satisfied, we therefore need (note that x ⊥r is invariant under at each n, which can be ensured by appropriate choices of γ n,in (ω, κ L , κ T ) at each order for n ≥ 1, since they determine both A n and a n . In turn, this justifies the bi-local thermal structure (5.19) and (5.20) which implies via (5.1) and (5.2) that the retarded correlation originates purely from the ingoing mode and the advanced correlation function purely from the outgoing mode at all orders. We emphasize that although γ n,out s are only functions of the frequency ω, γ n,in s are functions of both frequency and momenta. With γ n,in determined we know the RHS of (5.38) completely. Doing the frequency and momenta integrals shown in (5.16), we can obtain the first and similarly higher order corrections to all Schwinger-Keldysh propagators.

Consistency checks
The general result for the hidden bi-local thermal structure of the hydrodynamic correlation functions given by Eqs. (5.16)-(5.20) to all orders satisfies a simple consistency check. When J 1 = J 2 , i.e. when the sources are the same at the two boundaries, we have only the in-going solution which is analytic. The on-shell solution also vanishes. When J 1 = J 2 , we have on the forward arm, and thus we see that the response in the forward arm is indeed given by the retarded correlation function, and this also follows from the bi-local thermal structure of the hydrodynamic correlation functions resulting in (5.23) as shown above. We can also see that the analytic in-going mode is indeed always related to causal response from the full non-linear evolution (we just consider the forward arm here). There is a unique regular solution to Einstein's gravity minimally coupled to a scalar field corresponding to an initial condition for the bulk scalar field Φ(v, ζ, x ⊥ ) at a constant Eddington-Finkelstein time hypersurface when the boundary source J(σ, ζ, x ⊥ ) is specified for all time in the future. This solution can be obtained via a Chebyshev grid in the radial direction [9] which implies analyticity at the horizon, and thus the causal evolution is built out of analytic ingoing modes. Our result then should follow in the linearized limit. The way to implement this explicitly was shown earlier in [10] which is thus reproduced by the horizon cap prescription. 23 A more non-trivial consistency check involves the outgoing mode. Consider an outgoing mode with frequency ω and momenta k L and k T in the Bjorken flow background. Evidently, from the discussion before, to all orders in the late proper time expansion, we have J 2 (ω, k L , k T ) = e βω J 1 (ω, k L , k T ) with β defined by ε 0 as in (3.51). Therefore, on the forward arm (where we choose to do the k L and k T integrations first) Above, we have used the result from (5.19) that G 12 e βω = G 21 and also the identity (5.2). Thus indeed we see that the outgoing mode gives advanced response. Note for both J 2 = e βω J 1 and G 12 e βω = G 21 (used above) to hold, we need absence of log terms at the horizon, which follows if γ n,out s are functions of ω only. The out-going solution which satisfies the regularity condition (4.14) at the horizon cap indeed has this property with γ n,out determined from the cancellation of the double and single poles of the equation of motion being a linear function of ω and independent of κ L and κ T . Furthermore, let us consider transients which are ingoing solutions with vanishing sources at the boundary. To be simplistic, let us consider the homogeneous transients first. Since we can map the leading order solution to the static black brane background, we should have where ω Q corresponds to the homogeneous quasinormal mode frequencies of the static black brane, as noted by Janik and Peschanski. A non-trivial consistency test of our ansatz for the bulk scalar field is that the condition for the sources to vanish at the first subleading order is The details of the numerical verification are shown in Appendix E. These transients can be added to the solution with specific sources since the equation of the scalar field is linear and are needed for matching with arbitrary initial conditions for the bulk scalar field, as discussed below. The case of the inhomogeneous transients is complicated since the effective momenta κ L and κ T also depend on the proper time, and will be examined in a later work.
Our key result that the horizon cap is pinned to the non-equilibrium event horizon is consistent with the feature in non-equilibrium quantum field theory that the evolution of the Schwinger-Keldysh correlation functions, which can be written in the form of the Schwinger-Dyson equations for the commutator and anti-commutator via functional derivatives of the two-particle irreducible effective action, is causal and uniquely determined once we give the initial conditions for these in the initial state [1]. The bulk analogue is that the evolution of the field configuration on the full complexified geometry with the horizon cap should be uniquely and causally determined for given intial conditions. For this to hold, the horizon cap should cover the spacetime outside the event horizon, which is indeed the case.

On initial conditions and seeing behind the event horizon
The discussion on initial conditions is subtle. 24 We only discuss this briefly here and hope to address it fully in the future. For the full complexified geometry, in order to obtain a unique solution for the bulk scalar field on both arms, we need to specify Φ 1 (v, ζ, x ⊥ ) and Φ 2 (v, ζ, x ⊥ ) at an initial time σ 0 on both arms, and also the sources J 1 (σ, ζ, x ⊥ ) and J 2 (σ, ζ, x ⊥ ) at the two boundaries for all times in the future. Our horizon cap prescription leads to unique solutions for the in-going and out-going modes corresponding to J 1 (σ, ζ, x ⊥ ) and J 2 (σ, ζ, x ⊥ ) in a large proper time 24 We thank Kostas Skenderis for discussions on these issues. expansion, however it is not guaranteed we can match with arbitrary initial conditions Φ 1 (v, ζ, x ⊥ ) and Φ 2 (v, ζ, x ⊥ ) at σ = σ 0 . This is a conceptual problem although the matching with data for the initial density matrix is a separate issue from obtaining the correlation functions in the limit where the state hydrodynamizes and forgets most details of the initial conditions; we have concerned ourselves with the latter here. The two arms of the Schwinger-Keldysh contour specify the two in-states which have overlap, in principle, with arbitrary bulk field configurations at initial time, corresponding to Φ 1 (v, ζ, x ⊥ ) and Φ 2 (v, ζ, x ⊥ ) respectively. The full on-shell gravitational action beyond the hydrodynamic limit should yield the matrix element with T c denotes time-ordering in the closed time contour, and Φ 1 | and Φ 2 | denoting states in the dual field theory corresponding to the semi-classical bulk field configurations. 25 This issue can be partly addressed by allowing for the ingoing transients discussed above. These do not change the sources at the two boundaries but they modify the initial conditions at the two slices. However, this is not enough because we have a pair of initial conditions, one each for each arm. The transients are analytic at the horizon cap and do not affect the initial conditions on the two arms independently. Furthermore, we discussed that we need S out−out = 0 for consistency and this also constrains possible initial conditions on the two slices. Although this needs to be investigate further, presently we may conclude that one may not be able to obtain semiclassical solutions corresponding to arbitrary in-in states meaning that in the semi-classical gravity approximation decoherence which suppresses some off-diagonal matrix elements is built in. This issue is also relevant for the general approach of Skenderis and Balt van Rees [18]. We also note that there is possibility of adding other semi-classical complex saddles of the gravitational action which do not have a natural hydrodynamic limit.
There is another independent route for matching with initial conditions. The late proper time expansion of the correlation functions (5.33) is divergent, and would require a trans-series completion which would naively be of the form (like multi-instanton series) hydrodynamic relaxation of the Green's function for fixed separations of the reparametrized spacetime arguments x 1 and x 2 . It is possible that these functions are approximately constants. Crucially the Stokes data C α (σ r , ζ r , x ⊥r ) are determined by initial conditions. Instead of being constants as in the case of (τ), the Stokes data are now functions. The mathematical formulation of the trans-series is required to formulate a precise way to capture information about the initial state via Stokes data, and therefore the quantum fluctuations behind the event horizon. In order to match with the initial conditions, we should add the transients and therefore the Stokes data to hydrodynamic gravitational background as well. Furthermore, additional Stokes data for the evolving event horizon which can be captured by the time-dependent residual gauge transformation is needed too (recall that the residual gauge transformation itself is expressed in the large proper time expansion for which a trans-series completion is necessary like for (σ)). One could interpret the latter feature as the horizon hair (more precisely a proper residual gauge transformation) playing a crucial role in decoding the interior of the event horizon. 26

Conclusions and outlook
We believe that our method for computing the real time correlation functions of the hydrodynamic Bjorken flow can be generalized to generic situations where the state hydrodynamizes, meaning that the dynamics of one-point functions can be captured by an asymptotic series expansion which is generated by the hydrodynamic evolution of the temperature and velocity fields. The key would be to see if there exists Weyl transformations with non-trivial spacetime dependence which can map the flow at late time to a configuration with constant temperature and entropy density although time translation symmetry may not be present. In this case, as demonstrated here for the Bjorken flow, the dual black hole would attain a horizon with constant surface gravity and area at late time, although a time-like Killing vector may be eternally absent. The event horizon's shape would fluctuate in space and time keeping the total area fixed when entropy production ceases in this limit. Nevertheless, with appropriate spacetime reparametrizations, the horizon cap prescription can be made to work as in the case of the Bjorken flow provided at the leading order in the large (suitably reparametrized) time expansion the equation of motion of bulk fields can be mapped to that on a static black brane. We will pursue this direction in the future, particularly the Gubser flow [52].
Another important problem would be to understand the Borel resummation of the hydrodynamic series (5.33) and compute the Schwinger-Keldysh correlation functions of the holographic hydrodynamic attractor [20]. Furthermore, the late time thermal nature of the correlation functions in terms of the reparametrized spacetime arguments could be of phenomenological relevance in heavy ion collisions especially with regards to the dynamics of heavy quarks and bound states, and jets in the expanding quark-gluon plasma [53].
Our prescription can be used to construct the quantum generalization of classical stochastic hydrodynamics (see [54,55] for reviews) in holographic theories by considering backreaction of the fluctuations on the background geometry systematically. This is especially important in context of superfluid fluctuations since quantum dynamics is important at coherence time and length scales which are shorter than the scattering time and the mean free path respectively (see [14] for an excellent related discussion)-non-linearities can potentially cause novel non-trivial effects such as quantum corrections to the long time tails. 27 More generally, we would like use the horizon cap method to study holographic (evaporating) black holes interacting with heat baths or dynamical reservoirs, and understand the reconstruction of the islands [51,57] (which include the black hole interior) from Hawking quanta. In such cases, semi-holographic formulations for open quantum systems (see [58,59] and also [60] for instance) can provide useful models. 27 For a recent work on the role of classical fluctuations of the superfluid order parameter in modifying transport properties see [56].

Acknowledgments
We thank Kostas Skenderis and Balt van Rees for helpful discussions and sharing many insights on the issues discussed in Section 6. AM acknowledges support from the Ramanujan Fellowship and Early Career Research Award of the Science and Engineering Board of Department of Science and Technology of India, new faculty seed grant of IIT Madras, the institute of excellence schemes of the Ministry of Education of India and also IFCPAR/CEFIPRA funded project no. 6304-3. The research of AB is supported by IFCPAR/CEFIPRA funded project no. 6304-3.
A Details of the five dimensional metric dual to the Bjorken flow upto second order in the late-time expansion Consider the late-time expansion of Einstein's equations based on the ansatz (3.43) along with (3.46). At zeroth order, the six components of Einstein equations read, We solve these equations based on the conditions which is motivated from the fact that at very late times, we have an AdS 5 black hole with constant horizon area and surface gravity. Once these initial conditions are imposed, then at each subleading order, the six equations can be repackaged to the following three equations where S 1,i , S 2,i and S 3,i are the sources at the i-th order which depend on a j , l j , k j and their derivatives, for j < i. Now at each subleading order, the most general solution to these equations are given by, where the Greek letters with subscript i correspond to the integration constants at i-th order and PI 1,i , PI 2,i and PI 3,i are the particular solutions determined by the sources S 1,i , S 2,i and S 3,i respectively. The expressions for the PIs for i = 1 are simple and already mentioned in (3.47)-(3.48). For i = 2, we get − 16 arctan (x) log 1 + e 2i arctan (x) + 8 arctan (x) log 1 − ie 2i arctan (x) + 2π log 1 − ie 2i arctan (x) where Li 2 (x) is the polylogarithmic function of order 2. The expressions for the particular solution get increasingly complicated at higher orders. However, the simple structure of the homogeneous solutions remain the same at every order. The integration constants associated with the homogeneous solutions can be fixed at every order in the following way: • The fixation of α i has been detailed in (4.15) − (4.17). The same pattern repeats for the rest of α i>2 .
• The integration constants β i are fixed by mapping the coefficients of x 4 terms in (A.3) to a 4 (τ), l 4 (τ) and k 4 (τ) of (3.33) (or a 4 (σ), l 4 (σ) and k 4 (σ) in the Weyl transformed coordinate) and then solving the constraint equations (3.34). geometry given by (3.43), consider radial null geodesics in this spacetime whose equation is given by, For a static geometry, v E = cons. and the event horizon is simply given by the zero of g σσ . However, in a dynamical geometry, the event horizon will depend on σ. Consider the late-time expansion of the horizon of the form, Now to determine v Ei , we use the equation (B.1) along with the known late-time expansion of A(v, σ), At leading order, we have where α 1 and α 2 are the residual gauge parameters associated with the late-time expansion of A(v, σ). Now the gauge parameters can be uniquely fixed by demanding regularity of the horizon cap. In d = 4, this gives (recall (4.16) and (4.17)) α 1 = −3, These values, in turn fix the event horizon to 1/ε 1/d 0 upto the second subleading order. Thus, regularization of the horizon cap leads to vanishing of the sub-leading corrections v Ei (for i > 0) to the static event horizon. We expect this feature to remain true to all orders in the late-time expansion.

B.2 Apparent horizon:
The apparent horizon is a null hypersurface that acts as a boundary between the trapped and untrapped regions, whose location is given by the product of the expansion parameters θ ± , i.e. Θ = e f θ + θ − , where the factor e f is defined below. The trapped region is characterized by Θ > 0, where the light rays directed outward propagates inward, whereas in the un-trapped region the light rays directed outward propagates outward and is characterized by Θ < 0. Therefore Θ = 0 gives the location of the apparent horizon.
Here we will adopt the dual-null formalism [61][62][63][64] to study location of the apparent horizon, where one defines a pair of null hypersurfaces Σ ± parameterised by scalars ζ ± , with associated one-forms n ± = −dζ ± . The null normal vector associated with these hypersurfaces are given by, where e f is the normalization factor given by, e f = −g µν n − µ n + ν .
Next we define the expansion parameters θ ± , where µ is the spatial volume element of the geometry in which the hypersurfaces are defined and L ± is the lie derivative along the null normal vectors, l ± . Finally we introduce the invariant 29 quantity Θ = e f θ + θ − , whose zero gives the location of the apparent horizon.
In case of our geometry (3.43), we consider a pair of null hypersurface defined by constant values of the retarded and advanced radial null coordinates, whose normal one-forms are given by, where N + and N − are the overall normalizability factor that can be determined by the integrability condition d(dn ± ) = 0. This ensures that the one-forms (B.6) and (B.7) are exact. However there is no need of computing them explicitly, as when we compute Θ the contribution of N + and N − from e f will cancel with the ones coming from θ − and θ + . Now since the apparent horizon depends on σ, one can consider a late-time expansion of the apparent horizon similar to the event horizon as, To determine the coefficients v Ai we solve Θ = 0 at every order. In case of any spacetime dimensions d + 1, the leading order result is universal and turns out to be v A0 = 1/ε 1/d 0 , which coincides with the event horizon. However, the sub-leading corrections are dimension dependent. For d = 4, the first and second sub-leading corrections to the leading terms are Note that, for α 1 = −3 the first order correction to the apparent horizon vanishes similar to the case of the event horizon. However, at second order, the value of α 2 (given by (4.17)) for which the the correction to the event horizon vanishes, now renders v A2 > 0. So the apparent horizon will lie inside the event horizon (see Figure 4). Again, we expect this feature to remain true at higher orders as well.  Here ε 0 is set to 1 and the gauge parameters α 1 and α 2 are so chosen that the event horizon is fixed to 1. The apparent horizon always lies inside the event horizon (at greater v).
C More details of the bulk scalar field in the AdS 5 Bjorken flow background In this appendix, we will provide brief details of the solution of Klein-Gordon equation at the leading (4.5) and first subleading order (4.10) (for n = 1) and validate the near-horizon behaviours (4.9) and (4.13) for the same. To be specific, we will consider the example of d = 4.
At leading order, the homogenous Klein-Gordon equation (4.5) simply takes the form of that in a static black brane geometry. So we can expand the solution in the standard basis provided by the ingoing and outgoing modes. Near the horizon, these modes admit the following expansions Now, without loss of generality, we can set p 0 = 1 and q 0 = 1. Then, in the limit v → ε −1/4 0 we reproduce the behaviours given in (4.9).
E Transients and γ 0 The quasinormal frequencies and transients capture the causal response of a black hole to small perturbations. Typically, for a black hole these are determined using the spectral representation method. However, here we adapt a simpler approach for the computation of the same. For the homogeneous transients, consider the field ansatz (4.7) with only the ingoing modes of a massless scalar, i.e.
where φ n,in (v, ω) obeys the same equations (4.5) and (4.10) (and hence admit the solutions (4.9) and (4.12)) with κ L , κ T = 0. The residual gauge parameters α i appearing in the subleading solutions are fixed such that at every order, the event horizon is pinned at v h = ε −1/4 0 with associated Hawking temperature T = ε 1/4 0 /π. For concreteness, we will again consider d = 4. Recall that, at leading order, the ingoing mode(with zero spatial momentum) admits the the near-horizon expansion To compute the transients, we introduce the dimensionless decay rate λ = ω/πT and also scale the horizon to set v h = 1. In terms of the dimensionless parameter, the solution (E.2) when expanded near the boundary v = 0 reads, The QNM corresponds to a 0 (λ qnm ) = 0, implying that there is no contribution to the source from the leading order field. The lowest QNM frequency obtained in this method turns out to be λ qnm = −2.7668i ± 3.11945 , which agrees with the one obtained from spectral representation method upto order 10 −10 .
Similarly at subleading orders n > 0, the sources can be made to vanish, i.e. a n>0 (λ qnm ) = 0 by suitably fixing γ 0 (for n = 1) and γ n,in (for n > 1) as functions of λ qnm . For example, corresponding to this lowest λ qnm , γ 0 (λ qnm ) obtained by setting a 1 (λ qnm ) = 0 at the first subleading order is, γ 0 (λ qnm ) = −0.68669i ± 0.779863 , (E. 4) which turns out to be same as the one computed from the outgoing mode analysis (4.16), These two results agree upto order 10 −9 . See [65] for a recent relevant work in the large D limit.