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

Avik Banerjee, Toshali Mitra, Ayan Mukhopadhyay

► To cite this version:

Avik Banerjee, Toshali Mitra, Ayan Mukhopadhyay. Correlation functions of the Bjorken flow in the holographic Schwinger-Keldysh approach. Physical Review Research, 2023, 5 (4), pp.043230. 10.1103/physrevresearch.5.043230. hal-04783917

HAL Id: hal-04783917

https://hal.science/hal-04783917v1

Submitted on 14 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Correlation functions of the Bjorken flow in the holographic Schwinger-Keldysh approach

Avik Banerjee*

Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université,
Université de Paris, F-75005 Paris, France
and Center for Strings, Gravitation and Cosmology, Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India

Toshali Mitra

Asia Pacific Center for Theoretical Physics, Pohang 37673, Korea;
The Institute of Mathematical Sciences, Chennai 600113, India;
and Homi Bhabha National Institute, Training School Complex, Anushakti Nagar, Mumbai 400094, India

Ayan Mukhopadhyay

Instituto de Física, Pontificia Universidad Católica de Valparaíso, Avenida Universidad 330, Valparaíso, Chile
and Center for Strings, Gravitation and Cosmology, Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India

<figure><img src="image_2.png" /></figure>

(Received 22 April 2023; accepted 15 November 2023; published 11 December 2023)

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 expand 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.

DOI: 10.1103/PhysRevResearch.5.043230

## I. INTRODUCTION

### A. Motivation and aims

One of the outstanding issues in the holographic duality, which 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 (KMS) periodicity at an appropriate temperature, and how the occupation numbers of quasiparticles equilibrate or evolve to new fixed points [1–4]. These correlation functions are actually indispensable for understanding dynamics far from equilibrium in the strong-coupling limit where the system cannot be described by quasiparticles. 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–8] providing a primary example and numerical relativity [9] providing a powerful tool, an explicit computation of a

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>* avikphys02@gmail.com</footnote>

<footnote>† toshali.mitra@apctp.org</footnote>

<footnote>‡ ayan.mukhopadhyay@pucv.cl</footnote>
generating functional for hydrodynamic Schwinger-Keldysh (SK) correlation functions has not been achieved yet.<fnref n="1" />

The object of interest is the generating functional,

$$ \exp(iW[J_1, J_2]) = \mathrm{Tr} \left( \hat{\rho}_0 T_c \exp \left( -i \oint dt \int d^{d-1}x J(t, \mathbf{x}) \hat{O}(t, \mathbf{x}) \right) \right), \quad (1) $$

in the dual field theory, where $\hat{\rho}_0$ denotes the initial density matrix, $\oint$ the closed time Schwinger-Keldysh contour composed of the forward and backward arms—$\oint = \int_{-\infty}^{\infty} + \int_{\infty}^{\infty}$, where the source $J(t, \mathbf{x})$ is specified such that it is $J_1(t, \mathbf{x})$ and $J_2(t, \mathbf{x})$ on the forward and backward arms of the contour, respectively, and $T_c$ denotes contour ordering. Formally, we can rewrite the density matrix as

$$ \hat{\rho}_0 = \int [\mathcal{D}\phi_1][\mathcal{D}\phi_2]\rho_0(\phi_1, \phi_2)|\phi_1\rangle\langle\phi_2| \quad (2) $$

in terms of a basis $|\phi\rangle$ of field configurations and construct the kernel

$$ \mathcal{K}(\phi_1, \phi_2; J_1, J_2) = \langle \phi_2 | T_c \exp(-i \oint dt \int d^{d-1}x J(t, \mathbf{x}) \hat{O}(t, \mathbf{x})) | \phi_1 \rangle. \quad (3) $$

Then the functional $W[J_1, J_2]$ can be obtained from

$$ \exp(iW[J_1, J_2]) = \int [\mathcal{D}\phi_1][\mathcal{D}\phi_2]\rho_0(\phi_1, \phi_2)\mathcal{K}(\phi_1, \phi_2; J_1, J_2). \quad (4) $$

When $\hat{\rho}_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 of the evolution of an expanding system on the forward light cone (see Sec. III A 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_{\tau\tau} = \epsilon(\tau)$ via Ward identities, where $\tau = (t^2 - z^2)^{1/2}$ is the proper time of an observer comoving with the flow and $z$ being the longitudinal coordinate along which the expansion happens. (Here we use units where $c = 1$.) At late time, $\epsilon(\tau)$ is described by

hydrodynamics and thus it reaches a perfect fluid expansion, so that

$$ \epsilon(\tau) \approx \epsilon_0 \left(\frac{\tau_0}{\tau}\right)^{\frac{d}{d-1}}. \quad (5) $$

In the hydrodynamic regime, it can be described by a single constant parameter, namely

$$ \mu := \epsilon_0 \tau_0^{\frac{d}{d-1}}. \quad (6) $$

The full hydrodynamic series for $\epsilon(\tau)$ in powers of $\tau^{-\frac{d-2}{d-1}}$ (essentially a derivative expansion) is given in terms of $\mu$ (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^\mu(t, \mathbf{x})$, the energy density $\epsilon(t, \mathbf{x})$ [or equivalently the temperature $T(t, \mathbf{x})$], etc., and we would not require the knowledge of the detailed (off-diagonal) matrix elements $\rho(\phi_1, \phi_2)$ of the state or the kernel $\mathcal{K}(\phi_1, \phi_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 space-time configurations of conserved currents is also known as the large deviation functional [13] which can be computed in many models studied in classical nonequilibrium statistical mechanics.<fnref n="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 also present concrete steps for understanding how to go beyond the hydrodynamic limit and recover the initial state.

## B. A brief historical review

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

<footnote n="1">¹A limited number of observables can still be computed analytically or numerically. The equal-time two-point functions can be computed via the geodesic approximation when the operator has a large scaling dimension, even out of equilibrium. The out-of-equilibrium retarded correlation function can also be computed by implementing linear causal response appropriately—see Ref. [10] for a general prescription. Furthermore, equal-time Green’s function can be computed for operators with large anomalous scaling dimensions in the geodesic approximation and has been used to understand thermalization [11].</footnote>

<footnote n="2">²The large deviation functional gives the probability for a macroscopic space-time profile of a conserved charge or current density which does not necessarily satisfy the hydrodynamic equations. In a quantum system, the off-diagonal matrix elements in the basis of macroscopic field configurations for the conserved charges and currents eventually decohere, but the decoherence would be of interest. See Ref. [14] for a recent discussion on the possibility of a quantum generalization of large-deviation functional methods.</footnote>
to compute the out-of-equilibrium retarded correlation func-
tion 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\beta/2$ (note $\beta = 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 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 VI.

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 $\epsilon$ (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$ specified independently at the two boundaries.

C. Summary of the results

Our first result is the demonstration that the CGL horizon cap prescription reproduces both the KMS periodicity and the ingoing boundary-condition prescription of Son and Starinets (needed to obtain the retarded correlation function from causal response) in thermal equilibrium. Earlier in Ref. [19], both of these were demonstrated only up to quadratic orders in the frequency. Our proof relies on a simple matrix factorization of thermal correlation functions which is reproduced by the holographic method for extracting the correlation function of a scalar operator in the field theory via the on-shell action of the dual bulk scalar field whose mass is determined by the scaling dimension of the operator.<fnref n="3" />

Our primary tool for extending the CGL prescription into
the (conformal) hydrodynamic Bjorken regime is Weyl rescal-

ing that maps the asymptotic perfect fluid expansion to a
flow in which the temperature and entropy density become
*constant* at late time in a nontrivial background metric *without any timelike Killing vector*. This Weyl transformation can be
lifted to a bulk diffeomorphism [20,21] in which the dual
black hole's event and apparent horizon coincide at a fixed
radial location at a very late proper time. The area and surface
gravity of the horizons, and therefore the entropy density and
the temperature of the dual fluid, remain constant at late time,
but the horizons shrink in the directions transverse to the
flow and expand in the longitudinal direction. In the absence
of any timelike Killing vector, the late-time behavior of the
Weyl transformed Bjorken flow is *not* thermal although the
temperature and energy density become constant.

The requirement that the dual geometry has a regular future
black-hole horizon necessitates viscous corrections to the per-
fect fluid Bjorken flow, which implies the same for the Weyl
transformed version mentioned above since the bulk dual of
the latter is obtained simply via bulk diffeomorphism. Fur-
thermore, bulk regularity determines the precise values of *all*
transport coefficients of the holographic theory order by order
in the derivative expansion (equivalent to the large proper time
expansion) with the Bjorken flow giving a special case of
the fluid-gravity correspondence [6]. The large proper time
expansion of the holographic Bjorken flow has been worked
out to very high orders [22]. We discuss the bulk dual of the
Bjorken flow and its Weyl transformed version in Sec. III in
detail.

We want to emphasize that the Weyl rescaling (and hence
the dual bulk diffeomorphism) is determined purely by the late
perfect fluid flow regime that is mapped to that of a constant
temperature in a nontrivial background metric in the field
theory (which is also the boundary metric of the dual black
brane geometry after the bulk diffeomorphism) as mentioned
above. The Weyl transformation and consequently the Weyl
rescaled background metric of the field theory do *not* receive
any correction at first and higher orders in the derivative (i.e.,
large proper time) expansion.

As a result of the bulk diffeomorphism which implements
the Weyl transformation of the dual Bjorken flow, the Klein-
Gordon equation for a bulk scalar field with arbitrary mass can
be mapped to that of a static black brane at late proper time. In
this map, the frequency needs to be appropriately scaled, and
the momenta in the static black brane geometry are identified
with *comoving* longitudinal and transverse momenta.<fnref n="4" /> For
vanishing momenta, this reproduces the result of Janik and
Peschanski [23].

The utility of this map to the static black brane is to es-
tablish the horizon cap prescription in the asymptotic perfect
fluid limit. It follows that the Schwinger-Keldysh correlation
functions of the operator dual to the bulk scalar field can be
mapped to thermal correlation functions after suitable space-
time reparametrizations in the asymptotic perfect fluid limit

<footnote n="3">³The most important ingredient for realization of the matrix factor-
ization is that only the terms which are obtained from the product
of an ingoing mode with an outgoing mode can contribute to the
quadratic on-shell action of the free bulk scalar field minimally
coupled to gravity.</footnote>

<footnote n="4">⁴Note that the map to the static black brane holds only at the leading order in the large proper time expansion. The comoving momenta, however, are also defined at the leading order itself via the Weyl transformed background metric of the field theory (boundary metric of the bulk geometry).</footnote>
of the Bjorken flow (i.e., when both the proper time arguments of the two-point correlation functions are sufficiently large). The nonthermal nature of the Schwinger-Keldysh correlation functions can thus be absorbed into these spacetime reparametrizations up to overall (proper time-dependent) Weyl factors.

We also show that the horizon cap prescription extends to all orders in the late proper time expansion. First, at first and higher orders, the Klein-Gordon equation of the bulk scalar takes the same form as at the leading order but with source terms, and can be systematically solved at each order in the large proper time expansion such that the leading behavior of the ingoing and outgoing modes are exactly the same as at the zeroth order (that can be mapped to the static black brane) implying that the standard analytic continuation at the horizon cap can be performed at each order. Remarkably, this is possible only if the horizon cap is pinned to the time-dependent event horizon (and not the apparent or other dynamical horizons which do not coincide with the event horizon beyond the zeroth order).

Second, the requirement of the near-horizon behavior of the ingoing and outgoing modes at the time-dependent event horizon does not completely determine all the order-by-order corrections. The undetermined coefficients are those that give the leading behavior of the ingoing modes at the horizon at first and higher orders in the proper time expansion. However, when the on-shell action of the scalar field is used to extract the Schwinger-Keldysh correlation functions, we find that they are consistent with all field theory identities provided one of these identities which hold for arbitrary nonequilibrium states provided one of these identities is used to determine the unfixed coefficients that give the leading near-horizon behavior of the ingoing modes.

Therefore, we establish the horizon cap prescription can be unambiguously extended to determine the Schwinger-Keldysh correlation functions of the entire hydrodynamic tail of the Bjorken flow at any arbitrary order in the large proper time expansion. Aside from satisfying all field theory identities, we find that the retarded propagator given by the ingoing modes reproduces the normalizable bulk solutions at complex frequencies (which map to quasinormal modes of the static black brane [23]) with vanishing sources to all orders although the relevant phase factors at first and higher orders are determined as functions of the frequency via the near-horizon behavior of the outgoing modes at real frequencies. The latter implies nontrivially that one can obtain the normalizable bulk solutions dual to the transients (nonequilibrium generalization of the collective excitations of the system) to all orders in the large proper time expansion from the retarded correlation function generalizing how quasinormal modes are obtained from the poles of the retarded propagator in thermal equilibrium. To all orders, the Schwinger-Keldysh correlation functions satisfy a matrix factorized form implying a bilocal thermal structure.

Let us briefly provide some explicit expressions. The d-dimensional Bjorken flow in the field theory is conveniently expressed in the Milne coordinates, $\tau = (t^2 - z^2)^{1/2}$ (the proper time with z being the longitudinal Minkowski coordinate), the rapidity $\zeta = \arctanh(z/t)$, and the transverse coordinates $\vec{x}_{\perp}$.

It is convenient to define

$$ \sigma = \tau^{\frac{d-2}{d-1}} \tau_0^{\frac{1}{d-1}}, \quad \hat{\zeta} = \zeta \tau_0, \qquad (7) $$

with $\tau_0$ being an arbitrarily chosen fixed proper time. Furthermore, for the arguments $(\sigma_1, \zeta_1, \vec{x}_{\perp 1})$ and $(\sigma_2, \zeta_2, \vec{x}_{\perp 2})$ of the two-point correlation function let us define

$$ \bar{\sigma} = \frac{1}{2}(\sigma_1 + \sigma_2), \quad \sigma_r = \sigma_1 - \sigma_2, \quad \tilde{\zeta}_r = (\hat{\zeta}_1 - \hat{\zeta}_2) \frac{\bar{\sigma}}{\tau_0}, \\ \tilde{x}_{\perp r} = |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}| \left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{1}{d-2}}. \qquad (8) $$

As discussed above, at late proper time, the Schwinger-Keldysh correlation functions $G(\sigma_1, \sigma_2, \hat{\zeta}_1 - \hat{\zeta}_2, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|)$ of the Bjorken flow can be mapped to thermal correlation functions $G_{\beta}$ at a specific temperature $T = \beta^{-1}$ as follows:

$$ G(\sigma_1, \sigma_2, \hat{\zeta}_1 - \hat{\zeta}_2, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|) \\ \rightarrow \left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{2\Delta_O}{d-2}} G_{\beta} \left( \frac{d-1}{d-2} \sigma_r, \sqrt{\tilde{\zeta}_r^2 + \tilde{x}_{\perp r}^2} \right). \qquad (9) $$

This limit implies that both $\sigma_1$ and $\sigma_2$ are large. Above $\Delta_O$ is the scaling dimension of the operator (which is determined by the mass of the dual bulk scalar field). Also both $G$ and $G_{\beta}$ are $2 \times 2$ matrices with indices determining whether the arguments $(\sigma_1, \zeta_1, \vec{x}_{\perp 1})$ and $(\sigma_2, \zeta_2, \vec{x}_{\perp 2})$ are in the forward or backward legs of the Schwinger-Keldysh time-contour. Finally, the temperature $T$ of the thermal is determined by the (only free) parameter $\epsilon_0$ of the Bjorken flow [see Eq. (5)] as follows:

$$ \beta = T^{-1} = \frac{4\pi\epsilon_0^{1/d}}{d}, \quad \epsilon_0 = \frac{16\pi G_N}{d-1}\epsilon_0, \qquad (10) $$

where the $(d+1)$-dimensional gravitational constant $G_N$ is given by the rank of the gauge group of the dual theory [such as, for instance, in $N=4$ super Yang-Mills theory with SU(N) gauge group in 4 dimensions, $G_N^{-1} = 2N^2/\pi$].<fnref n="5" /> Remarkably, the map to the thermal form involving space-time reparametrization implies the emergence of SO(d-1) rotational symmetry, and also time translation symmetry although these are absent in the original coordinates $\tau$ and $\zeta$.

As mentioned above, we can systematically include viscous and higher-order corrections to the correlation functions and obtain

$$ G(\sigma_1, \sigma_2, \hat{\zeta}_1 - \hat{\zeta}_2, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|) \\ = \left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{2\Delta_O}{d-2}} \sum_{n=0}^{\infty} \frac{1}{\bar{\sigma}^n \epsilon_0^{n/d}} G_n(T_0 \sigma_r, T_0 \tilde{\zeta}_r, T_0 \tilde{x}_{\perp r}), \qquad (11) $$

where $G_0$ coincides with the thermal correlation function $G_{\beta}$ given by (9).

We also note that our result that the horizon cap for the holographic hydrodynamic Bjorken flow should be pinned to the nonequilibrium event horizon captures the causal nature

<footnote n="5">⁵Apparently, the thermal form on the right-hand side of Eq. (10) depends on the choice of $\tau_0$. However, a shift in $\tau_0$ also changes the parameter $\epsilon_0$ appropriately so that the thermal form in Eq. (10) is unambiguous. This is explained in Sec. VC explicitly.</footnote>
of the Schwinger-Dyson equations for the real-time (out-of-equilibrium) correlation functions in the dual field theory.<fnref n="6" />

The series in Eq. (11) is not expected to be convergent, and therefore requires a trans-series completion with appropriate Stokes data (distinct from the Stokes data for the expectation value of the energy density which defines the bulk geometry) which should be actually functions of $\sigma_r$, $\tilde{\zeta}_r$, and $\tilde{x}_{\perp 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.

### D. Plan

The paper is organized as follows: In Sec. II, 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 different matrix factorization of thermal correlation functions. In Sec. III, 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 Sec. IV, 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 nonequilibrium event horizon. In Sec. V, 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 space-time reparametrizations, we show how we systematically obtain the corrections in a proper time expansion. We also discuss many nontrivial consistency checks of our results. In Sec. VI, we present a discussion on how a trans-series completion of this expansion can lead to seeing the quantum fluctuations behind the nonequilibrium event horizon, and matching with initial data lost during hydrodynamization.

Finally, we conclude in Sec. VII with an outlook.

## II. THE CROSSLEY-GLORIOISO-LIU 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 in Ref. [19]. In Ref. [19], it has also been verified that the retarded correlation function is implied by the ingoing boundary condition, as demanded by the Son-Starinets prescription [15] up to the quadratic order in frequency. These were sufficient to obtain a rudimentary effective theory of diffusion and dissipative hydrodynamics from holography [19,24–34]. Here, we present an 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.<fnref n="7" /> For other approaches, see Ref. [33].

The generating functional for the thermal Schwinger-Keldysh correlation functions in a quantum field theory is<fnref n="8" />

$$e^{W[J_1, J_2]} = \mathrm{Tr}[\hat{\rho}_\beta T_c e^{i \int dt \int d^{d-1}x [\hat{O}_1(\mathbf{x},t) J_1(\mathbf{x},t) - \hat{O}_2(\mathbf{x},t) J_2(\mathbf{x},t)]}], \quad (12)$$

with $\rho_\beta$ 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 $[\mathrm{Tr}(\hat{\rho}_\beta \cdot)]$

$$\begin{aligned} G_{11}(t-t', \mathbf{x}-\mathbf{x}') &= \frac{\delta^2 W}{\delta J_1(t, \mathbf{x}) \delta J_1(t', \mathbf{x}')} \\ &= -i \langle T(\hat{O}(t, \mathbf{x}) \hat{O}(t', \mathbf{x}')) \rangle, \\ -G_{12}(t-t', \mathbf{x}-\mathbf{x}') &= \frac{\delta^2 W}{\delta J_1(t, \mathbf{x}) \delta J_2(t', \mathbf{x}')} \\ &= i \langle \hat{O}(t', \mathbf{x}') \hat{O}(t, \mathbf{x}) \rangle, \\ -G_{21}(t-t', \mathbf{x}-\mathbf{x}') &= \frac{\delta^2 W}{\delta J_2(t, \mathbf{x}) \delta J_1(t', \mathbf{x}')} \\ &= i \langle \hat{O}(t, \mathbf{x}) \hat{O}(t', \mathbf{x}') \rangle, \\ G_{22}(t-t', \mathbf{x}-\mathbf{x}') &= \frac{\delta^2 W}{\delta J_2(t, \mathbf{x}) \delta J_2(t', \mathbf{x}')} \\ &= -i \langle \bar{T}(\hat{O}(t, \mathbf{x}) \hat{O}(t', \mathbf{x}')) \rangle. \end{aligned} \quad (13)$$

Succinctly we can write the above as

$$G_{ij}(x_1, x_2) = \frac{\delta^2 W}{\delta J_i(x_1) \delta J_j(x_2)} (-)^{i+j}, \quad (14)$$

with $(i, j) = (1, 2)$. The above holds even out of equilibrium with $\rho_\beta$ in (12) replaced by an arbitrary initial state $\rho_0$.

<footnote n="7">⁷The real-time prescription [17,18] of van Rees and Skenderis leads to the ingoing boundary condition as shown in Ref. [35] in thermal equilibrium. Our methods discussed here discuss a natural generalization away from equilibrium for out of equilibrium, especially hydrodynamic states.</footnote>

<footnote n="8">⁸Unless specified, we always put the backward arm of the Schwinger-Keldysh time contour infinitesimally below the real axis. Also, we often omit explicit mention of the appendage of the contour along the imaginary axis, which creates the thermal state in the infinite past.</footnote>

<footnote n="6">⁶The latter is manifest when written in terms of the coupled evolution of the statistical and spectral functions [1].</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. (left) The general Schwinger-Keldysh contour for an arbitrary initial density matrix $\hat{\rho}_0$. Field configurations need to be specified at the two ends of the closed time contour which are represented in bold. (right) In thermal equilibrium this simplifies. One needs an appendage to the contour extending along the negative imaginary axis by $\beta = T^{-1}$ and then impose periodic boundary conditions for the full contour.</figcaption></figure>

It can readily be shown that the KMS periodicity [arising from the Tr in (12) after extending the contour along the negative imaginary axis by $\beta$ 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 (backward) arms of the contour) defined as

$$G_{ij}(\omega, \mathbf{k}) = \int dt d^{d-1}x e^{i\omega(t-t')} e^{-i\mathbf{k}\cdot(\mathbf{x}-\mathbf{x}')} G_{ij}(t-t', \mathbf{x}-\mathbf{x}') \quad (15)$$

assume the form

$$G_{ij} = \begin{pmatrix} G_R(\omega, \mathbf{k})[1 + n(\omega)] - G_A(\omega, \mathbf{k})n(\omega) & -G_R(\omega, \mathbf{k}) - G_A(\omega, \mathbf{k})n(\omega) \\ -G_R(\omega, \mathbf{k}) - G_A(\omega, \mathbf{k})[1 + n(\omega)] & G_R(\omega, \mathbf{k})n(\omega) - G_A(\omega, \mathbf{k})[1 + n(\omega)] \end{pmatrix}, \quad (16)$$

where

$$G_R(t - t', \mathbf{x} - \mathbf{x}') = -i\theta(t - t') \langle [\hat{O}(t, \mathbf{x}), \hat{O}(t', \mathbf{x}')] \rangle \quad (17)$$

is the retarded propagator,

$$G_A(t - t', \mathbf{x} - \mathbf{x}') = -i\theta(t' - t) \langle [\hat{O}(t', \mathbf{x}'), \hat{O}(t, \mathbf{x})] \rangle \quad (18)$$

is the advanced propagator, and $n(\omega) = 1/(e^{\beta\omega} - 1)$ is the Bose-Einstein distribution function. It is easy to see from these definitions that

$$G_A^*(\omega, \mathbf{k}) = G_R(\omega, \mathbf{k}). \quad (19)$$

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). The Schwinger-Keldysh thermal correlation functions (16) obtained by differentiating the real-time partition function at a temperature $T = \beta^{-1}$ can be factorized as shown below:

$$G_{ij}(\omega, \mathbf{k}) = \sigma_3 \cdot \begin{pmatrix} A(\omega, \mathbf{k}) & B(\omega, \mathbf{k}) \\ A(\omega, \mathbf{k}) & B(\omega, \mathbf{k})e^{\beta\omega} \end{pmatrix} \cdot \begin{pmatrix} a(\omega, \mathbf{k}) & b(\omega, \mathbf{k}) \\ a(\omega, \mathbf{k}) & b(\omega, \mathbf{k})e^{\beta\omega} \end{pmatrix}^{-1}, \quad (20)$$

in which $\sigma_3 = \text{diag}(1, -1)$ is the third Pauli matrix, and

$$G_R(\omega, \mathbf{k}) = \frac{A(\omega, \mathbf{k})}{a(\omega, \mathbf{k})} = \frac{B^*(\omega, \mathbf{k})}{b^*(\omega, \mathbf{k})} = G_A^*(\omega, \mathbf{k}). \quad (21)$$

Clearly $A \to \lambda A$, $B \to \tilde{\lambda}B$, $a \to \lambda a$, and $b \to \tilde{\lambda}b$ gives the same thermal matrix, so the factorization is unique up to the multiplicative complex functions $\lambda(\omega, \mathbf{k})$ and $\tilde{\lambda}(\omega, \mathbf{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 ingoing Eddington-Finkelstein (EF) coordinates. The ingoing 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 $\mp\epsilon$ (i.e., $r_1 \to r_1 - i\epsilon$ and $r_2 \to r_2 + i\epsilon$). 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(\epsilon)$ 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

$$ds^2 = -\frac{2}{r^2} dr dt - \frac{1}{r^2} (1 - r^d \epsilon_0) dt^2 + \frac{1}{r^2} (dx_1^2 + \cdots + dx_{d-1}^2). \quad (22)$$

where $r$ is the bulk radial coordinate, $t$ is the Eddington-Finkelstein time, and the horizon is at $r = r_h = \epsilon_0^{-1/d}$. 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

$$\Phi(r, t, \mathbf{x}) = \int \frac{d\omega}{2\pi} \frac{d^{d-1}k}{(2\pi)^{d-1}} e^{-i\omega t} e^{i\mathbf{k}\cdot\mathbf{x}} \Phi(r, \omega, \mathbf{k}). \quad (23)$$

On-shell, $\Phi(r, \omega, \mathbf{k})$ is a sum of two linearly independent solutions $\phi_{\text{in}}(r, \omega, \mathbf{k})$ and $\phi_{\text{out}}(r, \omega, \mathbf{k})$ which are ingoing and
<figure><img src="image_4.png" /><figcaption>FIG. 2. (top) The radial contour in the complexified two-sheeted black brane geometry on a constant time hypersurface. The radial coordinate goes around the horizon $r = r_h$ in the complex plane forming the horizon cap, and connects the two arms of the Schwinger-Keldysh contour at the two boundaries at a constant Lorentzian time. The analytic continuation along the horizon cap gives a well-defined Dirichlet value problem for the bulk fields in terms of their sources specified at the two boundaries. (bottom) Shows the time evolution of the radial contour (shown in blue lines) in the full complexified space-time ending at the two Lorentzian arms at the two boundaries. The reversal of the direction of time in the second copy implies that the full complexified space-time has a single orientation after the two copies are smoothly glued with the horizon cap.</figcaption></figure>

outgoing at the horizon respectively. Therefore,

$$ \Phi(r, \omega, \mathbf{k}) = \phi_{\text{in}}(r, \omega, \mathbf{k})p(\omega, \mathbf{k}) + \phi_{\text{out}}(r, \omega, \mathbf{k})q(\omega, \mathbf{k}) \quad (24) $$

generally with $p(\omega, \mathbf{k})$ and $q(\omega, \mathbf{k})$ representing the arbitrary Fourier coefficients of the solutions which are ingoing and outgoing at the horizon, respectively. The latter thus provide a basis of solutions for given $\omega$ and $\mathbf{k}$, and can be uniquely defined via the following conditions:

$$ \phi_{\text{in}}(r_h, \omega, \mathbf{k}) = 1, \quad \lim_{r \to r_h} \phi_{\text{out}}(r, \omega, \mathbf{k})(r_h - r)^{-\frac{i\beta\omega}{2\pi}} = 1, \quad (25) $$

where

$$ \beta = T^{-1} = \frac{4\pi r_h}{d} \quad (26) $$

is the inverse Hawking temperature of the black brane, and $r_h = \epsilon_0^{-1/d}$ is the radial location of the horizon. Indeed, near the horizon ($r \approx r_h$),

$$ \Phi(r, \omega, \mathbf{k}) \approx p(\omega, \mathbf{k}) + (r_h - r)^{\frac{i\beta\omega}{2\pi}} q(\omega, \mathbf{k}), \quad (27) $$

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 space-time to another then implies that the Fourier coefficients of the on-shell solutions in the two copies are related by

$$ p_2(\omega, \mathbf{k}) = p_1(\omega, \mathbf{k}), \quad q_2(\omega, \mathbf{k}) = e^{\beta\omega} q_1(\omega, \mathbf{k}), \quad (28) $$

with 1 and 2 denoting the copies ending on the forward and backward arms of the time contour respectively at their boundaries. The on-shell solution in the full geometry can therefore be written in the following matrix form:

$$ \begin{pmatrix} \Phi_1(r, t, \mathbf{x}) \\ \Phi_2(r, t, \mathbf{x}) \end{pmatrix} = \int \frac{d\omega}{2\pi} \frac{d^{d-1}k}{(2\pi)^{d-1}} e^{-i\omega t} e^{i\mathbf{k}\cdot\mathbf{x}} \mathcal{M}(r, \omega, \mathbf{k}) \cdot \begin{pmatrix} p(\omega, \mathbf{k}) \\ q(\omega, \mathbf{k}) \end{pmatrix}, \quad (29) $$

with the matrix

$$ \mathcal{M}(r, \omega, \mathbf{k}) = \begin{pmatrix} \phi_{\text{in}}(r, \omega, \mathbf{k}) & \phi_{\text{out}}(r, \omega, \mathbf{k}) \\ \phi_{\text{in}}(r, \omega, \mathbf{k}) & e^{\beta\omega} \phi_{\text{out}}(r, \omega, \mathbf{k}) \end{pmatrix} \quad (30) $$

providing a basis of solutions for the entire complexified space-time comprising of the two copies smoothly glued at the horizon. The sources $J_1(\omega, \mathbf{k})$ and $J_2(\omega, \mathbf{k})$ specified at the two boundaries (see below) implement the Dirichlet boundary conditions that determine $p(\omega, \mathbf{k})$ and $q(\omega, \mathbf{k})$ uniquely for real frequencies and momenta, and thus yielding a unique bulk field configuration in the full complexified space-time.

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 $\hat{\phi}$, on the full complexified space-time, i.e.,

$$ W[J_1, J_2] = iS_{\text{on-shell}}. \quad (31) $$

Assuming minimal coupling to gravity, the on-shell quadratic action for the bulk scalar field $\Phi$ dual to a scalar operator takes the form

$$ S_{\text{on-shell}} = S_{\text{in-in}} + S_{\text{in-out}} + S_{\text{out-out}}. \quad (32) $$

The first piece $S_{\text{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_{\text{in-in}} = 0$. Note if we keep the ingoing 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_{\text{in-in}} = 0$ is consistent with field theory.

The second piece $S_{\text{in-out}}$, which is the sum of cross terms between the in and outgoing 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_{\text{in-out}} = S_{\text{bdy1}} + S_{\text{bdy2}}. \quad (33) $$

The third piece $S_{\text{out-out}}$ has a 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 nonanalytic piece $(r_h - r)^{\frac{i\beta\omega}{2\pi}}$, which we denote collectively as $S_\epsilon$. Essentially $S_\epsilon$
gets contributions from the following two terms:

$$S_{\epsilon} \propto \int d\omega \int d^d k \oint_{\epsilon} dr \sqrt{-G} [G^{rr} \partial_r \phi_{\text{out}}^* \partial_r \phi_{\text{out}} + G^{rt} (\partial_r \phi_{\text{out}}^* \partial_t \phi_{\text{out}} + \partial_r \phi_{\text{out}} \partial_t \phi_{\text{out}}^*)]. \quad (34)$$

Remarkably, the poles originating from these two terms cancel each other out [note that $\beta$ is given by (26)] resulting in $S_{\epsilon} = 0$. The remaining terms quadratic in the outgoing mode are analytic, so that $S_{\text{out-out}}$ is also the sum of two boundary terms. These two boundary terms cancel each other out as in the ingoing case. On the gravitational side, the easy way to see this is by first writing the contributions from the forward and backward arms of the contour separately. The boundary contributions from one arm in the integrand would be proportional to

$$\phi_b(-\omega, k)(\partial_r \phi)_b(\omega, k) + \dots,$$

with $\phi_b$ and $(\partial_r \phi)_b$ denoting the boundary values of $\phi$ and its radial derivative, respectively, and $\dots$ include counter terms, too. The contributions from the forward and backward parts of the contour come with opposite signs. If we consider the terms quadratic in the outgoing mode, then $\phi_b(\omega, k)$ picks up a factor of $e^{\beta\omega}$ while $\phi_b(-\omega, k)$ picks up a factor of $e^{-\beta\omega}$ via analytic continuation through the horizon cap, and the product of these factors is unity. Therefore, the contributions from the forward and backward contours cancel out leading to $S_{\text{out-out}} = 0$.

It is useful to see this also from the field-theory perspective. If we keep the outgoing modes only, then $J_2(\omega, \mathbf{k}) = J_1(\omega, \mathbf{k})e^{\beta\omega}$. In any field theory [29],

$$W[J_1(t, \mathbf{x}), J_2(t, \mathbf{x})] = W_T[J_2(t - i\beta, \mathbf{x}), J_1(t, \mathbf{x})], \\ \text{i.e., } W[J_1(t, \mathbf{x}), J_2(t + i\beta, \mathbf{x})] = W_T[J_2(t, \mathbf{x}), J_1(t, \mathbf{x})],$$

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.<fnref n="9" /> For $J_1(t, \mathbf{x}) = J_2(t, \mathbf{x}) = J(\omega, \mathbf{k})$, this amounts to

$$W_T[J(\omega, \mathbf{k}), J(\omega, \mathbf{k})] = W[J(\omega, \mathbf{k}), J(\omega, \mathbf{k})e^{\beta\omega}].$$

The left-hand side (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 right-hand side (RHS) of the above equation should vanish too, implying that

$$W[J_2(\omega, \mathbf{k}) = J_1(\omega, \mathbf{k})e^{\beta\omega}] = 0. \quad (35)$$

Thus, $S_{\text{out-out}} = 0$ is consistent with field theory. See also footnote 12 for a more straightforward verification that the thermal correlators in the dual theory originate from $S_{\text{in-out}}$ alone.

The upshot is that we obtain only two boundary contributions from the cross term between the ingoing and outgoing

modes, so that

$$S_{\text{on-shell}} = S_{\text{in-out}} = S_{\text{bdy1}} + S_{\text{bdy2}}, \quad (36)$$

where $S_{\text{bdy1}}$ and $S_{\text{bdy2}}$ are the contributions from the two boundaries after taking into account counter-terms necessary for holographic renormalization [36]. This implies that

$$S_{\text{on-shell}}[J_1, J_2] = \int dt \int d^{d-1}x [\langle O_1(t, \mathbf{x}) \rangle J_1(t, \mathbf{x}) - \langle O_2(t, \mathbf{x}) \rangle J_2(t, \mathbf{x})], \quad (37)$$

with

$$\begin{align}
\begin{pmatrix} J_1(\omega, \mathbf{k}) \\ J_2(\omega, \mathbf{k}) \end{pmatrix} &= S(\omega, \mathbf{k}) \cdot \begin{pmatrix} p(\omega, \mathbf{k}) \\ q(\omega, \mathbf{k}) \end{pmatrix}, \\
\begin{pmatrix} \langle O_1(\omega, \mathbf{k}) \rangle \\ \langle O_2(\omega, \mathbf{k}) \rangle \end{pmatrix} &= \mathcal{R}(\omega, \mathbf{k}) \cdot \begin{pmatrix} p(\omega, \mathbf{k}) \\ q(\omega, \mathbf{k}) \end{pmatrix} \\
&= (\mathcal{R} \cdot S^{-1})(\omega, \mathbf{k}) \cdot \begin{pmatrix} J_1(\omega, \mathbf{k}) \\ J_2(\omega, \mathbf{k}) \end{pmatrix}. \tag{38}
\end{align}$$

The matrices $\mathcal{R}$ and $S$ are defined as follows: Let the asymptotic ($r \approx 0$) expansions of the ingoing and outgoing modes be<fnref n="10" />

$$\begin{align}
\phi_{\text{in}}(r, \omega, \mathbf{k}) &= r^{d-\Delta}[a_0(\omega, \mathbf{k}) + \dots] + r^{\Delta}[A_0(\omega, \mathbf{k}) + \dots], \\
\phi_{\text{out}}(r, \omega, \mathbf{k}) &= r^{d-\Delta}[b_0(\omega, \mathbf{k}) + \dots] + r^{\Delta}[B_0(\omega, \mathbf{k}) + \dots]. \tag{39}
\end{align}$$

Then,

$$S(\omega, \mathbf{k}) = \begin{pmatrix} a_0(\omega, \mathbf{k}) & b_0(\omega, \mathbf{k}) \\ a_0(\omega, \mathbf{k}) & e^{\beta\omega}b_0(\omega, \mathbf{k}) \end{pmatrix} = \lim_{r \to 0} r^{\Delta-d} \mathcal{M}(r, \omega, \mathbf{k}), \quad (40)$$

and<fnref n="11" />

$$\mathcal{R}(\omega, \mathbf{k}) = (2\Delta - d) \begin{pmatrix} A_0(\omega, \mathbf{k}) & B_0(\omega, \mathbf{k}) \\ A_0(\omega, \mathbf{k}) & e^{\beta\omega}B_0(\omega, \mathbf{k}) \end{pmatrix} + \dots. \quad (41)$$

The $\dots$ above stands for (state-independent) contact terms which we ignore. Denoting

$$\hat{G}(\omega, \mathbf{k}) = (\sigma_3 \cdot \mathcal{R} \cdot S^{-1})(\omega, \mathbf{k}) \quad (42)$$

[with $\sigma_3 = \text{diag}(1, -1)$], we find from (37), (38), and (40)–(42) that<fnref n="12" />

$$\begin{align}
S_{\text{on-shell}}[J_1, J_2] &= \int \frac{d\omega}{2\pi} \frac{d^{d-1}k}{(2\pi)^{d-1}} [J_1(-\omega, -\mathbf{k}) \langle O_1(\omega, \mathbf{k}) \rangle \\
&\qquad -J_2(-\omega, -\mathbf{k}) \langle O_2(\omega, \mathbf{k}) \rangle] \\
&= \int \frac{d\omega}{2\pi} \frac{d^{d-1}k}{(2\pi)^{d-1}} [J_1(-\omega, -\mathbf{k}) J_2(-\omega, -\mathbf{k})] \\
&\qquad \cdot \hat{G}(\omega, \mathbf{k}) \cdot \begin{pmatrix} J_1(\omega, \mathbf{k}) \\ J_2(\omega, \mathbf{k}) \end{pmatrix}. \tag{43}
\end{align}$$

<footnote n="10">¹⁰The scaling dimension $\Delta$ is related to the mass via $\Delta = \frac{d}{2} + (\frac{d^2}{4} + m^2 L^2)^{1/2}$ with $L$ being the anti-de Sitter (AdS) radius.</footnote>

<footnote n="11">¹¹Note that, asymptotically, $\mathcal{M}(r, \omega, \mathbf{k}) = r^{d-\Delta}[S(\omega, \mathbf{k}) + \dots] + r^{\Delta}[\mathcal{R}(\omega, \mathbf{k}) + \dots]$.</footnote>

<footnote n="12">¹²$J_i(\omega, \mathbf{k}) = \int dt \, d^{d-1}x e^{i\omega t} e^{-i\mathbf{k} \cdot \mathbf{x}} J_i(t, \mathbf{x})$.</footnote>

<footnote n="9">⁹Succinctly, $W[J_1, J_2] = \text{Tr}(\rho_0 \tilde{U}[J_2]U[J_1])$, where $U[J_1]$ is forward evolution with source $J_1$ and $\tilde{U}[J_2]$ is backward evolution with source $J_2$. Similarly, $W_T[J_1, J_2] = \text{Tr}(\tilde{U}[J_2]\rho_0 U[J_1]) = \text{Tr}(\rho_0 U[J_1]\tilde{U}[J_2])$.</footnote>
Therefore, the identification (31) together with (13) implies that<fnref n="13" />

$$G_{ij}(\omega, \mathbf{k}) = \frac{\partial^2 S_{\text{on-shell}}}{\partial J_i(-\omega, -\mathbf{k}) \partial J_j(\omega, \mathbf{k})} (2\pi)^d = \hat{G}(\omega, \mathbf{k}). \quad (44)$$

From the matrix factorization of thermal correlation functions given by (20), we readily find from (40)–(42) that the correlation functions obtained by differentiating the on-shell gravitational action are thermal, i.e., assume the form (16), provided<fnref n="14" />

$$\begin{aligned} G_R(\omega, \mathbf{k}) &= (2\Delta - d) \frac{A_0(\omega, \mathbf{k})}{a_0(\omega, \mathbf{k})}, \\ G_A(\omega, \mathbf{k}) &= (2\Delta - d) \frac{B_0(\omega, \mathbf{k})}{b_0(\omega, \mathbf{k})}. \end{aligned} \quad (45)$$

Remarkably, the above are exactly the Son-Starinets prescriptions [15] for the retarded and advanced propagators according to which they are obtained from the ingoing and outgoing boundary conditions at the horizon, respectively. Furthermore, since the outgoing mode is time reverse of the ingoing mode (which is not manifest in the Eddington-Finkelstein gauge but can be evident from transforming to Schwarzschild-like coordinates),<fnref n="15" /> we should have

$$\frac{B_0^*(\omega, \mathbf{k})}{b_0^*(\omega, \mathbf{k})} = \frac{A_0(\omega, \mathbf{k})}{a_0(\omega, \mathbf{k})}, \quad \text{i.e., } G_A^*(\omega, \mathbf{k}) = G_R(\omega, \mathbf{k}), \quad (46)$$

holds.

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 [37]. However, in this case, the backward part of the time contour needs to be shifted by $\beta/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 nonequilibrium 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

<figure><img src="image_5.png" /><figcaption>FIG. 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 are specified on a constant $\tau$ hyperboloid.</figcaption></figure>

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 space-time with the horizon cap. However, this is outside the scope of the present work, and therefore we do not further discuss this issue. Finally, we note that the arguments presented here are simpler compared with Ref. [19] since we do not employ an expansion about $\omega = 0$ which obscures the analytic continuation at the horizon cap by producing $(\ln \omega)^n$ terms.

### III. BJORKEN FLOW, WEYL RESCALING, AND THE HOLOGRAPHIC DUAL

#### A. Bjorken flow and its Weyl rescaling

Bjorken flow [38] is a simple model describing the expanding plasma produced by heavy ion collisions. This model is based on the assumptions of boost invariance, and translational and rotational symmetries in the transverse directions of an expanding system. The evolution occurs inside a forward light cone, as shown in Fig. 3. It is convenient to describe the Bjorken flow in the Milne proper time $\tau$ and the rapidity $\zeta$ which are related to the Minkowski (lab frame) time $t$ and the longitudinal coordinate $z$ (along which the system is expanding) as

$$t = \tau \cosh \zeta \quad \text{and} \quad z = \tau \sinh \zeta.$$

The transverse coordinates $x_\perp$ are the same in both Milne and Minkowski coordinate systems. In the Milne coordinates, the Minkowski metric takes the form

$$ds^2 = -d\tau^2 + \tau^2 d\zeta^2 + ds_\perp^2, \quad (47)$$

where $ds_\perp$ 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 $\tau$. 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

<footnote n="13">¹³The reader can check that substituting $p(\omega, \mathbf{k}) + q(\omega, \mathbf{k}) = J_1(\omega, \mathbf{k})$ and $p(\omega, \mathbf{k}) + q(\omega, \mathbf{k})e^{i\beta\omega} = J_2(\omega, \mathbf{k})$ in the on-shell action (43), and using the thermal form of the propagators below, that indeed only cross terms between $p$ and $q$, i.e., the in and outgoing modes appear. There are no contributions from terms quadratic in $p$ or in $q$, implying that $S_{\text{in-in}} = S_{\text{out-out}} = 0$ as claimed above, and also $S_{\text{on-shell}} = S_{\text{in-out}}$.</footnote>

<footnote n="14">¹⁴It is obvious that, to map to the factorization in (20), we have to set $A = (2\Delta - d)A_0$, $B = (2\Delta - d)B_0$, $a = a_0$, and $b = b_0$.</footnote>

<footnote n="15">¹⁵Note that the notion of in- and outgoing modes are gauge invariant up to overall multiplicative factors, but these cancel out in the ratio of the normalizable to the non-normalizable modes. This is why the Son-Starinets prescription is also gauge invariant.</footnote>
symmetries of the Bjorken flow, can take the form

$$ \langle T_{\mu\nu} \rangle = \mathrm{diag} \left( \epsilon(\tau), \tau^2 p_L(\tau), \underbrace{p_T(\tau), \dots, p_T(\tau)}_{(d-2) \text{ times}} \right) \quad (48) $$

in the Milne coordinates for $d > 2$. Clearly, $\epsilon, p_L$, and $p_T$ denote the energy density, longitudinal and transverse pressures, respectively. The local conservation of energy and momentum $\nabla_\mu T^{\mu\nu} = 0$ implies that

$$ p_L(\tau) = -\epsilon(\tau) - \tau \partial_{\tau} \epsilon(\tau), \quad (49) $$

and the conformal Ward identity $T^\mu_\mu = 0$ imposes

$$ p_T(\tau) = \frac{2}{d-2} \left[ \epsilon(\tau) + \frac{1}{2} \tau \partial_{\tau} \epsilon(\tau) \right]. \quad (50) $$

Consequently, the evolution of the energy-momentum tensor is determined by $\epsilon(\tau)$ in a conformal field theory.

At large proper time $\tau$, the Bjorken flow admits a hydrodynamic description [39]. To explicitly map the energy-momentum tensor (48) to that of a fluid, we need to set the flow velocity as

$$ u^\mu = (1, \vec{0}), $$

i.e., $d\tau$ is comoving with the flow in the Milne coordinates. In a conformal system, the large proper time expansion of $\epsilon(\tau)$ is given by a single parameter, namely

$$ \mu := \epsilon_0 \tau_0^{\frac{d}{d-1}}, \quad (51) $$

which is determined by the initial conditions —$\epsilon_0$ is a constant energy density, and $\tau_0$ can be chosen to be the value of $\tau$ where we initialize. The large proper time expansion of $\epsilon(\tau)$ takes the form

$$ T_{\tau}^{\tau} = -\epsilon(\tau) = -\mu \tau^{-\frac{d}{d-1}} \left( 1 + \sum_{n=1}^{\infty} \lambda_n \mu^{-\frac{n}{d}} \tau^{-\frac{n(d-2)}{d-1}} \right), \quad (52) $$

where $\lambda_n$ are (state-independent) constants that are determined by the transport coefficients of the microscopic theory. For instance, $\lambda_1$ is related to the shear viscosity $\eta$ as

$$ \lambda_1 = - \frac{\eta(\epsilon)}{\epsilon^{\frac{d-1}{d}}}, \quad (53) $$

which should indeed be a constant in a conformal theory. The leading term of the expansion $\propto \tau^{-d/(d-1)}$ 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 explicitly mention it later. Under a Weyl transformation which transforms the metric and the energy-momentum tensor as

$$ ds^2 \rightarrow \tilde{ds}^2 = \Omega(x)^2 ds^2, \quad T_{\mu\nu} \rightarrow \tilde{T}_{\mu\nu} = \Omega(x)^{-d+2} T_{\mu\nu}, \quad (54) $$

the new solutions of the hydrodynamic equations are given by

$$ u^{\mu} \rightarrow \tilde{u}^{\mu} = \Omega(x)^{-1} u^{\mu}, \quad \epsilon \rightarrow \tilde{\epsilon} = \Omega(x)^{-d} \epsilon, \quad (55) $$

in any conformal theory. Consider the combined operation of the time reparametrization

$$ \sigma = \tau^{\frac{d-2}{d-1}} \tau_0^{\frac{1}{d-1}} \quad (56) $$

and the Weyl scaling with

$$ \Omega(\sigma) = \left(\frac{\tau_0}{\sigma}\right)^{\frac{1}{d-2}}, \quad (57) $$

under which the Milne metric (47) transforms to (with $\hat{\xi} = \zeta \tau_0$)

$$ \tilde{ds}^2 = - \frac{(d-1)^2}{(d-2)^2} d\sigma^2 + \frac{\sigma^2}{\tau_0^2} d\hat{\xi}^2 + \left(\frac{\sigma}{\tau_0}\right)^{-\frac{2}{d-2}} ds_\perp^2, \quad (58) $$

and the energy-momentum tensor given by (48)–(50) transforms to

$$
\begin{align}
\tilde{T}_{\sigma\sigma} &= \Omega^{-d+2} \tau'(\sigma)^2 \epsilon(\tau(\sigma)) = \frac{(d-1)^2}{(d-2)^2} \left(\frac{\sigma}{\tau_0}\right)^{\frac{d}{d-2}} \epsilon(\tau(\sigma)) \nonumber \\
&= \frac{(d-1)^2}{(d-2)^2} \tilde{\epsilon}(\sigma) \quad \left[ \text{note } \tilde{\epsilon}(\sigma) := \left(\frac{\sigma}{\tau_0}\right)^{\frac{d}{d-2}} \epsilon(\tau(\sigma)) \right], \nonumber \\
\tilde{T}_{\hat{\xi}\hat{\xi}} &= \Omega^{-d+2} \left(\frac{\tau(\sigma)}{\tau_0}\right)^2 p_L(\tau(\sigma)) \nonumber \\
&= -\left(\frac{\sigma}{\tau_0}\right)^{\frac{3d-4}{d-2}} [\epsilon(\tau(\sigma)) + \tau(\sigma)\epsilon'(\tau(\sigma))] \nonumber \\
&= \left(\frac{\sigma}{\tau_0}\right)^2 \frac{1}{d-1} [\tilde{\epsilon}(\sigma) - (d-2)\sigma\tilde{\epsilon}'(\sigma)], \nonumber \\
\tilde{T}_{ii} &= \Omega^{-d+2} P_T(\tau(\sigma)) \nonumber \\
&= \frac{2}{d-2} \frac{\sigma}{\tau_0} \left[ \epsilon(\tau(\sigma)) + \frac{1}{2}\tau(\sigma)\epsilon'(\tau(\sigma)) \right] \nonumber \\
&= \left(\frac{\sigma}{\tau_0}\right)^{-\frac{2}{d-2}} \frac{1}{d-1} [\tilde{\epsilon}(\sigma) + \sigma\tilde{\epsilon}'(\sigma)], \quad (59)
\end{align}
$$

with $\tilde{T}_{ii}$ denoting the diagonal transverse components and $'$ denoting the derivative with respect to the argument of the corresponding function. It follows that in the hydrodynamic limit that the Bjorken expansion (52) takes the resultant form

$$
\begin{align}
\tilde{T}_{\sigma}^{\sigma} &\equiv -\tilde{\epsilon}(\sigma) = -\left(\frac{\sigma}{\tau_0}\right)^{\frac{d}{d-2}} \epsilon(\tau(\sigma)) \nonumber \\
&= -\epsilon_0 \left( 1 + \sum_{n=1}^{\infty} \lambda_n \epsilon_0^{-\frac{n}{d}} \sigma^{-n} \right). \tag{60}
\end{align}
$$

The Weyl scaled metric (58) has the property that

$$ \sqrt{-\tilde{g}} = \frac{d-1}{d-2} \quad (61) $$

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 (60), we have

$$ \lim_{\sigma \to \infty} \tilde{\epsilon}(\sigma) = \epsilon_0. \quad (62) $$
Therefore, instead of a perfect fluid expansion, the flow attains a constant temperature, energy and entropy densities at late time although no timelike Killing vector exists in the background metric. The latter feature leads to viscous and higher-order corrections. The large (reparametrized) proper time expansion is determined by $\epsilon_0$, the final thermal value of the energy density, while $\tau_0$ appears in the Weyl scaling factor $\Omega$ as explicit in Eq. (57).

The Weyl scaling depends explicitly on $\tau_0$. However, note that, for $\tau_0 \to \xi\tau_0$, we obtain from (56) that $\sigma \to \xi^{1/(d-1)}\sigma$. Thus the Weyl factor given by (57) scales as $\Omega \to \xi^{1/(d-1)}\Omega$, implying that $\tilde{d}s^2 \to \xi^{2/(d-1)}\tilde{d}s^2$ and $\epsilon_0 \to \xi^{-1/(d-1)}\epsilon_0$. Therefore, the dimensionless variables $\epsilon_0^{-1/d}\sigma^{-1}$ (which provides the proper time expansion parameter) and $\sigma^n\tilde{T}_\sigma^\sigma$ are invariant under $\tau_0 \to \xi\tau_0$ and are thus independent of $\tau_0$.

We also note that the Weyl transformation studied here is itself not corrected at first and higher orders in the large and proper time expansion.

## B. Gravitational setup

The gravitational dual of the Bjorken flow [22,40–44] has been extensively studied in the literature with the late time evolution providing a primary example of the fluid/gravity correspondence [5–8] where large order resummation of the hydrodynamic series [22] has been explicitly carried out revealing the hydrodynamization [45] of a far-from-equilibrium state. When a state hydrodynamizes, the energy-momentum tensor can be described as an optimally truncated (divergent and asymptotic) hydrodynamic series even when it is far from equilibrium [22,45]. In the context of the Bjorken flow, the evolution of the energy density approaches a hydrodynamic attractor [22,45]. This is a generic property of a many-body relativistic system irrespective of whether its degrees of freedom interact weakly or strongly (see Ref. [46] for a recent review).

Here we 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 section. 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 no timelike 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 section, 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 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 nonequilibrium event horizon for regularity, and therefore this gauge freedom will play an important role. We 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:

$$R_{MN} - \frac{1}{2}RG_{MN} - \frac{d(d-1)}{2L^2}G_{MN} = 0. \quad (63)$$

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 (63), lives at the boundary ($\rho = 0$) in the boundary metric, which is defined as

$$g_{\mu\nu}^{b} = \lim_{\rho \to 0} \rho^2 G_{\mu\nu}, \quad (64)$$

where $a$ and $b$ stand for the field-theory indices. Since we are considering the Bjorken flow in the Milne metric (47), $g_{\mu\nu}^b$ should coincide with it. Similarly, if we consider the Weyl scaled version of the Bjorken flow, the boundary metric should coincide with (58).

Before considering the Bjorken flow, it is useful to first understand the vacuum solution, which is pure (maximally symmetric) $AdS_{d+1}$ space-time with the desired boundary metric. In the ingoing Eddington-Finkelstein gauge, the vacuum state in the Milne metric (47) is thus dual to

$$ds^2 = -\frac{2}{r^2}drd\tau - \frac{1}{r^2}d\tau^2 + \left(1 + \frac{\tau}{r}\right)^2 d\xi^2 + \frac{1}{r^2}ds_{\perp}^2, \quad (65)$$

where $r$ is the radial coordinate. Similarly, the vacuum in Weyl scaled metric (58) is dual to

$$ds^2 = -\frac{2}{v^2} \frac{d-1}{d-2} dv d\sigma - \frac{1}{v^2} \left( \frac{(d-1)^2}{(d-2)^2} + \frac{2(d-1)v}{(d-2)^2 \sigma} \right) d\sigma^2 + \frac{1}{\tau_0^2} \left(1 + \frac{\sigma}{v}\right)^2 d\hat{\xi}^2 + \frac{1}{v^2} \left(\frac{\sigma}{\tau_0}\right)^{-\frac{2}{d-2}} ds_{\perp}^2, \quad (66)$$

where $v$ is the radial coordinate. These bulk metrics (65) and (66) are related by the diffeomorphism

$$\tau = \tau_0^{-\frac{1}{d-2}} \sigma^{\frac{d-1}{d-2}}, \quad r = v \left(\frac{\sigma}{\tau_0}\right)^{\frac{1}{d-2}}. \quad (67)$$

For both cases, (65) and (66), we obtain the boundary metrics (47) and (58) from (64), after replacing $\rho$ with $r$ and $v$, respectively. Any Weyl transformation at the boundary is dual to a bulk diffeomorphism. Since the boundary metrics (47) and (58) are related by a Weyl transformation, (67) is simply a specific instance of this general feature of holographic duality. Note that $\tau$ and $\sigma$ are related exactly by the time reparametrization (56) at the boundary.<fnref n="16" />

<footnote n="16">¹⁶Diffeomorphisms such as (67) which implement global transformations on the dual state are called improper diffeomorphisms which are always part of residual gauge freedom after gauge fixing in the bulk. The latter can also have additional proper diffeomorphisms which do not affect the dual physical quantities</footnote>
Holographic renormalization [20,21,36,47] provides the framework for extracting the $\langle T_{\mu\nu} \rangle$ corresponding to the state in the field theory dual to a specific asymptotically AdS$_{d+1}$ bulk geometry. The procedure essentially amounts to co-variantly regularizing the Brown-York tensor on a cutoff hypersurface with local counterterms built out of the induced metric, and then taking this surface to the boundary. For the bulk geometry (65), $\langle T_{\mu\nu} \rangle = 0$ in the dual vacuum state living in the flat Milne metric (47) at the boundary. On the other hand, for the vacuum state living on the Weyl transformed Milne metric (58) which is dual to the bulk geometry (66), $\langle T_{\mu\nu} \rangle = 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)

$$ \langle \tilde{T}_{\mu\nu} \rangle = \frac{1}{8\pi G_N} \mathcal{A}_{\mu\nu}, \qquad (68) $$

with

$$ \begin{aligned} \mathcal{A}_{\mu\nu} &= \frac{1}{16} \left[ \frac{4}{3} \tilde{R}_{\mu\nu} \tilde{R} - 2 \tilde{R}_{\mu\rho} \tilde{R}^{\rho}_{\nu} - \tilde{g}_{\mu\nu} \left( \frac{1}{2} \tilde{R}^2 - \tilde{R}_{\rho\sigma} \tilde{R}^{\rho\sigma} \right) \right] \\ &= \operatorname{diag} \left( -\frac{1}{32\sigma^4}, -\frac{11}{648\sigma^2\tau_0^2}, \frac{25\tau_0}{648\sigma^5}, \frac{25\tau_0}{648\sigma^5} \right), \end{aligned} \qquad (69) $$

where $\tilde{g}$ denotes the Weyl rescaled background metric (58), $\tilde{R}_{\mu\nu}$ is the Ricci tensor built out of it, etc. It is easy to verify that

$$ \tilde{\nabla}_{\mu} \tilde{T}^{\mu}_{\nu} = 0, \qquad (70) $$

i.e., energy and momentum are conserved in the Weyl rescaled background metric (58) (with $\tilde{\nabla}$ being the covariant derivative built out of it), and

$$ \tilde{T}^{\mu}_{\mu} = -\frac{1}{8\pi G_N} \left( \frac{1}{24} \tilde{R}^2 - \frac{1}{8} \tilde{R}_{\rho\sigma} \tilde{R}^{\rho\sigma} \right) = \frac{1}{8\pi G_N} \frac{2}{27\sigma^4}. \qquad (71) $$

Using

$$ \frac{1}{G_N} = \frac{2N^2}{\pi}, \qquad (72) $$

we can readily find that (71) reproduces the Weyl anomaly of SU(N) N = 4 supersymmetric Yang-Mills theory [20,47].

An asymptotically AdS$_{d+1}$ metric dual to a Bjorken flow on the flat Milne metric (47) at the boundary, is a solution to the vacuum Einstein's equations (63) which takes the form

$$ \begin{aligned} ds^2 = & -\frac{2}{r^2} dr d\tau - \frac{A(r, \tau)}{r^2} d\tau^2 + \left(1 + \frac{\tau}{r}\right)^2 e^{L(r, \tau)} d\zeta^2 \\ & + \frac{e^{K(r, \tau)}}{r^2} ds_{\perp}^2, \end{aligned} \qquad (73) $$

with the following Dirichlet asymptotic boundary conditions:

$$ A(r, \tau) \to 1, \quad K(r, \tau) \to 0, \quad \text{and} \quad L(r, \tau) \to 0 \text{ as } r \to 0. \qquad (74) $$

These boundary conditions ensure that the boundary metric (64) (with r being the radial coordinate in place of the generic $\rho$) coincides with the Milne metric (47).

The Einstein equations (63) can be readily solved in the late-time expansion, as functions of the scaling

variable

$$ s = r \left( \frac{\tau_0}{\tau} \right)^{\frac{1}{d-1}}, \qquad (75) $$

and with the expansion parameter being

$$ (\tilde{\mu}^{1/d} \tau^{\frac{d-2}{d-1}})^{-1}, \qquad (76) $$

where $\tilde{\mu} := \epsilon_0 \tau_0^{d/(d-1)}$ is a constant which will be related to the single parameter $\mu$ of the Bjorken flow defined in (51) (or equivalently to $\epsilon_0$) below. Explicitly,

$$ \begin{aligned} A(r, \tau) &= 1 - \epsilon_0 s^d + \sum_{i=1}^{\infty} (\tilde{\mu}^{1/d} \tau^{\frac{d-2}{d-1}})^{-i} a_{(i)}(\epsilon_0^{\frac{1}{d}} s), \\ L(r, \tau) &= \sum_{i=1}^{\infty} (\tilde{\mu}^{1/d} \tau^{\frac{d-2}{d-1}})^{-i} l_{(i)}(\epsilon_0^{\frac{1}{d}} s), \\ K(r, \tau) &= \sum_{i=1}^{\infty} (\tilde{\mu}^{1/d} \tau^{\frac{d-2}{d-1}})^{-i} k_{(i)}(\epsilon_0^{\frac{1}{d}} s). \end{aligned} \qquad (77) $$

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 $s_h = \epsilon_0^{-1/d}$, i.e., at

$$ r_h = \epsilon_0^{-\frac{1}{d}} \left( \frac{\tau}{\tau_0} \right)^{\frac{1}{d-1}}. \qquad (78) $$

Together with the Dirichlet boundary conditions (74), 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 [22]. 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 (78) 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-Finkelstein coordinates. We therefore keep this gauge freedom unfixed here and later show how it is fixed by the 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 (47)] and also the $\langle T_{\mu\nu} \rangle$ of the dual Bjorken flow extracted from holographic renormalization is invariant. It is useful to see this explicitly. For illustration, let us consider the AdS$_5$ case. The asymptotic
expansions take the form:

$$
\begin{align}
A(r, \tau) &= 1 + r a_1(\tau) + r^2 \left( \frac{a_1(\tau)^2}{4} - a'_1(\tau) \right) + r^4 a_4(\tau) \nonumber \\
&\quad + \dots, \nonumber \\
K(r, \tau) &= r a_1(\tau) - r^2 \frac{a_1(\tau)^2}{4} + r^3 \frac{a_1(\tau)^3}{12} + r^4 k_4(\tau) + \dots, \nonumber \\
L(r, \tau) &= r a_1(\tau) - r^2 \left( \frac{a_1(\tau)^2}{4} + \frac{a_1(\tau)}{\tau} \right) \nonumber \\
&\quad + r^3 \left( \frac{a_1(\tau)^3}{12} + \frac{a_1(\tau)^2}{2\tau} + \frac{a_1(\tau)}{\tau^2} \right) + r^4 l_4(\tau) \nonumber \\
&\quad + \dots. \tag{79}
\end{align}
$$

Above, the function $a_1(\tau)$ is related to the residual gauge freedom and can be chosen arbitrarily. Furthermore, the constraints of Einstein's equations (63) impose

$$
\begin{align}
a_4'(\tau) &= -\frac{a_1(\tau)^4}{12\tau} - \frac{4}{3\tau}a_4(\tau) - \frac{8}{3\tau}k_4(\tau), \nonumber \\
l_4(\tau) + 2k_4(\tau) &= -\frac{a_1(\tau)}{\tau^3} - \frac{3}{4}\frac{a_1(\tau)^2}{\tau^2} - \frac{1}{4}\frac{a_1(\tau)^3}{\tau} \nonumber \\
&\quad - \frac{3}{32}a_1(\tau)^4. \tag{80}
\end{align}
$$

Using these constraints, one can find via holographic renormalization that

$$
\langle T_{\mu\nu} \rangle = \frac{3}{16\pi G_N} \begin{pmatrix} -a_4(\tau) & 0 & 0 & 0 \\ 0 & \tau^2(a_4(\tau) + \tau a'_4(\tau)) & 0 & 0 \\ 0 & 0 & -a_4(\tau) - \frac{1}{2}\tau a'_4(\tau) & 0 \\ 0 & 0 & 0 & -a_4(\tau) - \frac{1}{2}\tau a'_4(\tau) \end{pmatrix}. \quad (81)
$$

First, the above result is exact to all orders in the late proper time expansion. Second, we readily find that $\langle T_{\mu\nu}(\tau) \rangle$ is determined by $a_4(\tau)$ alone and is independent of the arbitrary function $a_1(\tau)$ capturing the residual gauge freedom in the asymptotic expansion after utilizing the constraints (80) in the renormalized Brown-York stress tensor. Thus $\langle T_{\mu\nu}(\tau) \rangle$ is invariant under the residual gauge transformation. Furthermore, comparing (81) with (48)–(50) (for $d=4$) we find that $\langle T_{\mu\nu}(\tau) \rangle$ takes the general form of the energy-momentum tensor of Bjorken flow with the identification

$$ \epsilon(\tau) = - \frac{3}{16\pi G_N} a_4(\tau). \quad (82) $$

One can repeat the same exercise in arbitrary dimensions ($d > 2$) and show that $\langle T_{\mu\nu}(\tau) \rangle$ obtained from holographic renormalization takes the general form given by (48)–(50) with the identification

$$ \epsilon(\tau) = - \frac{d-1}{16\pi G_N} a_d(\tau), \quad (83) $$

where $a_d(\tau)$ is the coefficient of $r^d$ in the asymptotic expansions of $A(r, \tau)$.

Furthermore, extracting $a_d(\tau)$ from (77) we obtain that at the leading order in the late proper time expansion,

$$
\begin{align}
\epsilon(\tau) &\approx \frac{d-1}{16\pi G_N} \epsilon_0 \left(\frac{\tau_0}{\tau}\right)^{\frac{d}{d-1}} = \frac{d-1}{16\pi G_N} r_h(\tau)^{-d} \nonumber \\
&= \frac{d-1}{16\pi G_N} \left(\frac{4\pi T(\tau)}{d}\right)^d, \tag{84}
\end{align}
$$

where we have used (78), and also (26) to define an instantaneous Hawking temperature $T(\tau)$ given by

$$ T(\tau) = \frac{d}{4\pi} \epsilon_0^{\frac{1}{d}} \left( \frac{\tau_0}{\tau} \right)^{\frac{1}{d-1}}. \quad (85) $$

Once again comparing with the general (hydrodynamic) late proper time expansion (52), we find that

$$ \epsilon_0 = \frac{d-1}{16\pi G_N} \epsilon_0. \quad (86) $$

For the case of AdS$_5$, the identification (72) implies that at late proper time

$$
\begin{align}
T(\tau) &\approx \frac{1}{\pi} \epsilon_0^{1/4} \left(\frac{\tau_0}{\tau}\right)^{\frac{1}{4}}, \nonumber \\
\epsilon(\tau) &\approx \frac{3}{8} \pi^2 N^2 T(\tau)^4 = \frac{3}{8\pi^2} N^2 \epsilon_0 \left(\frac{\tau_0}{\tau}\right)^{\frac{4}{3}}. \tag{87}
\end{align}
$$

For any $d$, $\epsilon(\tau)$ is given by a perfect fluid expansion given by (84) at late time with

$$ p_L(\tau) \approx p_T(\tau) \approx \frac{1}{d-1} \epsilon(\tau) \approx \frac{1}{16\pi G_N} r_h(\tau)^{-d}. \quad (88) $$

Thus, at late proper time, the energy density $\epsilon(\tau)$ and the pressures $p_L(\tau)$ and $p_T(\tau)$ are given by the thermal equation of state obtained from a static black brane geometry, but with a time-dependent temperature (85) which satisfies the Euler equations.

To construct a regular horizon cap, it is useful to change coordinates from $r$ and $\tau$ to $v$ and $\sigma$ following (67), as in the case of the space-time dual to the vacuum state. Note, it follows from (75) that $v$ and $s$ are the same. In these new coordinates, the metric (73) takes the form

$$
\begin{align}
ds^2 = & -\frac{2}{v^2} \frac{d-1}{d-2} dv d\sigma \nonumber \\
& -\frac{1}{v^2} \left( \frac{(d-1)^2}{(d-2)^2} A(v, \sigma) + \frac{2(d-1)v}{(d-2)^2 \sigma} \right) d\sigma^2 \nonumber \\
& + \frac{1}{\tau_0^2} \left( 1 + \frac{\sigma}{v} \right)^2 e^{L(v, \sigma)} d\hat{\zeta}^2 + \frac{1}{v^2} \left( \frac{\sigma}{\tau_0} \right)^{-\frac{2}{d-2}} e^{K(v, \sigma)} ds_{\perp}^2, \tag{89}
\end{align}
$$
and is dual to the Bjorken flow on the Weyl rescaled background metric (58) at the boundary. We want to emphasize that, just like the Weyl transformation, the dual bulk diffeomorphism which achieves the above form of the metric is given by (67) exactly and is therefore not corrected at first and higher orders in the large proper time expansion.

Holographic renormalization and the constraints of the Einstein's equations (63) imply that the energy-momentum tensor of the dual Bjorken flow takes the form

$$ \langle \tilde{T}_{\mu\nu} \rangle = T_{\mu\nu}^{W} + \frac{1}{8\pi G_N} \mathcal{A}_{\mu\nu}, \qquad (90) $$

where $T_{\mu\nu}^W$ takes the general Bjorken form with nonvanishing components given by (59) in which

$$ \tilde{\epsilon}(\sigma) = -\frac{d-1}{16\pi G_N} \frac{(d-1)^2}{(d-2)^2} a_d(\sigma), \qquad (91) $$

and $\mathcal{A}_{\mu\nu}$ is the Weyl anomaly appearing for even $d$. Comparing with (83), we indeed verify that the bulk diffeomorphism (67) implements the Weyl transformation and time reparametrization in the dual theory [with the Weyl factor given by (57)] and also reproduces its Weyl anomaly. Particularly, for $d=4$, we recall that $\mathcal{A}_{\mu\nu}$ is simply given by (69). The anomalous term is state-independent (and is always the same as in the Weyl transformed vacuum state). We again note that $\langle \tilde{T}_{\mu\nu} \rangle$ is invariant under the residual gauge symmetry since it is independent of $a_1(\tau)$ after we implement the gravitational constraints following the previous discussion.

Obviously, the late proper time expansion (77) takes the form

$$ \begin{aligned} A(v, \sigma) &= 1 - \epsilon_0 v^d + \sum_{i=1}^{\infty} (\epsilon_0^{1/d} \sigma)^{-i} a_{(i)}(\epsilon_0^{1/d} v), \\ L(v, \sigma) &= \sum_{i=1}^{\infty} (\epsilon_0^{1/d} \sigma)^{-i} l_{(i)}(\epsilon_0^{1/d} v), \\ K(v, \sigma) &= \sum_{i=1}^{\infty} (\epsilon_0^{1/d} \sigma)^{-i} k_{(i)}(\epsilon_0^{1/d} v). \end{aligned} \qquad (92) $$

Explicitly, for $d=4$,

$$ \begin{aligned} a_{(1)}(x) &= \alpha_1 \frac{x(1+x^4)}{3} + \frac{2x^4(1+x)}{3}, \\ k_{(1)}(x) &= \alpha_1 \frac{x}{3} + \frac{1}{2}g(x), \\ l_{(1)}(x) &= \alpha_1 \frac{x}{3} - g(x), \end{aligned} \qquad (93) $$

where

$$ g(x) = \frac{4}{3}x - \frac{1}{3}\ln(x^2 + 1) - \frac{2}{3}\ln(x + 1) - \frac{2}{3}\arctan x. \qquad (94) $$

Above, $\alpha_1$ is the dimensionless parameter associated with the residual gauge freedom. At any order in the late proper time expansion, the terms multiplying $\alpha_1$ in (93) remain the same, however, we should replace $\alpha_1$ by $\alpha_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

$$ v_h = \epsilon_0^{-1/d} + O(\sigma^{-1}). \qquad (95) $$

At late proper time the dual black brane has a constant surface gravity and area although the directions longitudinal to the flow keep expanding and those transverse to the flow keep contracting.

The above metric reproduces the late proper time expansion of $\tilde{\epsilon}(\sigma)$ which takes the form (60) with $\epsilon_0$ given by (83) and $\lambda_n$ taking specific values for a given $d$. Particularly, for any $d$, we obtain

$$ \lambda_1 = - \frac{1}{(d-1)^{\frac{d-1}{d}} (16\pi G_N)^{1/d}}. \qquad (96) $$

It is easy to verify from (53) and the equation of state [see (84)]

$$ T = \beta^{-1} = \frac{d}{4\pi} \epsilon_0^{1/d}, \quad \epsilon = \frac{d-1}{16\pi G_N} \epsilon_0 = \frac{d-1}{16\pi G_N} \left( \frac{4\pi T}{d} \right)^d \qquad (97) $$

that (96) implies via (53) that

$$ \frac{\eta}{s} = \frac{1}{4\pi}, \qquad (98) $$

for any $d > 2$. More details of the perturbative expansion are in Appendix A.

## IV. 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 analog of the bulk Schwinger-Keldysh contour with the horizon cap for the gravitational background itself. This is straightforward. The metric dual to the Weyl scaled Bjorken flow [given by Eqs. (89) and (92)] 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 = \epsilon_0^{-1/d}$ 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 $\sigma$ 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 the absence of additional sources (and with the same boundary metric on the two arms) is exactly zero, as should be the case.<fnref n="17" />

There is a crucial subtlety to this rather simple construction. We should worry about the residual gauge symmetries

<footnote n="17">¹⁷In the field theory, this is equivalent to the statement that $W[J_1 = J_2 = 0] = 0$.</footnote>
at first and higher orders in the late proper time expansion. In what follows, we 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 = \epsilon_0^{-1/d}$ 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 = \epsilon_0^{1/d}$ up to third order in the proper time expansion in the case of AdS$_4$ 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.

At the outset, we repeat our emphasis that the Weyl rescaling (and hence the dual bulk diffeomorphism) is determined purely by the late proper time perfect flow regime that is

mapped to that of a constant temperature in a nontrivial background metric in the field theory (which is also the boundary metric of the dual black brane geometry after the bulk diffeomorphism). This should be evident from our construction of the bulk dual of the Weyl-transformed Bjorken flow (in the $v$ and $\sigma$ coordinates) in the previous section. Both these Weyl transformation and the background metric of the field theory do not receive any correction at first and higher orders in the derivative (i.e., large proper time) expansion. Nevertheless, with such a Weyl transformation (bulk diffeomorphism) we will be able to ensure that the Klein-Gordon equation of the bulk scalar takes the same form at first and higher orders as that at the leading order but with source terms, and that the leading near-horizon behavior is exactly the same as that in the static black brane geometry (in the $v$ and $\sigma$ coordinates) to all orders in the large proper time expansion when the horizon cap is pinned to the evolving event horizon. This will be instrumental in establishing the horizon cap prescription in the hydrodynamic tail of the Bjorken flow.

To see the main advantage of working in the $v$ and $\sigma$ coordinates, note that the explicit form of the leading order metric [as evident from Eqs. (89) and (92)] is

$$ds^2 = -\frac{2}{v^2} \frac{d-1}{d-2} dv d\sigma - \frac{1}{v^2} \left( \frac{(d-1)^2}{(d-2)^2} [1 - v^d \epsilon_0 + O(\sigma^{-1})] + \frac{2(d-1)v}{(d-2)^2 \sigma} \right) d\sigma^2 + \frac{1}{\tau_0^2} \left(1 + \frac{\sigma}{v}\right)^2 [1 + O(\sigma^{-1})] d\hat{\zeta}^2 + \frac{1}{v^2} \left(\frac{\sigma}{\tau_0}\right)^{-\frac{2}{d-2}} [1 + O(\sigma^{-1})] ds_{\perp}^2. \quad (99)$$

It is useful to define the comoving momenta which depend on $\sigma$:

$$\kappa_L = k_L \frac{\tau_0}{\sigma}, \quad \vec{k}_T = \vec{k}_T \left( \frac{\tau_0}{\sigma} \right)^{-\frac{1}{d-2}}, \quad (100)$$

such that it takes care of the longitudinal expansion and transverse contraction of the boundary metric. A natural ansatz for the scalar field consistent with boost invariance of the background geometry is

$$\Phi(v, \sigma, \hat{\zeta}, \vec{x}_{\perp}) \approx e^{i\left(-\frac{d-1}{d-2}\omega\sigma + k_L\hat{\zeta} + \vec{k}_T \cdot \vec{x}_{\perp}\right)} f(v, \omega, \kappa_L, \vec{k}_T), \quad (101)$$

where remarkably there is no explicit $\sigma$ dependence in $f$ (only implicitly through $\kappa_L$ and $\vec{k}_T$), while in the phase we have the usual fixed conjugate momenta $k_L$ and $\vec{k}_T$. This is exact at leading order and can be further corrected to obtain a systematic expansion of the equations of motion in powers of $\sigma^{-1}$, as shown below. Particularly, the Klein-Gordon equation ($\Box - m^2)\Phi = 0$ for the bulk scalar field at the leading order $\sigma^0$ is then$^{18}$

$$\mathcal{D}f = O(\sigma^{-1}), \quad (102)$$

where

$$\mathcal{D} = v^2(1 - v^d \epsilon_0) \partial_v^2 - v(d - 1 - 2iv\omega + v^d \epsilon_0) \partial_v - \{m^2 + v[v(\kappa_L^2 + \kappa_T^2) + i(d - 1)\omega]\}. \quad (103)$$

Remarkably, the left-hand side of (102) is just the Klein-Gordon equation for the massive bulk scalar in the AdS$_{d+1}$ static black brane geometry (22) at temperature $\beta^{-1}$ given by (97), with $v$ substituting for $r$, and $\omega$, $\kappa_L$, and $\vec{k}_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 = \vec{k}_T = 0$) and massless case, this result reduces to the observation made by Janik and Peschanski [23] in the context of the transients of the Bjorken flow.

Note that the dependence of $f$ in Eq. (101) on the comoving momenta implies that we do not have separation of variables, but this should not be expected as the background at late time corresponds to an expanding boost-invariant perfect fluid, which has no timelike Killing vector. Nevertheless, the map to the Laplacian of a static black hole geometry is possible at leading order because the boost invariant perfect fluid is given only by a time-dependent temperature.

The ansatz (101) can be corrected to incorporate the viscous and all higher-order corrections to the gravitational Bjorken flow background systematically while ensuring the regularity of the horizon cap after fixing the residual gauge

<footnote>$^{18}$To obtain this we should substitute $k_L$ and $\vec{k}_T$ in (101) by the re-defined momenta (100) and remember to differentiate $f$ with respect to $\sigma$ as the redefined momenta (100) depend on $\sigma$.</footnote>
freedom perturbatively in late proper time expansion. This modified ansatz involves an expansion in $\sigma^{-1}$ with coefficients which are functions of $v$, $\omega$, and the comoving momenta, such that the equations of motion can be solved systematically in $\sigma^{-1}$ expansion as well. Obviously, this implies that we also take into account the implicit dependence of the comoving momenta on $\sigma$ while obtaining the equations of motion order by order in $\sigma^{-1}$.

Explicitly, the ansatz for the bulk scalar field in the full complexified space-time with 1 and 2 labeling the sheets of the bulk space-time ending at the forward and backward arms of the Schwinger-Keldysh contour, respectively, at their boundaries is

$$ \begin{pmatrix} \Phi_1(v, \sigma, \hat{\zeta}, \vec{x}_{\perp}) \\ \Phi_2(v, \sigma, \hat{\zeta}, \vec{x}_{\perp}) \end{pmatrix} = \int d\omega dk_L d^{d-2}k_T e^{i(k_L \hat{\zeta} + \vec{k}_T \cdot \vec{x}_{\perp})} e^{-i \frac{d-1}{d-2} \omega \sigma} \\ \times \left(\frac{\sigma}{\tau_0}\right)^{i\gamma_0(\omega/\epsilon_0^{1/d})} \sum_{n=0}^{\infty} (\epsilon_0^{1/d} \sigma)^{-n} \\ \times \mathcal{M}_n(v, \omega, \kappa_L, \vec{k}_T) \cdot \begin{pmatrix} p(\omega, k_L, k_T) \\ q(\omega, k_L, k_T) \end{pmatrix}, \quad (104) $$

where

$$ \mathcal{M}_n(v, \omega, \kappa_L, \vec{k}_T) = \begin{pmatrix} \phi_{n,in}(v, \omega, \kappa_L, \vec{k}_T) & \phi_{n,out}(v, \omega, \kappa_L, \vec{k}_T) \\ \phi_{n,in}(v, \omega, \kappa_L, \vec{k}_T) & e^{\beta\omega} \phi_{n,out}(v, \omega, \kappa_L, \vec{k}_T) \end{pmatrix}. \quad (105) $$

Above, in (104), a nonanalytic term proportional to $\sigma^{i\gamma_0}$ has been introduced with $\gamma_0$ being a new dimensionless function of only $\omega/\epsilon_0^{1/d}$. This nonanalytic 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 $\omega$, (104) has the appearance of a trans-series in $\sigma^{-1}$ with $\omega$ characterizing the continuous instanton exponent. Note $p$ and $q$ are functions of $\omega$, $k_L$, and $k_T$ (and not the redefined momenta $\kappa_L$ and $\vec{k}_T$), since we are integrating over $\omega$, $k_L$, and $k_T$, and utilizing the superposition principle to obtain the general solution with right behavior at the horizon cap. However $\mathcal{M}_n$ depend on $\kappa_L$ and $\vec{k}_T$, and therefore they have nontrivial derivatives with respect to $\sigma$. The coefficients $p$ and $q$, which are functions of the ordinary momenta ($k_L$ and $\vec{k}_T$), and coefficients of the ingoing and outgoing modes respectively (see more below), are determined by imposing Dirichlet boundary conditions at the two boundaries.

The structure of $\mathcal{M}_n$ in (105) 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 ingoing and outgoing at the horizon, and satisfying the normalization conditions given by (25). Explicitly, near the perturbative horizon $v = \epsilon_0^{1/d}$ where the horizon cap is located, they take

the forms<fnref n="19" />

$$ \begin{align} \phi_{0,in} &= 1 + \sum_{k=1}^{\infty} p_{0,k} (\epsilon_0^{-1/d} - v)^k, \nonumber \\ \phi_{0,out} &= (\epsilon_0^{-1/d} - v)^{i \frac{2\omega}{d\epsilon_0^{1/d}}} \left[ 1 + \sum_{k=1}^{\infty} q_{0,k} (\epsilon_0^{-1/d} - v)^k \right], \tag{106} \end{align} $$

respectively. Note we have used (97) to set $\beta\omega/(2\pi) = 2\omega/(d\epsilon_0^{1/d})$. For $n \ge 1$, $\phi_{n,in}$ and $\phi_{n,out}$ represent corrections to these zeroth-order solutions, as discussed below. In (105), we have assumed that $\phi_{n,in}$ and $\phi_{n,out}$ have the same behavior at the horizon cap $v = v_h = \epsilon_0^{-1/d}$ for $n \ge 1$ as 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.<fnref n="20" />

At higher orders, the equations determining $\phi_{n,in}$ and $\phi_{n,out}$ can be obtained from expanding $(\Box - m^2)\Phi = 0$ in the late proper time expansion after substituting $k_L$ and $\vec{k}_T$ in (104) by the redefined momenta (100), and isolating the $\sigma^{-n}$ term. We obtain that $\phi_{n,in}$ and $\phi_{n,out}$ satisfies the linear inhomogeneous ordinary differential equations

$$ \mathcal{D}\phi_{n,in} = \mathcal{S}_{n,in}, \quad \mathcal{D}\phi_{n,out} = \mathcal{S}_{n,out}, \quad (107) $$

with $\mathcal{D}$ being the same operator (103), which is simply that corresponding to the Klein-Gordon equation for the static black brane at temperature $\beta^{-1}$ given by (97), at all orders. The sources $\mathcal{S}_{n,in}$ and $\mathcal{S}_{n,out}$ are functions of $v$, $\omega$, $\kappa_L$, and $\vec{k}_T$, and depend on $\phi_{m,in}$ and $\phi_{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 (92) of the background metric, with $i \le n$. Note that the source is linear in the bulk field $\Phi$, and therefore splits into $\mathcal{S}_{n,in}$ and $\mathcal{S}_{n,out}$ at each order in the $\sigma^{-1}$ expansion for $n \ge 1$. Since $\phi_{n,in}$ for $n \ge 1$ corrects $\phi_{0,in}$, we include the particular solution with only $\phi_{m,in}$ (and $m < n$) appearing in the source term $\mathcal{S}_{n,in}$ in it. The particular solutions sourced by $\phi_{m,out}$ with $0 \le m < n$ which appear in $\mathcal{S}_{n,out}$ are added by definition to $\phi_{n,out}$ similarly.

The regularity of the horizon cap implies that, at first and higher orders in the proper time expansion, we should have

$$ \begin{align} & \lim_{v \to \epsilon_0^{-1/d}} (\epsilon_0^{-1/d} - v)^{-i \frac{2\omega}{d\epsilon_0^{1/d}}} \phi_{n,out}(v, \omega, \kappa_L, \vec{k}_T) \nonumber \\ &= -i\gamma_{n,out}(\omega, \kappa_L, \vec{k}_T), \nonumber \\ & \lim_{v \to \epsilon_0^{-1/d}} \phi_{n,in}(v, \omega, \kappa_L, \vec{k}_T) = -i\gamma_{n,in}(\omega, \kappa_L, \vec{k}_T), \tag{108} \end{align} $$

<footnote n="19">¹⁹ $\phi_{0,in} \to 1$ and</footnote>

<footnote>$$(\epsilon_0^{-1/d} - v)^{-i \frac{2\omega}{d\epsilon_0^{1/d}}} \phi_{0,out} \to 1$$</footnote>

<footnote>as $v \to \epsilon_0^{-1/d}$ are just normalizations.</footnote>

<footnote n="20">²⁰ 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.</footnote>
with $\gamma_{n,out}$ and $\gamma_{n,in}$ as new functions of $\omega$, $\kappa_L$, and $\bar{\kappa}_T$ for $n \ge 1$. This will ensure that the analytical dependence of the field on $v$ at the horizon cap $v = v_h = \epsilon_0^{-1/d}$ is the same to all orders in the proper time expansion. Note that at the zeroth order, the analogous conditions for $\phi_{0,in}$ and $\phi_{0,out}$ are simply given by 1 on the RHS in both cases by choice [see (106)]. We will show that the $\gamma_{n,out}$ and $\gamma_{n,in}$ can be determined uniquely for $n \ge 1$ via horizon cap regularity and field theory identities.

First, note that (108) implies that $\phi_{n,in}$ and $\phi_{n,out}$ can be written in the form

$$ \phi_{n,in} = \phi_{n(p),in} - i\gamma_{n,in}\phi_{0,in}, \\ \phi_{n,out} = \phi_{n(p),out} - i\gamma_{n,out}\phi_{0,out} \quad (109) $$

for $n \ge 1$ such that $\phi_{n(p),in}$ and $\phi_{n(p),out}$ are the particular solutions of the inhomogeneous ordinary differential equations (107). Both $\phi_{n(p),in}$ and $\phi_{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 \approx \epsilon_0^{1/d}$, we explicitly find that $\phi_{n(p),in}$ and $\phi_{n(p),out}$ behave as

$$ \lim_{v \to \epsilon_0^{-1/d}} (\epsilon_0^{-1/d} - v)^{-i \frac{2\omega}{d\epsilon_0^{1/d}}} \phi_{n(p),out} = O(\epsilon_0^{-1/d} - v), \\ \lim_{v \to \epsilon_0^{-1/d}} \phi_{n(p),in} = O(\epsilon_0^{-1/d} - v), \quad (110) $$

respectively. Equivalently,

$$ \phi_{n(p),out} = (\epsilon_0^{-1/d} - v)^{i \frac{2\omega}{d\epsilon_0^{1/d}}} \sum_{k=1}^{\infty} q_{n,k} (\epsilon_0^{-1/d} - v)^k, \\ \phi_{n(p),in} = \sum_{k=1}^{\infty} p_{n,k} (\epsilon_0^{-1/d} - v)^k, \quad (111) $$

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 $\alpha_n$ appearing in the background metric and $\gamma_{n,out}$ are chosen appropriately to cancel double and single poles appearing in the equation of motion (107) at the horizon cap for each $n \ge 1$. Furthermore, $\alpha_n$ are simply numerical constants (as

they are defined to be), and $\gamma_{n,out}$ are linear functions of $\omega/\epsilon_0^{1/d}$ only. Thus the outgoing solutions appearing in $M_n$ in our ansatz (104) are determined uniquely. As discussed before and explicitly shown again in the next section, this completely determines the advanced propagator of the Bjorken flow. The nonequilibrium 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 $\gamma_{n,in}$ uniquely as well for $n \ge 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 (110) [equivalently (111)] implies that $\gamma_{n,in}$ and $\gamma_{n,out}$ are simply the coefficients of the homogeneous solutions of the equations of motion (107) for $n \ge 1$. Therefore, $\gamma_{n,in}$ and $\gamma_{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 (110) [equivalently (111)] at the $n$th order determines $\gamma_{n-1,out}$ along with the gauge parameter $\alpha_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 $\gamma_{n,in}$ for $n \ge 1$. However, $\gamma_0$ in the ansatz (104) appears in both the ingoing and outgoing modes. Our construction passes a significant consistency test that the same function $\gamma_0(\omega/\epsilon_0^{1/d})$ determines the homogeneous transients (sourceless solutions which are ingoing at the horizon) with the argument $\omega/\epsilon_0^{1/d}$ taking values corresponding to the appropriately rescaled quasinormal mode frequencies of the static black hole, as will be discussed in Sec. V E. This is remarkable as we determine the function $\gamma_0$ analytically by imposing the regularity condition (111) 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 $\gamma_0$ and $\alpha_1$ at the first order in the proper time expansion. This requires the first-order correction to the background metric given by (93) and (94). We find that the equation of motion (107) for $\phi_{1,out}$ explicitly takes the following form near the horizon cap $v \approx \epsilon_0^{1/4}$ up to overall proportionality factors:

$$ \frac{1}{(\epsilon_0^{1/4} v - 1)^2} \frac{\omega}{4\epsilon_0^{1/4}} (\alpha_1 + 3)(\omega\epsilon_0^{-1/4} + 2i) + \frac{1}{(\epsilon_0^{1/4} v - 1)} \frac{i\omega}{\epsilon_0^{1/4}} \times \left( \frac{2\kappa_L^2 + 2\kappa_T^2 - m^2 - 3\omega^2 + 240\omega\epsilon_0^{1/4}}{24\epsilon_0^{1/2}} (\alpha_1 + 3) + \frac{(16 + 5\alpha_1)\omega^2 - 4\gamma_0\omega\epsilon_0^{1/4}}{\epsilon_0^{1/2}} \right) + \dots, \quad (112) $$

with $\dots$ denoting terms which are regular at $v = \epsilon_0^{1/4}$ if $\phi_{1,out}$ is of the form (111). See Appendix C for more details. We readily see that in order to have a solution of the desired form (110) we must impose

$$ \alpha_1 = -3, \quad \gamma_0 = \frac{\omega}{4\epsilon_0^{1/4}} \quad (113) $$

so that the double- and single-pole terms of the equation of motion appearing at the horizon cap vanish. The double-pole term determines $\alpha_1$ and the single pole term determines $\gamma_0$. As claimed before, we find that $\gamma_0$ is indeed a simple linear function of $\omega/\epsilon_0^{1/4}$ (rather just proportional to it) while the gauge parameter $\alpha_1$ is a numerical constant as it should be.
At the second order in the proper time expansion, similarly $\alpha_2$ and $\gamma_1$ are determined by the vanishing of the double and single pole terms in the equation of motion for $\phi_{2,out}$, respectively. Here we have to utilize the explicit second-order correction to the background metric given in Appendix A. Explicitly,

$$ \alpha_2 = \frac{1}{72}(20 + 9\pi - 12 \ln 2), \\ \gamma_1 = \frac{1}{16}\left(4i - \frac{9\pi + 4 - 24 \ln 2}{9} \frac{\omega}{\varepsilon_0^{1/4}}\right). \quad (114) $$

We find once again $\alpha_2$ is just a numerical constant as it should be and $\gamma_1$ is a linear function of $\omega/\varepsilon_0^{1/4}$.

It is indeed crucial that the $\gamma_{n,out}$ for $n \ge 0$ are functions of $\omega$ only and are independent of $\kappa_L$ and $\bar{\kappa}_T$ which depend on $\sigma$. Otherwise, the central assumption (108) [and thus (111)] 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 (111) remains valid to all orders even if $\gamma_{n,in}$ depend on $\kappa_L$ and $\bar{\kappa}_T$. The coefficients $p_{n,k}$ in (111) also involve derivatives of $\gamma_{m,in}$ with respect to $\kappa_L$ and $\bar{\kappa}_T$ with $m < n$. See Appendix C for more details.

Finally, we have explicitly verified that the values of the residual gauge parameters $\alpha_1$ and $\alpha_2$ are such that the event horizon is pinned to the horizon cap $v = v_h = \varepsilon_0^{-1/d}$ 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.<fnref n="21" /> 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. V E.

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 ingoing and outgoing modes respectively, and the cross term $S_{in-out}$. As in the thermal case discussed in Sec. II, $S_{in-in} = 0$. The ingoing 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 ingoing mode by setting $q = 0$ in (104), 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 Sec. II. $S_{out-out}$ potentially has a single pole $(v_h - v)^{-1}$ divergence which we denote as $S_\epsilon$. Explicitly,

$$ S_\epsilon \propto \int d\omega \int dk_L \int d^{d-2}k_T \oint_\epsilon dv \sqrt{-G} (G^{vv} \partial_v \phi_{n,out}^* \partial_v \phi_{n,out} + G^{v\sigma} (\partial_v \phi_{n,out}^* \partial_\sigma \phi_{n,out} + \partial_v \phi_{n,out} \partial_\sigma \phi_{n,out}^*) + \dots) \quad (115) $$

at the $n$th order in the late proper time expansion. We have verified that $S_\epsilon = 0$ for the solution with the regular behavior at the horizon cap given by (108), obtained for the appropriate choices of $\alpha_n$ and $\gamma_{n-1,out}$ as discussed above. One can also check that terms like $\phi_{k,out}^* \phi_{l,out}$ with $k+l \le 2n$ also do not contribute to $S_\epsilon$. However, our arguments for $S_{out-out} = 0$ for the thermal case in Sec. II 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.

Since we preserve the horizon cap regularity at each order, the analytic continuation of the outgoing mode across the horizon cap at each order works exactly in the same way as in the case of thermal equilibrium, i.e., $\phi_{n,out}(\omega, v, \kappa_L, \kappa_T)$ picks up a factor of $e^{\beta\omega}$, as evident from (109) and (111). Repeating the argument in Sec. II, each boundary term in $S_{out-out}$ involves one $\phi_b(\omega, \kappa_L, \kappa_T)$ and another $\phi_b(-\omega, \kappa_L, \kappa_T)$ (with $b$ standing for the boundary value) or the corresponding boundary values of the radial derivatives. Also, the contribution from the forward and backward arms come with opposite signs. While $\phi_b(\omega, \kappa_L, \kappa_T)$ picks up a $e^{\beta\omega}$ factor via analytic continuation across the horizon cap, $\phi_b(-\omega, \kappa_L, \kappa_T)$ picks up a $e^{-\beta\omega}$ factor, and these multiply to unity. Therefore, the boundary term contributions from the two arms cancel out resulting in $S_{out-out} = 0$.

Finally, 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.

As shown in the following section, the boundary correlation functions obtained solely from $S_{in-out}$ ensure that the (nonequilibrium) retarded correlation function is given always by the linear causal response. Furthermore, a major consistency check is that we reproduce the homogeneous transients as poles of the retarded Green's function in complexified $\omega$ not only at the leading order as computed in Ref. [23], but also at the subleading orders as shown in Appendix E and discussed further below. We will discuss more consistency tests in Sec. V E.

## V. THE REAL-TIME OUT-OF-EQUILIBRIUM CORRELATION FUNCTIONS

In the previous section, we have shown that via the bulk diffeomorphism dual to the Weyl transformation we can obtain the solution of the bulk scalar field which has the same analytic thermal nature of the near-horizon modes to all orders in the (hydrodynamic) large proper time expansion. Furthermore, this solution is unique up to the coefficients that give the leading near-horizon behavior of the ingoing modes, and

<footnote n="21">²¹Although we do not have a rigorous proof to all orders, we believe that this follows from the general result that the event horizon is generated by null geodesics which also determine the singularities arising in the equation of motion at the horizon cap.</footnote>
the on-shell action to all orders is the product of the ingoing and outgoing modes, as the terms which are quadratic in the ingoing modes or the outgoing modes cancel exactly between the radial arms corresponding to the forward and backward parts of the time contour due to the thermal nature of their leading near-horizon behavior. Here we compute the Schwinger-Keldysh correlation functions from the on-shell action of the bulk scalar field in an appropriate late-time expansion. We show that the correlation functions are consistent with the field theory identities which hold even for out-of-equilibrium states and one of them can be used to fix the coefficients that give the leading near-horizon behavior of the ingoing modes and thus determine the bulk scalar field solution uniquely. Furthermore, we show that the all-order hydrodynamic correlation function has a matrix factorized form giving a bilocal generalization of the corresponding form of the thermal correlation functions.

In Sec. VA, we discuss the general identities satisfied by the Schwinger-Keldysh correlation functions and discuss how they transform under Weyl transformations. In Sec. VB, we obtain the Schwinger-Keldysh correlation functions from the on-shell action and establish its matrix factorized form which is consistent with the field theory identities. We also show that one of the latter identities can determine the bulk solution of the dual scalar field uniquely and therefore the Schwinger-Keldysh correlation functions. In Sec. VC, we obtain the correlation functions in the perfect fluid limit and show that it is thermal after space-time reparametrizations up to overall Weyl factors. In Sec. VD, we establish the systematic late-time expansion of the correlation functions. Finally, in Sec. VE, we discuss consistency checks, including the feature find that the retarded propagator given by the ingoing modes reproduces the normalizable bulk solutions at complex frequencies (which map to quasinormal modes of the static black brane [23]) with vanishing sources to all orders although the relevant phase factors at first and higher orders are determined as functions of the frequency via the near-horizon behavior of the outgoing modes at real frequencies. The latter implies nontrivially that one can indeed obtain the normalizable bulk solutions dual to the nonequilibrium transients of the dual field theory to all orders in the large proper time expansion from the retarded correlation function generalizing how we can obtain the quasinormal modes (thermal transients) as poles of the thermal retarded correlation function.

### A. Some useful relations and their consequences

Some crucial identities are valid for the Schwinger-Keldysh correlation functions even in out-of-equilibrium states. The first such identity of interest is

$$ 
\begin{aligned}
G_R(x_1, x_2) &= -i\theta(x_1^0 - x_2^0)\mathrm{Tr}(\rho[\hat{O}(x_1), \hat{O}(x_2)]) \\
&= -i[\mathrm{Tr}(\rho\mathrm{T}(\hat{O}(x_1)\hat{O}(x_2))) - \mathrm{Tr}(\rho\hat{O}(x_2)\hat{O}(x_1))] \\
&= G_{11}(x_1, x_2) - G_{12}(x_1, x_2).
\end{aligned}
\quad (116) $$

Similarly, we can arrive at the identity

$$ 
\begin{aligned}
G_A(x_1, x_2) &= i\theta(x_2^0 - x_1^0)\mathrm{Tr}(\rho[\hat{O}(x_1), \hat{O}(x_2)]) \\
&= -i[\mathrm{Tr}(\rho\mathrm{T}(\hat{O}(x_1)\hat{O}(x_2))) - \mathrm{Tr}(\rho\hat{O}(x_1)\hat{O}(x_2))] \\
&= G_{11}(x_1, x_2) - G_{21}(x_1, x_2).
\end{aligned}
\quad (117) $$

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

$$ 
\begin{aligned}
G_{11}(x_1, x_2) &= G_{11}(x_2, x_1), & G_{12}(x_1, x_2) &= G_{21}(x_2, x_1), \\
G_{22}(x_1, x_2) &= G_{22}(x_2, x_1).
\end{aligned}
\quad (118) $$

Finally, it is obvious also that in any arbitrary (nonequilibrium) state

$$ G_R(x_1, x_2) = G_A(x_2, x_1). \quad (119) $$

We show that these identities imply the following for the horizon cap of the Bjorken flow. First, we use (116) to show that the retarded correlation function is always obtained from the ingoing mode and (117) to show that the advanced correlation function is always obtained from the outgoing mode. Second, due to (119), we can uniquely determine $\gamma_{n,in}$ for $n \ge 1$. These coefficients cannot be determined by the regularity condition (108) of the horizon cap unlike $\gamma_{n,out}$, as discussed previously, but together with $\gamma_{n,out}$ they provide unique solution of the bulk scalar field for specified sources at the two boundaries via our ansatz (104). We also see that, by satisfying (119), we can also ensure the validity of (118).

Another important issue is the Weyl transformation. Consider a background metric $g_{\mu\nu}$ and its Weyl rescaled version $\Omega^2(x)g_{\mu\nu}$ such that $\Omega(x)$ is a nonvanishing function. Disregarding the Weyl anomaly, the correlations functions of a scalar primary operator $\hat{O}$ of conformal dimension $\Delta_O$ in a conformal field theory living in these two background metrics would be related by

$$ 
\begin{aligned}
\mathrm{Tr}(\tilde{\rho}\hat{O}(x_1) \cdots \hat{O}(x_n))_{\Omega^2 g_{\mu\nu}} \\
&= \Omega^{-\Delta_O}(x_1) \cdots \Omega^{-\Delta_O}(x_n) \mathrm{Tr}(\rho\hat{O}(x_1) \cdots \hat{O}(x_n))_{g_{\mu\nu}},
\end{aligned}
\quad (120) $$

where $\tilde{\rho} = U^\dagger \rho U$ with $U$ denoting the unitary operator implementing the Weyl transformation. In the context of the holographic Bjorken flow, the states $\rho$ and $\tilde{\rho}$ 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 Sec. III B.

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 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

$$ \Omega^{\Delta_O}(\tau_1)\Omega^{\Delta_O}(\tau_2)\tilde{G}(x_1, x_2), $$

with $\tilde{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 $\Omega$ given by (56) and (57). 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 $\hat{\xi}_1$ and $\hat{\xi}_2$ only through the relative separation $\hat{\xi}_1 - \hat{\xi}_2$ in rapidity due to the boost invariance of the background, and similarly only on $|\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|$ due to translation and rotation symmetries of the background in the transverse spatial plane.

We note that we can practically use the Weyl transformation to compute the correlation functions because we are in the boost-invariant hydrodynamic regime where the state is

characterized only by a proper time-dependent temperature. A Weyl rescaling makes the temperature constant at late time leading us to set up a perturbative derivative expansion controlled by a fixed dimensionful parameter. However, there will be contributions to the correlation functions in the form of a generalized trans-series which depend on the initial conditions and are not constructed from hydrodynamic data. Such contributions which represent transients leading to hydrodynamization of correlation functions are discussed in Sec. VI. In this section, we discuss the Schwinger-Keldysh correlation functions in the hydrodynamic regime using a well-defined late proper time expansion.

## B. General structure of hydrodynamic correlation functions

It is useful to define a new variable

$$s(\sigma) = \frac{d-1}{d-2}\sigma - \tilde{\gamma}_0\epsilon_0^{-1/d}\ln(\sigma/\tau_0), \qquad (121)$$

with $\tilde{\gamma}_0 = \gamma_0\epsilon_0^{1/d}/\omega$ being a numerical constant (independent of $\omega$). We have already seen that the smoothness of the horizon cap requires that, in the case of AdS$_5$, $\tilde{\gamma}_0 = 1/4$ [see (113)].$^{22}$ In terms of this variable, we can readily see from (104) that the non-normalizable mode of the scalar field in the gravitational background dual to the Bjorken flow takes the form

$$\begin{pmatrix} J_1(\sigma, \hat{\xi}, \vec{x}_{\perp}) \\ J_2(\sigma, \hat{\xi}, \vec{x}_{\perp}) \end{pmatrix} = \int d\omega dk_L d^{d-2}k_T e^{-i\omega s(\sigma)} e^{i(k_L \hat{\xi} + \vec{k}_T \cdot \vec{x}_{\perp})} \sum_{n=0}^{\infty} (\epsilon_0^{1/d} \sigma)^{-n} S_n(\omega, \kappa_L, \vec{k}_T) \cdot \begin{pmatrix} p(\omega, k_L, k_T) \\ q(\omega, k_L, k_T) \end{pmatrix}, \qquad (122)$$

and similarly the expectation value of the dual operator (up to state-independent contact terms) is

$$\begin{pmatrix} (O_1(\sigma, \hat{\xi}, \vec{x}_{\perp})) \\ (O_2)(\sigma, \hat{\xi}, \vec{x}_{\perp}) \end{pmatrix} = (2\Delta_o - d) \int d\omega dk_L d^{d-2}k_T e^{-i\omega s(\sigma)} e^{i(k_L \hat{\xi} + \vec{k}_T \cdot \vec{x}_{\perp})} \sum_{n=0}^{\infty} (\epsilon_0^{1/d} \sigma)^{-n} \mathcal{R}_n(\omega, \kappa_L, \vec{k}_T) \cdot \begin{pmatrix} p(\omega, k_L, k_T) \\ q(\omega, k_L, k_T) \end{pmatrix}, \qquad (123)$$

where $S_n$ and $\mathcal{R}_n$ can be defined from the asymptotic expansion of $M_n$:

$$M_n(v, \omega, \kappa_L, \vec{k}_T) = v^{d-\Delta_o} [S_n(\omega, \kappa_L, \vec{k}_T) + \dots] + v^{\Delta_o} [\mathcal{R}_n(\omega, \kappa_L, \vec{k}_T) + \dots], \qquad (124)$$

and the labels 1 and 2 stand for the sheets of the bulk space-time ending on the forward and backward arms of the Schwinger-Keldysh contour, respectively.

It is also useful to define

$$\mathcal{R}(\sigma, \omega, \kappa_L, \vec{k}_T) = \sum_{n=0}^{\infty} (\epsilon_0^{1/d} \sigma)^{-n} \mathcal{R}_n(\omega, \kappa_L, \vec{k}_T), \quad S(\sigma, \omega, \kappa_L, \vec{k}_T) = \sum_{n=0}^{\infty} (\epsilon_0^{1/d} \sigma)^{-n} S_n(\omega, \kappa_L, \vec{k}_T). \qquad (125)$$

Clearly (105) leads to

$$\begin{align}
S(\sigma, \omega, \kappa_L, \vec{k}_T) &= \begin{pmatrix} a(\sigma, \omega, \kappa_L, \vec{k}_T) & b(\sigma, \omega, \kappa_L, \vec{k}_T) \\ a(\sigma, \omega, \kappa_L, \vec{k}_T) & e^{\beta\omega}b(\sigma, \omega, \kappa_L, \vec{k}_T) \end{pmatrix}, \nonumber \\
\mathcal{R}(\sigma, \omega, \kappa_L, \vec{k}_T) &= \begin{pmatrix} A(\sigma, \omega, \kappa_L, \vec{k}_T) & B(\sigma, \omega, \kappa_L, \vec{k}_T) \\ A(\sigma, \omega, \kappa_L, \vec{k}_T) & e^{\beta\omega}B(\sigma, \omega, \kappa_L, \vec{k}_T) \end{pmatrix}, \tag{126}
\end{align}$$

where

$$\begin{align}
\phi_{\text{in}}(v, \sigma, \omega, \kappa_L, \vec{k}_T) &= \sum_{n=0}^{\infty} (\epsilon_0^{1/d} \sigma)^{-n} \phi_{n, \text{in}}(v, \sigma, \omega, \kappa_L, \vec{k}_T), \nonumber \\
\phi_{\text{out}}(v, \sigma, \omega, \kappa_L, \vec{k}_T) &= \sum_{n=0}^{\infty} (\epsilon_0^{1/d} \sigma)^{-n} \phi_{n, \text{out}}(v, \sigma, \omega, \kappa_L, \vec{k}_T) \tag{127}
\end{align}$$

<footnote>$^{22}$In the case of AdS$_4$, $\tilde{\gamma}_0 = 2/3$.</footnote>
gives us $a(\sigma, \omega, \kappa_L, \vec{k}_T)$, $A(\sigma, \omega, \kappa_L, \vec{k}_T)$, $b(\sigma, \omega, \kappa_L, \vec{k}_T)$, and $B(\sigma, \omega, \kappa_L, \vec{k}_T)$ from their following near boundary expansions:

$$ \begin{aligned} \phi_{\text{in}}(v, \sigma, \omega, \kappa_L, \vec{k}_T) &= v^{d-\Delta_o}[a(\sigma, \omega, \kappa_L, \vec{k}_T) + \dots] + v^{\Delta_o}[A(\sigma, \omega, \kappa_L, \vec{k}_T) + \dots], \\ \phi_{\text{out}}(v, \sigma, \omega, \kappa_L, \vec{k}_T) &= v^{d-\Delta_o}[b(\sigma, \omega, \kappa_L, \vec{k}_T) + \dots] + v^{\Delta_o}[B(\sigma, \omega, \kappa_L, \vec{k}_T) + \dots]. \end{aligned} \quad (128) $$

It should also be obvious that if we define $a_n, A_n, b_n$, and $B_n$ via

$$ \begin{aligned} \phi_{n,\text{in}}(v, \omega, \kappa_L, \vec{k}_T) &= v^{d-\Delta_o}[a_n(\omega, \kappa_L, \vec{k}_T) + \dots] + v^{\Delta_o}[A_n(\omega, \kappa_L, \vec{k}_T) + \dots], \\ \phi_{n,\text{out}}(v, \omega, \kappa_L, \vec{k}_T) &= v^{d-\Delta_o}[b_n(\omega, \kappa_L, \vec{k}_T) + \dots] + v^{\Delta_o}[B_n(\omega, \kappa_L, \vec{k}_T) + \dots], \end{aligned} \quad (129) $$

then

$$ \begin{aligned} a(\sigma, \omega, \kappa_L, \vec{k}_T) &= \sum_{n=0}^{\infty} (\varepsilon_0^{1/d} \sigma)^{-n} a_n(\omega, \kappa_L, \vec{k}_T), & A(\sigma, \omega, \kappa_L, \vec{k}_T) &= \sum_{n=0}^{\infty} (\varepsilon_0^{1/d} \sigma)^{-n} A_n(\omega, \kappa_L, \vec{k}_T), \\ b(\sigma, \omega, \kappa_L, \vec{k}_T) &= \sum_{n=0}^{\infty} (\varepsilon_0^{1/d} \sigma)^{-n} b_n(\omega, \kappa_L, \vec{k}_T), & B(\sigma, \omega, \kappa_L, \vec{k}_T) &= \sum_{n=0}^{\infty} (\varepsilon_0^{1/d} \sigma)^{-n} B_n(\omega, \kappa_L, \vec{k}_T). \end{aligned} \quad (130) $$

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_{\text{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 (61)] is (see Appendix D for more details)

$$ \begin{aligned} & \tilde{G}_{ab}(\sigma_1, \sigma_2, \hat{\zeta}_1 - \hat{\zeta}_2, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|) \\ &= \frac{1}{\sqrt{-\tilde{g}_1 \sqrt{-\tilde{g}_2}} \delta^2 S_{\text{on-shell}}[J_1, J_2]} (-)^{a+b} \\ &= \frac{(d-2)^2}{(d-1)^2} \int d\omega dk_L d^{d-2}k_T e^{-i\omega[s(\sigma_1) - s(\sigma_2)]} e^{ik_L(\hat{\zeta}_1 - \hat{\zeta}_2) + i\vec{k}_T \cdot (\vec{x}_{\perp 1} - \vec{x}_{\perp 2})} \tilde{G}_{ab}(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T), \end{aligned} \quad (131) $$

where [compare with the thermal case (42)—recall that $\sigma_3$ is the Pauli matrix]

$$ \begin{aligned} \hat{G}(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) &= \frac{2\Delta_o - d}{2} (s'(\sigma_2)\sigma_3 \cdot \mathcal{R}(\sigma_1, \omega, \kappa_{L1}, \vec{k}_{T1}) \cdot \mathcal{S}^{-1}(\sigma_2, \omega, \kappa_{L2}, \vec{k}_{T2}) \\ &\quad + (\text{transpose } \sigma_1 \leftrightarrow \sigma_2, \omega \rightarrow -\omega, \kappa_{L1} \leftrightarrow -\kappa_{L2}, \vec{k}_{T1} \leftrightarrow -\vec{k}_{T2})), \end{aligned} \quad (132) $$

with

$$ \kappa_{L1} = k_L \frac{\tau_0}{\sigma_1}, \quad \kappa_{L2} = k_L \frac{\tau_0}{\sigma_2}, \quad \vec{k}_{T1} = \vec{k}_T \left(\frac{\tau_0}{\sigma_1}\right)^{-1/(d-2)}, \quad \vec{k}_{T2} = \vec{k}_T \left(\frac{\tau_0}{\sigma_2}\right)^{-1/(d-2)}. \quad (133) $$

The transpose above denotes matrix transposition, $\leftrightarrow$ denotes exchange operation and $\rightarrow$ replacement. The second term in (132) is produced by the symmetrization due to the differentiation in (131).

Explicitly [compare with the thermal case given by Eq. (43)],

$$ \begin{aligned} \hat{G}_{11}(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) &\equiv \hat{G}_R(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T)[1+n(\omega)] - \hat{G}_A(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T)n(\omega), \\ \hat{G}_{12}(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) &\equiv [\hat{G}_R(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) - \hat{G}_A(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T)]n(\omega), \\ \hat{G}_{21}(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) &\equiv [\hat{G}_R(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) - \hat{G}_A(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T)][1+n(\omega)], \\ \hat{G}_{22}(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) &\equiv \hat{G}_R(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T)n(\omega) - \hat{G}_A(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T)[1+n(\omega)], \end{aligned} \quad (134) $$

with $n(\omega) = 1/(e^{\beta\omega} - 1)$ is the Bose-Einstein distribution and

$$ \begin{aligned} \hat{G}_R(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) &\equiv \frac{(d-2)^2}{(d-1)^2} (2\Delta_o - d) \frac{s'(\sigma_1) + s'(\sigma_2)}{2} \frac{A(\sigma_1, \omega, \kappa_{L1}, \vec{k}_{T1})}{a(\sigma_2, \omega, \kappa_{L2}, \vec{k}_{T2})}, \\ \hat{G}_A(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) &\equiv \frac{(d-2)^2}{(d-1)^2} (2\Delta_o - d) \frac{s'(\sigma_1) + s'(\sigma_2)}{2} \frac{B(\sigma_1, \omega, \kappa_{L1}, \vec{k}_{T1})}{b(\sigma_2, \omega, \kappa_{L2}, \vec{k}_{T2})}. \end{aligned} \quad (135) $$

Above $\equiv$ denotes equality up to terms which vanish after the frequency and momentum integration in (131) and $\beta$ is given by (97).

In order that (134) and (135) follow from (132), and also for the general identities (119) and (118) to be satisfied, we should have

$$ \frac{A(\sigma_1, \omega, \kappa_{L1}, \vec{k}_{T1})}{a(\sigma_2, \omega, \kappa_{L2}, \vec{k}_{T2})} \equiv \frac{B(\sigma_2, -\omega, -\kappa_{L1}, -\vec{k}_{T1})}{b(\sigma_1, -\omega, -\kappa_{L2}, -\vec{k}_{T2})}, \quad (136) $$
i.e.,

$$ \hat{G}_R(\sigma_1, \sigma_2, \omega, \kappa_{L1}, \vec{k}_{T1}) \equiv \hat{G}_A(\sigma_2, \sigma_1, -\omega, -\kappa_{L2}, -\vec{k}_{T2}). \quad (137) $$

As mentioned, we show in Sec. V D that the above can be satisfied by appropriate choices of $\gamma_{n,ins}$ at each order in the late proper time expansion where we go to large values of the average reparametrized proper time $\bar{\sigma} = (\sigma_1 + \sigma_2)/2$ with fixed difference $\sigma_r = \sigma_1 - \sigma_2$. Thus we see that the correlation functions of the hydrodynamic Bjorken flow has a hidden and simple bilocal thermal structure to all orders in the late proper time expansion. We show this satisfies crucial consistency tests.

Using (134) and (135), and the identities (116) and (117), which are valid out-of-equilibrium, we readily see that the actual retarded and advanced propagators of the Weyl-transformed Bjorken flow is

$$ \begin{aligned} \tilde{G}_R(\sigma_1, \sigma_2, \hat{\zeta}_1 - \hat{\zeta}_2, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|) &= \frac{(d-2)^2}{(d-1)^2} (2\Delta_o - d) \frac{s'(\sigma_1) + s'(\sigma_2)}{2} \int d\omega dk_L d^{d-2}k_T e^{-i\omega[s(\sigma_1) - s(\sigma_2)]} \\ & \quad \times e^{ik_L(\hat{\zeta}_1 - \hat{\zeta}_2) + i\vec{k}_T \cdot (\vec{x}_{\perp 1} - \vec{x}_{\perp 2})} \frac{A(\sigma_1, \omega, \kappa_{L1}, \vec{k}_{T1})}{a(\sigma_2, \omega, \kappa_{L2}, \vec{k}_{T2})}, \\ \tilde{G}_A(\sigma_1, \sigma_2, \hat{\zeta}_1 - \hat{\zeta}_2, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|) &= \frac{(d-2)^2}{(d-1)^2} (2\Delta_o - d) \frac{s'(\sigma_1) + s'(\sigma_2)}{2} \int d\omega dk_L d^{d-2}k_T e^{-i\omega[s(\sigma_1) - s(\sigma_2)]} \\ & \quad \times e^{ik_L(\hat{\zeta}_1 - \hat{\zeta}_2) + i\vec{k}_T \cdot (\vec{x}_{\perp 1} - \vec{x}_{\perp 2})} \frac{B(\sigma_1, \omega, \kappa_{L1}, \vec{k}_{T1})}{b(\sigma_2, \omega, \kappa_{L2}, \vec{k}_{T2})}. \end{aligned} \quad (138) $$

Therefore, it follows from (135) 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-of-equilibrium outgoing mode as claimed before. Note that at the zeroth order, these results are automatic as it reduces to the thermal case described earlier.<fnref n="23" />

Finally, to obtain the Schwinger-Keldysh correlation function of the Bjorken flow, we should implement the Weyl transformation. Using (57), the Weyl transformation finally yields the correlation function of the Bjorken flow:

$$ G(\sigma_1, \sigma_2, \hat{\zeta}_1 - \hat{\zeta}_2, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|) = \left(\frac{\tau_0}{\sigma_1}\right)^{\frac{\Delta_0}{d-2}} \left(\frac{\tau_0}{\sigma_2}\right)^{\frac{\Delta_0}{d-2}} \tilde{G}(\sigma_1, \sigma_2, \hat{\zeta}_1 - \hat{\zeta}_2, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|). \quad (139) $$

### C. In the limit of the perfect fluid expansion

At a 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 space-time 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 correlation functions via appropriate space-time-reparametrizations. After undoing the Weyl transformation, we should obtain the Schwinger-Keldysh correlation functions of the perfect fluid expansion.

To obtain the correlation function in the perfect fluid expansion, we first take the limit in $\tilde{G}$ given by Eqs. (131)–(135), where we take the average reparametrized proper time coordinate $\bar{\sigma} = (1/2)(\sigma_1 + \sigma_2)$ to infinity keeping the relative reparametrized proper time coordinate $\sigma_r = \sigma_1 - \sigma_2$ fixed, and also $\hat{\zeta}_1 - \hat{\zeta}_2$ and $\vec{x}_{\perp 1} - \vec{x}_{\perp 2}$ fixed. In this limit,

$$ \kappa_{L1}, \kappa_{L2} \rightarrow \bar{\kappa}_L, \quad \bar{\kappa}_{T1}, \bar{\kappa}_{T2} \rightarrow \bar{\kappa}_T, \quad (140) $$

where

$$ \bar{\kappa}_L = k_L \frac{\tau_0}{\bar{\sigma}}, \quad \bar{\kappa}_T = \bar{k}_T \left(\frac{\tau_0}{\bar{\sigma}}\right)^{-1/(d-2)}. \quad (141) $$

As shown in the previous section, both $\phi_{in}$ and $\phi_{out}$ assume the form in the static black brane geometry, and therefore,

$$ \hat{G}_R \rightarrow \frac{A_0(\omega, \bar{\kappa}_L, \bar{\kappa}_T)}{a_0(\omega, \bar{\kappa}_L, \bar{\kappa}_T)}, \quad \hat{G}_A \rightarrow \frac{B_0(\omega, \bar{\kappa}_L, \bar{\kappa}_T)}{b_0(\omega, \bar{\kappa}_L, \bar{\kappa}_T)}, \quad (142) $$

given by the corresponding static black brane results and therefore all Schwinger-Keldysh correlation functions after the momentum integrals shown in (131) should assume the thermal form, i.e.,

$$ \tilde{G} \rightarrow G_{\beta} \left[ \sigma_1 - \sigma_2, (\hat{\zeta}_1 - \hat{\zeta}_2) \frac{\bar{\sigma}}{\tau_0}, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}| \left( \frac{\tau_0}{\bar{\sigma}} \right)^{\frac{1}{d-2}} \right], \quad (143) $$

<footnote n="23">²³The bilocal thermal structure of hydrodynamic correlation functions in holography was argued earlier in Ref. [48] via the Wigner transform. However, the arguments here were less rigorous and relied on the possibility of obtaining correlation functions utilizing only one copy of the nonequilibrium background where the interior of the perturbative horizon could be removed.</footnote>
with $G_β$ denoting the thermal correlation functions. To obtain the above form, we change variables in the momentum integrals from $k_L$ and $\tilde{k}_T$ to $\bar{\kappa}_L$ and $\bar{\kappa}_T$. Note the Jacobian of this transformation is identity. We also need to change the integration variable $\omega$ to $\omega = \omega(d-1)/(d-2)$ which yields a Jacobian $(d-2)/(d-1)$. Furthermore, we have used

$$s'(\sigma_1), s'(\sigma_2) \rightarrow \frac{d-1}{d-2}. \qquad (144)$$

The factor $(d-2)^2/(d-1)^2$ in (131) is canceled 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

$$G \rightarrow \left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{2\Delta_0}{d-2}} G_\beta \left(\frac{d-1}{d-2}(\sigma_1 - \sigma_2), (\hat{\zeta}_1 - \hat{\zeta}_2)\frac{\bar{\sigma}}{\tau_0}, |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|\left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{1}{d-2}}\right), \qquad (145)$$

the Schwinger-Keldysh correlation functions of the perfect fluid expansion. Since the static black brane has full rotational invariance in terms of the reparametrized space-time coordinates, we may also write

$$G \rightarrow \left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{2\Delta_0}{d-2}} G_\beta \left( \frac{d-1}{d-2}(\sigma_1 - \sigma_2), \sqrt{(\hat{\zeta}_1 - \hat{\zeta}_2)^2 \frac{\bar{\sigma}^2}{\tau_0^2} + |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|^2 \left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{2}{d-2}}} \right) \qquad (146)$$

for the correlation functions in the limit of the late time perfect fluid expansion. The explicit analytic forms of the thermal holographic propagators is known, see Refs. [33,49], for instance.

It is obvious that we are actually resumming over the associated $\bar{\sigma}$ factors in the spatial factors $\hat{\zeta}_1 - \hat{\zeta}_2$ and $\vec{x}_{\perp 1} - \vec{x}_{\perp 2}$. For brevity we denote the variables as $\sigma_r = \sigma_1 - \sigma_2$ (as defined before) and

$$\tilde{\zeta}_r = (\hat{\zeta}_1 - \hat{\zeta}_2)\frac{\bar{\sigma}}{\tau_0}, \quad \tilde{x}_{\perp r} = |\vec{x}_{\perp 1} - \vec{x}_{\perp 2}|\left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{1}{d-2}}. \qquad (147)$$

The late proper time expansion of the Schwinger-Keldysh correlation function then amounts to the following series:

$$G(x_1, x_2) = \left(\frac{\tau_0}{\bar{\sigma}}\right)^{\frac{2\Delta_0}{d-2}} \sum_{n=0}^{\infty} \frac{1}{\bar{\sigma}^n \epsilon_0^{n/d}} G_n(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r}) + \dots, \qquad (148)$$

with $\dots$ denoting trans-series type completion about which we discuss more below.

Since the choice of $\tau_0$ determines the effective final temperature $T = \beta^{-1}$ via (97), let us consider a scaling $\tau_0 \rightarrow \xi \tau_0$ to see if there is any ambiguity in the above result. Note that $G_β$ in (146) depends only on $T\sigma_r$ and $T(\tilde{\zeta}_r^2 + \tilde{x}_{\perp r}^2)^{1/2}$ as it is a thermal correlation matrix of a conformal field theory. Under this scaling, it is evident from (7) that $\sigma \rightarrow \xi^{\frac{1}{d-1}}\sigma$, and therefore $\sigma_r \rightarrow \xi^{\frac{1}{d-1}}\sigma_r$, $\bar{\sigma} \rightarrow \xi^{\frac{1}{d-1}}\bar{\sigma}$, $\tilde{\zeta}_r \rightarrow \xi^{\frac{1}{d-1}}\tilde{\zeta}_r$, and $\tilde{x}_{\perp r} \rightarrow \xi^{\frac{1}{d-1}}\tilde{x}_{\perp r}$. Also $\epsilon_0 \rightarrow \xi^{-\frac{d}{d-1}}\epsilon_0$ and $T \rightarrow \xi^{-\frac{1}{d-1}}T$ as discussed earlier in Sec. III A. Together these imply that $G_β$ and also $G_n$ in (148) are invariant under $\tau_0 \rightarrow \xi \tau_0$ since $T\sigma_r$ and $T\tilde{\zeta}_r$ and $T\tilde{x}_{\perp r}$ are invariant this scaling. However, the Weyl factor in (9) scales as $\xi^{\frac{2\Delta_0}{d-1}}$ implying that the dimensionless correlation function $\bar{\sigma}^{-2\Delta_0}G$ is invariant under the scaling of $\tau_0$.

### D. First and higher orders in the late proper time expansion

From the form of the Schwinger-Keldysh correlation functions given by Eqs. (131)–(135) which are valid to all orders, we can readily go beyond the perfect fluid limit systematically and construct the late proper time expansion (148). The contribution to the first-order correction comes from the following terms:

(1) The phase factor

$$e^{-i\omega[s(\sigma_1) - s(\sigma_2)]} \rightarrow e^{-i\omega\sigma_r \frac{d-1}{d-2}} \left[ 1 + i\omega\tilde{\gamma}_0\epsilon_0^{-1/d} \frac{\sigma_r}{\bar{\sigma}} + O(\bar{\sigma}^{-2}) \right]. \qquad (149)$$

(2) The Jacobian

$$s'(\sigma_1) \text{ and } s'(\sigma_2) \rightarrow \frac{d-1}{d-2} - \frac{\tilde{\gamma}_0}{\bar{\sigma}\epsilon_0^{1/d}} + O(\bar{\sigma}^{-2}). \qquad (150)$$

(3) Late-time expansion of the matrices $\mathcal{R}$ and $\mathcal{S}$ given by (126) in $\bar{\sigma}^{-1}$. This also includes the $\kappa_L$ and $\bar{\kappa}_T$ dependence. For instance, consider any generic function $f(\kappa_{L1})$ and $f(\bar{\kappa}_{T1})$,

$$f(\kappa_{L1}) = f(\bar{\kappa}_L) - \partial_{\bar{\kappa}_L} f(\bar{\kappa}_L) \frac{\sigma_r \bar{\kappa}_L}{2\bar{\sigma}}, \quad f(\bar{\kappa}_{T1}) = f(\bar{\kappa}_T) + \partial_{\bar{\kappa}_T} f(\bar{\kappa}_T) \frac{\sigma_r \bar{\kappa}_T}{2(d-2)\bar{\sigma}}, \qquad (151)$$

and similarly for $f(\kappa_{L2})$ and $f(\bar{\kappa}_{T2})$.
From (132), we obtain

$$ \begin{aligned} \hat{G}_1 = \frac{(2\Delta_0 - d)}{2} \frac{d-2}{d-1} \sigma_3 \cdot & \left[ \mathcal{R}_0 \cdot \mathcal{S}_1^{-1} + \mathcal{R}_1 \cdot \mathcal{S}_0^{-1} + (\mathcal{R}_0 \cdot \partial_{\bar{\kappa}_L} \mathcal{S}_0^{-1} - \partial_{\bar{\kappa}_L} \mathcal{R}_0 \cdot \mathcal{S}_0^{-1}) \frac{\sigma_r \bar{\kappa}_L}{2} \right. \\ & \left. - (\mathcal{R}_0 \cdot \partial_{\bar{\kappa}_T} \mathcal{S}_0^{-1} - \partial_{\bar{\kappa}_T} \mathcal{R}_0 \cdot \mathcal{S}_0^{-1}) \frac{\sigma_r \bar{\kappa}_T}{2(d-2)} - \tilde{\gamma}_0 \epsilon_0^{-1/d} \frac{d-2}{d-1} \mathcal{R}_0 \cdot \mathcal{S}_0^{-1} \right] (\omega, \bar{\kappa}_L, \bar{\kappa}_T) \\ & + (\text{transpose } \omega \to -\omega, \bar{\kappa}_L \to -\bar{\kappa}_L, \bar{\kappa}_T \to -\bar{\kappa}_T). \end{aligned} \quad (152) $$

An easier way to read off the first-order correction is to use (134) and (135) and consider similar $\bar{\sigma}^{-1}$ expansion of $\hat{G}_R$ and $\hat{G}_A$. Collecting first-order terms from (135), for instance, we obtain

$$ \begin{aligned} \hat{G}_{R,1} = \frac{(2\Delta_0 - d)}{2} \frac{d-2}{d-1} \left[ \frac{2A_1}{a_0} - \frac{2A_0a_1}{a_0^2} + \frac{A_0\sigma_r}{a_0^2} \left( \frac{\bar{\kappa}_T}{d-2} \partial_{\bar{\kappa}_T} - \bar{\kappa}_L \partial_{\bar{\kappa}_L} \right) a_0 \right. \\ \left. + \frac{\sigma_r}{a_0} \left( \frac{\bar{\kappa}_T}{d-2} \partial_{\bar{\kappa}_T} - \bar{\kappa}_L \partial_{\bar{\kappa}_L} \right) A_0 - \frac{\tilde{\gamma}_0 \epsilon_0^{-1/d} A_0}{a_0} \frac{d-2}{d-1} \right] (\omega, \bar{\kappa}_L, \bar{\kappa}_T), \end{aligned} \quad (153) $$

and

$$ \begin{aligned} \hat{G}_{A,1} = \frac{(2\Delta_0 - d)}{2} \frac{d-2}{d-1} \left[ \frac{2B_1}{b_0} - \frac{2B_0b_1}{b_0^2} + \frac{B_0\sigma_r}{b_0^2} \left( \frac{\bar{\kappa}_T}{d-2} \partial_{\bar{\kappa}_T} - \bar{\kappa}_L \partial_{\bar{\kappa}_L} \right) b_0 \right. \\ \left. + \frac{\sigma_r}{b_0} \left( \frac{\bar{\kappa}_T}{d-2} \partial_{\bar{\kappa}_T} - \bar{\kappa}_L \partial_{\bar{\kappa}_L} \right) B_0 - \frac{\tilde{\gamma}_0 \epsilon_0^{-1/d} B_0}{b_0} \frac{d-2}{d-1} \right] (\omega, \bar{\kappa}_L, \bar{\kappa}_T). \end{aligned} \quad (154) $$

This allows us to finally determine $\gamma_{1,in}$ and similarly $\gamma_{n,in}$ for $n > 1$ as follows. After integrating over the frequency and momentum integrals (i.e., doing integrals over $\omega$, $k_L$, and $\bar{k}_T$) as in (138), we obtain the first-order correction in the series expansion of the type (148):

$$ G_R(x_1, x_2) = \sum_{n=0}^{\infty} \frac{1}{\sigma^n \epsilon_0^{n/d}} G_{R,n}(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r}), \quad (155) $$

and similarly for $G_A(x_1, x_2)$. For the identities (116) and (117) to be satisfied, we therefore need (note that $\tilde{x}_{\perp r}$ is invariant under $x_1 \leftrightarrow x_2$)

$$ G_{R,n}(-\sigma_r, -\tilde{\zeta}_r, \tilde{x}_{\perp r}) = G_{A,n}(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r}) \quad (156) $$

at each $n$, which can be ensured by appropriate choices of $\gamma_{n,in}(\omega, \kappa_L, \bar{k}_T)$ at each order for $n \ge 1$, since they determine both $A_n$ and $a_n$. In turn, this justifies the bilocal thermal structure (134) and (135) which implies via (116) and (117) 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 $\gamma_{n,out}$ are only functions of the frequency $\omega$, $\gamma_{n,in}$ are functions of both frequency and momenta.

With $\gamma_{n,in}$ determined we know the RHS of (153) completely. Doing the frequency and momenta integrals shown in (131), we can obtain the first and similarly higher-order corrections to all Schwinger-Keldysh propagators.

## E. Consistency checks

The general result for the hidden bilocal thermal structure of the hydrodynamic correlation functions given by Eqs. (131)–(135) 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 ingoing solution which is analytic. The on-shell solution also vanishes. When $J_1 = J_2$, we have

$$ \langle O_1(x_1) \rangle = \int_{x_2} [G_{11}(x_1, x_2)J_1(x_2) - G_{12}(x_1, x_2)J_2(x_2)] = \int_{x_2} G_R(x_1, x_2)J_1(x_2) $$

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 bilocal thermal structure of the hydrodynamic correlation functions resulting in (138) as shown above. We can also see that the analytic ingoing mode is indeed always related to the causal response from the full nonlinear 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 $\Phi(v, \zeta, \tilde{x}_{\perp})$ at a constant Eddington-Finkelstein time hypersurface when the boundary source $J(\sigma, \zeta, \tilde{x}_{\perp})$ 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 Ref. [10] which is thus reproduced by the horizon cap prescription.<fnref n="24" />

A more nontrivial consistency check involves the outgoing mode. Consider an outgoing mode with frequency $\omega$ and momenta $k_L$ and $\vec{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(\omega, k_L, \vec{k}_T) = e^{\beta\omega} J_1(\omega, k_L, \vec{k}_T),$$

with $\beta$ defined by $\epsilon_0$ as in (97). Therefore, on the forward arm (where we choose to do the $k_L$ and $\vec{k}_T$ integrations first)

$$\begin{align*}
\langle O_1(x_1) \rangle &= \int d\omega \int d^d x_2 e^{-i\omega[s(\sigma_1)-s(\sigma_2)]} [\hat{G}_{11}(\omega, x_1, x_2)J_1(\omega, \vec{x}_2) - \hat{G}_{12}(\omega, x_1, x_2)J_2(\omega, \vec{x}_2)] \\
&= \int d\omega \int d^d x_2 e^{-i\omega[s(\sigma_1)-s(\sigma_2)]} [\hat{G}_{11}(\omega, x_1, x_2) - \hat{G}_{12}(\omega, x_1, x_2)e^{\beta\omega}]J_1(\omega, \vec{x}_2) \\
&= \int d\omega \int d^d x_2 e^{-i\omega[s(\sigma_1)-s(\sigma_2)]} [\hat{G}_{11}(\omega, x_1, x_2) - \hat{G}_{21}(\omega, x_1, x_2)]J_1(\omega, \vec{x}_2) \\
&= \int d\omega \int d^d x_2 e^{-i\omega[s(\sigma_1)-s(\sigma_2)]} [\hat{G}_A(\omega, x_1, x_2)J_1(\omega, \vec{x}_2) \\
&= \int d^d x_2 G_A(x_1, x_2)J_1(x_2).
\end{align*}$$

Above, we have used the result from (134) that $\hat{G}_{12}e^{\beta\omega} = \hat{G}_{21}$ and also the identity (117). Thus indeed we see that the outgoing mode gives advanced response. Note for both $J_2 = e^{\beta\omega}J_1$ and $\hat{G}_{12}e^{\beta\omega} = \hat{G}_{21}$ (used above) to hold, we need absence of log terms at the horizon, which follows if $\gamma_{n,out}$ are functions of $\omega$ only. The outgoing solution which satisfies the regularity condition (111) at the horizon cap indeed has this property with $\gamma_{n,out}$ determined from the cancellation of the double and single poles of the equation of motion being a linear function of $\omega$ and independent of $\kappa_L$ and $\vec{k}_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

$$\omega = \omega_Q,$$

where $\omega_Q$ corresponds to the homogeneous quasinormal mode frequencies of the static black brane, as noted by Janik and Peschanski. A nontrivial 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

$$\gamma_0 = \gamma_0(\omega_Q),$$

i.e., for AdS$_5$ we should have

$$\gamma_0(\omega) = \frac{\omega_Q}{4\epsilon_0^{1/4}}.$$

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 $\kappa_L$ and $\vec{k}_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 nonequilibrium event horizon is consistent with the feature in nonequilibrium 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 anticommutator 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 analog 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 initial conditions. For this to hold, the horizon cap should cover the space-time outside the event horizon, which is indeed the case.

## VI. ON INITIAL CONDITIONS AND SEEING BEHIND THE EVENT HORIZON

The discussion on initial conditions is subtle.<fnref n="25" /> 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

<footnote n="24">²⁴See Ref. [50] for the nonequilibrium retarded correlation function in a hydrodynamic expansion, and [51,52] for AdS-Vaidya and states corresponding to quenches, etc. The latter uses the prescription of Ref. [10].</footnote>

<footnote n="25">²⁵We thank Skenderis for discussions on these issues.</footnote>
$\Phi_1(v, \zeta, \vec{x}_{\perp})$ and $\Phi_2(v, \zeta, \vec{x}_{\perp})$ at an initial time $\sigma_0$ on both arms, and also the sources $J_1(\sigma, \zeta, \vec{x}_{\perp})$ and $J_2(\sigma, \zeta, \vec{x}_{\perp})$ at the two boundaries for all times in the future. Our horizon cap prescription leads to unique solutions for the ingoing and outgoing modes corresponding to $J_1(\sigma, \zeta, \vec{x}_{\perp})$ and $J_2(\sigma, \zeta, \vec{x}_{\perp})$ in a large proper time expansion. However, it is not guaranteed we can match with arbitrary initial conditions $\Phi_1(v, \zeta, \vec{x}_{\perp})$ and $\Phi_2(v, \zeta, \vec{x}_{\perp})$ at $\sigma = \sigma_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 $\Phi_1(v, \zeta, \vec{x}_{\perp})$ and $\Phi_2(v, \zeta, \vec{x}_{\perp})$, respectively. The full on-shell gravitational action beyond the hydrodynamic limit should yield the matrix element

$$ \langle \Phi_1 | T_c \left( \exp \left[ -i \oint J(x) \hat{O}(x) \right] \right) | \Phi_2 \rangle \quad (157) $$

with $T_c$ denoting time-ordering in the closed time contour, and $\langle \Phi_1|$ and $\langle \Phi_2|$ denoting states in the dual field theory corresponding to the semiclassical bulk field configurations.<fnref n="26" />

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. Although this needs to be investigated further, presently we may conclude that one may not be able to obtain semiclassical solutions corresponding to arbitrary in-in states meaning that decoherence, which suppresses some off-diagonal matrix elements, is built in the semiclassical gravity approximation. This issue is also relevant for the general approach of Skenderis and van Rees [18]. We also note that there is possibility of adding other semiclassical 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 (148) is divergent and would require a trans-series completion which would naively be of the form (like multi-instanton series)

$$ G(x_1, x_2) = \prod_{\alpha=1}^{\infty} \left\{ \sum_{n_{\alpha}=0}^{\infty} \left[ C_{\alpha}(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r}) \bar{\sigma}^{\gamma_{\alpha}(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r})} e^{-\xi_{\alpha}(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r})} \bar{\sigma} \right]^{n_{\alpha}} \right\} \Gamma_{n_1, n_2, \dots} (\bar{\sigma}, \sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r}), $$

with

$$ \Gamma_{n_1, n_2, \dots} (\bar{\sigma}, \sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r}) = \sum_{k=0}^{\infty} \Gamma_{k; n_1, n_2, \dots} (\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r}) \bar{\sigma}^{-k}. \quad (158) $$

This is a generalization of the trans-series for $\epsilon(\sigma)$ for the Bjorken flow which completes the divergent asymptotic hydrodynamic series at large proper time [22,45,55]. Above $\Gamma_{n_{\alpha}=0}$ coincides with the matrix hydrodynamic power series (148) which we have explicitly computed here; while generally $\Gamma_{n_1, n_2, \dots}$ is a similar matrix power series in $\bar{\sigma}^{-1}$ whose coefficients are related to those of (148) for each fixed $\sigma_r$, $\tilde{\zeta}_r$, and $\tilde{x}_{\perp r}$. Physically, the functions $\xi_{\alpha}(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r})$ and $\gamma_{\alpha}(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r})$ govern the hydrodynamic relaxation of the Green's function for fixed separations of the reparametrized space-time arguments $x_1$ and $x_2$. It is possible that these functions are approximately constants. Crucially the Stokes data $C_{\alpha}(\sigma_r, \tilde{\zeta}_r, \tilde{x}_{\perp r})$ are determined by initial conditions. Instead of being constants as in the case of $\epsilon(\tau)$, 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. 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 $\epsilon(\sigma)$]. 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.<fnref n="27" />

## VII. 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 nontrivial space-time dependence which can map the

<footnote n="26">²⁶In the approach of Skenderis and van Rees [18], the in-in states are defined using Euclidean path integrals. Even in this approach, it is unclear if one can glue arbitrary in-in states by bulk geometries. See Refs. [53,54] for some discussions on this issue.</footnote>

<footnote n="27">²⁷See Refs. [56,57] for recent reviews on current progress in resolving black-hole information paradoxes with substantial discussion on the role played by hair degrees of freedom at the horizon.</footnote>
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 timelike 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 space-time 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 [58].

Another important problem would be to understand the Borel resummation of the hydrodynamic series (148) and compute the Schwinger-Keldysh correlation functions of the holographic hydrodynamic attractor [22]. Furthermore, the late time thermal nature of the correlation functions in terms of the reparametrized space-time 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 [59].

Our prescription can be used to construct the quantum generalization of classical stochastic hydrodynamics (see Refs. [60,61] for reviews) in holographic theories by considering the backreaction of the fluctuations on the background geometry systematically. This is especially important in the 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 Ref. [14] for an excellent related discussion)—nonlinearities can potentially cause nontrivial effects such as quantum corrections to the long time tails.<fnref n="28" />

A more general goal is to develop a quantum generalization of the large deviation function [13] in classical nonequilibrium statistical mechanics, which is essentially a generalization of the equilibrium free energy to the hydrodynamic regime, assigning, for instance, probability to a macroscopic hydrodynamic configuration which may not satisfy the hydrodynamic equations. The approach of Ref. [33] which constructs an effective action at finite temperature from the complexified bulk space-time will be very relevant for such developments if it can be generalized to out-of-equilibrium situations at least in the hydrodynamic limit.

More generally, we would like to 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 [57,63] (which include the black-hole interior) from Hawking quanta. In such cases, semiholographic formulations for open quantum systems (see Refs. [64,65] and also [33], for instance) can provide useful models.

ACKNOWLEDGMENTS

We thank Kostas Skenderis and Balt van Rees for helpful discussions and sharing many insights on the issues discussed in Sec. VI. A.M. 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 the IFCPAR/CEFIPRA funded Project no. 6304-3. The research of A.B. is supported by IFCPAR/CEFIPRA funded Project no. 6304-3.

APPENDIX A: DETAILS OF THE FIVE-DIMENSIONAL METRIC DUAL TO THE BJORKEN FLOW
UP TO SECOND ORDER IN THE LATE-TIME EXPANSION

Consider the late-time expansion of Einstein's equations based on the ansatz (89) along with (92). At zeroth order, the six components of Einstein equations read

$$E_{vv} : 2k_0'(x)^2 + 4k_0''(x) + l_0'(x)^2 + 2l_0''(x) = 0,$$

$$E_{v\sigma} : a_0(x)x^2[4k_0''(x) + l_0'(x)^2 + 2l_0''(x)] + 2xa_0(x)k_0'(x)[xl_0'(x) - 6] + 3x^2a_0(x)k_0'(x)^2 - 6xa_0(x)l_0'(x) + 24a_0(x) + xa_0'(x)[2xk_0'(x) + xl_0'(x) - 6] - 24 = 0,$$

$$E_{\sigma\sigma} : a_0(x)x^2[4k_0''(x) + l_0'(x)^2 + 2l_0''(x)] + 2xa_0(x)k_0'(x)[xl_0'(x) - 6] + 3x^2a_0(x)k_0'(x)^2 - 6xa_0(x)l_0'(x) + 24a_0(x) + xa_0'(x)[2xk_0'(x) + xl_0'(x) - 6] - 24 = 0,$$

$$E_{s_{\perp}s_{\perp}} : 2x(a_0'(x)[xk_0'(x) + l_0'(x)] - 6) + xa_0''(x)) + 24a_0(x) - 24 + a_0(x)x\{k_0'(x)[xl_0'(x) - 6] + 2x[k_0''(x) + l_0''(x)] + xa_0'(x)^2 + xl_0'(x)^2 - 6l_0'(x)\} = 0,$$

$$E_{\xi\zeta} : a_0(x)[3x^2k_0'(x)^2 + 4x^2k_0''(x) - 12xk_0'(x) + 24] + 4x[xk_0'(x) - 3]a_0'(x) + 2x^2a_0''(x) - 24 = 0, \quad (A1)$$

<footnote n="28">²⁸For a recent work on the role of classical fluctuations of the superfluid order parameter in modifying transport properties, see Ref. [62].</footnote>
where $x \equiv \epsilon_0^{1/d} v$. We solve these equations based on the conditions

$$a_0(x) = 1 - x^4, \quad l_0(x) = 0, \quad k_0(x) = 0,$$

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:

$$2k_i''(x) + l_i''(x) = S_{1,i}, \quad (x^4 + 3)[k_i'(x) - l_i'(x)] + x(x^4 - 1)[k_i''(x) - l_i''(x)] = S_{2,i}, \\ \frac{3}{2}x^2 a_i''(x) - 6xa_i'(x) + 6a_i(x) - 2x^5[2k_i'(x) + l_i'(x)] = S_{3,i}, \quad (A2)$$

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

$$2k_i(x) + l_i(x) = \rho_i + \alpha_i x + PI_{1,i}, \quad k_i(x) - l_i(x) = \xi_i - \frac{\kappa_i}{4} \ln(1 - x^4) + PI_{2,i}, \\ a_i(x) = \beta_i x^4 + \alpha_i \frac{x(3 + x^4)}{3} + PI_{3,i}. \quad (A3)$$

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 $PI$ for $i = 1$ are simple and already mentioned in (93) and (94). For $i = 2$, we get

$$PI_{1,2} = \frac{1}{12}[8x^2 + 6 \ln(x^2 + 1) - \alpha_1^2(x + 2)x - 4\alpha_1(x + 2)x + 16x + 4 \ln(x + 1) + 2(2 - 6x) \arctan(x) + 4],$$

$$PI_{2,2} = \text{Re} \left\{ \frac{1}{24} \left[ -24x^2 - \frac{16x}{x^3 + x^2 + x + 1} + 16x + (4 + 4i) \arctan(x)^2 - (1 - i) \ln^2(i - x) \right. \right. \\ + 4 \ln^2(x + 1) - 2 \ln^2(-i + x) - (3 + i) \ln^2(i + x) - 4 \ln^2(x^2 + 1) + (2 + 8i) \arctan \left( \frac{x+1}{1-x} \right) \\ + 2i\pi \arctan(x) - 8 \arctan(x) + (2 - 8i) \arctan \left( \frac{x-1}{x+1} \right) - \pi \ln(2) + 4\pi \ln(1 + e^{-2i \arctan(x)}) \\ - 16 \arctan(x) \ln(1 + e^{2i \arctan(x)}) + 8 \arctan(x) \ln(1 - ie^{2i \arctan(x)}) + 2\pi \ln(1 - ie^{2i \arctan(x)}) \\ + 6 \ln(1-x) - (4 + 4i) \ln \left[ \left( \frac{1}{2} + \frac{i}{2} \right) (ix + 1) \right] \ln(x+1) \\ + 2 \ln(x+1) + 4 \ln(x-1) \ln \left[ \left( \frac{1}{2} + \frac{i}{2} \right) (-i + x) \right] - (6 + 2i) \ln \left( \frac{1}{2} (ix + 1) \right) \ln(i+x) \\ + 4 \ln(x-1) \ln \left[ \left( \frac{1}{2} - \frac{i}{2} \right) (i+x) \right] + (4 + 4i) \ln(x+1) \ln \left[ \left( -\frac{1}{2} - \frac{i}{2} \right) (i+x) \right] \\ - (2 - 2i) \ln(i-x) \ln \left( -\frac{1}{2}i(i+x) \right) - 4 \ln(-i+x) \ln \left( -\frac{1}{2}i(i+x) \right) + (2 - 2i) \ln(i-x) \ln(x^2+1) \\ - 4 \ln(x-1) \ln(x^2+1) + 4 \ln(-i+x) \ln(x^2+1) + (6 + 2i) \ln(i+x) \ln(x^2+1) \\ + 2\pi \ln(x^2+1) - 4 \ln(x^2+1) + 8 \ln(x+1) \ln \left( 1 - \left( \frac{1}{2} - \frac{i}{2} \right) (x+1) \right) \\ - 4\pi \ln \left( \sin \left( \arctan(x) + \frac{\pi}{4} \right) \right) + 2\pi \ln \left( \frac{x+1}{\sqrt{x^2+1}} \right) - (6 - 2i) \text{Li}_2 \left( \frac{ix}{2} + \frac{1}{2} \right) + 8i \text{Li}_2(-e^{2i \arctan(x)}) \\ - 4i \text{Li}_2(ie^{2i \arctan(x)}) + 4 \text{Li}_2\left(\left(-\frac{1}{2} + \frac{i}{2}\right)(x-1)\right) + 4 \text{Li}_2\left(\left(-\frac{1}{2} - \frac{i}{2}\right)(x-1)\right) + (4 + 4i) \text{Li}_2\left(\left(\frac{1}{2} + \frac{i}{2}\right)(x+1)\right) \\ + (4 - 4i) \text{Li}_2\left(\left(\frac{1}{2} - \frac{i}{2}\right)(x+1)\right) - (6 + 2i) \text{Li}_2\left(-\frac{1}{2}i(i+x)\right) + \left[ -\frac{8x}{x^3 + x^2 + x + 1} + 8x \right. \\ \left. + (2 + 4i) \arctan \left( \frac{x+1}{1-x} \right) - 4 \arctan(x) + (2 - 4i) \arctan \left( \frac{x-1}{x+1} \right) \right] \alpha_1 - \frac{8}{x^3 + x^2 + x + 1} \right\}$$
$$PI_{3,2} = \frac{x}{108} \{ \alpha_1^2(-9x^5 - 6x^4 + 3x + 6) - 4\alpha_1(9x^5 + 12x^4 - 9x - 8) + 4[-9x^5 - 9(x^4 + 1) \arctan(x) + x^3 \\ + 9x^3 \ln(x^2 + 1) + 9x + 16] \},$$

where $Li_2(x)$ is the polylogarithmic function of order two. 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:

(i) The fixation of $\alpha_i$ has been detailed in (112)-(114). The same pattern repeats for the rest of $\alpha_{i>2}$.

(ii) The integration constants $\beta_i$ are fixed by mapping the coefficients of $x^4$ terms in (A3) to $a_4(\tau)$, $l_4(\tau)$, and $k_4(\tau)$ of (79) [or $a_4(\sigma)$, $l_4(\sigma)$, and $k_4(\sigma)$ in the Weyl transformed coordinate] and then solving the constraint equations (80).

(iii) The constants $\rho_i$ and $\xi_i$ are fixed by demanding that metric asymptotes to (58) at the boundary.

(iv) Finally, the integration constants $\kappa_i$ are fixed such that the dual geometry is regular at the horizon to all orders.

For $i = 1, 2$, this gives

$$\beta_1 = \frac{2}{3}, \quad \kappa_1 = 2, \quad \xi_1 = 0, \quad \rho_1 = 0,$$

$$\beta_2 = -\frac{1}{54}[11 + 6 \ln(2)], \quad \kappa_2 = -\frac{1}{3}[2 \ln(2) - 3], \quad \rho_2 = -\frac{1}{3},$$

$$\xi_2 = \frac{1}{288} i \left[ -24\pi\alpha_1 + 48iC - (48 - 96i)\mathrm{Li}_2\left(\frac{1}{2} + \frac{i}{2}\right) + (48 + 96i)\mathrm{Li}_2\left(\frac{1}{2} - \frac{i}{2}\right) + (7 + 6i)\pi^2 \right. \\ \left. - 96i - 48\pi + 72i \ln^2(2) - (24 - 60i)\pi \ln(2) + 48\pi \ln(2) + 24i\pi \ln(1-i) \right], \quad (A4)$$

where $C$ is the Catalan's constant $C = 0.915966$.

One can verify that these solutions reproduce the same late-time expansion of the stress tensor (setting $\tau_0 = 1$ and $\varepsilon_0 = 1$)

$$T_{\tau\tau} \approx \tau^{-4/3} - \frac{2}{3}\tau^{-2} + \frac{1 + 2 \ln(2)}{18}\tau^{-8/3} + \dots,$$

$$T_{\xi\xi} \approx \frac{1}{3}\tau^{2/3} - \frac{2}{3} + \frac{5[1 + 2 \ln(2)]}{54}\tau^{-2/3} + \dots, \quad (A5)$$

$$T_{x_{\perp}x_{\perp}} \approx \frac{1}{3}\tau^{-4/3} - \frac{[1 + 2 \ln(2)]}{54}\tau^{-8/3} + \dots.$$

## APPENDIX B: THE EVENT HORIZON AND THE APPARENT HORIZON OF THE BJORKEN FLOW

This Appendix is devoted to a comparative discussion between the event horizon and the apparent horizon in the dynamical geometry dual to the Bjorken flow.

### 1. Event horizon

Classically, the event horizon marks the null hypersurface from which no signal can come out, i.e., the outgoing null rays become tangential to the hypersurface. To determine its location in the geometry given by (89), consider radial null geodesics in this space-time whose equation is given by

$$\frac{\partial v}{\partial \sigma} + \frac{v^2}{2} \frac{d-2}{d-1} g_{\sigma\sigma} = 0 \quad \text{with} \quad g_{\sigma\sigma} = \frac{1}{v^2} \left( \frac{(d-1)^2}{(d-2)^2} A(v, \sigma) + \frac{2(d-1)v}{(d-2)^2 \sigma} \right). \quad (B1)$$

For a static geometry, $v_E = \text{const}$ and the event horizon is simply given by the zero of $g_{\sigma\sigma}$. However, in a dynamical geometry, the event horizon will depend on $\sigma$. Consider the late-time expansion of the horizon of the form

$$v(\sigma) \equiv v_E(\sigma) = \sum_{i=0} v_{Ei} \sigma^{-i}. \quad (B2)$$

Now to determine $v_{Ei}$, we use the equation (B1) along with the known late-time expansion of $A(v, \sigma)$,

$$A(v, \sigma) = 1 - \varepsilon_0 v^d + \frac{1}{\varepsilon_0^{1/d} \sigma} a_1(v) + \frac{1}{\varepsilon_0^{2/d} \sigma^2} a_2(v) + \dots. \quad (B3)$$

At leading order, we have

$$v_{E0} = 1/\varepsilon_0^{1/d}$$
in any space-time dimensions $d+1$. The subleading corrections to the event horizon depend on the residual gauge parameters. For instance, the first and second subleading corrections in $d=4$ are

$$v_{E1} = \frac{(3 + \alpha_1)}{6\epsilon_0^{1/2}}, \quad v_{E2} = \frac{72\alpha_2 - 9\pi - 20 + 12 \ln(2)}{216\epsilon_0^{3/4}}, \qquad (B4)$$

where $\alpha_1$ and $\alpha_2$ are the residual gauge parameters associated with the late-time expansion of $A(v, \sigma)$. Now the gauge parameters can be uniquely fixed by demanding regularity of the horizon cap. In $d=4$, this gives [recall (113) and (114)]

$$\alpha_1 = -3, \quad \alpha_2 = \frac{9\pi + 20 - 12 \ln(2)}{72}. \qquad (B5)$$

These values, in turn fix the event horizon to $1/\epsilon_0^{1/d}$ up to the second subleading order. Thus, regularization of the horizon cap leads to vanishing of the subleading 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.

## 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 $\theta_{\pm}$, i.e., $\Theta = e^f \theta_+ \theta_-$, where the factor $e^f$ is defined below. The trapped region is characterized by $\Theta > 0$, where the light rays directed outward propagates inward, whereas in the untrapped region the light rays directed outward propagate outward and are characterized by $\Theta < 0$. Therefore, $\Theta = 0$ gives the location of the apparent horizon.

Here we adopt the dual-null formalism [43,66–68] to study location of the apparent horizon, where one defines a pair of null hypersurfaces $\Sigma^{\pm}$ parameterized by scalars $\zeta^{\pm}$, with associated one-forms $n^{\pm} = -d\zeta^{\pm}$. The null normal vector associated with these hypersurfaces are given by

$$l_{\pm}^{\mu} = e^{-f} g^{\mu\nu} n_{\nu}^{\mp},$$

where $e^f$ is the normalization factor given by

$$e^f = -g^{\mu\nu} n_{\mu} n_{\nu}^{\mp}.$$

Next, we define the expansion parameters $\theta_{\pm}$,

$$\theta_{\pm} = \mathcal{L}_{\pm} \ln(\mu),$$

where $\mu$ is the spatial volume element of the geometry in which the hypersurfaces are defined and $\mathcal{L}_{\pm}$ is the Lie derivative along the null normal vectors, $l_{\pm}$. Finally, we introduce the invariant$^{30}$ quantity $\Theta = e^f \theta_+ \theta_-$, whose zero gives the location of the apparent horizon.

In case of our geometry (89), 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

$$n^{-} = -N^{-} d\sigma, \quad n^{+} = -N^{+} \left\{ \left[ \frac{2v}{\sigma} + (d-1)A(v, \sigma) \right] d\sigma + 2(d-2)dv \right\}, \qquad (B6)$$

and the null normal vectors $l_{\pm}^{\mu}$ along with the normalization factor $e^f$ and the volume element $\mu$ reads

$$\begin{align}
e^f &= \frac{2(d-2)^2}{d-1} N^+ N^- v^2, \notag \\
l_+^\mu &= \left( \frac{1}{2(d-2)N^+}, 0, 0, \vec{0} \right), \notag \\
l_-^\mu &= \left( -\frac{\sigma(d-1)A(v, \sigma) + 2v}{2(d-2)N^- \sigma}, 1/N^-, 0, \vec{0} \right), \notag \\
\mu &= \frac{\sigma + v}{v^3 \sigma} \sqrt{e^{2K(v, \sigma) + L(v, \sigma)}}. \tag{B7}
\end{align}$$

where $N^+$ and $N^-$ are the overall normalizability factor that can be determined by the integrability condition $d(dn^{\pm}) = 0$. This ensures that the one-forms (B6) and (B7) are exact. However there is no need of computing them explicitly, as when we compute $\Theta$ the contribution of $N^+$ and $N^-$ from $e^f$ will cancel with those coming from $\theta_-$ and $\theta_+$.

<footnote>$^{29}$For $d=3$, the first subleading correction is $v_{E1} = 5(5 + 3\alpha_1)/3\epsilon_0^{2/3}$, which again depends on the residual parameter. However, these corrections do not have a universal form like the leading order.</footnote>

<footnote>$^{30}$This parameter $\Theta$ is invariant under reparameterization of the scalar $\zeta^{\pm} \rightarrow \xi^{\pm}(\zeta^{\pm})$ or interchange of $\zeta^{\pm} \rightarrow \zeta^{\mp}$, as discussed in Ref. [66].</footnote>
<figure><img src="image_6.png" /><figcaption>FIG. 4. The figure shows the evolution of the event horizon (blue curve) and the apparent horizon (yellow curve) in the dynamical geometry dual to the Bjorken flow, up to second order in late-time expansion. Here $\epsilon_0$ is set to one and the gauge parameters $\alpha_1$ and $\alpha_2$ are so chosen that the event horizon is fixed at one. The apparent horizon always lies inside the event horizon (at greater $v$).</figcaption></figure>

Now since the apparent horizon depends on $\sigma$, one can consider a late-time expansion of the apparent horizon similar to the event horizon as

$$v_A(\sigma) = \sum_{i=0} v_{Ai} \sigma^{-i}.$$

To determine the coefficients $v_{Ai}$ we solve $\Theta = 0$ at every order. In case of any space-time dimensions $d+1$, the leading-order result is universal and turns out to be

$$v_{A0} = 1/\epsilon_0^{1/d},$$

which coincides with the event horizon. However, the subleading corrections are dimension dependent. For $d=4$, the first and second subleading corrections to the leading terms are

$$ \begin{aligned} v_{A1} &= \frac{3 + \alpha_1}{6\epsilon_0^{1/2}}, \\ v_{A2} &= \frac{72\alpha_2 - 9\pi - 8 + 12 \ln(2)}{216\epsilon_0^{3/4}}. \end{aligned} \qquad (B8) $$

Note that, for $\alpha_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 $\alpha_2$ [given by (114)] for which the correction to the event horizon vanishes, now renders $v_{A2} > 0$. So the apparent horizon will lie inside the event horizon (see Fig. 4). Again, we expect this feature to remain true at higher orders as well.

## APPENDIX C: MORE DETAILS OF THE BULK SCALAR FIELD IN THE AdS$_5$ BJORKEN FLOW BACKGROUND

In this Appendix, we provide brief details of the solution of Klein-Gordon equation at the leading (102) and first subleading order (107) (for $n=1$) and validate the near-horizon behaviors (106) and (110) for the same. To be specific, we will consider the example of $d=4$.

At leading order, the homogenous Klein-Gordon equation (102) 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:

$$ \phi_{0,\text{in}}(v, \omega, \kappa_L, \vec{\kappa}_T) = p_0 \left\{ 1 + (\epsilon_0^{-1/4} - v) \frac{(\kappa_L^2 + \vec{\kappa}_T^2 - m^2 \epsilon_0^{1/2} + 3i\epsilon_0^{-1/4}\omega)}{2(2\epsilon_0^{1/4} - i\omega)} + O\left[(\epsilon_0^{-1/4} - v)^2\right] + \dots \right\}, \quad (C1) $$

$$ (\epsilon_0^{-1/4} - v)^{-\frac{i\omega}{2\epsilon_0^{1/4}}} \phi_{0,\text{out}}(v, \omega, \kappa_L, \vec{\kappa}_T) = q_0 \left\{ 1 + (\epsilon_0^{-1/4} - v) \frac{(2\kappa_L^2 + 2\vec{\kappa}_T^2 - 2m^2\epsilon_0^{1/2} - 3i\omega)}{4(2\epsilon_0^{1/4} + i\omega)} + O\left[(\epsilon_0^{-1/4} - v)^2\right] + \dots \right\}. \quad (C2) $$

Now, without loss of generality, we can set $p_0 = 1$ and $q_0 = 1$. Then, in the limit $v \to \epsilon_0^{-1/4}$ we reproduce the behavior given in (106).

At the first subleading order ($n=1$), the Klein-Gordon equation is given by (107) with

$$ S_{1,\text{in}} = \tilde{S}\phi_{0,\text{in}}(v, \omega, \kappa_L, \vec{\kappa}_T), \quad S_{1,\text{out}} = \tilde{S}\phi_{0,\text{out}}(v, \omega, \kappa_L, \vec{\kappa}_T), \qquad (C3) $$

where

$$ \begin{aligned} \tilde{S} = & -\frac{v}{6\epsilon_0^{1/4}} \{ \alpha_1 \epsilon_0^{1/4} v (2\kappa_L^2 v + 2\vec{\kappa}_T^2 v + 3i\omega) + 2v(2\kappa_L^2 - \vec{\kappa}_T^2) \tan^{-1}(\epsilon_0^{1/4} v) + 4(\kappa_L^2 + \vec{\kappa}_T^2) \epsilon_0^{1/4} v^2 \\ & + (2\kappa_L^2 - \vec{\kappa}_T^2) v \left[ \ln(\epsilon_0^{1/2} v^2 + 1) + 2 \ln(\epsilon_0^{1/4} v + 1) \right] + 12i\gamma_0 \epsilon_0^{1/4} + 6i\epsilon_0^{1/4} v\omega \} \\ & -\frac{1}{6} v \left[ (4\epsilon_0^{3/4} v^5 + 2\alpha_1 \epsilon_0 v^6 + 4\epsilon_0 v^6 + 2\alpha_1 v^2 + 4v^2) \partial_v^2 + (2\epsilon_0 v^5 - \alpha_1 v - 8i\gamma_0 v + 4\epsilon_0^{3/4} v^4 + \alpha_1 \epsilon_0 v^5 - 2v) \partial_v \right. \\ & \left. - 12\kappa_L \partial_{\kappa_L} + 6\vec{\kappa}_T \partial_{\vec{\kappa}_T} - 4\vec{\kappa}_T v \partial_{\vec{\kappa}_T} \partial_v + 8\kappa_L v \partial_{\kappa_L} \partial_v \right]. \end{aligned} \qquad (C4) $$

The action of $\tilde{S}$ on $\phi_{0,\text{out}}$ gives rise to singular terms having poles of order one and two at the horizon, as mentioned in (112). Upon removing these poles by suitably fixing $\alpha_1$ and $\gamma_0$ as in (113), the particular solution at the first subleading order for the
outgoing mode admits the near-horizon expansion,

$$ (\epsilon_0^{-1/4} - v)^{-\frac{i\omega}{2\epsilon_0^{-1/4}}} \phi_{1(p),\text{out}}(\omega, v, \kappa_L, \vec{k}_T) = \frac{(\epsilon_0^{-1/4} - v)}{48\epsilon_0^{1/4} + 24i\omega} \left\{ -2\kappa_L^2[\pi - 7 + \ln(64)] + \vec{k}_T^2[\pi + 2 + \ln(64)] \right. \\ \left. - 2m^2\epsilon_0^{1/2} - 4i\epsilon_0^{1/4}\omega - \omega^2 \right\} + O((\epsilon_0^{-1/4} - v)^2) + \dots \quad (C5) $$

However, for the ingoing mode no such divergent terms appear and we simply have the following near-horizon expansion:

$$ \phi_{1(p),\text{in}}(\omega, v, \kappa_L, \vec{k}_T) = (\kappa_L^2 \{-4\epsilon_0^{1/4}[\pi - 7 + \ln(64)] + 2i\omega[\pi - 3 + \ln(64)]\} \\ + \vec{k}_T^2 \{2\epsilon_0^{1/4}[\pi + 2 + \ln(64)] - i\omega[\pi + 6 + \ln(64)]\} \\ - 4m^2\epsilon_0^{3/4} + 2im^2\epsilon_0^{1/2}\omega + 6\epsilon_0^{1/4}\omega^2 + 12i\epsilon_0^{1/2}\omega) \frac{(\epsilon_0^{-1/4} - v)}{24(2\epsilon_0^{1/4} - i\omega)^2} \\ + O((\epsilon_0^{-1/4} - v)^2) + \dots \quad (C6) $$

Clearly, in the limit $v \to \epsilon_0^{-1/4}$, we reproduce (110) for $d = 4$ and $n = 1$. Similar behavior persists at higher orders as well.

## APPENDIX D: DERIVATION OF EQS. (131)

The source and vev of the dual operator at the boundary are given by

$$ J(\sigma_2, \hat{\xi}_2, \vec{x}_{\perp 2}) = \int d\omega dk_L d^{d-2}k_T e^{-i\omega s(\sigma_2)} e^{i(k_L \hat{\xi}_2 + \vec{k}_T \cdot \vec{x}_{\perp 2})} S(\sigma_2, \omega, \kappa_L, \vec{k}_T) \cdot Q(\omega, k_L, k_T), \quad (D1) $$

$$ \langle O(\sigma_1, \hat{\xi}_1, \vec{x}_{\perp 1}) \rangle = (2\Delta_O - d) \int d\omega dk_L d^{d-2}k_T e^{-i\omega s(\sigma_1)} e^{i(k_L \hat{\xi}_1 + \vec{k}_T \cdot \vec{x}_{\perp 1})} \mathcal{R}(\sigma_1, \omega, \kappa_L, \vec{k}_T) \cdot Q(\omega, k_L, k_T), \quad (D2) $$

where

$$ J(\sigma, \hat{\xi}, \vec{x}_{\perp}) = \begin{pmatrix} J_1(\sigma, \hat{\xi}, \vec{x}_{\perp}) \\ J_2(\sigma, \hat{\xi}, \vec{x}_{\perp}) \end{pmatrix}, \quad \langle O(\sigma, \hat{\xi}, \vec{x}_{\perp}) \rangle = \begin{pmatrix} \langle O_1(\sigma, \hat{\xi}, \vec{x}_{\perp}) \rangle \\ \langle O(\sigma, \hat{\xi}, \vec{x}_{\perp}) \rangle \end{pmatrix}, \quad Q(\omega, k_L, k_T) = \begin{pmatrix} p(\omega, k_L, k_T) \\ q(\omega, k_L, k_T) \end{pmatrix}. $$

To extract the Sk Green's function, we first invert (D1) to get

$$ Q(\omega, k_L, k_T) = \int d\sigma_2 d\hat{\xi}_2 d^{d-2}x_{\perp 2} s'(\sigma_2) e^{i\omega s(\sigma_2)} e^{-i(k_L \hat{\xi}_2 + \vec{k}_T \cdot \vec{x}_{\perp 2})} S^{-1}(\sigma_2, \omega, \kappa_L, \vec{k}_T) \cdot J(\sigma_2, \hat{\xi}_2, \vec{x}_{\perp 2}). \quad (D3) $$

Finally, plugging it back to (D2) and using linear-response theory,

$$ \langle O(\sigma_1, \hat{\xi}_1, \vec{x}_{\perp 1}) \rangle = \int d\sigma_2 d\hat{\xi}_2 d^{d-2}\vec{x}_{\perp 2} \tilde{G}(\sigma_1, \sigma_2, \hat{\xi}_1, \hat{\xi}_2, \vec{x}_{\perp 1}, \vec{x}_{\perp 2}) J(\sigma_2, \hat{\xi}_2, \vec{x}_{\perp 2}), $$

we read off the SK Green's function

$$ \tilde{G}(\sigma_1, \sigma_2, \hat{\xi}_1 - \hat{\xi}_2, \vec{x}_{\perp 1} - \vec{x}_{\perp 2}) = \int d\omega dk_L d^{d-2}k_T e^{-i\omega(s(\sigma_1) - s(\sigma_2))} e^{ik_L(\hat{\xi}_1 - \hat{\xi}_2) + i\vec{k}_T \cdot (\vec{x}_{\perp 1} - \vec{x}_{\perp 2})} \hat{G}(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T), \quad (D4) $$

where

$$ \hat{G}(\sigma_1, \sigma_2, \omega, k_L, \vec{k}_T) = \frac{1}{2}(s'(\sigma_2)\sigma_3 \cdot \mathcal{R}(\sigma_1, \omega, \kappa_{L1}, \vec{k}_{T1}) \cdot S^{-1}(\sigma_2, \omega, \kappa_{L2}, \vec{k}_{T2}) \\ + (\text{transpose } \sigma_1 \leftrightarrow \sigma_2, \omega \to -\omega, \kappa_{L1} \leftrightarrow -\kappa_{L2}, \vec{k}_{T1} \leftrightarrow -\vec{k}_{T2})), \quad (D5) $$

with

$$ \kappa_{L1} = k_L \frac{\tau_0}{\sigma_1}, \quad \kappa_{L2} = k_L \frac{\tau_0}{\sigma_2}, \quad \vec{k}_{T1} = \vec{k}_T \left( \frac{\tau_0}{\sigma_1} \right)^{-1/(d-2)}, \quad \vec{k}_{T2} = \vec{k}_T \left( \frac{\tau_0}{\sigma_2} \right)^{-1/(d-2)}. $$

First, note that the variable conjugate to the momentum $\omega$ is the reparametrized time $s(\sigma)$ whereas the boundary quantities depend on $\sigma$. This gives rise to the factor $s'(\sigma_2)$ in (D3) as the Jacobian of the transformation $s(\sigma_2) \to \sigma_2$. Also note that, due to translational invariance along the spatial directions, the matrices $S$ and $\mathcal{R}$ are independent of $\hat{\xi}$ and $\vec{x}_{\perp}$. Hence, on carrying out the momentum integrals in (D4), we would get $\tilde{G}(\sigma_1, \sigma_2, \hat{\xi}_1, \hat{\xi}_2, \vec{x}_{\perp 1}, \vec{x}_{\perp 2}) \equiv \tilde{G}(\sigma_1, \sigma_2, \hat{\xi}_1 - \hat{\xi}_2, \vec{x}_{\perp 1} - \vec{x}_{\perp 2})$.
APPENDIX E: TRANSIENTS AND $\gamma_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 (104) with only the ingoing modes of a massless scalar, i.e.,

$$ \Phi(v, \sigma) = \int d\omega e^{-i\frac{d-1}{2}\omega\sigma} \left(\frac{\sigma}{\tau_0}\right)^{i\gamma_0(\omega/\epsilon_0^{1/d})} \sum_{n=0}^{\infty} (\epsilon_0^{1/d}\sigma)^{-n} \phi_{n,\text{in}}(v, \omega)p(\omega), \quad (E1) $$

where $\phi_{n,\text{in}}(v, \omega)$ obeys the same equations (102) and (107) [and hence admit the solutions (106) and (109)] with $\kappa_L, \bar{\kappa}_T = 0$. The residual gauge parameters $\alpha_i$ appearing in the subleading solutions are fixed such that at every order, the event horizon is pinned at $v_h = \epsilon_0^{1/4}$ with associated Hawking temperature $T = \epsilon_0^{1/4}/\pi$. For concreteness, we again consider $d = 4$.

Recall that, at leading order, the ingoing mode (with zero spatial momentum) admits the near-horizon expansion

$$ \phi_{0,\text{in}}(v, \omega) = 1 + \frac{3i\epsilon_0^{1/4}\omega}{2(2\epsilon_0^{1/4} - i\omega)}(\epsilon_0^{-1/4} - v) + O((\epsilon_0^{-1/4} - v)^2) + \dots \quad (E2) $$

To compute the transients, we introduce the dimensionless decay rate $\lambda = \omega/\pi T$ and also scale the horizon to set $v_h = 1$. In terms of the dimensionless parameter, the solution (E2) when expanded near the boundary $v = 0$ reads

$$ \phi_{0,\text{in}}(v = 0, \lambda) \equiv a_0(\lambda) = 1 - \frac{3\lambda}{2(2i + \lambda)} + \dots \quad (E3) $$

The quasinormal mode corresponds to $a_0(\lambda_{qnm}) = 0$, implying that there is no contribution to the source from the leading order field. The lowest quasinormal mode frequency obtained in this method turns out to be

$$ \lambda_{qnm} = -2.7668i \pm 3.11945, $$

which agrees with the one obtained from spectral representation method up to order $10^{-10}$.

Similarly at subleading orders $n > 0$, the sources can be made to vanish, i.e., $a_{n>0}(\lambda_{qnm}) = 0$ by suitably fixing $\gamma_0$ (for $n = 1$) and $\gamma_{n,\text{in}}$ (for $n > 1$) as functions of $\lambda_{qnm}$. For example, corresponding to this lowest $\lambda_{qnm}$, $\gamma_0(\lambda_{qnm})$ obtained by setting $a_1(\lambda_{qnm}) = 0$ at the first subleading order is

$$ \gamma_0(\lambda_{qnm}) = -0.68669i \pm 0.779863, \quad (E4) $$

which turns out to be same as the one computed from the outgoing mode analysis (113),

$$ \gamma_0 = \frac{\omega_Q}{4\epsilon_0^{1/4}} = \frac{\lambda_{qnm}}{4} = -0.68669i \pm 0.779863. \quad (E5) $$

These two results agree up to order $10^{-9}$.

See Ref. [69] for a recent relevant work in the large-$D$ limit.

[1] J. Berges, Introduction to nonequilibrium quantum field theory, AIP Conf. Proc. 739, 3 (2004).

[2] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh field theory for driven open quantum systems, Rep. Prog. Phys. 79, 096001 (2016).

[3] M. Heyl, Dynamical quantum phase transitions: A review, Rep. Prog. Phys. 81, 054001 (2018).

[4] T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, Thermalization and prethermalization in isolated quantum systems: A theoretical overview, J. Phys. B: At., Mol. Opt. Phys. 51, 112001 (2018).

[5] G. Policastro, D. T. Son, and A. O. Starinets, The Shear viscosity of strongly coupled $N = 4$ supersymmetric Yang-Mills plasma, Phys. Rev. Lett. 87, 081601 (2001).

[6] S. Bhattacharyya, V. E. Hubeny, S. Minwalla, and M. Rangamani, Nonlinear fluid dynamics from gravity, J. High Energy Phys. 02 (2008) 045.

[7] R. Baier, P. Romatschke, D. T. Son, A. O. Starinets, and M. A. Stephanov, Relativistic viscous hydrodynamics, conformal invariance, and holography, J. High Energy Phys. 04 (2008) 100.

[8] M. Rangamani, Gravity and hydrodynamics: Lectures on the fluid-gravity correspondence, Classical Quantum Gravity 26, 224003 (2009).

[9] P. M. Chesler and L. G. Yaffe, Numerical solution of gravitational dynamics in asymptotically anti-de Sitter spacetimes, J. High Energy Phys. 07 (2014) 086.

[10] S. Banerjee, T. Ishii, L. K. Joshi, A. Mukhopadhyay, and P. Ramadevi, Time-dependence of the holographic spectral function: Diverse routes to thermalisation, J. High Energy Phys. 08 (2016) 048.

[11] V. Balasubramanian, A. Bernamonti, J. de Boer, N. Copland, B. Craps, E. Keski-Vakkuri, B. Muller, A. Schafer, M. Shigemori, and W. Staessens, Holographic thermalization, Phys. Rev. D 84, 026010 (2011).
[12] M. Le Bellac, *Thermal Field Theory*, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 1996).

[13] H. Touchette, The large deviation approach to statistical mechanics, Phys. Rep. **478**, 1 (2009).

[14] D. Bernard, Can the macroscopic fluctuation theory be quantized? J. Phys. A: Math. Theor. **54**, 433001 (2021).

[15] D. T. Son and A. O. Starinets, Minkowski space correlators in AdS/CFT correspondence: Recipe and applications, J. High Energy Phys. 09 (2002) 042.

[16] C. P. Herzog and D. T. Son, Schwinger-Keldysh propagators from AdS/CFT correspondence, J. High Energy Phys. 03 (2003) 046.

[17] K. Skenderis and B. C. van Rees, Real-time gauge/gravity duality, Phys. Rev. Lett. **101**, 081601 (2008).

[18] K. Skenderis and B. C. van Rees, Real-time gauge/gravity duality: Prescription, renormalization and examples, J. High Energy Phys. 05 (2009) 085.

[19] P. Glorioso, M. Crossley, and H. Liu, A prescription for holographic Schwinger-Keldysh contour in non-equilibrium systems, arXiv:1812.08785.

[20] M. Henningson and K. Skenderis, The holographic Weyl anomaly, J. High Energy Phys. 07 (1998) 023.

[21] S. de Haro, S. N. Solodukhin, and K. Skenderis, Holographic reconstruction of space-time and renormalization in the AdS/CFT correspondence, Commun. Math. Phys. **217**, 595 (2001).

[22] M. P. Heller, R. A. Janik, and P. Witaszczyk, Hydrodynamic gradient expansion in gauge theory plasmas, Phys. Rev. Lett. **110**, 211602 (2013).

[23] R. A. Janik and R. B. Peschanski, Gauge/gravity duality and thermalization of a boost-invariant perfect fluid, Phys. Rev. D **74**, 046007 (2006).

[24] J. K. Ghosh, R. Loganayagam, S. G. Prabhu, M. Rangamani, A. Sivakumar, and V. Vishal, Effective field theory of stochastic diffusion from gravity, J. High Energy Phys. 05 (2021) 130.

[25] M. Crossley, P. Glorioso, and H. Liu, Effective field theory of dissipative fluids, J. High Energy Phys. 09 (2017) 095.

[26] P. Glorioso, M. Crossley, and H. Liu, Effective field theory of dissipative fluids (II): Classical limit, dynamical KMS symmetry and entropy current, J. High Energy Phys. 09 (2017) 096.

[27] P. Glorioso, H. Liu, and S. Rajagopal, Global anomalies, discrete symmetries, and hydrodynamic effective actions, J. High Energy Phys. 01 (2019) 043.

[28] K. Jensen, R. Marjieh, N. Pinzani-Fokeeva, and A. Yarom, A panoply of Schwinger-Keldysh transport, SciPost Phys. **5**, 053 (2018).

[29] H. Liu and P. Glorioso, Lectures on non-equilibrium effective field theories and fluctuating hydrodynamics, PoS TASI2017, 008 (2018).

[30] J. de Boer, M. P. Heller, and N. Pinzani-Fokeeva, Holographic Schwinger-Keldysh effective field theories, J. High Energy Phys. 05 (2019) 188.

[31] B. Chakrabarty, J. Chakravarty, S. Chaudhuri, C. Jana, R. Loganayagam, and A. Sivakumar, Nonlinear Langevin dynamics via holography, J. High Energy Phys. 05 (2020) 165.

[32] T. He, R. Loganayagam, M. Rangamani, and J. Virrueta, An effective description of momentum diffusion in a charged plasma from holography, J. High Energy Phys. 01 (2022) 145.

[33] C. Jana, R. Loganayagam, and M. Rangamani, Open quantum systems and Schwinger-Keldysh holograms, J. High Energy Phys. 07 (2020) 242.

[34] T. He, R. Loganayagam, M. Rangamani, A. Sivakumar, and J. Virrueta, The timbre of Hawking gravitons: an effective description of energy transport from holography, J. High Energy Phys. 09 (2022) 092.

[35] B. C. van Rees, Real-time gauge/gravity duality and ingoing boundary conditions, Nucl. Phys. B, Proc. Suppl. **192-193**, 193 (2009).

[36] K. Skenderis, Lecture notes on holographic renormalization, Classical Quantum Gravity **19**, 5849 (2002).

[37] D. T. Son and D. Teaney, Thermal noise and stochastic strings in AdS/CFT, J. High Energy Phys. 07 (2009) 021.

[38] J. D. Bjorken, Highly relativistic nucleus-nucleus collisions: The central rapidity region, Phys. Rev. D: Part. Fields **27**, 140 (1983).

[39] S. Jeon and U. Heinz, Introduction to hydrodynamics, Int. J. Mod. Phys. E **24**, 1530010 (2015).

[40] E. Shuryak, Sang-Jin Sin, and I. Zahed, A gravity dual of RHIC collisions, J. Korean Phys. Soc. **50**, 384 (2007).

[41] R. A. Janik and R. B. Peschanski, Asymptotic perfect fluid dynamics as a consequence of Ads/CFT, Phys. Rev. D **73**, 045013 (2006).

[42] R. A. Janik, Viscous plasma evolution from gravity using AdS/CFT, Phys. Rev. Lett. **98**, 022302 (2007).

[43] S. Kinoshita, S. Mukohyama, S. Nakamura, and K.-Y. Oda, A holographic dual of Bjorken flow, Prog. Theor. Phys. **121**, 121 (2009).

[44] G. Beuf, M. P. Heller, R. A. Janik, and R. Peschanski, Boost-invariant early time dynamics from AdS/CFT, J. High Energy Phys. 10 (2009) 043.

[45] M. P. Heller and M. Spalinski, Hydrodynamics beyond the gradient eExpansion: Resurgence and resummation, Phys. Rev. Lett. **115**, 072501 (2015).

[46] A. Soloviev, Hydrodynamic attractors in heavy ion collisions: A review, Eur. Phys. J. C **82**, 319 (2022).

[47] V. Balasubramanian and P. Kraus, A stress tensor for anti-de Sitter gravity, Commun. Math. Phys. **208**, 413 (1999).

[48] A. Mukhopadhyay, Nonequilibrium fluctuation-dissipation relation from holography, Phys. Rev. D **87**, 066004 (2013).

[49] M. Dodelson, A. Grassi, C. Iossa, D. Panea Lichtig, and A. Zhiboedov, Holographic thermal correlators from supersymmetric instantons, SciPost Phys. **14**, 116 (2023).

[50] S. Banerjee, R. Iyer, and A. Mukhopadhyay, The holographic spectral function in non-equilibrium states, Phys. Rev. D **85**, 106009 (2012).

[51] Lata Kh Joshi, A. Mukhopadhyay, F. Preis, and P. Ramadevi, Exact time dependence of causal correlations and nonequilibrium density matrices in holographic systems, Phys. Rev. D **96**, 106006 (2017).

[52] C. Cartwright and M. Kaminski, Correlations far from equilibrium in charged strongly coupled fluids subjected to a strong magnetic field, J. High Energy Phys. 09 (2019) 072.

[53] J. F. Pedraza, A. Russo, A. Svesko, and Z. Weller-Davies, Sewing spacetime with lorentzian threads: Complexity and the emergence of time in quantum gravity, J. High Energy Phys. 02 (2022) 093.

[54] P. J. Martínez and G. A. Silva, Thermalization of holographic excited states, J. High Energy Phys. 03 (2022) 003.
[55] G. Basar and G. V. Dunne, Hydrodynamics, resurgence, and transasymptotics, Phys. Rev. D 92, 125011 (2015).

[56] S. Raju, Lessons from the information paradox, Phys. Rep. 943, 1 (2022).

[57] T. Kibe, P. Mandayam, and A. Mukhopadhyay, Holographic spacetime, black holes and quantum error correcting codes: A review, Eur. Phys. J. C 82, 463 (2022).

[58] S. S. Gubser, Symmetry constraints on generalizations of Bjorken flow, Phys. Rev. D 82, 085027 (2010).

[59] W. Busza, K. Rajagopal, and W. van der Schee, Heavy ion collisions: The big picture, and the big questions, Annu. Rev. Nucl. Part. Sci. 68, 339 (2018).

[60] P. Kovtun, Lectures on hydrodynamic fluctuations in relativistic theories, J. Phys. A: Math. Theor. 45, 473001 (2012).

[61] H. Liu and J. Sonner, Holographic systems far from equilibrium: A review, arXiv:1810.02367.

[62] E. Grossi, A. Soloviev, D. Teaney, and F. Yan, Soft pions and transport near the chiral critical point, Phys. Rev. D 104, 034025 (2021).

[63] A. Almheiri, R. Mahajan, J. Maldacena, and Y. Zhao, The Page curve of Hawking radiation from semiclassical geometry, J. High Energy Phys. 03 (2020) 149.

[64] Lata Kh Joshi, A. Mukhopadhyay, and A. Soloviev, Time-dependent NAdS₂ holography with applications, Phys. Rev. D 101, 066001 (2020).

[65] S. Mondkar, A. Mukhopadhyay, A. Rebhan, and A. Soloviev, Quasinormal modes of a semi-holographic black brane and thermalization, J. High Energy Phys. 11 (2021) 080.

[66] S. A. Hayward, General laws of black hole dynamics, Phys. Rev. D: Part. Fields 49, 6467 (1994).

[67] S. A. Hayward, Energy and entropy conservation for dynamical black holes, Phys. Rev. D 70, 104027 (2004).

[68] P. Figueras, V. E. Hubeny, M. Rangamani, and S. F. Ross, Dynamical black holes and expanding plasmas, J. High Energy Phys. 04 (2009) 137.

[69] J. Casalderrey-Solana, C. P. Herzog, and B. Meiring, Holographic Bjorken flow at large-$D$, J. High Energy Phys. 01 (2019) 181.