Quantum non-linear evolution of inflationary tensor perturbations

We study the quantum mechanical evolution of the tensor perturbations during inflation with non-linear tensor interactions. We first obtain the Lindblad terms generated by non-linear interactions by tracing out unobservable sub-horizon modes. Then we calculate explicitly the reduced density matrix for the super-horizon modes, and show that the probability of maintaining the unitarity of the squeezed state decreases in time. The decreased probability is transferred to other elements of the reduced density matrix including off-diagonal ones, so the evolution of the reduced density matrix describes the quantum-to-classical transition of the tensor perturbations. This is different from the classicality accomplished by the squeezed state, the suppression of the non-commutative effect, which is originated from the quadratic, linear interaction, and also maintains the unitarity. The quantum-to-classical transition occurs within 5 - 10 e-folds, faster than the curvature perturbation.


Introduction
While inflation [1,2,3] is postulated to resolve initial condition problems of the early universe, it also gives a natural account of how to generate the initial perturbations on super-horizon scales as observed in the cosmic microwave background (CMB) [4]. That is, during inflation the primordial quantum fluctuations leave the horizon and become classical perturbations as we observe today such as the temperature anisotropies in the CMB and inhomogeneously distributed galaxies [5,6,7,8,9]. As for this aspect, the quantum-to-classical transition is known to arise as the states for the scalar and tensor perturbations take a form of specific coherent superposition of excitations, called the squeezed state [10]. More concretely, through time evolution, some quadratic terms reflecting curved space-time make the state "squeezed" (see Appendix A) such that the non-commutative effect between the perturbation field variable and its canonical conjugate momentum becomes negligibly small. The probability density of the squeezed state is not localized in both the field variable and its canonical momentum, but spreads over the phase space. Instead, the suppression of the quantum mechanical nature of non-commutativity allows us to interpret the probability distribution as a classical one, just like the semi-classicality of the Wentzel-Kramers-Brillouin (WKB) approximation. Such a squeezing of the state has been discussed in various contexts in physics (see, e.g. [11,12]). As an application to the inflationary cosmology, generic massless particle case was discussed [13] and analyzed in terms of tensor [11,14] and curvature [15] perturbations, respectively.
On the other hand, the presence of the horizon in (quasi) de Sitter (dS) space-time during inflation enables us to separate the quantum fluctuation modes into those shorter and longer than the horizon scale according to their wavelengths. When we focus on the modes relevant for observations, say, the modes with wavelengths longer than the horizon scale, they are evolving non-linearly in the time-dependent background due to gravity, interacting with subhorizon modes which are not very relevant observationally. In that sense, the system of the super-horizon modes is open. After tracing out unobservable modes, or environment, the nonlinear interactions between modes in system and environment introduce the so-called Lindblad operators [16] which violate the unitary evolution of the density matrix. Since the Lindblad terms are equivalent to adding new Hamiltonian terms with randomly varying source [17], they have a close connection to the stochastic inflation formalism [18]. Indeed, inclusion of the Lindblad terms provide the most generic evolution of the density matrix under several assumptions: the evolution equation of the density matrix is Markovian, i.e. the density matrix at some specific time depends only on that at an earlier time, rather than that over a range of times, and linear in the density matrix. Moreover, the resulting density matrix is Hermitian, and has a unit trace with a positive expectation value, satisfying probability interpretation for the "mixed state" [19].
It has been known that the time evolution regulated by the Lindblad opeartors leads to the loss of quantum nature during inflation in another way. One way to see this is to investigate the "decoherence" -the decay of the interference effect [20,21,22], mimicking collapse of the wave function into some specific eigenstate of the observable (for a review on decoherence in generic context, see e.g. [23,24]). Recent studies [25,26,27,28,29] on the Lindblad operators during inflation emphasized this aspect. Meanwhile, the Lindblad operators also reveal the classicality of the system by converting the pure state -for inflation, the squeezed state -into the mixed state. While the Lindblad operators introduce non-unitary evolution of the pure The evolution of the density operator, with interactions, can be written in a different manner in which picture one is using. One general remark however is that the density matrix is not an usual operator, in the sense that it follows the von Neumann equation where the sign is opposite to the standard Heisenberg equation. We basically work in the Schrödinger picture. The reason is, despite of the existence of higher-order interactions, the evolution of the whole system is contained in the states in the Schrödinger picture, from which we construct the density matrix of our interest. Since in general the Hamiltonian is time-dependent and not commuting at different times, the density matrix in the Schrödinger picture is given by, up to second order explicitly in interaction, where is the free evolution operator with T (T ) denoting (anti-) time-ordering operator, is the "free" density matrix in the Schrödinger picture at τ in the absence of interaction, and H int,I (τ ) is the interaction Hamiltonian in the interaction picture. We are eventually interested in the "reduced" density matrix where the effects of the fast modes, or the environment, are integrated out. So we take the trace over the environment state |E i at the time τ , not the state |E(τ ) = U 0 (τ ; τ 0 )|E(τ 0 ) which has evolved from the initial environment state |E(τ 0 ) , with i being an abstract index for the a possible choice of the basis states for environment at τ : Being a choice of the basis states, |E i is complete: The time evolution of this reduced density matrix can be found as: since the environmental states we trace out are fixed at the time of tracing out the environment, |E i does not carry time-dependence. Thus to obtain the evolution equation for the reduced density matrix, we only need to take the trace of that for the full density matrix (however, see the discussion below). To proceed, we note that for H 0 we can completely separate the system and environment since we are considering free evolution, so H 0 = H 0,S + H 0,E with H 0,S , H 0,E = 0. To trace over |E i , we note that since |E i is the environment state at a later time τ , it has no overlap with the initial system state |S(τ 0 ) , i.e. |E i and |S(τ 0 ) are orthogonal. Meanwhile, the initial environment state |E(τ 0 ) may still have ample overlap with |E i , so only |E(τ 0 ) responds to |E i . Then, after some calculations, we find the evolution equation for ρ red as (see [28,29] for derivation) H 0,S , being the free-evolved density matrix of the system states, in the first term in (2.9), we have defined (2) (τ )|E i . We can write this term purely in terms of the commutators with free Hamiltonian of the system so the evolution of the reduced density matrix described by this term is unitary. Now, with time-dependent background there is ambiguity how to choose the system and environment states. We may take the state |E i to be traced out as the environmental Fock space, i.e. the modes on sub-horizon scales fixed at the moment τ , viz. those with k > aH = −1/τ . This is quite natural, but there is no obvious sharp distinction between the initial system and environment states |S(τ 0 ) and |E(τ 0 ) . Different choices could well be justified for different reasons. One choice that is at least computationally appealing is such that we choose |E(τ 0 ) in such a way that this corresponds to the modes up to which becomes the horizon scale k = aH = −1/τ at the moment τ , and |S(τ 0 ) denotes all the modes with smaller k. Then with (2.7) and E(τ 0 )|E(τ 0 ) = 1, some expressions become very simple, e.g. ρ • In the second term in (2.9), we have defined where H int,S is the interaction Hamiltonian in the Schorödinger picture and L i 's are the socalled Lindblad operators defined below. While written purely in terms of the commutator form, one thing we can further note is the existence of the "shift" in the Hamiltonian, H int . This means even the unitary part of the evolution is described by a Hamiltonian not identical to the unperturbed free Hamiltonian of the system, because the environment perturbs the free Hamiltonian in such a way that the free Hamiltonian that describes the unitary evolution should be in general shifted. Note that as can be inferred, H (eff2) int and the Lindblad terms given below originate from the same term, • Lastly, in the third term in (2.9) we have defined the Lindblad operators as (2.14) They come from the pure commutator contributions between the density matrix and the interaction Hamiltonian, and are responsible for the non-unitary evolution of the reduced density matrix: since the third term in (2.9) is real, this part does not lead to a simple sinusoidal oscillations but gives rise to an exponential change in some elements of the reduced density matrix. Or, some information on the environment is integrated out, now the whole system is in the so-called mixed state, or decoherence occurred. This is described in several articles for the scalar cosmological perturbation, i.e. the curvature perturbation [25,26,27,28,29]. Also note that the trace of the Lindblad terms is obviously zero, meaning that the total classical probability is preserved. This leads to interesting consequences as we will see in Section 5.
Before closing, it is illustrative to consider what if the contributions are only from the system sector. Then, the environment state |E i commutes with H int , so (2.13) and (2.14) become Then we can easily find where ρ S (τ 0 ) ≡ S(τ 0 ) S(τ 0 ) . What happens is, not surprisingly, that we now have the unitary evolution of the pure system part in the precisely equivalent form to the standard von Neumann equation. That means, if we consider the contribution that purely belongs to the system sector, not including any from the environment sector, such a contribution can be absorbed into the unitary evolution, not contributing to the non-unitary Lindblad operators.
That is, any effect in the system sector that comes from integrating out the environment sector is visible only when we include the contribution from the environment sector in the interaction Hamiltonian.

Interaction Hamiltonian for tensor perturbations
Now we consider the interaction Hamiltonian for tensor perturbations. We consider as the leading higher-order interaction cubic order Hamiltonian that contains purely tensor perturbations h ij . The cubic action for the tensor perturbation is [32,33] where H ≡ a ′ /a and ǫ ≡ −Ḣ/H 2 . There are two things to notice: 1. Some terms are not slow-roll suppressed, and 2. some terms do not possess any derivatives.
Then, the interaction Hamiltonian we consider is As we are interested in the Fourier modes, going to the momentum space using the free solutions (A.4), (A.5) and (A.15) gives the following explicit expression for the interaction Hamiltonian: where and similar for c 22 and c 23 , (3.5) and similar for c 32 and c 33 . Before we make further steps forward, let us investigate the time dependence of (3.3). Here, the creation and annihilation operators a λ i k i (τ ) and a λ i † −k i (τ ) are coming from the mode function solution h λ (k i , τ ), and are time-dependent -this is clearly not the case for the Schrödinger picture. In fact, promoting the canonical conjugate pair π k and v k to time-dependent quantum operators as in (A.10) means we make use of the Heisenberg picture. The general transformation of an operator A from the Schrödinger picture to the interaction picture is determined only by the free Hamiltonian as A I (t) = U † 0 (t; t 0 )A S U 0 (t; t 0 ) 1 , also applies to the creation and annihilation operators as a k,I (τ ) = U † 0 (τ ; τ 0 )a k (τ 0 )U 0 (τ ; τ 0 ) (3.7) and the same for a † −k (τ ), and the explicit form is given in terms of the Bogoliubov transformation (A.17). That is, (A.17) tells us how the time-independent -in the sense that it is defined at τ 0 and is fixed there -Schrödinger picture operators a k (τ 0 ) and a † −k (τ 0 ) are related to the timedependent interaction picture operators a k (τ ) ≡ a k,I (τ ) and a † −k (τ ) ≡ a † −k,I (τ ). Returning to (3.3), as mentioned above the creation and annihilation operators are from h λ (k i , τ ), where the relations to the operators at the initial moments are given by (A.17). That is, the evolution of the creation and annihilation operators in (3.3

) is determined by the quadratic Hamiltonian
If we are in the Heisenberg picture, the evolution of an operator -or the transformation from the Schrödinger picture to the Heisenberg picture -is determined by the full Hamiltonian, H 0 , and thus (3.3) is in the interaction picture. Note that Here, H int,S (τ ) is the interaction Hamiltonian in the Schrödinger picture in the sense that the operators are those at τ 0 , while the coefficients are evaluated at τ and the relation between them is given by the standard form: H int,I (τ ) = U † (τ ; τ 0 )H int,S (τ )U(τ ; τ 0 ). Another benefit to make use of the initial creation and annihilation operators at τ 0 is that, after all, we wish to compute the time evolution of the reduced density matrix of the system ρ red , for which we need to find the Lindblad operators (2.13) and (2.14). As we can see, they are essentially the interaction Hamiltonian sandwiched between the time-evolved environment vacuum states. Thus, combined with time evolution operators accompanied by H int,I , the initial annihilation operators directly eliminate the initial environment vacuum state |E(τ 0 ) , reducing amounts of computational efforts for the Lindblad operators. While we have the formal relation between creation and annihilation operators at τ and those at τ 0 as (3.7), practically we can also use (A.17) conveniently to write Then we can write (3.3) straightforwardly in terms of the combination of the coefficients at τ and the operators at τ 0 as where for the last equality we have set q ≡ −k. Since the coefficients c i(j) in (3.3) are basically combinations of polarization tensor, we can see for example where for the last equality we have renamed the dummy integration momentum variables as k 3 ≡ k ′ 1 and k 1 = k ′ 3 . Similar relations holds for other coefficients and other operator combinations, e.g.
Then we can see that the Hermitian conjugate pair of each term of (3.3) can be written as one another term, e.g. (3.14) Furthermore, since the structure of c 2i and c 3i is of the form shuffling the dummy momentum variables shuffles these coefficients, e.g. k 2 → k 1 , k 3 → k 2 and k 1 → k 3 gives where h 0 and h 1 are the coefficients of a λ 1 k 1 a λ 2 k 2 a λ 3 k 3 and a λ 1 † −k 1 a λ 2 k 2 a λ 3 k 3 in (3.10) respectively. Note that in the last step we have written the half of H int,I as the Hermitian conjugates of the other half. As emphasized before and written explicitly above, the coefficients h 0 and h 1 are at τ yet the operators are at τ 0 , viz. h 0 = h 0 (τ ) and a λ 1 k 1 a λ 2 k 2 a λ 3 k 3 = a λ 1 k 1 a λ 2 k 2 a λ 3 k 3 (τ 0 ) and so on.
4 Structure of Lindblad operators

Contributions of the environment and system sectors
Having found the explicit form of the interaction Hamiltonian (3.18), now we can compute the Lindblad operators explicitly. Before then, however, it is helpful to overview their structure rather than getting involved detailed calculations. For this purpose, we need to go back to the point how the Lindblad operators are constructed. As we can see from (2.13) and (2.14), the Lindblad operators L 1 and L 2 are the interaction Hamiltonian sandwiched between the time-evolved initial environment states |E(τ 0 ) and the final environment states |E i which we want to trace out. For example, we may choose |E(τ 0 ) in such a way that |E(τ 0 ) denotes the sub-horizon modes at the moment τ . Then, as (A.31), starting from the standard Bunch-Davies vacuum |0 0 where there is no excited particle for all k, we may split |0 0 as a sum of the infrared and ultraviolet parts: That is, instead of writing the initial states for the system and environment respectively as |S(τ 0 ) and |E(τ 0 ) , we now write them as |0 S and |0 E since we take the standard assumption that the initial vacuum is the Bunch-Davies one |0 0 for all modes. The infrared and ultraviolet vacua |0 S and |0 E are annihilated respectively by the annihilation operators of the infrared and ultraviolet modes, a k S and a k E where Thus separating the infrared and ultraviolet sectors of the interaction Hamiltonian enables us to treat one independent from the other. Consider, for example, H int,I (τ ) is of the following form: Upon applying the ultraviolet vacuum state |0 E , the corresponding ultraviolet annihilation operators are operational and there are only left the annihilation operators for the infrared modes: More generally, based on the form (3.10), let us schematically write H int,I (τ ) as where the overall coefficient is at τ , and the system (or infrared) and environment (or ultraviolet) sectors of H int,I , which are combinations of the three initial creation and annihilation operators, are at τ 0 so only the environment sector H E,m (τ 0 ) is responding to the initial environment vacuum |0 E . It is very important to note that the right-hand side of (4.6) denotes a single representative term in H int,I , and is not Hermitian, as clear from (3.10) -while H int,I itself is Hermitian also as clear from the same equation. In fact, to maintain the Hermiticity of H int,I we had better always think of not a single term but a single Hermitian pair. But for brevity we keep writing a single term as (4.6) while keeping in mind that we should remember the right-hand side of (4.6) is not Hermitian. Then (4.7) While H int,I is Hermitian, to keep the Hermitian-conjugateness of the Lindblad operators we write the Hermitian conjugate symbol for H int,I as H † int,I . The part sandwiched between the environment vacuum states above is some number -exactly speaking, some function of momenta, including zero as well. Note that for the Lindblad operator terms in (2.9), as (4.6), we may decompose since the time integration only acts on the coefficients: Here, we have intentionally used different notations for the environment (system) sector as G E(S),m , despite of the same time dependence at τ 0 , to avoid confusion in ordering the operators in L 1 (corresponding to H) and L 2 (corresponding to G). After simple calculations, we can schematically write the Lindblad terms of (2.9) as To compute the environment parts of the Lindblad terms, naively thinking, it seems that only the operator combinations that consists of three annihilation operators first and three creation operators next survive, e.g.
But given that the environment and system sectors of an operator are decomposed as (4.4), it needs more care. There are surviving contributions in every term of (3.18). For example, the first term of (3.18) contains no creation operator, so it seems that if both operator sets are of this type there would be no remaining contribution. But we do have a surviving contribution: with O S being the system operators not sandwiched between |0 E like ρ S and U 0,S (τ ; τ 0 ), then even if both two operator sets consist of annihilation operators only, we have where now the momentum integrations are only restricted to the system domain. But as argued in (2.17), those contributions in which the interaction Hamiltonian belongs only to the system sector can be absorbed into the unitary evolution and do not contribute to the Lindblad operators, so we do not consider such terms. After taking into account all possible contractions between operators in the environment sector, the contributions appear as momentum-dependent (and also polarization-dependent) coefficients in front of the system operators. Thus the naive density matrix elements of system do change under the influence of the environment sector, as is what the Lindblad operators describe. But considering the environment sectors only does not seem to tell us any more, e.g. the structure of the time-dependent coefficients like h 0 (τ ), since what is constrained from the environment sectors are: • Some momenta in the interaction Hamiltonian in L 1 and L 2 are related, e.g. k 1 = −k 4 , and • certain momenta belong to certain domains, e.g. k 1 and k 2 are environmental (i.e. |k 1 | ≫ aH and |k 2 | ≫ aH at τ ) while k 3 belong to the system sector (i.e. |k 3 | ≪ aH at τ ) and so on.
Thus, to make steps forward, we have to consider explicitly the system sector contributions. As we have started for the environment sector with the Bunch-Davies vacuum |0 E , so is the system vacuum |0 S . With ρ S (τ ) given by (2.10), as we can see in (4.9) the evolution operator U † 0,S (U 0,S ) in front of (behind) ρ S (τ ) cancels the counterpart in ρ S (τ ), so the initial vacuum state |0 S ( 0| S ) is exposed. Then the initial system creation or annihilation operators at τ 0 from the Lindblad operators can act directly on |0 S .
During the computation in the system sector, it is convenient to arrange the Lindblad terms in their simplest form in the sense that there is no operator or are only creation ones on the left of |0 S , and likewise there is no operator or are only annihilation operators on the right of 0| S , i.e. each term in dρ red /dτ is of schematic form where ρ mn is the time and momentum dependent coefficient for the contribution with m creation and n annihilation operators, with m, n ≤ 6. This form of the Lindblad term is irreducible, and applying the basis vectors on the left and right should give us the matrix elements in such bases directly.

Matrix elements of reduced density matrix
Now, we need to specify the "matrix bases" and to compute the matrix elements of the reduced density matrix with respect to those bases, or the states relevant for observations. The real issue here is, rather than computing the matrix elements which is quite straightforward, identifying what should be the observationally relevant states. Usually, the coherent state as an eigenstate of the annihilation operator is considered to be relevant for two reasons. First, it is the state with which the expectation values of operators follow the classical equation of motion. Second, the uncertainty is minimized in the coherent state. Thus, the coherent state coincides with our intuition of the classical solution. Especially, by taking expectation value of the field operators in terms of the coherent state 2 , we can treat any operator in the Hamiltonian, say A i , as a function of classical fields: where we distinguish the field operators (hatted) and the classical fields (unhatted). But as we have seen, the Lindblad terms are after computing out the system sector written in terms of the squeezed states: That is, the Lindblad terms are schematically of the form (4.12). Thus the squeezed states are already written in terms of the basis 15) and in this sense the squeezed states are the natural bases so that the calculations become particularly easy if we adopt the squeezed states as the pointer bases 3 . But obviously, this never means that we should adopt the squeezed states as the pointer bases. Rather, it just means that by adopting the squeezed states calculations becomes especially trivial. One possible argument in favour of the squeezed states as the pointer bases is from noticing that we do not observe directly primordial perturbations. What we do observe in reality is the temperature fluctuations δT /T 0 in the CMB and the inhomogeneous distribution of galaxies. According to the standard cosmology, they are originated from the initial conditions. Being (believed to be) a Gaussian random field, which is in very good agreement with recent 2 Obviously, the coherent state is different from an eigenstate of the field operator. The reason we can take π φ = 0 is: we are dealing with real fields so the expectation value of the operator π φ vanishes, since π φ ∼ a−a † [see (A.15)] and if a has a real eigenvalue, a † in the expectation value is given by the same eigenvalue. Furthermore, uncertainty is still non-zero even if it is minimized, which means that as soon as we consider the eigenstate of the field operator, canonical momentum probability distribution spreads so we can no longer fix π φ = 0.
3 Unlike the basis for the environment discussed in Section 2, the squeezed state basis is time dependent. However, such a time dependence corresponds to the unitary evolution of the Fock space states, so it is irrelevant to the non-unitary evolution through the Lindblad terms. observations, what is important for δT /T 0 is their statistical properties on the whole observed CMB surface, rather than the (classical) evolution of δT (x)/T 0 at a certain spatial location back to the moment of generation. Given that the CMB anisotropy spectrum can be written as where the transfer function is completely fixed by the background hot big bang cosmological parameters, the initial scalar power spectrum P R (k) is all that matters. For inflationary cosmology, P R (k) is identified as the power spectrum of the comoving curvature perturbation. Likewise, the B-mode polarization power spectrum is supposed to be originated from that of the primordial tensor perturbations P h (k): And such a state is claimed to be the squeezed state. This indeed is how precisely the "classicality" is achieved by the squeezed state. Through time evolution, the uncertainty of measuring the field operators and their conjugate momenta is rather exponentially increasing: ∆φ∆π φ ∝ e 2t . That is, the probability distribution spreads over in both φ and π φ . However, at the same time, the expectation values φπ φ and π φ φ in terms of the squeezed state tend to converge to the exponentially large value and noncommutativity between φ and π φ becomes suppressed [15]. Since it is a good approximation to assign φ and π φ to certain values simultaneously, the probability distribution is a function of a set of values (φ, π φ ) called the Wigner function which is interpreted as a classical probability distribution.
Once we identify the pointer bases as the set of the squeezed states like U 0,S |0 S given by (4.15), the calculations are extremely straightforward. As we may read from (4.12), the Lindblad terms are outer products of various basis vectors (4.15) along with time-and momentumdependent coefficients. To extract the matrix elements of the Lindblad terms, say, dρ red /dτ | ab with a and b being the indices of the basis vectors, from (4.12) we simply compute Then the reduced density matrix equation is described by a 7 × 7 matrix, with each row and column being distinguished by the number of creation (or annihilation) operators acting on the vacuum state. After straight calculations and considering the contributions of the products of polarization tensors, we end up with the following evolution equation of the reduced density matrix: where with k 12 = k 1 + k 2 and likewise q ab = q a + q b , with q a and q b being the momenta of external inand out-states. Note that in all contributions the configuration is the product of two identical momentum triangles with two momenta in the environment sector and the remaining one in the system sector, (EES) × (EES).

Solution of reduced density matrix
To find explicitly the coefficients of dρ red /dτ , we first note that since e λ ij is symmetric under i ↔ j, the trace of any permutation is the same: Tr e 1 (k 1 )e 2 (k 2 )e 3 (k 3 ) = the same for any order of e 1 (k 1 ), e 2 (k 2 ), e 3 (k 3 ) .
(5.1) Moreover, from the cosine law we can write h 0 as where k 123 = k 1 + k 2 + k 3 . Here, we have arranged the terms following the power of τ , as eventually we will be interested in the moment near the end of inflation, τ → 0. Notice that h 0 is symmetric under the exchange of momenta. Integrating over τ , we find (5.5) To proceed further, we need to compute the product of the polarization tensors. From the homogeneity and isotropy of the background, without loss of generality by arranging k 1 and k 2 in such a way that k 1 is along the z-direction and k 2 is confined in the xz-plane with the polar angle θ, i.e. k 1 = (0, 0, k 1 ) , k 2 = (k 2 sin θ, 0, k 2 cos θ) , (5.6) we can explicitly compute the polarization tensors for each vector from (C.2) and (C.3). Summing over all possible combinations of polarization indices we find Then, E 00 becomes where κ 1 = k 1 /H = −k 1 τ and κ 2 = k 2 /H = −k 2 τ . Note that being in the system (environment) sector, the magnitude of κ 1 (κ 2 ) is bounded up (down) to the horizon scale, i.e. k 1 (k 2 ) = H = aH. In the last equality we have used H = −1/τ , and we have defined the integrations over κ 1 , κ 2 and cos θ as a coefficient C SE . Also notice that overall there would remain no dependence on k 1 or k 2 , but only with weak dependence on the infrared cutoff that can be inferred by counting the power of κ 1 and κ 2 . Indeed, we can perform the integral analytically to find C SE explicitly as where ε ≪ 1 is the infrared cutoff for κ 1 . We can proceed in a similar manner to find E 11 to find, with q a = −q b ≡ q, where the coefficient C E is given by, withq ≡ q/H < 1, Note that C SE and C E are related by as should be the case to maintain the traceless property of (4.19).
We pause to comment on the infrared logarithmic divergence in C SE . Even though E 00 is divergent as ε → 0, the tracelessness of the Lindblad terms, or equivalently unit trace of the density matrix, guarantees that the divergence E 00 is canceled by that in E 11 as we can see from (5.12). Indeed, when the external momentum q involved in E 11 is extremely soft, our observational apparatus does not distinguish E 11 from E 00 . That is, if the momentum involved in the one-particle to one-particle process represented by E 11 is too small well below threshold resolution to be detected by any means, it is observationally identical to the process involving no particle excitation, E 00 . In that sense, the infrared cutoff ε reveals our limit of detecting the soft graviton, which suggests that we need to sum up E 00 and E 11 with q < ε inclusively to obtain the probability of no effective detection of graviton at initial and final states. This is exactly how we obtain the finite infrared amplitude of the process in quantum electrodynamics [34], even though the probability in this case is quantum mechanical whereas that in our density matrix is classical.
Finally, E 20 component is given by which is just E 11 multiplied by e 2iqτ . Since E 11 or E 20 correspond to the matrix element with two external legs, the exponential factor e −iqτ (e iqτ ) can be understood as the phase factor of the graviton wave function for the in-(out-) state. So the phase factors cancel each other in E 11 , meanwhile they are squared in E 20 and E 02 . Thus, finally, the evolution equation of the reduced density matrix is given by with E 00 , E 11 and E 20 = E * 02 given respectively by (5.8), (5.10) and (5.13).

Decoherence rate
Now, we can discuss the non-unitary effect in the reduced density matrix. First of all, the 00 element of the reduced density matrix ρ red is given by (see Appendix D) This corresponds to the otherwise unitary evolution of the squeezed state, U 0,S |0 S , i.e. the evolution from the initial Bunch-Davies system vacuum with which we have started. That is, if we take the interaction Hamiltonian H int into account, 1. H int on one hand generates, as shown in (2.17), the unitary evolution in the system sector to produce many non-zero elements of the density matrix. We can calculate them explicitly by following similar steps presented in Appendix B but here we do not show them as they are sub-leading unitary evolution effects.
2. On the other hand, H int also generates the non-unitary evolution through the Lindblad terms, reducing the probability to maintain the pure squeezed state from unity. That means, the reduced matrix element is no longer restricted to the 00 element, but spreads out to produce other reduced density matrix elements, which are interpreted as the generation of classical probabilities for corresponding processes. Still, the trace of the reduced density matrix is unity, preserving the total classical probability. In our case, ρ red | 11 is produced and Tr ρ red | 00 + ρ red | 11 ] = 1, where the trace is taken over the momentum space and polarization, guarantees this fact. At the same time, ρ red | 20 and ρ red | 02 are also produced since they are equivalent to one-particle to one-particle process: just exchanges of the in and out states. The non-unitary evolution is of second order in H int as shown in (5.15). Hence, our calculation is valid provided that we can treat H int perturbatively.
We can further read from (5.15) the rate that the pure state evolves into the mixed state through the non-unitary Lindblad terms -the decoherence rate per unit volume for the whole system modes. By exponentiating (5.15), and noting that q is an almost vanishing momentum so (2π) 3 δ (3) (q) corresponds to the comoving volume of three-dimensional space V (see Footnote 6), we can write ρ red | 00 = exp(−V Γdτ ) where, from (5.8), Here r ≡ ∆ 2 h /∆ 2 R is the tensor-to-scalar ratio and are the amplitudes of the scalar and tensor perturbations respectively. The decoherence rate (5.16) grows as the physical volume a 3 , but is suppressed by r∆ 2 R , or equivalently H 2 /m 2 Pl . That is, the suppression of Γ is because the inflationary scale is well below the full quantum gravity regime, irrespective of the inflationary dynamics. Further, for an interval of ∆τ , we find the change in ρ red | 00 is where ∆N is the number of e-folds elapsed during ∆τ . Here, we have absorbed the factor V H −3 into the normalization of the scale factor so that only the difference between a(τ ) and a(τ + ∆τ ), viz. ∆N, is highlighted. In other words, we need for ρ red | 00 to change by e −1 , i.e. for decoherence to occur, where we have taken the central value from the Planck 2018 result ∆ 2 R ≈ 2.101 × 10 −9 [4]. For a wide range of C SE and r, typically 5 ∆N dec 10.
Finally, notice that (5.10) can be compared to the decoherence rate for the curvature perturbation R representing the reduced matrix element R q |ρ red |R q given by [27] with η ≡ǫ/(Hǫ). Here, the momentum configuration is the squeezed one where the system momentum q is much smaller than the environment one. The unity in the right-hand side indicates that the initial pure state is implemented by |R q : without interaction Hamiltonian, we expect R q |ρ red |R q = 1. In addition, we find that it has the similar structure as (5.10). That means there should be similar features in the density matrix -for example the q −3 dependence implies that taking the trace of the density matrix over one particle states results in a logarithmic infrared divergence, which should be cancelled by other diagonal density matrix elements, e.g. ρ red | 00 . On the other hand, there are two differences: 1) we have an additional delta function δ (3) (q ab ) representing that in-and out-state have the same momentum. This is a result of the contraction between creation (annihilation) operator in H int of L 1 and annihilation (creation) one in Hamiltonians of L 2 , and 2) the appearance of the slow-roll parameters as expected. This reflects that the Goldstone boson nature of the curvature perturbation resulting from the spontaneous breaking of dS isometry [35]. As a result, the number of e-folds for decoherence to occur has an additional slow-roll suppression to give ∆N dec typically two times bigger than (5.19), ∆N dec 10.

Conclusions
In this article, we have considered in the dS background the time evolution of the squeezed state for tensor perturbations under the influence of non-linear interaction Hamiltonian. From this, we have shown how tracing out the unobservable sub-horizon modes converts the pure state, i.e. the initial squeezed state, into the mixed state through the Lindblad operators from the non-linear cubic interactions, in which super-(system) and sub-(environment) horizon modes interact with each other. Since the mixed state is interpreted as an ensemble of the pure states with classical probabilities represented by the density matrix elements, the process we have described is interpreted as a quantum-to-classical transition that breaks the unitarity. This needs to be distinguished from the traditional quantum-to-classical transition accomplished by the squeezed state. That is, the squeezed state can be understood as the time evolution of the vacuum state by the unitary operator, coming from the quadratic sector in which the interaction between super-and sub-horizon modes is not taken into account. In this case, the classicality is coming from suppression of non-commutativity between the field variables and their conjugate momenta. The probability of maintaining the squeezed state for the tensor perturbations is given by (5.15). The fact that the reduced density matrix has the unit trace is reflected in the generation of other matrix elements, say, 11 and equivalent 20 and 02 elements in our case. While we have found that only the processes in which two particles are involved -one-to one-particle, zero-to two-particle and two-to zero-particle -considering higher order in the interaction Hamiltonian should excite the processes with more particles involved. Furthermore, the rate of decoherence, the non-unitary quantum-to-classical transition, can be computed explicitly.
This can be read from the rate of change in the 00 element of the reduced density matrix element and is given by (5.16). Typically, it takes 5 ∆N dec 10 e-folds for decoherence to occur rapidly. This delay is more or less half of that for the curvature perturbation, because it accompanies one more slow-roll suppression due to its Goldstone boson nature, so that the decoherence occurs slower accordingly.

A Quadratic evolution of tensor perturbations A.1 Solutions of operators
Considering the spatial metric as 4 where h ij is the pure tensor perturbations with h i i = 0 and ∂ i h i j = 0, the quadratic action of h ij becomes Here, dτ ≡ dt/a is the conformal time, and the indices of h ij are raised and lowered by δ ij so we do not sharply distinguish the upper and lower indices of h ij , e.g. h ′ ij 2 should be understood as h ij h ′ ij . Since there are two physical degrees of freedom for h ij , we introduce the polarization tensor e ij (λ), with λ being the polarization index, which satisfies 5 so that we can write so that the quadratic action becomes Thus there are two copies of the identical action of a canonically normalized scalar field v λ for each polarization state. So from now on we just consider only one polarization state and drop the subscript λ, since for the other polarization we can follow precisely the same steps to describe the physical evolution of the system. Now from (A.6) we define the conjugate momentum π as usual: Then in terms of the Fourier mode as and the same for π(x), the Hamiltonian H = d 3 x πv ′ − L becomes Since we are interested in the evolution of the state, we may work in the Schrödinger picture where the time evolution is taken by the state and the operators remain fixed. For the study of cosmological perturbations usually the Heisenberg picture is taken where the operators are time evolving while the state is fixed, usually the vacuum. But in fact, it does not matter sharply in which picture we are working, since the physical results should be the same, especially at quadratic level. To proceed, we first promote the Fourier modes π k and v k to the time-dependent operatorsπ k (τ ) andv k (τ ). In terms of the creation and annihilation operatorsπ k (τ ) andv k (τ ) are given byπ (A.10) Note that with the operators being time-dependent, we are in the Heisenberg picture. One important difference from the conventional approach in the Heisenberg picture is that, the time dependence is given to the creation and annihilation operators, a k = a k (τ ) and a † k = a † k (τ ). Instead, the mode functions are set up at an initial time, τ = τ 0 . Meanwhile, in the standard approach the creation and annihilation operators are defined initially and the mode functions carry the time dependence. This difference from the standard approach is to make the transfer from the Heisenberg picture to the Schrödinger picture more easily.
Being canonical, we have the following canonical commutation relations: Notice that these relations are non-zero only for the same polarization. If we revive the polarization index, we have a λ k (τ ), a λ ′ † q (τ ) = (2π) 3 δ λλ ′ δ (3) (k − q). Using the decompositions (A.10) and the commutation relations impose that, at a initial time τ 0 , the mode functions v k and u k satisfy To set up the initial conditions at τ = τ 0 , we assume that the modes are deep inside the horizon, or the expansion of the universe is neglected. That is, we may only consider the decoupled first two terms in (A.9), say H 0 , to set up the initial conditions. Then, we can determine the initial mode function solutions by demanding that the expectation value of the free Hamiltonian operator H 0 is minimized at τ = τ 0 : where the expectation value is taken with respect to the initial vacuum state, which is assumed to match the standard vacuum state in Minkowski space as the expansion of the universe can be ignored. We can easily find the solutions: 14) and the canonical variables (A.10) becomê Substituting the solutions (A.15) into (A.9), from the Hamiltonian equations for the operators we can find the following coupled differential equations that a k (τ ) and a † −k (τ ) should satisfy: The general solution of these equations are given by the linear combination of the initial ones: which is the so-called Bogoliubov transformation. The commutation relations (A.11) lead to the constraint α k (τ ) and β k (τ ) are always subject to: and we can parametrize them in terms of the hyperbolic functions as where without loss of generality we take the time-dependent functions r k (τ ), ϕ k (τ ) and θ k (τ ) real. Assuming a perfect de Sitter background so that a ′ /a = −1/τ , we can find analytically the solutions as [13,15]

A.2 Evolution of the vacuum state
Having found the canonical variables in terms of the time-dependent creation and annihilation operators, now we can write the Hamiltonian density operator as Here, we have changed the order of a k a † k using the commutation relation of a k . The reason is because after all we are interested in the evolution of the initial vacuum state. So the annihilation operator on the right first annihilate the vacuum state, so that the operation of this part is to keep the vacuum state intact as we will see right below. For clarity, now let us concentrate on a single k-mode. This means now the Hamiltonian we need to consider is H, not necessarily H = d 3 k/(2π) 3 H. Equivalently, we can interpret this as "discretization" of the momentum, d 3 k/(2π) 3 → L −3 k where L 3 is the volume into which the mode of our interest permeates. With a given volume L 3 , we may now isolate the volume from the canonical creation and annihilation operators (A.11) in such a way that, since a k has a mass dimension of M −3/2 , and likewise for a † k . That means, the new (dimensionless) operatorsâ k andâ † k now satisfies 6 Then the relation between H and H is also clear. From (A.23), we can see So H is the Hamiltonian density in the sense that it only accounts for the contribution of a single k-mode, but it includes the whole volume over which the quantum field is effective. Yet, H is the Hamiltonian per unit volume, but the contribution of all k-modes is considered. So, in the conventional sense that density means something divided by the total volume, H is Hamiltonian "density" because the effect of the total volume is singled out. To concentrate on a single k-mode [because we have already taken into account all delta functions in the commutation relations (A.11)] let us write H = k H k where R k and S k are called rotation and squeezing operators respectively. Since the evolution operator (A.28) is coming from the quadratic Hamiltonian, we can dump the whole evolution into the (initial vacuum) state in the Schrödinger picture. Since (A.28) is factorized into (A.29) and (A.30), the action of each piece is also separated. For the rotation operator (A.29), as we expand the exponential, the operator terms annihilate the vacuum since it is first multiplied by the annihilation operatorsâ k andâ −k . Thus only the first term, which originally accounts for the volume over which the field permeates, survives. That is, with the vacuum state defined in the standard manner as 7 Thus the rotation operator R k does on the initial vacuum state |0 0 nothing but producing an irrelevant phase θ k . This is natural, since R k is coming from the decoupled, free part of the Hamiltonian H 0 for which the expansion of the universe is neglected, such that H 0 should do in the same manner as what is done in the Minkowski space. That is, the vacuum remains vacuum. The squeezing operator (A.30), we can rewrite in as exp e −2iϕ k tanh r k 2 â −kâk .
(A.33) Using a x = 1 + x log a + x 2 (log a) 2 /2! + · · · and again expanding the exponential, we find the operator termâ −kâk does nothing as it annihilates the vacuum, so that is the two-mode (with momenta k and −k for momentum conservation) state with the same occupation number n, i.e. n-particle state. Thus the squeezing operator S k , originated from the interacting part of the Hamiltonian H i , is responsible for the creation of two quanta with momenta k and −k. In other words, the cosmological vacuum fluctuations are amplified.

B Full reduced density matrix
In this section we present the matrix for the Lindblad terms before any further consideration based on the configuration of momenta. The Lindblad terms for the evolution of the reduced density matrix may well be represented by a 7×7 matrix, depending on how many operators are involved, running from 0 to 6 as: where each term is given as follows: 1. E 00 : This is the simplest case where both the bra-and ket-basis vectors a| and |b are the vacuum-evolved state. To obtain such terms, we multiply from the left U 0,S |0 S † = 0| S U † 0,S and from the right U 0,S |0 S . To have non-zero contributions, there must be the factor of the form U 0,S |0 S 0| S U † 0,S . Two contributions are coming from the case, using the notations in (4.9), H = a † 3 a † 2 a 1 , a † 2 a 1 a † 3 and a 1 a † 3 a † 2 and G = a 4 a † 6 a † 5 , a † 5 a 4 a † 6 and a † 6 a † 5 a 4 . If either H † or G has one selfcontraction, which the product of one squeezed triangle with all the momenta in the system sector and another squeezed triangle with squeezed momenta in the environment sector, sharing a common momentum q ≈ 0: (SSS) sq × (EES) sq . Also note that the same momenta mean the same polarization indices, e.g. h 1 (k 1 , q, −k 1 )g * 1 (k 4 , −q, −k 4 ) means λ 1 = λ 3 , λ 2 = λ 5 and λ 4 = λ 6 . If both H † and G have one self-contraction each at the same time, which is a product of two squeezed triangles with two momenta in the environment sector respectively, sharing a common momentum q ≈ 0: (EES) sq × (EES) sq .
The other two contributions are coming from the case H = a † 3 a † 2 a † 1 and G = a † 6 a † 5 a † 4 . If there is one cross-contraction between H † and G, which is a product of two identical triangles with one momentum in the environment sector respectively, i.e. k 3 = −k 12 + q ∈ k S with q ≈ 0: (ESS) × (ESS). If there are two cross-contractions between H † and G, which is a product of two identical triangles with two momenta in the environment sector respectively, i.e. k 3 = −k 12 + q ∈ k E : (EES) × (EES).
2. E 02 : This corresponds to the case where the bra-basis vector a| is the vacuum-evolved state and the ket-basis vector |b corresponds to the vacuum-evolved two-excited state with momenta k a and k b . To obtain such terms, we multiply from the left U 0,S |0 S † = 0| S U † 0,S and from the right U 0,S a † a a † b |0 S . To have non-zero contributions, there must be the factor of the form U 0,S |0 S 0| S a a ′ a b ′ U † 0,S . Two contributions are coming from the case H = a † 3 a † 2 a 1 , a † 2 a 1 a † 3 and a 1 a † 3 a † 2 and G = a † 4 a 5 a 6 , a 6 a † 4 a 5 and a 5 a 6 a † 4 . If either H † or G has one self-contraction, Here, in the first (second) four terms, the two momenta sets constitute independent squeezed triangles (q a ≈ 0 and q b ≈ 0 but in general q a = q b ), in one of which all in the system sector and the other two in the environment sector for k 1 and k 4 (k 4 and k 1 ) respectively: (SSS) sq × (EES) sq . Meanwhile, the last two terms share the squeezed momentum denoted in common by q ab ∈ k S : (EES) sq × (EES) sq . If both H † and G have one self-contraction each at the same time, which is the product of two independent squeezed triangles with two momenta in the environment sector respectively: (EES) sq × (EES) sq .
The other three contributions are coming from the case H = a † 3 a † 2 a † 1 and G = a 4 a † 6 a † 5 , a † 5 a 4 a † 6 and a † 6 a † 5 a 4 . If there is one self-contraction in G, which is the product of two squeezed triangles with two momenta in the system sector in one triangle and in the environment sector in the other triangle, respectively, sharing q ab ≈ 0 (and q ab ∈ k S ): (SSS) sq × (EES) sq . If there is one cross-contraction between H † and G, which is the product of two identical triangles with one momentum being in the environment sector while the other two in the system sector and with k 1 − q b ∈ k S : (ESS) × (ESS). If there are two cross-contractions between H † and G, which is the product of two triangles sharing one common momentum k 1 in the environment sector with k 1 − q b ∈ k E : (EES) × (EES).
3. E 04 : This corresponds to the case where the bra-basis vector a| is the vacuum-evolved state and the ket-basis vector |b corresponds to the vacuum-evolved 4-excited state with momenta k a , k b , k c and k d . To obtain such terms, we multiply from the left U 0,S |0 S † = 0| S U † 0,S and from the right U 0,S a † a a † b a † c a † d |0 S . To have non-zero contributions, there must be the factor of the form U 0,S |0 S 0| S a a ′ a b ′ a c ′ a d ′ U † 0,S . One contribution is coming from the case H = a † 3 a † 2 a 1 , a † 2 a 1 a † 3 and a 1 a † 3 a † 2 and G = a 4 a 5 a 6 . Only one self-contraction in H † is possible to give which is the product of a squeezed triangle with two unsqueezed momenta in the environment sector and an arbitrary triangle in the system sector, not sharing any common momentum: (EES) sq × (SSS).
The other two contributions are coming from the case H = a † 3 a † 2 a † 1 and G = a † 4 a 5 a 6 , a 6 a † 4 a 5 and a 5 a 6 a † 4 . If G has one self-contraction, which is the product of a squeezed triangle with two unsqueezed momenta in the environment sector and an arbitrary triangle in the system sector, not sharing any common momentum: (EES) sq × (SSS). If there is one cross-contraction between H † and G, + 23 perm in arranging q a , q b , q c , q d in h 0 and g 1 , which is the product of two triangles sharing a common momentum q ab ∈ k E , yet the other momenta being in the system sector: (ESS) × (ESS).
4. E 11 : This corresponds to the case where both the bra-and ket-basis vectors a| and |b are the vacuum-evolved one-excited state. To obtain such terms, we multiply from the left U 0,S a † a |0 S † = 0| S a a U † 0,S and from the right U 0,S a † b |0 S . To have non-zero contributions, there must be the factor of the form U 0,S a † a ′ |0 S 0| S a b ′ U † 0,S . Two contributions are coming from the case H = a † 3 a † 2 a 1 , a † 2 a 1 a † 3 and a 1 a † 3 a † 2 and G = a 4 a † 6 a † 5 , a † 5 a 4 a † 6 and a † 6 a † 5 a 4 . If either H † or G has one self-contraction, This is the product of a squeezed triangle with all momenta in the system sector and another, independent squeezed triangle with two momenta in the environment sector: (SSS) sq × (EES) sq . If both H † and G have one self-contraction each at the same time, which is the product of two independent squeezed triangles with two momenta in the environment sector respectively: (EES) sq × (EES) sq .
The other one contribution is coming from the case H = a † 3 a † 2 a † 1 and G = a † 6 a † 5 a † 4 . If there are two cross-contractions between H † and G, This is the product of 2 identical (because q a = q b due to the delta function) triangles with k 1 ∈ k E and k 1 − q b ∈ k E : (EES) sq × (EES) sq .
5. E 13 : This corresponds to the case where the bra-basis vector a| is the vacuum-evolved single-excited state and the ket-basis vector |b corresponds to the vacuum-evolved 3excited state. To obtain such terms, we multiply from the left U 0,S a † a |0 S † = 0| S a a U † 0,S and from the right U 0,S a † b a † c a † d |0 S . To have non-zero contributions, there must be the factor of the form U 0,S a † a ′ |0 S 0| S a b ′ a c ′ a d ′ U † 0,S . One contribution is coming from the case H = a † 3 a † 2 a † 1 and G = a 4 a † 6 a † 5 , a † 5 a 4 a † 6 and a † 6 a † 5 a 4 . Only one self-contraction in G is possible to give which is the product of one squeezed triangle with two momenta in the environment sector and another independent triangle with an arbitrary shape in the system sector: (EES) sq × (SSS).
6. E 22 : This corresponds to the case where both the bra-and ket-basis vectors a| and |b are the vacuum-evolved two-excited states. To obtain such terms, we multiply from the left U 0,S a † a a † b |0 S † = 0| S a b a a U † 0,S and from the right U 0,S a † c a † d |0 S . To have non-zero contributions, there must be the factor of the form One contribution is coming from the case H = a † 3 a † 2 a † 1 and G = a † 6 a † 5 a † 4 . If there is one cross-contraction between H † and G, where q ab = −q cd ∈ k E . This is the product of two triangles sharing a common momentum q ab = −q cd in the environment sector, while the rest momenta in the system sector: The other non-zero elements in (B.1) not explicitly presented, E 20 , E 31 and E 40 , are the counterparts of E 02 , E 13 and E 04 respectively, and can be computed almost identically.

C Reduction of density matrix elements
In the previous section, many configurations we have found are not arbitrary but specific, e.g. squeezed. This would align the momenta in a specific manner to simplify the polarization tensors. For explicit check, we adopt the unit vectors in the spherical coordinate system in such a way that the radial unit vector denotes the directionk, and the polar unit vectorê 1 and the azimuthal unit vectorê 2 denote the orientation in terms of the polar and azimuthal angles θ and ψ orthogonal tok:k = (sin θ cos ψ, sin θ sin ψ, cos θ) , e 1 = (cos θ cos ψ, cos θ sin ψ, − sin θ) , e 2 = (− sin ψ, cos ψ, 0) .

(C.3)
To check, let us consider the well-known case that k is aligned along z-direction. This corresponds to θ = ψ = 0 so thatk = (0, 0, 1). Then the polarization tensors become (also one contribution in E 20 and E 40 each). Being in the system and environment sector, we demand that the corresponding mode should be smaller and greater than the horizon scale H.
To form a triangle, we also demand that the sum of the amplitudes of the two system modes be greater than that of one environment mode. So, if k 1 and k 2 are in the system sector and k 3 is in the environment sector, the following relation should holds: where θ 12 is the angle between k 1 and k 2 . To saturate the bounds, θ 12 should lie between 0 (for k 3 = 2H with k 1 = k 2 = H) and 2π/3 (for k 3 = H with k 1 = k 2 = H). Thus, in this case the distinction between system and environment modes is not very clear, in the sense that they are different within only a factor of 2 or so. Thus, to make the difference between system and environment modes as prominent as possible to keep the validity of our effective approach of distinguishing system and environment, we take the extreme case of the folded configuration θ 12 = 0 with (k 1 , k 2 , k 3 ) = (k, k, −2k), i.e. k 3 is the environment mode. Then, all polarization tensors are transverse to all momenta: e a ij k i b = 0 for a, b = 1, 2, 3. This is obvious if we write, say, k 1 = k 2 = (0, 0, k) and k 3 = (0, 0, −2k), and multiply them to e + ij and e × ij in (C.4). Then we are only left with the terms that contain the trace of the product of three polarization tensors, e λa ij (k a )e λ b jk (k b )e λc ki (k c ). Further, fork = (0, 0, −1), θ = π and ψ = 0 so that from (C.2) and (C.3) we can explicitly find e + ij is the same as that in (C.4) but e × ij has the opposite sign to that in (C.4). Thus, for all 8 combination of polarizations, the individual trace for each polarization vanishes: Tr e + (k 1 )e + (k 2 )e + (k 3 ) = Tr e × (k 1 )e + (k 2 )e + (k 3 ) = · · · = 0 . (C.7) Thus, in the folded configuration which should be most significant for ESS, c 1 = c 2i = c 3i = 0 so that h 0 = h 1 = g 0 = g 1 = 0.

D Direct computation of reduced density matrix
We consider in the same squeezed basis (4.15) the reduced density matrix given by (2.6), where the full density matrix before reduction ρ(τ ) is given by (2.3). With ρ(τ 0 ) = |0 E |0 S 0| S | 0| E and U 0 = U 0,S U 0,E , we find (D.1) Note that the first term is just trivial: since at free level the system and environment sectors commute, i E i U 0 (τ, τ 0 )ρ(τ 0 )U † 0 (τ, τ 0 ) E i = U 0,S |0 S 0| S U † 0,S = ρ As mentioned in (2.11), this is the free-evolved density matrix of the system states ρ S (τ ). Now we turn to find the matrix element of ρ red (τ ) in the squeezed state basis (4.15). For the first term of (D.1), only the U 0,S |0 S -row and U 0,S |0 S -column is non-zero and is given by with our choice i E i |U 0,E E = 1, and other components all vanishes. For the second term of (D.1), let us first write H int,I explicitly as (4.8). Since H int,I is Hermitian, we can use H int,I = H † int,I for the second term, which will help further calculations as before. This gives 2nd term of (D.1) = −i Since the interaction Hamiltonian contains three combinations of (initial) creation and annihilation operators, 0 E G E,m (τ 0 ) 0 E and 0 E G † E,p (τ 0 ) 0 E identically vanish. Thus, there is no surviving contribution from the second term of (D.1).

(D.6)
We see that the structure of (D.6) is very similar to (4.9). Three differences are 1) the timeand momentum-dependent coefficient is the product of two g(τ )'s, 2) there is no prefactor of 1/2, and 3) the terms we have encountered during the calculation of dρ red /dτ in which ρ red is on the leftmost do not appear but they appear as the Hermitian conjugate terms. Then, we find ρ red as: where E 00 + E * 00 =