Superhorizon entanglement entropy from particle decay in inflation

In inflationary cosmology all particle states decay as a consequence of the lack of kinematic thresholds. The decay of an initial single particle state yields an \emph{entangled quantum state of the product particles}. We generalize and extend a manifestly unitary field theoretical method to obtain the time evolution of the quantum state. We consider the decay of a light scalar field with mass $M\ll H$ with a cubic coupling in de Sitter space-time. Radiative corrections feature an infrared enhancement manifest as poles in $\Delta=M^2/3H^2$ and we obtain the quantum state in an expansion in $\Delta$. To leading order in $\Delta$ the pure state density matrix describing the decay of a particle with sub-horizon wavevector is dominated by the emission of superhorizon quanta, describing \emph{entanglement between superhorizon and subhorizon fluctuations and correlations across the horizon}. Tracing over the superhorizon degrees of freedom yields a mixed state density matrix from which we obtain the entanglement entropy. Asymptotically this entropy grows with the \emph{physical} volume as a consequence of more modes of the decay products crossing the Hubble radius. A generalization to localized wave packets is provided. The cascade decay of single particle states into many particle states is discussed. We conjecture on \emph{possible} impact of these results on non-gaussianity and on the ``low multipole anomalies'' of the CMB.


I. INTRODUCTION
Quantum fluctuations during inflation seed the inhomogeneities which are manifest as anisotropies in the cosmic microwave background and primordial gravitational waves. In its simplest inception the inflationary stage can be effectively described as a quasi-deSitter space time. Early studies [1][2][3][4][5][6] revealed that de Sitter space time features infrared instabilities and profuse particle production in interacting field theories. During inflation the rapid cosmological expansion modifies the energy-uncertainty relation allowing "virtual" excitations to persist longer, leading to remarkable phenomena, which is stronger in de Sitter space time as clarified in ref. [7]. Particle production in a de Sitter background has been argued to provide a dynamical"screening" mechanism that leads to relaxation of the cosmological constant [8][9][10] through back reaction, much like the production of particle-antiparticle pairs in a constant electric field. More recently this mechanism of profuse particle production has been argued to lead to the instability of de Sitter space time [11,12].
A particular aspect of the rapid cosmological expansion is the lack of a global time-like killing vector which leads to remarkable physical effects in de Sitter space time, as it implies the lack of kinematic thresholds (a direct consequence of energy-momentum conservation) and the decay of fields even in their own quanta [13][14][15] with the concomitant particle production. This result that was confirmed in ref. [16][17][18] and more recently in ref. [19] by a thorough analysis of the S-matrix in global de Sitter space.
The decay of an initial single particle state into many particle states results in a quantum state that is kinematically entangled in momentum space: consider the example of a scalar field theory with cubic self-interaction and an initial single particle state with spatial physical momentum k, namely |1 k , this state decays into a two-particle states of the form p C p (t) |1 p |1 k− p where C p (t) is the time dependent amplitude of the two particle state with momenta p and k − p respectively. This is an entangled state that features non-trivial correlations between the product particles. In ref. [13,15] it is argued that in de Sitter space time with Hubble constant H, the largest decay amplitude corresponds to the case when one of the product particles features physical momenta p ≪ H, therefore, if the initial particle has physical momenta k ≫ H and one of the product particles features a momentum p ≪ H (the other with | k − p| ≫ H) the quantum entangled state features correlations between the sub and superHubble daughter particles.
We refer to these correlated pairs produced from the decay of a parent particle as entangled across the Hubble radius, namely "superhorizon" entanglement, referring to the Hubble radius in de Sitter space time as the horizon as is customary in inflationary cosmology.
Correlations of quantum fluctuations during a de Sitter inflationary stage have been recently argued [20] to lead to remarkable Hanbury-Brown-Twiss interference phenomena with potential observational consequences.
Unitary time evolution of an initial single particle state is a pure quantum state in which the product particles are kinematically entangled.
If a pure quantum state describes an entangled state of several subsystems and if the degrees of freedom of one of the subsystems are not observed, tracing the pure state density matrix over these unobserved degrees of freedom leads to a mixed state reduced density matrix. The entanglement entropy is the Von Neumann entropy associated with this reduced density matrix; it reflects the loss of information that was originally present in the quantum correlations of the entangled state.
The main purpose of this article is to study the entanglement entropy in the case of an initial quantum state describing a single particle state with physical momentum k ≫ H decaying into a pair of particles one with p ≪ H (superhorizon), and the other with | k − p| ≫ H (subhorizon) by tracing over the super-Hubble ("superhorizon") degrees of freedom. This entanglement entropy is a measure of the loss of information contained in the pair correlations of the daughter particles.
The entanglement entropy has been the focus of several studies in condensed matter systems [21][22][23][24], statistical physics and quantum field theory [25][26][27][28][29], black hole physics [30][31][32] and in particle production in time dependent backgrounds [33]. Most of these studies focus on entanglement between spatially correlated regions across boundaries. The entanglement entropy in de Sitter space-time for a free, minimally coupled massive scalar field has been studied in ref. [34] with the goal of understanding superhorizon correlations, and ref. [35] studied the entropy from momentum space entanglement and renormalization in an interacting quantum field theory in Minkowski space-time.
Our study differs from these studies in many ways: we are not considering spatially correlated regions, and momentum space entanglement resulting from the kinematics of particle decay in states of the same quanta is different from the cases studied in ref. [35] which considered momentum space entanglement in the interacting ground state of two coupled theories or a finite density case, both in a stationary, equilibrium situation, whereas we are interested in the time evolution of the reduced density matrix and the concomitant increase of the entanglement entropy in an interacting theory in de Sitter space time.
More recently the entanglement entropy in the ubiquitous case of particle decay in Minkowski space-time from tracing over the degrees of freedom of an unobserved daughter particle has been studied in ref. [36] as a characterization of an "invisible" decay complementary to missing energy.
We focus on light scalar fields with mass M 2 ≪ H 2 , for which radiative corrections feature infrared divergences that are manifested as poles in ∆ = M 2 /3H 2 ≪ 1 [13,15] in the self-energy leading to a consistent expansion in ∆. A similar expansion was recognized in refs. [37][38][39][40].
The field theoretic method introduced in ref. [15,36,41] that describes the nonperturbative time evolution of quantum states is extended here and then generalized to inflationary cosmology (for other applications of this field theoretical method see refs. [42]) to obtain the entangled quantum state from single particle decay to leading order in a ∆ expansion. We show explicitly that unitarity is manifest in the time evolution of the quantum state. From this state we construct the (pure) density matrix and trace over the contribution from superhorizon modes and obtain the entanglement entropy to leading order in a ∆ expansion. Whereas in ref. [43] the entanglement between only two modes was studied in de Sitter space time, ours is a full quantum field theoretical treatment that includes coupling between all modes as befits a local quantum field theory and consistently trace over all the superhorizon degrees of freedom.
We find that the entanglement entropy asymptotically grows with the physical volume as more wavevectors cross the Hubble radius. The method is generalized to a wave packet description of single particle states and we study in detail the case of wave packets sharply localized in momentum around a wavevector k 0 ≫ H and localized in space on scales much smaller than the Hubble radius all throughout the near de Sitter inflationary state. We find that under these conditions, the entanglement entropy for wavepackets is approximately the same as that for plane waves and assess the corrections.
As mentioned above, the lack of kinematic thresholds implies that quanta can decay on many quanta of the same field, in particular for cubic interactions a single particle state can decay into two particles of the same field, however the decay process does not stop at the two particle level, but instead is a cascade decay 1 → 2 → 3 → · · · . We provide a non-perturbative framework to study this cascade decay process and argue that for weak (cubic) coupling λ there is a hierarchy of time scales and the cascade is controlled by this weak coupling. The probability of multiparticle states is suppressed by λ 2 for each extra particle in the final state, the time scales of production and decay of multiparticle states are also separated by 1/λ 2 . We comment on possible relationship with non-gaussianity, in particular pointing out the relationship between the quantum correlations between subhorizon and superhorizon quanta from particle decay and the bispectrum of scalar perturbations in the squeezed (local) limit. Furthermore, we speculate as to whether the information "lost" as modes cross the horizon is "recovered" when the modes re-enter the horizon during the matter dominated era. This study then bridges the main concepts of entanglement between spatial regions explored in ref. [34], with momentum space entanglement and coarse graining [35] and quantum entanglement via particle decay [36] in inflationary cosmology.

II. QUANTUM FIELD THEORETICAL WIGNER-WEISSKOPF TREATMENT OF THE DECAY WIDTH
The method developed in refs. [15,36,41,42] is a quantum field theoretical generalization of the Wigner-Weisskopf method used in quantum optics [44,45].
We consider a scalar field minimally coupled to gravity in a spatially flat de Sitter spacetime with scale factor a(t) = e Ht . In comoving coordinates, the action is given by It is convenient to pass to conformal time η = −e −Ht /H with dη = dt/a(t) and introduce a conformal rescaling of the fields a(t)φ( x, t) = χ( x, η). (II. 2) The action becomes (after discarding surface terms that will not change the equations of motion) with primes denoting derivatives with respect to conformal time η and where for de Sitter spacetime where H (1,2) ν (z) are Hankel functions. Expanding the field operator in this basis yields The Bunch-Davies vacuum is defined such that a k |0 = 0 , (II. 12) and the Fock space states are obtained in the usual manner, i.e. by applying creation operators a † k to the vacuum. In what follows we consider a light scalar field with M ≪ H and write For light scalar fields with ∆ ≪ 1 quantum loop corrections feature an infrared enhancement from the emission and absorption of superhorizon quanta that is manifest as poles in ∆ [13,15]. Below we exploit the expansion in ∆ implemented in ref. [13,15] to leading order, isolating the most infrared sensitive contributions to the entanglement entropy from these processes.
In the Schröedinger picture the quantum states |Ψ(η) obey where in an expanding cosmology the Hamiltonian H(η) is generally a function of η in marked contrast to the situation in Minkowski space-time, where it is constant. Introducing the time evolution operator U(η, η 0 ) obeying the solution of the Schröedinger equation is |Ψ(η) = U(η, η 0 ) |Ψ(η 0 ) . Now separate out the interaction Hamiltonian by writing H(η) = H 0 (η) + H i (η) with H 0 (η) the non-interacting Hamiltonian, and introduce the time evolution operator of the free theory U 0 (η, η 0 ) satisfying 16) the interaction picture states are defined as Here U I (η, η 0 ) is the time evolution operator in the interaction picture and obeys and where χ is the free field Heisenberg field operator in eq.(II.11).

A. Transition amplitudes and probability
Now consider a cubic interaction Hamiltonian for a scalar field which we label as χ( x, η) after the conformal rescaling described above: We can then use the expansion of the scalar field χ given by (II.11) to compute the transition amplitude for a one particle state to decay into two particles χ k → χ p + χ k− p as depicted in fig. (1): The total transition probability is where q = | k − p|. Note that this kernel has the property that Introducing the identity 1 = Θ(η 2 − η 1 ) + Θ(η 1 − η 2 ) in the (conformal) time integrals and using (II.24) we find from which we obtain the transition rate as In Minkowski space-time (η → t), if the kinematics of the transition is allowed, i.e. if energy-momentum conservation obtains, the transition is to on-shell states and the transition probability grows linearly in time, exhibiting secular growth. In the long time limit the transition rate becomes a constant. This is basically how the result from Fermi's Golden rule comes about. If, on the other hand energy-momentum conservation is not fulfilled, the probability becomes constant at asymptotically long times, with a vanishing transition rate, describing virtual processes that contribute to wave function renormalization. A true decay of the quantum state is therefore reflected in secular growth of the transition probability and a transition rate that either remains constant or grows at asymptotically long time. In de Sitter space time the lack of a global time-like Killing vector implies the lack of kinematic thresholds. As discussed earlier in ref. [13,15,41] and confirmed in ref. [19], quanta of a single field can decay into other quanta of the same field regardless of the mass of the field. In this subsection, we review the work in refs. [15,36,41,42], as the implementation of the quantum field theoretical Wigner-Weisskopf formulation is crucial in constructing states whose time evolution is manifestly unitary.
Expanding the interaction picture state |Ψ(η) I in Fock states |n obtained as usual by applying the creation operators on to the (bare) vacuum state as which in terms of the coefficients C n (η) become it is convenient to separate the diagonal matrix elements from those that represent transitions, writing Although this equation is exact, it provides an infinite hierarchy of simultaneous equations when the Hilbert space of states |n is infinite dimensional. The Wigner-Weisskopf method consists of two main ingredients: i) truncation of the hierarchy at a given order in the perturbative expansion, ii) a Markovian approximation that yields the long time asymptotics of the coefficients. In ref. [41] the equivalence between the Wigner-Weisskopf method, the time evolution obtained from the Dyson resummation of propagators in terms of the self-energy and the dynamical renormalization group was shown in Minkowski space time. Hence this method provides a non-perturbative resummation to obtain the real time dynamics of quantum states.
We begin by implementing this program to lowest order, and provide a roadmap for implementation at arbitrary higher order in section (VI) where we also study "cascade processes" that are available in de Sitter space time.
Thus, consider the case when a state |A , say, couples to a set of states |κ , which in turn couple back to |A via H I . Then to lowest order in the interaction, the system of equation closes in the form where the κ =A is over all the intermediate states coupled to |A via H I representing transitions. By including the diagonal terms n|H I (η)|n C n specifically, we can also consider mass counterterms [15], however, we will neglect these terms in the sequel since we are not concerned with either mass generation or renormalization in this article.
Consider the initial value problem in which at time η = η 0 the state of the system is given by |Ψ(η = η 0 ) = |A so that (II.33) We can then solve (II.32) and substitute the solution back into (II.31) to find This integro-differential equation with memory yields a non-perturbative solution for the time evolution of the amplitudes and probabilities. We can construct an approximation scheme to solve this equation as follows. First note that the time evolution of C A (η) as determined by eqn. (II.35) is slow in the sense that the relevant time scale is determined by a weak coupling kernel Σ. This allows us to introduce a Markovian approximation in terms of an expansion in derivatives of C A as follows: define Integrating by parts in eq.(II.35) we obtain The first term has "erased" the memory in the kernel by setting both time arguments to be the time of interest, while the second term on the right hand side is formally of fourth order in H I . Integrating by parts successively as discussed in ref. [41] a systematic approximation scheme can be developed. To leading order in the coupling (second order in H I ), we will neglect the second term on the right hand side of (II.39), in which case eqn. (II.35) becomes When the state A is a single particle state, radiative corrections to the mass are extracted from E A and is identified as a (conformal) time dependent decay rate. Comparing these expressions with the transition probability (II.25) we see from (II.45) that and that Γ(η) is exactly the same as expression (II.26).

C. Unitarity
One of our main goals is to study the entanglement entropy from tracing over superhorizon degrees of freedom. Thus it is important to make sure that the loss of information encoded in the entanglement entropy is a genuine effect of the tracing procedure and not a consequence of approximations in the evolution of the quantum state. Unitarity follows from the set of equations (II.28), combining these with their complex conjugates it is straightforward to therefore with the initial conditions (II.33) it follows that Although this is an exact statement, we now show that the Wigner-Weisskopf approximation and its Markovian implementation maintain unitary time evolution.

III. PARTICLE DECAY: ENTANGLEMENT ACROSS THE HORIZON:
In the scalar theory described by eq.(II.20) the cubic interaction allows a single particle state |1 k to decay into two particle states |1 k− p ; 1 p [13,15]. To lowest order in the coupling the matrix element for this process is given up to an overall phase by Consider an initial single particle state |1 k at time η 0 . Upon time evolution in the interaction picture this state evolves into This is an entangled state in which pairs of particles with momenta k − p, p are correlated.
In particular if k is subhorizon and p is superhorizon, the quantum state (III.2) describes entanglement and correlation of particles across the horizon. The coefficients in the state (III.2) are the solutions of the (WW) equations, namely where the matrix element is given by eq.(III.1). We will focus on the asymptotic limit where The self-energy (II.36) is given by 1 were the matrix elements are given by eq.(III.1) leading to the result given by (II.23).
As discussed in detail in ref. [15], as ∆ → 0 the integral features infrared divergences from regions in which the momenta are superhorizon, namely pη, pη ′ ≪ 1 and | k− p|η, | k− p|η ′ ≪ 1. Both of these momentum regions yield the same infrared contribution as a single pole in ∆ [15], as can be seen as follows. For superhorizon modes (−pη ≪ 1) the mode functions (II.9) behave (up to an overall phase) as and for subhorizon modes −kη ≫ 1 Therefore for p ≪ (−1/η) and k ≫ (−1/η) the matrix element (III.1) becomes (up to an overall phase) (III.8) The contribution to the self-energy from superhorizon modes with p ≤ µ (−1/η) (with µ an infrared cutoff) yields The processes that contributes to leading order in ∆ is the emission of superhorizon quanta, depicted in fig. ( The simple rules to extract the leading order contribution in ∆ are given in ref. [15], where the cancellation of the infrared regulator µ from the contributions of the subhorizon modes, for which one can safely set ∆ = 0, is also shown in detail. In particular, the appendix of ref. [15] shows how the contribution of the subhorizon modes replaces the term ln[µ 2 ηη ′ ] → ln[k 2 ηη ′ ] which to leading order in ∆ can be written as 1+∆ ln[ The contribution from the region | k − p| ≪ µ yields an overall factor 2 in the self-energy so that to leading order in ∆ (as can be seen by rerouting the loop momentum) Using the result in eq.(II.41) we finally find that to leading order, where we have approximated α/2 z 2−2∆ 0 → 0 since −kη 0 ≫ 1 as the physical wavevector of the initial particle is deep inside the Hubble radius at the initial time and it is assumed to remain inside the Hubble radius during the evolution.

IV. ENTANGLEMENT ENTROPY
The pure state density matrix corresponding to the entangled state of eq.(III.2) is Now let us trace over the superhorizon physical wavevectors − pη 1. This leads us to the mixed state density matrix for modes whose wavelengths are inside the horizon during the evolution where the factor 2 accounts for the two regions of superhorizon (physical) momenta −pη < 1 and −| k − p|η < 1 which yield the same contribution, as can be easily seen after a relabelling of momenta. The entanglement entropy is the Von-Neumann entropy for the reduced density matrix, we find where the occupation numbers of the initial and produced quanta are given by Note that the unitarity relation in eq.(II.53) implies that as expected on physical grounds. At this stage it is important to highlight how unitarity is manifest to leading order in the ∆ expansion as this feature simplifies the calculation of the entanglement entropy considerably. The main point is that the unitarity relation (II.53) implies that the contribution of superhorizon modes is the dominant one. This can be seen clearly from the following arguments: to leading order in ∆ we can neglect the term 2∆ in the exponent of z in the solution (III.11) and in the term −∆ in the exponent of (−η) in the matrix element (III.8) for superhorizon modes. Now the coefficient where we used the definition of α given by eqn. (III.10) and changed variables of integration to η 1 = √ α/k y. This expression clearly exhibits that the contribution to |C p (k; η)| 2 from superhorizon modes, to leading order in ∆ can be written in the following factorized form: The dependence on ∆ is a manifestation of unitarity to leading order; if we compute the integral in eq.(IV.7) over superhorizon modes the ∆ in the numerator in eq.(IV.7) cancels the single pole in ∆ from the integral giving an O(1) contribution, which is what is necessary to satisfy the unitarity condition (II.53) to leading order in ∆. This result is similar to that found in the case of particle decay in Minkowski space time [36]: in this case the particles produced from the decay of a parent particle feature a Lorentzian distribution in energy, with width Γ the decay width of the parent particle and amplitude 1/Γ, so that the energy integral over the distribution is O(1). In ref. [36] it is proven to leading order in the perturbative expansion O(Γ) that this narrow distribution of large amplitude is the main reason for the fulfillment of unitarity to leading order in the Wigner-Weisskopf approximation. In the case of de Sitter space time, the distribution function of the particles produced with superhorizon wavevectors is ∝ ∆/p 3−2∆ whose momentum integral over the region of superhorizon momenta is also of O(1).
Thus in the limit ∆ ≪ 1 the sum p |C p (η)| 2 is dominated by the superhorizon momenta and from the unitarity relation (II.53) we find Although the integral in F [k; η] can be written in terms of error functions, the unitarity relation (II.53) and the result (IV.9) furnish a more direct evaluation. Consider To leading order in ∆, the sum is dominated by the superhorizon contributions from both regions of integrations p (−1/η) , | k − p| (−1/η) contributing equally, hence Then the factorized form (IV.7) for superhorizon modes, combined with eqn. (IV.12) leads to (IV. 13) and for −kη ≫ 1 and −pη ≪ 1 we find to leading order in ∆ (IV.14) the same result is valid in the region −kη ≫ 1 with −| k − p|η ≪ 1 by replacing p ↔ | k − p|.
The long wavelength limit of eq.(IV.14) requires a careful treatment. Since |C p (η)| 2 = n p (η) is the distribution function of particles, for a fixed volume V there is an infrared divergence in the occupation as p → 0. However, our goal is to trace over the superhorizon quanta from the decay since the initial conformal time −η 0 up to conformal time η → 0 − . This entails that the lower momentum cutoff is determined by the mode that just becomes superhorizon at the initial time, namely Now the calculation of the entanglement entropy is straightforward: let us consider where we have introduced x m = η η 0 (IV. 18) and changing integration variable to x = −pη It is now clear that we can set x m → 0 safely in I 1 and in the terms that do not feature poles in ∆ in I 2 . The terms in I 2 that feature the ln[x m ] and the (single) pole in ∆, namely ] yield the leading contribution for ∆, x m ≪ 1. Therefore for ∆ ≪ 1 and x m ≪ 1 we find for the entanglement entropy to leading order α is given in eqn. (III.10) and we have set −η 0 = 1/(a i H i ) and V = 1/(a 0 H 0 ) with a i , a 0 the scale factor, and H i , H 0 the values of the Hubble parameter at the beginning of inflation (i) and today (0) respectively, taking the physical volume today to be the Hubble volume, therefore a i H i /a 0 H 0 ≃ 1. The function Z − 1 + e −Z is manifestly (semi) positive and monotonically increasing, behaving as ≃ Z 2 /2 for Z ≪ 1 and as ≃ Z for Z ≫ 1. As η → 0 the entanglement entropy grows monotonically during the time evolution. We can Z[η] in terms of the number of e-folds since the beginning of inflation as where N e (η) is the number of e-folds during inflation at (conformal) time η and N T ≃ 60 is the total number of e-folds of the inflationary stage.
V. WAVE PACKETS: The discussion above treated the initial and product particles in terms of plane waves. However, given the existence of a horizon and the intricacies that can give rise to for nonlocalized states, we now generalize the treatment to the case of wave packets. Quantization in a finite volume V is used throughout. Fock states describing single particle plane wave states of momentum k, |1 k , are normalized such that (V.1) Localized single particle states are constructed as linear superpositions where C ( k; k 0 ; x 0 ) is the amplitude, normalized so that For a monochromatic plane wave C( k; k 0 ; x 0 ) = δ k, k 0 . The spatial wave function corresponding to the wave packet is given by The normalization (V.3) implies For a monochromatic plane wave it follows that Υ( x) is a volume normalized plane wave.
The total number of particles and average momentum of the wave packet are given by and respectively, where a † k ; a k are the creation and annihilation operators. If k 0 is identified with the average momentum of the wave packet we assume that and the isotropy of |C( k; 0, x 0 )| 2 .
As a specific example we consider Gaussian wave packets, where σ is the localization in momentum space. The spatial wave function is The spatial wave function is localized at x 0 with localization length 1/σ and the momentum wave function is localized at k 0 which is the average momentum in the wave packet and the momentum localization scale is σ. The plane wave limit is obtained by formally identifying σ/ √ π → 1/V 1/3 ; V → ∞.
In terms of these wave functions the overlap of two wave packets with different momenta localized at different spatial points is (V.11) In the limit σ → 0 the overlap becomes a Kronecker delta, and in particular for k 0 , q 0 ≫ σ it follows that the wavepackets are nearly orthogonal since the overlap is non-vanishing for ∆k = k 0 − q 0 ∼ σ so that ∆k/k 0 ≪ 1. From the identity (V.8) we can infer the following important property of these wave packets which will be useful below: Although this result is evident with the Gaussian wave packets (V.9) it is quite general for localized functions of k − k 0 . The wave packet description is easily incorporated into the Wigner-Weisskopf approach to the description of the full time evolution of the quantum state of the decaying parent particle. The interaction picture quantum state (II.27) is generally written as where the states |κ are multiparticle states, with the initial conditions C( k; k 0 ; x 0 ; η 0 ) = C( k − k 0 ; x 0 ) ; C κ (t = 0) = 0 , (V.14) where C( k − k 0 ; x 0 ) describe the localized wave packet of the single particle state at the initial time, for example (V.9). Generalizing the state (III.2) describing the time evolved state to lowest order in λ, to a wave packet localized at the origin in space with the gaussian profile (V.9), we can write with the initial condition with C( k − k 0 ; 0) given by (V.9).
Recall that our goal in this article was to obtain the entanglement entropy associated with the decay of single particle states with sub-Hubble physical momenta all throughout the inflationary stage, assuming that near de Sitter inflation lasts a finite time. Namely the physical wavelength of the single particle state is always deep within the Hubble radius during the evolution. A wave packet description of single particle states, therefore must be in terms of wave packets whose physical spatial localization scale is always much smaller than the Hubble radius. Hence, we will consider wavepackets that are i) sharply localized in comoving momentum with an average momentum k 0 with k 0 ≫ H; k 0 ≫ σ, the latter condition ensuring a sharp localization around k 0 , and ii) with comoving spatial localization scale 1/σ 1/H so that the wavepacket is localized well within the Hubble radius. Namely the condition for the wavepacket to describe single particle states with a sharp localization in momentum and with spatial localization length scale smaller than or of the order of the Hubble radius implies the following constraint: Furthermore, consistency in tracing over degrees of freedom with super-Hubble physical wavelengths requires that the wavepacket is mainly composed of components with comoving momenta corresponding to physical wavelengths that are always inside the Hubble radius throughout the near de Sitter stage. This condition requires −k 0 η ≫ −ση ≫ 1 so that components of the wavepacket with super Hubble physical wavelengths are exponentially suppressed. The Wigner Weisskopf method follows the steps described in detail above. The interaction Hamiltonian connects the single particle plane wave states |1 k with the two-particle plane wave states |1 k− p ; 1 p with matrix elements given by (III.1) leading to the set of equations Implementing the Markovian approximation as in the plane wave case with the initial conditions (V.16) we find where C k (η); C p (k; η) are the solutions of the Wigner-Weisskopf equations for plane waves, given by (III.11,III.4).
To obtain the reduced density matrix we would need to carry out the integration over the wavepacket variable k. The wave packet profile (as function of comoving wavevectors) is chosen to be sharply peaked at k 0 with a width σ ≪ k 0 . Therefore upon integration we can Taylor expand the integrand around k = k 0 and integrate term by term in the Taylor expansion in k − k 0 , because the wavepacket profile is a function of | k − k 0 | it follows that the corrections are a series in σ 2 /k 2 0 ≪ 1. An example of a quantity that must be integrated in k are the matrix elements (III.1), which upon being integrated with the wavepacket profile can be simply written as M(p; k 0 ; η) + O(σ 2 /k 2 0 ) + · · · . The same argument applies to the coefficients Therefore, to leading order in σ 2 /k 2 0 the reduced density matrix becomes We emphasize that the trace over the superhorizon modes leading to the reduced density matrix (V.22) has been carried out in the orthonormal plane wave basis.
Using the definition of the wavepacket single particle states (V.2) and the property (V.12) we finally find to leading order in σ 2 /k 2 0 ≪ 1 For k 0 ≫ p, σ the wave-packet states | k 0 − p, 0 contain plane wave components with subhorizon momenta ≃ k 0 − p since components with wavevectors that are very different from this value are exponentially suppressed. Therefore these wavepacket states are very nearly plane wave states with subhorizon momenta k 0 ≫ −1/η. Therefore, to leading order in σ 2 /k 2 0 , the reduced density matrix in terms of the wave packet single particle states features the same form as for the plane wave case with the only modification being the replacement of the single particle Fock states by the localized wavepacket states of single particles. As a corollary, to leading order in σ 2 /k 2 0 the entanglement entropy is the same either for localized wavepackets or plane waves.
The logarithmic dependence of the entanglement entropy (IV.20) on the volume factor has a clear statistical interpretation independent of whether the description is in terms of localized wavepackets or plane wave states. Consider a dilute gas of particles whose statistical distribution or phase space density is f p . The total density of particles is 24) and the Von-Neumann entropy of this (dilute) gas is If the number of particles remains finite in the large volume limit, namely if the particle density scales ∝ 1/V in this limit, then it follows that f p ∝ 1/V . On the contrary, if f p is independent of the volume as in the cases of the Maxwell-Boltzmann, Bose-Einstein or Fermi-Dirac distributions, the total density is finite in the infinite volume limit and the entropy is extensive. For a finite number of particles (vanishing particle density in the infinite volume limit) f p ∝ 1/V and the Von-Neumann entropy is not extensive, This is precisely the origin of the logarithmic dependence on the volume of the entanglement entropy: the initial state has one particle within a Hubble volume and the final state has one (of each) daughter particle, the distribution function of the daughter particles at asymptotically long times after the decay of the parent particle is |C χψ (p, ∞)| 2 ∝ 1/V the inverse volume dependence is the statement that there is a finite number of particles distributed in phase space 2 . Obviously this volume dependence is independent of whether the states are described by plane waves or wave packets, but is a statement of the simple fact that the number of particles in the volume V = (−1/η 0 ) 3 is finite. The dependence on the scale factor reflects the fact that more modes are crossing the Hubble radius, but the total number of particles described by these modes is still finite.

VI. CASCADE PROCESSES: THE WAY FORWARD
In the previous section we implemented the Wigner Weisskopf method to lowest order in λ 2 , but the method itself is much more general. It relies on a perturbative expansion, a truncation of the hierarchy at a given order in this expansion, and a resummation of the resulting self-energy terms that yield the long time asymptotics. For example, in quantum optics it has been implemented to study the cascade decay of many level atoms [45,54]. As shown in [41] this resummation is a real time version of the Dyson resummation of selfenergies and is equivalent to a dynamical renormalization group resummation of secular terms.
In this section we set up a roadmap to study higher order processes and along the way we exhibit the relation between the Wigner-Weisskopf method and the resummation of selfenergy diagrams and a diagrammatic expansion. Given the discussion on wavepackets in the previous section, we will restrict ourselves to treating the plane wave case.
The lack of kinematic thresholds in inflationary cosmology implies that the decay of quanta occur in a cascade process. For example with a cubic interaction as studied above, a state with a single quanta can decay into a state with two other quanta, in turn each one of the quanta in this state can decay into two other quanta, therefore a single particle state will decay via a "cascade": 1 → 2 → 3 → 4 · · · depicted in fig.(3). Each branch of the cascade corresponds to an interaction vertex and another power of the coupling, showing that the branches of the cascade are suppressed in perturbation theory. For example, the amplitude for 3 particles is down by a factor of λ (trilinear coupling) with respect to the two particle one, the four particle state is suppressed by another power of λ, etc.
To simplify notation, let us define matrix elements that connect a state of i quanta with a state of j quanta via the interaction Hamiltonian H I , (VI.1) Here [i], [j] describes the set of i, j quanta with different values of momenta. As studied above, we see that a state with a single quanta of (comoving) momentum k is connected via H I to a state with two quanta, with momenta q, k − q respectively. The matrix element for this process is [1]|H I (η)|[2] = 1 k |H I (η)|1 q ; 1 k− q where the set of values q defines the twoquanta states [2]. Thus the generic matrix elements between single quanta states and two quanta states in this set are [1] In what follows we will only consider connected diagrams or processes, neglecting disconnected diagrams which do not describe transitions but rather a renormalization of the vacuum state (for discussions see [41]). Consider the cascade decay of a single particle state |1 k into three particles, along with their inverse processes, neglecting the disconnected (vacuum) diagrams the typical sequence is shown in fig. (3) and the quantum state is given by |Ψ( k, η) = C 1 (k, η)|1 k + p C 2 ( k, p; η)|1 k− p ; 1 p + p, q C 3 ( k, p, q; η)|1 k− p ; 1 q ; 1 p− q + · · · (VI.2) The set of Wigner-Weisskopf equations are obtained straightforwardly as in the previous section. An important aspect in obtaining these equations is that a particular state with n particles with a fixed set of momenta has branched out from one "ancestor state", whereas it branches forward into an n + 1 particle state where the new particle has an arbitrary momentum that is summed over. As an example of this pattern consider the 3 particle state |1 k− p ; 1 q ; 1 p− q for a fixed value of p and q (the value of k is fixed by the initial state). This state branched out from the two particle state |1 k− p ; 1 p (up to relabelling the momenta and indistinguishability of the particle states), therefore it only has one "ancestor" as a consequence of momentum conservation. However, it branches out to 4 particle states of the form |1 k− p ; 1 q ; 1 l ; 1 p− q− l where the wavector l must be summed over.
The labels without brackets in the coefficients C n correspond to a particular state of n − particles with a fixed set of momenta compatible with total momentum conservation, whereas the sums over [n] are over the n-particle states compatible with the set of wavenumbers determined by momentum conservation. The terms shown in the hierarchy (VI.3,VI.4,VI.5) are the ones depicted in fig. (3) and their inverse processes: if the Hamiltonian connect the states [i] with the states [j] it also connects [j] back with [i], these are the inverse processes depicted in fig. (3). The two terms in eqns. (VI.4,VI.5) have an illuminating interpretation. The first terms correspond to the "population gain" of the states with two and three particles from the decay of their ancestors states with one and two particles respectively, while the second terms represent the "loss" or decay of the amplitudes into states with one more particle. Because of the initial conditions C 1 (η 0 ) = 1; C n =1 (η 0 ) = 0, it follows that d|C 2 | 2 /dη ∝ λ 2 + λ 3 + · · · ; d|C 3 | 2 /dη ∝ λ 4 + λ 5 + · · · so that the (conformal) time dependence of the coefficients also follows a hierarchy: the three particle state "fills up" on time scales ∝ 1/λ 2 larger than the two particle state, the four particle state on time scales ∝ 1/λ 2 larger than the three particle state, etc.
Let us consider truncating the hierarchy beyond the three particles intermediate state, namely set C [4] = C [5] = · · · = 0 along with all the other higher terms in the hierarchy. We then proceed to solve the equations from the bottom up with the initial conditions C 1 (η 0 ) = 1; C [2] (η 0 ) = C [3] (η 0 ) = · · · = 0. We obtain The first term in (VI.7) describes the build-up of the two-particle amplitude from the decay of the initial single particle state, whereas the second term describes the decay of the two-particle state into three particles via the cascade decay. Since the matrix elements are ∝ λ we can solve eqn. (VI.7) iteratively in perturbation theory up to the order considered in the hierarchy, in order to understand the time scales, To make the arguments clear, let us consider Minkowski space time and early time scales so that C 1 ≃ 1. Then the two particle amplitude builds up ∝ λt (with rate ∝ λ), and from eqn. (VI.6) we see that the three particle state builds up ∝ λ 2 t 2 ≪ λt, clearly reflecting that the population of the three particle state builds up much slower than that of the two particle state etc.
The build-up and decay integrals feature secular growth as η → 0 (long cosmic time), and the second step in the Wigner-Weisskopf method provides a non-perturbative resummation of these processes: writing (VI.7) as an integro-differential equatioṅ (VI.9) and introducing the Markovian approximation as in eqn. (II.37-II.39) (the second approximation in the Wigner-Weisskopf method) we find (VI.10) This compact expression reveals at once the build-up of the amplitude from C 1 and the eventual decay of the two-particle state encoded in γ 2 (η). A simple perturbative expansion of this expression up to O(λ 4 ) reproduces (VI.8) consistently with the Markovian approximation.
The last step is to insert this solution into (VI.3), solve the integro-differential equation for C 1 and insert this solution into (VI.8) and (VI.6) respectively. Obviously this procedure leads to a very complicated expression that is not very illuminating. However progress can be made by introducing the perturbative solution (VI.8), leading to the following integrodifferential equation for C 1 : [3]   The dashed lines cut through multiparticle states and indicate similar rules to the Cutkosky rules of quantum field theory that relate the absorptive parts of self-energy diagrams to intermediate multiparticle states.
In order to make progress in the solution of (VI.11) the second part of the Wigner-Weisskopf method invokes a Markovian approximation, just as that described in section (II B) implemented to lowest order. This approximation is again justified in a weak coupling expansion and is the statement that η derivatives of the coefficients are "slow" and can be systematically expanded perturbatively. The procedure follows the steps described by eqns. (II.37-II.40), integrating by parts the kernels of the integrals and keeping consistently up to O(λ 4 ) we findĊ The second term in the bracket in the first line arises from the derivative expansion of the term with the one-loop self-energy (see eqn. (II.39)); in ref. [41] this term is identified as a contribution to wave function renormalization. Therefore This expression provides a non-perturbative resummation of self-energies in real time up to two loops and includes the decay of the initial state into intermediate states with two and three particles. In Minkowski space time, the initial state decays as ∝ e −Γt with Γ = λ 2 γ 2 + λ 4 γ 3 + · · · corresponding to the contribution to the self energy from the two particle intermediate states (one loop) three particle intermediate states (two loops) etc, highlighting that the probability of production of the two particle intermediate state occurs on a time scale ∝ 1/λ 2 , that of the three particle intermediate state on ∝ 1/λ 4 etc.
Clearly the decay into two particle states occurs on shorter time scales as this process corresponds to the one-loop diagram, whereas decay into three particles occurs on much slower scales at this process corresponds to the two-loop contributions.
It remains to insert this solution into (VI.8) and in turn insert the solution for C [2] ,into (VI.6). Because the matrix elements M ij ∝ λ it follows that if we take C 1 ∝ O(λ 0 ), then C [2] ∝ λ ; C [3] ∝ λ 2 · · · . The quantum state obtained from the decay of a quanta with momentum k is given by C [2] (η)|[2] + [3] C [3] (η)|[3] + · · · , (VI. 15) The states |[2] = |1 p ; 1 k− p and |[3] = |1 p 1 ; 1 p 2 : 1 k− p 1 − p 2 and the sums over [2], [3] are over p and p 1 , p 2 respectively. Thus the probability of a given two particle state is given by |C 2 (η)| 2 ∝ λ 2 , of a given three particle state is |C 3 (η)| 2 ∝ λ 4 , etc. This is exactly as in the case of multiphoton processes in quantum electrodynamics or of an atomic cascade of a multilevel atom. Each photon in the final state is associated a probability that is proportional to α em , so multiphoton processes are suppressed by powers of the fine structure constant. In this case multiparticle final states are suppressed by powers of λ 2 for each extra particle in the final state. This is also the case in atmospheric air showers where very energetic particles decay via a cascade process where each branch of the cascade is down by a power of the coupling to the respective channel.
In the case of cascade decay in Minkowski space time, the probability of finding particles from a particular decay channel is given by the branching ratio of such channel Γ c /Γ tot , namely ratios of different powers of the couplings. Our result obviously entails the same physics: the probability of a state with three quanta is suppressed by λ 2 with respect to that with only two quanta, etc.
Furthermore, the explicit form of W (η) (VI.13) clearly shows the separation of time scales: the decay into two particles involves time scales ∝ 1/λ 2 and is determined by the product of matrix elements M 21 M 12 whereas the time scales for decay into three particle states is determined by the last term in (VI.13) which implies time scales ∝ 1/λ 4 . Therefore there is a hierarchy both in the probability of multiparticle states and the time scales associated with their production from the decay of the parent particle. The cascade decay processes are controlled by the small coupling λ.
The entanglement entropy can now be calculated by obtaining the reduced density matrix by tracing over the superhubble degrees of freedom in the pure state density matrix |Ψ( k, η) Ψ( k, η)| and is a straightforward implementation of the steps described in the previous section with the technical complication of the integration over the super Hubble subset of momenta in the multiparticle contributions. This is only a technical difficulty but not a conceptual roadblock, since the contribution to the entanglement entropy from higher multiplicity states will be suppressed by high powers of the coupling λ. An illustrative example in Minkowski space time is the cascade decay π − → µ − ν µ ; µ − → e − ν e ν µ : whereas the pion decays on a time scale ≃ 2.8 × 10 −8 secs the muon decays on a time scale ≃ 2.2 × 10 −6 secs therefore during a long time interval 10 −8 secs ≤ t ≤ 10 −6 secs the two particle state |µ − , ν µ yields the largest contribution to the quantum state.
Furthermore, the unitarity relation (II.48) entails that |C [2] (η)| 2 + [3] |C [3] (η)| 2 + · · · = 1 , (VI. 16) which was confirmed in the previous section to leading order in the coupling and ∆. In summary: the cascade decay is controlled by the perturbative nature of the interaction, the probability for multiparticle states being suppressed by powers of the coupling constant and the time scales associated with the formation of multiparticle states widely separated by larger powers of 1/λ. Furthermore, in the case under consideration here, the physical momentum of the initial state is taken to remain deep inside the Hubble radius at all times during inflation. At any large but fixed (conformal time) the initial state maintains a small but non-vanishing population, a two particle state being populated with probability λ 2 , a given three particle state with probability ∝ λ 4 etc. Therefore if the quasi-de Sitter inflationary stage lasts a finite (say ≃ 60) number of e-folds, the quantum state will be a linear superposition of many particle states and unitarity implies that each state features a perturbatively small population. An interesting and conceptually puzzling situation arises in the case of eternal de Sitter, since in this case, at asymptotically long times all states would have decayed to vanishing probability in clear contradiction with unitarity, but in this case all physical momenta eventually also become superHubble. Perhaps this puzzling aspect is related to the intriguing results of ref. [11] and deserves to be studied further.
While we have established a roadmap and a "proof of principle" of the method, undoubtedly there are several aspects that merit a deeper study such as infrared enhancement from superhorizon modes, the issue of unitarity in eternal de Sitter, the detailed aspects of the (conformal) time dependence of the amplitudes of multiparticle states etc. We postpone the study of these more technical details of the higher order processes to a future article.

VII. POSSIBLE RELATION TO NON-GAUSSIANITY.
The cubic interaction vertex suggests a relation between the decay amplitude (see fig.(1)) and the non-gaussian bispectrum which is the three point function of the field. The relationship with the bispectrum becomes more clear by introducing G(k, η, η ′ ) = g * ν (k; η)g ν (k; η ′ ) from which it follows that The imaginary part of the η ′ -integral is proportional to the bispectrum of the scalar field [47,48]. The main difference is that the self-energy is the integral over one of the momenta.
In particular the leading order in ∆, namely the contribution from the infrared enhanced, superhorizon modes, is determined by the highly squeezed limit shown in fig. (5), which corresponds to the local limit of the non-gaussian correlator. This connection highlights that this local limit is describing correlations between subhorizon and superhorizon modes, these are the correlations that yield the entanglement entropy upon tracing over the superhorizon degrees of freedom.
There are important differences between the scalar field theory with cubic interaction studied here, and the cubic interactions of curvature perturbations in the theory of nongaussian fluctuations [47,48], the main difference being both spatial and (conformal) time derivatives in the interactions. However, the study of ref. [49] found that transition probabilities of curvature perturbations (in single field slow roll inflation) are suppressed by slow roll parameters but enhanced by infrared logarithms, which are similar to those emerging in the ∆ → 0 limit in our study (corresponding to massless fluctuations), thus suggesting that the results obtained above may apply to the decay of curvature perturbations and superhorizon entanglement and concomitant entanglement entropy.

VIII. CONCLUSIONS AND FURTHER QUESTIONS
In inflationary cosmology all particle states decay as a consequence of the lack of a global time-like Killing vector which would in turn enforce kinematic thresholds. In this article we have studied the entanglement entropy from the decay of single particle states during de Sitter inflation in a theory of a light scalar field with M 2 ≪ H 2 and cubic interactions. The quantum state that describes the single particle decay and the produced particles is a two-particle state entangled by momentum conservation. We have extended and generalized the Wigner-Weisskopf method used in the treatment of spontaneous decay of atomic states to the realm of quantum field theory in an expanding cosmology, and implemented this method to obtain the quantum state that describes the decay of the parent particle and the production of the daughter particles. We showed in detail that this non-perturbative approximation is manifestly unitary. The amplitudes for the two-particle entangled state features infrared enhancements that are manifest as poles in ∆ = M 2 /3H 2 as a consequence of the emission of superhorizon quanta and we implement a consistent expansion in ∆ to leading order to obtain the (pure state) density matrix that describes the decay of the parent and production of daughter particles. When the parent particle's wavelength is inside the horizon, the density matrix elements for the produced particles are dominated by the contribution of superhorizon momenta of one of the daughter particles, describing entanglement, correlation and coherences across the horizon. Tracing the pure state density matrix over the superhorizon modes we obtain a mixed state density matrix from which we calculate the Von Neumann entanglement entropy, which describes the loss of information from the correlations between sub and superhorizon modes due to the non-observation of these latter states. We find that the entanglement entropy is enhanced in the infrared by a factor of ln[1/∆] and grows logarithmically with the physical volume as a consequence of more modes crossing the Hubble radius during the inflationary stage.
The generalization to the description of single particle states in terms of wavepackets spatially localized within the Hubble radius but localized in momentum was provided. Under the conditions that the average wavevector of the wave packet be associated with subHubble wavelengths all throughout the near de Sitter stage, we showed the equivalence between the plane wave and wave packet description and assessed the corrections in terms of the ratio of the width of the wavepacket in momentum space and the average momentum associated with the single particle state.
The lack of kinematic thresholds implies that particle decay occurs in a cascade process, namely 1 → 2 → 3 · · · . We have extended the Wigner-Weisskopf method to establish a framework to study the cascade decay and analyzed in detail the process up to a three particle branching in the cascade, but the results are quite general. We showed that for weak coupling (here we considered a cubic coupling) the probability of multiparticle states is suppressed by powers of the coupling, for example in the case of cubic coupling the three particle state is suppressed by O(λ 2 ) with respect to the two particle state, the four particle O(λ 4 ) etc. We have established a relation between the different multiparticle processes and higher order loop contributions in the self-energy, just as in the case of Cutkosky rules in Minkowski space-time. This relation clearly shows that just as the probability of higher multiparticle states is suppressed by high powers of the coupling, the time scales for decay into higher multiplicity states are widely separated by inverse powers of λ 2 . Therefore the cascade decay is controlled by the weak coupling, just as multiphoton processes in QED.
Further questions: a): An important feature of inflationary cosmology is that physical wavelengths that cross the Hubble radius during inflation, re-enter the Hubble radius (now the particle horizon) during radiation or matter domination and these quantum fluctuations are the seeds of temperature anisotropies and inhomogeneities.
The entanglement entropy that we have studied is a measure of the correlations between the entangled subhorizon and superhorizon degrees of freedom as a consequence of interactions, which brings the question of whether upon re-entry the fluctuation modes that were superhorizon during inflation "bring back" the quantum correlations and if so how are these manifest in the power spectrum of the CMB? Furthermore, going from quantum fluctuations of the curvature (or gravitational potential) to temperature fluctuations entails replacing quantum averages by statistical averages. Thus it is a relevant question whether this statistical averaging includes the quantum correlations from entanglement. Last but not least, if the quantum states can decay, it is conceivable that the lack of power in the low multipoles which is present in the cosmological data and has been persistent in the statistical analysis of WMAP7 [50], WMAP9 [51] and Planck [52] which reports a power deficit at low multipole with 2.5 − 3 σ significance and a recent statistical analysis of the combined dataset [53], may be due to the decay of the quantum fluctuations during the inflationary stage. Our study applies to a scalar field in de Sitter space time and in order to answer this question the analysis presented here must be applied to the case of scalar perturbations.
b:) Furthermore, and in relation with the question above, it is a tantalizing possibility that the superhorizon correlations become manifest as intensity correlations leading to interference phenomena akin to the Hanbury-Brown-Twiss effect discussed in ref. [20]. c): how the infrared enhancements modify the higher order multiparticle processes in the cascade decay, and how unitarity is manifest in the (formal) case of eternal de Sitter inflation.
While at this stage we do not yet see a clear observational consequence of the entanglement entropy beyond the theoretical conceptual aspect of information loss from the correlations and superhorizon entanglement, the exploration of potential observational consequences along with the questions raised above are worthy of further and deeper study, on which we expect to report in the future.
This study of the superhorizon entanglement entropy from particle decay bridges two concepts previously explored in the literature: the entanglement between spatially separated but correlated regions, in our case the correlations between sub and superhorizon quanta of the daughter particles, akin to the superhorizon correlations studied in ref. [34], and the momentum-space entanglement studied in ref. [35,36]. In our study the entanglement entropy is a result of both types of concepts, linked together by the interactions but with the distinct aspect of being a non-equilibrium process as a consequence of the cosmological expansion.