∎

11institutetext:  Crete Center for Theoretical Physics, Department of Physics,
University of Crete, Heraklion, Greece. 22institutetext: Laboratoire de
Physique de l’École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne
Université, Université de Paris, F-75005 Paris, France 33institutetext: Asia
Pacific Center for Theoretical Physics, Postech, Pohang 37673, Korea.
44institutetext: The Institute of Mathematical Sciences, Chennai 600113,
India. 55institutetext: Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar, Mumbai 400094, India 66institutetext: Instituto de
Física, Pontificia Universidad Católica de Valparaíso, Avenida Universidad
330, Valparaíso, Chile. 77institutetext: Center for Strings, Gravitation and
Cosmology, Indian Institute of Technology Madras, Chennai 600036, India
88institutetext: Faculty of Mathematics and Physics, University of Ljubljana,
Jadranska ulica 19, SI-1000 Ljubljana, Slovenia

# How Gubser flow ends in a holographic conformal theory

Avik Banerjee addr2,addr3  Toshali Mitra addr4,addr5,addr6  Ayan Mukhopadhyay
addr7,addr8  Alexander Soloviev addr1

###### Abstract

Gubser flow is an axis-symmetric and boost-invariant evolution in a
relativistic quantum field theory which is best studied by mapping
$\mathbf{R}^{3,1}$ to $dS_{3}\times\mathbf{R}$ when the field theory has
conformal symmetry. We show that at late de-Sitter time, which corresponds to
large proper time and central region of the future wedge within
$\mathbf{R}^{3,1}$, the holographic conformal field theory plasma can reach a
state in which $\varepsilon=P_{T}=-P_{L}$, with $\varepsilon$, $P_{T}$ and
$P_{L}$ being the energy density, transverse and longitudinal pressures,
respectively. We further determine the full sub-leading behaviour of the
energy-momentum tensor at late time. Restricting to flows in which the energy
density decays at large transverse distance from the central axis in
$\mathbf{R}^{3,1}$, we show that this decay should be faster than any power
law. Furthermore, in this case the energy density also vanishes in
$\mathbf{R}^{3,1}$ faster than any power as we go back to early proper time.
Hydrodynamic behavior can appear in intermediate time.

††journal: Eur. Phys. J. C

##  1 Introduction

Gubser flow Gubser:2010ze ; Gubser:2010ui  is a time-dependent evolution of a
relativistic many-body system which is boost invariant and has a rotational
symmetry about an axis. The original context of study of the Gubser flow is
the evolution of QCD matter produced in central heavy ion collisions, in which
case the axis of rotational symmetry is the beam axis.

Gubser flow in conformal field theory can be studied by mapping Minkowski
space ($\mathbf{R}^{3,1}$) to a product of three dimensional de-Sitter space
and the real line ($dS_{3}\times\mathbf{R}$), which makes the symmetries
manifest Gubser:2010ui . The most remarkable feature of Gubser flow is that at
late de-Sitter time, $\rho$, the evolution cannot be described by relativistic
hydrodynamics Gubser:2010ui ; Denicol:2014xca ; Denicol:2014tha ;
Behtash:2017wqg ; Martinez:2017ibh ; Denicol:2018pak ; Chattopadhyay:2018fzy ;
Behtash:2019qtk ; Dash:2020zqx . This large $\rho$ regime is large proper time
and the central region of the flow in the future wedge of the collision in
Minkowski space.

Our study of the holographic Gubser flow sets a primary example in which one
can exactly determine how a system can evolve out of hydrodynamic regime due
to expansion of some directions of space in a quantum field theory, and is
relevant to understanding evolution of quantum matter in Kasner geometries
Berger:2002st  describing spacetime singularities, etc.

In kinetic theories without particle production in the Gubser flow, it has
been shown Denicol:2014xca ; Denicol:2014tha ; Behtash:2017wqg ;
Martinez:2017ibh ; Denicol:2018pak ; Chattopadhyay:2018fzy ; Behtash:2019qtk ;
Dash:2020zqx  that the large $\rho$ behavior is free-streaming where inter-
particle interactions disappear. However, such ultra dilute regime is better
studied in quantum field theory. Here, we show that in a holographic conformal
theory, the deconfined plasma can reach a phase in which the energy density
$\varepsilon$, the transverse pressure $P_{T}$ and the longitudinal pressure
$P_{L}$ satisfy $P_{L}=-P_{T}$, and $\varepsilon=P_{T}$ in the large $\rho$
limit independently of the details of the state of the system at an early
time. Remarkably, this phase is similar to the color glass condensate regime
described by perturbative physics of saturated gluons Epelbaum:2013ekf ;
Muller:2019bwd ; Berges:2020fwq .

We also establish a systematic late-time expansion which gives all sub-leading
corrections to the colour glass condensate like behavior at late time. As in
the fluid-gravity correspondence Rangamani:2009xk , this late time behavior is
deduced from demanding that the future horizon is regular at each order in the
expansion. However, this late-time expansion is not a gradient expansion which
breaks down at late time. It captures all systematic corrections in powers of
$\exp(-\rho/L)$, where $\rho$ is the de-Sitter time and $L$ is the radius of
the $S^{2}$ factor at $\rho=0$. The dimensionful parameter analogous to $L$ in
$R^{3,1}$ is $q^{-1}$ which is a scale set by the initial conditions, and
which gives a systematic expansion in negative powers of $q\tau$ (with $\tau$
being the proper time) in $R^{3,1}$.

In physically realizable Gubser flows in $R^{3,1}$, the energy density should
decay at large transverse distance from the central axis at any proper time.
We show that in such flows the asymptotic late de-Sitter time behavior
obtained from bulk regularity implies that the decay of the energy density at
large transverse distance from the central axis in $R^{3,1}$ faster than any
power. Furthermore, we also show that in the limit of early de-Sitter time
$\rho\rightarrow 0$ and early proper time $\tau\rightarrow 0$, the energy
density should vanish faster than any power of $\exp(\rho/L)$ and $\tau$,
respectively. The systematic late de-Sitter time expansion itself carries the
full information of the state (represented by the initial conditions of the
gravitational solution at a finite de-Sitter time) when the symmetries of the
Gubser flow are fully preserved.

We postpone discussions about the relevance of our results to heavy-ion
collisions in the concluding section.

The plan of this paper is as follows. In Section 2, we discuss the map from
the future wedge of the collision in Minkowski space to
$dS_{3}\times\mathbf{R}$, and review the breakdown of the gradient
(hydrodynamic) expansion at late time. In Section 3, we doscuss the basic
setup of how the Gubser flow can be studied via dual gravitational dynamics in
strongly coupled large $N$ holographic conformal theories. In Section 4, we
demonstrate how the late de-Sitter time behavior of the Gubser flow can be
determined systematically and present general perturbative solution in a
perturbative expansion. In Section 5, we show that by assuming that the energy
profile decays in transverse directions to the axis in Minkowski space, we can
also determine the early de-Sitter time behavior and how the late de-Sitter
time expansion itself determines the energy density in the full future wedge
of the collision in Minkowski space. We show that the energy density should
decay at large transverse directions to the axis in Minkowski space faster
than any power, and also vanish faster than any power as we approach the
proper time $\tau=0$ which is the moment of the collision. In Section 6, we
discuss different implications of our results and present an outlook.

##  2 Preliminaries

For a conformal system on $\mathbf{R}^{3,1}$, the full symmetry of the Gubser
flow is $SO(3)_{q}\otimes SO(1,1)\otimes Z_{2}$. These symmetries can be made
manifest by performing diffeomorphism and Weyl rescaling of $\mathbf{R}^{3,1}$
to $dS_{3}\times\mathbf{R}$ Gubser:2010ui . The $SO(1,1)$ boost symmetry acts
additively and $Z_{2}$, the reflection symmetry about the collision plane,
acts reflexively on the $\mathbf{R}$ factor, which is physically the rapidity.
Furthermore, $dS_{3}$ is a contracting and then expanding $S^{2}$ on which
$SO(3)$ has a natural action. The variable $q$ parametrizes inequivalent ways
of embedding $SO(3)\otimes SO(1,1)$ in the full conformal group $SO(4,2)$.
Also, $q^{-1}$ has the dimension of length and is essentially the size of the
colliding systems.

The explicit map from $\mathbf{R}^{3,1}$ to $dS_{3}\times\mathbf{R}$ is via
the Milne coordinates, which are the proper time $\tau=\sqrt{t^{2}-z^{2}}$,
the rapidity $\eta={\rm arctan}(z/t)$, the radial coordinate $x_{\perp}$ and
the angular coordinate $\phi$ of the plane transverse to the $z$-axis (the
beam axis). In Milne coordinates, the metric on $\mathbf{R}^{3,1}$ is

${\rm d}s^{2}=-{\rm d}\tau^{2}+\tau^{2}{\rm d}\eta^{2}+{\rm
d}x_{\perp}^{2}+x_{\perp}^{2}{\rm d}\phi^{2}.$ (1)

The $dS_{3}\times\mathbf{R}$ coordinates are $\rho$, $\theta$, $\phi$ and
$\eta$, where

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle-L\,\,{\rm
arcsinh}\left(\frac{1-q^{2}\tau^{2}+q^{2}x_{\perp}^{2}}{2q\tau}\right),$ (2)
$\displaystyle\theta$ $\displaystyle=$ $\displaystyle{\rm
arctan}\left(\frac{2qx_{\perp}}{1+q^{2}\tau^{2}-q^{2}x_{\perp}^{2}}\right).$
(3)

Crucially, the future wedge, which is the causal future of the plane of the
collision at $\tau=0$, is contained in $dS_{3}\times\mathbf{R}$. The
corresponding metric

$\displaystyle{\rm d}\hat{s}^{2}$ $\displaystyle=\frac{L^{2}}{\tau^{2}}{\rm
d}s^{2}$ (4) $\displaystyle=-{\rm
d}\rho^{2}+L^{2}\cosh^{2}\left(\frac{\rho}{L}\right)\left({\rm
d}\theta^{2}+\sin^{2}\theta{\rm d}\phi^{2}\right)+L^{2}{\rm d}\eta^{2},$

is Weyl equivalent to the flat metric. From our discussion above, it is clear
that $\rho$ is invariant under $SO(3)_{q}\otimes SO(1,1)\otimes Z_{2}$, and
therefore physical quantities like the energy density should depend only on
$\rho$ explicitly up to the Weyl rescaling. We will denote physical quantities
measured in $dS_{3}\times\mathbf{R}$ with a hat.

The conformal Ward identites for the energy-momentum tensor are:

$\hat{\nabla}_{\mu}\hat{T}^{\mu}_{\,\,\nu}=0,\quad\hat{T}^{\mu}_{\,\,\mu}=\mathcal{A},$
(5)

where $\mathcal{A}$ is the Weyl anomaly. The $SO(3)_{q}\otimes SO(1,1)\otimes
Z_{2}$ symmetries and these Ward identities imply that $\hat{T}^{\mu}_{\nu}$
should be of the form

$\hat{T}^{\mu}_{\,\,\nu}=\hat{t}^{\mu}_{\,\,\nu}+\hat{\mathcal{A}}^{\mu}_{\,\,\nu}$
(6)

where

$\displaystyle\hat{t}^{\rho}_{\,\,\rho}=-\hat{\varepsilon}\left(\frac{\rho}{L}\right),$
(7)
$\displaystyle\hat{t}^{\theta}_{\,\,\theta}=\hat{t}^{\phi}_{\,\,\phi}=\hat{P}_{T}(\rho)=-\frac{1}{2}\coth\left(\frac{\rho}{L}\right)\hat{\varepsilon}^{\prime}\left(\frac{\rho}{L}\right)-\hat{\varepsilon}\left(\frac{\rho}{L}\right),$
(8)
$\displaystyle\hat{t}^{\eta}_{\,\,\eta}=\hat{P}_{L}(\rho)=\coth\left(\frac{\rho}{L}\right)\hat{\varepsilon}^{\prime}\left(\frac{\rho}{L}\right)+3\hat{\varepsilon}\left(\frac{\rho}{L}\right).$
(9)

Here, ′ denotes a derivative w.r.t. to the argument. It is easy to check that
$\hat{t}^{\mu}_{\,\,\nu}$ is separately conserved (i.e.
$\hat{\nabla}_{\mu}\hat{t}^{\mu}_{\,\,\nu}=0$) and is traceless. The anomalous
term $\hat{\mathcal{A}}^{\mu}_{\,\,\nu}$ (which is state-independent) is
explicitly

$\displaystyle\hat{\mathcal{A}}^{\mu}_{\,\,\nu}=\frac{N^{2}}{32\pi^{2}}{\rm
diag}\left(1,1,1,-3\right)L^{-4}$ (10)

for $\mathcal{N}=4$ $SU(N)$ super Yang-Mills-theory (SYM) Henningson:1998gx ;
Balasubramanian:1999re . We note that
$\mathcal{A}=\hat{\mathcal{A}}^{\mu}_{\,\,\mu}=0$.

It is clear from the above that the only input from the microscopic dynamics
that is needed to determine the energy-momentum tensor is the evolution of the
energy density $\hat{\varepsilon}(\rho)$, because the latter determines the
transverse pressure $\hat{P}_{T}$ and the longitudinal pressure $\hat{P}_{L}$
via (8) and (9), respectively. Weyl transformation yields the energy-momentum
tensor in Minkowski space:

${T}^{\mu}_{\,\,\nu}=\frac{L^{4}}{\tau^{4}}\hat{t}^{\mu}_{\,\,\nu},\,\,{\rm
i.e.}\,\,\varepsilon=L^{4}\tau^{-4}\hat{\varepsilon}.$ (11)

Note that the anomalous term disappears.

In the hydrodynamic regime, the general form of the energy-momentum tensor
given by (7)-(9) implies that the fluid is static in the ${\rm
d}S_{3}\times\mathbf{R}$ frame, i.e. the velocity field is of the form
$\hat{u}^{\rho}=1$ with other components vanishing. The hydrodynamic equations
determine $\hat{\varepsilon}(\rho)$. For a perfect conformal fluid,
$\hat{\varepsilon}(\rho)=\hat{\varepsilon}_{0}\cosh^{-8/3}(\rho/L)$, whereas
for a viscous flow (see A for an explicit solution), $\hat{\varepsilon}$ goes
to a constant at large $\rho$ indicating the breakdown of the derivative
expansion. However, at $\rho\sim 0$, the Knudsen number is small, and the
derivative expansion is equivalent to

$\hat{\varepsilon}(\rho)=\hat{\varepsilon}_{0}-\frac{4}{3}\hat{\varepsilon}_{0}\frac{\rho^{2}}{L^{2}}+\frac{16}{27}H_{0}\hat{\varepsilon}_{0}^{3/4}\frac{\rho^{3}}{L^{4}}+\mathcal{O}(\rho^{4}),$
(12)

with $H_{0}$ being the dimensionless (and constant) shear viscosity (see A).
This hydrodynamic evolution is only transitory.

##  3 Holographic Gubser flow

Any state in the universal sector of a holographic conformal gauge theory,
such as $\mathcal{N}=4$ SYM with $SU(N)$ gauge group, has a dual description
as a solution of pure classical Einstein’s gravity with a negative
cosmological constant in one higher spacetime dimension in infinite ’t Hooft
coupling and large $N$ limits Maldacena:1997re ; Witten:1998qj ; Gubser:1998bc
. Such a solution should be regular, i.e. without any naked singularity, and
with boundary metric $g^{(b)}_{\mu\nu}$ defined in terms of the bulk metric
$G_{\mu\nu}$ via

$\lim_{r\rightarrow 0}\frac{r^{2}}{l^{2}}G_{\mu\nu}\equiv g^{(b)}_{\mu\nu},$
(13)

coinciding with the physical background metric on which the gauge theory
lives. Above, $r$ is the holographic radial coordinate, which has an
interpretation in terms of the energy scale of the dual theory Freedman:1999gp
; Heemskerk:2010hk ; Bredberg:2010ky ; Kuperstein:2011fn ; Kuperstein:2013hqa
; Behr:2015aat ; Behr:2015yna ; Mukhopadhyay:2016fre , while $\mu$ and $\nu$
indices represent the directions spanning the boundary of spacetime at $r=0$.
The bulk cosmological constant is $\Lambda=-6/l^{2}$.

Any five-dimensional gravitational solution dual to a Gubser flow of the dual
gauge theory in $dS_{3}\times\mathbf{R}$ can be written in the form

$\displaystyle{\rm d}s^{2}=-\frac{2l^{2}}{r^{2}}{\rm d}r{\rm
d}\rho-\frac{l^{2}}{r^{2}}\left(1-\frac{r^{2}}{L^{2}}+A\left(\frac{r}{L},\frac{\rho}{L}\right)\right){\rm
d}\rho^{2}$
$\displaystyle+\frac{l^{2}(L\cosh(\rho/L)+r\sinh(\rho/L))^{2}}{r^{2}}$
$\displaystyle\times e^{B\left(\frac{r}{L},\frac{\rho}{L}\right)}\,\Big{(}{\rm
d}\theta^{2}+\sin^{2}\theta{\rm d}\phi^{2}\Big{)}$
$\displaystyle+\frac{l^{2}}{r^{2}}e^{C\left(\frac{r}{L},\frac{\rho}{L}\right)-2B\left(\frac{r}{L},\frac{\rho}{L}\right)}\,L^{2}{\rm
d}\eta^{2}.$ (14)

We have chosen the ingoing Eddington-Finkelstein gauge in which the regularity
of the future horizon can be readily examined. For the boundary metric to be
the metric (4) on ${\rm d}S_{3}\times\mathbf{R}$, we need

$A(r=0,\rho/L)=B(r=0,\rho/L)=C(r=0,\rho/L)=0.$ (15)

The solution dual to the vacuum state corresponds to $A=B=C=0.$ This solution
is locally $AdS_{5}$ (it can be explicitly checked that it is the maximally
symmetric metric). Although the metric dual to the vacuum of the theory in
$dS_{3}\times\mathbf{R}$ has all the $AdS_{5}$ isometries, we can readily see
from (3) that the discrete symmetry $\rho\rightarrow-\rho$ of
$dS_{3}\times\mathbf{R}$ at the boundary is not preserved due to the
$\sinh(\rho/L)$ term in the $\rho\rho$-component of the metric. As a result,
the vacuum state itself has a monotonically growing entropy as discussed in D.
We note that the map from $R^{3,1}$ to $dS_{3}\times\mathbf{R}$ is also
defined only for the future wedge $\tau>0$, so it also prefers an arrow of
time. Therefore from the Minkowski point of view, the monotonically growing
entropy in the vacuum of $dS_{3}\times\mathbf{R}$ could capture something
similar to the entropy of the Rindler observer. We have presented some
discussions on this issue in D.

Generally, as detailed in B, $A$, $B$ and $C$ have radial expansions of the
form

$\displaystyle
A=\sum_{n=0}^{\infty}a_{(n)}(\rho/L)\frac{r^{n}}{L^{n}},\,\,B=\sum_{n=0}^{\infty}b_{(n)}(\rho/L)\frac{r^{n}}{L^{n}},$
$\displaystyle C=\sum_{n=0}^{\infty}c_{(n)}(\rho/L)\frac{r^{n}}{L^{n}}$ (16)

Plugging these in Einstein’s equation, we find that the entire solution is
determined just by two inputs: $a_{(1)}(\rho/L)$, which is related to a proper
residual diffeomorphism that does not affect the boundary data, and
$a_{(4)}(\rho/L)$ which is physical and should be chosen such that no naked
singularity is present.

Physically $a_{(4)}(\rho/L)$ determines with the energy-momentum tensor of the
dual gauge theory state. The latter can be obtained systematically via the
standard holographic dictionary (see B) Henningson:1998gx ;
Balasubramanian:1999re ; deHaro:2000vlm ; Skenderis:2002wp . Noting that
$l^{3}/G_{N}=(2/\pi)N^{2}$ Maldacena:1997re , we find that the dual energy-
momentum tensor has the same form as that given by (6)-(10) (including the
anomalous term) with the identification

$\hat{\varepsilon}\left(\frac{\rho}{L}\right)=-\frac{3l^{3}}{16\pi
G_{N}}a_{(4)}\left(\frac{\rho}{L}\right)=-\frac{3N^{2}}{8\pi^{2}}a_{(4)}\left(\frac{\rho}{L}\right).$
(17)

Since a Weyl transformation at the boundary can be lifted to a bulk
diffeomorphism Henningson:1998gx ; deHaro:2000vlm , we can readily obtain the
Minkowski boundary metric (1) and the energy-momentum tensor (11) via an
appropriate bulk diffeomorphism.

##  4 Finding the late time solution in ${\rm d}S_{3}\times\mathbf{R}$

The problem of finding generic late time behavior amounts to setting up a
suitable late time expansion and determining the set of conditions which lead
to a regular future horizon in this expansion perturbatively. Furthermore, the
solution has to be normalizable, i.e. satisfy (15). Due to the exponential
expansion of $S^{2}$, the energy density in ${\rm d}S_{3}\times\mathbf{R}$
should dilute and we should eventually reach the vacuum. Therefore, in the
dual gravitational solution (3), $A$, $B$ and $C$ must eventually vanish. The
latter solution with $A=B=C=0$ has a horizon at $r=L$ where $\partial_{\rho}$
has vanishing norm. (Note this is not a Killing horizon.) Therefore at large
$\rho$, the horizon should be at $r=L$ at leading order in the late time
expansion, and just like in the case of fluid-gravity correspondence
Rangamani:2009xk , we should require that at each order in the perturbative
expansion, the behavior of $A$, $B$ and $C$ are smooth at $r=L$ implying
regularity of the future horizon.

It is instructive to first discuss the example of a massless scalar field
$\Phi$ which satisfies the Klein-Gordon equation $\nabla^{2}\Phi=0$, where
$\nabla^{2}$ is the Laplacian operator in the five-dimensional metric (3) dual
to the vacuum and with $A=B=C=0$. Demanding $SO(3)_{q}\otimes SO(1,1)\otimes
Z_{2}$ symmetry amounts to requiring that $\Phi$ depends only on $r$ and
$\rho$. It is easy to see that the background five dimensional metric is a
rational function of $\sigma=\exp(\rho/L)$. Therefore, we will expect that at
late time $\Phi(r,\rho)$ should behave as
$\sigma^{-\alpha}\sim\exp(-\alpha\rho/L)$.The consistent ansatz which fits
this late time behavior (at large $\sigma$) is

$\Phi(r,\rho)=\sigma^{-\alpha}\sum_{n=0}^{\infty}\phi_{n}(v)\sigma^{-2n}$ (18)

where $v=r/L$. However, if the solution is analytic at the horizon $r=L$, we
should expect the behavior near $r=L$ (i.e. $v=1$) to be given by

$\phi_{n}(v)=\sum_{m=0}^{\infty}\lambda_{n,m}(1-v)^{m},$ (19)

where $\lambda_{n,m}$ are pure numbers. Substituting (18) and (19) in the
Klein-Gordon equation, we can readily solve all other $\lambda_{n,m}$ in terms
of $\lambda_{0,m}$ which can be chosen freely. Finally, we would require the
condition of normalizability, i.e.

$\phi_{n}(v=0)=0.$ (20)

Solving $\lambda_{n,0}$ explicitly in terms of $\lambda_{0,0}$, we would get

$\phi_{0}(v=0)=\lambda_{0,0}f_{0}(\alpha).$ (21)

The allowed values of $\alpha$ are simply the roots of $f_{0}(\alpha)$.
Solving for $\lambda_{n,0}$ to very high orders, we get stable roots, which
are $4+2k$, for $k=0,1,2,\ldots$ From the ansatz (18) and the allowed values
of $\alpha$, it follows that if $\phi_{0}$ is normalizable then so is
$\phi_{n}$ for $n\geq 1$ provided one of the two integration constants at each
order is chosen so that $\phi_{n}(0)=0$. Since $\phi_{0}$ determines the
source term in the linear equation for $\phi_{1}$ (and so on) for which a
normalizable particular solution exists, and the homogeneous normalizable
solution also exists since it coincides with one of the allowed values of
$\alpha$, namely $\alpha=4+2k$ for $k\geq 1$, it follows that the generic late
time behavior is simply (18) with $\alpha=4$. At each order in the expansion
we thus have one independent integration constant after ensuring
$\phi_{n}(0)=0$ giving a normalizable solution which is regular at the future
horizon $v=1$. Explicitly,

$\displaystyle\phi_{0}(v)=\hat{o}_{0}\frac{v^{4}}{(1+v)^{4}},$
$\displaystyle\phi_{1}(v)=\hat{o}_{1}\frac{v^{4}(3+v^{2})}{3(1+v)^{6}}+\hat{o}_{0}\frac{4v^{6}}{3(1+v)^{6}},\,\,{\rm
etc.}$ (22)

The expectation value (e.v.) of the marginal operator $\hat{O}$ in the gauge
theory dual to the massless scalar field is essentially the $r^{4}$
coefficient in the radial expansion of $\Phi$. Ignoring backreaction of the
dual bulk scalar field, the generic behavior of the vacuum expectation value
(v.e.v.) of this operator with Gubser flow symmetries is thus

$\hat{O}(\rho)\sim\sigma^{-4}\sum_{k=0}^{\infty}\hat{o}_{k}\sigma^{-2k},$ (23)

where $\hat{o}_{k}$ are arbitrary real numbers.

The solutions for $\phi_{n}(v)$ turn out to be rational functions.
Furthermore, with specific choices of integration constants for the subleading
$\phi_{n}$, we can sum over the entire late time expansion as well. For
instance, when the leading behaviour is $\sigma^{-4}$ (i.e. $\hat{o}_{4}\neq
0$), the full summation with a specific set of choices yields

$\Phi(r,\rho)=\Gamma_{4}\frac{\sigma^{-4}v^{4}}{(1+v)^{5}}\frac{(1+v)^{2}\sigma^{2}-3-v^{2}}{1+\sigma^{2}+v(\sigma^{2}-1)},$
(24)

an exact normalizable solution of the Klein-Gordon equation with an arbitrary
constant coefficient $\Gamma_{4}$. When the leading behaviour is $\sigma^{-6}$
(i.e. $\hat{o}_{4}=0$), we similarly obtain

$\Phi(r,\rho)=\Gamma_{6}\frac{\sigma^{-6}v^{4}}{(1+v)^{7}}\frac{(1+v)^{2}(3+v^{2})\sigma^{2}-6-8v^{2}-v^{4}}{1+\sigma^{2}+v(\sigma^{2}-1)},$
(25)

with an arbitrary constant coefficient $\Gamma_{6}$, and so on. Thus for each
leading $\sigma^{-(4+2k)}$ behavior with $k=0,1,2,\ldots$, we obtain an exact
normalizable solution with arbitrary constant coefficients $\Gamma_{4+2k}$.

We readily note that both (24) and (25) diverge in the infinite past
$\sigma=0$ (i.e. for $\rho\rightarrow-\infty$) at the horizon $v=1$, and the
corresponding $\hat{O}(\rho)$ also blows up. However, since we are setting up
initial conditions at an arbitrary but finite time and looking into the
future, this does not bother us, as in the case of fluid-gravity
correspondence where generically we get singularities on the past horizon of
the leading order static black brane geometry Rangamani:2009xk ; Gupta:2008th
even though the future horizon is regular for appropriate choice of transport
coefficients. This singularity in the infinite past can be cured by
appropriate initial conditions. Thus, (23) gives the generic late-time
behavior of $\hat{O}(\rho)$ in the vacuum.

We can replicate the same strategy in full non-linear pure gravity with the
following ansatz consistent with the gravitational equations

$A(r,\rho)=\sum_{n=1}^{\infty}\sum_{m=0}^{\infty}\sigma^{-n\alpha-2m}A_{nm}(v),$
(26)

and with similar expansions for $B(r,\rho)$ and $C(r,\rho)$ with coefficients
$B_{nm}(v)$ and $C_{nm}(v)$, respectively. Normalizability implies that
$A_{nm}(v)$, $B_{nm}(v)$ and $C_{nm}(v)$ should vanish at $v=0$. Compared to
(18), we get a double summation here because of the non-linearity in the
gravitational equations. As detailed in C, both normalizability and regularity
at the late-horizon $v=1$ are obtained when $\alpha=0,2,4,6,\ldots$. However,
we need to ensure that we get non-vanishing solutions which are not pure
gauge. We have been able to find such physical solutions for

$\alpha=4+2k,\,\,k=0,1,2,\ldots,$ (27)

exactly as in the case of the massless scalar field. This implies that we can
simplify the expansion (26) to

$A(r,\rho)=\sum_{k=0}^{\infty}\sigma^{-4-2k}a_{k}(v),$ (28)

with similar expansions for $B(r,\rho)$ and $C(r,\rho)$ with coefficients
$b_{k}(v)$ and $c_{k}(v)$ respectively. Remarkably, these are of the same
forms as (18) for the massless scalar case. It follows that the late de-Sitter
time behavior of energy density is

$\hat{\varepsilon}(\rho)=-\frac{3l^{3}}{16\pi
G_{N}}\sigma^{-4}L^{-4}\sum_{k=0}^{\infty}\hat{e}_{k}\sigma^{-2k},$ (29)

with $\hat{e}_{k}$ being arbitrary real numbers as $\hat{o}_{k}$ are in (23).
The normalizable solutions which are regular at the late-time horizon $v=1$
and correspond to such late-time behavior (29) are given by the following
explicit functions appearing in (28) up to $\mathcal{O}(\sigma^{-8})$:

$\displaystyle a_{0}(v)$
$\displaystyle=\frac{v^{4}}{(1+v)^{2}}\hat{e}_{0},\,\,b_{0}(v)=-\frac{v^{4}}{2(1+v)^{4}}\hat{e}_{0},\,\,$
$\displaystyle c_{0}(v)$ $\displaystyle=\frac{2v^{5}}{5(1+v)^{4}}\hat{e}_{0},$
(30) $\displaystyle a_{1}(v)$
$\displaystyle=\frac{12\hat{e}_{0}v^{6}}{5(v+1)^{4}}+\frac{\hat{e}_{1}\left(9v^{2}+10v+10\right)v^{4}}{10(v+1)^{4}},\,\,$
$\displaystyle b_{1}(v)$
$\displaystyle=\frac{\hat{e}_{0}\left(-8v^{2}+4v-15\right)v^{4}}{5(v+1)^{6}}+\frac{\hat{e}_{1}\left(-7v^{2}+4v-25\right)v^{4}}{20(v+1)^{6}},$
$\displaystyle c_{1}(v)$
$\displaystyle=\frac{4\hat{e}_{0}\left(6v^{2}-7v+14\right)v^{5}}{35(v+1)^{6}}+\frac{\hat{e}_{1}\left(8v^{2}+7v+35\right)v^{5}}{35(v+1)^{6}},$
(31) $\displaystyle a_{2}(v)$
$\displaystyle=-\frac{\hat{e}_{0}^{2}v^{8}}{14(v+1)^{6}}+\frac{2\hat{e}_{0}\left(7v^{2}-4v-6\right)v^{6}}{5(v+1)^{6}}$
$\displaystyle+\frac{6\hat{e}_{1}\left(13v^{2}+14v+21\right)v^{6}}{35(v+1)^{6}}$
$\displaystyle+\frac{\hat{e}_{2}\left(29v^{4}+70v^{3}+126v^{2}+70v+35\right)v^{4}}{35(v+1)^{6}},$
$\displaystyle b_{2}(v)$
$\displaystyle=\frac{5\hat{e}_{0}^{2}v^{8}}{28(v+1)^{8}}$
$\displaystyle+\frac{\hat{e}_{0}\left(-91v^{4}+68v^{3}-364v^{2}-28v-105\right)v^{4}}{35(v+1)^{8}}$
$\displaystyle+\frac{\hat{e}_{1}\left(-82v^{4}+84v^{3}-588v^{2}+56v-315\right)v^{4}}{70(v+1)^{8}}$
$\displaystyle+\frac{\hat{e}_{2}\left(-19v^{4}+20v^{3}-168v^{2}+28v-140\right)v^{4}}{70(v+1)^{8}},$
$\displaystyle c_{2}(v)$
$\displaystyle=\frac{\hat{e}_{0}^{2}\left(-42v^{2}-250v-225\right)v^{8}}{1050(v+1)^{8}}$
$\displaystyle+\frac{4\hat{e}_{0}\left(17v^{4}-63v^{3}+78v^{2}-84v-42\right)v^{5}}{105(v+1)^{8}}$
$\displaystyle+\frac{2\hat{e}_{1}\left(8v^{4}-2v^{3}+66v^{2}-21v+28\right)v^{5}}{35(v+1)^{8}}$
$\displaystyle+\frac{2\hat{e}_{2}\left(8v^{4}+15v^{3}+90v^{2}+42v+84\right)v^{5}}{105(v+1)^{8}}.$
(32)

It is obvious from the above expressions that the perturbative solution is
regular at $v=1$ for arbitrary values of $\hat{e}_{0}$, $\hat{e}_{1}$, etc.
When $\hat{e}_{0}\neq 0$, and (8) and (9) give
$\hat{\varepsilon}=\hat{P}_{T}=-\hat{P}_{L}$ in the limit
$\sigma\rightarrow\infty$, i.e. for $\rho\gg L$. We note that $\hat{e}_{n}$
can be of either sign. For the demonstration of monotonic entropy growth, see
D.

As detailed further in Appendix C, we have not been able to rule out that no
physical solution exists for the case $\alpha=2$. We have been able to find
only a pure gauge solution at the leading order in this case. If a physical
solution does exist, then it will give the leading late time behavior for
generic states implying that $P_{T}/\varepsilon\rightarrow 0$ and
$P_{L}/\varepsilon\rightarrow 1$ exactly like in kinetic theories.
Nevertheless, our results above shows that with fine-tuned initial conditions,
we can obtain a novel behavior at late time which is not realizable in kinetic
theories where the energy density and the pressures should be positive. (Note
even if $\alpha=2$ corresponds to physical solutions, we can get negative
$\varepsilon$ and $P_{L}$ of equal magnitude, which is also not realizable in
kinetic theories.)

##  5 Beyond late de-Sitter time: On the behaviour in the entire future wedge
of Minkowski space

To understand the evolution from the Minkowski point of view, we first recall
how the future wedge $\tau>0$ of $\mathbf{R}^{3,1}$ maps to
$dS_{3}\times\mathbf{R}$. For this purpose, it is useful to refer to Fig. 1
where $\rho/L$ is plotted as a function of $qx_{\perp}$ for fixed values of
$q\tau$ using (2). We note that for $q\tau<0.3$, the system is in the early
regime where $\rho/L$ is of large magnitude and negative. For $0.3<q\tau<3$,
we obtain $-1<\rho/L<1$ in the central region where hydrodynamic behavior may
be expected. Also, for $q\tau\gg 3$, there can be a hydrodynamic regime on any
constant $\tau$ slice, within the annulus $q\tau-1<qx_{\perp}<q\tau+1$ where
$-1<\rho/L<1$. The large $qx_{\perp}$ region is always in the early regime
from the de-Sitter point of view. In the domain $q\tau\gg 3$ and
$x_{\perp}<\tau$, we obtain $\rho/L\gg 1$, and here our late time behavior
(29) is valid.

Figure 1: $\rho/L$ plotted as a function of $qx_{\perp}$ for various fixed
values of $q\tau$. The gray dashed lines at $-1$ and 1 mark the boundary of
the applicability of a hydro-like region (shaded blue).

It is useful to first examine the case of the free bulk scalar field for which
we have exact solutions such as (24). In the latter solution, the dual
operator behaves as

$\hat{O}(\sigma)\sim\sigma^{-4}\frac{\sigma^{2}-3}{\sigma^{2}-1}$ (33)

in $dS_{3}\times\mathbf{R}$ (recall that $\sigma=\exp(\rho/L)$). The
corresponding behaviour of ${O}(\sigma)$ on the future wedge of
$\mathbf{R}^{3,1}$ can be obtained after the necessary Weyl scaling:

${O}(\tau,x_{\perp})\sim\tau^{-4}\hat{O}(\sigma(\tau,x_{\perp})).$

Utilizing (2), we readily see that at large $x_{\perp}$, we obtain that

$O\sim\frac{x_{\perp}^{8}}{\tau^{8}},$

implying that the vacuum expectation value (v.e.v.) diverges at large distance
from the central axis on the entire future wedge. Furthermore as
$\tau\rightarrow 0$, we obtain that

$O\sim-\frac{(1+x_{\perp}^{2})^{4}}{\tau^{8}},$

implying that the v.e.v. diverges as we approach the initial time. The lesson
is that although (24) is a valid solution for the theory defined on
$dS_{3}\times\mathbf{R}$, it does not render a solution which can be realized
with physical initial conditions for the theory on the future wedge of
Minkowski space. The same can be said about the exact solution (25) as well.

Clearly, we need to constrain the behaviour in $dS_{3}\times\mathbf{R}$ so
that we can map to physically realizable solutions in Minkowski space. This
can be readily achieved by a re-arrangement of the late-time expansion of
$\hat{O}$ given by (23) in a way such that we can ensure that $O$ decays in
Minkowski space at large $x_{\perp}$. Each term in the sum (23), which is of
the form $\exp(-(4+2n)\rho/L)$ with $n$ a non-negative integer, grows with
$x_{\perp}$ in Minkowski space as $(x_{\perp}^{2}/\tau)^{4+2n}$ at any fixed
$\tau>0$. It is easy to see that the following expansion

$\displaystyle\hat{O}(\rho)=\sum_{k=0}^{\infty}\left(\frac{o^{a}_{k}}{\cosh{(\rho/L)}^{4+2k}}+\frac{o^{b}_{k}}{\left(2\cosh{(\rho/L)}+\sinh{(\rho/L)}\right)^{4+2k}}\right)$
(34)

is equivalent to the general large $\rho$ expansion (23) for $\rho\gg L$, but
with each term that decays along the transverse directions in Minkowski space
at any fixed $\tau>0$. (Note ${\rm sech}(\rho/L)$ decays with large
$x_{\perp}$ as $x_{\perp}^{-2}$ at any fixed $\tau>0$, etc.) In the second set
of terms above, we have used $2\cosh(x)+\sinh(x)$ instead of $\sinh(x)$, since
then the denominator does not vanish at any value of $\rho$. Also note that
without the second set of terms we get an additional $\rho\rightarrow-\rho$
(de-Sitter time-reversal) symmetry which is not a feature of generic Gubser
flow. Of course, there is nothing unique about the choice of basis in (34). In
fact, instead of $2\cosh(x)+\sinh(x)$, we can choose $3\cosh(x)+\sinh(x)$,
etc.

We readily note that (34) has double the number of coefficients when compared
to the late time expansion (23), as for each $\hat{o}_{k}$ in (23), we have a
pair $o^{a}_{k}$ and $o^{b}_{k}$ in (34). By expanding at large $\rho$, it is
easy to see that

$\hat{o}_{0}=16\left(o^{a}_{0}+\frac{o^{b}_{0}}{81}\right),\quad\hat{o}_{1}=64\left(o^{a}_{1}+\frac{o^{b}_{1}}{729}-o^{a}_{0}-\frac{o^{b}_{0}}{243}\right),\,\,{\rm
etc.}$ (35)

Each $\hat{o}_{k}$ is a linear combination of $o^{a}_{k}$ and $o^{b}_{k}$, and
lower order terms. If we expand at early de-Sitter time
$\rho\rightarrow-\infty$ (i.e. $\sigma\rightarrow 0$), we obtain that

$\displaystyle\hat{O}(\rho)=\sigma^{4}\sum_{k=0}^{\infty}\hat{o}^{e}_{k}\sigma^{2k},$
(36)

with

$\hat{o}^{e}_{0}=16\left(o^{a}_{0}+o^{b}_{0}\right),\quad\hat{o}^{e}_{1}=64\left(o^{a}_{1}+o^{b}_{1}-o^{a}_{0}-3o^{b}_{0}\right),\,\,{\rm
etc.}$ (37)

Each $\hat{o}^{e}_{k}$ is a linear combination of $o^{a}_{k}$ and $o^{b}_{k}$,
and lower order terms. Therefore, the doubling of terms in (34) is due to
parametrization of both late-time and early-time behavior. Furthermore, the
form (34) allows for systematic expansions only when
$|\rho|\rightarrow\infty$, i.e. at early and late time; and therefore this
form is useful in only these limits. Our argument for considering the
expansion (34) shows that the early time expansion (36) for the v.e.v. at
$\sigma\sim 0$ is automatically imposed by demanding that the v.e.v. decays at
large $x_{\perp}$ in Minkowski space and behaves as (23) at late de-Sitter
time. In case of the exact solution (24), we see that the v.e.v. diverges as
$\sigma^{-4}$ as $\sigma\rightarrow 0$, and therefore the v.e.v. also blows up
at large $x_{\perp}$ in Minkowski space. If the v.e.v. has to decay in
Minkowski space at large $x_{\perp}$ and behave as (23) at late de-Sitter
time, then it has to vanish as $\sigma^{4}$ or faster when we go back to early
de-Sitter time as evident from (36).

We will soon see that bulk regularity implies that $\hat{o}^{e}_{k}$ is
related to $\hat{o}_{k}$, i.e. $o^{a}_{k}$ and $o^{b}_{k}$ are not independent
of each other, and therefore the systematic late-time expansion of the v.e.v.
(23) in absence of sources contains the full information of initial conditions
(given by the radial profile of the bulk scalar field at any fixed $\rho$).
However, this inference is possible, only when in addition to bulk regularity
we impose that the v.e.v. decays at large $x_{\perp}$ in Minkowski space so
that (34) is valid.

We readily see that in Minkowski space (34) implies that at late proper time
$\tau\sim\infty$ (with $x=qx_{\perp}$ and $t=q\tau$):

$\displaystyle
O(x_{\perp},\tau)=q^{4}L^{4}\Bigg{(}\frac{16(81o^{a}_{0}+o^{b}_{0})}{81t^{8}}$
$\displaystyle+\frac{64(729o^{a}_{1}+o^{b}_{1}+729o^{a}_{0}(-1+x^{2})+3o^{b}_{0}(-1+3x^{2}))}{729t^{10}}+\mathcal{O}(t^{-12})\Bigg{)}.$
(38)

and at early proper time $\tau\sim 0$ (with $x=qx_{\perp}$ and $t=q\tau$):

$\displaystyle
O(x_{\perp},\tau)=q^{4}L^{4}\Bigg{(}\frac{16(o^{a}_{0}+o^{b}_{0})}{(x^{2}+1)^{4}}$
$\displaystyle+\frac{64(o^{a}_{1}+o^{b}_{1}+o^{a}_{0}(-1+x^{2})+o^{b}_{0}(-3+x^{2}))}{(x^{2}+1)^{6}}t^{2}+\mathcal{O}(t^{3})\Bigg{)}.$
(39)

The latter is similar to (36) and we see that we get finite v.e.v. in the
limit $\tau\rightarrow 0$. It also follows that from (34) that we can
understand the large $x_{\perp}$ behavior of the v.e.v. in Minkowski space via
the general form (with $x=qx_{\perp}$ and $t=q\tau$):

$\displaystyle
O(x_{\perp},\tau)=q^{4}L^{4}\Bigg{(}\frac{16(o^{a}_{0}+o^{b}_{0})}{x^{8}}+\frac{64(o^{a}_{0}+o^{b}_{0})(-1+t^{2})}{x^{10}}$
$\displaystyle+\frac{32(5(o^{a}_{0}+o^{b}_{0})+2(-7o^{a}_{0}-9o^{b}_{0}+o^{a}_{1}+o^{b}_{1})t^{2}+5(o^{a}_{0}+o^{b}_{0})t^{4})}{x^{12}}$
$\displaystyle+\frac{64(-1+t^{2})(5(o^{a}_{0}+o^{b}_{0})+2(-11o^{a}_{0}-17o^{b}_{0}+3o^{a}_{1}+3o^{b}_{1})t^{2}+5(o^{a}_{0}+o^{b}_{0})t^{4})}{x^{14}}$
$\displaystyle+\mathcal{O}(x^{-16})\Bigg{)}.$ (40)

We note that both the small $\tau$ and large $x_{\perp}$ behavior of the
v.e.v. in Minkowski space is determined systematically by (34).

In order to obtain the bulk solution of the scalar field which yields the
v.e.v. which behaves as (34) in de-Sitter space, we proceed with the following
anstaz for the massless bulk scalar field:

$\displaystyle\Phi(r,\rho)=\sum_{n=0}^{\infty}\frac{\phi^{a}_{n}(v)}{\cosh{(\rho/L)}^{4+2n}}+\sum_{n=0}^{\infty}\frac{\phi^{b}_{n}(v)}{\left(2\cosh{(\rho/L)}+\sinh{(\rho/L)}\right)^{4+2n}}$
(41)

With the above ansatz we can solve the massless Klein-Gordon equation in the
background (locally) pure $AdS_{5}$ metric in the late time expansion
($\sigma\rightarrow\infty$) in powers of $\sigma^{-1}$ as in the previous
section, and also in the early time expansion ($\sigma\rightarrow 0$) in
powers of $\sigma$. At $n$-th order in these expansions, we obtained two
coupled second order ordinary differential equations for $\phi^{a}_{n}(v)$ and
$\phi^{b}_{n}(v)$. The general solutions of these equations has four
integration constants. We impose that

  * •

the bulk scalar field is normalizable so that
$\phi^{a}_{n}(0)=\phi^{b}_{n}(0)=0$, and

  * •

the solution of the bulk scalar field is regular at late time
($\sigma\rightarrow\infty$) at the future horizon $v=1$.

These give three conditions so that instead of four we have only one
integration constant at each order which we denote as $\Gamma_{n}$.
Remarkably, we find that imposing that the solution is regular at late time
$\sigma\rightarrow\infty$ automatically implies that at each order the
solution is regular at all time and even at very early time $\sigma\rightarrow
0$. We recall that this feature is absent for the exact solution (24) and
(25). Explicitly, the above procedure yields that

$\displaystyle\phi^{a}_{0}(v)=\Gamma_{0}\frac{v^{4}}{\left(1+v\right)^{4}},\,\,\,\,\,\,\,\,\,\,\phi^{b}_{0}(v)=-\Gamma_{0}\frac{v^{4}}{\left(1+v\right)^{4}}$
$\displaystyle\phi^{a}_{1}(v)=\frac{v^{4}\left(2\Gamma_{0}\left(v^{2}+728v-1085\right)+729\Gamma_{1}\left(v^{2}+3\right)\right)}{728(v+1)^{6}},$
$\displaystyle\phi^{b}_{1}(v)=-\frac{3v^{4}\left(2\Gamma_{0}\left(243v^{2}+728v-119\right)+243\Gamma_{1}\left(v^{2}+3\right)\right)}{728(v+1)^{6}},\,\,{\rm
etc.}$ (42)

The above results can be understood in remarkably simple terms. At late time,
we recover the expansion (18) with $\alpha=4$ and coefficients $\phi_{n}(v)$
given by (4) (see more below) from (41). At early time, we obtain from (41)
the expansion

$\Phi(r,\rho)=\sigma^{4}\sum_{n=0}^{\infty}\phi^{e}_{n}(v)\sigma^{2n}$ (43)

with

$\phi^{e}_{0}(v)=16\left(\phi^{a}_{0}(v)+\phi^{b}_{0}(v)\right),\quad\phi^{e}_{1}(v)=64\left(\phi^{a}_{1}(v)+\phi^{b}_{1}(v)-\phi^{a}_{0}(v)-3\phi^{b}_{0}(v)\right),\,\,{\rm
etc.}$ (44)

We can check that at each order $\phi^{b}_{n}(v)$ is related to
$\phi^{a}_{n}(v)$ and lower order terms such that

$\phi^{e}_{n}(v)=0,\,\,{\rm for}\,\,n=0,1,2,\cdots,$ (45)

i.e.

$\lim_{\sigma\rightarrow
0^{+}}\sigma^{-n}\Phi(r,\rho)=\lim_{\rho\rightarrow-\infty}\exp(-n(\rho/L))\Phi(r,\rho)=0,\,\,\,\,{\rm
for}\,\,n=0,1,2,\cdots.$ (46)

We have checked the above to very high orders. Thus, we can conclude the late
time expansion determines the entire evolution including the initial
conditions when we assume that the v.e.v. decays at large transverse distance
from the central axis in Minkowski space and further impose bulk regularity
which necessitates the ansatz (41). Note that at any finite value of $\rho$
(i.e. $0<\sigma<\infty$), we obtain a non-trivial profile of the bulk scalar
field.

The v.e.v. can be be readily extracted from the leading $v^{4}$ term of the
Taylor expansion about the boundary at $v=0$ of $\phi^{a}_{n}$ and
$\phi^{b}_{n}$ in (5), and it takes the form (34) with

$\displaystyle o^{a}_{0}=\Gamma_{0},\,\,\,\,\,\,\,o^{b}_{0}=-\Gamma_{0}$
$\displaystyle
o^{a}_{1}=\frac{1}{728}(2187\Gamma_{1}-2170\Gamma_{0}),\,\,\,\,\,\,\,o^{b}_{1}=-\frac{3(729\Gamma_{1}-238\Gamma_{0})}{728},\,\,{\rm
etc.}$ (47)

Using (35), we can extract the coefficients of the general late time expansion
(23) which turn out to be

$\displaystyle\hat{o}_{0}=\frac{1280\Gamma_{0}}{81},\quad\hat{o}_{1}=192\Gamma_{1}-\frac{20608\Gamma_{0}}{81},\,\,{\rm
etc.}$ (48)

We can similarly extract the coefficients of the general early time expansion
(23) using (37). From (45), it is obvious that we should get

$\displaystyle\hat{o}^{e}_{n}=0,\,\,\,\,{\rm for}\,\,n=0,1,2,\cdots.$ (49)

i.e.

$\displaystyle\lim_{\sigma\rightarrow
0^{+}}\sigma^{-n}\hat{O}(\sigma)=\lim_{\rho\rightarrow-\infty}\exp(-n(\rho/L))\hat{O}(\rho)=0,\,\,\,\,{\rm
for}\,\,n=0,1,2,\cdots.$ (50)

This has further remarkable consequences. To see these, we note from (37) that
this feature implies specific relations between $o^{b}_{n}$ and $o^{a}_{n}$,
such as

$o^{b}_{0}=-o^{a}_{0},\quad o^{b}_{1}=-o^{a}_{1}+o^{a}_{0}+3o^{b}_{0},\,\,{\rm
etc.}$ (51)

Substituting these in the early $\tau$ and large $x_{\perp}$ expansions of the
v.e.v. in Minkowski space given by (5) and (5), respectively, we find that

$\displaystyle\lim_{\tau\rightarrow
0^{+}}\tau^{-n}O(\tau,x_{\perp})=\lim_{x_{\perp}\rightarrow\infty}x_{\perp}^{n}O(\tau,x_{\perp})=0,\,\,\,\,{\rm
for}\,\,n=0,1,2,\cdots,$ (52)

implying that the v.e.v. vanishes as we go back to $\tau=0$ and at large
transverse distance from the central axis faster than any positive power of
$\tau$ and any negative power of $x_{\perp}$, respectively.

To summarize, we solved for the bulk scalar field imposing bulk regularity and
the requirement that the v.e.v. of the dual operator should decay at large
transverse distance from the central axis in Minkowski space. We find that
these requirements imply that the remarkable feature (49) that the v.e.v. and
all its derivatives vanish at early de-Sitter time. This further implies that
(60) must also hold. Unfortunately, we have not been able to find any explicit
exact solution to all orders by choosing integration constants such that (49)
is satisfied. So it is not possible to say anything about the radial profile
of the bulk scalar field at a finite value of $\rho$ because the forms (34)
(and (41)) of the v.e.v. (and the bulk scalar field) give systematic
expansions only at large $|\rho|$. We will come back to this issue later.

We can repeat this entire analysis for pure gravity. We assume the following
form for $\hat{\epsilon}(\rho)$

$\displaystyle\hat{\varepsilon}(\rho)=\sum_{k=0}^{\infty}\left(\frac{\varepsilon_{k}}{\cosh{(\rho/L)}^{4+2k}}+\frac{\tilde{\varepsilon}_{k}}{\left(2\cosh{(\rho/L)}+\sinh{(\rho/L)}\right)^{4+2k}}\right),$
(53)

which is similar to (34) so that the energy density decays at large
$x_{\perp}$ in Minkowski space and it is compatible with the late de-Sitter
time expansion given by (29). This implies that the functions $A$, $B$ and $C$
in the bulk metric (3) should take the forms

$\displaystyle
A(r,\rho)=\sum_{n=0}\left(\frac{a^{a}_{n}\left(v\right)}{\cosh^{2n+4}\left(\rho/L\right)}+\frac{a^{b}_{n}\left(v\right)}{\left(\sinh\left(\rho/L\right)+2\cosh\left(\rho/L\right)\right)^{2n+4}}\right),$
$\displaystyle
B(r,\rho)=\sum_{n=0}\left(\frac{b^{a}_{n}\left(v\right)}{\cosh^{2n+4}\left(\rho/L\right)}+\frac{b^{b}_{n}\left(v\right)}{\left(\sinh\left(\rho/L\right)+2\cosh\left(\rho/L\right)\right)^{2n+4}}\right),$
$\displaystyle
C(r,\rho)=\sum_{n=0}\left(\frac{c^{a}_{n}\left(v\right)}{\cosh^{2n+4}\left(\rho/L\right)}+\frac{c^{b}_{n}\left(v\right)}{\left(\sinh\left(\rho/L\right)+2\cosh\left(\rho/L\right)\right)^{2n+4}}\right).$
(54)

The above form yield systematic expansions at $\rho\rightarrow\pm\infty$. In
these limits, we can solve Einstein’s equations systematically obtaining
coupled second order ordinary differential equations at each order. Imposing
normalizability and regularity at late time, we obtain unique normalizable
solutions with only one integration constant at each order (after fixing the
residual gauge as detailed in C). Exactly, like in the case of the bulk scalar
imposing bulk regularity at late time automatically implies the same at all
time including $\rho\rightarrow-\infty$. Explicitly,

$\displaystyle
a^{a}_{0}(v)=\frac{\gamma_{0}v^{4}}{80(v+1)^{2}},\,\,\,\,\,\,a^{b}_{0}(v)=-\frac{\gamma_{0}v^{4}}{80(v+1)^{2}},\,\,\,\,\,\,b^{a}_{0}=-\frac{\gamma_{0}v^{4}}{160(v+1)^{4}},\,\,\,\,\,\,b^{b}_{0}(v)=\frac{\gamma_{0}v^{4}}{160(v+1)^{4}},$
$\displaystyle
c^{a}_{0}(v)=\frac{\gamma_{0}v^{5}}{200(v+1)^{4}},\,\,\,\,\,\,c^{b}_{0}(v)=-\frac{\gamma_{0}v^{5}}{200(v+1)^{4}},\,$
$\displaystyle
a^{a}_{1}(v)=\frac{\gamma_{0}\left(2533v^{2}+3650v+10\right)v^{4}}{291200(v+1)^{4}}+\frac{\gamma_{1}\left(9v^{2}+10v+10\right)v^{4}}{7280(v+1)^{4}},\,\,\,\,\,\,$
$\displaystyle
a^{b}_{1}(v)=-\frac{3\gamma_{0}\left(3271v^{2}+6070v+2430\right)v^{4}}{291200(v+1)^{4}}-\frac{\gamma_{1}\left(9v^{2}+10v+10\right)v^{4}}{7280(v+1)^{4}},\,\,\,\,\,\,$
$\displaystyle
b^{a}_{1}(v)=-\frac{\gamma_{0}\left(3979v^{2}+7292v-35\right)v^{4}}{582400(v+1)^{6}}-\frac{\gamma_{1}\left(7v^{2}-4v+25\right)v^{4}}{14560(v+1)^{6}},\,\,\,\,\,\,$
$\displaystyle
b^{b}_{1}(v)=\frac{9\gamma_{0}\left(1251v^{2}+2428v+805\right)v^{4}}{582400(v+1)^{6}}+\frac{\gamma_{1}\left(7v^{2}-4v+25\right)v^{4}}{14560(v+1)^{6}},$
$\displaystyle
c^{a}_{1}(v)=\frac{\gamma_{0}\left(4352v^{2}+5131v-2569\right)v^{5}}{1019200(v+1)^{6}}+\frac{\gamma_{1}\left(8v^{2}+7v+35\right)v^{5}}{25480(v+1)^{6}},\,\,$
$\displaystyle
c^{b}_{1}(v)=-\frac{9\gamma_{0}\left(1616v^{2}+2835v+847\right)v^{5}}{1019200(v+1)^{6}}-\frac{\gamma_{1}\left(8v^{2}+7v+35\right)v^{5}}{25480(v+1)^{6}},\,\,\,\,\,\,{\rm
etc.}$ (55)

which reproduces the functions in (30), (31) and (32) at late de-Sitter time.
At early time, we obtain exactly like in the case of the scalar field that

$\lim_{\sigma\rightarrow 0^{+}}\sigma^{-n}A(r,\sigma)=\lim_{\sigma\rightarrow
0^{+}}\sigma^{-n}B(r,\sigma)=\lim_{\sigma\rightarrow
0^{+}}\sigma^{-n}C(r,\sigma)=0,\,\,{\rm for}\,\,n=0,1,2,\cdots.$ (56)

From the above, we readily find that

$\displaystyle\varepsilon_{0}=\gamma_{0}/80,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\tilde{\varepsilon}_{0}=-\gamma_{0}/80$
$\displaystyle\varepsilon_{1}=\frac{\gamma_{0}+40\gamma_{1}}{29120},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\tilde{\varepsilon}_{1}=-\frac{729\gamma_{0}+40\gamma_{1}}{29120},\,\,{\rm
etc.}$ (57)

The coefficients of the late-time expansion (29) turn out to be

$\displaystyle\hat{e}_{0}=\frac{16\gamma_{0}}{81},\quad\hat{e}_{1}=-\frac{8\left(363\gamma_{0}-40\gamma_{1}\right)}{3645},\,\,{\rm
etc.}$ (58)

As in the case of the scalar field we obtain that

$\displaystyle\lim_{\sigma\rightarrow
0}\sigma^{-n}\hat{\varepsilon}(\sigma)=0,\,\,\,\,{\rm for}\,\,n=0,1,2,\cdots.$
(59)

which implies that in Minkowski space

$\displaystyle\lim_{\tau\rightarrow
0^{+}}\tau^{-n}\varepsilon(\tau,x_{\perp})=\lim_{x_{\perp}\rightarrow\infty}x_{\perp}^{n}\varepsilon(\tau,x_{\perp})=0,\,\,\,\,{\rm
for}\,\,n=0,1,2,\cdots.$ (60)

To summarize, the requirements of regularity of the bulk metric (at late time)
and that the dual energy density decays at large $x_{\perp}$ in Minkowski
space imply that the energy density should satisfy the remarkable feature that
it and all its derivatives vanish at early de-Sitter time. The latter feature
further implies that the energy density vanishes in Minkowski space both in
the limits $\tau\rightarrow 0^{+}$ and $x_{\perp}\rightarrow\infty$ faster
than any positive power of $\tau$ and any negative power of $x_{\perp}$,
respectively. We can also conclude that when the energy density decays at
large $x_{\perp}$ on the future wedge of Minkowski space, the late de-Sitter
time expansion itself contains the full information of the initial conditions
set at any finite value of $\rho$ or $\tau$ since the coefficients
$\tilde{\varepsilon}_{n}$ in (53) get related to $\varepsilon_{n}$, and
therefore there is no independent information needed from early time to
describe the energy density on the entire future wedge.

The feature of the Gubser flow that we can deduce information about the
initial conditions by constraining that the energy density (and expectation
values of other operators) should decay at large transverse distance from the
central axis is certainly distinctive as such an extrapolation is not possible
for the Bjorken flow in which there is homogeneity in the transverse
directions. Essentially, this feature is a consequence of the symmetries of
the Gubser flow which ties the $x_{\perp}$ dependence to $\tau$ dependence as
both arise via the dependence on the de-Sitter time $\rho$.

We note that the form (53) is compatible with the hydrodynamic expansion (12)
at $\rho\sim 0$. However, this form of the energy density and also the form
(5) of the functions in the bulk metric admit systematic expansions only at
large $|\rho|$, and therefore we cannot systematically match at $\rho=0$ or
any finite value of $\rho$. The same can also be said about any resummation of
the late-time expansion. For instance, one can readily note that for the exact
solutions of the massless scalar $\Phi(v,\rho)$ given by (24) and (25), the
late time expansion in $\sigma^{-1}$ converges only for
$\sigma^{2}>(1-v)/(1+v)$ for any fixed $v$. Thus, matching with hydrodynamic
evolution, which can be valid only around $\rho\sim 0$ (i.e. $\sigma\sim 1$),
is difficult. Irrespective of whether hydrodynamic evolution is realized
around $\rho\sim 0$, it is difficult to determine the nature of the metric at
a finite value of $\rho$ because unlike the case of the scalar field, it is
apriori even difficult to estimate where the late and early time expansions
converge.

Therefore, only by solving Einstein’s equations numerically (say following the
method of characteristics developed in Chesler:2010bi ), we would be able to
determine which initial conditions set at a finite value of $\rho$ or $\tau$
would be realizable physically, i.e. those in which the energy density does
decay at large $x_{\perp}$. We do not expect that for all such initial
conditions the hydrodynamic expansion (12) would actually be realized at
$\rho\sim 0$. Numerical simulations are also necessary to establish
convergence properties of the late time and early time expansions obtained
from (53) and contrast it with the trans-series forms discussed in
Behtash:2017wqg ; Denicol:2018pak ; Behtash:2019qtk ; Dash:2020zqx ;
Soloviev:2021lhs .

##  6 Discussion

We have provided some exact results on how relativistic quantum matter evolves
in Gubser flow in holographic strongly coupled large $N$ conformal gauge
theories. The most generic behaviour is yet to be established and further
progress can be made by employing numerical relativity in this context.
Numerical analyses are needed to confirm many of the extrapolations which we
have made for the behavior on the entire future wedge in Minkowski space from
the late de-Sitter time behavior. However, the late de-Sitter time behavior
can be rigorously established just from bulk regularity alone.

The main implications of our results is that in flows with such special
symmetries, a relativistic quantum system can evolve out of hydrodynamics into
a new regime that is independent of the initial conditions. Such evolution can
therefore reveal many fundamental features of the underlying microscopic
theory. Particularly, we show that when the symmetries of the Gubser flow are
preserved, the late de-Sitter time evolution can lead to a color glass
condensate-like phase which is independent of the initial conditions. Of
course, this phase is very different from the perfect fluid state which is
realized for generic initial conditions which violate the symmetries of the
Gubser flow.

Another implication of our results is that the initial conditions which
realize such flows with special symmetries are further constrained beyond
these symmetries. We have shown that in Gubser flows in holographic strongly
coupled large $N$ conformal gauge theories, the energy density should vanish
faster than any power of the proper time $\tau$ as we go back to $\tau=0$, the
boundary of the future wedge. It is very unlikely that such a behavior can be
realized by colliding gravitational shock waves (which represent collision of
energy lumps in the dual gauge theory) as studied in Chesler:2010bi . The bulk
radial profiles of these shocks should be supported sufficiently far away from
the boundary so that the dual energy density (given by the fourth radial
derivative of the $A$ function in the bulk metric (3)) and its proper time
derivatives should vanish at the moment of collisions, and this could as well
imply that the radial profiles of these shocks actually vanish. Therefore,
instead of being generated by “asymptotic states,” the Gubser flow naturally
could be smoothly glued to pure vacuum.111Since $\varepsilon(x_{\perp},\tau)$
and all its $\tau$-derivatives vanish at $\tau=0$, one can argue via a
Fefferman-Graham type expansion as well that the bulk metric at $\tau=0$
should be close to pure $AdS_{5}$. This could be another alternative way to
prove that the bulk solution can be smoothly glued to pure $AdS_{5}$ outside
the future wedge. The non-trivial behavior in the future wedge evolving back
to the vacuum arises from self-consistent pressure gradients. Of course, we
need to numerically confirm the extrapolation we have made to the entire
future wedge from the late de-Sitter time behavior.

It would be also interesting to investigate bulk solutions which behave as a
Gubser flow in the vicinity of a specific azimuthal direction but is otherwise
homogeneous in transverse directions, so that it corresponds to a Bjorken flow
approximately. Such a solution can be interpreted as a jet embedded in an
expanding medium and also created at early proper time with non-trivial
initial conditions. The existence of such solutions could indicate that Gubser
flow can be applied to collective flow arising in jets if not the flow of the
full quark gluon plasma. Gubser flow indeed has recently been successfully
applied for understanding two and four particle correlations in jets with
large charge multiplicity arising in pp-collisions Taghavi:2019mqz ;
ALICE:2023lyr  (in Taghavi:2019mqz , it has been argued that four particle
cumulants can be related to fluctuations in initial conditions).

It would be interesting to understand these issues in QCD:

  * •

Is there a generic nature of the final phase of the evolution of Gubser flow
in de-Sitter time which is independent of the initial conditions, and if so,
then what are the characteristics of this phase?

  * •

What kind of constraints on the initial conditions lead to such flow with
special symmetries like Gubser flow (beyond the necessary symmetries)?

  * •

Can such initial conditions be realized perhaps with greater likelihood in
jets and other substructures within the full bulk evolution of the system?

These questions are relevant for the phenomenology of heavy-ion collisions.
Therefore, we plan to study Gubser flow in holographic non-conformal/confining
gauge theories Gursoy:2015nza ; Chattopadhyay:2021ive ; Jaiswal:2021uvv ;
Chen:2021wwh  and in semi-holographic Iancu:2014ava ; Banerjee:2017ozx ;
Kurkela:2018dku ; Mitra:2020mei ; Mitra:2022xtb  scenarios and see how both,
the initial conditions and the late de-Sitter regime for the Gubser flow
distinguishes confining behaviour from an emergent infrared critical point.

The advantage of studying Gubser flow in a holographic theory is that we can
understand the quantum information theoretic aspects of such a scenario in
which a many-body system can escape hydrodynamization, especially the
fundamental reason why and how despite lack of hydrodynamization the evolution
can reach a phase independent of the initial conditions. Many novel aspects of
quantum thermodynamics RevModPhys.93.035008  can be understood via explicit
computations of entanglement measures (see for example Kibe:2021qjy ;
Banerjee:2022dgv ). A preliminary discussion on the entropy production
captured by the growth of the area of the horizons, and puzzles regarding its
interpretation has been presented in D.

###### Acknowledgements.

It is a pleasure to thank Jorge Casalderrey-Solana, Suresh Govindarajan, Sašo
Grozdanov, Matti Järvinen, David Mateos, Giuseppe Policastro, Alexandre
Serantes and Amitabh Virmani for discussions. We thank Sašo Grozdanov, Edmond
Iancu, Matti Järvinen, David Müller, Mike Strickland, Victor Roy and Amitabh
Virmani for comments on the manuscript. Both AB and AM acknowledge support
from IFCPAR/CEFIPRA funded project no 6304-3. Research of AB is partially
supported by the European MSCA grant HORIZONMSCA-2022-PF-01-01 and by the
H.F.R.I call “Basic research Financing (Horizontal support of all Sciences)”
under the National Recovery and Resilience Plan “Greece 2.0” funded by the
European Union– NextGenerationEU (H.F.R.I. Project Number: 15384). The
research of AM is also supported by the Center of Excellence initiative of the
Ministry of Education of India, and the new faculty seed grant of IIT Madras.
AS is supported by funding from Horizon Europe research and innovation
programme under the Marie Skłodowska-Curie grant agreement No. 101103006 and
the project N1-0245 of Slovenian Research Agency (ARIS). TM has been supported
by an appointment to the JRG Program at the APCTP through the Science and
Technology Promotion Fund and Lottery Fund of the Korean Government. TM has
also been supported by the Korean Local Governments – Gyeongsangbuk-do
Province and Pohang City – and by the National Research Foundation of Korea
(NRF) funded by the Korean government (MSIT) (grant number 2021R1A2C1010834).

##  Appendix A Viscous Hydrodynamics

The conformal viscous Gubser flow has been studied in Gubser:2010ui . The
conformal perfect fluid energy-momentum tensor is
$\hat{t}^{\mu\nu}_{0}=\varepsilon(T)(4u^{\mu}u^{\nu}+g^{\mu\nu})$, where
$\varepsilon(T)\propto T^{4}$ is determined by the equation of state in terms
of the temperature $T$. Including the first order viscous correction, we have

$\displaystyle\hat{t}^{\mu\nu}=\hat{t}^{\mu\nu}_{0}-2\eta(T)\sigma^{\mu\nu},$
(61)

where $\sigma^{\mu\nu}$ is the shear stress tensor defined as

$\displaystyle\sigma^{\mu\nu}=\frac{1}{2}\Delta^{\alpha\mu}\Delta^{\beta\nu}\left(\nabla_{\alpha}u_{\beta}+\nabla_{\beta}u_{\alpha}\right)-\frac{1}{3}\nabla_{\lambda}u^{\lambda}\Delta^{\mu\nu},$
(62)

and $\eta(T)\propto T^{3}$ is the shear viscosity. It is useful to define
$H_{0}=\eta(T)\hat{\varepsilon}^{-3/4}(T)$ which is a constant in a conformal
field theory and is determined by microscopic dynamics. Above,
$\Delta^{\mu\nu}=u^{\mu}u^{\nu}+g^{\mu\nu}$. The conservation of the energy-
momentum tensor (i.e. the hydrodynamic equations) for the Gubser flow in ${\rm
d}S_{3}\times\mathbf{R}$ frame is simply

$\displaystyle\hat{\varepsilon}^{\prime}({\tilde{\rho}})+\frac{8}{3}\tanh({\tilde{\rho}})\hat{\varepsilon}({\tilde{\rho}})-\frac{16}{9}H_{0}\tanh^{2}({\tilde{\rho}})\hat{\varepsilon}({\tilde{\rho}})^{3/4}=0,$
(63)

with ${\tilde{\rho}}=\rho/L$. The solution of this equation is

$\displaystyle\varepsilon({\tilde{\rho}})$
$\displaystyle=\frac{\text{sech}^{4}({\tilde{\rho}})}{531441}\Big{(}4H_{0}L^{-1}\sinh^{3}({\tilde{\rho}}){\cosh^{1/3}({\tilde{\rho}})}$
${}_{2}F_{1}\left(\frac{7}{6},\frac{3}{2};\frac{5}{2};-\sinh^{2}({\tilde{\rho}})\right)+27\hat{\varepsilon}_{0}^{1/4}{\cosh^{1/3}({\tilde{\rho}})}\Big{)}{}^{4}.$
(64)

The viscous energy density $\hat{\varepsilon}$ in the limit
${\tilde{\rho}}\rightarrow 0$ gives the derivative expansion (12) as mentioned
in the main text. However, in the limit ${\tilde{\rho}}\rightarrow\infty$ the
energy density has the following expansion dominated by viscous term at the
leading order:

$\displaystyle\hat{\varepsilon}({\tilde{\rho}})=\frac{16}{81}H_{0}^{4}L^{-4}+\mathcal{O}(\cosh({\tilde{\rho}})^{-2/3})+\ldots$
(65)

This indicates the breakdown of derivative expansion in the large
${\tilde{\rho}}$ limit.

##  Appendix B Holographic renormalization

For compactness, we define ${\tilde{\rho}}=\rho/L$ and also use the variable
$v$ with $v=r/L$ as in the main text. We can solve for the metric functions
$A(v,{\tilde{\rho}}),B(v,{\tilde{\rho}})$ and $C(v,{\tilde{\rho}})$ in the
radial expansion near the boundary $v=0$ as indicated in (3) and obtain,

$\displaystyle
A(v,{\tilde{\rho}})=a_{(1)}\left({\tilde{\rho}}\right)v+\frac{1}{4}\left(a_{(1)}^{2}\left({\tilde{\rho}}\right)-4a_{(1)}^{\prime}\left({\tilde{\rho}}\right)\right)v^{2}$
$\displaystyle+a_{(4)}\left({\tilde{\rho}}\right)v^{4}+\ldots,$ (66)
$\displaystyle
B(v,{\tilde{\rho}})=a_{(1)}({\tilde{\rho}})\cosh^{2}({\tilde{\rho}})~{}v$
$\displaystyle-\frac{1}{4}a_{(1)}({\tilde{\rho}})\left(a_{(1)}({\tilde{\rho}})+4\tanh({\tilde{\rho}})\right)~{}v^{2}$
$\displaystyle+\frac{1}{12}a_{(1)}({\tilde{\rho}})\Big{(}\left(a_{(1)}({\tilde{\rho}})+3\tanh\left({\tilde{\rho}}\right)\right)^{2}+3\tanh^{2}\left({\tilde{\rho}}\right)\Big{)}~{}v^{3}$
$\displaystyle+b_{(4)}({\tilde{\rho}})~{}v^{4}+\ldots,$ (67) $\displaystyle
C(v,{\tilde{\rho}})=3a_{(1)}({\tilde{\rho}})~{}v-\frac{1}{4}a_{(1)}({\tilde{\rho}})\Big{(}8\tanh({\tilde{\rho}})+3a_{(1)}({\tilde{\rho}})\Big{)}~{}v^{2}$
$\displaystyle+\frac{1}{4}a_{(1)}({\tilde{\rho}})\left(\left(a_{(1)}({\tilde{\rho}})+2\tanh({\tilde{\rho}})\right)^{2}+4\tanh^{2}({\tilde{\rho}})\right)~{}v^{3}$
$\displaystyle+s_{(4)}({\tilde{\rho}})~{}v^{4}+\ldots,$ (68)

where $a_{(1)}\left({\tilde{\rho}}\right)$ is an arbitrary function. The
functions $b_{(4)}\left({\tilde{\rho}}\right)$ and
$s_{(4)}\left({\tilde{\rho}}\right)$ are determined by
$a_{(4)}\left({\tilde{\rho}}\right)$ and $a_{(1)}\left({\tilde{\rho}}\right)$
by the constraints of Einstein’s equations:

$\displaystyle
s_{(4)}\left({\tilde{\rho}}\right)=-\frac{1}{32}a_{(1)}\left({\tilde{\rho}}\right)\Big{(}a_{(1)}^{2}\left({\tilde{\rho}}\right)+4\tanh^{2}\left({\tilde{\rho}}\right)\Big{)}$
$\displaystyle\Big{(}3a_{(1)}\left({\tilde{\rho}}\right)+16\tanh\left({\tilde{\rho}}\right)\Big{)}+\frac{9}{8}a_{(1)}^{2}\left({\tilde{\rho}}\right)\tanh^{2}\left({\tilde{\rho}}\right),$
$\displaystyle
a_{(4)}^{\prime}\left({\tilde{\rho}}\right)=\frac{8}{3}\tanh\left({\tilde{\rho}}\right)\Big{(}b_{(4)}\left({\tilde{\rho}}\right)-a_{(4)}\left({\tilde{\rho}}\right)\Big{)}$
$\displaystyle+\frac{1}{12}\Big{(}a_{(1)}^{3}\left({\tilde{\rho}}\right)+4a_{(1)}\left({\tilde{\rho}}\right)\tanh^{2}\left(\frac{\rho}{L}\right)\Big{)}$
$\displaystyle\Big{(}a_{(1)}\left({\tilde{\rho}}\right)+8\tanh\left({\tilde{\rho}}\right)\Big{)}+\frac{5}{6}a_{(1)}^{2}\left({\tilde{\rho}}\right)\tanh^{2}\left({\tilde{\rho}}\right).$
(69)

We find that the radial expansions and thus $A$, $B$ and $C$ are entirely
determined by the two functions $a_{(4)}\left({\tilde{\rho}}\right)$ and
$a_{(1)}\left({\tilde{\rho}}\right)$.

However, $a_{(1)}\left({\tilde{\rho}}\right)$ is a residual gauge freedom and
can be set to zero by the diffeomorphism $v\rightarrow v+f({\tilde{\rho}})$
which preserves the ingoing Eddington-Finkelstein gauge. Furthermore, the
latter is a proper diffeomorphism, meaning that
$a_{(1)}\left({\tilde{\rho}}\right)$ does not affect boundary data. Indeed
extracting the dual energy-momentum tensor $\langle T^{\mu}_{\,\nu}\rangle$ at
the boundary from the renormalized on-shell gravitational action following
Henningson:1998gx ; Balasubramanian:1999re ; deHaro:2000vlm  and using (B), we
find that $a_{(1)}\left({\tilde{\rho}}\right)$ disappears. Explicitly, we
obtain

$\displaystyle\langle
T^{\mu}_{\,\nu}\rangle=\begin{pmatrix}\hat{t}^{\rho}_{\,\rho}&0&0&0\\\
0&\hat{t}^{\theta}_{\,\theta}&0&0\\\ 0&0&\hat{t}^{\phi}_{\,\phi}&0\\\
0&0&0&\hat{t}^{\eta}_{\,\eta}\end{pmatrix}+\mathcal{A}^{\mu}_{\,\nu}$ (70)

where

$\displaystyle\hat{t}^{\rho}_{\,\rho}\equiv-\hat{\varepsilon}\left({\tilde{\rho}}\right)=\frac{3l^{3}}{16\pi
G_{N}}a_{(4)}\left({\tilde{\rho}}\right),$
$\displaystyle\hat{t}^{\theta}_{\,\theta}\equiv\hat{P_{T}}\left({\tilde{\rho}}\right)=\frac{3l^{3}}{16\pi
G_{N}}\left(\frac{1}{2}\coth\left({\tilde{\rho}}\right)a_{(4)}^{\prime}\left({\tilde{\rho}}\right)+a_{(4)}\left({\tilde{\rho}}\right)\right),$
$\displaystyle\hat{t}^{\phi}_{\,\phi}\equiv\hat{P_{T}}\left({\tilde{\rho}}\right)=\hat{t}^{\theta}_{\,\theta},$
$\displaystyle\hat{t}^{\eta}_{\,\eta}\equiv\hat{P_{L}}\left({\tilde{\rho}}\right)=-\frac{3l^{3}}{16\pi
G_{N}}\Big{(}\coth\left({\tilde{\rho}}\right)a_{(4)}^{\prime}\left({\tilde{\rho}}\right)+3a_{(4)}\left({\tilde{\rho}}\right)\Big{)}.$
(71)

and $\mathcal{A}^{\mu}_{\,\nu}$ is the Weyl anomaly given by

$\displaystyle\hat{\mathcal{A}}^{\mu}_{\,\nu}$ $\displaystyle=$
$\displaystyle\frac{l^{3}}{128\pi
G_{N}}\Big{(}\frac{4}{3}R^{\mu}_{\,\nu}R-2R^{\mu}_{\,\kappa}R^{\kappa}_{\,\nu}-\delta^{\mu}_{\,\nu}\Big{(}\frac{1}{2}R^{2}-R_{\kappa\sigma}R^{\kappa\sigma}\Big{)}\Big{)}$
(72) $\displaystyle=$ $\displaystyle\frac{l^{3}}{64\pi
G_{N}}diag\left(1,1,1,-3\right)L^{-4},$

where the curvatures refer to those of the background metric (4) for the dual
theory (which is also the boundary metric of the five-dimensional bulk
spacetime). With the identification $l^{3}/G_{N}=2N^{2}/\pi$, we obtain (10).

##  Appendix C Finding values of $\alpha$ for pure gravity

Here we give details of how $\alpha$ has been determined replicating the
strategy for the massless scalar field discussed in the main text. The ansatz
(26) is a double expansion in $\sigma$. Assuming analytic behavior at the
future horizon $r=L$ (i.e. $v=1$), we obtain the near horizon expansion (with
numerical coefficients $a_{nm,p}$)

$\displaystyle A_{nm}(v)=\sum_{p=0}^{\infty}a_{nm,p}(1-v)^{p},$ (73)

and similarly for $B_{nm}(v)$ and $C_{nm}(v)$ in terms of the numerical
coefficients $b_{nm,p}$ and $c_{nm,p}$, respectively. Solving the
gravitational equations in this near-horizon expansion, we find that
$A_{10}(v)$, $B_{10}(v)$ and $C_{10}(v)$ are determined by the three
integration constants $a_{10,0}$, $b_{10,0}$ and $c_{10,0}$. Explicitly,

$\displaystyle
A_{10}(v)=a_{10,0}\left(1-(\alpha+1)(1-v)+\mathcal{O}((1-v)^{2})\ldots\right)$
$\displaystyle+b_{10,0}\left(-2\alpha(1-v)+\mathcal{O}((1-v)^{2})+\ldots\right)$
$\displaystyle+c_{10,0}\left(\frac{1}{2}\alpha(1+\alpha)(1-v)+\mathcal{O}((1-v)^{2})+\ldots\right),$
(74)

and also similarly for $B_{10}(v)$ and $C_{10}(v)$.

Normalizable solutions can be matched to the near-boundary expansion (B)-(B).
However, the solution of $A_{10}(v)$ given by (C) when expanded near the
boundary $v=0$ reads as

$\displaystyle A_{10}(v)$
$\displaystyle=P(\alpha)(a_{10,0}+(1+\alpha)(2b_{10,0}-c_{10,0}))$ (75)
$\displaystyle+\mathcal{O}(v)+\ldots$

where $P(\alpha)$ is polynomial in $\alpha$. Therefore, for normalizability,
$\alpha$ should be the roots of $P(\alpha)$ and these are $\alpha=2k$, with
$k=0,1,2,\ldots$. This can be found by performing the expansion (C) to high
orders and then re-expanding it as a Taylor series about $v=0$ as in the case
of the massless scalar.

For $k=0$, $A_{10}$, $B_{10}$ and $C_{10}$ vanish. For the case $k=1$ we have
been able to find a solution which corresponds to the boundary expansion
(B)-(B) with $a_{(1)}({\tilde{\rho}})\neq 0$ but $a_{(4)}({\tilde{\rho}})=0$.
Since $a_{(1)}({\tilde{\rho}})$ corresponds to a pure residual (proper) gauge
freedom, this solution for the case $k=1$ yields just a pure gauge
deformation. However, we cannot rule out this case completely as another
physical solution may exist. For $k=2,3,\ldots$, we are able to find physical
solutions as reported in the main text.

To see how we reach these conclusions, it is useful to use the residual gauge
freedom to set $a_{(1)}({\tilde{\rho}})=0$. In this case, the near boundary
expansions (B)-(B) start from $v^{4}$. For the comparison with near-horizon
expansions for $A_{10}$, $B_{10}$ and $C_{10}$ we should further impose a
late-time expansion of $a_{4}({\tilde{\rho}})$ of the form

$a_{(4)}({\tilde{\rho}})=\sum_{n=1}^{\infty}\sum_{m=0}^{\infty}\sigma^{-n\alpha-2m}a_{(4)nm}$
(76)

that follows from (26) with constant coefficients $a_{(4)nm}$. Setting the
coefficients of $v$ and $v^{2}$ to zero in (75), we get two relations between
the three integration constants $a_{10,0}$, $b_{10,0}$ and $c_{10,0}$. For
$k=0$, these simply set $a_{10.0}=0$ and then requiring that the boundary
values of $B_{10}$ and $C_{10}$ to vanish (so that we get normalizable
solutions) we obtain $b_{10.0}=c_{10.0}=0$. Then $A_{10}$, $B_{10}$ and
$C_{10}$ have to vanish. For $k\geq 2$, setting the coefficients of $v$ and
$v^{2}$ to zero in (75) give two relations to determine $b_{10,0}$ and
$c_{10,0}$ in terms of $a_{10,0}$. Using these values of $b_{10,0}$ and
$c_{10,0}$, we find a perfect agreement with the near boundary expansions
(B)-(B) corresponding to normalizable solutions where the residual gauge
freedom has been fixed with $a_{(1)}=0$. The remaining integration constant
$a_{10,0}$ simply determines the leading term in $a_{(4)}({\tilde{\rho}})$ in
(76), i.e. $a_{(4)10}$. For the allowed values of $\alpha$ which are $4+2n$,
with $n=0,1,2,\ldots$ the ansatz (26) can be simplified to (28), as mentioned
before. (Then $a_{(4)10}$ identified with $\hat{e}_{0}$ up to a numerical
constant, etc.) The case $k=1$ is somewhat tricky because the series expansion
of the equations near the horizon themselves do not have unique solutions. One
way of solving it leads to a pure gauge solution. We have not been able to
show that physical solutions do not exist for this case.

At higher orders in the late de-Sitter time expansion in case of $\alpha=4+2n$
and $n=0,1,2,\ldots$, we obtain normalizable solutions because the lower
orders source only normalizable particular solutions, while the homogeneous
solutions are also normalizable with a new arbitrary integration constant at
each order since $\alpha=6,8,\ldots$ are allowed to be the leading order
behavior also (when the lower order coefficients vanish). This argument is
similar to the case of the massless bulk scalar. We thus establish that (29)
gives realizable late de-Sitter time expansion for the dual energy density.

##  Appendix D Examining the entropy

Let us first examine the entropy $S$ in the de-Sitter frame. The entropy $S$
of the black hole given by (3) is given by the area $\Sigma_{H}$ of the
apparent or event horizon located at $v=v_{h}({\tilde{\rho}})$

$\displaystyle S=\frac{\Sigma_{H}}{4G_{N}}$ (77)

where $G_{N}$ is the Newton constant. Here we will examine the event horizon.
The area can be computed via
$\Sigma_{H}=\int_{v=v_{h}}\sqrt{\gamma}~{}d\theta~{}d\phi~{}d\eta$ with
$\gamma$ being the induced metric on the spatial sections of the event horizon
(generated by a congruence of null geodesics). Explicitly,

$\displaystyle\gamma=l^{6}\frac{e^{C(v_{h}({\tilde{\rho}}),{\tilde{\rho}})}\sin^{2}\theta\left(\cosh\left({\tilde{\rho}}\right)+v_{h}({\tilde{\rho}})\sinh\left({\tilde{\rho}}\right)\right)^{4}}{v_{h}^{6}({\tilde{\rho}})}.$
(78)

Here $\theta\in[0,\pi]$, $\phi\in[0,2\pi]$ and $\eta\in[-\infty,\infty]$. Due
to the infinite extent of $\eta$, it is better to define the entropy per unit
rapidity which is

$\displaystyle\frac{dS}{d\eta}=\pi\frac{l^{3}}{G_{N}}\frac{e^{\frac{C(v_{h}({\tilde{\rho}}),{\tilde{\rho}})}{2}}\left(\cosh\left({\tilde{\rho}}\right)+v_{h}({\tilde{\rho}})\sinh\left({\tilde{\rho}}\right)\right)^{2}}{v_{h}^{3}({\tilde{\rho}})}.$
(79)

The corresponding entropy density per unit rapidity reads

$\displaystyle\frac{ds}{d\eta}=\frac{1}{4\pi
L^{3}\cosh^{2}\left({\tilde{\rho}}\right)}\frac{dS}{d\eta}.$ (80)

The factor of $4\pi\cosh^{2}\left({\tilde{\rho}}\right)$ is the area of the
sphere at the boundary (measured by the background metric (4)). The location
of the event horizon $v_{h}({\tilde{\rho}})$ can be determined by the radial
null geodesic equation, i.e.

$\displaystyle\frac{{\rm d}v_{h}({\tilde{\rho}})}{{\rm
d}{\tilde{\rho}}}+\frac{1-v_{h}^{2}({\tilde{\rho}})+A(v_{h}({\tilde{\rho}}),{\tilde{\rho}})}{2}=0$
(81)

with the condition that $v_{h}({\tilde{\rho}}=\infty)=1$, since the horizon
coincides with that of the solution dual to the vacuum in the limit
${\tilde{\rho}}\rightarrow\infty$.

In the state dual to the vacuum, we simply have

$\displaystyle\frac{ds}{d\eta}=\frac{l^{3}}{4L^{3}\cosh^{2}\left({\tilde{\rho}}\right)G_{N}}\left(\cosh\left({\tilde{\rho}}\right)+\sinh\left({\tilde{\rho}}\right)\right)^{2},$
$\displaystyle=\frac{l^{3}}{L^{3}G_{N}}\times\frac{1}{4}\left(1+\tanh\left({\tilde{\rho}}\right)\right)^{2}.$
(82)

This implies that the entropy density per unit rapidity monotonically
increases from zero to a constant value with the de-Sitter time as shown in
Fig. 2.

Figure 2: Vacuum entropy density per unit rapidity plotted as a function of
${\tilde{\rho}}$. The figure shows that at late-time $ds/d\eta$ goes to a
constant.

At this point, it would seem strange to associate an entropy to the vacuum
state in the Gubser flow both from the de-Sitter and Minkowski future wedge
points of view. In the latter case, this is the entropy associated with an
accelerated observer like the Rindler observer. The appearance of entropy in
the de-Sitter vacuum is more subtle and has been discussed recently in
Chandrasekaran:2022cip . Although the bulk metric dual to the vacuum in de-
Sitter space is locally $AdS_{5}$, note that bulk diffeomorphism can produce a
non-trivial entropy as in the case of the map from pure three-dimensional
anti-de-Sitter space ($AdS_{3}$) to the Banados Teitelboim Zanelli (BTZ) black
hole Banados:1992wn  (however the physical interpretation of the entropy in
our case is more similar to that in Emparan:1999gf .)

The bulk dual of aGubser flow admits the late-time expansion (28), and thus
the location of the event horizon which is a solution of (81) with the
condition $v_{h}({\tilde{\rho}}=\infty)=1$, admits the following expansion

$\displaystyle
v_{h}({\tilde{\rho}})=1+\sum_{m=0}^{\infty}v_{h,m}e^{-(4+2m){\tilde{\rho}}}$
(83)

at late time. Using the explicit perturbative solution for the metric given by
(30)-(32), we obtain

$\displaystyle v_{h,0}=\frac{\hat{e}_{0}}{40},\quad
v_{h,1}=\frac{24\hat{e}_{0}+29\hat{e}_{1}}{2240},\quad\ldots$ (84)

and

$\displaystyle
L^{3}\frac{G_{N}}{l^{3}}\frac{ds}{d\eta}=1-2e^{-2{\tilde{\rho}}}+3\left(1-\frac{\hat{e}_{0}}{80}\right)e^{-4{\tilde{\rho}}}+\ldots$
(85)

In the case of the vacuum with $0=\hat{e}_{0}=\ldots$, the above series has
alternating signs. However, the function is monotonically increasing (and
assuming the value $1$ in the limit ${\tilde{\rho}}\rightarrow\infty$) as
evident from the exact expression (D) (that is plotted in Fig. 2). The
monotonic growth of the entropy density per unit rapidity is not affected
perturbatively and it reaches the constant vacuum value at late time. It has
been shown in Maldacena:2012xp  that the entanglement entropy can scale as the
volume in de-Sitter space with a numerical factor which has a maximal value.
If $\hat{e}_{0}$ is positive, then the energy density at leading order is
negative, and $ds/d\eta$ is less than in the vacuum. On the other hand if
$\hat{e}_{0}$ is negative, then the energy density at leading order is
positive, and $ds/d\eta$ is larger than that in the vacuum. Here, it is also
useful to note that it has been argued that the entropy of an excited state in
de-Sitter space is expected to be less than that of the vacuum
Chandrasekaran:2022cip . It will be useful to understand the case of the
holographic Gubser flow better with a clear algebraic interpretation of
$ds/d\eta$ computed here.

## References

  * (1) S.S. Gubser, Phys. Rev. D 82, 085027 (2010).  DOI 10.1103/PhysRevD.82.085027
  * (2) S.S. Gubser, A. Yarom, Nucl. Phys. B 846, 469 (2011).  DOI 10.1016/j.nuclphysb.2011.01.012
  * (3) G.S. Denicol, U.W. Heinz, M. Martinez, J. Noronha, M. Strickland, Phys. Rev. Lett. 113(20), 202301 (2014).  DOI 10.1103/PhysRevLett.113.202301
  * (4) G.S. Denicol, U.W. Heinz, M. Martinez, J. Noronha, M. Strickland, Phys. Rev. D 90(12), 125026 (2014).  DOI 10.1103/PhysRevD.90.125026
  * (5) A. Behtash, C.N. Cruz-Camacho, M. Martinez, Phys. Rev. D 97(4), 044041 (2018).  DOI 10.1103/PhysRevD.97.044041
  * (6) M. Martinez, M. McNelis, U. Heinz, Phys. Rev. C 95(5), 054907 (2017).  DOI 10.1103/PhysRevC.95.054907
  * (7) G.S. Denicol, J. Noronha, Phys. Rev. D 99(11), 116004 (2019).  DOI 10.1103/PhysRevD.99.116004
  * (8) C. Chattopadhyay, U. Heinz, S. Pal, G. Vujanovic, Nucl. Phys. A 982, 287 (2019).  DOI 10.1016/j.nuclphysa.2018.11.005
  * (9) A. Behtash, S. Kamata, M. Martinez, H. Shi, JHEP 07, 226 (2020).  DOI 10.1007/JHEP07(2020)226
  * (10) A. Dash, V. Roy, Phys. Lett. B 806, 135481 (2020).  DOI 10.1016/j.physletb.2020.135481
  * (11) B.K. Berger, Living Rev. Rel. 5, 1 (2002).  DOI 10.12942/lrr-2002-1
  * (12) T. Epelbaum, F. Gelis, Phys. Rev. Lett. 111, 232301 (2013).  DOI 10.1103/PhysRevLett.111.232301
  * (13) D. Müller, Simulations of the Glasma in 3+1D.  Ph.D. thesis, Vienna, Tech. U. (2019)
  * (14) J. Berges, M.P. Heller, A. Mazeliauskas, R. Venugopalan, Rev. Mod. Phys. 93(3), 035003 (2021).  DOI 10.1103/RevModPhys.93.035003
  * (15) M. Rangamani, Class. Quant. Grav. 26, 224003 (2009).  DOI 10.1088/0264-9381/26/22/224003
  * (16) M. Henningson, K. Skenderis, JHEP 07, 023 (1998).  DOI 10.1088/1126-6708/1998/07/023
  * (17) V. Balasubramanian, P. Kraus, Commun. Math. Phys. 208, 413 (1999).  DOI 10.1007/s002200050764
  * (18) J.M. Maldacena, Adv. Theor. Math. Phys. 2, 231 (1998).  DOI 10.4310/ATMP.1998.v2.n2.a1
  * (19) E. Witten, Adv. Theor. Math. Phys. 2, 253 (1998).  DOI 10.4310/ATMP.1998.v2.n2.a2
  * (20) S.S. Gubser, I.R. Klebanov, A.M. Polyakov, Phys. Lett. B 428, 105 (1998).  DOI 10.1016/S0370-2693(98)00377-3
  * (21) D.Z. Freedman, S.S. Gubser, K. Pilch, N.P. Warner, Adv. Theor. Math. Phys. 3, 363 (1999).  DOI 10.4310/ATMP.1999.v3.n2.a7
  * (22) I. Heemskerk, J. Polchinski, JHEP 06, 031 (2011).  DOI 10.1007/JHEP06(2011)031
  * (23) I. Bredberg, C. Keeler, V. Lysov, A. Strominger, JHEP 03, 141 (2011).  DOI 10.1007/JHEP03(2011)141
  * (24) S. Kuperstein, A. Mukhopadhyay, JHEP 11, 130 (2011).  DOI 10.1007/JHEP11(2011)130
  * (25) S. Kuperstein, A. Mukhopadhyay, JHEP 11, 086 (2013).  DOI 10.1007/JHEP11(2013)086
  * (26) N. Behr, A. Mukhopadhyay, Phys. Rev. D 94(2), 026002 (2016).  DOI 10.1103/PhysRevD.94.026002
  * (27) N. Behr, S. Kuperstein, A. Mukhopadhyay, Phys. Rev. D 94(2), 026001 (2016).  DOI 10.1103/PhysRevD.94.026001
  * (28) A. Mukhopadhyay, Int. J. Mod. Phys. A 31(34), 1630059 (2016).  DOI 10.1142/S0217751X16300593
  * (29) S. de Haro, S.N. Solodukhin, K. Skenderis, Commun. Math. Phys. 217, 595 (2001).  DOI 10.1007/s002200100381
  * (30) K. Skenderis, Class. Quant. Grav. 19, 5849 (2002).  DOI 10.1088/0264-9381/19/22/306
  * (31) R.K. Gupta, A. Mukhopadhyay, JHEP 03, 067 (2009).  DOI 10.1088/1126-6708/2009/03/067
  * (32) P.M. Chesler, L.G. Yaffe, Phys. Rev. Lett. 106, 021601 (2011).  DOI 10.1103/PhysRevLett.106.021601
  * (33) A. Soloviev, Eur. Phys. J. C 82(4), 319 (2022).  DOI 10.1140/epjc/s10052-022-10282-4
  * (34) S.F. Taghavi, Phys. Rev. C 104(5), 054906 (2021).  DOI 10.1103/PhysRevC.104.054906
  * (35) S. Acharya, et al., JHEP 03, 092 (2024).  DOI 10.1007/JHEP03(2024)092
  * (36) U. Gursoy, M. Jarvinen, G. Policastro, JHEP 01, 134 (2016).  DOI 10.1007/JHEP01(2016)134
  * (37) C. Chattopadhyay, S. Jaiswal, L. Du, U. Heinz, S. Pal, Phys. Lett. B 824, 136820 (2022).  DOI 10.1016/j.physletb.2021.136820
  * (38) S. Jaiswal, C. Chattopadhyay, L. Du, U. Heinz, S. Pal, Phys. Rev. C 105(2), 024911 (2022).  DOI 10.1103/PhysRevC.105.024911
  * (39) Z. Chen, L. Yan, Phys. Rev. C 105(2), 024910 (2022).  DOI 10.1103/PhysRevC.105.024910
  * (40) E. Iancu, A. Mukhopadhyay, JHEP 06, 003 (2015).  DOI 10.1007/JHEP06(2015)003
  * (41) S. Banerjee, N. Gaddam, A. Mukhopadhyay, Phys. Rev. D 95(6), 066017 (2017).  DOI 10.1103/PhysRevD.95.066017
  * (42) A. Kurkela, A. Mukhopadhyay, F. Preis, A. Rebhan, A. Soloviev, JHEP 08, 054 (2018).  DOI 10.1007/JHEP08(2018)054
  * (43) T. Mitra, S. Mondkar, A. Mukhopadhyay, A. Rebhan, A. Soloviev, Phys. Rev. Res. 2(4), 043320 (2020).  DOI 10.1103/PhysRevResearch.2.043320
  * (44) T. Mitra, S. Mondkar, A. Mukhopadhyay, A. Rebhan, A. Soloviev, JHEP 04, 041 (2024).  DOI 10.1007/JHEP04(2024)041
  * (45) G.T. Landi, M. Paternostro, Rev. Mod. Phys. 93, 035008 (2021).  DOI 10.1103/RevModPhys.93.035008.  URL https://link.aps.org/doi/10.1103/RevModPhys.93.035008
  * (46) T. Kibe, A. Mukhopadhyay, P. Roy, Phys. Rev. Lett. 128(19), 191602 (2022).  DOI 10.1103/PhysRevLett.128.191602
  * (47) A. Banerjee, T. Kibe, N. Mittal, A. Mukhopadhyay, P. Roy, Phys. Rev. Lett. 129(19), 191601 (2022).  DOI 10.1103/PhysRevLett.129.191601
  * (48) V. Chandrasekaran, R. Longo, G. Penington, E. Witten, JHEP 02, 082 (2023).  DOI 10.1007/JHEP02(2023)082
  * (49) M. Banados, C. Teitelboim, J. Zanelli, Phys. Rev. Lett. 69, 1849 (1992).  DOI 10.1103/PhysRevLett.69.1849
  * (50) R. Emparan, JHEP 06, 036 (1999).  DOI 10.1088/1126-6708/1999/06/036
  * (51) J. Maldacena, G.L. Pimentel, JHEP 02, 038 (2013).  DOI 10.1007/JHEP02(2013)038