Entanglement and out-of-equilibrium dynamics in holographic models of de Sitter QFTs

In this paper we study various aspects of entanglement entropy in strongly-coupled de Sitter quantum field theories in various dimensions. We find gravity solutions that are dual to field theories in a fixed de Sitter background, both in equilibrium and out-of-equilibrium configurations. The latter corresponds to the Vaidya generalization of the AdS black hole solutions with hyperbolic topology. We compute analytically the entanglement entropy of spherical regions and show that there is a transition when the sphere is as big as the horizon. We also explore thermalization in time-dependent situations in which the system evolves from a non-equilibrium state to the Bunch-Davies state. We find that the saturation time is equal to the light-crossing time of the sphere. This behavior is faster than random walk and suggests the existence of free light-like degrees of freedom.


Introduction
Quantum field theory in curved spacetimes [1,2] is a subject of great relevance that has lead to many interesting areas of research in the past few decades. Although it is believed to be a good physical description in circumstances where quantum gravitational effects do not play a dominant role, the reader must also be aware of its limitations. For instance, current research suggests that effective QFT might break down in places where the curvature of spacetime is not so large, leading to various well-known puzzles and paradoxes [3][4][5]. Albeit gravity is treated classically, QFT in a fixed background has provided us historically with a tool to investigate qualitatively some of the most fundamental questions in quantum gravity, e.g. understanding various aspects of black hole thermodynamics and deriving the physical consequences of inflation and modern cosmology.
QFT in curved spacetimes is often discussed in terms of free field theory but it has been extended to perturbative studies of weakly coupled field theories. While most of the qualitative features are visible at this level, understanding the strong-coupling and nonperturbative dynamics is of interest. Indeed, in the context of cosmology it is plausible that strong dynamics might have played an important role in the early universe.
In recent years, the discovery of the AdS/CFT correspondence [6][7][8] has provided tools for the study of a large class of strongly-coupled QFTs. To date, there have been a number of proposals for holographic theories living in de Sitter space and other cosmological backgrounds [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. The purpose of the present paper is to study entanglement generated during the cosmological evolution [29][30][31][32][33][34][35][36][37][38]. Entanglement entropy (or geometric entropy) is an important concept and a useful tool in quantum field theories and quantum many body systems, and serves as a probe to characterize states of matter with long range correlations. As a first step we will consider strongly coupled gauge theories in de Sitter spacetime in the large-N limit and we will analytically compute entanglement entropies of spherical regions in the conformally flat patch and the static patch using the holographic prescription. We will then use entanglement entropy as a tool to study thermalization of out of equilibrium configurations in de Sitter space.
It is well known that entanglement entropy of a spatial region in a local field theory is UV-divergent The remarkable feature of the entanglement entropy in de Sitter is that S div is fully regulated by subtracting the flat space result. 1 On the other hand, S f inite contains information about the long range entanglement and it is expected to be more sensitive to the curvature of space-time. Another local way to deal with the UV divergences was introduced in [39] which amounts to compare entanglement entropies of spheres of similar radii.
Let us now focus on the finite part S f inite . When the size of the sphere is much smaller than the de Sitter horizon, R << H −1 , S f inite is expected to be the same as in flat space-time. However, for a particular conformal field theory 2 (CFT) which has a gravity dual we will show that S f inite in (1 + 1) and (2 + 1) dimensions is exactly the same as in flat space-time for R < H −1 .
In (3 + 1)-dimensions, S f inite is more sensitive to the curvature and we find for where c is a constant and A is the proper area of the sphere. We will also show that the entanglement entropy of a sphere of radius R in the static patch of de Sitter is the same as the entanglement entropy of a sphere of proper radius R in the conformally flat patch of de Sitter when the size of the sphere is smaller than the size of the event horizon. Behaviors of the extremal surfaces for radii smaller and larger than the horizon are entirely different and as a consequence the entanglement entropy undergoes a phase transition at R = H −1 . This phase transition is a signal of a drastic change in correlations at distance R = H −1 . It is important to note that only a super-observer can "see" this phase transition of the entanglement entropy.
The above feature is more prominent in the so-called renormalized entanglement entropy, introduced in [39], which is a derived quantity that has some advantages over entanglement entropy. Entanglement entropy of a region of size R is sensitive to the physics from scale R all the way down to the short distance cut-off , whereas renormalized entanglement entropy S d (R) is expected to be most sensitive to degrees of freedom at scale R and it somewhat naturally describes the RG flow of the entanglement entropy. In particular, for the vacuum of Lorentz invariant, unitary QFTs S d (R) in (1+1) and (2+1) dimensions is monotonically decreasing and non-negative, providing a central function out of the entanglement entropy. In de Sitter we will show that S d (R) = constant when R < H −1 . In (1 + 1) and (2 + 1) dimensions S d (R) = 0 when R > H −1 because regions separated by a distance R > 1/H are causally disconnected. In (3 + 1)-dimensions, the renormalized entanglement entropy is more complicated and it is neither monotonic nor positive for R > H −1 . It is perhaps an indication that the definition of S 4 (R) should be modified in order to construct a central function out of the entanglement entropy.
Entanglement is also a useful probe in out-of-equilibrium configurations. For instance it has been shown that, in comparison to other non-local observables (two-point functions and Wilson loops), entanglement entropy equilibrates the latest when the system undergoes a global quench [40][41][42], thus setting the relevant time-scale for the approach to thermal equilibrium. We will explore the issue of thermalization in dS space and we will use entanglement entropy to characterize the time evolution of the system. The time-dependent configurations relevant in the present context are achieved by turning on for a short time interval a uniform density of sources. The relevant interval of time δt is taken to be much smaller than any other scale in the system, and then is turned off. The work done by the sources take the system to an excited state which subsequently equilibrates under the evolution of the same Hamiltonian before the quench.
In the context of QFT, the study of such non-equilibrium situations is a serious challenge and a topic of current research. In a recent paper [43], Calabrese and Cardy showed that for a variety of (1 + 1)-dimensional CFTs as well as for some lattice models, entanglement entropy for a segment of length grows linearly in time as and then reaches saturation at some t = t sat . Here, ∆S represents the difference of the entanglement entropy with respect to the initial state, and S sat is the equilibrium value after it reaches saturation. They argued that this remarkably simple behavior can be understood from a simple model of entanglement propagation using free-streaming quasiparticles traveling at the speed of light. This kind of non-equilibrium configurations has also been considered holographically by studying the gravitational collapse of a shell of null dust in AdS, i.e. by means of the so-called Vaidya geometries, starting with the seminal works [44][45][46] and continuing with a large body of work that includes the recent additions . Of particular relevance for the present context are the results of [69,74]. In these papers Liu and Suh showed that, for holographic theories with a gravity dual, entanglement entropy also undergoes a series of regimes resembling those in phase transitions: pre-local-equilibration quadratic growth, post-local-equilibration linear growth, memory loss, and saturation. These results apply for strongly-coupled QFTs in Minkowski space.
A natural question that arises here is that if these universal features also show up in theories living in different backgrounds, and in de Sitter space in particular. The study of out-of-equilibrium physics in the theories we are considering in the present paper is also interesting on its own right. On one hand, it is known that the early universe went through a series of epochs with varying temperature and intricate stronglycoupled dynamics. While there are known techniques in field theory that allows us to handle certain non-equilibrium problems [75][76][77], the fact that the theories under consideration are strongly-coupled render the perturbative methods unreliable. Here is where holography plays a useful role.
This rest of the paper is organized as follows: in section 2, we study gravity solutions that are dual to field theories in a fixed de Sitter background. We show explicitly two possible slicings, one that is dual to the static patch and a second one that is dual to the conformally flat patch. In section 3 we give a brief review of entanglement entropy, the Ryu-Takayanagi prescription and two useful derived quantities that we shall study in the present context: mutual information and renormalized entanglement entropy. In section 4 and 5 we compute the entanglement entropy in the Bunch-Davies vacuum of de Sitter QFTs in various dimensions. Section 6 is devoted to study the evolution of entanglement entropy in time-dependent configurations. Finally, in section 7 we make some comments about our results and close with conclusions.

Gravity solutions with dS slices
The purpose of this section is to study gravity solutions to Einstein's equations with a foliation such that the boundary metric corresponds to de Sitter space in a given coordinate system. In (d + 1)-dimensions, the Einstein-Hilbert action with negative cosmological constant is given by which gives the following equations of motion Here Any asymptotically AdS metric can be written in the Fefferman-Graham form [78] from which the dual CFT metric ds 2 = g µν (x)dx µ dx ν can be directly read off as g µν (x) = g µν (0, x). The full function g µν (z, x) also encodes data dual to the expectation value of the CFT stress-energy tensor T µν (x). More specifically, in terms of the near-boundary expansion the standard GKPW recipe for correlation functions [7,8] leads after appropriate holographic renormalization to [79][80][81] where X (d) µν = 0 ∀ odd d (reflecting the fact that for odd boundary dimensions, there are no gravitational conformal anomalies), and X (2d) µν for d ≥ 3 given by similar but longer expressions that we will not transcribe here. In (2.6) it is understood that the indices of the tensors g µν (x) are raised with the inverse boundary metric g µν (x).

The static patch
The de Sitter spacetime in d-dimensions has an isometry group of SO(d, 1). For free field theory, there is a family of de-Sitter invariant vacuum states and it is known as the α-vacua [82,83]. However, among these states only the Bunch-Davies (or Euclidean) vacuum [84] reduces to the standard Minkowski vacuum in the limit H → 0.
This state is well defined on the entire manifold but, for concreteness, in this section we will center our discussion in the static patch of de Sitter, which covers the causal diamond associated with a single geodesic observer, The name "static" comes from the fact that it has a killing vector ∂ t associated with the isometry of time translations. Therefore, energy as well as entropy are well defined quantities. For such an observer, the Bunch-Davies vacuum is characterized by a temperature T dS = H/2π that is associated to the presence of a cosmological horizon at r = 1/H, where H denotes the value of the Hubble constant [85]. This can be obtained by continuing to Euclidian space and imposing regularity on the horizon.
To obtain a solution dual to a QFT in a given background we start by writing the d + 1-dimensional metric of the bulk in the Fefferman-Graham form. In particular, for static dS we can assume that the bulk is independent of time, The next step is to write the functions f (r, z), j(r, z) and h(r, z) as a series expansion in z and solve order by order the Einstein equations (2.2).
f (r, z) = f 0 (r) + f 2 (r)z 2 + f 4 (r)z 4 + · · · , j(r, z) = j 0 (r) + j 2 (r)z 2 + j 4 (r)z 4 + · · · , (2.9) The first coefficients in the above expansions are related to the boundary metric, so if we want the boundary theory to be defined on the static patch of de Sitter we have to impose that f 0 (r) = 1 − H 2 r 2 , j 0 (r) = 1 1 − H 2 r 2 , and h 0 (r) = r 2 . (2.10) After plugging the ansatz (2.10) in Einstein equations (2.2), we find that for all d the solutions are truncated at order O(z 4 ), For these metrics, the curvature scalars are given by which show that these backgrounds are completely smooth and absent of singularities. The solutions have a regular Killing horizon at z = 2/H with constant surface gravity, and this is in fact related to the temperature of the field theory T dS = H/2π. In addition, there is also the expected horizon at r = 1/H ∀ z.
It is interesting to note that the authors of [24] derived a family of solutions dual to QFTs in de Sitter space for arbitrary temperature (not necessarily on the de Sitter invariant vacuum). These solutions are found to be related by a bulk diffeomorphism (that acts as a boundary conformal transformation) to the so-called hyperbolic (or topological) black holes described in [86][87][88]. We will come back to these solutions in section 6.
The energy-momentum tensor of the boundary theory can be obtained from (2.5), and turns out to be T µν (x) = 0 ∀ odd d, This result captures the correct conformal anomaly present in even dimensions.

The conformally flat patch
To obtain a gravity solution dual to a QFT living in another dS chart we can proceed as we did in the previous section. However, note that for our previous background, the and g dS µν (x) being the d-dimensional de Sitter metric. This fact provide us with a convenient shortcut: the bulk metric for any other dS chart can be obtained just by a coordinate transformation of (2.8). We can thus express g dS µν (x) in any coordinate system and it will immediately related to our previous solution via a trivial diffeomorphism that does not mix the z-coordinate (and thus preserving the Fefferman-Graham form).
Here we are interested in the planar patch of de Sitter. The planar patch is conformally related to Minkowski space and it covers one half of global de Sitter, Here η is the so-called conformal time and is related to the usual FRW time by the For this coordinate system then, the corresponding bulk metric reads where the function f (z) is given in (2.17). The analysis of curvature invariants is the same as in equations (2.12)-(2.14). The boundary stress-energy tensor is related to the one in the static patch by a trivial diffeomorphism: and T µν (x) = 0 ∀ odd d.

Entanglement entropy, mutual information and number of degrees of freedom
When we consider an arbitrary quantum field theory with many degrees of freedom, we can ask about the entanglement of the system. To describe the system, we define a density matrix, ρ, which is a self-adjoint, positive semi-definite, trace class operator. Now, on a constant time Cauchy surface let us imagine dividing the system into two subsystems A and A c , where A c is the complement of A. The total Hilbert space then factorizes as H total = H A ⊗ H A c . For an observer who has access only to the subsystem A, the relevant quantity is the reduced density matrix defined as The entanglement between A and A c is measured by the entanglement entropy, which is defined as the von Neuman entropy using this reduced density matrix In AdS/CFT, Ryu and Takayanagi [89] conjectured a formula to compute the entangle- ment entropy of quantum systems with a static (d + 1) dimensional gravitational dual. The entanglement entropy of a region A of the quantum field theory is given by Figure 3. The two disjoint sub-systems A and B, each of length l along X-direction and separated by a distance x. The schematic diagram on the right shows the possible candidates for minimal area surfaces which is relevant for computing S A∪B . See [92] for a detailed discussion.
where G N is the bulk Newton's constant, and γ A is the (d − 1)-dimensional area surface such that ∂γ A = ∂A. For background with time dependence, this proposal has been successfully generalized to [90] where now minimal area surface is replaced by extremal surface. It is well known that entanglement entropy of a spatial region in a local field theory is UV-divergent Only local physics contributes to the UV-divergent piece S div . On the other hand S f inite contains information about the long range entanglement. Mutual information is a quantity that is derived from entanglement entropy. Mutual information between two disjoint sub-systems A and B is defined as where S A , S B and S A∪B denote entanglement entropy of the region A, B and A ∪ B respectively with the rest of the system (see fig. 3 for an example). Mutual information is a UV-finite quantity and hence it does not depend on regularization scheme. Moreover, as showed in [91], given an operator O A in the region A and O B in the region B, mutual information sets an upper bound and thus measures the total correlation between the two sub-systems: including both classical and quantum correlations. Renormalized entanglement entropy, introduced in [39], is another derived quantity which is UV-finite. For a region A with a single length scale R, it is defined as is the entanglement entropy of the region A. Renormalized entanglement entropy S d (R) has some nice properties: (i) for a CFT S d (R) is a constant, (ii) for a renormalizable quantum field theory both S d (R → 0) and S d (R → ∞) are constants and as R is increased from zero to infinity S d (R) interpolates between them, (iii) S d (R) is expected to be most sensitive to degrees of freedom at scale R. S d (R) somewhat naturally describes the RG flow of the entanglement entropy with distance scale [39,93]. Particularly, for d = 2, 3 and 4, S d (R) is given by

Entanglement entropy in (1+1) dimensions
In this section we will compute the entanglement entropy of an interval x ∈ − a 2 , a 2 at time η = η 0 for a strongly coupled QFT living in the conformally flat patch of de Sitter in (1 + 1) dimensions, using the holographic prescription (3.4). We choose the standard Bunch-Davies vacuum state of the field theory and the dual bulk metric is given by (2.20). The minimal area surface can be parametrized by two functions: x(z) and η(z) with the boundary conditions The area functional is given by, Before, we proceed, let us perform a change of variable that will make our life easier 3) x(z) = − Hη(z)r(z). Note that ∆r(0) = ∆x(0)/H|η 0 | is the proper length of the interval ∆x(0). In terms of these new functions, the action becomes with boundary conditions: where τ 0 is related to η 0 and l is the proper length of the interval a. The action is independent of η, however now we are restricted to lH < 2, i.e. the interval is smaller than the size of the horizon. Equations (4.2, 4.5) lead to two conserved quantities: where, c 1 , c 2 are constants and (4.9) Therefore, for lH < 2, we obtain dτ dx = e Hτ c 1 c 2 . (4.10)

Connected solution
When lH < 2, we can use action (4.5) to obtain U-shaped solutions. From the action (4.5) it is clear that for these solutions τ (z) = τ 0 and only r(z) has a nontrivial profile. Let us introduce rH = sin θ and in terms of θ(z), equation of motion is obtained to be and the area is given by where, z c is the closest approach point obtained from the boundary condition (4.13)  Therefore, the entanglement entropy is given by, where, c = 3L

Disconnected solution
The disconnected solution is given by : x(z) = ± a 2 , η(z) = η 0 and the corresponding entanglement entropy is given by Therefore, there is a transition from the connected solution to the disconnected solution at lH = 2. It is important to note that l < 2/H is the region accessible to a single observer and hence only a "super-observer" can "see" this transition of entanglement entropy.
Finally the entanglement entropy is given by, where, Θ is the Heaviside theta function. Few comments are in order: note that in (1+1) dimensions, entanglement entropy of an interval of length l in flat space is the same as entanglement entropy of an interval of proper length l in the conformally flat patch of de Sitter, when the interval is smaller than the size of the horizon for a single observer. However, a "super-observer" will see a difference in the behavior of the entanglement entropy.

Mutual information
Let us now calculate mutual information of two intervals A and B of length a (proper length l) separated by a distance b (proper length x). Mutual information undergoes an interesting entanglement/disentanglement "phase-transition" (see figure 5) which is qualitatively different from what has been discussed in [92]. We will consider three separate cases: In this case, mutual information is independent of H: Let us now comment on a few key features of mutual information in this case. When the intervals are comparable to the horizon size (i.e. lH ∼ O(1)), entanglement between them increases. However, mutual information between any two intervals separated by a proper distance xH ≥ 2 is identically zero. It implies that for any two intervals A and B, separated by a distance larger than the horizon size, ρ A∪B = ρ A ⊗ ρ B and hence they are completely disentangled.

Renormalized entanglement entropy
For d = 2, S 2 (R) is monotonically decreasing for all Lorentz-invariant, unitary QFTs and hence it indeed describes RG flow of the entanglement entropy. It is interesting to investigate how S 2 (R) behaves for de Sitter QFTs in the Bunch-Davies vacuum state. Particularly, when lH < 2, renormalized entanglement entropy is non-zero and independent of lH S 2 (l) = c 3 . (4.20) However, for lH > 2 it vanishes S 2 (l) = 0. (4.21) Thus, the renormalized entanglement entropy undergoes a phase transition. This is not very surprising since S 2 (l) measures entanglement correlations of a system at distance scale l and two points separated by a distance l > 2/H are causally disconnected.

Entanglement entropy of a sphere in d−dimensions
In this section, we will compute the entanglement entropy of a spherical region d−1 i=1 x 2 i ≤ a 2 for a field theory living in the conformally flat patch of d−dimensional de Sitter spacetime. We choose the Bunch-Davies vacuum state and the corresponding dual (d + 1)−dimensional bulk metric is given by (2.20). We will show that the entanglement entropy of a sphere of radius R in the static patch of de Sitter is the same as entanglement entropy of a sphere of proper radius R in the conformally flat patch of de Sitter when the size of the sphere is smaller than the size of the horizon.
Because of the spherical symmetry, the extremal surface can be parametrized by two functions: ρ(z) and η(z) and the area functional is given by where, Ω d−1 is the area of a unit sphere in d − 1 dimensions. The boundary conditions are where, is the short distance cut-off that we need to introduce in order to regularize the area of the extremal surface. The proper radius of the sphere is given by, R = a/(|η|H).

RH < 1: Connected solution
When ρ(z)/|η(z)| < 1, we can again perform a change of variables: 3 where, Hr(z) < 1. In terms of these new functions, the action becomes with boundary conditions: where τ 0 is related to η 0 by equation (5.3) and R is the proper radius. Before proceeding further, a few comments are in order: First, note that the action (5.5) is independent of τ and hence it is clear that τ (z) = τ 0 is a solution. Also note that the action (5.5) is identical to that for a spherical region in the static patch of de Sitter, where the bulk metric is given by (2.8). Therefore, entanglement entropy of a sphere of radius R in the static patch of de Sitter with the observer at the center of the sphere, is the same as entanglement entropy of a sphere of proper radius R in the conformally flat patch of de Sitter, provided the size of the sphere is smaller than the size of the horizon. However, it is not very clear how to generalize this result for non-spherical regions because boundary conditions are different in different coordinates. In the conformally flat patch, temporal boundary condition is given by η =constant, whereas, in the static patch it is given by τ =constant. These two boundary conditions are identical only for spherical regions. Now we will obtain the equation of motion for the radial profile using the action (5.5) It has a simple solution: where R 0 is a constant which can be computed using the boundary condition r( ) = R, yielding Now the area is given by, where, such that r(z c ) = 0.

RH > 1: Disconnected solution
Now we will do the following change of variable, This change of variable is valid only for Hr > 1. In terms of these new functions, the action becomes  with RH > 1 and τ 0 is related to η 0 and R is the proper radius. Note that the action is again independent of τ . Solutions of the equations of motion are τ (z) =τ 0 , (5.16) .
Because the area of the extremal surface is divergent, we will again introduce a short distance cut-off z = such that r( ) = R. Therefore, Behaviors of the extremal surfaces for RH < 1 and RH > 1 are entirely different and for RH > 1, extremal surfaces do not penetrate the region inside RH = 1; this feature is schematically shown in figure 6. 4 Also note that the behavior of the extremal surfaces near the Killing horizon is rather unique, since the near horizon region of the bulk contributes insignificantly to the area of the surfaces. 5 Presence of the Killing horizon is the bulk reflection of the boundary causal structure and it is responsible for the phase transition of the entanglement entropy.

Entanglement entropy in (2 + 1)-dimensions
For H = 0 and d = 3, we recover the flat space results: For non-zero H and RH < 1, we obtain where, c 3 = L 2 /4G 3+1 N . The fact that S H (R) = S 0 (R), for regions smaller than the size of horizon, is probably related to the absence of conformal anomaly in odd dimensions.
For RH > 1, entanglement entropy is an area-worth quantity There is a "phase transition" of the entanglement entropy at RH = 1 where the first derivative of the entanglement entropy is discontinuous:

Entanglement entropy in (3 + 1)-dimensions
For H = 0 in d = 4, we have For non-zero H and RH < 1, we obtain 4 A similar behavior of the extremal surfaces has also been observed for Kasner-AdS soliton background [35]. 5 For a discussion on behavior of the extremal surfaces near an event horizon see [94,95]. Hence, where, c 4 = L 3 /4G 4+1 N . Note that in d = 4, energy momentum tensor of the boundary theory T µν (x) = 0. It is also interesting to note that at finite temperature in flat spacetime S T (R) − S 0 (R) ∼ R 4 T 4 . Whereas in de Sitter at strong coupling S H (R) − S 0 (R) ∼ −R 2 T 2 dS and hence the behavior is completely different from what one would obtain for the same field theory in flat space-time at T = T dS .
For RH > 1, the entanglement entropy is given by Again there is a phase transition at RH = 1. This phase transition is a signal of a drastic change in correlations at distance R = 1/H. It is important to note that R < 1/H is the region accessible to a single observer and hence only a super-observer can "see" this phase transition of the entanglement entropy.
In the limit R >> 1/H, we obtain and hence the finite piece is proportional to the proper area of the sphere. Note that in the limit R >> 1/H, we do not have a term proportional to the number of e-foldings mainly because our extremal surfaces are disconnected in this limit. If one imposes smoothness of the solution as an additional criterion for it to be a good Ryu-Takayanagi surface, then in the limit R >> 1/H one will get an extra term which is proportional to the number of e-foldings [32].

Renormalized entanglement entropy
The entanglement entropy of a sphere is continuous but not smooth at R = 1/H. As a consequence, the renormalized entanglement entropy undergoes a phase transition at R = 1/H where it is discontinuous; however a single observer will never see this phase transition.

d = 3
For the vacuum of Lorentz invariant, unitary QFTs S 3 (R) is monotonically decreasing and non-negative [96], providing a central function for the F-theorem. It is interesting to investigate how S 3 (R) behaves for de Sitter QFTs in the Bunch-Davies vacuum state. When RH < 1, renormalized entanglement entropy is given by, For RH > 1, Thus, the renormalized entanglement entropy undergoes a phase transition. This is expected since S 3 (R) measures entanglement correlations of a system at distance scale R and two regions separated by a distance R > 1/H are causally disconnected.

d = 4
In d ≥ 4 dimensions, the renormalized entanglement entropy is more complicated, since there are indications that it is neither monotonic nor non-negative [39]. 6 However it is still expected that S 4 (R → 0) > S 4 (R → ∞) because of the a-theorem [97,98]. A similar behavior is observed for QFTs in de Sitter space-time. For RH < 1, renormalized entanglement entropy is a constant The renormalized entanglement entropy undergoes a phase transition at R = 1/H and for RH > 1, S 4 (R) is neither positive nor monotonically decreasing function of R However, S 4 (R → 0) > S 4 (R → ∞) still holds. One perhaps can argue that the behavior of S 4 (R) is still acceptable because non-monotonic, non-positive regions are not accessible to a single observer. It is also possible that non-positive, non-monotonic behavior of S 4 (R) is an indication that the definition of S 4 (R) should be modified in order to construct a central function out of the entanglement entropy.

Thermalization of dS QFTs
In this section we are interested in studying entanglement entropy in a thermalizing state of the same dS QFTs. We start by revisiting some generalities of the solutions found in [24] for holographic theories of dS with T = T dS and then we construct the time-dependent generalization of these geometries. We can think of these bulk solutions as the Vaidya version of the so-called hyperbolic black holes described in [86][87][88].

Hyperbolic black holes and dS QFTs with T = T dS
The basic idea goes as follows: the static patch of dS d is conformally related to the Lorentzian hyperbolic cylinder R × H d−1 where H d−1 is the Euclidean hyperboloid, (6.1) 6 It is known that S 4 (R) is non-monotonic and negative for several systems; for example GPPZ flow which describes the flow of N = 4 SYM to a confining theory under a mass deformation (which has UV-dimensions ∆ = 3) [39].
Moreover, holographic duals for theories living in the hyperbolic cylinder with Euclidean time period given by β = 1/T are know. These are the so-called hyperbolic (or topological) black holes, where dΣ 2 d−1 = dξ 2 + sinh 2 ξ dΩ 2 d−2 is the metric of the Euclidean hyperboloid. Alternatively, the mass can be given as where ρ + is the radius of the event horizon, given by the largest positive positive root of f (ρ). To find the Hawking temperature T for this solution, we must impose that the Euclidian time t E be periodically identified with appropriate period, t E ∼ t E + β. A simple calculation shows that This relation can be inverted to find which allows us to take T as the parameter that determines the solution. The case µ = 0 is isometric to AdS and is not properly a black hole as it is completely nonsingular. However, it covers a smaller portion of the entire manifold and the coordinate patch breaks down at ρ = ρ + = L. The Killing horizon in this case is analogous to a Rindler horizon, with associated inverse temperature β = 2πL, and non-vanishing area. Solutions with µ = 0 possess a true singularity at ρ = 0. In contrast with the regular AdS black holes, the zero temperature solution of (6.2) is different from the one that is isometric to AdS. In fact, there is a range of negative values for µ such that the solutions still possess regular horizons and have sensible thermodynamics. The minimum values of µ and ρ + allowed, for which the horizon degenerates, are For these values the corresponding black hole is extremal. In general, the Hawking temperature is monotonically increasing with respect to the black hole mass. In Figure  7 we can see this behaviour for black holes of various dimensions. Interestingly, the Penrose diagram for a hyperbolic black hole with negative µ is like that of a Reissner-Nordström-AdS black hole. For positive µ it is instead like that of a Schwarzschild-AdS black hole [99][100][101][102][103][104][105]. Now, the field theory dual to the hyperbolic black hole (6.2) lives in the hyperbolic cylinder which, according to (6.1), is conformal to the static patch of de Sitter. The fact that we are dealing with a CFT presents us with an interesting possibility: among the infinite number of conformal frames available to us, we can choose the one related to R×H d−1 via the specific Weyl transformation ds 2 → (1−H 2 r 2 )ds 2 . As first explained in [106], Weyl transformations in the boundary are dual to specific bulk diffeomorphism. 7 In our case, the desired transformation can be achieved by defining a new coordinate which brings the metric (6.2) into the form This bulk geometry is now dual to a CFT on the static patch of dS d at a temperature given by (6.5), which does not have to coincide with the de Sitter temperature T dS . For the particular case T = T dS , the background (6.9) differs from (2.8) only by a simple coordinate transformation. From the boundary theory perspective, we are taking the period of the Euclidean time to be 1/T instead of 2π/H. Of course, unless T = T dS , the Euclidean manifold has a conical singularity which map in the Lorentzian signature to the de Sitter horizon. Nevertheless, as long as we restrict our attention to study the physics inside the horizon, there is nothing wrong about the corresponding Lorentzian spacetime. Physically, we can think of this as having an external heat bath sitting at the de Sitter horizon which keeps the system at a temperature different from T dS . Therefore, the geometry does not smoothly continue through the horizon because we run into the thermal bath, which is visible as the conical singularity in the Euclidean space. The nature of this external bath becomes clear upon inspection of the holographic stress tensor, which implies that it can be identified as a radiation source. More specifically, separating the trace-free and trace parts, the energy-momentum tensor of the boundary theory is found to be [24]: The function F d (T, H) has the property that it vanishes for T = T dS . Thus, as we can see, the holographic stress tensor provides a smooth interpolation between the theory on de Sitter space at temperature T dS and finite temperature physics in Minkowski space. The extra piece that appears for even dimensions comes naturally upon integration of the gravitational conformal anomaly.

Gravitational collapse and black hole formation
Our goal is to study the approach to thermalization of entanglement entropy in a timedependent setup. The time-dependent configuration we are looking for here should capture the physics of gravitational collapse in the bulk which translates into the study of black hole formation. Thus, we will construct a generalized version of the Vaidya spacetimes for the so-called hyperbolic black holes.
To find the corresponding background, we have to couple the action in (2.1) with an external source S = S 0 + α S ext , (6.13) where α is a constant. For the physics we want to study in the present context we do not need to specify the form of S ext . The equations of motion in this case will take the following form (6.14) These field equations lead to the well-known AdS-Vaidya solution. For black holes with planar horizons, the relevant physics is restricted to the Poincaré patch and the metric takes the form This background has been used extensively in the program of holographic thermalization [40][41][42]. The metric (6.15) is written in terms of Eddington-Finkelstein coordinates (so that v labels ingoing null trajectories) and represents a (d + 1)-dimensional infalling shell of null dust. 8 To see this directly, let us analize the matter contents that leads to the above solution. The function m(v) is arbitrary and captures the information of the black hole formation. Quite generally, m(v) is chosen to interpolate between zero as v → −∞ (corresponding to pure AdS) and a constant value as v → ∞ (corresponding to AdS-Schwarzschild). A particular example of such a function is where v 0 is a parameter that fixes the thickness of the shell. With this choice, the external source must yield the following energy-momentum tensor If we identify k µ = δ µv , then we get [40] T ext µν ∼ k µ k ν , with k 2 = 0 , (6.19) which is characteristic of null dust. We can proceed in the same way and generalize the hyperbolic black holes (6.9) to include a time-dependent mass. We find it convenient to define an inverse radius z = L 2 /ρ, so that the boundary of AdS now lies at z → 0. Going to Eddington-Finkelstein coordinates, where the coordinate v is defined as we obtain 9,10 Similar to the planar AdS-Vaidya (6.15), the energy-momentum tensor of the external source yields which implies that the infalling shell is made of null dust. In order to avoid possible issues in the boundary theory, we want to constrain the mass function in order to satisfy the null energy condition in the bulk. 11 This implies in particular that m(v) should be an increasing function of v. In the context of the hyperbolic black holes we are considering here, however, there is a novel possibility that is not present in the traditional AdS-Vaidya backgrounds: we can start in a state with zero mass and end up with a finite positive mass or we can start with a negative mass (equal or above the extremal value) and finish with an increased mass (positive or negative). For concreteness, and with the aim of studying a situation of physical relevance, we will focus in the following mass function: where M ext ≤ M 1 < 0. In the boundary theory, this choice is equivalent to prepare the system in a state with T < T dS and then letting it evolve to the Bunch-Davies vacuum; in other words, the system undergoes a period of particle production and ends as a thermal bath of quanta at temperature T = T dS . To be even more specific, we will set so that the initial state is at zero temperature. We will see later that the relevant physics we are interested here is not affected by the specific value of M 1 . We can also consider a mass function of the type, with M 2 > 0. This corresponds to a situation in which we start in the Bunch-Davies vacuum and then evolve to a state with T > T dS . We will briefly comment on this possibility. In Figure 8 we plot the behavior of the mass functions m 1 (v)/|M 1 | and m 2 (v)/M 2 with respect to v for various values of v 0 . It is clear that in the thin shell limit, i.e. as v 0 → 0, the variation of the mass functions is sharply localized around v = 0 and approximates to a step function. Also, we have to bear in mind that the value of M 1 will depend on the dimension d according to (6.25). In the remainder of this section, we will study the evolution of entanglement entropy in these time-dependent backgrounds in the thin-shell approximation.

Holographic thermalization
We are now ready to compute the entanglement entropy in our dynamical background. For the ease of the computation, we will take as our entangling surface a spherical region of radius R inside the static patch, so that RH ≤ 1. We define a rescaled radius,r = rH, and parametrize the extremal surface by functions z(r) and v(r). The area functional in this case is given by Here the mass function is time-dependent and is given by either (6.24) or (6.26). The equations of motion resulting from the above area functional are quite involved and therefore we do not present them explicitly. Also note that we do not have any conservation equation since the action depends explicitly on r.
To solve these equations we impose the following boundary conditions z( ) = z * + corrections , z ( ) = 0 + corrections , where z * and v * are two free parameters, and is a small number. The "corrections" in the above expressions are obtained in the following way: we first write z(r) and v(r) as power series inr and then we expand the equations of motion aroundr = 0. Finally, we solve the resulting equations order by order and then we evaluate the solutions atr = .
So far z * and v * are two free parameters that generate the numerical solutions for z(r) and v(r). The boundary data can be obtained from where z 0 is an UV cut-off andt ≡ tH is the boundary time. In Figure 9 we plot sample numerical solutions for z(r) for fixed v * and various values of z * . Some of them cross the shell located at v = 0 and refract. For a fixedR this always happens for arbitrary early times (arbitrary negative v * ). However, the converse is not always true. In flat space the extremal surfaces always cross the shell for arbitrary large volumes. However, in dS space we need RH ≤ 1 in order to have connected solutions. Hence, if we fix v * , not always it is possible to find solutions that go deep enough into the bulk to cross the shell. To generate the corresponding thermalization curves S(t) for a fixedR, we do the following: for a fixed v * we vary z * until the read-off value ofR hits the desired value. This generates profiles for z(r) and v(x), which we can use to compute the area functional as well as the boundary time. Then, we vary v * and repeat the process again.
We are interested in the finite contribution to the entanglement entropy S(t) for a fixed radiusR. In order to study this quantity numerically, we fix the UV cutoff to be a small number, typically of the order of 10 −3 , and we define a rescaled entanglement entropy ∆S(t) by subtracting the entropy of the initial state ∆S(t) = S(t)−S(−∞). 12,13 Exploring the behavior of ∆S(t) as we change the radiusR, we note some general properties. Qualitatively, our results agree with those of [69,74] for Vaidya geometries with planar horizons. First, the evolution for very early times appears to be weakly dependent on the sphere size. The pre-local-equilibration regime is almost quadratic in time and it is followed by a post-local-equilibration linear growth phase at intermediate times. The entropy grows faster asR is increased. Finally, there is a period of memory loss prior to equilibration and then the entropy abruptly flattens out at late times after it reaches saturation. In Figure 10 (left panel) we plot this behavior for different values of the sphere radius. The saturation timet sat , on the other hand, shows a strong dependence on the sphere radius. It first increases linearly asR is increased, with a slope of order unity, and then blows up logarithmically asR approaches the horizon. This is surprising, it suggests the existence of light-like degrees of freedom that do not random walk, in agreement with the results of [43]. To see this, note that a radially outward light rays in de Sitter space obey the geodesic equation which can be solved and inverted to obtaint = tanh −1 (r). Forr 1 we gett =r+O(r 2 ), whereas forr ∼ 1 the leading contribution behaves liket = 1 2 (log(2) − log(1 − r)) + O(1 −r). This result holds true for arbitrary dimension. We verified numerically this behavior for different values of d = 2, 3, 4 and we found agreement to high accuracy. The results are plotted in Figure 10 (right panel). It is important to mention that this behavior might be particular for spherical regions. If this holds true for more general situations, one might conjecture that the saturation time is given by the time it takes for a light-ray to reach the farthest point of a particular entangling surface with respect to the observer. However, to prove this conjecture one would need a microscopic description of the thermalization process, which is beyond the scope of this paper.
Finally, we also explored the possibility of varying the initial and final masses of the black hole, corresponding to quenches evolving from and to different thermal states. One of such examples was the mass function given by equation (6.26). While the equilibration value of S sat showed a significant variation depending on the situation (and on the value of d), the behavior oft sat was found to be a robust feature, and quite insensitive to the details of the quench.

Conclusions and future directions
In this paper we have studied several properties of entanglement entropy in QFTs in de Sitter space, both in the static patch and the conformally flat patch. The theories in consideration have bulk duals coming from the standard Einstein-Hilbert action with negative cosmological constant and hence can be thought as models that belong to a universality class of strongly-coupled CFTs in the large-N limit. We started by choosing the standard vacuum state and obtained analytically extremal surfaces in the bulk with boundary condition taken as a spherical region of definite radius. According to the Ryu-Takayanagi prescription, the area of these solutions is interpreted as entanglement entropy for a spherical region in the boundary theory. Behaviors of extremal surfaces for R < H −1 and R > H −1 are qualitatively different. This implies that the entanglement entropy and renormalized entanglement entropy undergo phase transitions at R = H −1 . When R < H −1 , extremal surfaces are connected and U-shaped; whereas for R > H −1 , extremal surfaces are disconnected. We believe this is a bulk refection of a drastic decrease in entanglement at distances R > H −1 in the boundary theory. For realistic cosmological scenarios where a period of accelerated expansion is followed by decelerated expansion, regions outside the horizon will re-enter the horizon. It is reasonable to speculate that in that case, the extremal surfaces will go through transitions from disconnected to connected shapes and consequently the entanglement entropy will also undergo a reverse phase transition. We will investigate this transition in more detail in the future.
We have also studied the rich phase structure of entanglement/disentanglement phase-transition of mutual information in (1 + 1) dimensions. We found that mutual information between any two intervals separated by a proper distance x ≥ 2/H is identically zero, implying that for any two intervals A and B, separated by a distance larger than the horizon size, the density matrix ρ A∪B = ρ A ⊗ ρ B and hence they are completely disentangled. It is difficult to compute mutual information of spherical regions in higher dimensions, however, we expect qualitatively a similar behavior of mutual information in higher dimensions. A detailed calculation of mutual informations in higher dimensions would undoubtedly be useful to make this conclusion more robust.
We also considered out-of-equilibrium configurations in which the theory is undergoing a thermal quench. In order to do so, we constructed the Vaidya version of the well-known hyperbolic (or topological) black holes. In the static case, these solutions represent other states of the boundary theory (non-Bunch-Davies states) that are maintained externally in thermal equilibrium at a temperature T = T dS . When time dependence is introduced, we claim that these bulk backgrounds describe the dynamics of dS QFTs in a thermalizing, out-of-equilibrium state. We gave particular attention to situations in which the initial state is taken to be at zero temperature and ends as a thermal bath at temperature T = T dS , but we showed that our results hold for more general situations. For a fixed spherical region, the entanglement entropy follows a series of stages similar to what was found in [69,74] for theories living in flat space. The saturation time t sat , on the other hand, is found to depend on the sphere radius R: for RH 1 it increases linearly but then it blows up as RH → 1. This result accounts for the fact that, for a static observer, any signal that is sent radially outwards takes an infinitely amount of time to reach the horizon due to an infinite blueshift. The behavior of t sat is independent of the number of dimensions. We argue that this behavior could be accounted for if we think of the process of thermalization in terms of a streaming of light-like degrees of freedom.
Last but not least, it is worth emphasizing that the theories we are considering in the present paper are conformal field theories. In addition, we took spherical regions for the entangling surfaces which preserve the symmetries of the static patch. It would be interesting to investigate entanglement entropy for other shapes of the entangling surface, although it would be computationally more challenging. Another interesting possibility to address in the future would be to investigate the behavior of entanglement entropy in non-conformal theories on de Sitter space (see for instance [9,10,18,28]) and study the interplay between H and the different phase transitions.
Center, which is supported by the College of Natural Sciences and the Department of Astronomy at the University of Texas at Austin and the McDonald Observatory.