Entanglement Entropy of Cosmological Perturbations

We show that the entropy of cosmological perturbations originating as quantum vacuum fluctuations in the very early universe, including the contribution of the leading nonlinear interactions, can be viewed as momentum space entanglement entropy between sub- and super-Hubble modes. The interactions between these modes causes decoherence of the super-Hubble fluctuations which, in turn, leads to a non-vanishing entropy of the reduced density matrix corresponding to the super-Hubble inhomogeneities. In particular, applying this to inflationary cosmology reveals that the entanglement entropy produced by leading order nonlinearities dominates over that coming from the squeezing of the vacuum state unless inflation lasts for a very short period. Furthermore, demanding that this entanglement entropy be smaller than the thermal entropy at the beginning of the radiation phase of standard cosmology leads to an upper bound on the duration of inflation which is similar to what is obtained from the Trans-Planckian Censorship Conjecture.


I. INTRODUCTION
There has recently been a lot of interest in entanglement entropy in the context of quantum field theory and gravity (see e.g., [1] for reviews). In particular, the entanglement entropy of a conformal field theory is holographically related to properties of the bulk in the context of the AdS/CFT (anti-de-Sitter bulk/conformal field theory on the boundary [2]) correspondence (see e.g., [3]). In the same context, entanglement entropy can be related to properties of black holes in the AdS bulk [4]. The relationship between the bulk Einstein equations and properties of entanglement of the boundary CFT was explored in [5]. Entanglement entropy considerations have also been applied directly to black holes physics (see [6] for a review), and to de Sitter space in [7,8]. There are also attempts to build up space-time itself from quantum entanglement [9].
Most considerations of entanglement are based on a position space separation of the domain; for example, the separation between the inside of a black hole and the outside. However, in cosmology it is more natural to work in momentum space because it is the properties of the momentum modes of cosmological fluctuations which are generally probed (such as the power spectrum). Momentum space entanglement has been considered in [10] (see also [11]), and we will use methods from that work extensively.
Entanglement is a crucial, and rather essential, feature of quantum mechanical systems. In many early universe scenarios, the cosmological fluctuations which we measure today are postulated to emerge from quantum vacuum perturbations. This is the case not only in inflationary cosmology [12], but also in the Ekpyrotic sce-nario [13] and in the matter bounce scenario [14]. Cosmological perturbations (see e.g., [15,16] for reviews) are small amplitude fluctuations about the homogeneous and isotropic cosmological background. Because of their small amplitude, the inhomogeneities are generally described in Fourier space. To leading order, each Fourier mode evolves independently, and each mode obeys a harmonic oscillator equation with a time-dependent mass. The Hubble radius H −1 (t) (where H is the Hubble expansion rate) plays a key role in the dynamics of the modes: on sub-Hubble scales the canonical fluctuation variable oscillates, while it is squeezed on super-Hubble scales.
Successful early universe scenarios have the common feature that the fluctuation modes which are probed today in cosmological observations were sub-Hubble in the early universe phase, thus allowing a causal generation mechanism. In the classes of models we consider here, the initial state for the fluctuations is taken to be the quantum vacuum state 1 . When the fluctuation modes exit the Hubble radius, their state becomes a squeezed vacuum state. The Hilbert space of states thus naturally divides into two parts -the super-Hubble mode space H A (t) and the sub-Hubble mode space H B (t): where H k is the harmonic oscillator Hilbert space of the k'th mode and H −1 c (t) stands for the comoving Hubble radius. It is natural to consider the space of super-Hubble modes to be the system we consider, and the space of sub-Hubble modes to be the bath which we integrate over. Note that the comoving Hubble radius decreases as a function of time in the early universe phase of the models which we consider. This means that modes exit the Hubble radius. Hence, the boundary between the two Hilbert spaces H A and H B depends on time: the dimension of the system Hilbert space is increasing. This is a specific feature of a system on a dynamically expanding background. Furthermore, although not explicitly stated above, we shall assume an ultraviolet (UV) cutoff (M Pl ) for the bath modes so that there is always a constant supply of modes which we integrate over. We assume that some underlying UV theory is able to provide the details of the dynamics of the modes lying in the range k > M Pl and shall not consider them in our work.
As mentioned above, in this paper, we consider the entropy of the space of super-Hubble modes which results from the entanglement with the bath of sub-Hubble modes. The question of entropy of cosmological perturbations has been considered previously. For example, in [18][19][20] the entropy of a classical field was studied, and the results were applied to compute the entropy of cosmological perturbations and gravitational waves in an inflationary universe. In [18][19][20], the source of entropy can be traced back to the loss of information about the phases of the fluctuations for super-Hubble modes, while a similar calculation for the coherent state basis was shown in [21]. In [22], the issue of entropy of cosmological perturbations was reconsidered, taking the loss of information which leads to entropy generation to be the loss of information due to the spreading of the wave function of the super-Hubble modes which results from squeezing. Entropy generation as a consequence of coupling to an environment was studied in [23]. In [24], entropy generation of cosmological fluctuations as a consequence of a truncation of the hierarchy of Green's functions was considered.
What was not considered in these previous works on entropy generation is the role of nonlinearities. Because of the nonlinear nature of the Einstein equations, there is always a mixing of modes for cosmological perturbations. In particular, there is a mixing between the suband super-Hubble modes. As discussed in [25][26][27][28], this leads to decoherence of the reduced density matrix of super-Hubble modes 2 . This decoherence is crucial in order to explain why the cosmological perturbations become classical even though they have a quantum origin. The resulting density matrix of the super-Hubble modes is no longer that of a pure state, and hence leads to a non-vanishing entropy which we compute in this paper. We stress that, as shall become apparent later on, we calculate a lower bound on the amount of entanglement entropy of scalar density perturbations, produced in any model of inflation, due to the minimal gravitational nonlinearities which must always be present. Additional couplings or fields, or considering interactions between scalar and tensor modes, would lead to enhanced amounts of entropy production.
There are some similarities between our work and that of [33], where decoherence through neglecting observationally inaccessible correlators was considered, and that of [34] where decoherence via entropy field loops was studied (decoherence of fluctuations through entropy loops was considered earlier in [35]). There is also a connection with the work of [36] where super-Hubble entanglement through inflaton decay was considered.
Our notation is as follows: We use natural units in which the speed of light, Planck's constant and Boltzmann's constant are set to one. We consider a spatially flat background cosmology such that the metric can be written as where η is the conformal time which is related to the physical time t via dt = adη, and x are the comoving spatial coordinates. The Hubble parameter is given in terms of the scale factor a(t) by where the overdot represents the derivative with respect to t. We emphasize that the Hubble radius plays a crucial role in our analysis. Sub-Hubble modes of the canonical fluctuation variable oscillate while those on super-Hubble scales are squeezed [15,16]. We denote the Planck mass by M Pl . In the next section, we give a first pass at arriving at the entropy of cosmological perturbations due to the squeezing of super-Hubble modes during inflation. In Sec-III, we review the well-known argument that interaction between the perturbation modes, arising from minimal gravitational nonlinearities, leads to a suppression of the off-diagonal terms in the density matrix for the super-Hubble modes. This justifies an assumption used in Sec-II for calculating the entropy due to the squeezed state. Finally, having set up our dominant interaction term in Sec-III, we go on to calculate the entanglement entropy density for our system (super-Hubble) modes in Sec-IV. We estimate an order of magnitude for the upper bound of this quantity and show that it is greater than the entropy for the squeezed vacuum, as calculated in Sec-II. In Sec-V, interestingly we find an upper bound on the duration of inflation by requiring that this entanglement entropy remains smaller than the thermal entropy produced at the end of inflation 3 . We discuss our main findings in Sec-VI.

A. The Squeezed Vacuum
We consider linear scalar cosmological perturbations about the background metric (2). Assuming that the matter source of the fluctuations has no anisotropic stress, the perturbations are described by a single field ζ(x, t), the curvature perturbation in comoving gauge. The metric including these fluctuations is The action for cosmological perturbations has a canonical kinetic term if we use the rescaled field (we are following the notation of [39]) with where H is the first "slow-roll" parameter defined via and c 2 s is the speed of sound squared of the matter source. Although, later on, we shall only consider models of single-field inflation with no derivative self-couplings, we are keeping c s = 1 at this stage so that our expressions remain as general as possible 4 .
The linear cosmological perturbations about the classical background geometry can be canonically quantized [12]. We insert the ansatz for the fluctuating metric and matter into the total action (joint gravitational and matter action) and expand to quadratic order. Since at linear order each Fourier mode evolves independently, we can reduce the quantization to the standard quantization of a set of harmonic oscillators, the oscillators having a time dependent mass coming from the time dependence of the background. In terms of the usual ladder operators, the quadratic Hamiltonian H 2 corresponding to scalar cosmological perturbations takes the form This quadratic Hamiltonian generates the following equation of motion for the ladder operators Given an initial condition at an instant of time, η 0 , we can solve for this as In the above, r k and φ k are the squeezing parameter and the squeezing angle, whereas θ k denotes the action of the rotation operator. The number of particles in a given mode k is proportional to the squeezing parameter n k ∼ sinh 2 r k . For inflation, the leading order time-dependence of these parameters is given by [40] Given the quadratic Hamiltonian, the evolution operator U 0 (η) can be written as where S k (r k , φ k ) and R k (θ k ) are the two-mode squeezing and rotation operators, respectively, which are defined as [40] S k := exp At the level of the quadratic Hamiltonian, the U 0 (η) is unitary. However, once interaction terms are introduced, the evolution becomes necessarily non-unitary in the presence of bath modes [39]. The effect of the rotation operator is only to change the phase and would be of no consequence to us, and hence we drop it from hereon. The effect of the two-mode squeezing operator on the vacuum leads to the squeezed vacuum, which is defined as For a given mode k, it is easy to see that this state is normalized, as follows: as required. The squeezed vacuum of all the modes can be obtained in a straightforward manner as the tensor product state

B. The Reduced Density Matrix
The straightforward definition of the density matrix, corresponding to the squeezed state given in (20), is If we calculate the entropy corresponding to this state, naturally this is going to be zero since it is a pure state, given by the evolution of the vacuum under the quadratic Hamiltonian (8). More concretely, the density matrix expressed in terms of the two-mode occupation number basis reads which is still a pure density matrix. Let us show this more explicitly, as follows. Our state can be written as a product state where |ψ A is the product state of all the super-Hubble modes, and |ψ B over the sub-Hubble modes. Since we are focusing on the super-Hubble modes, our reduced density matrix is obtained by tracing over the the sub-Hubble mode Hilbert space.
where the sum is over the basis states of the Hilbert space of sub-Hubble modes. In the absence of entanglement between the sub-and super-Hubble modes, and given that the states of both subsystems are pure, the reduced density matrix ρ A also corresponds to that of a pure state and hence has vanishing entropy. So far, however, we have neglected any coarse graining or nonlinear effects. In particular, we have neglected entanglement effects between sub-and super-Hubble modes which are inevitably present because the equations of gravity are nonlinear. In the following we will take a first look at the entropy of cosmological perturbations after loss of some information about the state. In the following section we then show that this loss of information is an inevitable consequence of the entanglement between sub-and super-Hubble modes.

C. First View on Entanglement Entropy of Cosmological Pertubations
In order to get a non-vanishing von-Neumann entropy of the reduced density matrix ρ A , we need to coarse-grain it in a suitable way to derive a mixed density matrix. In [18,19], it was observed that the phase associated with the squeezing angle is sensitively dependent on the density perturbation whereas the amplitude is not. As a consequence, the coarse-grained entropy in [18,19] was defined by averaging over the squeezing angle, which also leads to decoherence. In our setup, a similar "averaging" over the squeezing angles would lead to setting the offdiagonal elements to zero in the number basis, leading to a reduced density matrix of the form A different perspective of arriving at the above form for the reduced density matrix would be to consider only the diagonal entries of (22), whereas assuming that the off-diagonal elements quickly fall-off to zero. The usefulness of this perspective lies in the fact that one does not have to refer to the phase in order to derive the reduced density matrix. However, now we need to justify our choice of ignoring the off-diagonal elements for the density matrix. One way to argue would be to consider that there are a lot of particles created for a given mode, with opposite momenta, with φ k being the phase of each of these particle pairs. But if we want to use the destruc-tive interference, while group averaging these phases, as being responsible for suppressing the off-diagonal terms, then we are back to our previous argument. Instead, one might follow the arguments of [22,41] to justify the reduction of the density matrix as a result of assuming a distribution of coherent states as our initial state -instead of the usual vacuum -as a manifestation of our ignorance regarding initial conditions. If one assumes this as the starting point, it can be shown that the offdiagonal terms are naturally suppressed as long as one invokes equipartition of probabilities for the initial states in the ensemble [41]. We are neither advising this approach nor suggesting that it is better than considering the averaging procedure over random phases, but just pointing out that there have been different justifications for considering the above form of the reduced density matrix (25). In the next section, we will give an improved analysis and explain the decay of the off-diagonal elements as a consequence of decoherence resulting from entanglement between the modes.
The von-Neumann entropy associated with this reduced density matrix is given by First, we expand the product in the logarithm as a sum of logs, i.e.
Using this in (26), we can rewrite the entropy density (per comoving volume) as Using the normalization (19), the term in the first parentheses is equal to 1. The entropy gets simplified to In the large occupation number limit, n k = sinh 2 r k 1, we get back the same expression for the entropy density s ≈ k ln sinh 2 r k , as derived in [18,19]. However, we derived this result from the von-Neumann entropy formula for a quantum density matrix instead of using the Shannon entropy for a classical field. Note that one should expect that our expression matches that for the classical calculation, done earlier, only in the large squeezing limit. In this sense, one should view k ln sinh 2 r k as the classical limit of the von-Neumann entropy calculated here within a quantum field theoretic approach, and it is thus compatible with previous results [18,19] of considering the entropy of a classical field. Our result also matches with previous works as presented in [22].
In the case of slow-roll inflation with an approximately constant Hubble constant we can estimate the resulting entropy density by integrating over all super-Hubble modes and apply a infrared cutoff: we do not consider modes with wavelengths larger than the Hubble radius H −1 at the beginning of inflation. With the convention that the scale factor is set to one at the beginning of inflation, this implies that in (29) we need to integrate over all values of k with H < k < aH. At any time, this integral is dominated by the modes exiting the Hubble radius at that time, and we thus obtain 5 To obtain the entropy density per physical volume element, we have to divide the above by a 3 , and we hence get Before moving on, let us note that the entropy calculated in this section is not quite an entanglement entropy as it arises from the squeezing of the cosmological perturbations. The way we manage to get a nonzero result for a density matrix arising from a quadratic Hamiltonian (8) is by employing some yet-to-be-specified coarse-graining, due to which the pure density matrix in (22) is reduced to a mixed one (25), by ignoring the off-diagonal terms. In the next section, we shall give a rigorous argument as to how gravitational nonlinearities, responsible for decohering the quantum fluctuations into classical perturbations, necessarily render the density matrix diagonal. In this way, the entanglement between sub-and super-Hubble modes, due to mode-mixing arising from gravitational non-linearities, is also responsible for the entropy of cosmological perturbations calculated above 6 .

III. NONLINEARITIES, DECOHERENCE AND ENTROPY GENERATION
Here we review the analysis of [28] which shows how the purely gravitational interactions which are inevitably present because of the nonlinearity of General Relativity lead to a decoherence of the reduced density matrix of the super-Hubble modes as a consequence of the interaction with the sub-Hubble fluctuations. For our purposes, we will focus on the case of inflation.
We shall now take into account the effects of the cubic Hamiltonian in addition to the quadratic Hamiltonian discussed in the previous section. This is the leading term which generates entanglement between the sub-and super-Hubble modes. We are considering the full cubic action for the density perturbations in the presence of a single matter field. If the matter is a canonically normalized scalar field, then the speed of sound c s = 1. In more general models, c 2 s can be smaller than one, and this can significantly increase the size of the cubic interaction terms, resulting in a significant contribution to the equilateral-shape non-Gaussianity parameter f N L . However, as a first pass, let us only consider vanilla matter models with c s = 1, which should be sufficient to estimate a lower bound on the entanglement entropy for models of inflation.
We take the form of the cubic contribution to the Hamiltonian from [42], from now on setting c s = 1, which is a generalization of the results from [43].
We have also introduced the second "slow-roll" parameter: We shall ignore the non-local terms that containχ since those are not the dominant terms in the action. Additionally, there are also terms which would get cancelled with each other (such as theη H term in the second line would get cancelled by a similar term from the third line of (32)). Since, in the case of inflation, the dominant mode of ζ has frozen out on super-Hubble scales, we will neglect interaction terms which containζ. Furthermore, we shall restrict our analyses only to the leading order terms in the slow-roll parameters, and would thus be left with the second term in the first line of (32) (the other terms being higher orders in H and η H , or contain ȧ ζ). Hence, the dominant term in the interaction Hamiltonian is (after integration by parts, and recalling that The H int we are considering arises purely from gravitational non-linearities, originating from the cubic Lagrangian given in (32). As discussed above, in a model of single-field slow roll inflation without any derivative selfinteraction, this would be the dominant term. However, for a nontrivial speed of sound model, there can a different term which significantly enhances the cubic interaction. This would lead to both a faster rate of decoherence as well as a greater amount of entanglement entropy. In this sense, our calculation should be understood to yield the minimum amount of entanglement entropy that must be produced in any inflationary model; multiple fields or more complicated interactions would only enhance our results. Note here that there is an additional term, not shown above in the cubic Lagrangian, that is of the exact same form, ζ 2 (∂ 2 ζ), but with a pre-factor H η H a [44]. This term is part of a large number of terms which are typically removed by a field redefinition [43] and do not affect the correlation functions for calculating the bispectrum.
Strictly speaking, we should keep this term if we are interested in calculating the entropy corresponding to the ζ field (and not for the redefined one). However, we drop it here to avoid additional clutter since it is straightforward to include its effects at the end by adding a factor of H η H , in addition to the 2 H in (34), to our results. Having setup our interaction terms, we begin the evolution at the conformal time η 0 , in a pure Gaussian product state of all of the modes, which has the wave function where in this case we have indicated which variables the individual states depend on. As a consequence of the interactions, the state evolves into at a later time η, where the third factor is a consequence of the interaction Lagrangian.
The interaction contribution to the wave function is given by where k, k stand for sub-Hubble modes, and q stands for a super-Hubble mode, and the kernel function F(k, k , q; η) is given by an integration over time of the interaction Hamiltonian in momentum space (see [28] for details) with the property that its imaginary part blows up as η → 0. In the above, the integration runs over all momenta with the property that k + k + q = 0 (momentum conservation).
The reduced density matrix of the super-Hubble modes can be obtained by integrating over the sub-Hubble ones. In the field representation we have where D B stands for the integration over the sub-Hubble modes B. Eq. (38) yields where D[ζ,ζ] is the decoherence factor. Focusing on a single super-Hubble mode q, the decoherence factor is where the time dependence of the factors has been suppressed, where P G is a property of the Gaussian wavefunction, and As is clear from (40), the decoherence factor decays in time on super-Hubble scales since the imaginary part of F blows up. Note that the decoherence effect is dominated by the Hubble scale modes. There is no UV divergence in the loop diagram which produces the interaction. This is a consequence of the specific form of our interaction Lagrangian.
To conclude this section, we have reviewed how the interaction with the sub-Hubble modes leads to decoherence of the super-Hubble ones. For a particular mode, decoherence happens after Hubble radius crossing. The important thing for us is the fact that decoherence leads to the damping of the off-diagonal terms such that the reduced density matrix of the super-Hubble modes become diagonal very quickly. In the previous section, we had calculated the entropy corresponding to the squeezing of the super-Hubble modes, assuming that their density matrix turns diagonal which we have given a concrete justification for in this section. As promised, we show that even for the entropy of the cosmological perturbations which solely arises from the quadratic Hamiltonian, nonlinearities play a crucial role by making the density matrix diagonal. In the following section, we will compute the entanglement entropy which the non-Gaussianities directly generate.

A. Setup
Having set up our interaction terms, let us discuss how one can calculate the entanglement entropy of the cosmological perturbations due to the effects of these coupling terms. To calculate the entanglement entropy, we shall follow the prescription of [10], and generalize their results for flat spacetime to inflationary backgrounds.
Given our breakup of the Hilbert space (1) H = H A ⊗H B into system and environment modes, our Hamiltonian can be expressed as where H A,B denote the free part of the Hamiltonian and λ is a time-dependent constant. The ground state of the free theory, neglecting the interactions, is denoted by |0, 0 = |0 ⊗ |0 , and one can write the interacting vacuum of the entangled system as where |n denotes an n-particle state of the system (in fact, a product state over all super-Hubble k modes), and |N is the corresponding state for the bath.
Following the analyses of [10], one finds that the leading order contribution to the entanglement entropy for such a system can be written as where we can express the matrix element C n,N in terms of standard perturbation theory as Note that our definition of C n,N differs slightly from that of [10] for later convenience. Before going on to calculate this matrix element, and the corresponding entanglement entropy for our cosmological system, let us review the flat space calculation first through an explicit example.

B. Calculation for flat space
Considering a cubic interaction term, one can write the action for a massive scalar field as For a flat (3 + 1)-dimensional spacetime, the field can be decomposed in terms of the usual ladder operators as where ω k = √ m 2 + k 2 . Here, instead of putting the fields in a box as in [10], we choose to work with continuous field variables, as would be more appropriate for cosmological perturbations later on. However, we still have a scale µ which separates our system from the environment, using the same convention as in [10]. In other words, we are interested in calculating the entanglement entropy between the modes with momenta k above and below µ. In this case, the only nontrivial contribution to the matrix element would be from an excited state of a 3-particle one which can be written as Recalling that the interaction Hamiltonian is (λ/3!) ϕ 3 , λ having dimension of mass, the required matrix element (45) can be written as In the second line above, we only keep the creation operators as required, whereas in the third line we have used the orthonormality property of the inner product to eliminate the integrals over (k 1 , k 2 , k 3 ). In the final step, we used the integration over the spatial coordinate, and the remaining delta function implies that at least one of the spatial momenta must be above, and at least one below, the scale demarcating the system and the environment. The entanglement entropy for this system can be then evaluated by plugging in the above expression into (44) s flat ent = −λ 2 log λ 2 where the integrals are over a set of momenta such that there can only be two configurations of interest -either one of (p 1 , p 2 , p 3 ) is greater than µ while the rest are below µ, or two of them are above while one is below µ. We have also divided the total entanglement entropy by the (infinite) volume to express it as an entanglement entropy density (≡ S flat ent /Vol).

C. Vacuum & Interaction Hamiltonian
Let us first outline the differences we anticipate between the flat space calculation above and our case for cosmological perturbations. Firstly, the interaction parameter λ = λ(η) will now be time-dependent. Secondly, the vacuum for the system modes is now given by the squeezed vacuum, and the mode functions corresponding to the vacuum in curved spacetime will have a different form of their momentum dependence. Since the vacuum of the super-Hubble modes will now be the squeezed vacuum, there now are contributions of terms with both creation and annihilation operators in our case. Finally, a major conceptual difference arises from the fact that the scale separating our system from the bath is given by the (comoving) Hubble scale which is time-dependent since we are working with comoving coordinates (and, in addition, by itself has a weak time-dependence of its physical value during inflation), and is not some arbitrary, tunable parameter µ as in the flat space case. With this in mind, let us begin by factoring the Hamiltonian for the overall system as where the H sys and H bath is the quadratic Hamiltonian, for the super and sub-Hubble modes respectively, as given in (8). Next, we write down the vacuum modes for the unperturbed systems, ignoring nonlinearities, as The |0, 0 is the vacuum state for both the system as well as the bath modes. For the super-Hubble modes, the vacuum is given by the squeezed state as given in (20). On the other hand, we have the usual Minkowski vacuum for the sub-Hubble modes, denoted by |0 . The explicit form of the interaction Hamiltonian naturally depends on the choice of the interaction term we choose between the perturbation modes. As mentioned earlier, for this paper, we shall restrict ourselves to only cubic perturbation terms which arise naturally from gravitational nonlinearities in any model of inflation, as captured by our interaction Lagrangian given in (32). We emphasize once again that considering more complicated interactions or more fields can lead in a different term dominating H int , which would end up producing enhanced amounts of entanglement entropy. In this precise sense, we give a lower bound on the amount of entropy production coming from scalar modes during inflation.
For our dominant interaction term of the form we can write down the interaction Hamiltonian by converting the ζ field to our canonical field χ, and then expanding in terms of the creation and annihilation operators in momentum space. We find the following expression [45]: where all the terms in the parentheses (. . . ) are the same and include all possible (momentum-conserving) combinations of the ladder operators. We have also defined . The difference in the momenta dependence of our choice of H int from, say, one with time-derivatives such as L 3 ∼ ζ(ζ ) 2 , would be that some of the terms in the expression above would come with a minus sign since, in that case, the interaction term couples the field with its conjugate momentum [39]. The prefactor is given by (keeping in mind that we go from cosmic time to conformal time) where, as anticipated, we get a time-dependent interaction parameter. We now have all the ingredients -the vacuum state and the interaction Hamiltonian -to calculate the matrix element given in (45).

D. Matrix element
Let us revisit our calculation of the matrix element for the cubic Lagrangian in Minkowski space. The crucial difference between that calculation and the one for inflation would be that instead of only keeping the term which solely involves creation operators from the interacting Hamiltonian, we shall also have to consider terms of the form c k1 c † −k2 c † −k3 and c k1 c k2 c † −k3 . This is easy to understand since for the case of flat spacetime, the only nonzero contribution for the matrix element between the Minkowski vacuum and an excited state (with, say, three particles for a cubic interaction) can come if we sandwich a term consisting of three creation operators in between. If there exists any annihilation operator, it would simply annihilate the vacuum, resulting in zero. On the other hand, for inflation, we have a tensor product of the Minkowski vacuum for the sub-Hubble modes and the squeezed vacuum for thr super-Hubble ones (52). In this case, the ladder operator(s) corresponding to the sub-Hubble modes must be creation ones c † −k whereas the one(s) corresponding to the super-Hubble modes can be either c † −k or c k . This is so because an annihilation operator c k does not annihilate the squeezed vacuum |SQ(k, η) . One can see this explicitly from the form of the two-mode squeezed vacuum, as given in (17).
Having said this, let us list all the possible choices of interaction terms which can appear in the matrix elements: • Terms of the form c † −k c † −k c † −k : There can be either two system (super-Hubble) modes and one bath (sub-Hubble) mode or vice-versa.
• Terms of the form c k c † −k c † −k : There can be either two system modes and one bath mode or vice-versa. However, the annihilation operator must always correspond to the super-Hubble mode.
• Terms of the form c k c k c † −k : There must be two system modes, corresponding to the two annihilation operators, and can, therefore, only be one bath mode.
• The terms proportional to c k c k c k necessarily yield zero for the matrix element since the annihilation operator corresponding to any of the bath modes annihilates the Minkowski vacuum.
Let us consider the first case in detail in the following calculation while we leave the details of the other terms for the Appendix. Therefore, the term of interest for us from the H int (54), for calculating (45), is the following: Next, we need to find the explicit action of a creation operator on the squeezed vacuum. Using the definition of the two-mode squeezed state from (17), we can formally express the action of a creation operator on it as Schematically, it implies that we are considering an excited state with a particle of energy p over our squeezed vacuum. A similar iteration would create higher order excited states over the squeezed vacuum. However, recall that for a cubic interaction term, the only non-zero contribution to the matrix element comes from having the first excited state over both the squeezed and the Minkowski vacuum. Also, since we are only considering cubic interactions, there can be only two choices -either one of the modes is in the system and two are in the bath or two of them are in the system while one is in the bath. However, it will be clear from the following that the dominant contribution to the entanglement entropy comes from having two of the modes in the bath and one in the system. This is not at all surprising keeping in mind that the decoherence rate is also dominated by having two short-wavelength modes and one long-wavelength one. Let us consider the former option first, i.e. p 1 , p 2 > aH while p 3 < aH. The appropriate excited state to consider is of the form The only other novelty for our calculation is the effect of the squeezed vacuum on the inner product. Recall the standard result where we have written things schematically to avoid clutter. To explicitly see how this result comes about, one should write down the unitary transformation of the creation and annihilation operator under the squeezing operator, i.e. S † cS and S † c † S as linear combinations of c, c † , dropping all momenta indices. Also, note that S † = S −1 . See the Appendix for more details. The rest of the calculation follows exactly that of flat space, and it is easy to evaluate the matrix element as (c † c † c † ) C sq n,N = (2π) 3 1 + sinh 2 r p3 It is clear that for our choice of p 1 , p 2 ∈ bath while p 3 ∈ system, the dominant term in the above comes from the third term C n,N ∝ p1p2 p3 . It is also evident from the above calculation that if we had two modes in the system and one in the bath, then the dominant term in the matrix element would have the form (c † c † c † ) C fold n,N = (2π) 3 1 + sinh 2 r p2 1 + sinh 2 r p3 where we have chosen p 1 > aH and p 2 , p 3 < aH. Already at this stage we can see that the entanglement entropy for cosmological perturbations, during inflation, peaks in the "squeezed" limit p 3 p 1 ≈ p 2 , given the momentum structure of the matrix element, for (c † c † c † ) C sq n,N whereas it gets its maximum contribution in the "folded" limit p 3 + p 2 ≈ p 1 for the other case (c † c † c † ) C fold n,N .

E. Entanglement entropy
Let us recall the formula for the leading order term in the entanglement entropy where a sum is implied on both types of C n,N calculated in (59) and (60). Note our slight difference in convention of defining the matrix element C n,N as in (45), with that of [10]. In order to calculate the entanglement entropy, we need to reinstate factors of the coupling parameter λ(η) = √ H / 2 √ 2 a(η) M Pl as well as the energy corresponding to the ground and excited states, both for the Minkowski and the squeezed vacuum, considered above. However, what we really need in the above formula is the energy difference between the first excited state and the ground state, for both the Minkowski and the squeezed vacua. This is the same for both the system and the bath modes and is given by ω k := k for (nearly) massless scalar excitations.
Note that the sum over (n, N ) translates into integrals over all the momentum modes in the formula (44). Recall that there was a similar integral over all momentum modes also in the expression of the entropy arising from the squeezing part of the quadratic Hamiltonian, as shown in (29). However, unlike in that case, we would have the integrals over all momentum conserving configurations involving (p 1 , p 2 , p 3 ) and not over individual modes as is expected for an entanglement entropy coming from cubic interactions. Keeping this is mind, the entanglement entropy (per unit comoving volume) is given by where we have only kept the dominant terms from the matrix elements (59) and (60). It is important to discuss the limits of the above integral first: We have introduced M Pl as the natural physical UV cutoff and the comoving wavenumber at the beginning of inflation as the infrared cutoff. We set a i = 1 for the scale factor at the beginning of inflation (and therefore, in our convention, a is always > 1). We also assume that the Hubble parameter, H, remains constant during inflation. Furthermore, the UV cutoff for the comoving momenta is given by aM Pl which signifies the fact that the integration of the environment is over a fixed number of bath modes, even though we are considering an accelerating background. This is so because although the environment is continuously depleted by modes getting redshifted into the system, there is also a constant supply of modes from the UV into the bath 7 . However, the system has an increasing phase space of modes as more and more modes become super-Hubble as time goes on, and given our infrared cutoff which states that there were no comoving modes which were super-Hubble before inflation started. Naturally, we have to assume that inflation starts at a finite time in the past which reinforces the need of having an UV cutoff for the perturbation modes.
Let us now estimate the integrals I 1 and I 2 given in (62). For I 1 , when we have two bath modes and one system mode, the integrand would naturally have its largest contribution coming from the squeezed limit, as shown below: In the second line, we have killed the p 1 integral using the delta function, introducing the angle Θ between p 2 and p 3 . In the next line, we introduce the crucial approximation that the integrand peaks in the limit Θ → π/2 and p 2 p 3 , i.e. the squeezed limit. This would help us in getting an upper bound on the entanglement entropy corresponding to the I 1 term. We have also used the expression for the squeezing parameter from (11) and used the approximation that 1 + sinh r k ≈ sinh r k , for large squeezing, in this step. It is then easy to see that the integration over the bath modes is dominated by the upper limit (the UV cutoff scale), while the integral over the system mode p 3 is dominated by the lowest value of p 3 , i.e. by the infrared (IR) cutoff scale. We have only kept the leading terms in the integrals in the same spirit to arrive at our lower estimate for the entropy density, ignoring numerical factors. We note that a factor of a 3 should be divided from the final result in order to account for the entanglement entropy density (total entropy per unit physical volume). We are then left with a factor of (a/a i ) 2 (recall, we have set a i = 1) and this reflects the fact that the phase space of the system modes is growing, and the contribution to the p 3 integral is dominated by the IR cutoff. Collecting everything, the estimate 8 of the entanglement entropy per unit physical volume coming from I 1 is given by where a > 1 is such that the number of e-foldings of inflation is given by N := ln a in our convention.
Let us now first show that the contribution coming from I 2 to the entanglement entropy density would be subdominant to the above result. In this case of having two system and one bath mode, the largest contribution to the integrand would come from the folded limit p 1 ≈ p 2 + p 3 . Following the calculation as in the previous case, we can arrive at an upper bound for the estimate of this term in a similar way. However, note that once we eliminate the integral over the bath mode p 1 using the delta function, none of the system mode integrals which are left have any dependence on M Pl . The other difference lies in the additional squeezing terms leading to an extra factor of the IR cutoff in the final result, namely, Once again, we have expressed this final result in terms of the entanglement entropy per unit physical volume and have only given a rough estimate of the upper bound.
Thus, we find which means that s I2 ent shall always remain subdominant to s I1 ent , provided f > 1 ⇒ N < 3 ln(M Pl /H). In the next section, we shall show that combining the observed scalar power spectrum with the fact that the entanglement entropy of cosmological perturbations during inflation remain smaller than the thermal entropy produced during (p)reheating leads to this condition being always satisfied. Therefore, we can always ignore the entanglement entropy corresponding to having two system and one bath mode when compared to that of having two sub-and one super-Hubble mode.
Note that the above estimates were calculated using the approximations of squeezed and folded shapes, in which the integrands reach their peak values. The full integrals do not lend themselves to having simple analytic forms and we have thus avoided writing them down explicitly. The effect of removing these approximations would result in some small numerical factors appearing in front of our estimates, as in (64). However, recall that we have only shown here the result of the calculation of the entanglement entropy coming from the terms of the form c † −k c † −k c † −k , arising from the interaction Hamiltonian in (54). As mentioned earlier, there are other terms, proportional to c k c † −k c † −k and c k c k c † −k , which also contribute to the entanglement entropy. As shown in the Appendix, in the limit of large squeezing, r k 1, the contribution of all of these terms are either proportional to s I1 ent or to s I2 ent . Naturally, we neglect the terms proportional to s I2 ent since they are sub-dominant. And the terms which are proportional to s I1 ent shall add to our estimate for the entanglement entropy density (64). All of this is to say that in our order of magnitude estimate for the entanglement entropy density of cosmological perturbations during inflation, there should be some O(1) numerical factor appearing, namely There are two sources which contribute to this O(1) number -one from the additional terms, as shown in the Appendix, and the other coming from the fact that we are estimating the integral by its upper bound. From now on, we shall drop this number as well as the logarithmic factor in our upcoming discussions. Now that we have an estimate for the entanglement entropy due to the gravitational nonlinearities, let us compare this with the contribution coming from the squeezing part of the quadratic Hamiltonian, as in (29). As mentioned earlier, for large r k 1, the entropy density (per physical volume), coming from (29), is given by (31) where we have, once again, ignored some small numerical factors. Although s ent corresponding to cubic interactions arising from gravitational nonlinearities is suppressed by a factor of H (as it should be), it is still greater than s sq . One way to easily see this is to approximate the value of the observed scalar power spectrum as As we shall see from the bounds on N that we will derive in the next section, this quantity t > 1 and thus the entanglement entropy from non-Gaussianities would be larger than that corresponding to the squeezed vacuum, provided inflation lasts a reasonable amount of time and is not fine-tuned to be extremely small. This is quite a remarkable result since this implies that the entanglement entropy due to (cubic) gravitational nonlinearities are larger than that due to the (squeezing part of the) quadratic action!

V. UPPER BOUND ON THE DURATION OF INFLATION
We have seen that the entanglement entropy density of cosmological perturbations produced by nonlinearities builds up during a period of inflation as where N is the number of e-foldings of inflation, and a i is the value of the scale factor at the beginning of inflation (which we had set equal to 1 in the last section, for simplicity). In order to allow a graceful exit from inflation consistent with the second law of thermodynamics, it is important to make sure that the entropy due to these interactions remain subdominant to the entropy in the thermal radiation state after inflation. This thermal entropy density is given by where T R is the initial temperature of the radiation bath, and g * is the number of spin degrees of freedom in the radiation bath. Assuming rapid thermalization after inflation, and nearly constant Hubble parameter during inflation, this yields which results from the TCC [38]. Note that this bound on the duration of the inflationary phase is the same as derived in [46], where it was argued that beyond that time the de Sitter phase cannot be given a well-defined classical background interpretation due to the buildup of entanglement 9 . We are thus led to speculate the the TCC may have a derivation based on entropy considerations and the second law of thermodynamics. It is already known that entropy considerations have also proven useful [49] to derive the de Sitter swampland conjecture [50,51], one of the various constraints on effective field theories to be consistent with string theory (see e.g. [52,53] for reviews).
Note that we have derived a lower bound on the entanglement entropy due to the minimal gravitational nonlinearities (ignoring those due to tensor perturbations). We might speculate that if we were to do a more detailed calculation, our entropy bound on N might turn out to be in even closer agreement with the bound from the TCC. Note that the bound (76) can be relaxed if we consider H to be decreasing substantially during inflation, or if the thermal history of the universe after inflation is nonstandard. However, as shown in [54,55], in these cases the TCC bound is also relaxed. Note, also, that if we take into account entanglement entropy due to modes which were already super-Hubble at the beginning of inflation, the bound can be strengthened, in the same way that the TCC bound is strengthened if we consider pre-inflation evolution [55,56]. Finally, it has also been pointed out that deriving the TCC from different quantum gravity arguments can, by itself, lead to a refinement of it [57] and can bring it closer to our bound.
Returning to the discussion at the end of the previous section concerning the ratio of the entropies produced by nonlinear entanglement effects on one hand, and by pure decoherence of the linear modes on the other, we see that if the duration of inflation saturates the above bound (76), then the entanglement entropy dominates by a factor of (M Pl /H) 3/2 , the result we promised to derive earlier. In other words, unless inflation lasts for a very short period of time, s ent would always dominate over s sq .
Note that a related bound on the duration of inflation based on entanglement considerations was given in [58], where it was argued that, interpreting the current horizon entropy of the Universe as entanglement entropy, there is a number of e-foldings of inflation before which there is no entropy and we cannot talk about a de Sitter background.

VI. CONCLUSIONS AND DISCUSSION
In this work, we have derived the entanglement entropy of inflationary scalar perturbations, corresponding to nonlinearities arising from gravity. Although entropy of cosmological perturbations is a rich subject by itself, what is novel to our work is that we calculate the entanglement entropy to the leading order of cubic interactions, going beyond the calculation of entropy corresponding to the squeezing of the super-Hubble vacuum state. Remarkably, we show that this cubic (and higher order) interactions are essential even to calculate the entropy corresponding to the quadratic Hamiltonian. This is so because decoherence arising from these terms is what is responsible for reducing the pure density matrix to a mixed one, by suppressing the off-diagonal terms. These higher order interaction Hamiltonians themselves lead to mode-couplings such that there is an entanglement between the super-and sub-Hubble modes which is a direct manifestation of the quantum origin of these vacuum fluctuations 10 .The entanglement entropy corresponding to these interactions is what we have calculated for the first time by treating the super-Hubble modes as our system and the sub-Hubble ones as a bath.
Our result shows that the entanglement entropy density scales as H 2 M Pl (a/a i ) 2 , where a i is the scale factor at the beginning of inflation. In order to allow for a graceful exit from inflation consistent with the second law of thermodynamics, this entropy must be smaller than the thermal entropy after inflation. This leads to an upper bound on the duration of inflation which is very close to the bound obtained from the TCC. Interestingly, the nonlinearities produce the dominant contribution to the entropy of cosmological perturbations, surpassing the one for the squeezed vacuum, provided > (H/M Pl ) (a i /a) done in momentum space. It is easy to appreciate this properly if one compares our result with that for determining the full non-unitary evolution of the density matrix of the system modes as has been done, for instance, in [39] (see [45] for the case of tensor modes). The time evolution of the reduced density matrix involves non-Hamiltonian terms, and might even contain non-Markovian terms, which depend on the so-called Lindblad operator. If one were to try and calculate the solution of the time-dependent reduced density matrix and then evaluate the von Neumann entropy associated with it, the calculation would become much harder and rather intractable. In this paper, we give a complementary way of calculating the entanglement entropy without having to deal with the full dynamics since, as emphasized earlier, we only require to calculate certain matrix elements for our purposes. The fact that these two seemingly different methods yield the same result for the entanglement entropy has been shown in [59] for any quantum field theory. In addition, going to momentum space makes it easy to impose a UV cutoff for the bath modes, as has been done in this case.
The natural next step for us would be to calculate the entanglement entropy corresponding to primordial gravitational waves. Once again, assuming the simplest model of inflation, nonlinearities would arise from gravitational interactions which would lead to decoherence and entropy production. Therefore, this calculation would also give an improved lower bound on the amount of entropy which must be produced in any model of inflation. Furthermore, the leading interactions between the tensor perturbations are not slow-roll suppressed which typically lead them to decohere faster than their scalar counterpart [45]. Anticipating along similar lines, we expect that the entanglement entropy of tensor modes would be some-what enhanced, and this will be studied in future work. The cubic interactions coupling tensor and scalar modes also need to be taken into account which will result in enhancing both the entanglement entropy density of the scalar as well as the tensor perturbations.
Finally, we note that our analysis has been done in the context of inflationary cosmology, but the methods also apply to other early universe scenarios in which the primordial fluctuations are quantum in origin, in particular to the matter bounce and to the Ekpyrotic scenarios.