Entanglement in cosmology

We compute the evolution of the entanglement entropy for a massless field within a spherical region throughout the inflationary period and the subsequent era of radiation domination, starting from the Bunch-Davies vacuum. In order to focus on the entanglement of modes that are directly accessible to observations, we impose an ultraviolet cutoff set by the wavelength of the last mode that exited the horizon at the end of inflation. The transition of each mode towards a squeezed state upon horizon exit during inflation and the additional squeezing when radiation domination sets in enhance the entanglement entropy. Shortly after the transition to the radiation-dominated era, a volume term develops and becomes the leading contribution to the entropy at late times, as is common for systems lying in squeezed states. We estimate the magnitude of the entropy and discuss its interpretation in the light of the quantum to classical transition for modes exiting the horizon during inflation. Our results raise the possibility that the quantum nature of weakly interacting fields, such as gravitational waves resulting from tensor modes during inflation, may be detectable in today's universe. On the other hand, an observer with no knowledge of the degrees of freedom beyond the horizon would interpret the entropy as thermal. From this point of view, the reheating after inflation would be a result of quantum entanglement.


Introduction
According to the inflationary paradigm, the large-scale structure of our universe originated in vacuum fluctuations during inflation [1][2][3][4][5].Quantum field fluctuations that were stretched beyond the horizon by the expansion were transformed into classical stochastic fluctuations.After the end of inflation, they re-entered the horizon at successive times depending on their wavelengths, and generated the observed structure through the process of gravitational collapse of overdense regions.
A rather obscure point in the above scenario is the quantum to classical transition upon horizon exit.In inflationary cosmology a field fluctuation can be expressed in terms of momentum modes whose mode function involves a growing and a decaying term.After horizon exit, the mode function becomes dominated by the first term and loses its oscillatory form (it freezes) [6].Moreover, the dominance of the growing contribution causes the field and its conjugate momentum to commute.As a result, the field after horizon crossing is viewed as a classical stochastic field, and its quantum expectation value is considered as the classical stochastic average.Despite the simplicity of the above picture, it is important to keep in mind that the full quantum field and its conjugate momentum always obey the canonical commutation relation.This is guaranteed by the presence of the decaying term in the mode function.In this sense, the field never loses its quantum nature.It has been argued that the observation of quantum properties of the field that may survive until today is very difficult because of the enormous difference in the amplitudes of the growing and the decaying term, which would require a precision of 90 orders of magnitude in the measurement of the momentum [7].Devising experiments that could look for a violation of Bell's inequalities in the cosmological context seems like a formidable task [7][8][9].
At the conceptual level, the question persists as to whether the quantum origin of the fluctuations can be mirrored in certain quantities that probe beyond the classical level.A very interesting example of such a quantity is the entanglement entropy.In the simplest approach, one may consider the momentum-space entanglement between high and lowmomentum modes, such as between modes with physical momenta below and above the Hubble scale H [10][11][12][13].However, this would vanish for a free field, as long as the initial state can be written as a tensor product of states, one for each momentum mode, as in the Minkowski vacuum.Since each mode evolves independently, the reduced density matrix, resulting from some modes being traced over, would be that of a pure state.A nontrivial result is obtained only in the presence of an effective coupling between different momentum modes, as in [14].
The entropy associated with the entanglement between degrees of freedom localized within two spatial regions separated by an entangling surface is more promising.The calculation is more difficult, as one now has to trace over the degrees of freedom in the interior or the exterior.On the other hand, the reduced density matrix would not correspond to a pure state and the entanglement entropy would be nonvanishing even in free field theory.Explicit calculations of real-space entanglement entropy in flat space have been carried out for non-interacting or highly symmetric quantum field theories, and mostly in lower dimensions [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33].The entropy has interesting properties, such as a dependence on the area of the entangling surface, which indicates a similarity with black hole entropy [15,16].
There are few results on the behaviour of real-space entanglement entropy in a cosmological setting [27,34,35].The explicit calculation of [27] focused on the subleading logarithmic term, whose coefficient is universal.We are interested instead in the leading term, for which we seek a physical interpretation.Our approach follows closely the calculation of the entanglement entropy in flat (3+1)-dimensional spacetime in the seminal work by Srednicki [16].The original calculation focused on the leading term in the entropy, while the logarithmic correction was computed in [24].The basic formalism for the generalization of the calculation to the case of a Friedmann-Robertson-Walker (FRW) background was developed in [35,36] and was applied to a toy (1+1)-dimensional model.In this work we carry out the analysis for a (3+1)-dimensional cosmological background that undergoes a transition from an inflationary era to a period of radiation domination.
It must be mentioned also that the entanglement entropy can be derived through the Ryu-Takayanagi proposal [37][38][39] in the context of the AdS/CFT correspondence [40][41][42], and results have been obtained in time-dependent backgrounds for theories that have gravitational duals [43,44].Our aim here is to look directly at the details of the mechanism of entanglement in a cosmological setting for fields minimally coupled to gravity.For this, we perform an explicit calculation of the reduced density matrix in the case of a single massless scalar field, from which we deduce the entropy.The effect of a non-minimal coupling has been explored in [45,46].
The main effect of the expansion on the quantum mechanics of the scalar field is that the field modes evolve from a simple oscillator ground state to a squeezed state [47].There is a long history of studies on the connection between squeezing and entropy [48][49][50][51][52][53][54][55].As we discussed above, the dominance of the growing term usually leads to the conclusion that the quantum properties of the field are not visible in late-time observations that focus on classical local quantities [56][57][58][59][60][61].However, the entanglement entropy may provide an exception, as it is a purely quantum, non-local quantity that does not have a classical analogue.More specifically, the squeezing of canonical modes generically increases the entanglement between local degrees of freedom and is expected to also increase the entropy.
A conceptual question that arises when trying to interpret the entanglement entropy is how to control the ultraviolet (UV) divergence that it displays.In flat (3+1)-dimensional spacetime the entropy scales ∼ 1/ϵ 2 , with ϵ a short-distance cutoff.For the vacuum state this cutoff is intrinsically connected with the radius of the entangling surface R, as they appear together in the leading term ∼ R 2 /ϵ 2 that realizes the area law.The standard interpretation is that this term quantifies the very strong entanglement between the short-distance modes on either side of the entangling surface.
It is difficult to remove the divergent term through some kind of renormalization procedure.On the other hand, if it is viewed as physical there is an ambiguity in the identification of the length scale ϵ.In a theory that includes gravity, one may adopt the logic that ϵ must be of the order of the Planck length.However, there are fundamental difficulties with such an assumption.The huge difference between the Planck scale and macroscopic scales of interest would mean that the entanglement entropy would always be dominated by the UV and would not carry any interesting information about the long-distance physics.For an expanding background the situation is more problematic.The continuous stretching of physical length scales and the corresponding redshifting of frequencies imply that the UV is continuously replenished by new modes that emerge from sub-Planckian distances.One can only make arbitrary assumptions about the state that such modes occupy, as well as their entanglement with the longer ones.
The quantum to classical transition upon horizon crossing during inflation can give a natural way to bypass the above fundamental issues.There is a certain mode of comoving wavenumber k s which crossed the horizon at the end of inflation and immediately re-entered.Modes with wavenumbers k > k s remained subhorizon at all times and never went through the process of freezing and the dominance of the growing term in the mode function.In this sense, they have always constituted vacuum fluctuations.Of course, the value of k s is not fixed precisely, but the freezing of adiabatic modes occurs sufficiently fast for k s to be determined up to a factor of order 1.The modes with k < k s are the ones directly accessible to experiment and constitute the observable universe.Even though they appear mostly classical, their quantum nature may still be visible in quantities such as the entanglement entropy.The advantage of this logic is that the entanglement of interest is due to modes with wavelengths above a UV cutoff ϵ ∼ 1/k s .
The absence of modes with k > k s can be justified only for very weakly interacting fields, for which mode-mode coupling can be ignored.For interacting fields, all modes are eventually excited.This does not result in an additional limitation on the type of fields that are possibly entangled today, as interacting fields are expected to thermalize and lose quantum coherence anyway.Our analysis is relevant for very weakly interacting fields as, for example, gravitational waves resulting from tensor modes during inflation.We treat them as freely evolving fields on the cosmological background, and we make the underlying assumption throughout the paper that their quantum coherence is not lost through some secondary process during the whole evolution of the universe until today.
In order to obtain a feeling of the magnitude of k s , we may use the fact that a gravitational wave generated N e-foldings before the end of inflation, with frequency set by the approximately constant Hubble parameter H, is redshifted to a frequency f today [62] where N CMB ≃ 60 indicates the number of e-foldings at which the modes of the cosmic microwave background (CMB) exit the horizon.Modes that exit the horizon at the end of inflation (N ≃ 0) would have a frequency today f ∼ 10 8 Hz, which sets the cutoff in the spectrum of gravitational waves generated by inflation [62].The corresponding wavelength is λ s ∼ 1m.This is a scale much longer than any conceivable fundamental UV cutoff, but still much shorter than the Hubble radius today.It is interesting also to note that a possible infrared cutoff k l may exist for the modes contributing to the entanglement entropy.Let us assume that inflation started at some initial time, after some unspecified pre-inflationary era.Even if modes that were superhorizon at that time went through the process of freezing later on, their enhancement will be smaller than if they had grown throughout inflation starting from below the horizon.The power spectrum for such modes will be suppressed relatively to its value in the scale-invariant range.It is then natural to expect that these modes do not contribute significantly to the entanglement entropy.This, admittedly rather speculative, argument suggests that k l can be identified with the wavenumber of the mode that exited the horizon at the beginning of inflation.
In the following we focus on the scenario in which the UV cutoff is identified with the last mode that crossed the horizon during inflation.Its wavenumber satisfies k s /a t ≃ H, where H is the constant value of the Hubble parameter during inflation, and a t is the value of the scale factor at the transition from inflation to the era of radiation domination.The UV regularization is implemented by considering a discretized version of the free scalar theory on a lattice of comoving spacing ϵ.The highest mode has a comoving wavenumber ∼ 1/ϵ.Without loss of generality we can set a t = 1.This results in the condition k s ≃ H, or ϵ ∼ 1/H, where H is the constant value of the Hubble parameter during inflation.With these assumptions, the physical lattice spacing ϵa never exceeds the physical Hubble radius during the whole evolution from the de Sitter (dS) phase to the era of radiation domination (RD).
The paper is organized as follows: In section 2 we summarize the formalism that is needed in order to compute the entanglement entropy in a time-dependent background.We describe the discretized version of the theory of a massless scalar field and the form of the wave function of the canonical modes.We also provide a concise summary of the method developed in refs.[35,36] for the calculation of the reduced density matrix and its eigenvalues, from which the entanglement entropy is deduced.In section 3 we present the results of a numerical calculation of the entanglement entropy in a time-dependent background that evolves from an inflationary era to radiation domination.In section 4 we give our conclusions.

Formalism
In this section we summarize the formalism that we employ for the calculation of the entanglement entropy.We follow the approach of Srednicki [16], which was generalized to the case of an expanding background in [35,36].We refer the reader to these publications for the details.

The discretized Hamiltonian
We consider a free real scalar field in a FRW background, described by the metric in terms of the conformal time τ .Through the definition ϕ(τ, x) = f (τ, x)/a(τ ), the action of the field can be written as where the prime denotes differentiation with respect to the conformal time τ .The field f (τ, x) has a canonically normalized kinetic term.As we consider spherical entangling surfaces, it is convenient to define the spherical moments of the field and its momentum as where Y lm are real spherical harmonics.The radial coordinate can be discretized by introducing a lattice of concentric spherical shells with radii r j = jϵ, where j is an integer obeying 1 ≤ j ≤ N .The radial distance between successive shells introduces an UV cutoff equal to 1/ϵ, while the total size of the lattice L = N ϵ sets an IR cutoff equal to 1/L.By defining the discretized degrees of freedom as so they are canonically commuting, we arrive at the Hamiltonian All dimensionful quantities in the problem can be expressed in units of the comoving lattice spacing ϵ.This is also the case for the time parameter τ , as the combination Hτ that determines the time evolution becomes a function of τ /ϵ, as is apparent from eq. (2.5).This normalization corresponds to setting ϵ = 1.There is a correspondence between ϵ and the comoving scale k s we discussed in the introduction, i.e. ϵ ∼ 1/k s .Similarly, we have that We shall analyze the scenario in which the evolution of the field starts during inflation in a de Sitter (dS) phase and continues in the era of radiation domination (RD).We approximate the transition as instantaneous, neglecting the complications of reheating.It is also possible to consider the transition to an era of matter domination (MD), either directly from the dS phase [35] or from the RD era.We have considered such scenarios as well, but they lead to results that are very similar to the scenario we discuss in detail, albeit with much more complicated analytical expressions.Assuming that a transition from a dS to a RD background takes place at τ = 0, the scale factor takes the form so that The Hubble parameter is continuous through τ = 0.However, the effective mass term −a ′′ /a in eq.(2.5) is discontinuous.
The entanglement entropy between the interior and exterior of a sphere of radius R can be obtained from the density matrix of the harmonic system that corresponds to the discretized field theory, via tracing out the oscillators with jϵ < R in order to compute the reduced density matrix for the degrees of freedom outside the entangling surface (or vice versa).For the vacuum of the theory, the state of the system of oscillators is the product of the 'ground states' of the modes that diagonalize the Hamiltonian.The assumption of a Bunch-Davies vacuum implies that as 'ground state' of a mode we must define the solution of the time-dependent Schrödinger equation for this mode that reduces to the usual simple harmonic oscillator ground state as τ → −∞.The wave function of each mode depends on a linear combination of the various f lm,j , i.e. the corresponding canonical coordinate.Since modes with different l and m indices do not mix, each eigenfunction actually involves one set of (l, m).The effect of the expanding background is encoded in the term ∼ (a ′′ /a)f 2 lm,j , which is identical for all (l, m).For a background that evolves from a dS to a RD era, the discretized Hamiltonian for the free field takes the form where flm,j are the canonical coordinates (linear combinations of f lm,j ), ω lm,j the corresponding eigenfrequencies when the time dependence is neglected, and At this point, we need to mention a subtle issue concerning the cutoff imposed on the eigenfrequencies of the canonical modes of the above Hamiltonian.In the original calculation by Srednicki [16], all discrete values of l = 0, ..., ∞ were taken into account.For a given value of ϵ, the sum over m and l appearing in eq.(2.8) converges in the calculation of the entanglement entropy, so that the UV divergence is well controlled by ϵ.In our approach, however, the fundamental issue is not simply one of regularization, but of the exclusion of all physical modes that remain subhorizon during the whole evolution.Allowing l to take arbitrarily large values incorporates fluctuations in the tangential directions that should have been excluded.In order to resolve this issue we have adopted the following strategy: We first determine the maximal eigenfrequency ω max in the l = 0 sector, which is of order 1/ϵ.In all other l sectors, we exclude the eigenfrequencies that exceed ω max , and the corresponding canonical modes.The resulting entanglement entropy still has a UV cutoff set by ϵ, but a difference may appear in nonuniversal contributions.In particular, the leading term ∼ R 2 /ϵ 2 for a static background will have a different coefficient than the one computed in [16].Actually, the coefficient turns out to be smaller than in [16], as an infinite number of modes are excluded.The procedure may be viewed as a different regularization of the entanglement entropy.However, our interpretation is that it provides a result that incorporates only the modes that went through the physical process of freezing upon horizon exit.
The cutoff procedure is also related to the presence of a new scale in the calculation of the entropy, set by the Hubble constant H during inflation.In the introduction we argued that the relevant scale of highest frequency for the problem is the one with comoving wavenumber k s , such that the mode exits the horizon at the end of inflation and immediately re-enters in the RD era.For this mode k s /a t = H, where a t is the value of the scale factor at the transition from the dS to the RD era.We have set a t = 1 without loss of generality, while the shortest mode that we keep in the calculation has an eigenfrequency set by ϵ.This leads to the identification ϵ ∼ 1/H, which we employ in the next section.
Of course, one could ignore the issue of mode freezing and treat Hϵ as a free parameter.In this case, the scale 1/ϵ can be viewed as a coarse-graining scale in a Wilsonian procedure of integrating out high-energy modes.As the momentum modes of a free field are decoupled, the integration would result again in a free theory, albeit with fewer modes.Even though this logic seems appealing in momentum space, it is less transparent in real space.The degrees of freedom in this case are the field values at every point in the discretized space (the lattice), which are coupled through the kinetic term.Varying the lattice spacing essentially defines a new entropy that depends on an arbitrary UV scale.Moreover, Hϵ ≫ 1 results in a discretized version of the field theory such that the characteristic length scale of the gravitational background can become smaller than the physical lattice spacing during part of the evolution.On the other hand, Hϵ ≪ 1 results in the dominance of the UV fluctuations and a result for the entanglement entropy identical to that in static space at all times.The natural choice Hϵ ∼ 1 that we advocate avoids the above complications.

The wave function of the canonical modes
The quantization procedure treats the field as a collection of quantum harmonic oscillators.In order to make the analogy with the language of quantum mechanics clearer, we follow the notation of previous work [16,35,36] and use the variable x instead of f .It must be kept in mind that it is the field value that is treated as a harmonic oscillator.In this section the variable x has nothing to do with spatial coordinates.The Hamiltonian (2.8) becomes diagonal when expressed in terms of canonical modes, which must be put at their 'ground state' in the vacuum of the theory.The form of the Hamiltonian implies that we need the 'ground state' eigenfunction of the harmonic oscillator with a time-dependent frequency of the form (2.10) We set ϵ = 1 for simplicity in the expressions.Since the Hamiltonian has explicit time dependence, this 'ground state' is not an eigenstate of the Hamiltonian, but rather a solution of the time-dependent Schrödinger equation that reduces to the ground state of the simple harmonic oscillator with constant frequency ω 0 , as τ → −∞.
A solution of the time-dependent Schrödinger equation for an oscillator with frequency given by eq.(2.10) can be obtained in several steps [63], following the Lewis-Riesenfeld method [64,65].First one must find a solution b(τ ) of the Ermakov equation with boundary conditions appropriate to the problem at hand.In terms of b(τ ), the solution of the time-dependent Schrödinger equation is given by where F 0 (τ, x) is a solution of the standard simple harmonic oscillator with constant frequency ω 0 , namely a linear combination of the wave functions with H n (x), n = 0, 1, 2, ... the Hermite polynomials.The solutions (2.12) are not energy eigenstates, as the problem has an explicit time dependence.The 'ground state' must be selected through appropriate boundary conditions.We require that the solutions are reduced to the standard wave functions of the harmonic oscillator for τ → −∞, consistently with the selection of the Bunch-Davies vacuum for the field theory.For the 'ground state', we have F 0 (τ, x) = F 0 0 (τ, x), while the function b(τ ) must satisfy b(τ ) → 1 for τ → −∞.At the time τ = 0 corresponding to the transition to the RD era, we impose the continuity of the wave function.This fixes both b(0) and b ′ (0) and determines uniquely the form of b(τ ) for τ > 0. We refer the reader to [35] for the details of this straightforward calculation.We list here the expressions for the function b(τ ) in the dS and RD eras, which determine the 'ground state' for our problem: It can be seen easily that b dS tends to 1 for τ → −∞, while b dS and b RD , as well as their derivatives, are continuous at τ = 0.

The calculation of the entanglement entropy
In this work we calculate the entanglement entropy following the original approach by Srednicki [16].A huge advantage of this approach is that it provides the reduced density matrix as an intermediate result.This matrix contains all the information about entanglement, and the entropy can be derived from its spectrum.It must be kept in mind that the inverse process is not possible, since the entanglement entropy provides only partial information about entanglement.
Srednicki's method treats the field as a collection of coupled harmonic oscillators.When the state of these oscillators is specified, the determination of the density matrix is straightforward.The next step is to trace out some of the degrees of freedom, corresponding to the field values on one side of an entangling surface, in order to derive the reduced density matrix.The difficulty of this task depends on the state in which the overall system lies.The diagonalization of the Hamiltonian decouples this system by describing it in terms of the field canonical modes.In the original calculation [16] the field is considered in flat spacetime and is assumed to be at its ground state, i.e. each decoupled canonical mode is at its ground state.This greatly simplifies the calculation, as the ground state is Gaussian, and tracing out degrees of freedom can be carried out via the evaluation of Gaussian integrals.In our case the field theory is considered in an expanding background.There is an effective time-dependent mass term, which deforms the ground state of canonical modes so that it corresponds to the Bunch-Davis vacuum during the dS phase.The continuity of the wave function uniquely determines the state during the subsequent RD phase as well.The state is Gaussian at all times, as is evident from eq. (2.12), even though the coefficient of the quadratic term in the exponent is complex.This implies that the tracing out of degrees of freedom can be carried out via the evaluation of Gaussian integrals, as in the original calculation.
The state of the overall system reads where the vector x denotes collectively the field values f lm,j , following the notation of [16,35,36], as we explained in the previous subsection.The matrix W is given by (2.17) The matrix Ω is the eigenfrequency matrix that corresponds to the time-independent part of the Hamiltonian of the overall system.It is the positive square root of the matrix of couplings K between the degrees of freedom.In other words, if x i are the coordinates and π i the conjugate momenta, the time-independent part of the Hamiltonian reads As is evident by equation (2.5), the Hamiltonian can be divided into angular momentum sectors that do not interact, allowing for the writing of the matrix K as K = l,m K l , where The matrices K l have positive eigenvalues k li .Each N × N matrix K l has 2 N square roots, that have eigenvalues ± √ k li .As positive square root Ω l of the matrix K l we define its square root that has only positive eigenvalues.Then, the matrix Ω reads Ω = l,m Ω l .In eq.(2.17) we have indicated explicitly the dependence of the function b(τ ) of the previous subsection on the eigenvalues of this matrix.The matrix W is in general a complex symmetric matrix.The overall system is described by the density matrix which is Gaussian.We consider as subsystem 1 the set of n degrees of freedom described by the coordinates x j , where j ≤ n.The rest of the degrees of freedom comprise the complementary subsystem 2. We trace out subsystem 2. It is convenient to introduce the block notation where, in an obvious manner, the matrix A is an n×n matrix, the vector x 1 is an n-dimensional vector, and so on.The matrices A and C are complex symmetric matrices like W , whereas the matrix B has no specific symmetry property and is not even a square matrix.It is a matter of simple algebra with Gaussian integrals to show that the reduced density matrix describing subsystem 1 assumes the form where In order to calculate the entanglement entropy, one needs to specify the eigenvalues of the above reduced density matrix.In order to do so, we must solve for its eigenfunctions, which satisfy There are two differences between our case and the simpler case of a field theory at its ground state in flat space [16].The first one is that the matrix δ does not vanish.This is a very innocent change, as it can be shown that the matrix δ does not alter the eigenvalues of the reduced density matrix, but only its eigenfunctions [35,36].Therefore, for the purpose of the specification of the spectrum of the reduced density matrix, the matrix δ can be set to zero by hand.
The second difference induces several complications.The matrix β is a complex Hermitian matrix, whereas in the original calculation [16] it is real and symmetric.If these properties persisted in our case, the matrices γ and β would be simultaneously diagonalizable via an appropriate coordinate transformation, resulting in a reduced density matrix that could be factored into the tensor product of density matrices describing a single degree of freedom each.Each of these density matrices would describe a simple harmonic oscillator mode at a thermal state, with temperatures given as functions of the eigenvalues of the matrix γ −1 β.In our case this is not possible and the calculation is more complicated.Nevertheless, it can be shown that the general structure of the eigenfunctions remains identical.There is still a 'ground' eigenstate, n 'first excited' eigenstates and so on.However, the reduced density matrix cannot be factored into density matrices of a single degree of freedom each.
The complete analysis of the problem is presented in refs.[35,36], to which the reader is referred for the details.As it turns out, the eigenvalues of the reduced density matrix can be specified in terms of the matrix W that satisfies the quadratic matrix equation where This matrix appears in the exponent of the Gaussian part of the eigenfunctions of the reduced density matrix.For example, the 'ground' Gaussian eigenstate reads Then, the eigenvalues of the reduced density matrix assume the form where m i , i = 1, 2, . . ., n, take any non-negative integer value and the parameters ξ i are the eigenvalues of the matrix Ξ = βT (I + W) −1 . (2.30) Notice that the matrix Ξ is neither real or Hermitian, yet it turns out that it possesses real eigenvalues.The reason is that it is similar to a Hermitian matrix, as shown in section 4.4.2 of ref. [36].
A complication that appears in this calculation is that the matrix equation (2.26) has more than one solutions.Only one of them gives rise to normalizable eigenfunctions of the reduced density matrix.This is the only one which, through eq.(2.30), corresponds to a matrix Ξ that has no eigenvalues larger than 1, and thus to an appropriately normalized spectrum of the reduced density matrix (2.29).The problem of selecting the correct solution of eq.(2.26) can be bypassed by upgrading it to the eigenvalue problem of a matrix of higher dimension.More specifically, it can be shown that the eigenvalues of the matrix Ξ are a subset of the eigenvalues of the matrix These eigenvalues λ are the roots of the characteristic polynomial Because of its structure, the eigenvalues come in pairs, with one element being the inverse of the other.Thus, half of the eigenvalues are smaller than 1 and the rest are larger than 1.The eigenvalues of the matrix Ξ that corresponds to the correct W are exactly the eigenvalues of M that are smaller than 1.This allows the specification of the spectrum of the reduced density matrix in a way that is amenable to a numerical calculation.The entanglement entropy is directly given by the formula The reader is referred to [35,36] for the proof of the above statements and more details.
3 Numerical results

General expectations
The formalism presented in the previous section provides the basis for the numerical calculation of the entanglement entropy associated with a massless free field in an expanding background.As we discussed in the introduction, we assume the presence of an UV cutoff k s for the wavenumbers of the modes that contribute to the entropy.For a discretized theory on a lattice, this cutoff is set by the lattice spacing ϵ, so that k s ∼ 1/ϵ.An IR cutoff k l may also be present, even though we find that it does not have an effect on the entropy.For the discretized theory, the value of k l is set by the total size of the lattice.The structure of the theory implies that the result for the entropy depends on the combination Hϵ, with H the constant value of the Hubble parameter during the dS era, which also sets the initial value of this parameter in the subsequent RD era.Our assumption that the relevant physical degrees of freedom are the ones that exit the horizon during inflation and re-enter during the RD era implies that Hϵ ∼ 1 (for a scale factor set equal to 1 at the dS-RD transition).Before presenting the numerical results, it is instructive to consider the expected form of the entanglement entropy and its dependence on the various scales of the theory.In the dS phase the leading contributions to the entropy, apart from a constant, are expected to be We have used physical parameters and followed the notation of [27], neglecting a term involving the mass of the field.Here A p = 4πa 2 (τ )R 2 is the proper area of the entangling surface and ϵ p = a(τ )ϵ the physical UV cutoff.The terms proportional to c 1 and c 2 are present in a flat background as well.The coefficient c 1 is regularization-scheme dependent and was first The entanglement entropy for a spherical region as a function of the radius of the entangling surface and time for Hϵ = 1.The radius of the total spherical lattice is L = N ϵ, and we depict the results for N = 200 (brown surface), N = 100 (red surface), and N = 50 (green surface).We also indicate the entropy at the dS to RD transition (black curve) and the location of the comoving horizon (dashed, red curve).computed in the seminal work of Srednicki [16], while the universal constant c 2 = −1/90 was computed in [23,24].The term proportional to c 4 is an additional UV divergence appearing in a curved background.The terms proportional to c 5 and c 6 are genuine IR effects related to the expansion.The value c 6 = 1/90 was found in [27].
The numerical accuracy of our calculation allows the precise determination of only the leading area terms.The logarithmic terms proportional to c 2 and c 6 cannot be reliably specified within the accuracy of our numerical calculations.The logarithmic dependence of the term proportional to c 4 affects the coefficient of the area law and is tractable.We postpone its precise analysis for a future publication.It must be pointed out that no volume term is expected to develop in the dS phase, despite the fact that the modes lie in squeezed states.The technical reasons in the context of our approach are explained in section 5 of [35].
The form of the entanglement entropy in the RD phase is more complicated.Several curvature invariants may contribute, which are reduced to the terms involving H in the dS phase.Similarly to before, an area term is present.However, as we shall see, a new feature appears: a volume term develops and becomes the dominant contribution to the entropy.A particular question, which we address in subsection 3.3, is whether this term arises from UV contributions, or it is a low-energy effect.

The evolution of the entropy
In fig. 1 we depict the entanglement entropy for a spherical region as a function of the radius of the entangling surface and time for Hϵ = 1.The entropy and radius axes are logarithmic.
The discretization of the theory and the implementation of the UV cutoff were described in subsection 2.1.The radius of the total spherical lattice is L = N ϵ, and we depict the results for N = 200 (brown surface), N = 100 (red surface), N = 50 (green surface).We have displaced slightly the three surfaces in the vertical direction in order to make them more visible.The radius of the entangling surface is R = nϵ, with n = 1, ..., 25, so that it is always smaller that half of the radius of the overall lattice and the interior region contains manifestly fewer degrees of freedom than the exterior.The transition from the dS to the RD era takes place at τ = 0, in a region characterized by a strong increase of the entanglement entropy.The thick black line depicts the entropy as a function of n at τ = 0. We have also indicated the location of the comoving horizon at various times through a dashed, red line.For τ < 0 in the dS era, the comoving horizon shrinks and a diminishing part of the total system remains subhorizon.For τ > 0 in the RD era, the comoving horizon grows and the subhorizon part of the system grows.At the dS to RD transition only a small part, with size of the order of the lattice spacing, is subhorizon.
We observe that, for τ → −∞, the entropy becomes independent of time.In the logarithmic plot, ln S is a linear function of ln n with slope equal to 2 to a very good approximation, apart from expected deviations for small n, resulting from subleading contributions to the entropy (constant and logarithmic terms).We have ignored such corrections, as the precision of our calculation is not sufficient to determine them reliably, and we have focused on large values of n for which terms of quadratic or of higher order dominate.For τ → −∞, we have performed a fit of the entropy with a quadratic function in the region n ≥ 15 in order to determine the coefficient s of the quadratic term.We obtain s ≃ 0.09 for all three values of N .This is to be compared with the value s ≃ 0.3 quoted by Srednicki [16].The difference arises from the different regularization scheme that we use, as we discussed in detail in subsection 2.1.In order to confirm the reliability of our numerics, we have also computed the entropy using the regularization of [16], reproducing the value s ≃ 0.3.The form of the entropy is consistent with the expectation that it should coincide with the entropy in a static background when the physical radius of the entangling surface is much smaller than the Hubble radius.This is the case for τ → −∞, when the scale factor approaches zero and the term proportional to c 1 in eq.(3.1) dominates.This result is also consistent with the assumption of a Bunch-Davies vacuum.
When the time approaches τ = 0, there is a strong increase of the entanglement entropy while the background is still in the dS phase.The dependence on n remains strictly quadratic until the transition to the RD era.This feature is expected for a free theory, as was discussed in [35].The increase of the entropy arises because the terms proportional to c 4 and c 5 in eq.(3.1) become dominant.In particular, we have verified that the entropy grows ∼ a 2 (τ ), following the growth of the physical radius of the entangling surface when the comoving radius is kept fixed, as in our calculation.The transition to the RD phase is followed by a further increase of the entropy, associated with the additional squeezing of the wave function of the canonical modes.For τ → +∞ the entanglement entropy develops an almost constant form, with a profile that indicates deviations from a purely quadratic dependence on n.We investigate this form in more detail below.What is apparent in fig. 1 is a weak oscillatory behaviour that depends on N .In particular, the wavelength of the oscillatory pattern is comparable to that of the longest mode allowed by the finite lattice.This pattern is a finitesize effect that would be absent in an infinite system.It is also apparent that the oscillations decay with time, with the resulting asymptotic form of the entropy being independent of N .We have analyzed numerically the behaviour of the entropy at late times, but we have not identified a significant residual effect associated with an IR cutoff k l ∼ 1/L.If such an effect exists, it is not visible with the numerical precision of our calculation.
The next issue that we would like to analyze is the dependence of the entanglement entropy on the Hubble scale H during inflation.As we have already mentioned, the relevant dimensionless parameter is the product Hϵ.The value of s we quoted corresponds to the coefficient of the term n 2 = R 2 /ϵ 2 = (aR) 2 /(aϵ) 2 that involves the ratio of the physical radius to the physical UV cutoff.In fig. 2 we depict the form of the entropy as a function of the radius and time for Hϵ = 0.1 (green surface), Hϵ = 1 (red surface) and Hϵ = 10 (brown surface).The size of the lattice is N = 100.
During the dS era, if the UV cutoff ϵ is much shorter than the Hubble radius 1/H, as in the case Hϵ = 0.1, the entropy is dominated by the short-distance entanglement between degrees of freedom on either side of the entangling surface.The expansion of the background induces only a subleading effect even when the physical entangling surface is larger than the horizon.During the subsequent RD era, the physical cutoff shrinks even further relatively to the Hubble radius, so that the entropy continues to be dominated by the UV, similarly to the situation in static space.
For larger values of Hϵ during inflation, the asymptotic entanglement entropy of the discretized theory at late times becomes larger than that for a static background.As we discussed in the introduction, we focus on the case Hϵ = 1, for which modes that did not cross the horizon during inflation and never froze are not taken into account in the calculation of the entropy.For this choice of UV cutoff, the enhancement of the entropy due to the expansion of the background can be significant relatively to the static-space case, as can be seen in fig.2, in which the vertical axis is logarithmic.Notice that the physical UV cutoff today, which is of the order of 1m as we discussed in the introduction, is much shorter than the Hubble radius today.However, the effect of the squeezing of the wave function during inflation and the transition to the RD era results in an increase of the entropy by at least two orders of magnitude relatively to what one would expect for a static background.Since the determination of the shortest wavelength that froze at the end of inflation is rather imprecise, we have also considered the case Hϵ = 10, for which the enhancement of the entropy relatively to the static case is higher by an additional two orders of magnitude.

The volume term
Apart from the enhancement of the entanglement entropy, a very interesting feature that is apparent in figs. 1 and 2 is the dependence of the entropy on the radius of the entangling surface at late times.As long as the background is in the dS phase, the entropy is dominated by a quadratic term.However, after the transition to the RD era, a different pattern emerges.
In fig. 3 we depict the entanglement entropy as a function of n = R/ϵ at an early time τ /ϵ = −100 in the dS era and a late time τ /ϵ = +200 in the RD era.The size of the lattice is L = N ϵ with N = 100.We fit the shape of the curve for n ≥ 15 with a polynomial that includes a quadratic and a cubic term.The assumed form of the entropy is where we have indicated explicitly the use of the UV cutoff in the parameterization.The coefficients s and c are in general functions of τ /ϵ.Their values at specific times are given within the inset.At early times c vanishes, while at late times a good fit is obtained only if it takes a nonvanishing value.The red, dashed curve and the blue, dotted one correspond to fits with only a quadratic or cubic term, respectively.It is clear that at late times they both are unsatisfactory.
In fig. 4 we depict the coefficents s (left plot) and c (right plot) of the quadratic and the cubic term, respectively, as a function of time.At early times we have s ≃ 0.09 and c = 0.While the background is still in the dS phase, s grows, whereas c remains zero.The cubic term emerges only after the time τ = 0 of the dS to RD transition.The oscillatory patterns arise because of the finite-size effect that we discussed earlier in this section, which is associated with the longest mode that can exist within the finite lattice.The oscillations decay with time and the two parameters asymptotically settle to constant values: s ≃ 2.5 and c ≃ 0.13.The growth of the cubic term is a very interesting, novel feature.It indicates that strong entanglement is not confined only to narrow regions on either side of the entangling surface, but spreads throughout the whole bulk of the system.We discuss the implications of this result in the following section.
An important question is whether the appearance of the volume term is an effect related to the UV modes of the system.The parameterization (3.2) gives the impression that the dependence on the UV cutoff is strong.However, adopting the logic that led to eq. (3.1) suggests a parameterization of the form with R p = a(τ )R the physical radius.Apart from the effects related to the expansion, the area term receives contributions from the UV range.However, this is not necessarily true for the cubic term.The explicit UV dependence has now been shifted into the definition of c.However, this could be counterbalanced by the dependence of c on Hϵ.In fig. 5 we depict c and c/(Hϵ) 3  for Hϵ = 0.1 (blue curve), Hϵ = 0.5 (orange curve), and Hϵ = 1 (green curve).Notice that the size of the lattice is L = N ϵ with N = 200, i.e. twice that in fig. 4.This explains the larger period of the oscillations associated with the longest mode of the system (a finite-size effect discussed in the previous subsection).The slightly negative values of c within a short time interval are caused by the limitations of the numerically calculation and the fit to the results.
The important conclusion that can be drawn from fig. 5 is that the UV dependence of c is such that c is largely independent of ϵ.Even though c varies by 3 orders of magnitude for the values of Hϵ we considered, the average value of c/(Hϵ) 3 varies by less than a factor of 2. The other important point is that the volume term is determined by the comoving radius R = R p /a(τ ).When this term dominates, the entanglement entropy within a spherical region that follows the expansion is roughly constant.
The qualitative picture that emerges is that the entanglement entropy within a spherical region of fixed comoving radius R grows during the dS era, even for a constant comoving UV cutoff ϵ.This growth is caused by the squeezing of the wave function of the various field modes by the expansion, and is additive to the standard entanglement entropy in a static background.However, the dependence on the radius is quadratic, so that the effect satisfies an area law.The transition to the RD phase induces the emergence of a volume term in the entropy.The most plausible interpretation is that the entanglement is spread by the expansion over larger distances, until it encompasses the degrees of freedom of the whole system.The volume contribution to the entropy is roughly constant for a fixed comoving radius.In this sense, the expansion during the RD phase does not generate additional entropy, but redistributes over the whole system the entropy produced during the dS phase.

Conclusions
We have followed the evolution of the entanglement entropy of a real, massless, scalar field throughout the inflationary period and the subsequent era of radiation domination.We have assumed that during inflation the field is in the Bunch-Davies vacuum.As a result, the entanglement of the short-distance modes is very similar to that at the ground state in a static background.For longer modes, the transition towards a squeezed state upon horizon exit of each mode during inflation and the additional squeezing when radiation domination sets in enhance the entanglement entropy.Even though we did not analyse it explicitly, a similar behaviour is expected during matter domination, as has been verified in a previous study [35].
The entanglement entropy, even in a static background, is an UV-divergent quantity that requires the introduction of an UV cutoff.We have identified this cutoff with the wavelength λ s of the last mode that exited the horizon at the end of inflation and re-entered immediately in the RD era.Modes with shorter wavelengths did not go through the process of freezing and the dominance of the growing term in the mode function.We have adopted the view that they can be considered as vacuum fluctuations even today, while we focused on the entanglement due to modes that are directly accessible to observations.A question of particular interest concerns the dependence of the entropy on the size of the system.The characteristic pattern in a static (3 + 1)-dimensional background is that the entropy is proportional to the area of the entangling surface, similarly to the black hole entropy [16].This feature is attributed to the strong entanglement of short-distance modes on the two sides of the entangling surface.One would expect that the rapid expansion during the dS era would stretch the entangled modes sufficiently far from the entangling surface, so that a volume effect might appear.However, this expectation is not realized.The subtle technical reasons are explained in section 5 of [35].On the other hand, the expansion during the RD era induces a drastic modification.A volume term develops and becomes the dominant contribution to the entropy at late times.It must be noted that the presence of a volume term is not unexpected for systems lying in squeezed states, such as those developing through the expansion of the background [36].It is rather a peculiarity of the cosmological evolution that the volume term appears only after the universe has evolved into the RD era.
The behaviour that we outlined above must be contrasted with the standard picture of quantum to classical transition upon horizon exit for the cosmological fluctuations produced during inflation.In our analysis we have fully accounted for the quantum properties of the field, without discarding the decaying mode, no matter how smaller it may become than the growing one.In this sense, the entropy that we have computed should be attributed to the quantum entanglement of the degrees of freedom.For the short-distance modes it coincides with the entanglement entropy in a static background.
At late times during the RD era, the volume contribution to the entanglement entropy within a fixed comoving radius approaches a constant value.A similar behaviour is expected during matter domination.The possibility that the quantum properties of a massless field can be traced even today seems very exciting.Even though we studied a free scalar field, our analysis could be relevant for all very weakly interacting fields as, for example, gravitational waves resulting from tensor modes during inflation, under the assumption that they evolve almost freely and their quantum coherence is not lost during the whole evolution of the universe until today through some secondary process.Finding appropriate experimental signatures of the quantum origin of such fields seems a difficult task [7][8][9].However, the question is of fundamental importance [66,67].
Even though the picture outlined above seems the most natural, the final form of the entropy, and in particular the appearance of a volume term, points towards a possibly different interpretation.It is known that the effective decoherence underlying the quantum to classical transition is achieved only if the decaying mode is neglected.On the other hand, it has also been pointed out that even an exponentially small decaying mode may be important for the calculation of the entropy of the perturbations [57].Moreover, there is no clear criterion that can characterize the system as classical when the decaying mode is not neglected and the state remains pure [68][69][70].It is then possible that the smallness of the decaying term results in a system with significant entropy that must be interpreted in classical terms.Our result emerges through the non-trivial form of the density matrix induced by a rapidly expanding background that contains a horizon.A possible interpretation would attribute thermal characteristics to the reduced density matrix resulting from tracing out degrees of freedom beyond a certain radius, which are classically inaccessible to an observer.The regions beyond the horizon provide the typical example for such a situation.The presence of the volume term would then support the interpretation of the entropy of the observable universe as thermal, even though its origin would lie in the presence of the horizon.
Let us recall the simple example of two coupled harmonic oscillators lying at their ground state [16].The reduced density matrix that describes either of the two is the density matrix of a single harmonic oscillator 1 , lying at a mixed state, which is thermal with a temperature that depends on the coupling between the two oscillators.An observer with knowledge of the existence of both oscillators concludes that the form of this reduced density matrix results from the entanglement between them and the corresponding entropy is entanglement entropy.However, an observer who has access only to one of the two oscillators cannot reach such a conclusion.According to this observer, the only way to interpret the entropy is as thermal.
In our study we calculated the entanglement entropy in the expanding universe splitting the system into two through a spherical entangling surface.A realization of this surface is provided by the cosmological horizon, which physically prohibits an observer from measuring the observables associated to the degrees of freedom in the exterior.A hypothetical observer with knowledge of the entire universe would attribute the entropy to the entanglement between the interior and the exterior.However, an observer confined within the horizon can only perceive the entropy as thermal.The presence of the volume term that we deduced implies that this constrained observer does not just see entropy emanating from the cosmological horizon, similarly to a black hole, but concludes that the bulk of the interior has been heated.This conclusion is reached despite the fact that the whole system still lies at a pure state.
One must bear in mind that in our case the reduced density matrix does not correspond to a system in exact thermal equilibrium.When a subsystem of a harmonic system contains many degrees of freedom, the reduced density matrix is organized through effective canonical modes, each of which lies at a different temperature.When squeezing is present this picture is not exact, since these canonical modes cannot be associated with a real combination of the local degrees of freedom.However, the appearance of many temperatures (as many as the degrees of freedom of the subsystem) is still valid.This fact is a consequence of the integrability of our toy model, which is a free field theory.The thermalization hypothesis suggests that in a realistic theory the reduced density matrix would actually be thermal.
We note that there is no distinct point in the evolution that we observed which can be associated with a qualitative change in the nature of the entropy.The evolution is smooth, as can be seen in figs. 1 and 2. The structure of the reduced density matrix is complicated and does not provide any hints for a qualitative transition.On the other hand, the thermodynamic interpretation is appealing because it is consistent with the quantum to classical transition.The quantum effects are tied to the decaying mode, which can be many orders of magnitude smaller than the growing one for realistic scenarios of inflation.A thermal entropy is consistent with the picture of a stochastic classical field, as a result of the dominance of the growing mode.
The thermodynamic interpretation could provide a quantum-mechanical realization of reheating after inflation.In our analysis the background is introduced by hand and not as a solution of the Einstein equations for an appropriate equation of state.In this sense, the RD era does not arise as a result of the field properties.On the other hand, it is intriguing that the appearance of a volume term, a characteristic feature of thermal entropy, is connected to the transition to the RD era.
It is usually assumed that the reheating occurs through the strong interactions between the products of the inflaton decay.Interpreting our result as thermal entropy would imply that all fields (interacting or not) that exist at the end of inflation can be treated in the same way, since they would all become part of the thermal environment during reheating.Of course, there is no indication of a common temperature for all the constituents of the universe.However, our understanding of the situation is still very incomplete.
Our numerical results indicate that the magnitude of the entropy can be significant.If we estimate it through the volume term that develops during the era of radiation domination, we get a value for the observable universe ∼ (Hλ s ) −3 , with a cutoff λ s ∼ 1m, as we discussed in the introduction, while the current value of the Hubble radius is 1/H ∼ 10 26 m.This gives a value ∼ 10 78 for the entropy, which is to be compared with the standard thermal entropy ∼ 10 88 associated with the plasma in the early universe, transferred to the photons and neutrinos today.There is a discrepancy of 10 orders of magnitude, which can be understood if we recall that the typical wavelength of the cosmic background photons is ∼ 10 −3 m.The origin of the difference lies in that the scale of the standard thermal entropy is connected to the energy density of the inflaton background, while the scale of the entropy we computed for a massless field is set by the last mode that exited the horizon at the end of inflation.The two would be comparable if the energy scale of inflation was close to the Planck scale.However, the observations indicate an energy scale roughly three orders of magnitude lower.Despite the quantitative discrepancy, the interpretation of the entropy as thermal does not seem completely off the mark.
The basic conclusion that can be reached through our analysis is that the rapid expansion of the background during inflation and the transition to the RD era generate a significant amount of entropy even for a system of a non-interacting field.This entropy is carried by the long-distance modes that go through the process of horizon exit during inflation and are accessible to experiments today.It can be defined and computed as entanglement entropy, and is intrinsically linked to the presence of a horizon.It seems unlikely that the possible loss of quantum coherence, which we did not address in this work, will result in the suppression of the entropy.It seems more likely that the reduced density matrix, employed by an observer with no access to the regions beyond the horizon, will remain nontrivial, with a thermodynamic interpretation.This speculation is consistent with the standard picture of the quantum to classical transition during inflation.More work is needed in order to obtain a deeper understanding of these issues.The notion that the entropy of the universe can be attributed to the presence of the cosmological horizon merits further exploration.

Figure 2 .
Figure 2. The entropy as a function of the radius of the entangling surface and time for Hϵ = 0.1 (green surface), Hϵ = 1 (red surface) and Hϵ = 10 (brown surface).The size of the lattice is N = 100.

Figure 3 .
Figure 3.The entanglement entropy as a function of n = R/ϵ at an early time τ /ϵ = −100 in the dS era and a late time τ /ϵ = +200 in the RD era.The size of the lattice is L = N ϵ with N = 100.The red, dashed curve and the blue, dotted one correspond to fits with only a quadratic or cubic term.

Figure 4 . 5 .
Figure 4.The coefficients s (left plot) and c (right plot) of the quadratic and the cubic term, respectively, as a function of time, for the theory of fig. 3.