aainstitutetext: Departament de Física Quàntica i Astrofísica and Institut de
Ciències del Cosmos (ICC),
Universitat de Barcelona, Martí i Franquès 1, ES-08028, Barcelona,
Spainbbinstitutetext: Institut für Theoretische Physik, Goethe Universität,
Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germanyccinstitutetext:
Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís
Companys 23,
ES-08010, Barcelona, Spain.ddinstitutetext: Theoretical Physics Department,
CERN, CH-1211 Genève 23, Switzerland

# Strong-coupling dynamics and entanglement in de Sitter space

Jorge Casalderrey-Solana b Christian Ecker a,c David Mateos d and Wilke van
der Schee<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>frankfurt.de<EMAIL_ADDRESS>
###### Abstract

We use holography to study the dynamics of a strongly-coupled gauge theory in
four-dimensional de Sitter space with Hubble rate $H$. The gauge theory is
non-conformal with a characteristic mass scale $M$. We solve Einstein’s
equations numerically and determine the time evolution of homogeneous gauge
theory states. If their initial energy density is high compared with $H^{4}$
then the early-time evolution is well described by viscous hydrodynamics with
a non-zero bulk viscosity. At late times the dynamics is always far from
equilibrium. The asymptotic late-time state preserves the full de Sitter
symmetry group and its dual geometry is a domain-wall in AdS5. The approach to
this state is characterised by an emergent relation of the form
$\mathcal{P}=w\,\mathcal{E}$ that is different from the equilibrium equation
of state in flat space. The constant $w$ does not depend on the initial
conditions but only on $H/M$ and is negative if the ratio $H/M$ is close to
unity. The event and the apparent horizons of the late-time solution do not
coincide with one another, reflecting its non-equilibrium nature. In between
them lies an “entanglement horizon” that cannot be penetrated by extremal
surfaces anchored at the boundary, which we use to compute the entanglement
entropy of boundary regions. If the entangling region equals the observable
universe then the extremal surface coincides with a bulk cosmological horizon
that just touches the event horizon, while for larger regions the extremal
surface probes behind the event horizon.

††preprint: ICCUB-20-026, CERN-TH-2020-195

##  1 Introduction

Understanding the dynamics of non-Abelian gauge theories beyond the weak-
coupling limit is an important challenge. In this regime a quasi-particle
description is likely not applicable and one must resort to a different
intuition in order to understand the physics, especially out of equilibrium.
Holography provides a powerful framework with which a variety of theories can
be analysed from first principles in this regime. In this context the quasi-
particle intuition is replaced by intuition based on higher-dimensional
gravity, black hole horizons, etc. In the case in which the gauge theory
dynamics takes place in flat space, holography has provided valuable
qualitative insights into the properties of Quantum Chromodynamics (QCD),
especially into the far-from-equilibrium dynamics of its deconfined phase (see
e.g. CasalderreySolana:2011us  and references therein). These properties are
explored experimentally via the small drops of deconfined QCD matter that are
created in Heavy Ion Collisions (HIC) (for a recent review, see Busza:2018rrf
). Since these little fireballs are violently produced at an initial
temperature just a few times the QCD deconfinement temperature, the physics
immediately after the collision is non-weakly coupled and far from
equilibrium. One of the insights provided by holography is that the system
becomes well described by hydrodynamics at a time at which the viscous
corrections are still very large Heller:2011ju ; Chesler:2010bi ; Casalderrey-
Solana:2013aba . During its subsequent evolution the fireball expands, cools
down and eventually hadronises. In QCD in equilibrium, this transition is
realised as a smooth crossover Aoki:2006we .

Extending the analysis to the dynamics of gauge theories in curved spacetime
is interesting from several viewpoints. At a theoretical level, the spacetime
curvature may lead to new effects and hence richer dynamics.
Phenomenologically, one motivation comes from Cosmology, where the dynamics of
the gauge theory is coupled to an expanding spacetime. This situation was
certainly realised about one microsecond after the Big Bang, when the
decreasing temperature of the Universe crossed the QCD critical temperature
and quarks and gluons became bound into hadrons. Another interesting scenario
comes from the possibility that the physics beyond the Standard Model might be
completed at some high-energy scale by a Grand Unified non-Abelian gauge
Theory (GUT). In this scenario there may be implications for the early
Universe, and out-of-equilibrium effects may arise if, for example, the GUT
theory undergoes a phase transition (see for example Schwaller:2015tja ;
Caprini:2015zlo ). Finally, it has recently been suggested that dark matter
may be strongly self-interacting (recent reviews include Kribs:2016cew ;
Tulin:2017ara ), in which case a complete understanding of the dark sector
would require going beyond perturbative methods.

This paper is an exploratory investigation aimed at understanding the dynamics
of strongly coupled matter in a cosmological context via holography. We
emphasize that, unlike in e.g. Strominger:2001pn ; McFadden:2009fg , our goal
is not to provide a dual holographic description of the cosmological
gravitational field but only of the strongly coupled matter that lives in this
background, as in e.g. Koyama:2001rf ; Marolf:2010tg ; Ghoroku:2012vi ;
Buchel:2016cbj ; Buchel:2017lhu ; Buchel:2017pto ; Buchel:2019qcq ;
Buchel:2019pjb . We will therefore assume that the four-dimensional
gravitational field is prescribed a priori and use five-dimensional gravity to
describe only the dynamics of the four-dimensional gauge theory. We also
stress that, at this very early stage, we are mainly motivated by theoretical
curiosity, with the phenomenological motivation being mostly inspirational.
For this reason we will not be guided by an attempt to describe a realistic
scenario but rather make a number of simplifying assumptions.

The first one is that we will ignore the backreaction of the matter on the
expanding metric. The second one is that we will consider the simplest
possible expanding geometry, namely de Sitter (dS) space. The third
simplification is that we will restrict our attention to spatially homogeneous
states. And the fourth one concerns the gauge theory that we will study. This
will be defined by the condition that it be a four-dimensional, non-conformal
theory with the simplest possible gravity dual. The non-conformal nature of
the theory is absolutely crucial in order to uncover the physics that we are
interested in. The reason is that de Sitter space is conformal to Minkowski
space. Roughly speaking this means that, up to the effect of the conformal
anomaly, the physics of a conformal theory in dS is the same as in flat space
Apostolopoulos_2009 .

Despite these simplifications, we will still be able to capture several novel
effects. These include the fact that at late times the apparent and the event
horizons do not coincide, reflecting the non-equilibrium nature of the state,
or the existence of an entanglement horizon in between them that cannot be
penetrated by extremal surface anchored at the boundary.

The rest of the paper is structured as follows. In section 2 we introduce the
holographic model. In section 3 we review the thermodynamic and transport
properties of the model on flat space and discuss how we fix ambiguities due
to anomalies in the hydrodynamic approximation. In section 4 we introduce the
numerical algorithm we use to solve the dual gravity problem, explain how we
construct initial states and analyse our results for the time evolution of the
model on de Sitter space, the properties of late time states and how they are
approached. In section 5 we discuss entanglement and horizon entropies. We end
with a summary and discussion in section 6.

##  2 Holographic model

We follow a bottom-up approach and use five-dimensional Einstein-dilaton
gravity with non-trivial potential to model the dynamics of a strongly coupled
field theory with broken conformal symmetry in four dimensions. We will
consider the same holographic model as in Attems:2016ugt . That reference
explored the thermodynamics and the transport properties of the dual gauge
theory in flat space. We will review these properties in Sec. 3. In the
current section we will focus on the extension that is needed on the
holographic side in order to describe the dynamics in de Sitter space.

The action of the holographic model is given by

$S=\frac{2}{8\pi
G}\int_{\mathcal{M}}d^{5}x\sqrt{-g}\left(\frac{1}{4}R[g]-\frac{1}{2}(\partial\phi)^{2}-V(\phi)\right)+\frac{1}{8\pi
G}\int_{\partial\mathcal{M}}d^{4}x\sqrt{-\gamma}K+S_{ct}\,.$ (1)

Here $G$ is the five-dimensional Newton’s constant, $R[g]$ is the Ricci scalar
associated to the five-dimensional bulk metric $g_{\mu\nu}$ on $\mathcal{M}$,
$\gamma_{ij}$ is the metric induced on a four-dimensional slice near the
boundary $\partial\mathcal{M}$, and

$K=\gamma^{ij}K_{ij}=\gamma^{ij}\nabla_{i}n_{j}$ (2)

is the trace of the extrinsic curvature $K_{ij}$ associated to this slice. The
second term on the right-hand side of (1) is the familiar Gibbons–Hawking
term. The third term in (1) will be described below. The equations of motion
take the form

$\displaystyle R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}$
$\displaystyle=2\partial_{\mu}\phi\partial_{\nu}\phi-2V(\phi)g_{\mu\nu}-(\partial\phi)^{2}g_{\mu\nu}\,,$
(3a) $\displaystyle\nabla^{2}\phi$ $\displaystyle=\frac{\partial
V}{\partial\phi}\,.$ (3b)

The potential $V(\phi)$ encodes the properties of the dual gauge theory. We
wish to choose the simplest possible potential with the following two
properties: (i) it describes a _non_ -conformal theory, and (ii) the vacuum of
the theory in flat space is described by a completely regular solution on the
gravity side. Following Ref. Attems:2016ugt  we therefore choose the potential

$L^{2}V(\phi)=-3-\frac{3}{2}\phi^{2}-\frac{1}{3}\phi^{4}+\left(\frac{1}{3\phi_{M}^{2}}+\frac{1}{2\phi_{M}^{4}}\right)\phi^{6}-\frac{1}{12\phi_{M}^{4}}\phi^{8}\,,$
(4)

which can be derived from the superpotential

$L\,W\left(\phi\right)=-\frac{3}{2}-\frac{\phi^{2}}{2}+\frac{\phi^{4}}{4\phi^{2}_{M}}$
(5)

via the relation

$V(\phi)=-\frac{4}{3}W\left(\phi\right)^{2}+\frac{1}{2}W^{\prime}\left(\phi\right)^{2}\,.$
(6)

$L$ is a length scale. The dimensionless constant $\phi_{M}$ is a free
parameter that controls the degree of non-conformality of the model, for
example the maximum value of the bulk viscosity. For concreteness, in this
paper we will choose

$\phi_{M}=2\,.$ (7)

Both $V(\phi)$ and $W(\phi)$ have a maximum at $\phi=0$ and a minimum at
$\phi=\phi_{M}$. Each of these extrema yields an AdS solution of the equations
of motion with constant $\phi$ and radius $L^{2}=-3/V(\phi)$. In the gauge
theory each of these solutions is dual to a fixed point of the Renormalisation
Group (RG) with a number of degrees of freedom $N^{2}$ proportional to
$L^{3}/G$. In top-down models this relation is known precisely. For example,
in the case in which the gauge theory is $\mathcal{N}=4$ SYM with $N$ colours
we would have

$\frac{L^{3}}{8\pi G}=\frac{N^{2}}{4\pi^{2}}\,.$ (8)

In our bottom-up model we will take this as a definition of the number of
degrees of freedom in the gauge theory, $N$, at each fixed point.

The potential (4) leads to three important properties of the model: First, the
resulting geometry is asymptotically AdS5 in the UV with radius $L$, since
$V(0)=-3/L^{2}$. Second, the second derivative of the potential at $\phi=0$
implies that, in this asymptotic region, the scalar field has mass
$m^{2}=-3/L^{2}$. Following the standard quantisation analysis this means
that, in the UV, this field is dual to an operator in the gauge theory,
$\hat{{\mathcal{O}}}$, with scaling dimension $\Delta_{\textrm{UV}}=3$. The
value of the source $M$ of this operator introduces a scale responsible for
the breaking of conformal invariance. Third, the solution near $\phi=\phi_{M}$
is again AdS5 with a different radius

$L_{\rm
IR}=\sqrt{-\frac{3}{V\left(\phi_{M}\right)}}=\frac{1}{1+\frac{1}{6}\phi_{M}^{2}}L\,.$
(9)

In this region the effective mass of the scalar field differs from its UV
value and it is given by

${m^{2}_{\rm
IR}}=\frac{12}{L^{2}}\left(1+\frac{1}{9}\phi_{M}^{2}\right)=\frac{12}{L^{2}_{IR}}\frac{\left(1+\frac{1}{9}\phi_{M}^{2}\right)}{\left(1+\frac{1}{6}\phi_{M}^{2}\right)^{2}}\,.$
(10)

As a consequence, the operator $\hat{{\mathcal{O}}}$ at the IR fixed point has
dimension

$\Delta_{\textrm{IR}}=2+2\sqrt{1+\frac{{m^{2}_{\rm IR}}L_{\rm
IR}^{2}}{4}}=6\,\left(1+\frac{\phi_{M}^{2}}{9}\right)\left(1+\frac{\phi_{M}^{2}}{6}\right)^{-1}\,.$
(11)

As we will review in Sec. 3, when the gauge theory is placed in flat space
there exists an RG flow between the UV and the IR fixed points. The crossover
takes place at the scale $M$ and the geometry dual to the entire flow is
completely regular. In most expressions below we will fix the radius of the UV
AdS solution to unity, i.e., we will set $L=1$.

In order to understand the UV properties of the theory, such as anomalies, UV
divergences, etc. we will solve the Einstein’s equations near the boundary in
a power expansion in the so-called Fefferman–Graham (FG) coordinate $\rho$.
This has dimensions of (length)2 and in terms of it the near-boundary metric
takes the form

$ds^{2}=\frac{{\rm d}\rho^{2}}{4\rho^{2}}+\gamma_{ij}(\rho,x){\rm d}x^{i}{\rm
d}x^{j}\,.$ (12)

The boundary is located at $\rho=0$ and is parametrised by the coordinates
$x^{i}$ with $i=0,\ldots,3$. Near the boundary the metric and the scalar field
take the form

$\displaystyle\gamma_{ij}(\rho,x)$
$\displaystyle=\frac{1}{\rho}\Bigg{\\{}g_{(0)ij}(x)+\rho\,g_{(2)ij}(x)+\rho^{2}\Big{[}g_{(4)ij}(x)+h_{(4)ij}(x)\log\rho\Big{]}+O(\rho^{3})\Bigg{\\}}\,,$
(13a) $\displaystyle\phi(\rho,x)$
$\displaystyle=\rho^{1/2}\Bigg{\\{}\phi_{(0)}(x)+\rho\,\Big{[}\phi_{(2)}(x)+\psi_{(2)}(x)\log\rho\Big{]}+O(\rho^{2})\Bigg{\\}}\,.$
(13b)

As we will see below, the logarithmic terms are related to the presence of
anomalies. The first term $g_{(0)ij}(x)$ is the boundary metric. We will be
interested in the case in which this is a maximally symmetric, four-
dimensional spacetime with constant positive curvature $R=12H^{2}$, namely a
dS4 metric with Hubble rate $H$:

$ds_{\mathrm{b}}^{2}=g_{(0)ij}{\rm d}x^{i}{\rm d}x^{j}=-{\rm
d}t^{2}+e^{2Ht}{\rm d}\vec{x}^{2}\,.$ (14)

Similarly, we will assume that the first term in the expansion of the scalar
field is a constant that defines the characteristic mass scale in the gauge
theory:

$\phi_{(0)}=M\,.$ (15)

The first term of the action (1) suffers from large-volume divergences, as can
be verified by substituting the expansions (13) into the action. These
divergences can be regularised and renormalised by a procedure called
holographic renormalisation (see e.g. deHaro:2000vlm ; Bianchi:2001de ;
Bianchi:2001kw ), which makes the action finite and the variational principle
well-defined. This procedure is implemented by including in (1) the
counterterm action

$S_{\mathrm{ct}}=\frac{1}{8\pi G}\int{\rm
d}^{4}x\sqrt{-\gamma}\Bigg{[}\left(-\frac{1}{8}R[\gamma]-\frac{3}{2}-\frac{1}{2}\phi^{2}\right)+\frac{1}{2}\left(\log\rho\right)\mathcal{A}+\left(\alpha\mathcal{A}+\beta\phi^{4}\right)\Bigg{]}\,,$
(16)

defined on a timelike, constant-$\rho$ hypersurface near the boundary. The
induced metric on this hypersurface is denoted $\gamma_{ij}$ and $R[\gamma]$
is the associated Ricci scalar. The second term of (1) is also understood to
be evaluated on this slice at $\rho$, the first term of (1) is understood to
be evaluated by integrating down to this slice, and the limit $\rho\to 0$ is
understood to be taken at the end of the calculation.

In (16), $\mathcal{A}(\gamma_{ij},\phi)$ is the so-called conformal anomaly,
which in our case is given by

$\mathcal{A}=\mathcal{A}_{g}+\mathcal{A}_{\phi}$ (17)

where

$\mathcal{A}_{g}=\frac{1}{16}(R^{ij}R_{ij}-\frac{1}{3}R^{2})$ (18)

is the holographic gravitational conformal anomaly and

$\mathcal{A}_{\phi}=-\frac{\phi^{2}}{12}R$ (19)

is the conformal anomaly due to matter. In these equations all the terms are
functionals of the metric $\gamma_{ij}$ and of the scalar field $\phi$ induced
on the $\rho$-hypersurface. However, making use of the expansions (13) we see
that the product with the determinant of the induced metric yields a finite
contribution in the limit in which the cut-off is removed, since

$\lim_{\rho\to
0}\,\sqrt{-\gamma}\,\mathcal{A}\left(\gamma_{ij},\phi\right)=\lim_{\rho\to
0}\,\left[\frac{1}{\rho^{4}}\,\sqrt{-g_{(0)}}\right]\Big{[}\rho^{4}\mathcal{A}\left(g_{(0)ij},\phi_{(0)}\right)\Big{]}=\sqrt{-g_{(0)}}\,\mathcal{A}\left(g_{(0)ij},\phi_{(0)}\right)\,.$
(20)

For this reason we will often think of the anomaly as evaluated on the
boundary values of the fields, in which case (18) and (19) reduce to

$\displaystyle\mathcal{A}_{g}$ $\displaystyle=$
$\displaystyle-\frac{3}{4}H^{4}\,,$ (21) $\displaystyle\mathcal{A}_{\phi}$
$\displaystyle=$ $\displaystyle-M^{2}H^{2}\,.$ (22)

The fact that $\sqrt{-\gamma}\mathcal{A}$ yields a finite result has two
consequences. First, it means that the logarithmic term in (16) cancels a
purely logarithmic divergence from the bulk action. The requirement that this
cancellation takes place fixes uniquely the form of the anomaly, including the
values of all the numerical coefficients in (18) and (19). The presence of
this logarithmic term on the gravity side breaks diffeomorphism invariance and
is dual to the presence of the conformal anomaly in the dual gauge theory.

The second consequence is that the anomaly itself (without the $\log$) can be
added to the counterterm action with an arbitrary coefficient, which we named
$\alpha$ in (16). It is important to note that not just the anomaly but any
local, finite term that is invariant under the symmetries of the theory can be
added to the counterterm action with an arbitrary coefficient. These terms can
be constructed out of non-negative powers of the scalar field and of curvature
invariants of the induced metric $\gamma_{ij}$ in such a way that their
overall mass dimension is four. The $\beta\phi^{4}$ term is an example of such
a term. Other possible terms include combinations of $R_{ij}R^{ij}$, $R^{2}$
and $\phi^{2}$ that are linearly independent of $\mathcal{A}$. Therefore we
could replace the last term of (16) by

$\alpha\mathcal{A}+\beta\phi^{4}+\delta_{1}\left(R_{ij}R^{ij}+R^{2}\right)+\delta_{2}\,\phi^{2}R+\cdots\,.$
(23)

The freedom to add these terms with arbitrary coefficients $\delta_{i}$ is
part of the general freedom in the choice of renormalisation scheme. The
coefficient $\alpha$ plays a special role since it can be shifted by a scale
transformation, which is implemented via the following rescaling of the
coordinates

$x_{i}=\lambda x_{i}^{\prime}\,,\qquad\rho=\lambda^{2}\rho^{\prime}\,,$ (24)

where $\lambda$ is a positive real number. It is easy to see that the effect
of this transformation is to shift the counterterm action by a term of the
form $(\log\lambda)\mathcal{A}$, which in turn can be absorbed through the
redefinition $\alpha\to\alpha+\log\lambda$. The freedom to rescale $\rho$, or
equivalently to shift $\alpha$, is the freedom to choose a renormalisation
scale. We thus see that the freedom to choose a renormalisation scheme
includes, but is larger than, the freedom to choose a renormalisation scale.
This statement is well known on the gauge theory side. In order to renormalise
the theory it is not enough to choose a renormalisation scale since finite
parts must also be fixed. For example, the difference between the MS and
$\overline{\mbox{MS}}$ schemes is precisely the choice of the finite parts.
Below we will discuss the effect of the anomaly on the gauge theory
observables of interest to us, namely the expectation values of the stress
tensor and of the scalar operator.

The value $\beta=1/4\phi_{M}^{2}$ is special because in this case the
$\beta\phi^{4}$ term combines with the second and the third summands in the
first term of (16) to give precisely the superpotential (5). This means that,
if the theory (1) is the bosonic truncation of a supersymmetric theory with
superpotential $W$, then in flat space this choice of $\beta$ leads to a
supersymmetric renormalisation scheme. Motivated by this discussion, in this
paper we will set to zero all the coefficients in (23) but $\alpha$ and
$\beta$. This implies no loss of generality since physically meaningful
quantities are scheme-independent.

Substitution of the expansions (13) in the equations of motion (3) determines
several coefficients Bianchi:2001kw . The Klein–Gordon equation for the scalar
field fixes the logarithmic coefficient $\psi_{(2)}$ in terms of $g_{(0)ij}$
and $\phi_{(0)}$ as

$\psi_{(2)}=\frac{1}{24}\phi_{(0)}R=\frac{1}{2}MH^{2}\,.$ (25)

Unless otherwise indicated, in this and in subsequent equations it is
understood that the curvature tensors are those associated to the boundary
metric $g_{(0)ij}$. At leading order Einstein’s equations determine

$g_{(2)ij}=-\frac{1}{2}\left(R_{ij}-\frac{1}{6}R\,g_{(0)ij}\right)-\frac{\phi_{(0)}^{2}}{3}g_{(0)ij}=-\left(\frac{1}{2}H^{2}+\frac{1}{3}M^{2}\right)g_{(0)ij}\,.\\\
$ (26)

The logarithmic part at subleading order fixes

$h_{(4)ij}=h_{(4)ij}^{\mbox{\tiny{grav}}}-\frac{1}{12}R_{ij}\phi_{(0)}^{2}\,,$
(27)

where

$\displaystyle h_{(4)ij}^{\mbox{\tiny{grav}}}=$
$\displaystyle\frac{1}{8}R_{ikjl}R^{kl}-\frac{1}{48}\nabla_{i}\nabla_{j}R+\frac{1}{16}\nabla^{2}R_{ij}-\frac{1}{24}RR_{ij}$
$\displaystyle+\left(\frac{1}{96}R^{2}-\frac{1}{96}\nabla^{2}R-\frac{1}{32}R_{kl}R^{kl}\right)g_{(0)ij}$
(28)

is the purely gravitational part. For conformally flat metrics, such as (14),
this part vanishes, hence

$h_{(4)ij}=-\frac{1}{12}R_{ij}\phi_{(0)}^{2}=-\frac{1}{4}H^{2}M^{2}g_{(0)ij}\,.$
(29)

The subleading non-logarithmic part of Einstein’s equations fixes the trace of
$g_{(4)ij}$:

$\displaystyle\mathrm{Tr}g_{(4)ij}$
$\displaystyle=-2\phi_{(0)}{\phi}_{(2)}+\frac{5}{72}R\phi_{(0)}^{2}+\frac{1}{16}(R_{ij}R^{ij}-\frac{2}{9}R^{2})+\frac{2}{9}\phi_{(0)}^{4}$
(30a)
$\displaystyle=-2M{\phi}_{(2)}+\frac{5}{6}M^{2}H^{2}+\frac{1}{4}H^{4}+\frac{2}{9}M^{4}\,,$
(30b)

as well as its covariant divergence

$\displaystyle\nabla^{j}g_{(4)ij}\,=\,\,$
$\displaystyle\nabla^{j}\Bigg{\\{}-\frac{1}{8}\Big{[}\mathrm{Tr}\left(g_{(2)}^{\,\,\,\,\,2}\right)-\left(\mathrm{Tr}g_{(2)}\right)^{2}\Big{]}g_{(0)ij}+\frac{1}{2}\left(g_{(2)}^{\,\,\,\,\,2}\right)_{ij}-\frac{1}{4}g_{(2)ij}\,\mathrm{Tr}\,g_{(2)}$
$\displaystyle\qquad\,-\frac{3}{2}h_{(4)ij}-g_{(0)ij}\phi_{(0)}\Big{(}\phi_{(2)}+\psi_{(2)}\Big{)}\Bigg{\\}}$
(31a) $\displaystyle\,=\,\,$
$\displaystyle\nabla^{j}\Bigg{\\{}\left(\frac{1}{4}H^{4}+\frac{5}{24}H^{2}M^{2}+\frac{1}{9}M^{4}-M{\phi}_{(2)}\right)g_{(0)ij}\Bigg{\\}}\,.$
(31b)

The coefficients in (13) determine the holographic stress tensor as

$\displaystyle\langle\hat{T}_{ij}\rangle$ $\displaystyle=\lim\limits_{\rho\to
0}\frac{2\rho^{-2}}{\sqrt{-\gamma}}\frac{\delta S}{\delta\gamma^{ij}}$
$\displaystyle=2\left(\frac{N^{2}}{4\pi^{2}}\right)\Bigg{\\{}g_{(4)ij}+\frac{1}{8}\left[\mathrm{Tr}g_{(2)}^{2}-(\mathrm{Tr}g_{(2)})^{2}\right]g_{(0)ij}-\frac{1}{2}g_{(2)}^{2}+\frac{1}{4}g_{(2)ij}\mathrm{Tr}g_{(2)}$
$\displaystyle\qquad\,\,\,+\phi_{(0)}\left(\phi_{(2)}-\frac{1}{2}\psi_{(2)}\right)g_{(0)ij}+\alpha\left(\mathcal{T}^{g}_{ij}+\mathcal{T}_{ij}^{\phi}\right)+\left(\frac{1}{18}+\beta\right)\phi_{(0)}^{4}g_{(0)ij}\Bigg{\\}}\,.$
(32)

In this and in subsequent equations we have made use of (8) to replace $G$ in
favour of $N$, which makes the expected $N^{2}$-scaling of the stress tensor
manifest. The contributions $\mathcal{T}^{g}_{ij}$ and
$\mathcal{T}_{ij}^{\phi}$ come from the variation of $\mathcal{A}_{g}$ and
$\mathcal{A}_{\phi}$ in (16), respectively, and are given by:

$\displaystyle\frac{1}{2}\mathcal{T}^{g}_{ij}\,=\,\,$ $\displaystyle
h_{(4)ij}^{\mbox{\tiny{grav}}}=0\,,$ (33)
$\displaystyle\mathcal{T}^{\phi}_{ij}\,=\,\,$
$\displaystyle-\frac{1}{6}\phi_{(0)}^{2}\left(R_{ij}-\frac{1}{2}Rg_{(0)ij}\right)=\frac{1}{2}M^{2}H^{2}g_{(0)ij}\,.$
(34)

As we mentioned above, the first equation follows from the conformal flatness
of the dS metric (14). The expectation value of the scalar operator in the
field theory is given by

$\langle\hat{{\mathcal{O}}}\rangle=\lim\limits_{\rho\to
0}\frac{\rho^{-\Delta_{\mathrm{UV}}/2}}{\sqrt{-\gamma}}\frac{\delta
S}{\delta\phi}=2\left(\frac{N^{2}}{4\pi^{2}}\right)\left\\{-2\phi_{(2)}+(1-4\alpha)\psi_{(2)}-4\beta\phi_{(0)}^{3}\right\\}\,.$
(35)

In the presence of external sources the holographic stress tensor satisfies
anomaly-corrected Ward identities. These can be obtained from the variation of
the renormalised on-shell action

$\delta S[\delta g_{(0)},\delta\phi_{(0)}]=\int{\rm
d}^{4}x\sqrt{g_{(0)}}\left(\frac{1}{2}\langle\hat{T}_{ij}\rangle\delta
g_{(0)}^{ij}+\langle\hat{{\mathcal{O}}}\rangle\delta\phi_{(0)}\right)\,.$ (36)

Invariance of the action under diffeomorphisms

$\delta
g_{(0)}^{ij}=-(\nabla^{i}\xi^{j}+\nabla^{j}\xi^{i})\,,\quad\delta\phi_{(0)}=\xi^{i}\nabla_{i}\phi_{(0)}\,,$
(37)

leads to the diffeomorphism Ward identity

$\nabla^{i}\langle\hat{T}_{ij}\rangle=-\langle\hat{{\mathcal{O}}}\rangle\nabla_{j}\phi_{(0)}\,.$
(38)

Weyl transformations

$\delta g_{(0)}^{ij}=-2\sigma
g^{ij}\,,\quad\delta\phi_{(0)}=-(d-\Delta_{\mathrm{UV}})\sigma\phi_{(0)}\,,$
(39)

give the anomaly-corrected conformal Ward identity

$\langle\hat{T}^{i}_{i}\rangle=-(d-\Delta_{\mathrm{UV}})\phi_{(0)}\langle\hat{{\mathcal{O}}}\rangle+\left(\frac{N^{2}}{4\pi^{2}}\right)\left(\mathcal{A}_{g}+\mathcal{A}_{\phi}\right)\,,$
(40)

where $d$ is the spacetime dimension of the boundary theory and
$\mathcal{A}_{g}$ and $\mathcal{A}_{\phi}$ are given in (18) and (19),
respectively. In our case the Ward identities reduce to

$\nabla^{i}\langle\hat{T}_{ij}\rangle=0\,,\qquad\langle\hat{T}^{i}_{i}\rangle=-M\langle\hat{{\mathcal{O}}}\rangle-\left(\frac{N^{2}}{4\pi^{2}}\right)\left(\frac{3}{4}H^{4}+M^{2}H^{2}\right)\,.$
(41)

We are now ready to discuss the effect of the anomaly on physical observables
such as the stress tensor and the scalar operator. To see this consider again
the rescaling (24). In the gauge theory this is equivalent to rescaling $H$
and $M$ as

$H^{\prime}=\lambda H\,,\qquad M^{\prime}=\lambda M\,.$ (42)

Following Bianchi:2001de , we note that the rescaling above leaves the FG form
of the metric (12) invariant and transforms all the expansion coefficients
homogeneously,

$g_{(0)ij}^{\prime}=g_{(0)ij}\,,\qquad
g_{(2)ij}^{\prime}=\lambda^{2}g_{(2)ij}\,,\qquad
h_{(4)ij}^{\prime}=\lambda^{4}h_{(4)ij}\,,$ (43)

except for $g_{(4)ij}$, which acquires an inhomogeneous piece due to the
logarithmic term in (13):

$g_{(4)ij}^{\prime}=\lambda^{4}g_{(4)ij}+2\lambda^{4}\log\lambda\,h_{(4)ij}\,.$
(44)

Similarly, the coefficients $\phi_{(0)}$ and $\psi_{(2)}$ in the expansion of
the scalar field transform homogeneously, whereas $\phi_{(2)}$ acquires an
inhomogeneous piece:

$\phi_{(2)}^{\prime}=\lambda^{3}\phi_{(2)}+2\lambda^{3}\log\lambda\,\psi_{(2)}\,.$
(45)

It follows that the stress tensor and the scalar expectation value transform
as

$\displaystyle\langle\hat{T}_{ij}^{\prime}\rangle$ $\displaystyle=$
$\displaystyle\lambda^{4}\langle\hat{T}_{ij}\rangle+4\left(\frac{N^{2}}{4\pi^{2}}\right)\lambda^{4}\log\lambda\,h_{(4)ij}\,,$
(46) $\displaystyle\langle\hat{{\mathcal{O}}}^{\prime}\rangle$
$\displaystyle=$
$\displaystyle\lambda^{3}\langle\hat{{\mathcal{O}}}\rangle-8\left(\frac{N^{2}}{4\pi^{2}}\right)\lambda^{3}\log\lambda\,\psi_{(2)}\,,$
(47)

namely

$\displaystyle\langle\hat{T}_{ij}(\lambda H,\lambda M)\rangle$
$\displaystyle=$
$\displaystyle\lambda^{4}\langle\hat{T}_{ij}(H,M)\rangle-\lambda^{4}\log\lambda\,\left(\frac{N^{2}}{4\pi^{2}}\right)\,H^{2}M^{2}\,g_{(0)ij}\,,$
(48) $\displaystyle\langle\hat{{\mathcal{O}}}(\lambda H,\lambda M)\rangle$
$\displaystyle=$
$\displaystyle\lambda^{3}\langle\hat{{\mathcal{O}}}(H,M)\rangle-\lambda^{3}\log\lambda\,\left(\frac{N^{2}}{4\pi^{2}}\right)\,4MH^{2}\,,$
(49)

where we have made use of (25) and (29). This immediately implies that these
expectation values must take the form

$\displaystyle\langle\hat{T}_{ij}(H,M)\rangle$ $\displaystyle=$ $\displaystyle
H^{4}\,t_{ij}\left(\frac{H}{M}\right)-\log\left(\frac{H}{\mu}\right)\,\left(\frac{N^{2}}{4\pi^{2}}\right)\,H^{2}M^{2}\,g_{(0)ij}\,,$
(50) $\displaystyle\langle\hat{{\mathcal{O}}}(H,M)\rangle$ $\displaystyle=$
$\displaystyle
H^{3}\,f\left(\frac{H}{M}\right)-\log\left(\frac{H}{\mu}\right)\,\left(\frac{N^{2}}{4\pi^{2}}\right)\,4MH^{2}\,,$
(51)

where $\mu$ is some arbitrary reference scale, a remnant of the
renormalisation process much like the renormalisation scale in QFT. The first
and second terms on the right-hand sides transform homogeneously and
inhomogeneously under the rescaling (42), respectively. Needless to say, one
could rewrite the first terms in a variety of forms, for example as
$M^{4}t_{ij}(H/M)$ for the stress tensor, etc. Also, one could replace
$\log(H/\mu)$ by $\log(H/M)+\log(M/\mu)$, thus redefining

$t_{ij}\to
t_{ij}+\log\left(\frac{H}{M}\right)\,\left(\frac{N^{2}}{4\pi^{2}}\right)\,\frac{M^{2}}{H^{2}}\,g_{(0)ij}\,.$
(52)

Note also that there is no loss of generality in assuming that the scale $\mu$
is the same in both equations, since the difference can again be absorbed in a
redefinition of the homogeneous terms.

The key conclusion is that, because of the anomaly, expectation values in the
field theory do not only depend on the ratio $H/M$, but on the two independent
dimensionless ratios that can be built from $M,H$ and $\mu$. Put differently,
in order to specify the theory it is not enough to specify the ratio between
$H$ and $M$, but instead both scales must be specified independently with
respect to some arbitrary reference scale $\mu$. The freedom in the choice of
this scale is part of a bigger freedom in the choice of renormalisation
scheme, as we discussed around (23). Throughout this paper we will measure all
dimensionful quantities in units of $M$ and, when necessary, we will fix the
renormalisation scheme by specifying $\alpha$ and $\beta$.

In the flat-space limit, namely if $H=0$, both the anomaly (21) and its
contributions (33) and (34) to the stress tensor vanish identically. This
means that in this case both the stress tensor and the scalar operator
transform covariantly under scale transformations. In other words, the non-
homogeneous terms in the equations above vanish. Moreover, the only non-zero
finite term among all the possible ones in (23) is the $\beta\phi^{4}$ term.
This produces a contribution to the stress tensor (2) that shifts its value by
a term proportional to the boundary metric $g_{(0)ij}=\eta_{ij}$, namely it
shifts the energy density and the pressure by opposite amounts. Therefore the
choice of scheme in the flat-space case reduces entirely to fixing the energy
or the pressure of some reference state, for example that of the vacuum. We
will come back to this point in the next section.

In this work we will only consider states that are homogeneous and isotropic,
for which the associated energy momentum tensor takes the diagonal form

$\langle\hat{T}^{i}_{j}\rangle=\mathrm{diag}\left\\{-\epsilon(t),p(t),p(t),p(t)\right\\}\,.$
(53)

When plotting numerical results we will often use “reduced” quantities such as
reduced energy density, reduced pressure and reduced expectation value of the
scalar operator defined as

$\mathcal{E}(t)\equiv\frac{2\pi^{2}}{N^{2}}\,\epsilon(t)\,,\qquad\mathcal{P}(t)\equiv\frac{2\pi^{2}}{N^{2}}\,p(t)\,,\qquad{\mathcal{O}}(t)\equiv\frac{2\pi^{2}}{N^{2}}\,\langle\hat{\mathcal{O}}(t)\rangle\,.$
(54)

In terms of these variables, the trace Ward identity (41) takes the form

$\mathcal{E}(t)-3\mathcal{P}(t)=M{\mathcal{O}}(t)+\frac{1}{2}\left(\frac{3}{4}H^{2}+M^{2}H^{2}\right)\,.$
(55)

##  3 Dynamics in flat space

###  3.1 Thermodynamics and transport

In this section we review the most salient thermodynamic and transport
properties of the holographic model on flat space, studied in detail in
Attems:2016ugt . This is useful because later we will use thermal equilibrium
states on flat space to initialise the time evolution of non-equilibrium
states on dS4 and compare their evolution to viscous hydrodynamics with
transport coefficients presented in this section.

The gauge/gravity correspondence maps thermodynamic equilibrium states on the
gauge theory side to equilibrium black brane geometries on the gravity side.
In our case these are homogeneous and isotropic solutions of the equations of
motion (3) with a regular horizon and asymptotically AdS boundary conditions
for the metric and appropriate asymptotic scaling for the scalar field. A
convenient gauge to construct these solutions is one where the holographic
coordinate is identified with the scalar field111The function $H(\phi)$
appearing in this section should not be confused with the Hubble rate $H$
appearing throughout the whole paper.

$ds^{2}=e^{2A(\phi)}\left(-H(\phi){\rm d}\tau^{2}+{\rm
d}\vec{x}^{2}\right)-2e^{A(\phi)+B(\phi)}\,{\rm d}\tau{\rm d}\phi\,.$ (56)

In this gauge the boundary is located at $\phi=0$ and the value of the scalar
field at the horizon $\phi_{\mathrm{h}}$ is determined by the condition
$H(\phi_{\mathrm{h}})=0$. After introducing a master field

$G(\phi)=\frac{{\rm d}}{{\rm d}\phi}A(\phi)\,,$ (57)

the equations of motion (3) can be rewritten in terms of a single master
equation222As in Attems:2016ugt , we normalise the scalar field differently
than in Gubser:2008ny , which is the reason why some of the coefficients in
(58) differ from those in the corresponding master equation in Gubser:2008ny .
Gubser:2008ny ; Attems:2016ugt

$\frac{G^{\prime}(\phi)}{G(\phi)+\frac{4V(\phi)}{3V^{\prime}(\phi)}}=\frac{{\rm
d}}{{\rm
d}\phi}\log\left(\frac{1}{3G(\phi)}-2G(\phi)+\frac{G^{\prime}(\phi)}{2G(\phi)}-\frac{G^{\prime}(\phi)}{2\left(G(\phi)+\frac{4V(\phi)}{3V^{\prime}(\phi)}\right)}\right)\,.$
(58)

Close to the horizon a solution to the master equation can be expressed as a
power series:

$G(\phi)=-\frac{4V(\phi_{\mathrm{h}})}{3V^{\prime}(\phi_{\mathrm{h}})}+\frac{2}{3}(\phi-\phi_{\mathrm{h}})\left(\frac{V(\phi_{\mathrm{h}})V^{\prime\prime}(\phi_{\mathrm{h}})}{V^{\prime}(\phi_{\mathrm{h}})^{2}}-1\right)+O\left(\left(\phi-\phi_{\mathrm{h}}\right)^{2}\right)\,.$
(59)

Close to the boundary, $\phi\to 0$, the master field can be expanded as

$G(\phi)=-\frac{1}{\phi}+\cdots\,.$ (60)

In practice we obtain a one–parameter family of solutions for $G(\phi)$,
parametrised by the value of $\phi_{\mathrm{h}}$, by numerically integrating
(58) from a value of $\phi$ close to the horizon to a value close to the
boundary using boundary conditions for $G$ and $G^{\prime}$ constructed from
(59). The metric functions in (56) can then be obtained through the relations

$\displaystyle A(\phi)$ $\displaystyle=$
$\displaystyle-\log\left(\frac{\phi}{M}\right)+\int_{0}^{\phi}{\rm
d}\tilde{\phi}\left(G(\tilde{\phi})+\frac{1}{\tilde{\phi}}\right)\,,$ (61)
$\displaystyle B(\phi)$ $\displaystyle=$
$\displaystyle\log\left(\left|G(\phi)\right|\right)+\int_{0}^{\phi}{\rm
d}\tilde{\phi}\frac{2}{3G(\tilde{\phi})}\,,$ (62) $\displaystyle H(\phi)$
$\displaystyle=$
$\displaystyle-\frac{e^{2B(\phi)}\left(4V(\phi)+3G(\phi)V^{\prime}(\phi)\right)}{3G^{\prime}(\phi)}\,.$
(63)

The temperature and the entropy density of field theory states dual to these
numerically constructed geometries can be expressed in terms of (61) and (62)
evaluated at the horizon (see Gubser:2008ny ; Attems:2016ugt  for details):

$T=\frac{A(\phi_{\mathrm{h}})-B(\phi_{\mathrm{h}})}{4\pi}\,,\quad\quad
s=2\pi\left(\frac{N^{2}}{4\pi^{2}}\right)e^{3A(\phi_{\mathrm{h}})}\,.$ (64)

In Fig. 1 (left) we show the reduced entropy density $\mathcal{S}\equiv
2\pi^{2}s/N^{2}$ divided by $T^{3}$ as a function of $T/M$ for $\phi_{M}=2$.
The dotted and the dashed black lines indicate, respectively, the infinite-
and the zero-temperature limits (recall that we are setting $L=1$)

$\lim\limits_{T\to\infty}\,\frac{\mathcal{S}}{\pi^{4}T^{3}}=1\,,\qquad\lim\limits_{T\to
0}\,\frac{\mathcal{S}}{\pi^{4}T^{3}}=L_{\mathrm{IR}}^{3}=\frac{27}{125}\,,$
(65)

with $L_{\mathrm{IR}}$ given by (9). As explored in detail in Attems:2016ugt ,
for real values of $\phi_{M}$ the model has a smooth crossover between the IR
and UV fixed point.

Figure 1: (Left) Entropy density as function of temperature for $\phi_{M}=2$.
The black dashed and dotted lines indicate respectively the IR- and the UV-
limits (65). (Center) Ratio of reduced pressure to energy density as a
function of the reduced energy density. (Right) Ratio of bulk viscosity to
energy density as a function of the reduced energy density.

The energy density and the pressure can be extracted from the stress tensor
discussed in Sec. 2. In thermal equilibrium they can be equivalently obtained
from the knowledge of the entropy density as a function of the temperature via
the thermodynamic relations

$\displaystyle p_{\mathrm{eq}}$ $\displaystyle\equiv p_{0}+\int_{0}^{T}{\rm
d}T^{\prime}\,s(T^{\prime})\,,$ (66a) $\displaystyle h$
$\displaystyle\equiv\epsilon+p_{\mathrm{eq}}=Ts\,,$ (66b)

where $p_{0}$ is the pressure of the vacuum state in the $T\rightarrow 0$
limit and $h$ is the enthalpy density. The value of $p_{0}$ is not fixed by
thermodynamic considerations or by the equations of motion. In contrast, $T$
and $s$ are uniquely defined via (64), and hence so is the enthalpy density.
It follows that the energy density is also defined up to an arbitrary constant
equal to $-p_{0}$. As we saw towards the end of Sec. 2, the choice of
renormalisation scheme in flat space boils down precisely to the choice of
this constant. Thus one way to fix the scheme in flat space and uniquely
determine the stress tensor is to impose that the vacuum energy density and
pressure vanish. As discussed in Sec. 2, this corresponds to the choice
$\beta=1/4\phi_{M}^{2}$ and this is what we implicitly assume in the rest of
this section.

In Fig. 1 (center) we show the ratio of reduced pressure and energy density.
The equilibrium values for $p$ and $\epsilon$ constitute the equation of state
(EoS) $p_{\mathrm{eq}}(\epsilon)$. We use the deviation of the ratio $w\equiv
p/\epsilon={\mathcal{P}}/{\mathcal{E}}$ from 1/3 as a measure of the amount of
conformal symmetry breaking. Similarly, we could measure it via the deviation
of the speed of sound squared $c_{\mathrm{s}}^{2}\equiv{\rm
d}p_{\mathrm{eq}}/{\rm d}\epsilon$ from its conformal value
$c_{\mathrm{s,CFT}}^{2}\equiv 1/(d-1)=1/3$. One advantage of $w$ over
$c_{s}^{2}$ is that the former can be computed without taking derivatives, and
hence it is defined instantaneously, which will be useful when we study the
model in de Sitter space. In flat space, $w$ asymptotes to the conformal value
$w=1/3$ in the low- and high-energy density regimes. In between, at energy
densities comparable to the scale of the theory, ${\mathcal{E}}\approx M^{4}$,
$w$ deviates significantly from its conformal value.

At leading order in the hydrodynamic expansion (see below) the transport
properties in flat space are determined by two coefficients: the shear
viscosity and the bulk viscosity. The ratio of shear viscosity over entropy
density, $\eta/s=1/4\pi$, is universal in all holographic theories with an
Einstein gravity dual Kovtun:2004de . This means that knowledge of the entropy
density (66) is sufficient to determine the shear viscosity in our case. The
bulk viscosity can be obtained from the logarithmic derivative of the entropy
density with respect to the value of the scalar field at the horizon
Eling:2011ms :

$\frac{\zeta}{s}=\frac{1}{\pi}\left(\frac{{\rm d}\log s}{{\rm
d}\phi_{\mathrm{h}}}\right)^{-2}\,.$ (67)

In Fig. (1) (right) we plot the ratio of bulk viscosity and energy density as
a function of the reduced energy density. As we will see in Sec. 4.4, the
specific combination $9M\zeta/\epsilon$ measures the viscous contribution to
$w$ in dS4. The bulk viscosity vanishes at small and large energy densities
where the model is conformal. In between, at energy densities comparable to
the scale of the theory, ${\mathcal{E}}\approx M^{4}$, $\zeta$ is non-zero.

###  3.2 Hydrodynamics

The thermodynamic analysis of the previous section only applies to static
equilibrium states. It serves as a starting point for describing the long
wavelength dynamics of the system in a hydrodynamic approximation. The modern
interpretation of hydrodynamics as an effective theory in terms of a gradient
expansion has been reviewed many times (see e.g. romatschke_romatschke_2019
and references therein). The purpose of this section is to review some of the
basic definitions in order to fix notation and to show how scheme dependence
manifests itself in the hydrodynamic expansion.

In the absence of other conserved charges, long-wavelength excitations in the
gauge theory are solely controlled by the dynamics of the energy-momentum
tensor. In this limit the stress tensor can be approximated in terms of a
derivative expansion

$\displaystyle\langle\hat{T}^{ij}\rangle$ $\displaystyle=$
$\displaystyle\epsilon\,u^{i}u^{j}+p_{\mathrm{eq}}(\epsilon)\Delta^{ij}-\eta(\epsilon)\,\sigma^{ij}-\zeta(\epsilon)\Delta^{ij}\overline{\nabla}^{k}u_{k}+O(\overline{\nabla}^{2})\,,$
(68) $\displaystyle\sigma^{ij}$ $\displaystyle=$
$\displaystyle\Delta^{ik}\Delta^{jl}(\overline{\nabla}_{k}u_{l}+\overline{\nabla}_{l}u_{k})-\frac{2}{3}\Delta^{ij}\overline{\nabla}_{k}u^{k}\,,$
(69)

where $u^{i}$ is the fluid velocity and
$\overline{\nabla}^{i}\equiv\Delta^{ij}\nabla_{j}$, with
$\Delta^{ij}=g^{ij}_{(0)}+u^{i}u^{j}$, is the projection of the covariant
derivative to the spatial components in the local rest frame of the fluid. The
EoS $p_{\mathrm{eq}}(\epsilon)$ and the transport coefficients
$\eta(\epsilon)$ and $\zeta(\epsilon)$ are functions of the energy density
that depend on the microscopic details of the theory. The hydrodynamic
approximation involves a choice of hydrodynamic variables and specifying these
variables is called a choice of frame. We choose the Landau frame, in which
the velocity $u^{i}$ and energy density $\epsilon$ are defined as the timelike
eigenvector and the eigenvalue, respectively, of the stress tensor, namely
$\langle\hat{T}^{ij}\rangle u_{j}=-\epsilon u^{i}$.

The leading term in the expansion (68) is called the ideal hydrodynamic part.
It describes the flow of energy and momentum in terms of a locally
equilibrated ensemble locally boosted to non-vanishing velocity $u^{i}$.
Higher-order terms are expressed in terms of gradients of energy density and
fluid velocity. In the following we will neglect terms of
$O(\overline{\nabla}^{2})$ and only consider the leading ideal and sub-leading
viscous part of (68). We will write this term as

$\langle\hat{T}^{ij}\rangle=\Big{(}\epsilon+p_{\mathrm{eq}}(\epsilon)\Big{)}u^{i}u^{j}+p_{\mathrm{eq}}(\epsilon)\,g_{(0)}^{ij}+\Pi^{ij}\,,$
(70)

where the viscous tensor $\Pi^{ij}$ in general contains contributions due to
bulk and shear viscosity.333The appearance of shear and bulk viscosity as
leading contributions in $\Pi^{ij}$ is specific to our choice of Landau frame.
In other frames these contributions are different Kovtun:2019hdm .  In this
work we will only consider homogeneous and isotropic flows without shear
stresses, in which case the viscous tensor simplifies to

$\Pi^{ij}=-\left(g_{(0)}^{ij}+u^{i}u^{j}\right)\,\zeta(\epsilon)\overline{\nabla}_{k}u^{k}\,.$
(71)

The possible scheme dependence of the microscopic energy-momentum tensor also
manifests itself in the hydrodynamic approximation. As discussed in Sec. 3.1,
$p_{\mathrm{eq}}$ given in (66) is not uniquely defined but contains an
arbitrary contribution $p_{0}$ identified as the vacuum pressure. As explained
above, in flat space the entire freedom in the choice of renormalisation
scheme reduces to the choice of this constant. Thus one way to proceed is to
make an explicit choice and perform all calculations in that scheme.
Alternatively, we may work with manifestly scheme-independent quantities as
follows. We first define the excess pressure and the excess energy density
over the vacuum as

$\Delta\epsilon\equiv\epsilon+p_{0}\,,\qquad\Delta
p_{\mathrm{eq}}=p_{\mathrm{eq}}-p_{0}\,.$ (72)

These are scheme-independent, and we may then view the EoS as a relation of
the form $\Delta p_{\mathrm{eq}}=\Delta p_{\mathrm{eq}}(\Delta\epsilon)$. Next
we rewrite the second term in (70) as

$p_{\mathrm{eq}}g_{(0)}^{ij}=\Delta
p_{\mathrm{eq}}g_{(0)}^{ij}+\langle\hat{T}^{ij}_{\mathrm{vac}}\rangle\,,\quad\langle\hat{T}^{ij}_{\mathrm{vac}}\rangle\equiv
p_{0}g_{(0)}^{ij}\,.$ (73)

This separates (70) into a scheme-dependent vacuum contribution
$\langle\hat{T}^{ij}_{\mathrm{vac}}\rangle$ and a scheme-independent
contribution

$\langle\Delta\hat{T}^{ij}\rangle\equiv\langle\hat{T}^{ij}\rangle-\langle\hat{T}^{ij}_{\mathrm{vac}}\rangle\,.$
(74)

Note that the fact that $\langle\hat{T}^{ij}_{\mathrm{vac}}\rangle$ is
proportional to the background metric $g_{(0)}^{ij}$ is consistent with the
expectation that the vacuum must respect the symmetries of this background.
Since $\langle\Delta\hat{T}^{ij}\rangle$ and $\Delta\epsilon$ are scheme-
independent, we can now define a scheme-independent velocity field defined
through the relation

$\langle\Delta\hat{T}^{ij}\rangle u_{j}=-\Delta\epsilon\,u^{i}\,.$ (75)

The scheme-independent part of (70) is then given by

$\langle\Delta\hat{T}^{ij}\rangle=\Big{(}\Delta\epsilon+\Delta
p_{\mathrm{eq}}(\Delta\epsilon)\Big{)}u^{i}u^{j}+\Delta
p_{\mathrm{eq}}(\Delta\epsilon)g_{(0)}^{ij}-\left(g_{(0)}^{ij}+u^{i}u^{j}\right)\,\zeta(\Delta\epsilon)\overline{\nabla}_{k}u^{k}\,.$
(76)

##  4 Dynamics in de Sitter space

###  4.1 The dual gravity problem

We now turn to the main subject of this work: the far-from-equilibrium
dynamics of our strongly coupled non-conformal gauge theory on a time-
dependent background geometry. For this we numerically solve the fully non-
linear equations of motion (3) of the dual gravity theory, following the
method reviewed in e.g. Chesler:2013lia ; vanderSchee:2014qwa ; Ecker:2018jgh
, and extract the time evolution of the expectation values of the energy
momentum tensor and the scalar operator from the solution near the boundary.

We are ultimately interested in the evolution of field theory observables on
dS4. However, it is useful to set up the problem with a slightly more general
boundary metric of Friedmann–Lemaître–Robertson–Walker type

$ds_{\mathrm{b}}^{2}=g_{(0)ij}{\rm d}x^{i}{\rm d}x^{j}=-{\rm
d}t^{2}+S_{0}(t)^{2}{\rm d}\vec{x}^{2}\,.$ (77)

For $S_{0}(t)=e^{Ht}$ the boundary metric (77) becomes the dS4 (14) with
curvature scalar $R=12H^{2}$.

We use generalised Eddington–Finkelstein (EF) coordinates to parametrise the
bulk geometry and the scalar field

$ds^{2}=-A(r,t){\rm d}t^{2}+2{\rm d}r{\rm d}t+S(r,t)^{2}{\rm
d}\vec{x}^{2}\,,\quad\phi=\phi(r,t)\,,$ (78)

where the asymptotic boundary is located at $r=\infty$. The line element (78)
has a residual gauge freedom in the radial coordinate

$\displaystyle r$ $\displaystyle\rightarrow\bar{r}\equiv r+\xi(t)\,,$ (79a)
$\displaystyle A(r,t)$ $\displaystyle\rightarrow\bar{A}(\bar{r},t)\equiv
A(\bar{r}-\xi(t),t)+2\partial_{t}\xi(t)\,,$ (79b) $\displaystyle S(r,t)$
$\displaystyle\rightarrow\bar{S}(\bar{r},t)\equiv S(\bar{r}-\xi(t),t)\,,$
(79c)

which we exploit in our numerical scheme to fix the coordinate value
$r_{\mathrm{AH}}$ at the apparent horizon, defined by the condition
$\dot{S}(r_{\mathrm{AH}},t)=0$, to a constant.

Using (78) the equations of motion (3) result in the following set of
equations

$\displaystyle S^{\prime\prime}$
$\displaystyle=-\frac{2}{3}S\left(\phi^{\prime}\right)^{2}\,,$ (80a)
$\displaystyle\dot{S}^{\prime}$
$\displaystyle=-\frac{2\dot{S}S^{\prime}}{S}-\frac{2SV}{3}\,,$ (80b)
$\displaystyle\dot{\phi}^{\prime}$
$\displaystyle=\frac{V^{\prime}}{2}-\frac{3\dot{S}\phi^{\prime}}{2S}-\frac{3S^{\prime}\dot{\phi}}{2S}\,,$
(80c) $\displaystyle A^{\prime\prime}$
$\displaystyle=\frac{12\dot{S}S^{\prime}}{S^{2}}+\frac{4V}{3}-4\dot{\phi}\phi^{\prime}\,,$
(80d) $\displaystyle\ddot{S}$
$\displaystyle=\frac{\dot{S}A^{\prime}}{2}-\frac{2S\dot{\phi}^{2}}{3}\,,$
(80e)

where a prime denotes a radial derivative, $f^{\prime}\equiv\partial_{r}f$,
and an overdot is short-hand for the modified derivative
$\dot{f}\equiv\partial_{t}f+\frac{1}{2}A\partial_{r}f$.

Imposing the metric (77) and the asymptotic behaviour of the scalar field

$\lim\limits_{r\to\infty}\phi(r,t)=\frac{M}{r^{d-\Delta_{\mathrm{UV}}}}=\frac{M}{r}$
(81)

as boundary conditions, solutions of (80) can be expressed near the boundary
as generalised power series:

$\displaystyle A(r,t)=$
$\displaystyle\,\,r^{2}+2r\xi(t)+\xi(t)^{2}-2\xi^{\prime}(t)+\frac{S_{0}^{\prime}(t)^{2}-2S_{0}(t)S_{0}^{\prime\prime}(t)}{S_{0}(t)^{2}}-\frac{2M^{2}}{3}+\frac{a_{(4)}(t)}{r^{2}}$
$\displaystyle+\frac{2M^{2}S_{0}^{\prime\prime}(t)}{3S_{0}(t)}\frac{\log(r)}{r^{2}}+O(r^{-3})\,,$
(82a) $\displaystyle S(r,t)=$
$\displaystyle\,\,S_{0}(t)r+S_{0}^{\prime}(t)+\xi(t)S_{0}(t)-\frac{M^{2}S_{0}(t)}{3r}+\frac{M^{2}(3\xi(t)S_{0}(t)-S_{0}^{\prime}(t))}{9r^{2}}$
$\displaystyle+\frac{M(4M^{3}S_{0}(t)^{2}-72\bar{\phi}_{(2)}(t)S_{0}(t)^{2}+48M\xi(t)S_{0}(t)S_{0}^{\prime}(t)+9M(S_{0}^{\prime}(t)^{2}+S_{0}(t)S_{0}^{\prime\prime}(t)))}{216S_{0}(t)r^{3}}$
$\displaystyle+\frac{M^{2}(S_{0}^{\prime}(t)^{2}+S_{0}(t)S_{0}^{\prime\prime}(t))}{6S_{0}(t)}\frac{\log(r)}{r^{3}}+O(r^{-4})\,,$
(82b) $\displaystyle\phi(r,t)=$
$\displaystyle\,\,\frac{M}{r}-\frac{M\xi(t)}{r^{2}}+\frac{\bar{\phi}_{(2)}(t)}{r^{3}}-\frac{M\left(S_{0}^{\prime}(t)^{2}+S_{0}(t)S_{0}^{\prime\prime}(t)\right)}{2S_{0}(t)^{2}}\frac{\log(r)}{r^{3}}+O(r^{-4})\,.$
(82c)

Note that the fall-off coefficient $\bar{\phi}_{(2)}$ in EF coordinates is
generically different from the one in FG coordinates ${\phi}_{(2)}$. The
coefficients $a_{(4)}(t)$ and $\bar{\phi}_{(2)}(t)$ in these series cannot be
determined from the near-boundary analysis but need to be extracted from the
full bulk solution. The equations of motion impose the following relation on
these coefficients

$\displaystyle
a_{(4)}^{\prime}(t)+4a_{(4)}(t)\frac{S_{0}^{\prime}(t)}{S_{0}(t)}$
$\displaystyle=$
$\displaystyle-\frac{4}{3}M\left(2\bar{\phi}_{(2)}(t)\frac{S_{0}^{\prime}(t)}{S_{0}(t)}+\bar{\phi}_{(2)}^{\prime}(t)\right)+\frac{16M^{4}S_{0}^{\prime}(t)}{27S_{0}(t)}$
$\displaystyle+\frac{2}{3}M^{2}\left(\frac{S_{0}^{\prime}(t)S_{0}^{\prime\prime}(t)}{S_{0}(t){}^{2}}+\frac{4S_{0}{}^{(3)}(t)}{3S_{0}(t)}+\frac{4\xi(t)^{2}S_{0}^{\prime}(t)}{S_{0}(t)}+4\xi(t)\xi^{\prime}(t)\right)\,.$

This relation follows from the momentum constraint (80e) and implies covariant
conservation of the holographic stress tensor in the boundary theory, namely
the first Ward identity in (41).

The EF coordinate system (78) is useful to obtain time-dependent solutions of
the equations of motion numerically. However, our expressions for the
expectation values of the stress tensor (2) and the scalar operator (35)
assume the FG coordinate system (12). Although we could recompute the
corresponding expressions in EF gauge, it will prove more convenient to relate
the EF and the FG coefficients. For this purpose we need to find the
asymptotic coordinate transformation between the EF and the FG coordinate
systems. We first write a series ansatz for the EF coordinates $r_{EF}$ and
$t_{EF}$ in powers of the radial FG coordinate $r_{FG}$, related to $\rho$ in
(12) through $\rho=1/r_{FG}^{2}$:

$\displaystyle r_{EF}(r_{FG},t_{FG})$
$\displaystyle=\sum_{n=1}^{\infty}\Big{[}r_{n}(t_{FG})+\rho_{n}(t_{FG})\log(r_{FG})+\cdots\Big{]}(r_{FG})^{n}\,,$
(84a) $\displaystyle t_{EF}(r_{FG},t_{FG})$
$\displaystyle=t_{FG}+\sum_{n=1}^{\infty}\Big{[}t_{n}(t_{FG})+\tau_{n}(t_{FG})\log(r_{FG})+\cdots\Big{]}(r_{FG})^{-n}\,,$
(84b)

where dots stand for therms with higher powers of $\log(r_{FG})$. All these
logarithmic terms appear because we are working with a curved boundary metric.
The metric transforms as follows

$g^{FG}_{\mu\nu}=\frac{\partial x_{EF}^{\alpha}}{\partial
x_{FG}^{\mu}}\frac{\partial x_{EF}^{\beta}}{\partial
x_{FG}^{\nu}}g_{\alpha\beta}^{EF}\,,$ (85)

where $x_{EF}^{\mu}=(r_{EF},t_{EF})$ and $x_{FG}^{\mu}=(r_{FG},t_{FG})$.
Substituting the general expressions for metrics in EF and FG coordinates

$g^{EF}_{\mu\nu}=\begin{pmatrix}0&1\\\ 1&g^{EF}_{1,1}\end{pmatrix}\,,\qquad
g^{FG}_{\mu\nu}=\begin{pmatrix}r_{FG}^{2}&0\\\ 0&g^{FG}_{1,1}\end{pmatrix}$
(86)

into the transformation law (85) leads to a set of two equations

$\displaystyle 0$
$\displaystyle=\frac{\partial{t_{EF}}}{\partial{r_{FG}}}\frac{\partial{t_{EF}}}{\partial{t_{FG}}}g^{EF}_{1,1}+\frac{\partial{r_{EF}}}{\partial{r_{FG}}}\frac{\partial{t_{EF}}}{\partial{t_{FG}}}+\frac{\partial{t_{EF}}}{\partial{r_{FG}}}\frac{\partial{r_{EF}}}{\partial{t_{FG}}}\,,$
(87a) $\displaystyle 0$
$\displaystyle=r_{FG}^{2}-2\frac{\partial{r_{EF}}}{\partial{r_{FG}}}\frac{\partial{t_{EF}}}{\partial{r_{FG}}}-\left(\frac{\partial{t_{EF}}}{\partial{r_{FG}}}\right)^{2}g^{EF}_{1,1}\,,$
(87b)

which can be solved order by order in $r_{FG}$ and $r_{FG}\log(r_{FG})$ for
the expansion coefficients in (84). For the leading behaviour of the metric in
FG coordinates we find the following expressions

$\displaystyle g^{FG}_{1,1}$
$\displaystyle=-r^{2}+\frac{\bar{\phi}_{(0)}^{2}}{3}-\frac{S_{0}^{\prime}(t)^{2}-2S_{0}(t)S_{0}^{\prime\prime}(t)}{2S_{0}(t)^{2}}-\frac{\bar{\phi}_{(0)}^{2}S_{0}^{\prime\prime}(t)}{2S_{0}(t)}\frac{\log(r)}{r^{2}}-\left(\frac{\bar{\phi}_{(0)}^{4}}{36}+\frac{3a_{(4)}(t)}{4}\right)\frac{1}{r^{2}}$
$\displaystyle+\frac{2\bar{\phi}_{(0)}^{2}S_{0}(t)^{2}\big{(}2S_{0}^{\prime}(t)^{2}-3S_{0}(t)S_{0}^{\prime\prime}(t)\big{)}-3(S_{0}^{\prime}(t)^{2}-2S_{0}(t)S_{0}^{\prime\prime}(t))^{2})}{48S_{0}(t)^{4}}\frac{1}{r^{2}}+O(r^{-4})\,,$
(88a) $\displaystyle g^{FG}_{2,2}$
$\displaystyle=r^{2}S_{0}(t)^{2}-\frac{1}{3}\bar{\phi}_{(0)}^{2}S_{0}(t)^{2}-\frac{1}{2}S_{0}^{\prime}(t)^{2}+\frac{\bar{\phi}_{(0)}^{2}}{6}\left(2S_{0}^{\prime}(t)^{2}+S_{0}(t)S_{0}^{\prime\prime}(t)\right)\frac{\log(r)}{r^{2}}$
$\displaystyle+S_{0}(t)^{2}\left(\frac{19}{108}\bar{\phi}_{(0)}^{4}-\frac{1}{4}a_{(4)}(t)-\frac{2}{3}\bar{\phi}_{(0)}\bar{\phi}_{(2)}(t)\right)\frac{1}{r^{2}}$
$\displaystyle+\left(\frac{27S_{0}^{\prime}(t)^{4}+186\bar{\phi}_{(0)}^{2}S_{0}(t)^{3}S_{0}^{\prime\prime}(t)+288\bar{\phi}_{(0)}^{2}S_{0}(t)^{4}\xi(t)^{2}}{432S_{0}(t)^{2}}\right)\frac{1}{r^{2}}+O(r^{-4})\,,$
(88b)

where we have set $r\equiv r_{FG}$ and $t\equiv t_{FG}$ to shorten notation.
After replacing the radial coordinate by $\rho\equiv 1/r^{2}$ we obtain
explicit expressions for the non-vanishing metric components in FG coordinates
(12) in terms of the coefficients in EF gauge

$g_{(0)tt}=-1,\quad g_{(0)xx}=S_{0}(t)^{2}\,,$ (89)
$g_{(2)tt}=\frac{\bar{\phi}_{(0)}^{2}}{3}-\frac{S_{0}^{\prime}(t){}^{2}-2S_{0}(t)S_{0}^{\prime\prime}(t)}{2S_{0}(t){}^{2}},\quad
g_{(2)xx}=-\frac{2\bar{\phi}_{(0)}^{2}S_{0}(t)^{2}+3S_{0}^{\prime}(t)^{2}}{6}\,,$
(90) $\displaystyle g_{(4)tt}$
$\displaystyle=-\frac{3a_{(4)}}{4}-\frac{\bar{\phi}_{(0)}^{4}}{36}+\frac{\bar{\phi}_{(0)}^{2}\left(2S_{0}^{\prime}(t){}^{2}-3S_{0}(t)S_{0}^{\prime\prime}(t)\right)}{24S_{0}(t){}^{2}}-\frac{\left(S_{0}^{\prime}(t){}^{2}-2S_{0}(t)S_{0}^{\prime\prime}(t)\right){}^{2}}{16S_{0}(t){}^{4}}\,,$
(91a) $\displaystyle g_{(4)xx}$
$\displaystyle=-\frac{a_{(4)}(t)S_{0}(v){}^{2}}{4}-\frac{2\bar{\phi}_{(0)}\bar{\phi}_{(2)}(t)S_{0}(t){}^{2}}{3}+\frac{2}{3}\bar{\phi}_{(0)}^{2}\xi(v)^{2}S_{0}(v)^{2}$
$\displaystyle+\frac{31\phi_{(0)}^{2}S_{0}(t)S_{0}^{\prime\prime}(t)}{72}+\frac{19\bar{\phi}_{(0)}^{4}S_{0}(v){}^{2}}{108}+\frac{S_{0}^{\prime}(v){}^{4}}{16S_{0}(v){}^{2}}\,,$
(91b)
$h_{(4)tt}=\frac{\bar{\phi}_{(0)}^{2}S_{0}^{\prime\prime}(t)}{4S_{0}(t)},\quad
h_{(4)xx}=-\frac{\bar{\phi}_{(0)}^{2}\left(2S_{0}^{\prime}(v){}^{2}+S_{0}(t)S_{0}^{\prime\prime}(t)\right)}{12}\,.$
(92)

The corresponding near-boundary expansion of the scalar field reads

$\displaystyle\phi(r,t)$
$\displaystyle=\bar{\phi}_{(0)}\rho^{1/2}+\left(\bar{\phi}_{(2)}(t)-\frac{\bar{\phi}_{(0)}^{3}}{6}+\frac{\bar{\phi}_{(0)}\left(S_{0}^{\prime}(t){}^{2}-2S_{0}(t)S_{0}^{\prime\prime}(t)\right)}{4S_{0}(t){}^{2}}-\bar{\phi}_{(0)}\xi(t)^{2}\right)\rho^{3/2}$
$\displaystyle+\frac{\bar{\phi}_{(0)}\left(S_{0}^{\prime}(t)^{2}+S_{0}(t)S_{0}^{\prime\prime}(t)\right)}{4S_{0}(t)^{2}}\rho^{3/2}\log(\rho)+O(\rho^{5/2})\,.$
(93a)

From the above expression one can now read off the relations between the
expansion coefficients in FG and EF coordinates

$\displaystyle\phi_{(0)}$ $\displaystyle=$
$\displaystyle\bar{\phi}_{(0)}=M\,,$ (94) $\displaystyle\psi_{(2)}(t)$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\bar{\psi}_{(2)}(t)\,,$ (95)
$\displaystyle\phi_{(2)}(t)$ $\displaystyle=$
$\displaystyle\bar{\phi}_{(2)}(t)-\frac{\bar{\phi}_{(0)}^{3}}{6}+\frac{\bar{\phi}_{(0)}\left(S_{0}^{\prime}(t){}^{2}-2S_{0}(t)S_{0}^{\prime\prime}(t)-4\xi(t)^{2}S_{0}(t)^{2}\right)}{4S_{0}(t){}^{2}}\,.$
(96)

Using these relations we can express the non-vanishing components of the
holographic stress tensor in terms of the coefficients in the near boundary
expansion in EF gauge:

$\displaystyle\mathcal{E}(t)=$
$\displaystyle\,\,-\frac{3a_{(4)}(t)}{4}-M\bar{\phi}_{(2)}(t)+\frac{3S_{0}^{\prime}(t)^{4}}{16S_{0}(t)^{4}}+M^{2}\left(\xi(t)^{2}+\frac{S_{0}^{\prime}(t)^{2}}{8S_{0}(t)^{2}}+\frac{2S_{0}^{\prime\prime}(t)}{3S_{0}(t)}\right)$
$\displaystyle-M^{2}\alpha\frac{S_{0}^{\prime}(t)^{2}}{2S_{0}(t)^{2}}-M^{4}\left(\beta-\frac{7}{36}\right)\,,$
(97a) $\displaystyle\mathcal{P}(t)=$
$\displaystyle\,\,-\frac{a_{(4)}(t)}{4}+\frac{1}{3}M\bar{\phi}_{(2)}(t)+\frac{S_{0}^{\prime}(t)^{2}\left(S_{0}^{\prime}(t)^{2}-4S_{0}(t)S_{0}^{\prime\prime}(t)\right)}{16S_{0}(t)^{4}}$
$\displaystyle-\frac{M^{2}}{3}\left(\xi(t)^{2}+\frac{S_{0}^{\prime}(t)^{2}}{8S_{0}(t)^{2}}+\frac{13S_{0}^{\prime\prime}(t)}{12S_{0}(t)}\right)+M^{2}\alpha\left(\frac{S_{0}^{\prime}(t)^{2}}{6S_{0}(t)^{2}}+\frac{S_{0}^{\prime\prime}(t)}{3S_{0}(t)}\right)+M^{4}\left(\beta-\frac{5}{108}\right)\,.$
(97b)

Similarly, the expectation value of the operator $\hat{{\mathcal{O}}}$ is
given by

$\displaystyle{\mathcal{O}}(t)=$
$\displaystyle\,\,-2\bar{\phi}_{(2)}(t)+M\left(2\xi(t)^{2}-\frac{S_{0}^{\prime}(t)^{2}}{4S_{0}(t)^{2}}+\frac{5S_{0}^{\prime\prime}(t)}{4S_{0}(t)}\right)$
$\displaystyle-M\alpha\left(\frac{S_{0}^{\prime}(t)^{2}}{S_{0}(t)^{2}}+\frac{S_{0}^{\prime\prime}(t)}{S_{0}(t)}\right)-M^{3}\left(4\beta-\frac{1}{3}\right)\,.$
(98)

As mentioned above, the coefficients $\alpha$ and $\beta$ in (97) and (4.1)
encode the scheme-dependence of $\mathcal{E},\mathcal{P}$ and $\mathcal{O}$.

###  4.2 Numerical procedure

Our main interest is to compute the time evolution of $\mathcal{P}$,
$\mathcal{E}$ and ${\mathcal{O}}$ in a dynamical background geometry. For this
purpose we have to solve the set of equations (80) with consistent initial and
boundary conditions for the metric and the scalar field. The preparation of
initial states and our specific choice of time-dependent boundary conditions
for the metric will be discussed in the next section. Here we concentrate on
the evolution algorithm assuming that the initial data and the boundary metric
are known.

The set of equations (80) has a nested structure that allows us to treat them
on every slice of constant EF time as ordinary differential equations in the
radial coordinate for the functions $S$, $\dot{S}$, $\phi$, $\dot{\phi}$ and
$A$. Knowing these functions on a single time slice is sufficient to evolve
them to the next slice. In practice we do not directly solve (80) but rather a
set of equations for a new set of functions
$\\{\tilde{S},\tilde{\dot{S}},\tilde{\phi},\tilde{\dot{\phi}},\tilde{A}\\}$.
These are obtained from the original variables by subtracting all the
divergent terms, as well as some finite ones, that are explicitly known from
the near-boundary analysis (93). The new functions are therefore manifestly
finite and take the form

$\displaystyle A(r,t)=$
$\displaystyle\,\,r^{2}+2r\xi(t)+\xi(t)^{2}-2\xi^{\prime}(t)+\frac{S_{0}^{\prime}(t)^{2}-2S_{0}(t)S_{0}^{\prime\prime}(t)}{S_{0}(t)^{2}}-\frac{2M^{2}}{3}$
$\displaystyle+\log(r)\left(\frac{A_{\log,1}(t)}{r^{2}}+\frac{A_{\log,2}(t)}{r^{3}}\right)+\frac{\tilde{A}(r,t)}{r^{2}}\,,$
(99a) $\displaystyle S(r,t)=$
$\displaystyle\,\,S_{0}(t)r+S_{0}^{\prime}(t)+\xi(t)S_{0}(t)-\frac{M^{2}S_{0}(t)}{3r}+\frac{M^{2}(3\xi(t)S_{0}(t)-S_{0}^{\prime}(t))}{9r^{2}}$
$\displaystyle+\log(r)\left(\frac{S_{\log,1}(t)}{r^{3}}+\frac{S_{\log,2}(t)}{r^{4}}+\frac{S_{\log,3}(t)}{r^{5}}\right)+\frac{\log(r)^{2}}{r^{5}}S_{\log,4}(t)+\frac{\tilde{S}(r,t)}{r^{3}}\,,$
(99b) $\displaystyle\phi(r,t)=$
$\displaystyle\,\,\frac{M}{r}-\frac{M\xi(t)}{r^{2}}+\log(r)\left(\frac{\phi_{\log,1}(t)}{r^{3}}+\frac{\phi_{\log,2}(t)}{r^{4}}\right)+\frac{\tilde{\phi}(r,t)}{r^{3}}\,,$
(99c) $\displaystyle\dot{S}(r,t)=$
$\displaystyle\,\,\frac{1}{2}r^{2}S_{0}(t)+r\left(S_{0}(t)\xi(t)+S_{0}^{\prime}(t)\right)+\frac{3\left(S_{0}(t)\xi(t)+S_{0}^{\prime}(t)\right)^{2}-M^{2}S_{0}(t)^{2}}{6S_{0}(t)}-\frac{2M^{2}S_{0}^{\prime}(t)}{9r}$
$\displaystyle+\log(r)\left(\frac{\dot{S}_{\log,1}(t)}{r^{2}}+\frac{\dot{S}_{\log,2}(t)}{r^{3}}\right)+\frac{\tilde{\dot{S}}(r,t)}{r^{2}}\,,$
(99d) $\displaystyle\dot{\phi}(r,t)=$
$\displaystyle\,\,-\frac{M}{2}+\log(r)\left(\frac{\dot{\phi}_{\log,1}(t)}{r^{2}}+\frac{\dot{\phi}_{\log,2}(t)}{r^{3}}\right)+\frac{\tilde{\dot{\phi}}(r,t)}{r^{2}}\,.$
(99e)

The functions $A_{\log,1}(t)$, $A_{\log,2}(t)$, etc. are explicitly known from
the near-boundary analysis but are too long to be displayed here. Furthermore,
we switch to the inverse radial coordinate $z\equiv 1/r$ where the boundary is
located a $z=0$.

We now implement the following procedure to solve the initial value problem:

  1. 1.

For a given radial profile of the scalar field $\tilde{\phi}(z,t_{0})$ at some
initial time $t_{0}$ we first solve the second-order Hamiltonian constraint
equation (80a) for $\tilde{S}(z,t_{0})$. In principle this differential
equation requires boundary conditions, but in our (subtracted) formulation
this equation (and the following except the one for $\tilde{\dot{S}}$) has
regular singular points, and demanding regularity of the solution fixes the
boundary condition. The simplest way to achieve this is to use spectral
methods (see below), which are by construction regular.

  2. 2.

Next we use $\tilde{\phi}(z,t_{0})$ and $\tilde{S}(z,t_{0})$ in (80b) and
solve for $\tilde{\dot{S}}(z,t_{0})$. The boundary condition for this function
reads

$\displaystyle\tilde{\dot{S}}(z=0,t_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{36}S_{0}(t_{0})\Big{[}18M\left(\bar{\phi}_{2}(t_{0})-M\xi(t_{0})^{2}\right)+18a_{4}(t_{0})-5M^{4}\Big{]}$
(100)
$\displaystyle+\frac{1}{144}M^{2}\Big{[}32\xi(t_{0})S_{0}^{\prime}(t_{0})-61S_{0}^{\prime\prime}(t_{0})\Big{]}+\frac{3M^{2}S_{0}^{\prime}(t_{0})^{2}}{16S_{0}(t_{0})}\,.$

Evaluating this requires knowledge of $a_{4}(t_{0})$ and
$\bar{\phi}_{2}(t_{0})$, whereby $a_{4}(t_{0})$ is required as a separate
initial condition and $\bar{\phi}_{2}(t_{0})$ needs to be extracted from the
initial data for the scalar field.

  3. 3.

Once $\tilde{\phi}(z,t_{0})$, $\tilde{S}(z,t_{0})$ and
$\tilde{\dot{S}}(z,t_{0})$ are known we solve (80c) for
$\tilde{\dot{\phi}}(z,t_{0})$.

  4. 4.

Once $\tilde{\phi}(z,t_{0})$, $\tilde{S}(z,t_{0})$, $\tilde{\dot{S}}(z,t_{0})$
and $\tilde{\dot{\phi}}(z,t_{0})$ are known we can solve (80d) for
$\tilde{A}(z,t_{0})$.

  5. 5.

Once $\tilde{A}(z,t_{0})$ and $\tilde{\dot{\phi}}(z,t_{0})$ are known, we use
the definition of the dot-derivative

$\tilde{\dot{\phi}}(z,t_{0})=\partial_{t}\tilde{\phi}(z,t_{0})+\frac{1}{2}z^{2}\tilde{A}(z,t_{0})\partial_{z}\tilde{\phi}(z,t_{0})$
(101)

to solve for time derivative of the initial data,
$\partial_{t}\tilde{\phi}(z,t_{0})$. This equation depends on the gauge choice
$\xi^{\prime}(t_{0})$, which we fix by demanding that the apparent horizon
stays at a constant value of the radial coordinate. This is done by solving
$\partial_{t}\dot{S}(z_{\rm AH},t)=0$ for $\xi^{\prime}(t_{0})$.

  6. 6.

Finally, we compute $a_{4}^{\prime}(t_{0})$ from (4.1), whereby we obtain
$\phi_{2}^{\prime}(t_{0})$ from the near-boundary expansion of
$\partial_{t}\tilde{\phi}(z,t_{0})$. We subsequently evolve the initial data
to the next time slice $t_{1}=t_{0}+\Delta t$ for example via

$a_{4}(t_{1})=a_{4}(t_{0})+a_{4}^{\prime}(t_{0})\Delta
t\,,\qquad\tilde{\phi}(z,t_{1})=\tilde{\phi}(z,t_{0})+\partial_{t}\tilde{\phi}(z,t_{0})\Delta
t\,,$ (102)

and the start over at the first entry in this list.

The only element missing in the above discussion is the initial value of the
gauge function $\xi(t)$. We use this freedom to fix the initial apparent
horizon (AH) at $z_{\rm AH}=1/2$. This is done by computing $\dot{S}$ for the
gauge $\xi_{0}(t_{0})=0$ first, and then numerically solving the apparent
horizon equation $\dot{S}(z_{\rm hor},\,t_{0})=0$ (see also (121)). We then
update

$\xi_{1}(t_{0})=r_{\rm AH}-r_{\rm hor}=\xi_{0}(t_{0})+2-1/z_{\rm hor}\,,$
(103)

where $z_{\rm hor}=1/r_{\rm hor}$ is the current location of the apparent
horizon. To increase precision this procedure is repeated a few times, each
time using an updated $\xi(t_{0})$.

We solve the resulting equations numerically with a Chebyshev pseudo-spectral
method (see e.g. boyd01 ) using typically $N=60$ grid points in the radial
direction. For the stepping between time slices (sixth entry in the list
above) we use a fourth-order Adams–Bashforth method press2007numerical  with a
time step $\Delta t=1/(10N^{2})$. After each time step we evaluate the
momentum constraint (80e) to monitor the accuracy of the numerical evolution.

###  4.3 Initial states and time evolution

We are interested in studying the time evolution of gauge theory states in
dS4. There are of course many ways to construct initial states for this
evolution. The strategy that we will follow is to start with thermal
equilibrium states in flat space. We will then smoothly increase the value of
$H$ at the boundary from $H=0$ to its desired value in each case. This leads
to a transient period of time in which ${\rm d}H/{\rm d}t\neq 0$ and the
boundary geometry interpolates between flat space and dS4. After this time $H$
becomes constant and we are in the desired situation of studying the dynamics
of an excited gauge theory state in de Sitter space. A natural question that
arises is whether the resulting state at late times is sensitive to the
specific way in which we prepare the initial state. One of our main results is
that the answer to this question is negative. In this sense the way in which
we initialise the evolution implies no loss of generality.

The thermal equilibrium states in flat space could be constructed in FG
coordinates with the procedure presented in Sec. 3.1 and then numerically
transformed to the EF gauge (78), which is better suited for time evolution.
However, in practice it is simpler to construct the solutions directly in EF
coordinates by relaxation. To do this, we start with some initial guess for
$a_{4}(t_{0})$ and for the scalar field profile $\tilde{\phi}(z,t_{0})$. This
is in general an excited state, which we then evolve with flat boundary
conditions for the metric ($S_{0}=1$) using the algorithm outlined in the
previous section. After a few units of $tM$ the state relaxes to thermal
equilibrium.

Figure 2: Construction of initial states in flat space. All dimensionful
quantities are given in units of $M$ and we choose the renormalisation scheme
$\beta=1/16$. (Left) Initial guess (red) and resulting thermalised equilibrium
configuration (blue) for the scalar field profile. (Center) Evolution of
${\mathcal{E}}(t)$, ${\mathcal{P}}(t)$ and ${\mathcal{O}}(t)$. (Right) Ring-
down of $|{\mathcal{O}}(t)-{\mathcal{O}}(5)|$ (blue) and
$|{\mathcal{P}}(t)-{\mathcal{P}}(5)|$ (red) together with an exponential fit
for the lowest quasinomal mode $\omega_{1}\approx 6.00+4.64i$ (doted black).

An example of this procedure is shown in Fig. 2, where all dimensionful
quantities are given in units of $M$ and we have chosen $\beta=1/16$. In the
left plot we show the initial guess $\tilde{\phi}(z,t_{0})=0$ (red) together
with the thermalised equilibrium result $\tilde{\phi}(z,t=5)$ (blue). Dots
indicate the non-equidistant distribution of points on the Chebyshev grid in
the radial direction. For the subsequent evolution with time-dependent
boundary metric at least 60 grid points together with 70-digit-accurate
arithmetic is used. In the middle plot we show how ${\mathcal{E}}$,
${\mathcal{P}}$ and ${\mathcal{O}}$ evolve towards their equilibrium values.
Energy conservation in flat space and homogeneity of the state imply that the
energy density $\mathcal{E}$ remains constant during the evolution. The
evolutions of ${\mathcal{P}}$ and ${\mathcal{O}}$ are not independent but are
constrained by the Ward identity (55). Close to the equilibrium state both
oscillate according to the quasinormal ring-down shown in the plot on the
right. Our numeric evolution allows us to extract an estimate for this lowest
quasinormal mode at zero momentum given by $\omega_{1}\approx 6.00+4.64i$.

We now turn on the Hubble rate at the boundary so that the boundary metric
changes smoothly from Minkowski to dS4. We implement this by imposing the
following relation on the function $S_{0}(t)$ in the boundary metric (77)

$\frac{S_{0}^{\prime}(t)}{S_{0}(t)}=H\,\frac{\Big{[}1+\tanh\Big{(}\Omega\,(t-t_{\ast})\Big{)}\Big{]}}{2}.$
(104)

This relation mimics a time-dependent Hubble rate that changes from zero at
$t\ll t_{\ast}$ to $H$ at $t\gg t_{\ast}$ in a time of order $1/\Omega$. We
will refer to this transient period as “the quench”. The corresponding form of
$S_{0}(t)$ follows from integrating (104) subject to the initial condition
$S_{0}(t=0)=1$

$S_{0}(t)=e^{\frac{Ht}{2}}\left[\frac{\cosh\Big{(}\Omega(t-t_{\ast})\Big{)}}{\cosh(\Omega
t_{\ast})}\right]^{\frac{H}{2\Omega}}\,.$ (105)

In Fig. 3 we illustrate this procedure for several choices of parameters. We
will refer to each of these choices as a “protocol”.

Figure 3: (Left) Different protocols for $S_{0}^{\prime}(t)/S_{0}(t)$, all
using our maximal value of $H=3$. (Right) All protocols lead to the same late-
time value of $\mathcal{E}$.

We choose a theory with $H/M=3$, we measure all dimensionful quantities in the
figure in units of $M$ and we fix the renormalisation scheme by choosing
$\alpha=0,\beta=1/16$. On the left plot we show the ratio
$S_{0}^{\prime}(t)/S_{0}(t)$ and on the right the evolution of the energy
density. We vary both the “quench parameter” $\Omega$, which controls the
length of the transient period, and the parameter $a_{4}$, which controls the
initial energy density. We find that the state at $t\gg t_{\ast}$ does not
depend on the values of these parameters. In particular, we see in Fig. 3
(right) that the energy density at late times always approaches the same
value. This convergence to the same late-time state is remarkable in view of
the vastly different values of the initial energy density, and it means that
the late-time state only depends on $H$ and $M$. We will explore this
dependence in our subsequent simulations. Given the independence of the
initial conditions, in most simulations we will use the same values
$\Omega=4$, $t_{\ast}=1$ and $a_{4}(t=0)=-100$. The only exception will be
when $H\geq 2$, in which case we will use $\Omega=8$ in order to have a quench
time $\Omega^{-1}$ that is sufficiently shorter than the expansion rate
$H^{-1}$.

After a time $t-t_{\ast}\gg\Omega^{-1}$ the boundary geometry settles down to
dS4 and we can analyse the evolution of field theory states on this expanding
background. In Fig. 4 we show the evolution of $\mathcal{E}$, $\mathcal{P}$,
$\mathcal{O}$ and $\mathcal{E}+\mathcal{P}$ (which in equilibrium would be
enthalpy) for a number of different values for the expansion rate $H$.

Figure 4: (Top) Energy density (left) and pressure (right) as functions of
time for different values of $H$. (Bottom) Expectation value of the scalar
operator (left) and the combination $\mathcal{E}+\mathcal{P}$ as functions of
time for different values of $H$. The black dashed line is the expansion
factor $e^{-3Ht}$ of a fluid with zero pressure. All dimensionful quantities
are measured in units of $M$ and we set $\alpha=0,\beta=1/16$.

Shortly after the quench, the energy density and the pressure start to
decrease rapidly due to the expansion of space. At some later time the initial
energy of the plasma is almost entirely depleted. From this point onward the
energy and the pressure are dominated by their (scheme-dependent, see below)
vacuum contributions equivalent to a pure cosmological constant, meaning that
at late-times $\mathcal{E}=-\mathcal{P}$. The evolution of $\mathcal{O}$ is
entirely fixed in terms of the energy density and the pressure by the Ward
identity (40). The scheme-dependence of the late-time limits

$\mathcal{E}_{\infty}\equiv\lim\limits_{t\to\infty}\mathcal{E}(t)\,,\qquad\mathcal{P}_{\infty}\equiv\lim\limits_{t\to\infty}\mathcal{P}(t)\,,\qquad\mathcal{O}_{\infty}\equiv\lim\limits_{t\to\infty}\mathcal{O}(t)$
(106)

is encoded in the dependence on $\alpha$ and $\beta$ of Eqs. (97). As
mentioned above, in order to fix the scheme we choose
$\beta=1/4\phi_{M}^{2}=1/16$ throughout the paper. Moreover, in this section
(but not elsewhere) we set $\alpha=0$. Note that the enthalpy
$\mathcal{E}+\mathcal{P}$ is scheme-independent because the dependence on
$\alpha$ and $\beta$ cancels out in this combination. In Fig. 4 (bottom right)
we see that $\mathcal{E}+\mathcal{P}$ decays for any value of $H$, where to
guide the eye we have added a black dashed line proportional to $e^{-3Ht}$
(appropriate for a pressureless fluid).

As mentioned above, we set $\alpha=0$ only in this section. In the following
we will use variables which have their late time values subtracted. In other
words, we will measure energy, pressure, etc with respect to their late-time
asymptotic values. These subtracted quantities have the advantage that they
are scheme-independent because the dependence on $\alpha,\beta$ cancels out.
Conceptually, the reason for this is that the scheme dependence is the same
for any state. For any given values of $H,M$, working with subtracted
variables is equivalent to choosing a scheme in which

$\alpha=-\frac{2}{M^{2}H^{2}}\mathcal{E}_{\infty}|_{\alpha=0},$ (107)

with $\mathcal{E}_{\infty}|_{\alpha=0}$ the final numerical value of
$\mathcal{E}$ in Fig. 4. Note that we treat different values of $H,M$ as
corresponding to different theories for each of which a separate
renormalization scheme can be chosen. None of our later results will hence
directly compare stress tensors for different values of $H$ or $M$. The logic
behind the choice (107) is similar to the logic behind the choice
$\beta=1/4\phi_{M}^{2}$. In the latter case we require the energy and pressure
of the vacuum state to vanish. In the case of (107) we require the energy and
the pressure of the late-time, asymptotic state to vanish. In this scheme, the
values of the energy density and the pressure during the evolution can be
interpreted as those of excitations on top of the late-time state. In the next
subsection we will study the dynamics of these excitations.

###  4.4 Hydrodynamic regime

One of the lessons of holographic studies is that hydrodynamics becomes a good
description of a plasma once the quasi-normal modes of the system have
decayed. At strong coupling, the decay time is of order of the inverse
temperature $1/T$ in both conformal Heller:2011ju ; Chesler:2009cy ;
Chesler:2010bi  and non-conformal Buchel:2015saa ; Attems:2016ugt  theories.
Since the expansion rate is of order $1/H$, we expect that hydrodynamics will
provide a good description provided that $T>H$. In this regime the expansion
can then be seen as an almost-adiabatic process in which local properties of
matter are close to those in equilibrium. Under these circumstances the
expectation value of the energy-momentum tensor of the system can be
approximated in terms of the gradient expansion discussed in Sec. 3.2. If
$H\gg M$ or $H\ll M$ then the energy density at which the hydrodynamic
description ceases to be valid is close to the UV or to the IR fixed point of
the gauge theory. In these regions the dynamics is quasi-conformal, the bulk
viscosity is close to zero, and the relation between energy and pressure is
essentially determined by symmetry. It follows that the most interesting range
of parameters is $H\sim M$, which in our units means $H\sim 1$. Therefore we
will focus on the evolution of states with initial energy densities of order 1
and we will vary the value of $H$ from a few times larger to a few times
smaller than 1. As suggested by the discussion above, we will find a
qualitative change in the applicability of hydrodynamics around $H=1$.

In order to compare the evolution of the holographic energy-momentum tensor to
the hydrodynamic approximation it is convenient to work with hydrodynamic
variables that are manifestly scheme-independent, namely independent of
$\alpha$. To do this we start with the general form of the energy-momentum
tensor for homogeneous and isotropic states

$\langle\hat{T}^{ij}\rangle=\frac{N^{2}}{2\pi^{2}}\left\\{\left({\mathcal{E}}(t)+{\mathcal{P}}(t)\right)u^{i}u^{j}+{\mathcal{P}}(t)g^{ij}_{(0)}\right\\}\,,$
(108)

where $u^{i}=\left(1,\,\vec{0}\right)$ is a future pointing time-like vector
and ${\mathcal{E}}(t)$ and ${\mathcal{P}}(t)$ are reduced energy density and
pressure in the rest frame defined by $u^{i}$, respectively. This general form
is independent of the hydrodynamic approximation and applies in particular to
states preserving the symmetries of dS4, which in addition to (108) also
satisfy the relation ${\mathcal{E}}(t)=-{\mathcal{P}}(t)$. As shown in Fig. 4
(bottom right) and studied in more detail in Sec. 4.6, all our states are
attracted to states preserving the symmetries of dS4 such that the energy-
momentum tensor of these late-time states can be written as

$\langle\hat{T}^{ij}_{\infty}\rangle=\frac{N^{2}}{2\pi^{2}}\,{\mathcal{P}_{\infty}}g^{ij}_{(0)}\,,$
(109)

where ${\mathcal{P}_{\infty}}$ is scheme-dependent because it depends on
$\alpha$. In the spirit of Sec. 3.2 we build scheme-independent combinations
of the energy-momentum tensor components in which all dependencies on $\alpha$
cancel out:

$\langle\Delta\hat{T}^{ij}\rangle=\langle\hat{T}^{ij}\rangle-\langle\hat{T}^{ij}_{\infty}\rangle\,,\quad\Delta{\mathcal{E}}={\mathcal{E}}+{\mathcal{P}_{\infty}}\,,\quad\Delta{\mathcal{P}}={\mathcal{P}}-{\mathcal{P}_{\infty}}\,.$
(110)

These are the variables we use to build the hydrodynamic benchmark we compare
to the evolution obtained from solving the time-dependent dual gravity
problem. Using (110) in (76) we obtain the hydrodynamic approximation for the
pressure in the Landau frame444Both bulk viscosities in (111) and (68) are
manifestly scheme independent. A subtle question regarding scheme dependence
is however if these two functions are truly the same functions. Even though
they have both been obtained by the same prescription (defining energy
densities with respect to the maximally symmetric state), these backgrounds
are not precisely equal. Nevertheless, since the scheme dependence has the
form $\alpha H^{2}$ any difference will scale as $H^{2}$ and can hence be
attributed to higher order viscous effects that we do not take into account in
this paper.

$\Delta{\mathcal{P}^{{\rm
hydro}}}(t)\equiv\frac{N^{2}}{2\pi^{2}}\,\Big{\\{}\Delta
p_{\mathrm{eq}}(\Delta{\mathcal{E}}(t))-3H\zeta(\Delta{\mathcal{E}}(t))\Big{\\}}+O(H^{2})\,,$
(111)

where we have used $\overline{\nabla}_{i}u^{i}=3H$. This expression shows that
$H$ controls the size of spatial gradients in the velocity field, where $H=0$
gives the leading term in the gradient expansion which only depends on the
equilibrium pressure on flat space. In the following we study a number of
specific examples to see how well this approximation with and without the
viscous contribution agrees with the exact solution obtained from the full
holographic simulation.

Figure 5: Scheme-independent ratio of the excess pressure over the excess
energy density defined in (110) as a function of time for the same initial
state but different values of $H$.${}^{\ref{more}}$ The vertical grey bands
indicate the duration of the quench. The solid black curves labelled “strong
coupling” correspond to the exact result obtained holographically. Dotted blue
and dashed red curves correspond to the ideal and viscous hydrodynamic
approximations, respectively. In the plots with $H<1$ viscous hydrodynamics
provides a good approximation to the exact result for some time after the
quench. During part of this period the difference between ideal and viscous
hydrodynamics is of order unity, meaning that gradient corrections are as
large as the ideal terms. In the plots with $H\geq 1$ viscous hydrodynamics
never provides a good approximation to the exact result.

Fig. 5 shows the time evolution of the excess pressure divided by the excess
energy density, which in equilibrium is determined by the equation of state,
for six different values of $H$.555More precisely, what is shown is actually
$3\mathcal{P}/\mathcal{E}$ for the value of $\alpha$ specified in (107) which,
after the quench, coincides with $3\Delta\mathcal{P}/\Delta\mathcal{E}$.
Recall that all dimensionful quantites are measured in units of $M$. The black
lines show the exact, strongly coupled evolutions obtained holographically,
the red dashed lines show the viscous hydrodynamic approximation according to
(111) and the blue dotted lines show the ideal hydrodynamic approximation
according to the equation of state $\Delta{\mathcal{P}_{\rm
eq}}(\Delta{\mathcal{E}}(t))$. Note that the latter approaches unity at very
early and very late times, as expected from the existence of the UV and IR
fixed points in the gauge theory (see Fig. 1), but deviates from this value in
between due to the non-conformal nature of the gauge theory. Before the quench
all states are in thermal equilibrium on Minkowski and as expected all curves
agree. Although after the quench the gradients due to the dS expansion become
important, for $H<1$ there is a period of time during which the results are
well described by viscous hydrodynamics. For example, for $H=1/3,1/2$ this
agreement extends to times around $H(t-t_{*})\sim 1.4$. It is remarkable that
around these times the difference between ideal and viscous hydrodynamics is
of order one, indicating that the first-order viscous corrections are as large
as the ideal terms. This provides another example of hydrodynamics working
with large gradients Heller:2011ju ; Chesler:2010bi ; Casalderrey-
Solana:2013aba , in this case in a dynamical spacetime. For $H>1$ the
evolution after the quench is never well described by hydrodynamics and we
conclude that the evolution is far from equilibrium.

It is interesting to ask how the results of Fig. 5 would change when using a
higher initial energy density, which for $H=1,2$ and a 19.3 times higher
initial energy density is shown in Fig. 6. In that case the quench has a more
moderate effect on the state and the energy density will take longer to cool
down under the exponential expansion. Both these effects can at late times be
incorporated by an extra shift of $H\Delta t\approx 0.76$, after which the
$\Delta\mathcal{E}$ evolutions agree and hence also the viscous hydrodynamic
results. Before the quench the energy densities are different and
correspondingly the equation of state, as visible in the figure. It is also
visible that the lower energy density has a larger effect from the quench
(small wiggle around $t=0$), even though the effect is modest. This comparison
clearly shows that the statement that viscous hydrodynamics provides a good
approximation at early times for small expansion rates $H$ is valid
independent of the initial energy density, but it has to be interpreted at a
time where viscous corrections are sizeable, whereby this time can depend on
the initial energy density.

Figure 6: For $H=1$ (left) and $H=2$ (right) we compare the result of Fig. 5
with an equivalent run that starts with an energy density 19.3 times higher.
The original run is displayed in lighter colors and gray, whereas the high-
energy run is displayed in dark colors and black. The high-energy run is
shifted by an extra time of $H\Delta t=0.76$ for both figures, in order to
guarantee that the energy evolution of both simulations agree at late times.
After this shift in time there is virtually no difference between the results,
except small differences before the quenches (displayed gray and dark gray).

Importantly, as time passes the energy density in the expanding background
decreases. At some point the energy density will become of the same scale as
the gradients, and hydrodynamics need no longer apply. Indeed, we find that at
later time as $\Delta{\mathcal{E}}\rightarrow 0$ the full evolution always
deviates from the hydrodynamic evolution and can characteristically develop a
negative pressure excess for $H\gtrsim 1$. Note that since we compare the
pressure excess the sign of this pressure is unrelated to the negative
pressure of a cosmological constant contribution, rather it means that the
final pressure in Fig. 4 is approached from below. We stress that this
negative pressure is a completely out-of-equilibrium effect which may not be
inferred directly from the equilibrium dynamics of the plasma in flat space.
The next subsection will describe this late time evolution in more detail.

###  4.5 Effective equation of state and quasinormal modes

In the previous section we have shown that the late-time evolution of our
system is not well described by hydrodynamics. In this section we will
demonstrate that it is instead possible to describe the dynamics in terms of
small perturbations around the maximally symmetric late-time state. Because
the late-time state has a dual description in terms of a geometry with a
horizon, the dynamics of small perturbations in the gauge theory is determined
by “quasi-normal modes” (QNM) of the dual black brane. We have used quotation
marks in the previous sentence to emphasize that we are slightly abusing the
nomenclature, because the term QNM is normally used in the context of
stationary solutions, whereas the late-time, black brane solution of to us is
not stationary. Nevertheless, the fact that the time dependence of our
solution is only “along the spatial gauge theory directions” leads to many
familiar properties for excitations that are homogeneous along these
directions. The spectrum of such excitations for a different non-conformal
theory on de Sitter space was computed holographically in Buchel:2017lhu  via
a perturbative expansion around the conformal limit. The low-lying QNMs were
found to be purely imaginary, with the mysterious exception of the third mode,
which was found to have non-vanishing real and imaginary part. While we will
not attempt a direct calculation of the QNM spectrum of our model in this
paper, we will extract an estimate for the lowest QNM from the numerical
solution and confirm its purely dissipative nature.

Once the excess energy density as a function of time is known we can obtain
the pressure from the covariant conservation of the energy-momentum tensor
$\nabla^{i}\langle\Delta\hat{T}_{ij}\rangle=0$, which for the dS4 metric (14)
evaluates to

$\Delta\mathcal{P}(t)=-\Delta\mathcal{E}(t)-\frac{\Delta\mathcal{E}^{\prime}(t)}{3H}\,.$
(112)

At late times $\Delta{\mathcal{E}}(t)$ is well described by

$\Delta{\mathcal{E}}(t)=A\,e^{-i\omega t}\,,$ (113)

where $A$ is the amplitude of the fluctuation and $\omega$ is a purely
imaginary quasinormal mode with $\mathrm{Im}\omega<0$. By (112), this implies
that $\Delta{\mathcal{P}}(t)$ and $\Delta{\mathcal{E}}(t)$ satisfy at
sufficiently late times an EoS-like relation of the form

$\Delta{\mathcal{P}}(t)=w_{\rm eff}\Delta{\mathcal{E}}(t)$ (114)

with

$w_{\rm
eff}=-1+\frac{i\omega}{3H}=-1+\frac{\left(-\mathrm{Im}\omega\right)}{3H}\,.$
(115)

This effective EoS explains the constant late-time ratio
$3\Delta{\mathcal{P}}(t)/\Delta{\mathcal{E}}(t)$ in Fig. 5. We can extract an
estimate for $w_{\rm eff}$ from our numerical results at late times. In Fig. 7
we show two examples for this.

Figure 7: Decay of ${\mathcal{E}}+{\mathcal{P}}$ (solid blue) for $H=1/2$
(left) and $H=2$ (right), together with the QNM fit (dashed black) and their
absolute difference (dashed red).

We find $\omega=-3.398H\,i$ $(-2.782H\,i)$ for $H=1/2$ $(2)$ which gives
$3w_{\rm eff}=0.398$ $(-0.218)$ in precise agreement with the final values of
$3\Delta{\mathcal{P}}(t)/\Delta{\mathcal{E}}(t)$ shown in the middle plot in
the top (bottom) row of Fig. 5.

Since the dS expansion dilutes the energy density of the initial state, one
can think of the time evolution in dS as a dynamical implementation of the
gauge theory RG group. It may therefore seem surprising that, in general,
$w_{\rm eff}\neq 1/3$ at late times, given that our gauge theory is conformal
in the infrared when formulated in flat space. The reason is that, when the
theory is placed in a dS spacetime, the Hubble rate $H$ acts as an IR cut-off
Ghosh:2017big , in a way similar to the effect of placing the theory at finite
temperature. As a consequence, the IR behaviour of the gauge theory on dS is
approximately conformal only in the limits $H\ll M$ or $H\gg 1$. In the first
case the RG flow of the gauge theory in dS is similar to that in flat space,
in the sense that it explores almost all possible energy scales and is only
cut-off very close to the IR fixed point of the theory. This behaviour is
consistent with the top row of Fig. 5 (recall that $M=1$), where we see from
the late-time behaviour that $w_{\rm eff}$ grows towards 1 as we move from the
right plot ($H=1$) to the left plot ($H=1/3$). In the second case the IR cut-
off is so close to the UV fixed point that the entire RG flow is approximately
conformal. This is supported by the perturbative analysis of Buchel:2017lhu ,
which suggests that conformal symmetry is restored at late times in the limit
of large expansion rates $H\gg M$. We indeed find evidence for this behaviour
in Fig. 5: $w_{\rm eff}$ obtains its minimal value $w_{\rm eff}^{\rm
min}=-0.298$ for $H=1$ (top-right plot) and grows monotonically for $H>1$
(bottom row). Although we were not able to obtain numerical results for $H\gg
3$ because the numerics becomes increasingly challenging, we expect $w_{\rm
eff}$ to approach the value $1/3$ in the limit of large $H$.

To strengthen our conclusion that the late-time geometry is determined by
quasi-normal modes we show in Fig. 8 three evolutions for $H=1$ in which we
quench the boundary metric function $S_{0}(t)$ by a factor
$1+Ae^{-32(t-t_{q})^{2}}$ with amplitudes $A=0.05,0.01$ and $0.005$ at
quenching times $t_{q}=3.5,5$ and $7$, respectively.

Figure 8: Evolution of ${\mathcal{E}}+{\mathcal{P}}$ (left) and
$3\Delta{\mathcal{P}}/{\mathcal{E}}$ (right) for $H=1$ quenched with different
amplitude at different times. The quenches are realised by exciting the
boundary metric with Gaussian perturbations whose widths are indicated by the
red, green and blue bands. The quench performs work on the system, which leads
to an increase of ${\mathcal{E}}$ and explains the increase of
${\mathcal{E}}+{\mathcal{P}}$. After the quench the ratio of pressure and
energy density quickly returns to the late-time value of the unquenced
evolution (the small oscillations are verified to be numerical artefacts).

The quenches shown perform work on the system, which leads to an increase of
energy and, as a consequence, also to an increase of
${\mathcal{E}}+{\mathcal{P}}$. Quenching the late-time dynamics of the system
excites higher QNMs. The curves in Fig. 8 show the decay of these modes within
a time of order $1/H$, after which the system returns to a state whose late-
time dynamics is entirely characterized by $w_{\rm eff}$, i.e. by the first
QNM. In the next section we discuss the maximally symmetric late-time state in
more detail.

###  4.6 Late-time solution

At late times the geometry and the scalar field approach the following form
(see also Buchel:2017pto )

$ds_{\infty}^{2}=-A_{\infty}(r){\rm d}t^{2}+2{\rm d}r{\rm
d}t+e^{2Ht}S_{\infty}(r)^{2}{\rm
d}\vec{x}^{2}\,,\quad\phi=\phi_{\infty}(r)\,,$ (116)

where $A_{\infty}$, $S_{\infty}$ and $\phi_{\infty}$ are the time independent
late-time limits of the metric functions and the scalar field defined as

${A_{\infty}}(r)\equiv\lim\limits_{t\to\infty}A(t,r)\,,\quad{S_{\infty}}(r)\equiv\lim\limits_{t\to\infty}e^{-Ht}S(t,r)\,,\quad\phi_{\infty}(r)\equiv\lim\limits_{t\to\infty}\phi(t,r)\,.$
(117)

In the following we will drop the index “$\infty$” and implicitly assume this
limit in all appearances of $A$, $S$ and $\phi$. Under these conditions the
equations of motion simplify to

$\displaystyle 0$
$\displaystyle=A\,\phi^{\prime\prime}+\left(A^{\prime}+3\left(\frac{A\,S^{\prime}}{S}+H\right)\right)\phi^{\prime}-V^{\prime}\,,$
(118a) $\displaystyle 0$ $\displaystyle=\phi^{\prime
2}+\frac{3S^{\prime\prime}}{2S}\,,$ (118b) $\displaystyle 0$
$\displaystyle=A^{\prime\prime}+\frac{4\left(S^{\prime}\left(A^{\prime}+3H\right)+A\,S^{\prime\prime}\right)}{S}+\frac{2A\,S^{\prime
2}}{S^{2}}+2A\,\phi^{\prime 2}+4V\,,$ (118c) $\displaystyle 0$
$\displaystyle=A^{2}\left(\frac{2\left(S\,S^{\prime\prime}+S^{\prime
2}\right)}{S^{2}}+\frac{2\phi^{\prime
2}}{3}\right)+\frac{A\,S^{\prime}\left(A^{\prime}+8H\right)}{S}+H\left(2H-A^{\prime}\right)+\frac{4A\,V}{3}\,,$
(118d) $\displaystyle 0$
$\displaystyle=\frac{3\left(S^{\prime}\left(A^{\prime}+6H\right)+2A\,S^{\prime\prime}\right)}{2S}+\frac{3A\,S^{\prime
2}}{S^{2}}+A\phi^{\prime 2}+2V\,.$ (118e)

The late time geometry (116) possesses an event horizon located at
$r=r_{\mathrm{EH}}$ defined by the condition

$A(r_{\mathrm{EH}})=0\,.$ (119)

Using this condition in (118d) one finds that the surface gravity at the event
horizon equals the Hubble rate

$\kappa_{\mathrm{EH}}=\frac{1}{2}A^{\prime}(r_{\mathrm{EH}})=H\,.$ (120)

Because the geometry depends on time, the apparent horizon, i.e. the outermost
trapped lightlike surface, does not coincide with the event horizon. The
radial position of the apparent horizon $r_{\mathrm{AH}}$ is determined by
$\dot{S}|_{r_{\mathrm{AH}}}=0$, which in the coordinate system (117) gives the
condition

$0=S(r_{\mathrm{AH}})\frac{{\rm d}}{{\rm
d}t}e^{Ht}+\frac{1}{2}e^{Ht}A(r_{\mathrm{AH}})S^{\prime}(r_{\mathrm{AH}})\,.$
(121)

Using this condition in (118d) together with (118e) allows us to express
surface gravity at the apparent horizon

$\kappa_{\mathrm{AH}}=\frac{1}{2}A^{\prime}(r_{\mathrm{AH}})=-H\,.$ (122)

We arrive at the conclusion that surface gravity at the apparent horizon
equals minus surface gravity at the event horizon

$\kappa_{\mathrm{AH}}=-\kappa_{\mathrm{EH}}\,.$ (123)

Above, we used the ansatz (116) and the equations of motion to derive the
surface gravity of the event and of the apparent horizon. This ansatz can be
seen as educated guess motivated by our numerical results which at late times
agree with (116) very accurately. However, it is possible to arrive at (120)
without invoking the equations of motion or taking guidance from numerical
analysis. For this we use as starting point the assumption that the system
evolves towards a vacuum state that by definition obeys the symmetries of the
background geometry. On de Sitter space this state is known as Bunch–Davis
vacuum Bunch:1978yq . By the holographic duality the bulk metric dual to this
vacuum state has to satisfy the isometries of dS4 as well. A domain wall
parametrization in FG coordinates makes these isometries manifest666Solutions
of this type have been constructed in Ghosh:2017big  for a variety of gauge
theories.

$ds_{\mathrm{DW}}^{2}={A_{\infty}}(\rho)\left(-{\rm
d}t^{2}_{\mathrm{FG}}+e^{2Ht_{\mathrm{FG}}}{\rm d}\vec{x}^{2}\right)+{\rm
d}\rho^{2}\,,$ (124)

where the time and radial coordinate are related to those in (116) by

$\rho(r)=-\int^{\infty}_{r}{\rm
d}r^{\prime}{A_{\infty}}(r^{\prime})^{-1/2}\,,\quad
t_{\mathrm{FG}}(r,t)=t-\int^{\infty}_{r}{\rm
d}r^{\prime}{A_{\infty}}(r^{\prime})^{-1}\,.$ (125)

Demanding that (116) can be transformed to the manifestly dS4-isometric domain
wall form (124) gives the relation

${S_{\infty}}(r)=\sqrt{{A_{\infty}}(r)}\,e^{-H\int_{r}^{\infty}{\rm
d}r^{\prime}\,{A_{\infty}}(r^{\prime})^{-1}}\,.$ (126)

Beyond the event horizon ($r<r_{\mathrm{EH}}$) the metric function
$A_{\infty}(r)$ is negative and naively one could expect the right-hand side
of (126) to acquire an imaginary part. The exponent, however, diverges as $r$
approaches the event horizon, which leads to a term $a\log(r-r_{\rm EH})$ when
assuming a regular near horizon expansion of the form
$A_{\infty}(r)=a(r-r_{\rm EH})+{O\left((r-r_{\rm EH})^{2}\right)}$. The
imaginary parts exactly cancel and ${S_{\infty}}$ remains real inside the
event horizon if

${A_{\infty}}(r)=2H\left(r-r_{\mathrm{EH}}\right)+O\left(\left(r-r_{\mathrm{EH}}\right)^{2}\right)\,.$
(127)

This gives precisely the same result for
$\kappa_{\mathrm{EH}}=\lim\limits_{r\to
r_{\mathrm{EH}}}\frac{1}{2}A^{\prime}(r)$ as (120), only by demanding dS4
symmetry and regularity of the solution at the horizon. It is possible to
derive (122) in an analogous way as well.

In Fig. 9 we show surface gravity (left) and area density (right) of event and
apparent horizon for $H=1$ and $H=2.5$.

Figure 9: Surface gravity $\kappa$ (left) and area density (right) of event
and apparent horizons for two different values of the Hubble rate.

Initially, when the geometry is static and dual to a thermal state on flat
space, the event and apparent horizon coincide and therefore have equal
surface gravity and area density that depends on the choice of the initial
temperature. On the expanding background the locations of apparent and event
horizon deviate and so do the respective surface gravities and area densities.
In accordance with the analytic analysis, the numeric evolution evolves
towards a solution where surface gravity of event and apparent horizon
precisely satisfy $\kappa_{\mathrm{EH}}=-\kappa_{\mathrm{AH}}=H$. For the
gravity dual of $\mathcal{N}=4$ SYM theory on de Sitter space apparent and
event horizon densities can be straightforwardly computed because the bulk
geometry is known explicitly Apostolopoulos:2008ru ; Buchel:2016cbj . This
geometry has no apparent horizon but an event horizon whose area density is
given by

$S(r_{\rm
EH},\,t)^{3}=8H^{3}S_{0}(t)^{3}=8\kappa_{\mathrm{EH}}^{3}S_{0}(t)^{3}\,.$
(128)

For comparison, a thermal state in $\mathcal{N}=4$ SYM on flat space has
$S(r_{\rm EH})^{3}=\kappa_{\mathrm{EH}}^{3}/8$.

Indeed, we show in Fig. 9 (right) that as $H$ increases from $1$ to $2.5$ and
thereby approaching the conformal UV fixed point ($M/H\to 0$), the apparent
horizon area density ($S(r_{\rm hor},t)^{3}/S_{0}(t)^{3}$) decreases and the
normalised event horizon area density approaches the conformal value at late
times. The properties of the late time solution are solely determined by the
value of $H$ and $M$, but the way they are approached depends on the initial
conditions. Thermal initial states with large entropy, and therefore with
large initial event and apparent horizons area, lead to decreasing area
densities of both, apparent and event horizons. An example for this is the
simulation for $H=1$ shown in Fig. 9 (right). However, if the initial entropy
is such that the corresponding horizon area is smaller than the area density
of the late-time horizon, the area density grows. This is precisely what we
find for example for the event horizon area density for $H=2.5$ (dashed black)
shown in Fig. 9 (right). It may seem surprising that the event horizon area
density increases even though the surface gravity decreases substantially.
This, however, can qualitatively be explained by comparing the analytic area
densities of $\mathcal{N}=4$ SYM in flat space to de Sitter space, where
indeed at fixed surface gravity the area density is much smaller on flat
space. The area density of the apparent horizon (blue) decreases monotonically
in all cases shown.

This is a clear indication that the holographic interpretation of bulk
horizons areas as entropy in the dual field theory is subtle when the boundary
theory is expanding. We will elaborate on this in the discussion. However, we
emphasize that the comoving area density of the apparent horizon $S(r_{\rm
hor},t)^{3}$ is a monotonously growing function of time in accordance with
Hawinkgs area theorem hawking1972  and the more recent discussion777We are
grateful to Alex Buchel for bringing this to our attention. Buchel:2017pto ;
Buchel:2019pjb  in the context of a holographic model similar to the one we
use here.

In the next section we analyse the entanglement properties of the de Sitter
vacua constructed in this section and comment on the subtleties involved in
the assignment of an effective entropy to their horizon area densities.

##  5 Entanglement and horizon entropies

The holographic duality maps the area of event horizons in time independent
bulk geometries to the thermodynamic entropy of thermal states in the dual
field theory. In time-dependent geometries this mapping is obscured, firstly
because the dual field theory is not in thermal equilibrium and thermodynamic
concepts of entropy and temperature do not apply and secondly because the
mapping of the horizon to the boundary is not necessarily unique. To study the
time evolution of non-equilibrium states in the field theory it can
nevertheless be useful to define an effective temperature in terms of surface
gravity of the event horizon888We will also look at the apparent horizon, even
though the physical interpretation of temperatures defined in terms of
apparent horizons is problematic, because it in general depends on the slicing
of spacetime..

In time independent geometries event and apparent horizon coincide and the
entropy density in the dual field theory is uniquely defined in terms of the
area density of the event horizon in the bulk. Because de Sitter space is
expanding it is non-trivial to map points in the boundary to points on the
horizon in the bulk whose area density changes due to the exponentially
growing scale factor (see also Buchel:2019qcq  for related difficulties in
interpreting the apparent horizon area).

A robust and gauge independent notion of entropy is given by entanglement
entropy Holzhey:1994we  of spatial subregions $\mathcal{R}$ in the QFT defined
as

$S_{\mathcal{R}}=-\mathrm{Tr}_{\mathcal{R}}\hat{\rho}_{\mathcal{R}}\log\hat{\rho}_{\mathcal{R}}\,,$
(129)

where $\hat{\rho}_{\mathcal{R}}=\mathrm{Tr}_{\bar{\mathcal{R}}}\hat{\rho}$
denotes the reduced density matrix obtained by performing on the full density
matrix $\hat{\rho}$ a partial trace over the degrees of freedom outside
$\mathcal{R}$. For simplicity we will assume these subregions to be spatial
balls at some fixed time $t=t_{0}$ with radius $\ell$

$\mathcal{R}=\\{t=t_{0},0\leq
r\leq\ell,0\leq\theta\leq\pi,0\leq\varphi<2\pi\\}\,.$ (130)

The coordinate $r$ in this section is the radial coordinate at the boundary
and should not be confused with the holographic coordinate. Since the density
matrix $\hat{\rho}$ can be time-dependent, i.e. defined in terms of time-
dependent states, entanglement entropy is also well defined for states that
are not in thermal equilibrium. For $H=0$ and $\ell\to\infty$ the entangling
region $\mathcal{R}$ covers an entire spacelike slice of Minkowski space. In
this limit $\hat{\rho}_{\mathcal{R}}=\hat{\rho}$ and (129) equals the von
Neumann entropy of the full density matrix, i.e., the thermodynamic entropy of
a quantum state in thermal equilibrium. For $H\neq 0$ the QFT inherits the
causal structure of de Sitter background which has a cosmological event
horizon located at $r=H^{-1}$. Although spatial regions of size $\ell>H^{-1}$
are causally disconnected, quantum states on such regions can be entangled
Maldacena:2012xp .

Direct field theory computations of entanglement entropy are only possible in
exceptional cases like for example in two-dimensional CFTs Calabrese:2004eu
and in dimensions higher than two only in non-interacting QFTs
Srednicki:1993im . The holographic duality replaces the field theory
computation of entanglement entropy by a much simpler extremisation problem
for the surface area $\mathcal{A}_{\mathcal{R}}$ of a codimension two surface,
homologous to $\mathcal{R}$, in the bulk theory Ryu:2006bv ; Hubeny:2007xt

$S_{\mathcal{R}}=\frac{\mathcal{A}_{\mathcal{R}}}{4G}\,.$ (131)

The boundary of the relevant surface coincides with the boundary of the
entangling region $\mathcal{R}$ in the field theory and extremises the area
functional in the bulk theory

$\mathcal{A}_{\mathcal{R}}[X]=\int{\rm
d}^{3}\sigma\sqrt{\mathrm{Det}\left(\partial_{a}X^{\mu}\partial_{b}X^{\nu}g_{\mu\nu}\right)}\,.$
(132)

The surface embedding $X^{\mu}=X^{\mu}(\sigma^{a})$ is parametrised by three
intrinsic coordinates for which we choose $\sigma^{a}=\\{r,\theta,\varphi\\}$.
The entangling regions (130) do not break spherical symmetry in the boundary
theory. We can then parametrise the bulk surface with

$X^{\mu}(z)=\\{z(r),t(r),r,\theta,\varphi\\}\,.$ (133)

This choice simplifies the area functional considerably. Integration over the
angular coordinates $\theta$ and $\phi$ can be performed explicitly and the
remaining expression takes the form of a geodesic action

$\mathcal{A}_{\mathcal{R}}[X]=4\pi\int{\rm
d}r\sqrt{\bar{g}_{\alpha\beta}(z(r),t(r))\frac{{\rm d}X^{\alpha}}{{\rm
d}r}\frac{{\rm d}X^{\beta}}{{\rm d}r}}\quad\mathrm{s.t.}\quad
X^{\mu}(0)=\\{0,t_{0},\ell,\theta,\varphi\\}\,,$ (134)

where the metric $\bar{g}_{\alpha\beta}$ is related by a conformal factor to a
three dimensional subspace ($\alpha,\beta=\\{t,z,r\\}$) of the bulk metric
(116) (see also Ecker:2015kna )

$\bar{ds}^{2}=\bar{g}_{\alpha\beta}{\rm d}x^{\alpha}{\rm
d}x^{\beta}=\left(r\,S(z)\,a(t)\right)^{4}g_{\alpha\beta}{\rm d}x^{\alpha}{\rm
d}x^{\beta}\,.$ (135)

where we use $a(t)\equiv e^{Ht}$ in the following to simplify notation. The
equations of motion that follow from $\delta\mathcal{A}_{\mathcal{R}}=0$ take
the form of a non-affine geodesic equation

$\frac{{\rm d}^{2}X^{\alpha}}{{\rm
d}r^{2}}+\Gamma^{\alpha}_{\beta\gamma}\frac{{\rm d}X^{\beta}}{{\rm
d}r}\frac{{\rm d}X^{\gamma}}{{\rm d}r}=J\frac{{\rm d}X^{\alpha}}{{\rm d}r}\,,$
(136)

where $\Gamma^{\alpha}_{\beta\gamma}$ is the Levi-Civitá connection associated
to $\bar{g}_{\alpha\beta}$ and is meant to be evaluated at the location of the
surface $X^{\alpha}(r)$; the viscous friction term on the right hand side
includes the Jacobian $J=\frac{{\rm d}^{2}\tau(r)}{{\rm d}r^{2}}/\frac{{\rm
d}^{2}\tau(r)}{{\rm d}r^{2}}$ that originates from transforming from the
affine parameter $\tau$ defined by $\frac{{\rm d}X^{\alpha}(\tau)}{{\rm
d}\tau}\frac{{\rm d}X^{\beta}(\tau)}{{\rm d}\tau}\bar{g}_{\alpha\beta}=1$ to
the non-affine parameter $r$.

For $\phi=0$ the field theory is $\mathcal{N}=4$ SYM theory and (136) has a
simple analytic solution $z(r)=\sqrt{\ell^{2}-r^{2}}$ and
$t(r)=-\log\left(1+z(r)\right)$ in a gauge where $A(z)=z^{-2}-1$. In this case
the de Sitter vacuum can be mapped by a conformal transformation to the
Minkowski vacuum Buchel:2017pto  and the holographic entanglement entropy is
equal to the area of a hemisphere in AdS5 Ryu:2006bv .

For the non-conformal case with $\phi\neq 0$ the geodesic equations are too
long to display here and no closed solutions are available. At the turning
point of the bulk surface located at $z_{\ast}\equiv z(r=0)$ the equations
simplify by symmetry $z^{\prime}(r=0)=t^{\prime}(r=0)=0$ to

$\displaystyle t^{\prime\prime}(r)$
$\displaystyle=-3z(r)^{2}a(t(r))^{2}S(z(r))S^{\prime}(z(r))$ (137a)
$\displaystyle z^{\prime\prime}(r)$
$\displaystyle=3z(r)^{2}a(t(r))S(z(r))\left[z(r)^{2}a(t(r))A(z(r))S^{\prime}(z(r))-S(z(r))a^{\prime}(t(r))\right]\,.$
(137b)

The sign of $z^{\prime\prime}(r=0)$ determines if the bulk surface can reach
the boundary ($z^{\prime\prime}(r=0)<0$) or not ($z^{\prime\prime}(r=0)>0$).
The condition $z^{\prime\prime}(r)=0$ defines the ‘entanglement horizon’,
i.e., a barrier in the bulk that extremal surfaces attached to the boundary do
not cross whose location is determined by

$0=z^{2}a(t)A(z)S^{\prime}(z)-S(z)a^{\prime}(t)=z^{2}A(z)S^{\prime}(z)-HS(z)\,.$
(138)

Plugging this relation into the definition of $\dot{S}$ gives

$\dot{S}+\frac{1}{2}z(r)^{2}A(z(r))a(t(r))S^{\prime}(z(r))=0\,.$ (139)

Except for the second term, this relation is equal to equation (121) that
determines the location of the apparent horizon. Since
$z(r)^{2}A(z(r))S^{\prime}(z(r))$ is positive and monotonic in $z(r)$ the
entanglement horizon is located between event and apparent horizon. For
boundary geometries that are de Sitter the entanglement horizon is a
Lagrangian surface of the bulk geometry, or in other words a surface with zero
surface gravity (see also Engelhardt:2013tra  for similar results on extremal
surface barriers). This can be shown by combining (138) with the Einstein
equations (118b) to (118e) which gives

$\kappa_{\rm Ent}=-\frac{z^{2}}{2}A^{\prime}(z)=0\,.$ (140)

Furthermore, when the geometry is static ($H=0$), we know from (121) that the
apparent horizon ($\dot{S}=0$) and event horizon ($A=0$) coincide. From (139)
we then see that this also solves the entanglement horizon equation, so that
all three horizons coincide as expected Hubeny:2012ry .

We compute the entanglement entropy numerically by shooting, i.e., by
integrating (136) from the turning point $z_{\ast}$ of the geodesics to a
cutoff-value for the holographic coordinate
$z_{\mathrm{cut}}=\epsilon/(1-\xi\epsilon)$ close to the boundary located at
$z=0$. Note that the gauge freedom $\xi$ in $S(r)=r+\xi+\mathcal{O}(1/r)$ near
the boundary leads to an additional $1/\epsilon$ dependent divergence in the
entanglement entropy. We eliminate the $1/\epsilon$ divergence with the gauge
transformation $r\rightarrow r-\xi$ that ensures $S(r)=r+\mathcal{O}(1/r)$ and
in addition fixes the residual gauge freedom ($r\rightarrow r+\xi$) in the
radial coordinate. The cutoff regulated value for the entanglement entropy can
then be obtained by numerically solving the following integral

$S_{\mathcal{R}}^{\mathrm{cut}}=\frac{\pi}{G}\int_{z_{\mathrm{cut}}}^{z_{\ast}}{\rm
d}z\left(S(z)a(t(z))r(z)\right)^{2}\sqrt{-A(z)t^{\prime}(z)^{2}-\frac{2t^{\prime}(z)}{z^{2}}+\left(S(z)a(t(z))r^{\prime}(z)\right)^{2}}\,.$
(141)

In Fig. 10 we show extremal surfaces in the gauge where the apparent horizon
is fixed at $z_{\rm AH}=1/2$, obtained by shooting from different values of
$z_{\ast}$, for several spherical entangling regions of different radius
$\ell$.

Figure 10: Extremal surfaces with various values of $z_{\ast}$ in the $z-r$
plane (left) and $t-r$ plane (right) for $H=1$ (top) and $H=2.5$ (bottom).
Solid, dotted and dashed black lines indicate the radial location of event,
entangling and apparent horizon, respectively. Red dashed curves are extremal
surfaces with $z_{\ast}=z_{\mathrm{EH}}$ which end on the boundary precisely
at the location ($r=1/H$) of the cosmological horizon.

Interestingly, surfaces with $z_{\ast}=z_{\mathrm{EH}}$ end on the boundary
precisely at $r=1/H$ where the cosmological horizon is located. The red dashed
curves in Fig. 10 are two examples for such surfaces with
$z_{\ast}=0.3597\,(0.1961)$ for $H=1\,(2.5)$. Extremal surfaces of entangling
regions larger than the cosmological horizon ($\ell>1/H$) probe regions behind
the event horizon (solid black), but never beyond the entanglement horizon
(black dotted). In our gauge, where the apparent horizon (black dashed) is
located at $z_{\mathrm{AH}}=1/2$, the entanglement horizon is located at
$z=0.4554\,(0.3434)$ and $z_{\rm EH}=0.3597\,(0.1961)$ for $H=1\,(2.5)$. For
comparison, in the gravity dual of $\mathcal{N}=4$ SYM theory on de Sitter the
event horizon is located at $z_{\mathrm{EH}}=1/(2H-\xi)$ and extremal surfaces
can probe until $z=1/(H-\xi)$ in the gauge where $z_{\rm AH}=-1/\xi$.

In time independent geometries the apparent, event and entanglement horizon
coincide and bound the region that is causally connected to the boundary, as
well as the region that can be probed by extremal surfaces anchored at the
boundary. For de Sitter we find that extremal surfaces can penetrate the event
horizon, but only if their entangling region is super-horizon, i.e. larger
than the cosmological horizon. This implies that only a ‘super-observer’ with
access to information in a region larger than the observable universe could in
principle be able to reconstruct the dual spacetime behind the event horizon
from field theory data.

A physical explanation for this phenomenon is that the extremal surface
corresponding to the observable universe in the boundary is itself a
cosmological horizon in the bulk, this time of an observer at the origin on
the boundary. This is illustrated in Fig. 11, where we shoot a family of null
geodesics (blue) from a generic point on the bulk extremal surface (dashed
red). This family is indeed just able to reach the origin at the boundary,
which asserts our statement that this point is an element of the (bulk)
cosmological horizon of an observer at the origin. The green curves on the
other hand represent a family of null geodesics that originate from a generic
point beyond the cosmological horizon and are therefore unable to reach the
origin. The starting point of this second family of geodesics is then element
of a cosmological horizon in the bulk for an observer located at $Hr=0.4$. We
stress that the origin is not a special place, but only defined as the origin
of the entangling region used in this example.

Fig. 12 shows the areas associated to the extremal surfaces for various values
of the cut-off $\epsilon$. Asymptotically the area equals Maldacena:2012xp

$\mathcal{A}=4\pi\ell^{2}\left(\frac{1}{2\epsilon^{2}}+\frac{1}{3}\log(\epsilon)+O(1)\right)+O(\log(\ell)).$
(142)

The leading divergence is clear from the left figure, whereas the subtracted
version (right) shows the difference with the leading order divergence.
Unfortunately it is numerically difficult to extract higher order coefficients
that depend on the non-conformality of the theory. Naively from Fig. 10 it may
seem that for large $\ell$ we should obtain a volume law scaling, by the usual
argument that the extremal surface hovers at the entanglement horizon and
hence gets a contribution proportional to $S(z_{\rm Ent},\,t)^{3}$ times the
volume of the region. For time independent geometries this is indeed the case,
but interestingly our extremal surfaces have a non-trivial time dependence
through $S(z,\,t)=S_{0}(z)a(t(r))$ and $t(r)$. In fact, as illustrated in Fig.
10 and more explicitly in Fig. 13 the time as the surface hovers the horizon
is for $H=1$ approximately given by $t_{\rm Ent}\approx-0.3-\log(\ell)$. In
combination with $a(t)=e^{Ht}$ and $S_{0}(z_{Ent})\approx 0.69$ this implies
that the part of the surface999In the following we assume that this is the
entire surface, which is valid for $\ell\gg 1/H$. hovering the entanglement
horizon has an area $\frac{4}{3}\pi\ell^{3}0.134/\ell^{3}\approx 0.559$ which
is only a constant term as opposed to the usual volume scaling of the time
independent setting.

Figure 11:  Blue and green curves are two families of null geodesics,
originating on a generic point of the extremal surface with
$z_{\ast}=z_{\mathrm{EH}}$ and a generic point in the bulk that is not
enclosed by the extremal surface. The fact that this family just reaches the
observer at the origin implies that the extremal surface is the (bulk)
cosmological horizon of an observer at the origin ($r=0$).

Figure 12: Entanglement entropy for our $H=1.0$ (left), rescaled by the area
and cut-off $\epsilon^{2}$, together with the leading order divergence
subtracted (right), for different values of the cut-off. Figure 13:  Minimum
time reached on the extremal surface (see also Fig. 10) as a function of
length. At large $\ell$ this time grows logarithmically with $\ell$.

##  6 Discussion

We have used a holographic model to study the dynamics of a strongly coupled
non-conformal gauge theory in four-dimensional de Sitter space. The four-
dimensional dS metric is prescribed a priori and is probed by the strongly
coupled matter. In other words, the five-dimensional gravitational model
provides a dual description of the dynamics of the gauge theory matter but not
of the four-dimensional gravitational field on which this matter propagates.

We have carefully explained the holographic renormalisation of the model and
the anomalies that arise in the dual field theory due to the curved boundary
geometry and the scalar field. After reviewing thermodynamic and transport
properties of the model in flat space, we have presented numerical results for
the fully non-linear time evolution of finite-temperature states towards the
Bunch-Davis vacuum at late times. We have shown that the approach to the de
Sitter vacuum is characterised by an effective relation of the form
${\mathcal{P}}=w{\mathcal{E}}$. This is different from the equilibrium
equation of state in flat space, so much so that actually $w$ is negative if
the ratio $H/M$ is close to unity.

We have studied in detail the properties of the de Sitter vacua of the
holographic model and we have analyzed the different horizons that arise in
the bulk geometry. The connection between event horizons and thermodynamics
found for black holes Hawking:1974sw  also applies to cosmological horizons
PhysRevD.15.2738 . Therefore, in analogy with the Bekenstein-Hawking
temperature of black holes PhysRevD.7.2333 ; Hawking:1974sw , an observer
living at the boundary would associate a temperature to the cosmological
horizon of de Sitter space

$T_{\mathrm{dS}}=\frac{\kappa_{\mathrm{dS}}}{2\pi}=\frac{H}{2\pi}\,,$ (143)

where the surface gravity $\kappa_{\mathrm{dS}}$ is evaluated at the horizon
and equals the Hubble rate $H$. Interestingly, we found two temperatures in
our five-dimensional dual description: one at the bulk event horizon, equal to
the result by Hawking and Gibbons of $H/2\pi$, and another temperature at the
deeper apparent horizon, equal to $-H/2\pi$ (also found in Buchel:2019qcq ).
In the literature there are several works Klemm:2004mb ; Visser:2019muv ;
Jacobson:2019gco  suggesting such a negative temperature based on the first
law. In our case the apparent horizon is however causally disconnected from
the boundary as well as time-slicing dependent, which suggests that indeed the
positive temperature of the event horizon is the physical temperature.

In analogy with the Bekenstein-Hawking law, which relates the entropy of a
black hole to the area of its event horizon, one can also associate a
gravitational entropy to de Sitter space

$S_{\mathrm{dS}}=\frac{A_{\mathrm{dS}}}{4G_{4}}=\frac{\pi}{H^{2}G_{4}}\,,$
(144)

where $A_{\mathrm{dS}}=4\pi/H^{2}$ is the area of the cosmological horizon and
$G_{4}$ is the Newton’s constant of the four-dimensional boundary theory. Note
that this entropy is formally infinite in our case since the boundary metric
is non-dynamical and hence implicitly we are setting $G_{4}=0$. Therefore the
gravitational entropy (144) should not be confused with the entropy that one
may want to assign to the area densities of the event, entanglement and
apparent horizons that we studied in sections (4.6) and (5). In principle
these would be related to the entropy of the matter in de Sitter space, but
this relation is not straight-forward. For a stationary geometry the three
areas agree and can be identified with the entropy density of the boundary
gauge theory. Firstly, in our expanding geometry the horizons do not coincide
with one another, and furthermore the mapping between points at the horizon
and at the boundary is ambiguous. Secondly, such an entropy density
interpretation suggests a volume law, whereas we showed in Fig. 12 in
combination with Fig. 13 that for large regions the entanglement horizon
contribution to the entanglement entropy is just a constant term. The
divergent piece of the entanglement entropy satisfies an area law, which
prohibits a direct extraction of the IR part of the entropy (see also
Maldacena:2012xp ). It is hence difficult to have a direct interpretation of
the entropy of de Sitter itself, but we note that in a theory with dynamical
gravity it is conjectured that this entropy is limited by a Bekenstein-Hawking
term of $A/4G_{4}$. This entropy, or part thereof, can then potentially be
identified with the entanglement entropy whereby $1/G_{4}$ plays the role of
the UV cut-off Susskind:1994sm .

A further result of our analysis of the entanglement entropy of boundary
regions was that the extremal surfaces corresponding to entangling regions
that coincide with the boundary observable universe exactly touch the bulk
event horizon and are in fact a bulk cosmological horizon. It would be
interesting to understand analytically why the extremal surface associated
with the boundary cosmological horizon is itself a bulk cosmological horizon.

Our analysis of perturbations around the late-time state showed that, after a
quench, the state relaxes within a time $1/T\sim 1/H$ (see Fig. 7 and 8), in
agreement with Buchel:2017lhu . In our case this time can be a parametrically
different from $\Delta\mathcal{E}^{1/4}$. This hence gives further credibility
that the de Sitter temperature provides a physical temperature. On the
holographic side this can be understood by the fact that the relaxation time
is determined by the distance between the boundary and the event horizon,
which is indeed proportional to $1/H$.

For our late-time solution the energy density excess over the asymptotic late-
time solution decreases exponentially. As a consequence, at the time when the
negative excess pressure becomes relevant the energy density excess
$\Delta\mathcal{E}$ will quickly become smaller than $T^{4}$, with $T\sim H$
the background de Sitter temperature (note that this temperature is the
minimal temperature, accelerating observers will see an even higher
temperature Narnhofer:1996zk ). This raises the question of whether the energy
and pressure excesses in this regime can be measured by an actual observer,
since the relevant modes will have wavelengths larger than the observable
Universe.

###### Acknowledgements.

We thank Alex Buchel and Michal Heller for their comments on the manuscript.
JCS and DM are supported by grants FPA2016-76005-C2-1-P, FPA2016-76005-C2-2-P,
2014-SGR-104, 2014-SGR-1474, SGR-2017-754, MDM-2014-0369,
PID2019-105614GB-C21, PID2019-105614GB-C22. They also acknowledge financial
support from the State Agency for Research of the Spanish Ministry of Science
and Innovation through the “Unit of Excellence María de Maeztu 2020-2023”
award to the Institute of Cosmos Sciences (CEX2019-000918-M).

## References

  * (1) J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal and U. A. Wiedemann, _Gauge/String Duality, Hot QCD and Heavy Ion Collisions_. Cambridge University Press, 2014, 10.1017/CBO9781139136747, [1101.0618].
  * (2) W. Busza, K. Rajagopal and W. van der Schee, _Heavy Ion Collisions: The Big Picture, and the Big Questions_ , _Ann. Rev. Nucl. Part. Sci._ 68 (2018) 339 [1802.04801].
  * (3) M. P. Heller, R. A. Janik and P. Witaszczyk, _The characteristics of thermalization of boost-invariant plasma from holography_ , _Phys. Rev. Lett._ 108 (2012) 201602 [1103.3452].
  * (4) P. M. Chesler and L. G. Yaffe, _Holography and colliding gravitational shock waves in asymptotically AdS 5 spacetime_, _Phys. Rev. Lett._ 106 (2011) 021601 [1011.3562].
  * (5) J. Casalderrey-Solana, M. P. Heller, D. Mateos and W. van der Schee, _From full stopping to transparency in a holographic model of heavy ion collisions_ , _Phys. Rev. Lett._ 111 (2013) 181601 [1305.4919].
  * (6) Y. Aoki, G. Endrodi, Z. Fodor, S. Katz and K. Szabo, _The Order of the quantum chromodynamics transition predicted by the standard model of particle physics_ , _Nature_ 443 (2006) 675 [hep-lat/0611014].
  * (7) P. Schwaller, _Gravitational Waves from a Dark Phase Transition_ , _Phys. Rev. Lett._ 115 (2015) 181101 [1504.07263].
  * (8) C. Caprini et al., _Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions_ , _JCAP_ 04 (2016) 001 [1512.06239].
  * (9) G. D. Kribs and E. T. Neil, _Review of strongly-coupled composite dark matter models and lattice simulations_ , _Int. J. Mod. Phys. A_ 31 (2016) 1643004 [1604.04627].
  * (10) S. Tulin and H.-B. Yu, _Dark Matter Self-interactions and Small Scale Structure_ , _Phys. Rept._ 730 (2018) 1 [1705.02358].
  * (11) A. Strominger, _The dS / CFT correspondence_ , _JHEP_ 10 (2001) 034 [hep-th/0106113].
  * (12) P. McFadden and K. Skenderis, _Holography for Cosmology_ , _Phys. Rev. D_ 81 (2010) 021301 [0907.5542].
  * (13) K. Koyama and J. Soda, _Strongly coupled CFT in FRW universe from AdS / CFT correspondence_ , _JHEP_ 05 (2001) 027 [hep-th/0101164].
  * (14) D. Marolf, M. Rangamani and M. Van Raamsdonk, _Holographic models of de Sitter QFTs_ , _Class. Quant. Grav._ 28 (2011) 105015 [1007.3996].
  * (15) K. Ghoroku and A. Nakamura, _Holographic Friedmann equation and N=4 supersymmetric Yang-Mills theory_ , _Phys. Rev. D_ 87 (2013) 063507 [1212.2304].
  * (16) A. Buchel, M. P. Heller and J. Noronha, _Entropy Production, Hydrodynamics, and Resurgence in the Primordial Quark-Gluon Plasma from Holography_ , _Phys. Rev. D_ 94 (2016) 106011 [1603.05344].
  * (17) A. Buchel, _Ringing in de Sitter spacetime_ , _Nucl. Phys. B_ 928 (2018) 307 [1707.01030].
  * (18) A. Buchel and A. Karapetyan, _de Sitter Vacua of Strongly Interacting QFT_ , _JHEP_ 03 (2017) 114 [1702.01320].
  * (19) A. Buchel, _Entanglement entropy of ${\cal N}=2^{*}$ de Sitter vacuum_, _Nucl. Phys. B_ 948 (2019) 114769 [1904.09968].
  * (20) A. Buchel, _$\chi\rm{SB}$ of cascading gauge theory in de Sitter_, _JHEP_ 05 (2020) 035 [1912.03566].
  * (21) P. S. Apostolopoulos, G. Siopsis and N. Tetradis, _Cosmology from an anti-de sitter–schwarzschild black hole via holography_ , _Physical Review Letters_ 102 (2009) .
  * (22) M. Attems, J. Casalderrey-Solana, D. Mateos, I. Papadimitriou, D. Santos-Oliván, C. F. Sopuerta et al., _Thermodynamics, transport and relaxation in non-conformal theories_ , _JHEP_ 10 (2016) 155 [1603.01254].
  * (23) 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 (2001) 595 [hep-th/0002230].
  * (24) M. Bianchi, D. Z. Freedman and K. Skenderis, _How to go with an RG flow_ , _JHEP_ 08 (2001) 041 [hep-th/0105276].
  * (25) M. Bianchi, D. Z. Freedman and K. Skenderis, _Holographic renormalization_ , _Nucl. Phys. B_ 631 (2002) 159 [hep-th/0112119].
  * (26) S. S. Gubser and A. Nellore, _Mimicking the QCD equation of state with a dual black hole_ , _Phys. Rev._ D78 (2008) 086007 [0804.0434].