The first heat: production of entanglement entropy in the early universe

TL;DR

This paper investigates how entanglement entropy evolves in the early universe, showing rapid growth and saturation after transitioning from de Sitter to flat space, with implications for cosmological decoherence.

Contribution

It numerically demonstrates the transition of entanglement entropy from area law to volume law after the universe's transition, revealing a partially thermalized state.

Findings

Entanglement entropy grows rapidly after the de Sitter phase ends.

The entropy saturates at volume law scaling.

The final state resembles a generalized Gibbs ensemble.

Abstract

Entanglement entropy (EE) of a spatial region quantifies correlations between the region and its surroundings. For a free scalar in the adiabatic vacuum in de Sitter space the EE is known to remain low, scaling as the surface area of the region. Here, we study the evolution of entanglement after the universe transitions from de Sitter to flat space. We concentrate on the case of a massless minimally coupled scalar. We find numerically that, after the de Sitter stage ends, the EE and the R\'enyi entropy rapidly grow and saturate at values obeying the volume law. The final state of the subsystem (region) is a partially thermalized state reminiscent of a generalized Gibbs ensemble. We comment on application of our results to the question of when and how cosmological perturbations decohere.

Equations

$${\tilde{n}}_{\ell }=(e^{f_{\ell }}-1{)}^{-1}$$

$$S=\int dt\left(\frac{1}{2}{\tilde{a}}^2(t)\sum _i{\dot{\varphi }}_i^2-\frac{1}{2}\sum _{ij}{\varphi }_i{\tilde{K}}_{ij}(t){\varphi }_j\right),$$

$${\tilde{K}}_{ij}(t)=K_{ij}+m^2{\tilde{a}}^2(t){\delta }_{ij},$$

$${\mathcal{H}}_{BH}=\sum _i\left[-J\sum _e(a_i^{†}{a}_{i+e}+a_{i+e}^{†}{a}_i)+\frac{1}{2}U{n}_i^2-(\mu +\frac{1}{2}U)n_i\right],$$

$$n_i=\bar{n}+\delta n_i$$

$${\mathcal{H}}_{BH,2}=J\bar{n}\sum _{i,e}\left[({\theta }_{i+e}-{\theta }_i{)}^2+\frac{1}{4{\bar{n}}^2}(\delta n_{i+e}-\delta n_i{)}^2\right]+\frac{1}{2}U\sum _i\delta n_i^2.$$

$$\sum _{ij}{\theta }_i{L}_{ij}{\theta }_j=\sum _{i,e}({\theta }_i-{\theta }_{i+e}{)}^2.$$

$$Z_{ij}=\frac{1}{2J\bar{n}}\left({\delta }_{ij}+\frac{J}{2\bar{n}U}{L}_{ij}\right).$$

$${\mathcal{H}}_{BH,2}=\frac{1}{2}\varphi \left(\hat{L}+\frac{J}{2\bar{n}U}{\hat{L}}^2\right)\varphi +J\bar{n}U{\pi }^2.$$

$$\frac{d}{dt}\sum _i\delta n_i=-2J\bar{n}\sum _{ij}{L}_{ij}{\theta }_j=0.$$

$$S=\int d\eta a(\eta )\left(\frac{1}{2}\sum _i{{\varphi }_i^{\prime }}^2-\frac{1}{2}\sum _{ij}{\varphi }_i{\tilde{K}}_{ij}{\varphi }_j\right),$$

$${\pi }_i={\tilde{a}}^2(t){\dot{\varphi }}_i=a(\eta ){\varphi }_i^{\prime }$$

$$a(\eta)=\left\{\begin{array}[]{ll}-1/H\eta,&\eta\leq\eta_{0}\\ -1/H\eta_{0},&\eta>\eta_{0}\,,\end{array}\right.$$

$${\varphi }_i=\sum _k{q}_k{O}_{ki},$$

$$q_k(\eta )=\frac{1}{2}\sqrt{\pi H}\left[b_k\eta H_{\nu }^{(1)}(-{\omega }_k\eta )+b_k^{†}\eta H_{\nu }^{(2)}(-{\omega }_k\eta )\right],$$

$${\nu }^2=1-\frac{m^2}{H^2}.$$

$$q_k(\eta )=\frac{1}{\sqrt{2{\Omega }_{k}a({\eta }_0)}}\left(c_k{e}^{-i{\Omega }_k\eta }+c_k^{†}{e}^{i{\Omega }_k\eta }\right),$$

$${\Omega }_k^2={\omega }_k^2+m^2{a}^2({\eta }_0),$$

$${\omega }_{k,\mathrm{ph}}={\omega }_k/a({\eta }_0),$$

$${\omega }_k|{\eta }_0|={\omega }_{k,\mathrm{ph}}/H\ll 1.$$

$$t=H\eta /|{\eta }_0|+\text{mbox}const.$$

$$S_E$$

$$S_R(q)$$

$$C_{IJ}(\eta )=\frac{1}{2}⟨\{{\xi }_I(\eta ),{\xi }_J(\eta )\}⟩,$$

$$\xi =({\varphi }_1,\dots ,{\varphi }_n,{\pi }_1,\dots ,{\pi }_n),$$

$${\rho }_A(\eta )=\prod _{\ell =1}^n(1-e^{-{\beta }_{\ell }(\eta )})e^{-{\beta }_{\ell }(\eta )a_{\ell }^{†}(\eta )a_{\ell }(\eta )},$$

$${\beta }_{\ell }=\ln \frac{{\gamma }_{\ell }+1/2}{{\gamma }_{\ell }-1/2}.$$

$$s_E(\gamma )$$

$$s_R(q;\gamma )$$

$${\omega }_k^2=4-2\cos k_x-2\cos k_y,$$

$$b_E({\eta }_0)[2\mathrm{d}]$$

Full text

1 Introduction

Entanglement entropy (EE) offers a way to characterize thermalization in isolated quantum systems. At late times, for a generic nonlinear system, the EE of a subsystem is expected to approach the thermodynamic entropy but, as recent work has shown, surprisingly large amounts of it can be produced already at the early stages, well before the nonlinear effects have had time to set in. For example, in a lattice gas escaping from a small container into a larger one, the initial growth of the EE can be attributed to free streaming [1]. The growth of the EE after a global quench in an integrable field theory has also been ascribed to free streaming, in that case, of quasiparticles produced by the quench [2, 3, 4].

The examples recounted above correspond to particular choices of the initial out-of-equilibrium state. In the present paper, we look at yet another type of the initial state—a highly squeezed quantum state inherited from an earlier stage of the evolution. In cosmology, perhaps the best known example of this is the state of super-horizon perturbations that have been amplified (or, rather, have frozen in) during inflation [5, 6, 7, 8]. The analog of the post-quench dynamics in this case is the evolution of the perturbations after inflation, including the time of the second horizon crossing, when a given mode of perturbation reenters the horizon. We will assume that these large-wavelength perturbations remain in the linear regime, so they can be described by a free field theory.

The simplest model of the scale factor of the universe in this context is to connect, at some value $\eta_{0}$ of the conformal time, a de Sitter expansion, with the scale factor $a(\eta)=-1/H\eta$, where $H$ is a constant Hubble parameter, to a static flat universe with the scale factor $a(\eta_{0})$. One motivation for considering this problem is to understand better the “thermal” properties of cosmological perturbations. Earlier work [9, 10, 11, 12] has used various coarse-graining procedures to define the entropy of those, resulting in fairly large values—for a massless minimally coupled scalar in de Sitter, of order $H^{3}$ per unit physical volume. These results have been questioned later [13] as implying too much decoherence at the time of the second horizon crossing, compared to what is allowed by observations. Computations of the EE, on the other hand, do not rely on coarse graining. We therefore ask: can entanglement of a region with its surroundings lead to a finite density of the EE, and how large can that density be?

At first glance, the answer to the first question may seem an obvious no: a finite density implies an entropy scaling as the volume of an entangling region, while it is known that the EE of a free scalar in the adiabatic vacuum in de Sitter only scales as the surface area [14]. Suppose, however, that we think not of the EE in de Sitter itself but rather the one obtaining in a post-quench fashion after the de Sitter stage has ended. Studies of quenches in integrable theories [2, 3, 4] have shown that the entropy after a quench saturates at a value obeying the volume law, and we might expect the same for the present case. We start by verifying that expectation, using a numerical method.

Numerical calculations of the EE in Gaussian field theories are by now entirely standard (see, for example, Ref. [4] and references therein). The novel element in our case is the dependence of the EE on the degree of squeezing, as represented by the parameter $\eta_{0}$ defined above. The quantity we ultimately want is the entropy per degree of freedom (DOF) or, more precisely, per DOF that has been amplified during inflation. This can be computed on a lattice without a need to take the continuum limit. The reason is that the frequency of the lattice modes is bounded from above, so at a sufficiently small $|\eta_{0}|$ all the modes are amplified. To estimate the entropy per amplified mode, we can simply divide the entropy of a region by the number of lattice sites in it.

Fortuitously, a massless scalar in an expanding universe with two spatial dimensions is isomorphic in the infrared to the Bogoliubov phonon in an atomic superfluid with a variable coupling between the atoms. The latter system can then be considered as a “quantum simulator” for the former. In particular, a comparison of experimental data for the superfluid to the results presented here may allow one to isolate the effects of non-Gaussianity that were written into the quantum state during the time when the coupling was large. The question of how much information can be encoded in such non-Gaussianities is of interest for both the early universe and the physics of black holes. Because of this potential application of our results, we have chosen to tailor our presentation to the two-dimensional case, with an occasional summary of results for three dimensions.

In two spatial dimensions (2d), we compute the EE and the second Rényi entropy directly for rectangular regions of various sizes. In three dimensions (3d), we use the following modification. We consider, as the entangling region, a right rectangular column maximally extended along the third coordinate $z$ and compute the entropy of the modes independent of $z$, that is, having momentum $k_{z}=0$ in that direction. The number of these, for a column of square cross section with side length $L$ (in lattice units), is $L^{2}$. The volume law corresponds to the entropy per DOF being $L$-independent, i.e., the entropy of the column scaling as $L^{2}$. We find that, in both dimensionalities, the entropies of the regions grow rapidly after the end of the de Sitter stage, until they saturate at values obeying the volume law. The saturation entropies depend logarithmically on the parameter $\eta_{0}$. We discuss the significance of these logarithms in Sec. 4.

For a more precise characterization of the final state, we compute the distribution of the thermal parameters $\beta_{\ell}$ [15] corresponding to the Williamson normal modes. These are the modes for which the density matrix of the subsystem (region) is diagonal. Those parameters can be compared to an a priori different set of thermal parameters, $f_{\ell}$, which determine the occupation numbers

[TABLE]

of the Bogoliubov transformed operators corresponding to the subsystem’s usual normal modes. $f_{\ell}$ can be computed directly in de Sitter under the assumption that the subsystem is (quasi)isolated. They can be written as $f_{\ell}=f(\widetilde{\omega}_{\ell},\eta_{0})$ where $f$ is a known function of the normal frequency $\widetilde{\omega}_{\ell}$ and the degree of squeezing. We find that, at large squeezing, the distributions of $\beta_{\ell}$ (in the final state) and $f_{\ell}$ match fairly well. This suggests that the system undergoes a restricted version of thermalization: at large enough times, the Williamson normal modes of a subsystem and its usual normal modes become more or less the same.

The matching of the distributions of $\beta_{\ell}$ and $f_{\ell}$ lends credence to one of the coarse-graining proposal of Refs. [9, 10, 12], namely, coarse-graining in the occupation number basis, provided that it is applied not in de Sitter itself but rather after that stage has ended and the modes have been already oscillating for some time. The distributions do not match at earlier times. In addition, we recall that, unlike for coarse-graining, the density matrix in our case is that of a subsystem; the system as a whole remains in a pure state.

Matching of the distributions is also reminiscent of relaxation towards a generalized Gibbs ensemble (GGE) in integrable models [16, 17]. The difference is that, in our case, the particle numbers $N_{\ell}$, to which the expectation values (1) correspond, are not strictly conserved, because the subsystem is in fact not isolated. Our results are in keeping with the general principle of statistical physics, according to which the equilibrium density matrix of a subsystem should depend only on additive quasi-conserved quantities, i.e., those that would be conserved if it were not for exchange of them with the surroundings.

2 The lattice model

The main object of our study is a Gaussian scalar field $\phi_{i}$ defined on a two-dimensional (2d) square lattice (although we also discuss the 3d case). We consider two ways of defining dynamics for $\phi_{i}$. The first is by discretizing the action of a continuum relativistic scalar minimally coupled to gravity in a spatially flat expanding universe. The resulting action is

[TABLE]

where $\tilde{a}(t)$ is the scale factor of the universe,

[TABLE]

$K_{ij}$ is a time-independent symmetric positive semidefinite matrix, and $m$ is the mass of the scalar. Note that we have absorbed the area of the unit cell into $\tilde{a}^{2}$, so $\tilde{a}$ has dimension of time. We will present numerical results only for the massless case, when $\tilde{K}_{ij}=K_{ij}$, but some of the intermediate expressions apply more generally.

The second approach to defining dynamics is to start with the 2d Bose-Hubbard model in the superfluid phase,†††For the phase diagram of the Bose-Hubbard model, see Ref. [18]. where the scalar is the Bogoliubov phonon. On a square lattice, the Bose-Hubbard Hamiltonian is

[TABLE]

where $e$ goes over the positive lattice directions, $n_{i}=a_{i}^{\dagger}a_{i}$, and $J,U,\mu$ are parameters. We take $J>0$ and $U>0$, the latter corresponding to a repulsive interaction. In the superfluid phase, at low energies, the number density $n_{i}$ is close to its average $\bar{n}$, i.e.,

[TABLE]

where $\delta n_{i}$ is small. Defining the order parameter phase $\theta_{i}$ via $a_{i}=\sqrt{n_{i}}e^{-i\theta}$ and expanding in small $\delta n_{i}$ and $\theta_{i+e}-\theta_{i}$, one obtains the quadratic Hamiltonian

[TABLE]

Define symmetric matrix $L_{ij}$ such that

[TABLE]

To identify the Bogoliubov phonon, one rescales $\theta_{i}$ and $\delta n_{i}$ into a new canonical pair $\phi_{i}$, $\pi_{i}$, as follows: $\phi=\hat{Z}^{-1/2}\theta$, $\pi=\hat{Z}^{1/2}\delta n$ (in matrix notation), where $\hat{Z}$ is made of the following matrix elements:

[TABLE]

Suppose for a moment that the parameters in (8) are time-independent. Then, the new quadratic Hamiltonian is simply Eq. (6) rewritten in the new variables:

[TABLE]

To establish a connection with the cosmological model (2), we consider the case when $\bar{n}$ is large enough for the second term in the brackets in (8) and (9) to be always negligible. If we neglect this term and switch time-dependence on in $U$ but not in $J$ or $\bar{n}$, the matrix (8) will remain time-independent, and the expression (9) will still apply. The theory then becomes equivalent to one with the action (2) with zero mass and the following identifications: $\hat{K}=\hat{L}$ and $\tilde{a}^{2}(t)=[2J\bar{n}U(t)]^{-1}$.

For a possible experimental realization, it is important that a simulation of the theory (2) with atoms in an optical lattice requires only time-dependence of $U$ and not, for instance, of the hopping matrix element $J$. This is specific to two spatial dimensions.

The sum $\delta N=\sum_{i}\delta n_{i}$ represents a variation of the total particle number and, if the latter is conserved, must be zero. The Hamiltonian (6) guarantees that, if imposed initially, this condition will hold at all times:

[TABLE]

This is because $\sum_{i}L_{ij}=0$, meaning that the symmetric matrix $\hat{L}$ has a constant as its zero mode. In the continuum limit, the particle number is conserved by the Neumann ($\nabla\theta=0$) boundary condition on $\theta$ but not by the Dirichlet (fixed $\theta$) one. Eq. (2) with $\hat{K}=\hat{L}$ can be viewed as a discretization of a continuum action with the field subject to the Neumann condition.

In what follows, we remove the constant mode from the consideration entirely by assuming that we work with the ensemble with a fixed total particle number.

Upon transition to the conformal time $\eta$ via $d\eta=dt/\tilde{a}(t)$, the action becomes

[TABLE]

where $a(\eta)\equiv\tilde{a}[t(\eta)]$, and a prime denotes a derivative with respect to $\eta$. Recall that $\tilde{a}$ has dimension of time, so $\eta$ is dimensionless.

The canonical momentum conjugate to $\phi_{i}$ is

[TABLE]

and can be obtained from either (2) or (11).

We adopt the following form of the scale factor

[TABLE]

corresponding to a de Sitter expansion that connects, at some $\eta=\eta_{0}<0$, to a flat spacetime.

To obtain the mode expansion of the Gaussian field in this spacetime, we first diagonalize the matrix $\hat{K}$ by an orthogonal transformation

[TABLE]

resulting in a set of eigenvalues $\lambda_{k}=\omega_{k}^{2}$ (as noted above, we only keep the nonzero ones) and the corresponding eigenvectors (the normal modes), which are the rows of the orthogonal matrix $\hat{O}$; $q_{k}$ are the new canonical coordinates.

The time-dependence of the amplitudes $q_{k}$ is obtained by solving a linear differential equation. The result is as follows. At $\eta<\eta_{0}$,

[TABLE]

where $b_{k}^{\dagger}$ are time-independent creation operators, $H_{\nu}$ are the Hankel functions, and

[TABLE]

The latter is for the general massive case (3). The normalization coefficient in (15) is chosen so that $q_{k}$ and $p_{k}=a(\eta)q^{\prime}_{k}$ are correct canonical pairs.

At $\eta>\eta_{0}$,

[TABLE]

where

[TABLE]

and $c_{k}^{\dagger}$ is another set of time-independent operators. By matching the expansions (15) and (17) at $\eta=\eta_{0}$, one expresses one set of operators in terms of the other.

At fixed $\nu$, the only dependence of Eqs. (15) and (17) on the Hubble parameter is in the prefactors: $q_{k}$ scales as $\sqrt{H}$ and its conjugate momentum $p_{k}$ as $1/\sqrt{H}$. This dependence can be removed by a canonical transformation without affecting the value of the entropy. So, in computations of the entropy we can set $H=1$. When we wish to express results in terms of the physical frequency

[TABLE]

we will restore $H$. For example, the condition that a given mode is significantly amplified in de Sitter can be written as

[TABLE]

Another place where $H$ plays a role is the relation between the conformal time $\eta$ and the laboratory time $t$ during the flat space segment of the evolution:

[TABLE]

Note that, for a given $H$, the same interval of $\eta$ corresponds to different intervals of $t$ for different values of $\eta_{0}$.

3 Entanglement entropy

At $\eta\to-\infty$, all modes (on a finite lattice, with the constant mode excluded) evolve adiabatically, and the vacuum of the operators $b_{k}$ from (15) is the usual adiabatic vacuum of a free scalar in de Sitter. This is the state for which we will be studying correlations as represented by the EE and the Rényi entropy for different spatial regions. The regions are defined by selecting a set of lattice sites and are time independent. We refer to these regions as subsystems.

To obtain the entropies, we first construct the density matrix $\rho_{A}$ of a region $A$ by tracing out the local degrees of freedom $(\phi_{i},\pi_{i})$ corresponding to lattice sites outside the region. The EE and the $q$-th Rényi entropy of the subsystem are then defined by the usual formulas

[TABLE]

A convenient way to define $\rho_{A}$ for a Gaussian state is to require that it reproduces the covariance matrix

[TABLE]

where the braces denote anticommutator, and the vector

[TABLE]

whose components $\xi_{I}$ appear in (24), combines the coordinates and momenta for the lattice sites inside the region $A$. The averaging is over the adiabatic vacuum. For a Gaussian state with zero $\langle\phi_{i}\rangle$ and $\langle\pi_{i}\rangle$, such as the one here, the covariance matrix carries all the information there is about the state. The averages in (24) are then reinterpreted as those computed in a mixed state in a Gaussian theory that contains only the variables in $A$, with the help of a density matrix $\rho_{A}\equiv\rho_{A}(\eta)$.

The $2n\times 2n$ matrix $\hat{C}$ that has $C_{IJ}$ as its matrix elements has the following normal form, called the Williamson normal form [19]: $\hat{C}=\hat{{\cal S}}^{T}\mbox{diag}(\bm{\gamma},\bm{\gamma})\hat{{\cal S}}$, where $\hat{{\cal S}}$ is a symplectic matrix. The $n$ entries $\gamma_{\ell}$ that comprise the vector $\bm{\gamma}$ are called the symplectic eigenvalues of $\hat{C}$. The density matrix $\rho_{A}(\eta)$ can be written as a product of thermal density matrices for a set of oscillators $a_{\ell}(\eta)$, $a_{\ell}^{\dagger}(\eta)$ corresponding to the Williamson normal modes, as follows [15]:

[TABLE]

where

[TABLE]

The entropies (22) and (23) are then the sums of those of the individual oscillators, that is, $S_{E}=\sum_{\ell}s_{E}(\gamma_{\ell})$ and $S_{R}(q)=\sum_{\ell}s_{R}(q;\gamma_{\ell})$, where

[TABLE]

and $\gamma_{\pm}=\gamma\pm{1\over 2}$. This method was used for a calculation of the entropy growth after a quench in a Gaussian theory in flat spacetime in Ref. [4].

We now present numerical results for the model (11) on a square 2d lattice with $K_{ij}$ set equal to $L_{ij}$, Eq. (7), and $a(\eta)$ of the form (13). The eigenvalues of this $\hat{K}$, for an $M\times N$ lattice, are

[TABLE]

where $k_{x}$ is quantized in units of $\pi/M$, and $k_{y}$ in units of $\pi/N$. Note that $\omega_{k}$ is bounded from above by $2\sqrt{2}$, so for a sufficiently small $|\eta_{0}|$ all lattice modes are amplified.

With the Hubble parameter set to one as described at the end of Sec. 2, the only parameters left in the model are the mass parameter $\nu$, Eq. (16), and the transition time $\eta_{0}$ in Eq. (13). We can also vary the total size of the system and the size and shape of the entangling region. Here, we present results for the massless case $\nu=1$ and the total lattice size of $31\times 33$.

Fig. 1 shows the EE and the second ($q=2$) Rényi entropy for a $10\times 10$ square in the middle of the system for two different values of $\eta_{0}$. We see that, after the expansion is switched off, the entropies grow rapidly, until they saturate at values that are large compared to the entropy the region had in de Sitter.

The dependence of the saturation entropy on the subsystem linear size $L$ is shown in Fig. 2. One finds a good volume scaling with a coefficient of $L^{2}$ that depends on $\eta_{0}$. That dependence, for $|\eta_{0}|\leq 0.01$, is well fit by

[TABLE]

for the EE and the second Rényi entropy, respectively. Because for these small $|\eta_{0}|$ all modes we have on the lattice are amplified, we can view (31) and (32) as estimates of the entropies per degree of freedom (DOF) of the subsystem.

For a parallel computation in three dimensions, we consider an entangling region in the form of a right vertical column that extends for all values of the third ($z$) coordinate and has a constant cross-section of $L\times L$ sites in the $(x,y)$ plane. To estimate the entropy per an amplified DOF, we restrict attention to the $k_{z}=0$ mode of the field. Then, the total number of DOF in the subsystem is $L^{2}$, just as in the 2d case. The difference with that case is in the form of the mode functions at $\eta<\eta_{0}$: instead of (15) we now have

[TABLE]

where

[TABLE]

The results are similar to those for two dimensions. For the massless case ($\nu=3/2$), the counterparts of Eqs. (31) and (32) are

[TABLE]

for $|\eta_{0}|\leq 0.01$.

Finally, we record an estimate for the entropies per unit physical volume, corresponding to the saturation values (35) and (36). In terms of the physical frequency (19), the logarithms in (35)–(36) become those of $\omega_{k,{\rm ph}}/H$, so

[TABLE]

and similarly for $S_{R}$. The step function $\theta$ restricts the integral to the amplified modes only. In the continuum limit, this subtracts the vacuum contribution of the high-frequency modes. In that case, $\omega_{k,{\rm ph}}=k_{{\rm ph}}$, and the entropy densities are of order $H^{3}$ (in units where the speed of quasiparticles is equal to 1).

4 Thermal interpretation

The coefficients of $\ln|\eta_{0}|$ in the fits of Sec. 3, for both dimensionalities, are noteworthy. Earlier papers [9, 10, 11, 12] have sought to define entropy of a scalar field in de Sitter by means of a coarse-graining procedure, i.e., by neglecting the off-diagonal elements of the density matrix in a preferred basis. The density matrix in question is that of the entire system; in our case, that would involve the full $M\times N$ lattice sites. The results of such coarse-graining depend on which basis is chosen. One of the choices amounts to using the ideal gas expression

[TABLE]

where $k$ goes over the normal modes of the full system (as before, we exclude the constant), $n_{k}=\sinh^{2}r_{k}$, and $r_{k}$ is the squeezing parameter. According to (38), a mode with a large $|r_{k}|$ contributes $s_{k}\approx\ln n_{k}$ of entropy to the total.

We list here expressions for $n_{k}$ for a massless minimally coupled scalar in two and three spatial dimensions. In 3d, $n_{k}=1/(2\omega_{k}\eta)^{2}$ [20], so in the limit of large squeezing $s_{k}\approx-2\ln(\omega_{k}|\eta|)$. In two dimensions, the general expression for $n_{k}$ is more complicated, but at large $r_{k}$ it simplifies into $n_{k}\approx 1/|2\pi\omega_{k}\eta|$, so $s_{k}\approx-\ln(\omega_{k}|\eta|)$. (The coefficient of the logarithm is half of that in 3d.) These results are obtained directly in de Sitter and make no reference to the time $\eta_{0}$ at which it ends. If, however, we substitute $\eta=\eta_{0}$, they match Eqs. (35) and (31).‡‡‡Curiously, the expression for $n_{k}$ in two dimensions, when written in terms of the physical frequency (19) as $n_{k}\approx H/(2\pi\omega_{k,{\rm ph}})$, coincides with the low-frequency tail of the Planck distribution with the Gibbons-Hawking [21] temperature $T_{GH}=H/2\pi$.

We find this correspondence nontrivial. For one thing, our computations do not involve any coarse graining. For another, Eq. (38) refers to the normal modes of the full system, while our results to a subsystem. Finally, our expressions (31) and (35) are for the saturation entropy, which obtains after both the de Sitter stage and the subsequent growth of the EE have already ended. During the de Sitter stage itself, the EE is much smaller (cf. Fig. 1). Indeed, at that time, it does not even scale as the volume of the region, only as the surface area [14]. If we think of entropy as representing decoherence (i.e., the loss of information on the relative phases of the basis states), we have to conclude that fluctuations of the field do not decohere significantly while they are outside the horizon but only do so after they start oscillating. This, in particular, may help resolve the tension [13] between the large value of $S_{c.g.}$ and the amount of coherence between the first and second horizon crossings that is required to explain the observational data.

The occupation number $n_{k}$ at the end of the de Sitter stage is a definite function of the mode frequency $\omega_{k}$ and the parameter $\eta_{0}$. If the subsystem were isolated, the occupation numbers $\tilde{n}_{\ell}$ of its normal modes would be given by the same function, but now of the subsystem’s own normal frequencies $\widetilde{\omega}_{\ell}$. For a more detailed characterization of the final state, we compare the spectrum of the parameter $f_{\ell}$ obtained from these $\tilde{n}_{\ell}$ via Eq. (1) to that of the thermal parameter $\beta_{\ell}$, Eq. (27). We concentrate on the limit of large squeezing, when $f_{\ell}\approx 1/\tilde{n}_{\ell}$. If $\tilde{n}_{\ell}$ obeyed the Planck distribution at temperature $T$, we would have $f_{\ell}=\widetilde{\omega}_{\ell}/T$. As it is, according to the expressions listed after Eq. (38), $f_{\ell}\approx 2\pi\widetilde{\omega}_{\ell}|\eta_{0}|$ in two dimensions and $f_{\ell}\approx 4\widetilde{\omega}_{\ell}^{2}\eta_{0}^{2}$ in three.

Even if were assured of thermalization, we could not expect a perfect correspondence between the individual $\beta_{\ell}$ and $f_{\ell}$. Indeed, $f_{\ell}$ are computed under the assumption that the subsystem is perfectly isolated, while in practice it is not. So, for a comparison, we bin $\beta_{\ell}$ and $f_{\ell}$ into distributions $D_{\beta}(\eta)$ and $D_{f}$, which count the number of times the corresponding parameter ($\beta_{\ell}$ or $f_{\ell}$) falls into a bin centered on a given value. Note that $D_{\beta}$ depends on time, while $D_{f}$, being fully determined by $\eta_{0}$ and the spectrum of the normal modes, does not. A sample result for two dimensions is shown in Fig. 3. Results for three dimensions are similar. We wish to stress that a good agreement between the two distributions, seen in Fig. 3, develops only at sufficiently large times. At $\eta=0$, $\beta_{\ell}$ are distributed over a much broader range. As discussed in the Introduction, we interpret the agreement at large times as a restricted version of thermalization.

5 Conclusion

Our main conclusion is that a thermal interpretation (along the lines of Sec. 4) of the entanglement entropy of a Gaussian scalar in a region of the de Sitter spacetime looks plausible, provided one applies it not to fluctuations in that spacetime itself but rather to those that evolve after the de Sitter stage has ended. It is as if the spacetime stores and then releases latent heat. If experimental results for the quantum simulator (the Bose-Hubbard model) described in Sec. 2 become available, comparing them to the results obtained here for the Gaussian theory may elucidate the role of nonlinearities.

Acknowledgments

We thank Chen-Lung Hung and Martin Kruczenski for discussions. This work was supported in part by the U.S. Department of Energy (grant DE-SC0019202) and by the W. M. Keck Foundation.