Strong-coupling dynamics and entanglement in de Sitter space

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 AdS$_5$. 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.


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. [1] 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 [2]). 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 [3][4][5]. During its subsequent evolution the fireball expands, cools down and eventually hadronises. In QCD in equilibrium, this transition is realised as a smooth crossover [6].
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 [7,8]). Finally, it has recently been suggested that dark matter may be strongly self-interacting (recent reviews include [9,10]), 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. [11,12], 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. [13][14][15][16][17][18][19][20]. 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 [21].
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.

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 [22]. 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 Here G is the five-dimensional Newton's constant, R[g] is the Ricci scalar associated to the five-dimensional bulk metric g µν on M, γ ij is the metric induced on a four-dimensional slice near the boundary ∂M, and is the trace of the extrinsic curvature K ij associated to this slice. The second term on the right-hand side of (2.1) is the familiar Gibbons-Hawking term. The third term in (2.1) will be described below. The equations of motion take the form Rg µν = 2∂ µ φ∂ ν φ − 2V (φ)g µν − (∂φ) 2 g µν , (2.3a) The potential V (φ) 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. [22] we therefore choose the potential which can be derived from the superpotential L is a length scale. The dimensionless constant φ 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 φ M = 2 . (2.7) Both V (φ) and W (φ) have a maximum at φ = 0 and a minimum at φ = φ M . Each of these extrema yields an AdS solution of the equations of motion with constant φ and radius L 2 = −3/V (φ). 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 N = 4 SYM with N colours we would have 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 (2.4) leads to three important properties of the model: First, the resulting geometry is asymptotically AdS 5 in the UV with radius L, since V (0) = −3/L 2 . Second, the second derivative of the potential at φ = 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,Ô, with scaling dimension ∆ 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 φ = φ M is again AdS 5 with a different radius In this region the effective mass of the scalar field differs from its UV value and it is given by 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 ρ. This has dimensions of (length) 2 and in terms of it the near-boundary metric takes the form The boundary is located at ρ = 0 and is parametrised by the coordinates x i with i = 0, . . . , 3.
Near the boundary the metric and the scalar field take the form 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 dS 4 metric with Hubble rate H: (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: The first term of the action (2.1) suffers from large-volume divergences, as can be verified by substituting the expansions (2.13) into the action. These divergences can be regularised and renormalised by a procedure called holographic renormalisation (see e.g. [23][24][25]), which makes the action finite and the variational principle well-defined. This procedure is implemented by including in (2.1) the counterterm action defined on a timelike, constant-ρ hypersurface near the boundary. The induced metric on this hypersurface is denoted γ ij and R[γ] is the associated Ricci scalar. The second term of (2.1) is also understood to be evaluated on this slice at ρ, the first term of (2.1) is understood to be evaluated by integrating down to this slice, and the limit ρ → 0 is understood to be taken at the end of the calculation. In (2.16), A(γ ij , φ) is the so-called conformal anomaly, which in our case is given by is the holographic gravitational conformal anomaly and is the conformal anomaly due to matter. In these equations all the terms are functionals of the metric γ ij and of the scalar field φ induced on the ρ-hypersurface. However, making use of the expansions (2.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 (2.20) For this reason we will often think of the anomaly as evaluated on the boundary values of the fields, in which case (2.18) and (2.19) reduce to The fact that √ −γA yields a finite result has two consequences. First, it means that the logarithmic term in (2.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 (2.18) and (2.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 α in (2.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 γ ij in such a way that their overall mass dimension is four. The βφ 4 term is an example of such a term. Other possible terms include combinations of R ij R ij , R 2 and φ 2 that are linearly independent of A. Therefore we could replace the last term of (2.16) by (2.23) The freedom to add these terms with arbitrary coefficients δ i is part of the general freedom in the choice of renormalisation scheme. The coefficient α plays a special role since it can be shifted by a scale transformation, which is implemented via the following rescaling of the coordinates where λ 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 λ)A, which in turn can be absorbed through the redefinition α → α + log λ. The freedom to rescale ρ, or equivalently to shift α, 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 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 β = 1/4φ 2 M is special because in this case the βφ 4 term combines with the second and the third summands in the first term of (2.16) to give precisely the superpotential (2.5). This means that, if the theory (2.1) is the bosonic truncation of a supersymmetric theory with superpotential W , then in flat space this choice of β leads to a supersymmetric renormalisation scheme. Motivated by this discussion, in this paper we will set to zero all the coefficients in (2.23) but α and β. This implies no loss of generality since physically meaningful quantities are scheme-independent.
Substitution of the expansions (2.13) in the equations of motion (2.3) determines several coefficients [25]. The Klein-Gordon equation for the scalar field fixes the logarithmic coefficient ψ (2) in terms of g (0)ij and φ (0) as 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 The logarithmic part at subleading order fixes is the purely gravitational part. For conformally flat metrics, such as (2.14), this part vanishes, hence The subleading non-logarithmic part of Einstein's equations fixes the trace of g (4)ij : as well as its covariant divergence The coefficients in (2.13) determine the holographic stress tensor as In this and in subsequent equations we have made use of (2.8) to replace G in favour of N , which makes the expected N 2 -scaling of the stress tensor manifest. The contributions T g ij and T φ ij come from the variation of A g and A φ in (2.16), respectively, and are given by: As we mentioned above, the first equation follows from the conformal flatness of the dS metric (2.14). The expectation value of the scalar operator in the field theory is given by In the presence of external sources the holographic stress tensor satisfies anomalycorrected Ward identities. These can be obtained from the variation of the renormalised on-shell action Invariance of the action under diffeomorphisms leads to the diffeomorphism Ward identity give the anomaly-corrected conformal Ward identity where d is the spacetime dimension of the boundary theory and A g and A φ are given in (2.18) and (2.19), respectively. In our case the Ward identities reduce to 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 (2.24). In the gauge theory this is equivalent to rescaling H and M as Following [24], we note that the rescaling above leaves the FG form of the metric (2.12) invariant and transforms all the expansion coefficients homogeneously, except for g (4)ij , which acquires an inhomogeneous piece due to the logarithmic term in (2.13): g (4)ij = λ 4 g (4)ij + 2λ 4 log λ h (4)ij . (2.44) Similarly, the coefficients φ (0) and ψ (2) in the expansion of the scalar field transform homogeneously, whereas φ (2) acquires an inhomogeneous piece: It follows that the stress tensor and the scalar expectation value transform as where we have made use of (2.25) and (2.29). This immediately implies that these expectation values must take the form where µ 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 (2.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/µ) by log(H/M ) + log(M/µ), thus redefining Note also that there is no loss of generality in assuming that the scale µ 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 µ. 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 µ. The freedom in the choice of this scale is part of a bigger freedom in the choice of renormalisation scheme, as we discussed around (2.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 α and β.
In the flat-space limit, namely if H = 0, both the anomaly (2.21) and its contributions (2.33) and (2.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 (2.23) is the βφ 4 term. This produces a contribution to the stress tensor (2.32) that shifts its value by a term proportional to the boundary metric g (0)ij = η 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 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 In terms of these variables, the trace Ward identity (2.41) takes the form 3 Dynamics in flat space

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 [22]. This is useful because later we will use thermal equilibrium states on flat space to initialise the time evolution of non-equilibrium states on dS 4 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 (2.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 field 1 In this gauge the boundary is located at φ = 0 and the value of the scalar field at the horizon φ h is determined by the condition H(φ h ) = 0. After introducing a master field The function H(φ) appearing in this section should not be confused with the Hubble rate H appearing throughout the whole paper. 2 As in [22], we normalise the scalar field differently than in [26], which is the reason why some of the coefficients in (3.3) differ from those in the corresponding master equation in [26].
Close to the horizon a solution to the master equation can be expressed as a power series: Close to the boundary, φ → 0, the master field can be expanded as In practice we obtain a one-parameter family of solutions for G(φ), parametrised by the value of φ h , by numerically integrating (3.3) from a value of φ close to the horizon to a value close to the boundary using boundary conditions for G and G constructed from (3.4). The metric functions in (3.1) can then be obtained through the relations The temperature and the entropy density of field theory states dual to these numerically constructed geometries can be expressed in terms of (3.6) and (3.7) evaluated at the horizon (see [26,22] for details): In Fig. 1 (left) we show the reduced entropy density S ≡ 2π 2 s/N 2 divided by T 3 as a function of T /M for φ 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) with L IR given by (2.9). As explored in detail in [22], for real values of φ M the model has a smooth crossover between the IR and UV fixed point. 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 where p 0 is the pressure of the vacuum state in the T → 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 (3.9), 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 β = 1/4φ 2 M 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 constitute the equation of state (EoS) p eq ( ). We use the deviation of the ratio w ≡ p/ = P/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 2 s ≡ dp eq /d from its conformal value c 2 s,CFT ≡ 1/(d − 1) = 1/3. One advantage of w over c 2 s 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, E ≈ 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, η/s = 1/4π, is universal in all holographic theories with an Einstein gravity dual [27]. This means that knowledge of the entropy density (3.11) 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 [28]: (3.12) 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 ζ/ measures the viscous contribution to w in dS 4 . 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, E ≈ M 4 , ζ is non-zero.

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. [29] 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 where u i is the fluid velocity and is the projection of the covariant derivative to the spatial components in the local rest frame of the fluid. The EoS p eq ( ) and the transport coefficients η( ) and ζ( ) 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 are defined as the timelike eigenvector and the eigenvalue, respectively, of the stress tensor, namely T ij u j = − u i .
The leading term in the expansion (3.13) 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(∇ 2 ) and only consider the leading ideal and sub-leading viscous part of (3.13). We will write this term as where the viscous tensor Π ij in general contains contributions due to bulk and shear viscosity. 3 In this work we will only consider homogeneous and isotropic flows without shear stresses, in which case the viscous tensor simplifies to The possible scheme dependence of the microscopic energy-momentum tensor also manifests itself in the hydrodynamic approximation. As discussed in Sec. 3.1, p eq given in (3.11) 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 These are scheme-independent, and we may then view the EoS as a relation of the form ∆p eq = ∆p eq (∆ ). Next we rewrite the second term in (3.15) as This separates (3.15) into a scheme-dependent vacuum contribution T ij vac and a scheme- Note that the fact that T ij vac is proportional to the background metric g ij (0) is consistent with the expectation that the vacuum must respect the symmetries of this background. Since ∆T ij and ∆ are scheme-independent, we can now define a scheme-independent velocity field defined through the relation The scheme-independent part of (3.15) is then given by 4 Dynamics in de Sitter space

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 (2.3) of the dual gravity theory, following the method reviewed in e.g. [31][32][33], 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 dS 4 . However, it is useful to set up the problem with a slightly more general boundary metric of Friedmann-Lemaître-Robertson-Walker type For S 0 (t) = e Ht the boundary metric (4.1) becomes the dS 4 (2.14) with curvature scalar R = 12H 2 .
We use generalised Eddington-Finkelstein (EF) coordinates to parametrise the bulk geometry and the scalar field where the asymptotic boundary is located at r = ∞. The line element (4.2) has a residual gauge freedom in the radial coordinate which we exploit in our numerical scheme to fix the coordinate value r AH at the apparent horizon, defined by the conditionṠ(r AH , t) = 0, to a constant. Using (4.2) the equations of motion (2.3) result in the following set of equations where a prime denotes a radial derivative, f ≡ ∂ r f , and an overdot is short-hand for the modified derivativeḟ ≡ ∂ t f + 1 2 A∂ r f . Imposing the metric (4.1) and the asymptotic behaviour of the scalar field as boundary conditions, solutions of (4.4) can be expressed near the boundary as generalised power series: Note that the fall-off coefficientφ (2) in EF coordinates is generically different from the one in FG coordinates φ (2) . The coefficients a (4) (t) andφ (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 This relation follows from the momentum constraint (4.4e) and implies covariant conservation of the holographic stress tensor in the boundary theory, namely the first Ward identity in (2.41). The EF coordinate system (4.2) 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.32) and the scalar operator (2.35) assume the FG coordinate system (2.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 F G , related to ρ in (2.12) through ρ = 1/r 2 F G : where dots stand for therms with higher powers of log(r F G ). All these logarithmic terms appear because we are working with a curved boundary metric. The metric transforms as follows Substituting the general expressions for metrics in EF and FG coordinates into the transformation law (4.9) leads to a set of two equations which can be solved order by order in r F G and r F G log(r F G ) for the expansion coefficients in (4.8). For the leading behaviour of the metric in FG coordinates we find the following expressions where we have set r ≡ r F G and t ≡ t F G to shorten notation. After replacing the radial coordinate by ρ ≡ 1/r 2 we obtain explicit expressions for the non-vanishing metric components in FG coordinates (2.12) in terms of the coefficients in EF gauge The corresponding near-boundary expansion of the scalar field reads From the above expression one can now read off the relations between the expansion coefficients in FG and EF coordinates (4.20) 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: .
Similarly, the expectation value of the operatorÔ is given by As mentioned above, the coefficients α and β in (4.21) and (4.22) encode the schemedependence of E, P and O.

Numerical procedure
Our main interest is to compute the time evolution of P, E and O in a dynamical background geometry. For this purpose we have to solve the set of equations (4.4) 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 (4.4) 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,Ṡ, φ,φ 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 (4.4) but rather a set of equations for a new set of functions {S,S,φ,φ,Ã}. 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 (4.17). The new functions are therefore manifestly finite and take the form 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 ≡ 1/r where the boundary is located a z = 0. We now implement the following procedure to solve the initial value problem: 1. For a given radial profile of the scalar fieldφ(z, t 0 ) at some initial time t 0 we first solve the second-order Hamiltonian constraint equation (4.4a) forS(z, t 0 ). In principle this differential equation requires boundary conditions, but in our (subtracted) formulation this equation (and the following except the one forS) 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. Next we useφ(z, t 0 ) andS(z, t 0 ) in (4.4b) and solve forS(z, t 0 ). The boundary condition for this function readṡ Evaluating this requires knowledge of a 4 (t 0 ) andφ 2 (t 0 ), whereby a 4 (t 0 ) is required as a separate initial condition andφ 2 (t 0 ) needs to be extracted from the initial data for the scalar field.
6. Finally, we compute a 4 (t 0 ) from (4.7), whereby we obtain φ 2 (t 0 ) from the nearboundary expansion of ∂ tφ (z, t 0 ). We subsequently evolve the initial data to the next time slice t 1 = t 0 + ∆t for example via 26) 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 ξ(t). We use this freedom to fix the initial apparent horizon (AH) at z AH = 1/2. This is done by computingṠ for the gauge ξ 0 (t 0 ) = 0 first, and then numerically solving the apparent horizon equationṠ(z hor , t 0 ) = 0 (see also (4.45)). We then update where z hor = 1/r hor is the current location of the apparent horizon. To increase precision this procedure is repeated a few times, each time using an updated ξ(t 0 ). We solve the resulting equations numerically with a Chebyshev pseudo-spectral method (see e.g. [34]) 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 [35] with a time step ∆t = 1/(10N 2 ). After each time step we evaluate the momentum constraint (4.4e) to monitor the accuracy of the numerical evolution.

Initial states and time evolution
We are interested in studying the time evolution of gauge theory states in dS 4 . 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 dH/dt = 0 and the boundary geometry interpolates between flat space and dS 4 . 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 (4.2), 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φ(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. An example of this procedure is shown in Fig. 2, where all dimensionful quantities are given in units of M and we have chosen β = 1/16. In the left plot we show the initial guessφ(z, t 0 ) = 0 (red) together with the thermalised equilibrium resultφ(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 E, P and O evolve towards their equilibrium values. Energy conservation in flat space and homogeneity of the state imply that the energy density E remains constant during the evolution. The evolutions of P and O are not independent but are constrained by the Ward identity (2.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 ω 1 ≈ 6.00 + 4.64i.
We now turn on the Hubble rate at the boundary so that the boundary metric changes smoothly from Minkowski to dS 4 . We implement this by imposing the following relation on the function S 0 (t) in the boundary metric (4.1) This relation mimics a time-dependent Hubble rate that changes from zero at t t * to H at t t * in a time of order 1/Ω. We will refer to this transient period as "the quench". The corresponding form of S 0 (t) follows from integrating (4.28) subject to the initial condition S 0 (t = 0) = 1 (4.29) In Fig. 3 we illustrate this procedure for several choices of parameters. We will refer to each of these choices as a "protocol". 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 α = 0, β = 1/16. On the left plot we show the ratio S 0 (t)/S 0 (t) and on the right the evolution of the energy density. We vary both the "quench parameter" Ω, 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 t * 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 Ω = 4, t * = 1 and a 4 (t = 0) = −100. The only exception will be when H ≥ 2, in which case we will use Ω = 8 in order to have a quench time Ω −1 that is sufficiently shorter than the expansion rate H −1 .
After a time t − t * Ω −1 the boundary geometry settles down to dS 4 and we can analyse the evolution of field theory states on this expanding background. In Fig. 4 we show the evolution of E, P, O and E + P (which in equilibrium would be enthalpy) for a number of different values for the expansion rate H. 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- is encoded in the dependence on α and β of Eqs. (4.21). As mentioned above, in order to fix the scheme we choose β = 1/4φ 2 M = 1/16 throughout the paper. Moreover, in this section (but not elsewhere) we set α = 0. Note that the enthalpy E + P is scheme-independent because the dependence on α and β cancels out in this combination. In Fig. 4 (bottom right) we see that E + 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 α = 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 α, β 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 In the latter case we require the energy and pressure of the vacuum state to vanish. In the case of (4.31) 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.

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 [3,36,4] and non-conformal [37,22] 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 M or H 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 ∼ M , which in our units means H ∼ 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 α. To do this we start with the general form of the energy-momentum tensor for homogeneous and isotropic states where u i = 1, 0 is a future pointing time-like vector and E(t) and 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 dS 4 , which in addition to (4.32) also satisfy the relation E(t) = −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 dS 4 such that the energy-momentum tensor of these late-time states can be written as where P ∞ is scheme-dependent because it depends on α. In the spirit of Sec. 3.2 we build scheme-independent combinations of the energy-momentum tensor components in which all dependencies on α cancel out: 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 (4.34) in (3.21) we obtain the hydrodynamic approximation for the pressure in the Landau frame 4 where we have used ∇ 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. 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. 5 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 (4.35) and the blue dotted lines show the ideal hydrodynamic approximation according to the equation of state ∆P eq (∆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 * ) ∼ 1.4. It is remarkable that around these times the difference between 4 Both bulk viscosities in (4.35) and (3.13) 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 α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. 5 More precisely, what is shown is actually 3P/E for the value of α specified in (4.31) which, after the quench, coincides with 3∆P/∆E. 3 Δ/Δℰ Figure 5. Scheme-independent ratio of the excess pressure over the excess energy density defined in (4.34) as a function of time for the same initial state but different values of H. 5 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 ≥ 1 viscous hydrodynamics never provides a good approximation to the exact result.
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 [3][4][5], 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∆t ≈ 0.76, after which the ∆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. 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 ∆E → 0 the full evolution always deviates from the hydrodynamic evolution and can characteristically develop a negative pressure excess for H 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.

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 latetime, 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 [17] 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 ∇ i ∆T ij = 0, which for the dS 4 metric (2.14) evaluates to At late times ∆E(t) is well described by where A is the amplitude of the fluctuation and ω is a purely imaginary quasinormal mode with Imω < 0. By (4.36), this implies that ∆P(t) and ∆E(t) satisfy at sufficiently late times an EoS-like relation of the form This effective EoS explains the constant late-time ratio 3∆P(t)/∆E(t) in Fig. 5. We can extract an estimate for w eff from our numerical results at late times. In Fig. 7  It may therefore seem surprising that, in general, w eff = 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 [38], 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 M or H 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 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 [17], which suggests that conformal symmetry is restored at late times in the limit of large expansion rates H M . We indeed find evidence for this behaviour in Fig. 5: w eff obtains its minimal value w min eff = −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 3 because the numerics becomes increasingly challenging, we expect w 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 quasinormal 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−tq) 2 with amplitudes A = 0.05, 0.01 and 0.005 at quenching times t q = 3.5, 5 and 7, respectively. The quenches shown perform work on the system, which leads to an increase of energy and, as a consequence, also to an increase of E + 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 eff , i.e. by the first QNM. In the next section we discuss the maximally symmetric late-time state in more detail. 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 E and explains the increase of E + 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).

Late-time solution
At late times the geometry and the scalar field approach the following form (see also [18]) where A ∞ , S ∞ and φ ∞ are the time independent late-time limits of the metric functions and the scalar field defined as  In the following we will drop the index "∞" and implicitly assume this limit in all appearances of A, S and φ. Under these conditions the equations of motion simplify to The late time geometry (4.40) possesses an event horizon located at r = r EH defined by the condition A(r EH ) = 0 . Using this condition in (4.42d) one finds that the surface gravity at the event horizon equals the Hubble rate 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 AH is determined byṠ| r AH = 0, which in the coordinate system (4.41) gives the condition 0 = S(r AH ) d dt e Ht + 1 2 e Ht A(r AH )S (r AH ) . (4.45) Using this condition in (4.42d) together with (4.42e) allows us to express surface gravity at the apparent horizon We arrive at the conclusion that surface gravity at the apparent horizon equals minus surface gravity at the event horizon Above, we used the ansatz (4.40) 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 (4.40) very accurately. However, it is possible to arrive at (4.44) 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 [39]. By the holographic duality the bulk metric dual to this vacuum state has to satisfy the isometries of dS 4 as well. A domain wall parametrization in FG coordinates makes these isometries manifest 6 where the time and radial coordinate are related to those in (4.40) by Demanding that (4.40) can be transformed to the manifestly dS 4 -isometric domain wall form (4.48) gives the relation Beyond the event horizon (r < r EH ) the metric function A ∞ (r) is negative and naively one could expect the right-hand side of (4.50) to acquire an imaginary part. The exponent, however, diverges as r approaches the event horizon, which leads to a term a log(r − This gives precisely the same result for κ EH = lim r→r EH 1 2 A (r) as (4.44), only by demanding dS 4 symmetry and regularity of the solution at the horizon. It is possible to derive (4.46) 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. 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 κ EH = −κ AH = H. For the gravity dual of N = 4 SYM theory on de Sitter space apparent and event horizon densities can be straightforwardly computed because the bulk geometry is known explicitly [40,16]. This geometry has no apparent horizon but an event horizon whose area density is given by For comparison, a thermal state in N = 4 SYM on flat space has S(r EH ) 3 = κ 3 EH /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 → 0), the apparent horizon area density (S(r 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 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 hor , t) 3 is a monotonously growing function of time in accordance with Hawinkgs area theorem [41] and the more recent discussion 7 [18,20] 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.

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 horizon 8 .
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 [19] for related difficulties in interpreting the apparent horizon area).
A robust and gauge independent notion of entropy is given by entanglement entropy [42] of spatial subregions R in the QFT defined as S R = −Tr RρR logρ R , (5.1) 7 We are grateful to Alex Buchel for bringing this to our attention. 8 We 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.
whereρ R = TrRρ denotes the reduced density matrix obtained by performing on the full density matrixρ a partial trace over the degrees of freedom outside R. For simplicity we will assume these subregions to be spatial balls at some fixed time t = t 0 with radius 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ρ can be timedependent, 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 → ∞ the entangling region R covers an entire spacelike slice of Minkowski space. In this limit ρ R =ρ and (5.1) equals the von Neumann entropy of the full density matrix, i.e., the thermodynamic entropy of a quantum state in thermal equilibrium. For H = 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 > H −1 are causally disconnected, quantum states on such regions can be entangled [43]. Direct field theory computations of entanglement entropy are only possible in exceptional cases like for example in two-dimensional CFTs [44] and in dimensions higher than two only in non-interacting QFTs [45]. The holographic duality replaces the field theory computation of entanglement entropy by a much simpler extremisation problem for the surface area A R of a codimension two surface, homologous to R, in the bulk theory [46,47] The boundary of the relevant surface coincides with the boundary of the entangling region R in the field theory and extremises the area functional in the bulk theory The surface embedding X µ = X µ (σ a ) is parametrised by three intrinsic coordinates for which we choose σ a = {r, θ, ϕ}. The entangling regions (5.2) do not break spherical symmetry in the boundary theory. We can then parametrise the bulk surface with X µ (z) = {z(r), t(r), r, θ, ϕ} . where the metricḡ αβ is related by a conformal factor to a three dimensional subspace (α, β = {t, z, r}) of the bulk metric (4.40) (see also [48]) ds 2 =ḡ αβ dx α dx β = (r S(z) a(t)) 4 g αβ dx α dx β .
where we use a(t) ≡ e Ht in the following to simplify notation. The equations of motion that follow from δA R = 0 take the form of a non-affine geodesic equation where Γ α βγ is the Levi-Civitá connection associated toḡ αβ and is meant to be evaluated at the location of the surface X α (r); the viscous friction term on the right hand side includes the Jacobian J = d 2 τ (r) dr 2 / d 2 τ (r) dr 2 that originates from transforming from the affine parameter τ defined by dX α (τ ) dτ dX β (τ ) dτḡ αβ = 1 to the non-affine parameter r. For φ = 0 the field theory is N = 4 SYM theory and (5.8) has a simple analytic solution z(r) = √ 2 − r 2 and t(r) = − log (1 + z(r)) 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 [18] and the holographic entanglement entropy is equal to the area of a hemisphere in AdS 5 [46].
The sign of z (r = 0) determines if the bulk surface can reach the boundary (z (r = 0) < 0) or not (z (r = 0) > 0). The condition z (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 (z) − S(z)a (t) = z 2 A(z)S (z) − HS(z) . (5.10) Plugging this relation into the definition ofṠ giveṡ S + 1 2 z(r) 2 A(z(r))a(t(r))S (z(r)) = 0 . (5.11) Except for the second term, this relation is equal to equation (4.45) that determines the location of the apparent horizon. Since z(r) 2 A(z(r))S (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 [49] for similar results on extremal surface barriers). This can be shown by combining (5.10) with the Einstein equations (4.42b) to (4.42e) which gives Furthermore, when the geometry is static (H = 0), we know from (4.45) that the apparent horizon (Ṡ = 0) and event horizon (A = 0) coincide. From (5.11) we then see that this also solves the entanglement horizon equation, so that all three horizons coincide as expected [50]. We compute the entanglement entropy numerically by shooting, i.e., by integrating (5.8) from the turning point z * of the geodesics to a cutoff-value for the holographic coordinate z cut = /(1 − ξ ) close to the boundary located at z = 0. Note that the gauge freedom ξ in S(r) = r + ξ + O(1/r) near the boundary leads to an additional 1/ dependent divergence in the entanglement entropy. We eliminate the 1/ divergence with the gauge transformation r → r − ξ that ensures S(r) = r + O(1/r) and in addition fixes the residual gauge freedom (r → r + ξ) in the radial coordinate. The cutoff regulated value for the entanglement entropy can then be obtained by numerically solving the following integral In Fig. 10 we show extremal surfaces in the gauge where the apparent horizon is fixed at z AH = 1/2, obtained by shooting from different values of z * , for several spherical entangling regions of different radius . Interestingly, surfaces with z * = z 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 * = 0.3597 (0.1961) for H = 1 (2.5). Extremal surfaces of entangling regions larger than the cosmological horizon ( > 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 AH = 1/2, the entanglement horizon is located at z = 0.4554 (0.3434) and z EH = 0.3597 (0.1961) for H = 1 (2.5). For comparison, in the gravity dual of N = 4 SYM theory on de Sitter the event horizon is located at z EH = 1/(2H − ξ) and extremal surfaces can probe until z = 1/(H − ξ) in the gauge where z AH = −1/ξ.
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 'superobserver' 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 . Asymptotically the area equals [43] A = 4π 2 1 2 2 + 1 3 log( ) + O(1) + O(log( )). (5.14) 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 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 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  Figure 11. Blue and green curves are two families of null geodesics, originating on a generic point of the extremal surface with z * = z 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). the horizon is for H = 1 approximately given by t Ent ≈ −0.3 − log( ). In combination with a(t) = e Ht and S 0 (z Ent ) ≈ 0.69 this implies that the part of the surface 9 hovering the entanglement horizon has an area 4 3 π 3 0.134/ 3 ≈ 0.559 which is only a constant term as opposed to the usual volume scaling of the time independent setting.

Discussion
We have used a holographic model to study the dynamics of a strongly coupled nonconformal 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 P = wE. 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 [51] also applies to cosmological horizons [52]. Therefore, in analogy with the Bekenstein-Hawking temperature of black holes [53,51], an observer living at the boundary would associate a temperature to the cosmological horizon of de Sitter space where the surface gravity κ 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π, and another temperature at the deeper apparent horizon, equal to −H/2π (also found in [19]). In the literature there are several works [54][55][56] 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 where A dS = 4π/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 (6.2) 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 [43]). 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 [57].
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 ∼ 1/H (see Fig. 7 and 8), in agreement with [17]. In our case this time can be a parametrically different from ∆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 ∆E will quickly become smaller than T 4 , with T ∼ H the background de Sitter temperature (note that this temperature is the minimal temperature, accelerating observers will see an even higher temperature [58]). 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.