{'fraction_non_alphanumeric': 0.04822665062904327, 'fraction_numerical': 0.018343031760996067, 'mean_word_length': 4.01835311391166, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Observational constraints on quantum decoherence during inflation

Since inflationary perturbations must generically couple to all degrees of freedom present in the early Universe, it is more realistic to view these fluctuations as an open quantum system interacting with an environment. Then, on very general grounds, their evolution can be modelled with a Lindblad equation. This modified evolution leads to quantum decoherence of the system, as well as to corrections to observables such as the power spectrum of curvature fluctuations. On one hand, current cosmological observations constrain the properties of possible environments and place upper bounds on the interaction strengths. On the other hand, imposing that decoherence completes by the end of inflation implies lower bounds on the interaction strengths. Therefore, the question arises of whether successful decoherence can occur without altering the power spectrum. In this paper, we systematically identify all scenarios in which this is possible. As an illustration, we discuss the case in which the environment consists of a heavy test scalar field. We show that this realises the very peculiar configuration where the correction to the power spectrum is quasi scale invariant. In that case, the presence of the environment improves the fit to the data for some inflationary models but deteriorates it for others. This clearly demonstrates that decoherence is not only of theoretical importance but can also be crucial for astrophysical observations.


Introduction
One of the deepest insights of modern cosmology is that all structures in our Universe (galaxies, clusters of galaxies, Cosmic Microwave Background -CMB -anisotropies etc.) originate from vacuum quantum fluctuations stretched by the expansion and amplified by gravitational instability [1][2][3][4][5][6] during an early epoch of accelerated expansion named inflation [7][8][9][10][11][12]. This idea is strongly supported by the data [13,14], in particular by the fact that we observe an almost scale invariant power spectrum of curvature perturbations.
However, the quantum origins of the perturbations also raise new issues. Clearly, the structures observed today are classical objects and the question of how the quantumto-classical transition occurred in a cosmological context remains unanswered and has been the subject of many investigations . It is widely believed that decoherence [37][38][39] could have played an important role in that process [40][41][42][43][44][45][46][47][48][49][50][51][52]. On more general grounds, cosmological fluctuations have to couple (at least gravitationally) to the other degrees of freedom present in the Universe. They should thus be treated as an open quantum system rather than an isolated one, and studying decoherence is an effective way to investigate the role played by these additional degrees of freedom. As a consequence, decoherence is important not only for theoretical considerations but also for observational reasons. In other words, even if one denies its relevance in foundational issues of quantum mechanics, on the practical side, it must be taken into account since the presence of these interacting extra degrees of freedom appears to be unavoidable.
Under some very general conditions (to be discussed in this article), the evolution of an open quantum system can be modelled through a Lindblad equation, which describes how the interaction with the environment modifies the evolution of the system. This modification is such that the off-diagonal terms of the system density matrix go to zero in a preferred basis selected by the form of the interaction. Although this does not solve the quantum measurement problem [53,54] since it does not explain how a definite outcome is obtained, it explains how a preferred basis can be selected and why we do not see superpositions of macroscopic objects.
Decoherence is not an instantaneous process but proceeds over a finite time scale, controlled by the strength of the interaction between the system and its environment. Over the duration of that interaction, the environment can also change the diagonal elements of the density operator, i.e. the probabilities associated to the possible final results. This is why in general, the environment does not only suppress interferences between possible outcomes of a measurement, it also changes their predicted probabilities of occurrence.
In a cosmological setting, this means that if decoherence occurs in the early Universe, the statistical properties of cosmological perturbations may be modified. Since the later are very well constrained, in particular by measurements of the CMB temperature and polarisation anisotropies [13,14], this opens up the possibility to observationally constrain cosmic decoherence. The question we investigate in the present work is therefore the following: Which interactions with which environments allow sufficient decoherence to take place in the early Universe while preserving the standard statistical properties of primordial cosmological fluctuations as predicted from inflation and confirmed by observations?
This article is organised as follows. In Sec. 2, we present the general Lindblad equation formalism for cosmological perturbations during inflation. We pay special attention to how correlation functions of the environment enter this equation, and to the condi-tions under which it is valid. We also relate its parameters to microphysical quantities in the case where the environment is made of a heavy scalar field. In Sec. 3, we calculate environment induced corrections to the power spectrum of cosmological perturbations in the case where the interaction is linear in the Mukhanov-Sasaki variable (see Sec. 3.1, where an exact solution for the full density matrix is obtained), and in the case where the interaction is quadratic (see Sec. 3.2, where mode coupling renders the analysis more involved but the modified power spectrum can still be calculated exactly). Requiring that quasi scale invariance is preserved, we derive an upper bound on the effective interaction strength. In Sec. 4, we calculate the level of decoherence that cosmological perturbations undergo during inflation, in the case of linear interactions in Sec. 4.1 and for quadratic interactions in Sec. 4.2. Requiring that decoherence has occurred by the end of inflation, we derive a lower bound on the effective interaction strength that we compare with the above-mentioned upper bound. This allows us to identify the class of viable scenarios. In Sec. 5, we generalise our approach to arbitrary order interactions (i.e. cubic and beyond), and we show how the power spectrum and the decoherence rate can still be calculated exactly. In Sec. 6, we summarise our results and conclude by mentioning a few possible extensions. Finally, the paper ends with a series of technical appendices. Appendix A presents a detailed derivation of the Lindblad equation, focusing on the assumptions needed to obtain it, in order to determine whether and when they are satisfied in a cosmological context. In Appendix B, we discuss the case where the environment is made of a heavy scalar field, which allows us to relate the phenomenological parameters appearing in the Lindblad equation to microphysical quantities. In Appendix C, we explain how an exact solution to the Lindblad equation can be found if the system is linearly coupled to the environment, and Appendix D provides details for the calculation of the power spectrum if the coupling is quadratic.

Lindblad equation for inflationary perturbations in interaction with an environment
During inflation curvature perturbations are described by the Mukhanov-Sasaki variable [1,55] v (η, x), where η denotes conformal time and x is the conformal spatial coordinate. This variable is a combination of the perturbed inflaton field and of the Bardeen potential, the latter being a generalisation of the gravitational Newtonian potential [56]. The free evolution Hamiltonian for cosmological perturbations,Ĥ v , can be expressed as [57]Ĥ andp k =v k is the Fourier transform of its conjugate momentum. Here a prime denotes a derivative with respect to conformal time η. The Hamiltonian (2.1) describes a collection of parametric oscillators, the time-dependent frequency of which is given by where a is the Friedman-Lemaître-Robertson-Walker scale factor, 1 = 1 − H /H 2 the first slow-roll parameter and H = a /a = aH, H being the Hubble parameter (not to be confused with the Hamiltonian).

Quantising inflationary perturbations
Sincev(η, x) is real, one hasv −k =v † k . Decomposingv k andp k into real and imaginary parts according tov k = v R k + iv I k / √ 2 andp k = p R k + ip I k / √ 2 , this gives rise tô v R −k =v R k andv I −k = −v I k , and similar expressions forp R,I ±k . This shows that not all v k 's are independent degrees of freedom and that only the variables v R k and v I k for k ∈ R 3+ must be quantised, i.e. k runs on half of the Fourier space. This is done through the canonical commutation relations [v k ,p q ] = iδ (k + q) , (2.4) which also imply that v † k ,p q = v k ,p † q = iδ (k − q). The quantum state of the perturbations is represented by a wavefunctional Ψ[v(η, x)]. In Fourier space, it reads Ψ [v(η, x)] = Ψ {v s k (η)} k∈R 3+ , s∈{R,I} (2.5) i.e. it is a function of the infinite number of Fourier components in R 3+ . Since the free Hamiltonian (2.1) can be written as a sum of independent Hamiltonians on each component of the Fock space,Ĥ v = R 3+ d 3 k s=R,IĤ s k , withĤ s k = (p s k ) 2 /2+ω 2 (v s k ) 2 /2, if the wavefunctional can be factorised initially it remains so at later times, (2.6) Here, the dependence in η has been dropped in the two latest expressions for display convenience, and the product has to be understood as a tensorial one (sometimes noted ⊗). As will be shown below, in the presence of non-linear interaction (i.e. in the presence of mode coupling), this factorisation is no longer possible and one has to work with Eq. (2.5).
In order to include non-pure states in the analysis, one has to work in terms of the density matrixρ v = |Ψ [v] Ψ [v]|. In the free theory (2.1), the factorisation (2.6) gives rise to a similar one for the density matrix, (2.7) As mentioned above, when non-linear interactions are introduced, this no longer holds. The evolution of the system is controlled by the Schrödinger equation d |Ψ [v] /dt = −iĤ v |Ψ [v] or, equivalently, by the Liouville-von Neumann equation If the state is factorisable, this can also be written in Fourier space as follows. Time differentiating Eq. (2.7), one first has Then, using the fact thatĤ v = R 3+ d 3 k s=R,IĤ s k , the commutator in Eq. (2.8) can be expressed as This confirms that, in the absence of non-linear interactions, each Fourier subspace can be treated independently from the others.

Including the interaction with an environment
The previous considerations assume that the cosmological perturbations can be modelled as an isolated system. In practice however, there are other degrees of freedom in the Universe that, on generic grounds, interact with the perturbations. This is why cosmological perturbations, here described by the set of variables v k (η) or, equivalently, v(η, x), should rather be modelled as an open system interacting with the "environment" comprising all other degrees of freedom associated with other fields, cosmological perturbations outside our causal horizon, physics beyond the UV or IR cutoffs of the theory, etc. The total Hamiltonian can be written aŝ whereĤ v is the Hamiltonian (2.1),Ĥ env is the free evolution Hamiltonian for the environment that we will not need to specify, g is a dimensionless coupling constant andĤ int is the interaction Hamiltonian. Requiring that the system and the environment couple through local interactions only, it can be expressed aŝ H int (η) = d 3 xÂ (η, x) ⊗R (η, x) , (2.12) whereÂ belongs to the system sector andR belongs to the environment sector.

Lindblad equation
Let us notice that the interaction term (2.12) is, strictly speaking, not of the form usually required to derive a Lindblad equation. However, in Appendix A, we show that the usual treatment can be generalised to an interaction term of the form (2.12). In that Appendix, it is shown that, even if the full system starts off being described by a density matrixρ for which there are no initial correlations between the system and the environment,ρ (η in ) =ρ v (η in ) ⊗ρ env (η in ), then, at latter times, the system and the environment become entangled. Since we are only interested in tracking the evolution of the cosmological perturbations, let us introduce the reduced density matrix where the environment degrees of freedom have been traced out. Under the assumption that the autocorrelation time ofR in the environment, η c , is much shorter than the time scale over which the system evolves, in Appendix A the reduced density matrix (2.13) is shown to follow the non-unitary evolution equation where C R is the same-time correlation function ofR in the environment, C R (x, y) = R (η, x)R(η, y) , and the coefficient γ is related to the coupling constant g and to the autocorrelation (conformal) time η c ofR in the environment, according to γ = 2g 2 η c . (2.15) This parameter is, in general, time-dependent, and in what follows we assume that it is given by a power law in the scale factor γ = γ * a a * p , (2.16) where p is a free index and a star refers to a reference time. For convenience, we take it to be the time when the pivot scale k * = 0.05 Mpc −1 crosses the Hubble radius during inflation. The time dependence of γ comes from the fact that g and η c are found to depend on time when it comes to concrete models as will be exemplified below. The Lindblad operator (2.14) is also a generator of all quantum dynamical semigroups [58], i.e. of the transformations F η (ρ) of the density matrix indexed by the time parameter η that satisfy the Markovian property F η [F η (ρ)] = F η+η (ρ). It therefore allows one to investigate the dynamics of observable cosmological fluctuations as an open quantum system on very generic grounds.

Correlation function of the environment
If the environment is in a statistically homogeneous configuration, C R depends on x − y only, and if statistical isotropy is satisfied too, it simply depends on |x − y|. Assuming also that a single physical length scale E is involved, it is a function of a|x − y|/ E , and in what follows, for simplicity, we assume it to be a top-hat function where Θ(x) is 1 if x < 1 and 0 otherwise andC R is a constant. Notice that the appearance of the scale factor a in the argument of the top-hat function is due to the fact that x and y are comoving coordinates while the correlation length E is a physical scale. With a given model for the environment, the correlation function can in principle be calculated more precisely. In the following, the form of C R (x, y) will be left unspecified as much as possible. In any case, one expects the above approach to be a good approximation, since only when physical scales are of the order of the correlation length of the environment can our modelling be slightly inaccurate.

A heavy scalar field as the environment
In Appendix B, the case where the environment consists of a scalar field ψ with mass M H is investigated. In practice, the coupling betweenv andψ is assumed to be of the formĤ where µ is a fixed mass scale parameter, √ −g is the square root of minus the determinant of the metric andφ =v/a. The correlation function ofψ m can be calculated using renormalisation techniques for heavy scalar fields on de-Sitter space-times [59][60][61] and one findsC (2. 19) In this expression, σ(m) is 1 or 0 depending on whether m is even or odd, and "!!" denotes the double factorial. Defining the Lindblad operator aŝ the ansatz (2.16) is realised with to be confused with the determinant of the metric) is time dependent and effectively scales as g ∝ λa 4−n . Since the correlation cosmic time of the environment t c is constant, η c = t c /a scales as the inverse of the scale factor. 1 We conclude, using Eq. (2.15), that γ ∝ g 2 η c ∝ a 7−2n . Finally, since H ∝ a − 1 at first order in slow roll,C R is not strictly constant but scales as a −6m 1 . This slow time dependence can be absorbed in the definition of γ, by shifting p → p − 6m 1 * , and in replacing H with H * in Eq. (2.19). For linear interactions, n = 1, p 5, for quadratic interactions, n = 2, p 3, and in Sec. 3 we will show why these behaviours are in fact very remarkable.
The typical time scale over which the system evolves is of order the Hubble time, so the assumption that it is much longer than the environment autocorrelation time, which is necessary in order to derive the Lindblad equation, amounts to M H. In Appendix B, it is shown that this condition also guarantees that ψ is a test field (i.e. does not substantially contribute to the energy budget of the Universe). Another requirement for the validity of the Lindblad approach is that the interaction between the system and the environment does not affect much the behaviour of the environment and only perturbatively affects the system. In Appendix B this is shown to be the case if λ H 6−3m M 2m−2 /(µ 4−n−m φ n ). One concludes that, if the additional field ψ is sufficiently heavy, and if the coupling constant λ is sufficiently small, the influence of ψ on the dynamics of the Mukhanov-Sasaki variable associated with φ can be studied with the Lindblad equation (2.14), with the parameters given in Eqs. (2.16)-(2.21).

Evolving quantum mean values
In order to extract observable predictions from the quantum state described by the density matrixρ v , quantum expectation values Ô = Tr ρ vÔ (2.22) have to be calculated, whereÔ is an arbitrary operator acting in the Hilbert space of the system. When the Lindblad equation (2.14) cannot be fully solved, it may also be convenient to restrict the analysis to (a subset of) such expectation values, as will be shown below. Differentiating Eq. (2.22) with respect to time and plugging Eq. (2.14) in leads to (2.23) In this expression, ∂Ô/∂η accounts for a possible explicit time dependence of the operator O. The term describing the interaction between the system and the environment can be written in Fourier space, and one obtains (2.24) To derive this expression, we have assumed that the environment is placed in a statistically homogeneous configuration such that, as stated above, C R (x, y) = C R (x − y), and the functions C R (x − y) andÂ(x) have been Fourier expanded in a similar way as in Eq. (2.2). Using the top-hat ansatz (2.17), one obtains  26) where the amplitude at the origin has been matched.

Power spectrum
In the previous section we have shown how interactions with the environment can be modelled through the addition of a non-unitary term in the evolution equation of the density matrix for the system, which becomes of the Lindblad type (2.14). As explained in Sec. 1, this new term leads to the dynamical suppression of the off-diagonal elements of the density matrix when written in the basis of the eigenstates of the operator through which the system couples to the environment [Â in the notations of Eq. (2.14), which here we take to be some power of the Mukhanov-Sasaki variablev, see Eq. (2.20)]. This will be explicitly shown in Sec. 4 and is at the basis of the "decoherence" mechanism. However, Eq. (2.14) also has a unitary term which comes from the free Hamiltonian of the system. If that term couples the evolution of the diagonal elements of the density matrix to the non-diagonal ones, the Lindblad term also induces modifications of the diagonal elements of the density matrix, hence of the expected probabilities of observing given values ofv, that is to say of the observable predictions for measurements of the system. In particular, the power spectrum of cosmological curvature perturbations is altered by the presence of the Lindblad term, and in this section we calculate the corrected power spectrum. We then determine how observations constrain the size of this correction, and thus place bounds on the strength of interactions with the environment.

Linear interaction
Let us first consider the case where the system couples linearly to the environment, A(x) =v(x), i.e. n = 1 in Eq. (2.20). We will show that the Lindblad equation can be solved exactly in that case, i.e. all elements of the density matrix will be given explicitly.
In particular, we will show how the power spectrum can be extracted and how the correction it receives from the Lindblad term can be studied.

Lindblad equation in Fourier space
Let us first show that if n = 1, the Lindblad equation (2.14) decouples into a set of independent Lindblad equations in each Fourier subspace. Since we have shown that this is already the case in the free theory, see Eq. (2.10), it is enough to consider the interaction term only and to show that it has the same property. From Eq. (2.14), the interaction term is given by 1) where in the first equality we have Fourier expanded C R andv and in the second equality we have used the decompositionv p = v R p + iv I p / √ 2 introduced in Sec. 2.1. In Eq. (3.1), the two last terms vanish. Indeed, one can split the integral over R 3 into two pieces, namely over R 3+ and R 3− , and, in the second piece, perform the change of variable p → −p. Using the symmetry relationv R −p =v R p andv I −p = −v I p together with the fact that the correlation function of the environment only depends on the modulus of the wavevector (hence is independent of its sign), one obtains that the two pieces of the integral cancel out each other. If the state is factorisable as in Eq. (2.7), one then finds that the interaction term can be similarly factorised, The fact that the interaction term is linear inv(η, x) thus preserves the property that each Fourier mode evolves separately, and combining Eqs. (2.9), (2.10), and (3.2), one obtains dρ s Let us also notice that a particular comoving scale appears in the interaction term. Indeed, in order for Eq. (3.3) to have the correct dimension, γC R (k) must be homogeneous to the square of a comoving wavenumber. In what follows, we denote this scale k γ , and using the form (2.26), it can be written as , (3.7) where the quantities I k , J k and K k are defined by and v k (η) is the solution of the Mukhanov-Sasaki mode equation v k + ω 2 (k)v k = 0 with initial conditions set in the Bunch-Davies vacuum. With Eq. (2.7), this provides an exact solution to the full Lindblad equation (2.14). Because of the linearity of the interaction term, the state is still Gaussian, and one can check that it is properly normalised, Tr(ρ v ) = 1. If one sets γ = 0, i.e. if one switches off the interaction with the environment, one can also check that . The standard two-mode squeezed state, which is a pure state, is therefore recovered in that limit.
Notice that the ability to set initial conditions for the perturbations in the Bunch-Davies vacuum, one of the most attractive features of inflation, is preserved by the Lindblad equation, thanks to the presence of the environment correlation function. Indeed, as long as the mode has not crossed out the environment correlation length, it is unaffected by the presence of the environment, see Eq. (2.26).

Two-point correlation function
In the solution (3.7), one can check that the diagonal elements of the density matrix, i.e. the coefficients obtained by setting v , are affected by the presence of the environment since they involve the quantity J k defined in Eq. (3.9), which depends on γ. This confirms that the observable predictions one can draw from the state (3.7) are modified by the interaction with the environment. Since this state is still Gaussian, this modification is entirely captured by the change in the two-point correlation function, i.e. the power spectrum, that we now calculate.
The quantum mean value of (v s k ) 2 can be expressed as Making use of Eq. (3.7), this integral is Gaussian and can be performed easily, leading to In the absence of interaction with the environment, J k = 0 and one recovers the standard result. The power spectrum of curvature perturbations can be directly obtained from the relation ζ k = v k /(a √ 2 1 M Pl ), and this leads to (3.15)

Alternative derivation of the power spectrum
Before proceeding to the "concrete" calculation of the modified power spectrum, i.e. of the quantity J k /|v k | 2 , let us notice that the above result can also be obtained without solving for the Lindblad equation (2.14) entirely, but by restricting the analysis to two-point correlators. This technique will be of particular convenience in the case of quadratic interaction with the environment since there, no explicit solution to the Lindblad equation can be found, see Sec. 3.2. Let us first consider the case of one-point correlators. Making use of Eq. (2.24) withÔ =v k andÔ =p k , one has This is nothing but the Ehrenfest theorem since no correction due to the interaction with the environment appears in these equations. Combined together, they lead to v k + ω 2 (k) v k = 0, i.e. the classical Mukhanov-Sasaki equation for v k . Since the state is initially symmetric, v k = p k = 0, it remains so at all times.
If one then considers two-point correlators of the form Ô (3.17) One can see that the environment induces a new term proportional to γ only in the evolution equation for p k 1p k 2 , but since the previous equations are all coupled together, it affects the evolution of all two-point correlators. The appearance of a Dirac function δ(k 1 + k 2 ) in the modification to the last equation is important since it means that the interaction with the environment preserves statistical homogeneity, i.e. solutions of the form can be found. 2 If the environment correlator also preserves statistical isotropy,C R (k) = C R (k), as is the case of the ansatz (2.26) adopted in this work, the above system admits statistically homogeneous and isotropic solutions where P OO in Eq. (3.20) These equations can be combined into a single third-order equation for P vv only, which reads P vv + 4ω 2 P vv + 4ω ωP vv = S 1 , (3.21) with S 1 being defined as which is a function of x1 − x2 only and is invariant by translating both x1 and x2 with a constant displacement vector u. Conversely, one can show that if Ô (x1)Ô (x2) = Ô (x1 + u)Ô (x2 + u) for all u, then Ô k 1Ô k 2 = Ô k 1Ô k 2 e iu(k 1 +k 2 ) for all u, hence Eq. (3.18) must be true.
The consistency check is then to verify that the solution obtained in Eq. (3.13) is indeed a solution of Eq. (3.21). Using the explicit form of J k given by Eq. (3.9), we see that differentiating Eq. (3.13) requires to differentiate a function of the form η −∞ dη g(η , η), which gives g(η, η) + η −∞ dη ∂[g(η , η)]/∂η. This leads to the term corresponding to g(η, η) being absent since proportional to Im 2 [v k (η) v * k (η)] = 0. Similar considerations lead to the third derivative of P vv , which can be expressed as (3.24) Using the mode equation for v k , v k + ω 2 v k = 0, it is then easy to show that P vv given in Eq. (3.13) is indeed a solution of Eq. (3.21).

Power spectrum in the slow-roll approximation
Our final move is to use the slow-roll approximation to calculate ∆P k = J k /|v k | 2 , where J k is defined by Eq. (3.9) and v k is the solution to the Mukhanov-Sasaki equation that is normalised to the Bunch-Davies vacuum in the sub-Hubble limit. This is done in details in Appendix C.2 and here we simply quote the results. Defining ν ≡ 3/2 + 1 * + 2 * /2, The result can then be expanded in two limits. The first one consists in using the requirement that the environment autocorrelation time t c is much shorter than the typical time scale over which the system evolves, H −1 . If the environment correlation time and length are similar, t c ∼ E , as is the case for the example discussed in Appendix B where the degrees of freedom contained in a heavy scalar field play the role of the environment, this amounts to H E 1. The second limit consists in evaluating the power spectrum at the end of inflation, where all modes of astrophysical interest today are far outside the Hubble radius, i.e. are such that −kη 1. In these limits, the dominant contribution to the power spectrum depends on the value of p and three cases need to be distinguished. If p > 3 + (2 + 2ν)/(1 + 1 * ), referred to as case one if what follows, then . The second case is when 3 + 1/(1 + 1 * ) < p < 3 + (2 + 2ν)/(1 + 1 * ), for which the quantity ∆P k reads (3.26) Finally remains the third case where p < 3 + 1/(1 + 1 * ) and one obtains For complete consistency, these expressions must also be expanded in slow roll since we have used this approximation before. At first order, the form of the result for the three cases (i = 1, 2, 3) is given by where A i (k), B i , C i and D i can be calculated from Eqs. (3.25), (3.26) and (3.27) and are given by , (3.29) , (3.32) where γ E 0.577 is the Euler-Masheroni constant and ψ(z) is the digamma function.
Let us note that in the first and second cases, the result is independent of E , since the main contribution to the integral (3.9) comes from the neighbourhood of its upper bound. It means that the details of the shape of the environment correlator C R are irrelevant in these cases, and the top-hat ansatz (2.17) we have employed does not lead to any loss of generality. In the third case however, the result is directly dependent on E . The slow-roll corrections given in Eq. (3.36) may therefore lie beyond the accuracy level of the present calculation where the environment correlator is approximated with a top-hat function. In any case, the observational constraints derived below make use of the expression (3.35) for the amplitude A 3 only in this case, and are therefore robust.
Another remark is that in the second and third cases, the power spectrum settles to a stationary value at late time since none of the expressions (3.26) and (3.27) depends on time. In the first case however, the power spectrum is not frozen on large scales and continues to evolve as is revealed by the amplitude A 1 in Eq. (3.29), which can also be written as (3.37) In this expression, N − N * is the number of e-folds elapsed since the pivot scale crosses the Hubble radius. The exponential dependence on N − N * explains why we have not expanded this term, and the time-dependent term of Eq. (3.29), in slow roll. Let us note that the power of e N −N * is positive since the condition p > 3 + (2 + 2ν)/(1 + 1 * ) is precisely what defines case number one. The correction to the standard result is therefore enhanced by a very large factor in this case.

Observational constraints
In order to check the validity of the above calculations, we have numerically integrated the power spectrum for different values of p, and compared the result to our analytical approximations. The comparison is displayed in Fig. 1, where one can check that the analytical approximations (3.29)-(3.36) match very well the numerical result (because of the perfect overlap they can in fact not even be distinguished).
The structure of Eq. (3.14) implies that the power spectrum in the presence of decoherence is made of two branches, the standard, almost scale-invariant one, and a new branch which strongly deviates from scale invariance (except for the case p 5 discussed below). The scale k t at which the transition between the two branches occurs is such that A i (k t ) ∼ 1, and Eqs. (3.37), (3.32) and (3.35) give rise to   where ∆N * corresponds to N −N * evaluated at the end of inflation, i.e. ∆N * = N end −N * . More precisely, in the first case, i.e. when p > 8 − 3 1 * + 2 * , Eq. (3.37) shows that the correction to the standard power spectrum scales as k 3 regardless of the value of p. This can be checked on the bottom-right panel of Fig. 1 where, for p = 9, 10 and 11, the corrections grow for large values of k (small scales) and have the same slope. Since we observe an almost scale-invariant power spectrum, the part k > k t of the power spectrum has to be outside the observable window, k t k * . Together with Eq. (3.38), it gives rise to Through Eq. (3.4), this directly constrains the interaction strength with the environment, to very small values. The second case corresponds to 4 − 1 * < p < 8 − 3 1 * + 2 * , and Eq. (3.32) implies that the correction scales as k p−5 , such that it modifies the power spectrum at k > k t if p > 5 and k < k t if p < 5 (the case p = 5 is singular and is discussed separately below), with a p-dependent slope. This can be checked on the top-right and bottom-left panels of Fig. 1. For instance, for p = 4 (top-right panel), A 2 (k) ∝ k −1 and, indeed, the correction grows on large scales. On the contrary, if we consider, say p = 6 and p = 7 (bottom-right panel), one has A 2 (k) ∝ k 1 and A 2 (k) ∝ k 2 respectively, and, indeed, the corrections now grow on small scales and with different slopes. In order for observed modes to lie outside the modified, non scale-invariant part of the power spectrum, one needs to have k t k * for p < 5 and k t k * for p > 5. In both cases, with Eq. (3.38) it leads to Finally, the third case is defined by p < 4 − 1 * . The scaling of A 3 (k) is the same as for A 2 (k) as can be checked on the top-left panel in Fig. 1. For instance, p = 2 corresponds to A 3 (k) ∝ k −3 and p = 3 to A 2 (k) ∝ k −2 . Observational constraints on scale invariance impose k t k * , which, with Eq. (3.38), translates into Notice that the constraint k t k * (when p < 5) is conceptually different from the constraint k t k * (when p > 5). Indeed, in the former case, one removes the corrections to the power spectrum outside our observational window, while in the later case, one pushes the corrections to scales that are smaller than the ones probed in the CMB but that could still be of astrophysical interest.

Case of a heavy scalar field as the environment
As can be seen in Fig. 1, of particular interest is the case p = 5, see the pink line in the bottom-left panel, where the power spectrum is almost scale invariant even in the modified branch. A crucial remark is that the microphysical example studied in Appendix B, where the environment is made of the degrees of freedom contained in a heavy test scalar field, precisely gives p 5. More precisely, from Eq. (2.21), one has p = 5 − 6m 1 * in that case. Let us study the observational predictions of this model in more details. Combining the standard expression of the power spectrum, namely 7296 is a numerical constant, with Eqs. (3.14) and (3.26), one has In this expression, recall that k γ /k * is given by Eq. (3.5). This formula differs in two ways from the standard one (3.42). First, the amplitude is no longer a function of the inflationary energy scale H * and the first slow-roll parameter 1 * only, but now also depends on the ratio k γ /k * . If one assumes that tensor perturbations are not affected by the presence of the environment and that their power spectrum is still given by the standard formula P h 2H 2 * /(π 2 M 2 Pl ), the tensor-to-scalar ratio r ≡ P h /P ζ , which is the standard case in given by r| standard = 16 1 * , now reads When k γ /k * 1, one recovers the standard result, otherwise r is reduced compared to the free theory. Second, the spectral index n S ≡ 1 + d ln P ζ /d ln k is also modified, and instead of the standard formula n S | standard = 1 − 2 1 * − 2 * , we now have When k γ /k * 1, one recovers the standard result n S n S | standard but if k γ /k * 1, one obtains n S n S | standard − (6m − 2) 1 * . The shift in the spectral index ∆n S is negative (at least for m > 1/3) and proportional to 1 * , so it is still compatible with quasi scale invariance but it may have consequences for particular models of inflation.
To illustrate this result, in Fig. 2 we have displayed the spectral index and the tensor-to-scalar ratio predicted by three single-field models of inflation: Higgs inflation [63] (or the Starobinsky model [7], HI) where the potential is given by , and natural inflation [65] (NI) where V (φ) ∝ 1+cos(φ/f ). The predictions are calculated with the ASPIC library [66,67], where the free parameters appearing in the potential   Figure 2. Spectral index n S and tensor-to-scalar ratio r for three single-field models of inflation, Higgs inflation (or the Starobinsky model) "HI" and power-law inflation "PLI" in the left panel, and natural inflation "NI" in the right panel. The black lines correspond to the one-and twosigma contours obtained from the Planck 2015 T T data combined with the high-C T E + C EE likelihood and the low-temperature plus polarisation likelihood (PlanckTT,TE,EE+lowTEB in the notations of Ref. [62]) together with the BICEP2-Keck/Planck likelihood described in Ref. [14]. The blue disks correspond to the standard results in absence of decoherence, while the green, red and pink disks correspond to the modified power spectrum (3.43) with k γ /k * = 1, for m = 1, 2 and 3 respectively.
are varied together with the reheating temperature 3 according to the priors of Ref. [69], and one assumes the averaged equation-of-state parameter during reheating to vanish. The black lines are the one-and two-sigma contours of the BICEP2-Keck/Planck 2015 likelihood [14,62]. The blue disks correspond to the standard predictions (3.42), while the green, red and pink disks correspond to the modified predictions (3.43) k γ /k * = 1, for m = 1, 2 and 3 respectively. For the plateau model of Higgs inflation, because the predicted tensor-to-scalar ratio in the free theory is small, 1 * is small and the shift in the spectral index is also small, such that the model still provides a good fit to the data in the presence of decoherence. For power-law inflation, in the absence of decoherence the model is strongly disfavoured, since it predicts a too large value for n S when r is sufficiently small. In the presence of decoherence however, the decrease in the value of 3 In the case where kγ/k * does not vanish, the amplitude of the power spectrum is modified and the constraint on the inflationary energy scale H * arising from the measured amplitude [68] of the power spectrum, P ζ 2 × 10 −9 , is modified too. This implies that the upper bound on the reheating temperature is decreased and that the uncertainty coming from the reheating expansion history is reduced. This effect is however minimal in Fig. 2 where kγ/k * = 1 (it would be larger for kγ/k * 1, except for powerlaw inflation where the potential is conformally invariant and gives predictions that are independent of the reheating expansion history).
both n S and r is such that the model becomes viable for some values of its free parameter α. For natural inflation, the standard version of the model is disfavoured since it predicts a too small value for the spectral index. By further decreasing the spectral index, decoherence makes it worse and the model is even more disfavoured.
One concludes that for the interaction model proposed in Appendix B, the observational constraint on the interaction strength, here parametrised by k γ /k * , depends on the specific model of inflation. For models predicting the right value of the spectral index and very low values for the tensor-to-scalar ratio, such as Higgs inflation, decoherence does not change much the predictions and as a result there is no strong constraint on k γ /k * . For models predicting too large values for the spectral index, such as power-law inflation, the model becomes viable only in the presence of decoherence and observations place a lower bound on the interaction strength. On the contrary, for models predicting too small values for the spectral index, such as natural inflation, or models predicting values for the spectral index that are in agreement with the data and values for the tensor-to-scalar ratio that are not far from the current upper bound, decoherence can only make the model worse, and observations impose an upper bound on the interaction strength.

Quadratic interaction
Let us now consider the case where the system couples quadratically to the environment, A(x) =v 2 (x), i.e. n = 2 in Eq. (2.20). Contrary to the case of linear interactions, because of mode coupling, the Lindblad equation (2.14) does not decouple into a set of independent Lindblad equations for each Fourier mode, and cannot be solved entirely. However, the power spectrum can still be calculated by making use of the technique presented in Sec. 3.1.4.

Two-point correlation functions
Let us start by deriving the equation of motion for the two-point correlation functions. Using the fact that the Fourier transform of a squared function is the convolution product of its Fourier transform, namelŷ As explained in Sec. 3.1.4, one can derive the equations of motion for one-point correlators and, as before, find that, in agreement with Ehrenfest theorem, the standard equations are unmodified, see Eqs. (3.16). For two-point correlators, as in Sec. 3.1.4, only the equation for p k 1p k 2 is changed and one finds An important remark is that the term involving γ is not explicitly proportional to δ(k 1 + k 2 ) as was the case for linear interactions, where the presence of the Dirac function guaranteed that the Lindblad term preserved statistical homogeneity. One may therefore be concerned that the above system generates violations to statistical homogeneity. However, if the system is solved through a perturbative expansion in γ, during the first iteration the Lindblad term contains the correlator v k+k 1v −k+k 2 evaluated in the free theory, which is proportional to δ(k 1 + k 2 ). This guarantees that the solution that is obtained at the first iteration is statistically homogeneous. Since it sources the equation at the second iteration, the solution is again statistically homogeneous, so on and so forth. The result is therefore statistically homogeneous, 4 and Eq. (3.48) admits a solution of the form given by Eq. (3.18). IfC R is also isotropic,C R (k) =C R (k), the system for isotropic solutions then reduces to As was done in the linear interaction case, one can combine the above equations in order to get a third order equation for P vv only, Compared with the corresponding equation (3.21) for linear interactions, where the power spectrum is sourced by the Fourier transform of the environment correlation function, in the present case it is sourced by the convolution product of the Fourier transform of the environment correlation function with the power spectrum itself. This makes clear that, as announced, quadratic interactions yield mode coupling, since the power spectrum at a given mode is sourced by its value at all other modes.

Solving the equation for the power spectrum
The third-order differential equation for the power spectrum P vv , Eq. (3.51), is of the form where the source S 2 is a function of time that involves the power spectrum P vv itself evaluated at all scales, namely This is therefore an integro-differential equation, that couples all modes together, and that is very difficult to solve in full generality. However, at leading order in γ, one can use the free theory to calculate S 2 , which becomes a fixed function of time. Then, Eq. (3.52) is turned into an ordinary differential equation that can be solved as we now explain. If one wants to go to higher orders in γ, one can repeat the procedure and embed it in a recursive expansion in γ, but that would be inconsistent with the standard derivation of the Lindblad equation, see Appendix A, which is valid at linear order in γ only. Inspired by the fact that Eq. (3.9) provides a solution to Eq. (3.21), let us introduce the function S k defined by where v k (η) is a mode function, that is to say a solution of the Mukhanov-Sasaki equation v k + ω 2 v k = 0. Since the Wronskian of v k is conserved through the Mukhanov-Sasaki equation, the factor in front the integral in Eq. (3.54) is a constant. Using similar techniques as the ones used below Eq. (3.21), it is straightforward to check that S k is a particular solution of the equation of motion (3.52) for P vv .
In fact, one can show that this particular solution is independent of the mode function v k (η) one has chosen. Indeed, S k (η) is the solution of Eq. (3.52) that, as can be shown from Eq. (3.54), satisfies S k (η 0 ) = S k (η 0 ) = S k (η 0 ) = 0. It is therefore unique. In practice, to evaluate it, one can use the Bunch-Davies normalised mode but the result is independent of that choice, and S k (η) therefore carries a single integration constant, namely η 0 . The full solution can be obtained by adding the solution of the homogeneous equation (i.e. without the source term) P vv + 4ω 2 P vv + 4ωω P vv = 0. As already mentioned for the linear case, one can If one sets the mode function appearing in the first term of the right-hand side of Eq. (3.55) to the Bunch-Davies normalised one, P vv matches the Bunch-Davies result in the infinite past if η 0 = −∞, and this leads to

Calculation of the source term
The next step is to calculate the source term, that is to say the convolution product between the power spectrum in the free theory and the Fourier transform of the environment correlator. This is done in details in Appendix D.1 and here, we simply quote the result. By performing the angular integration, one can first show that 3). The integral over p contains a UV part (p > aH, sub-Hubble scales) and an IR part (p < aH, super-Hubble scales). The UV part is finite and subdominant at late time. In any case, as usually done, it is regularised away by adiabatic subtraction [70,71], which amounts to setting the upper bound of the integral over p in Eq. (3.57) to the comoving Hubble scale −1/η.
The IR part is divergent, and can be regularised by imposing an IR cutoff, that corresponds to the comoving Hubble scale at the onset of inflation and that we call −1/η IR . In the integral over z,C R is such that it vanishes when its argument is larger than the comoving correlation wavenumber of the environment, a/ E , and goes to a constant value in the opposite limit. Since the integral over p is now restricted to super-Hubble modes, p < aH, and given that we assumed E H −1 when deriving the Lindblad equation (at least if E ∼ t c , see Appendix A), one has p a/ E . Two cases can then be distinguished. If k a/ E , both |k − p| and k + p are much smaller than a/ E , andC R ( √ z ) is constant over its integration range. If k a/ E , both |k − p| and k + p are much larger than a/ E , andC R ( √ z ) vanishes over its full integration range. This means that in Eq. (3.57), one can simply replaceC R ( √ z ) byC R (k), and one obtains In Appendix D.1, it is shown how this formula arises as a leading order expansion in H E starting from the ansatz (2.17), here we simply gave a heuristic argument. Finally, on super-Hubble scales and neglecting slow-roll corrections, the power spectrum in the free theory is given by P vv (p) = (2p) −1 (−pη) −2 . Plugging this formula into Eq. (3.58), one is led to Notice that this expression was obtained neglecting all slow-roll corrections. Indeed, since the integral of the power spectrum over all modes appears in the source function, a slow-roll expansion around the pivot scale k * cannot be used consistently to describe the entire set of modes, and one would have to specify the details of the inflationary potential over the entire inflating domain [72] in order to calculate the source term beyond the de-Sitter limit. For this reason, we will ignore all following slow-roll corrections since their inclusion would be inconsistent. A more refined calculation would have to be carried out on a model-by-model basis. From Eqs. (3.51) and (3.59), the source function is then given by (3.60)

Power spectrum and observational constraints
The final step consists in inserting this source with Eq. , where the equivalent of Eq. (3.15) for the quadratic interaction reads ∆P k ≡ S k /|v k | 2 . One finds that three regimes have to be distinguished. The first regime corresponds to when p > 6, and gives (3.61) The second case is when 2 < p < 6, for which one has Finally, the third regime is when p < 2. In that case, the modification to the power spectrum reads (3.63) In these expressions, N − N IR ≡ ln(η IR /η) denotes the number of e-folds elapsed since the onset of inflation, and we have introduced the dimensionless coefficient that characterises the strength of the interaction with the environment. One notices that the cases p = 2 and p = 6 are singular and must be treated separately, giving rise to These analytical formulas are superimposed with a numerical integration of the power spectrum in Fig. 3, where one can check that the agreement is excellent (because of the perfect overlap, they are even hard to distinguish). As for the linear case, the power spectrum is made of two branches, one which matches the standard formula, and one where scale invariance is strongly broken (with the exception of the case p = 3 that is discussed below separately). The transition between these two branches occurs at the scale k t such that ∆P kt = 1, and expressions similar to Eq. (3.38) can be derived from Eqs. (3.61)-(3.65) (we do not write them down here for display convenience but they are straightforward). A case-by-case analysis reveals several other similarities with the linear interaction results.
In the first case indeed, i.e. for large values of p, the power spectrum is not frozen on large scales and continues to increase, leading to a very large enhancement of the correction to the standard power spectrum at late time. The modified branch of the power spectrum scales as k 3 , similarly to what was seen in Eq. (3.37). The requirement that observed modes are scale invariant, k t k * , leads to the constraint   where N T corresponds to N −N IR evaluated at the end of inflation, i.e. N T = N end −N IR and stands for the total duration of inflation. Let us notice that the 1/N T term originates from the last term in Eq. (3.61), which dominates when N T is very large. In the second case, i.e. for intermediate values of p, the power spectrum is frozen on large scales and independent of the shape and correlation length of the environment correlator. This is again similar to the linear case. The power spectrum is modified on small scales, k > k t , if p > 3, and on large scales, k < k t , if p < 3. The requirement that the modified, non-scale invariant branch of the power spectrum is unobserved leads to The exception evading this constraint is the case p = 3, where the power spectrum is almost scale invariant even in the modified branch. Strikingly, p 3 again corresponds to the scaling expected in the model proposed in Appendix B, where the environment is made of a heavy test scalar field. In this case, σ γ is related to the microphysical parameters of the model according to  (3.64). Similarly to what was done in Sec. 3.1.7, constraints on γ could be placed from the current observational bounds on n S and r within specific models of inflation, but this would require to include slow-roll corrections to the calculation of the source term, as explained above.
In the third case finally, i.e. for small values of p, the power spectrum freezes on large scales. The amplitude of the correction depends explicitly on the environment correlation length E , and it scales with the wavenumber in the same manner as for intermediate values of p. This is again similar to the linear case. The power spectrum is modified on large scales k < k t , and requiring that the non scale-invariant part is unobserved leads to the constraint . (3.69)

Decoherence
In this section, we show how the addition of a non-unitary term in the evolution equation (2.14) of the density matrix of the system, that models the interaction with environmental degrees of freedoms, leads to the dynamical suppression of its off-diagonal elements in the basis of the eigenstates of the interaction operator (here the Mukhanov-Sasaki variablev). Since this decoherence mechanism is thought to play a role in the quantum-to-classical transition of primordial cosmological perturbations, as explained in Sec. 1, we calculate the required strength of interaction that leads to decoherence at the end of inflation. We then compare this lower bound with the upper bound derived in the previous section from the requirement that quasi scale invariance is preserved. We thus identify the models for which successful decoherence occurs without spoiling standard observational predictions.

Linear interaction
In the case of linear interactions with the environment, in Sec. 3.1.2 (see also Appendix C) we have shown how Eq. (2.14) can be solved exactly and the density matrix was given in Eq. (3.7). Since it is written in the basis of the eigenstates ofv k , through which the system couples to the environment, the suppression of its off-diagonal elements directly allows us to study decoherence.

Decoherence criterion
Let us consider an off-diagonal element of the density matrixρ s k , away from the diagonal by a distance ∆v k , that we write v s where Eq. (3.13) has been used for the power spectrum P vv (k), and where we have introduced the parameter In Eq. (4.1), the factor 1/4 has been separated out from the definition of δ k since it is present even in the absence of interactions with the environment (contrary to δ k ), such that δ k characterises the additional decrease of the off-diagonal elements produced by the environment. Successful decoherence is characterised by the condition δ k 1, which can be justified by either of the three following reasons.
First, the typical distance between two realisations v s,(1) k and v s,(2) k of the Mukhanov-Sasaki variable v k is by definition given by the square root of its expected second moment, P vv (k) . If one replaces ∆v k by P vv (k) in Eq. (4.1), one can see that a large suppression of the off-diagonal element corresponds to a large value of δ k .
Second, δ k appears in Eq. (4.1) added to 1/4, which corresponds to the standard decrease of the off-diagonal term in the absence of an environment. It seems therefore natural to compare the environment-induced suppression of the off-diagonal element with the standard one, and to define "decoherence" by the requirement that the former is larger than the later, δ k 1. Third, decoherence is sometimes characterised by the so-called purity of the state, defined as Tr ρ s k 2 (for other measures of coherence see e.g. Ref. [73]). Making use of Eq. (3.7), this can be expressed as a Gaussian integral and one obtains When δ k 1, Tr ρ s k 2 1 and the state remains pure, while when δ k 1, Tr ρ s k 2 1 and the state becomes highly mixed. This is another reason why we define decoherence with the criterion δ k 1.

Calculation of the decoherence parameter
At linear order in γ, which is the order at which the Lindblad equation has been established in Appendix A, Eq. (4.2) gives rise to since one recalls that I k , J k and K k all carry a factor γ, see Eqs.  many cancellations occur (the reason behind all these cancellations will be made explicit in Sec. 4.1.5) and the following expression is obtained (4.5) The corresponding time evolution of δ k is displayed on the left panel of Fig. 4 for different values of p and a fixed value of k (k = k * ), and on the right panel for a fixed value of p (p = 3) and different values of k. On the left panel, one can see that δ k takes off very rapidly as soon as the mode under consideration crosses out the correlation length of the environment, and either settles to a stationary value afterwards (if p ≤ 2) or continues to grow (if p > 2). The case p = 3 for different values of k is displayed on the right panel, where one can see that after crossing out the environment correlation length, δ k remains stationary for some transient period of time and starts to increase again after a few e-folds. Since δ k always increases and since it takes off later for smaller scales (larger values of k), it is larger on larger scales (smaller values of k).
This behaviour can be analytically understood by computing δ k at the end of inflation, when the modes of astrophysical interest today are well outside the Hubble radius and Eq.  least if E ∼ t c ), one obtains (4.6) In this equation, one can see that if the coefficient in the argument of the exponential is positive, namely p > 3 + 2(1 − ν)/(1 + 1 ) (or, if one neglects slow-roll corrections, p 2), then δ k grows on large scales. If, on the contrary, p 2, then the exponential becomes very quickly negligible and one is left with the first term, which is constant. This is agreement with the above discussion about Fig. 4. In Fig. 5, we have represented δ k calculated at the end of inflation using the exact result (4.5) and the analytical approximation (4.6). In the left panel, δ k (η end ) is plotted as a function of p and for a few values of k, and in the right panel it is plotted as a function of k and for a few values of p. The coloured lines correspond to the exact result (4.5) while the black lines stand for the analytical approximation (4.6). Evidently, they match very well (and are in fact hard to distinguish).
Since decoherence at observable scales is characterised by the condition δ k * 1, Eq. (4.6) allows us to calculate the minimum interaction strength that is required for decoherence to complete before the end of inflation. One obtains that decoherence occurs when (4.8) In particular, one can see that when p < 3 + (2 − 2ν)/(1 + 1 * ) 2, decoherence and quasi scale invariance cannot be achieved simultaneously, since Eq. (3.41) and the first of Eq. (4.7) directly contradict each other. This is because, from Eqs. (3.27) and (4.6), one can see that ∆P k = 2δ k in that case.
This allows one to understand Fig. 6 where the situation is summarised. This figure represents the regions in parameter space (p, k γ /k * ) where quasi scale invariance and decoherence can or cannot be realised. The light grey region corresponds to situations where the interaction strength with the environment is too small to yield decoherence. The medium grey region is where it is too large to preserve quasi scale invariance, and the dark grey region is where both problems occur (no decoherence and scale invariance breaking). For p 2, one can see that decoherence cannot be realised without spoiling the quasi scale invariance of the power spectrum, in agreement with the above discussion. For p 2, there are intermediate values of k γ /k * for which decoherence is obtained while preserving quasi scale invariance. This region is coloured in Fig. 6, and the colour code quantifies how many e-folds since Hubble crossing are needed in order to complete decoherence. One can check that the larger the interaction strength, the fewer e-folds it takes to reach decoherence. In general, decoherence happens after Hubble crossing (light blue to red regions), but it can also occur soon after crossing out the environment correlation length and before crossing out the Hubble radius (dark blue regions).
The most striking feature of the plot is probably the thin vertical line centred at p = 5. As shown in Sec. 3.1.7, the origin of this line is the fact that for p = 5, the correction to the power spectrum caused by the Lindblad term is itself scale invariant. As a consequence, even a very large value of the coupling constant can lead to quasi N decoherence − N * Figure 6. Regions in parameter space (p, k γ /k * ) where decoherence and quasi scale invariance can or cannot be realised. The light grey region corresponds to values of p and k γ /k * where the interaction strength with the environment, parametrised by k γ /k * , see Eq. (3.4), is too small to lead to decoherence. The medium grey region is where it is too large to preserve quasi scale invariance, and the dark grey region is where both problems occur (no decoherence and scale invariance breaking). The coloured region corresponds to parameters where perturbations decohere and scale invariance is preserved. The colour code, indicated by the vertical bar, quantifies how many e-folds since Hubble crossing it takes before complete decoherence is reached. Here, decoherence is characterised by the condition δ k * > 10, and quasi scale invariance by the condition |n S −n S | < 5σ nS , wheren S 0.96 and σ nS 0.006 are the mean value and standard deviation of the Planck spectral index measurement [68] respectively. The spectral index and the decoherence parameter are computed numerically without approximations. scale invariance. Let us also recall that p 5 is precisely the case that corresponds to the model described in Appendix B where inflationary perturbations couple to heavy scalar degrees of freedom. Fig. 6 highlights again the remarkable property of this type of environment, which allows cosmological perturbations to decohere during inflation without spoiling their scale invariance.

Can we neglect the free Hamiltonian?
It is sometimes argued that decoherence can be estimated without taking the free Hamiltonian into account, since no decoherence occurs in the free theory. In principle, this is true only if decoherence is so rapid that the free evolution of the system can be neglected while it occurs. Here we re-calculate the decoherence parameter in the absence of free evolution in order to determine when this is indeed the case.
Neglecting the free Hamiltonian terms in Eq. (3.6), one has Using the Bunch-Davies vacuum as the initial state, this can be readily integrated and gives rise to Eq. (4.1) if one replaces δ k with δ no free k , defined as (4.11) This expression needs to be compared with Eq. (4.6). Like Eq. (4.6), it is made of two terms and which one dominates at late time depends on the value of p. By comparing the exponentials of each term, one finds that if p < 3 + 1/(1 + 1 * ) 4, the first term dominates. In Eq. (4.6), the dominant term depends on whether p is smaller or larger than 3+2(1−ν)/(1+ 1 * ) 2 but since neither term matches the first term of Eq. (4.11), the two expressions strongly disagree in that case. If p > 3 + 1/(1 + 1 * ) 4, however, the second term in Eq. (4.11) dominates, which matches the second term of Eq. (4.6), up to a factor [(p − 3) (1 that is of order one, and the two expressions agree in that case. Interestingly, p = 3 + 1/(1 + 1 * ) corresponds to the limiting value between cases 2 and 3 in the classification introduced in Sec. 3.1.5, where it was shown that p > 3 + 1/(1 + 1 * ) is the condition for the power spectrum to be independent of the environment correlation shape and length. We conclude that only in this case can the free Hamiltonian be neglected when computing the decoherence parameter. Further insight into this result is provided at the end of Sec. 4.1.5.

Alternative derivation of the decoherence parameter
Before moving on to study decoherence in the presence of quadratic interactions, let us show how the above result can be obtained without solving for the Lindblad equation (2.14) entirely. Indeed, in the case of quadratic interactions, a full solution to Eq. (2.14) is not available and we will need an alternative technique.
The starting point is to use Eq.
where in the second equality, the Lindblad equation (3.3) written in Fourier subspaces has been used. In this expression, using the cyclicity of the trace operator, one finds that the first term vanishes and the second one can be written as (4.13) where in the second equality, Tr(ρ s k [ρ s k ,v s k ]v s k ) has been written explicitly in terms of the elements of the density matrixρ s k . The next step is to notice that at linear order in γ, it is enough to evaluate the right-hand side of the above equation in the free theory, where the density matrix is given by Eq. (3.11). In this case, the integrals appearing in Eq. (4.13) are Gaussian and can be performed explicitly, and one obtains Tr(ρ s k [ρ s k ,v s k ]v s k ) = P vv (k). The relation (4.13) can then be readily integrated, and expanding Eq. (4.3) at leading order in γ, Tr(ρ s where in the second line we have recast the result in terms of the source function S 1 defined in Eq. (3.22). Several remarks about this expression are in order.
First, let us notice that Eq. (4.14) is valid beyond the linear expansion in γ. Indeed, if one uses the exact solution (3.7) to the Lindblad equation to calculate the integrals appearing in the right-hand side of Eq. (4.13), one obtains Tr(ρ s k [ρ s k ,v s k ]ρ s k ) = P vv (k)(1 + 4δ k ) −3/2 = P vv (k)Tr 3 (ρ s k 2 ). This allows one to write Eq. (4.13) as a linear differential equation for Tr(ρ s k 2 ), the solution of which gives rise to Eq. (4.14) when combined with Eq. (4.3). The formula (4.14) is therefore exact.
Second, at linear order in γ, the right-hand side of Eq. (4.14) can be evaluated in the free theory where P vv = |v k | 2 , where the mode function v k is given by Eq. (C.24) in terms of Bessel functions. Making use of Eqs. (2.26) and (C.25), one obtains Eq. (4.5). This elucidates the numerous cancellations that appeared when going from Eq. (4.4) to Eq. (4.5), and which result from the equivalence between Eq. (4.4) and Eq. (4.14). In fact, without expanding in γ, the right-hand side of Eq. (4.14) can be evaluated with P vv = |v k | 2 + J k , see Eq. (3.13), and this gives rise to Eq. (4.2). This explains why, since P vv is linear in γ, δ k is quadratic in γ.
Third, it is interesting to compare Eq. (4.14) with Eq. (4.10), which was obtained by neglecting the influence of the free Hamiltonian. One can see that the only difference between these two expressions is that in Eq. (4.10), the power spectrum is taken out of the integral and evaluated at the time η. This explains why only the contribution from the upper bound of the integral of Eq. (4.10) is correctly computed, up to a prefactor of order one.

Quadratic interaction
Let us now study decoherence as produced by quadratic interactions with an environment. As already explained, because of mode coupling, the Lindblad equation (2.14) does not decouple into a set of independent Lindblad equations for each Fourier mode, and cannot be solved entirely. This means that the calculation of Sec. 4.1.2 cannot be reproduced here, and that the alternative technique presented in Sec. 4.1.5 must instead be employed. Since the state is not factorisable into Fourier subspaces in the presence of quadratic interactions, i.e. Eq. (2.7) does not apply, the effective density matrix on the space s k has first to be defined by tracing over all other degrees of freedom,

Decoherence criterion
Making use of the linearity and ciclycity of the trace operator, one has (4. 16) In this expression, dρ v /dη needs to be replaced by the Lindblad equation (2.14) (witĥ A =v 2 ), which contains two terms. The first one involves the free Hamiltonian (2.1) and is, in practice, difficult to incorporate in the following calculation. However, the contribution from this term to the decoherence parameter only reflects the correlations that develop between the Fourier subspace s k under consideration and the other Fourier subspaces, over which we have traced over, see Eq. (4.15). The reason for this partial trace is not that we do not "observe" the other Fourier degrees of freedom, but is because we want to define a decoherence parameter for each Fourier subspace, by analogy with the linear case. If the state were factorisable, the Hamiltonian contribution to Eq. (4.16) would vanish, which is why it is simply discarded in what follows. The second term coming from Eq. (2.14) is the Lindblad term. Fourier expandingv and C R , it gives rise to Let us consider the case where k / ∈ {±k 2 , ±(k 1 +k 2 ), ±k 3 , ±(k 1 −k 3 )}. By removinĝ ρ s k from the trace in the second line of Eq. (4.16) (since it commutes with allv operators), one is left with a full (as opposed to partial) trace, that vanishes. This means that k must be equal, up to a sign (recall that k and k live in R 3+ while k 1 , k 2 and k 3 live in R 3 ), to one of the wavenumbers that index thev operators in Eq. (4.17). If it is equal to one such wavenumber only and the other three are different, one can show that the trace vanishes again, such that k must be equal to 2, 3 and all 4 wavenumbers that index thê v operators in Eq. (4.17). Let us discuss these three possibilities separately.
We first examine the situation where k is equal to two of the wavenumbers that index thev operators in Eq.
where we have used the decompositionv k = (v R k + iv I k )/ √ 2 . The same result is obtained with k = k 2 = k 3 −k 1 , k = −k 1 −k 2 = −k 3 or k = −k 1 −k 2 = −k 1 +k 3 , if the condition k 1 = −2k is enforced. In the same manner, if k = −k 2 = k 3 , k = −k 2 = k 1 − k 3 , k = k 1 + k 2 = k 3 or k = k 1 + k 2 = k 1 − k 3 , and if the condition k 1 = 2k is enforced, the same result as in Eq. (4.18) is obtained, except that the power spectrum is evaluated at k − k 1 instead of k + k 1 . One can check that all other configurations give a vanishing result.
Then, one can show that the case where k is equal to three of the wavenumbers that index thev operators in Eq. (4.17) gives contributions that are always proportional to the quantum mean value of a single mode function operator and therefore vanish, see the discussion below Eq. (3.16). Finally remains the situation where all wavenumbers are, up a sign, equal. For instance, let us consider the case where k 1 = −2k, k 2 = k and k 3 = −k, for which one has The only other non-vanishing configuration of this type is when k 1 = 2k, k 2 = −k and k 3 = k, which gives the same result. These contributions are however suppressed by a volume factor with respect to the ones giving Eq. (4.18), since they do not vanish only for a single configuration of the wavenumbers k 1 , k 2 and k 3 , while the configurations leading to Eq. (4.18) leave one wavenumber free. This is why, in the factorisable state, one obtains that Tr 1 − 2δ k , and this gives rise to which is analogous to Eq. (4.14). This suggests that the above formula is in fact generic, in the same way that Eq. (3.56) for the power spectrum applies for any source function.

Calculation of the decoherence parameter
Since the calculation is being performed at linear order in γ, the right-hand side of Eq. (4.22) must be evaluated in the free theory where P vv = |v k | 2 , and the mode function where σ γ has been defined in Eq. The corresponding time evolution of δ k is displayed on the left panel of Fig. 7 for different values of p and at the pivot scale k * . The main difference with the case of linear interactions shown in Fig. 4 is that, here, δ k continues to increase at late time for all values of p that are shown. This can be understood analytically by computing δ k close to the end of inflation, when the modes of astrophysical interest today are well outside the Hubble radius and Eq. (4.5) can be expanded in the limit −kη 1.
Further assuming that E H −1 * , the integrals I 3 and I 4 can be expanded according to Eqs. (D.22) and (D.26) respectively, and one obtains (4.24) In the right panel of Fig. 7, we have represented δ k calculated at the end of inflation using the exact result (4.23) (coloured lines) and the analytical approximation (4.24) (black lines), and they match very well.
In Eq. (4.24), one can see that if p > 0, the first term dominates at late time and δ k grows on large scales, in agreement with what can be seen on the left panel of Fig. 7. If, on the contrary, p < 0, the exponential in the first term becomes very quickly negligible and one is left with the second term, which is constant. Recalling that decoherence at observable scales is characterised by the condition δ k * 1, Eq. (4.24) allows us to calculate the minimum interaction strength that is required for decoherence to complete before the end of inflation, (4.25)

Combining with observational constraints
As for the case of linear interactions, these lower bounds can be combined with the upper bounds derived in Sec. 3.2.4 from the requirement that the quasi scale invariance of the power spectrum is preserved. One obtains that the range of values for the interaction strength, here parametrised by σ γ , such that decoherence occurs without spoiling scale invariance is given by (4.26) In particular, one can see that when p < 0, the correction to the power spectrum given in Eq. (3.63) and the decoherence parameter given by the first line of Eq. (4.24) are directly related, ∆P k | 3 = 2δ k , which explains why decoherence cannot occur without spoiling the quasi scale invariance of the power spectrum. The situation is summarised in Fig. 8, which represents the regions in parameter space where quasi scale invariance and decoherence can or cannot be realised. The conventions are the same as the ones used in Fig. 6 in the case of linear interactions, and all remarks made about that figure apply here too. As for linear interactions, the , is too small to lead to decoherence. The medium grey region is where it is too large to preserve quasi scale invariance, and the dark grey region is where both problems occur (no decoherence and scale invariance breaking). The coloured region corresponds to parameters where perturbations decohere and scale invariance is preserved. The colour code, indicated by the vertical bar, quantifies how many e-folds since Hubble crossing it takes before complete decoherence is reached. This map is obtained with similar conventions as in Fig. 6. striking feature of Fig. 8 is the presence of a thin vertical line centred at p = 3, for which the correction to the power spectrum caused by the Lindblad term is scale invariant, and observations do not constrain the interaction strength. Let us stress again that p 3 precisely corresponds to the model described in Appendix B where inflationary perturbations couple to heavy scalar degrees of freedom.

Generalisation to higher-order interactions
So far, we have shown that the effective inclusion of environmental degrees of freedom, through a Lindblad equation, gives rise to a modified power spectrum and to a decoherence parameter that can be calculated from a source function, see Eq. (3.56) and Eq. (4.22) respectively. In the case of linear interactions, n = 1, the source function is given by Eq. (3.22), namely For quadratic interactions, n = 2, one can rewrite Eq. (3.53) as where the power spectrum has to be evaluated in the free theory at leading order in γ. As explained in Sec. 3.2.3 (see also Appendix D.1), the convolution product between the power spectrum and the environment correlator, which appears in Eq. (5.2), is dominated by its IR contribution where p 1 is super Hubble (if the UV contribution is not simply removed by adiabatic subtraction). Since the environment correlatorC R selects out modes that have crossed out the environment correlation length E , which is much smaller than the Hubble radius, the source function is well approximated by the environment correlator times the integrated power spectrum, Before concluding this work, let us try to generalise our approach to higher-order interaction terms. For instance, let us consider a cubic interaction,Â =v 3 . In that case, a long but straightforward calculation, similar to the one presented at the beginning of Sec. 3.2.1, gives the same system as in Eq. (3.48), except that the second line of the last entry now reads (2π) 9/2 δ (k 1 + k 2 ) dp 1 dp 2 dp 3CR (p 1 + p 2 + k 1 ) T 4 (p 1 , p 2 , p 3 ) .

(5.4)
In this expression, the trispectrum has been defined according to v p 1v p 2v p 3v p 4 = T 4 (p 1 , p 2 , p 3 )δ (p 1 + p 2 + p 3 + p 4 ). One then obtains the differential equation (3.52) for the power spectrum, where the source function is given by S 3 , which is such that the quantity written in Eq. (5.4) equals S 3 (k 1 )δ(k 1 + k 2 )/2. At leading order in γ, it can be evaluated in the free theory where the state is Gaussian and one has v p 1v p 2vk3v Figure 9. Feynman diagrams representation of the source function. Straight and wiggly lines stand for propagators of the Mukhanov-Sasaki variablev, and of the environment operatorR it couples to, respectively. The first diagram corresponds to linear interactions ∝vR, and the second diagram is for quadratic interactions ∝v 2R . p 4 ) + P vv (p 1 ) P vv (p 2 ) δ(p 1 + p 4 )δ(p 2 + p 3 ) according to Wick theorem. This gives rise to +C R (k) dp 1 dp 2 P vv (p 1 ) P vv (p 2 ) .

(5.5)
As for quadratic interactions, in the limit E H −1 and keeping the IR component of these integrals only, this reduces to (5.6)

Diagrammatic calculation of the source
The formulas obtained for the source function for linear interactions in Eq. (5.1), for quadratic interactions in Eq. (5.2), and now for cubic interactions in Eq. (5.5), can be understood with the diagrammatic representation shown in Figs. 9 and 10. In these Feynman diagrams, straight and wiggly lines represent propagators of the Mukhanov-Sasaki variablev, and of the environment operatorR it couples to, respectively. The first diagram in Fig. 9 stands for linear interactions of the formvR. This is why one straight line and one wiggly line are attached to each vertex. Momentum conservation imposes that k 1 = −k 2 = k and that q = k. This is why the source function is simply proportional to g 2C R (k) ∝ γC R (k), in agreement with Eq. (5.1).
The second diagram in Fig. 9 stands for quadratic interactions of the formv 2R , which is why each vertex has two straight lines and one wiggly line. Momentum conservation imposes that k 1 = −k 2 = k and that q = k − p 1 . The loop integral then gives rise to dp 1 P vv (p 1 )C R (k − p 1 ), which indeed corresponds to Eq. (5.2).
The two left diagrams in Fig. 10 stand for cubic interactions of the formv 3R and correspond to the two ways one can have three straight lines and one wiggly line per vertex, while having a single wiggly line in the diagram. In the left top diagram, momentum conservation imposes that k 1 = −k 2 = k and that q = k − p 1 − p 2 . The loop integral then gives rise to dp 1 dp 2 P vv (p 1 )P vv (p 2 )C R (k − p 1 − p 2 ), which indeed corresponds to the first term of Eq. (5.5). In the left bottom diagram, momentum Figure 10. Feynman diagrams representation of the source function, as in Fig. 9, for cubic interactions ∝v 3R (left diagrams), and for quartic interactions ∝v 4R (right diagrams).
conservation imposes that k 1 = −k 2 = k and that q = k, so the loop integral is given by dp 1 dp 2 P vv (p 1 )P vv (p 2 )C R (k), which indeed corresponds to the second term in Eq. (5.5). The multiplicity of the top diagram is 2 since the lines labeled by p 1 and p 2 are indistinguishable, again in agreement with Eq. (5.5). This diagrammatic representation is useful since it allows one to easily guess the form of the source function if the coupling is of arbitrary order n, In this expression, j is the number of system propagators that are involved in the loop with the environment propagator, such that n − 1 − j is the number of tadpole system propagators, and a j is the corresponding multiplicity factor. Let us discuss in more detail the above expression using again the case n = 3, see Fig. 10 (left diagrams). We have seen that there are two diagrams, one with no tadpole (top-left diagram) and the other one with two tadpoles (bottom-left diagram). For the first one, one has j = 2 and, indeed, one checks that n−1−j = 0. As already mentioned, the corresponding multilplicity factor is a 2 = 2. The argument of the correlation function is given by k − 2 m=1 p m = k − p 1 − p 2 in agreement with Eq. (5.5). For the second diagram, j = 0 since we see that no system propagator is involved in a loop with the environment propagator. As already mentioned, this implies that the number of tadpoles is n−1−j = 2. The multiplicity factor is a 0 = 1. Let us also notice that the sum appearing in the argument of the correlation function is such that the upper bound (zero) is less that the lower bound (one). In that case, this simply means that this sum vanishes by notational convention. There is no diagram with one tadpole (the number of tadpoles being always even) and, hence, a 1 = 0. The total number of terms is therefore n−1 j=0 = a 0 + a 2 = 3 which is indeed the number of terms obtained by expanding a four-point correlation function by means of the Wick theorem.
This discussion can be repeated for higher-order cases. For instance, for n = 4, see Fig. 10 (right diagrams), one still has two terms, one with no tadpole and multiplicity a 3 = 6 (right top diagram) and another with two tadpoles and multiplicity a 1 = 9 (right bottom diagram), the total number of terms being 15 (using a 0 = a 2 = 0), again in agreement with Wick theorem.
In practice, we do not need to specify the value of a j for each diagram since each term in the sum over j yields the same contribution. Indeed, as for quadratic interactions, in the limit E H −1 , the IR component of the integrals appearing in Eq. (5.7) all reduce to the environment correlator evaluated at mode k times the integrated power spectrum, and one obtains Here, (2n − 3)!! = n−1 j=0 a j is, by definition, the number of possibilities to partition a set of 2(n − 1) elements into n − 1 pairs and follows from Wick's theorem, and the power spectrum is integrated on super-Hubble scales only. One can check that when n = 1, Eqs. (5.7) and (5.8) reduce to Eq. (5.1), when n = 2, they reduce to Eqs. (5.2) and (5.3) respectively, and when n = 3, they reduce to Eqs. (5.5) and (5.6) respectively.
As explained in Sec. 3.2.3, the calculation of the integrated power spectrum cannot be performed generically in the slow-roll approximation since it involves a wide range of modes over which slow-roll expansions may break down. This is why a slow-roll analysis can be carried out only on an inflationary model-by-model basis. Otherwise, neglecting slow-roll corrections, one has where we recall that η IR is an IR cutoff, that e.g. corresponds to the time at which inflation started. Making use of Eqs. (2.16) and (2.26), one obtains where k γ is a comoving scale that characterises the strength of the interaction with the environment and that is defined by (5.11) When n = 1, this is consistent with Eq. (3.4). When n = 2, the above expression is singular and the interaction strength is quantified by a dimensionless parameter σ γ instead of a comoving scale k γ , see Eq. (3.64). In Eq. (5.10), one can see that the source function varies in time as S n ∝ η 5−2n−p , if one neglects the logarithmic term and slow-roll corrections (which are in fact of the same order). The crucial remark is that if the environment consists of a heavy scalar field, we have shown that p = 7 − 2n, see Eq. (2.21), so S n ∝ η −2 , which does not depend on n anymore, and which is precisely the behaviour that produces a scale-invariant correction to the power spectrum, as we have shown when n = 1 in Sec. 3.1.7 and when n = 2 in Sec. 3.2.4. This implies that the remarkable property found in this work, namely that an environmental heavy test scalar field ψ does not spoil the observed quasi scale invariance of the power spectrum for linear and quadratic interactions, is in fact fully generic and does not depend on the order of the interaction. By linearity, it is true for any coupling of the type f (φ)g(ψ), as soon as f and g can be Taylor expanded.

Power Spectrum
As shown in Sec. 3.2.2, the form (3.56) for the power spectrum is valid for any source function. In the case of higher-order interactions, one therefore has Plugging in the result from the diagrammatic calculation of the source function detailed above, see Eq. (5.10), one obtains where the integrals I 2n−2 (ν) and I 2n−1 (ν) are defined by  19) and (D.23), respectively. When n ≥ 3 this is no longer the case, but these integrals can still be approximated in the following way. Depending on the value of α n , the integrals I 2n−2 and I 2n−1 receive their dominant contribution from the neighbourhood of their lower bound z ∼ −kη 1 (case 1 in the language of Sec. 3.1.5), from the neighbourhood of the intermediate value z ∼ 1 (case 2 in the language of Sec. 3.1.5), or from the neighbourhood of their upper bound z ∼ (H * E ) −1 1 (case 3 in the language of Sec. 3.1.5). In each case, the logarithm term does not vary much over the integration range that provides the main contribution, and so it can be approximated as a constant over that range. This corresponds to a leading-order saddle-point, or steepest-descent, approximation, and in practice, it boils down to taking Eq. (C.33) and to multiplying each term by ln n−1 (−kη IR /z), where z needs to be replaced by either −kη, 1, or (H * E ) −1 , depending on where the contribution comes from. This gives rise to Obviously, when n = 1 this reduces to Eq. (C.33). When n = 2 however, Eqs. (D. 22) and (D. 26) are recovered only at leading order in the logarithms. In principle, one could carry out the saddle-point approximation at higher orders, by Taylor expanding the logarithmic terms, but in practice, Eq. (5.15) already provides reliable estimates. Together with Eq. (5.13), it leads to the following relative correction to the power spectrum The three cases are similar to the ones discussed in the situation of linear or quadratic interactions. If p < 6 − 2n, the power spectrum freezes out on large scales and the amplitude of the correction to the standard result depends on the environment correlation length E . If 6 − 2n < p < 10 − 2n, the power spectrum still freezes out on large scales but its amplitude no longer depends on E . Since the correction scales as k p+2n−7 , the specific case p 7 − 2n preserves the quasi scale invariance of the power spectrum, and this precisely corresponds to the situation where the environment consists of a heavy test scalar field, as already pointed out. If p > 10 − 2n, the power spectrum continues to increase on large scales. Except in the case p 7 − 2n, the observed quasi scale invariance of primordial cosmological fluctuations thus places an upper bound on the interaction strength with the environment, here parametrised by k γ /k * , that will be given below.

Decoherence
For now let us determine whether decoherence proceeds before the end of inflation or not. In the case of linear and quadratic couplings, we have shown that the decoherence parameter is given in terms of the source function by Eq. (4.22). Strictly speaking, we have not shown that this formula generalises to higher orders. One can however easily convince oneself that a calculation similar to the one presented in Sec. 4.2.1 yields the same structure for the decoherence parameter at higher values of n, namely The only unknown is the numerical prefactor in Eq. (5.17), but since it was found to be the same in the cases n = 1 and n = 2, here we assume that it does not depend on n (even if it did, it would not change the observational constraints derived below, since they are based on order-of-magnitude considerations only) and is therefore still given by 1/2. This leads to δ k (N ) = π 8 sin 2 (πν) (5.19) If p < 4 − 2n, the first term in the braces dominates at late time, which implies that δ k reaches a constant value that depends explicitly on the environment correlation length E . If p > 4 − 2n, the second term dominates and δ k continues to increases on large scales, leading to more efficient decoherence.

Observational constraints
Requiring that decoherence occurs by the end of inflation at observable scales [δ k * (N end ) 1] but that the power spectrum remains unaltered at observable scales [∆P k * (N end ) 1] finally leads to the constraint (5.20) In particular, one can see that when p < 4 − 2n, decoherence can never occur without spoiling the quasi scale invariance of the power spectrum, since from Eqs. (5.16) and (5.18), one has ∆P k 2δ k in that case. When n = 1 and n = 2, Eq. (5.20) reduces to Eq. (4.8) and Eq. (4.26), respectively. This implies that the structure of Figs. 6 and 8 generalises to higher-order couplings.

Case of a massive scalar field as the environment
The only exception evading these constraints is the case p 7 − 2n, for which the correction to the power spectrum is itself quasi scale invariant and only the lower bound on k γ /k * applies. More precisely, the quasi scale-invariant correction to the power spectrum coming from the environment may improve or deteriorate how a given model fits the data, but one has to study each model separately, as was done in Sec. 3.1.7 for the case of linear interactions, and there is no model-independent conclusion to be drawn. Let us also mention that this requires to incorporate slow-roll corrections into the calculation, which we did for linear interactions, but which otherwise involves the calculation of the source function beyond the de-Sitter limit, that a priori depends on the field dynamics over the entire inflationary period.
As already pointed out, it is remarkable that p 7 − 2n precisely corresponds to the model proposed in Appendix B where the environment is made of the degrees of freedom contained in a heavy test scalar field. In this case, combining Eqs. (2.21) and (2.19) with Eq. (5.11), one obtains In terms of the physical parameters of the model, we therefore find that, using Eq. (5.20), decoherence occurs if the coupling constant satisfies λ > λ decoherence , where (5.22) For instance, if we take ∆N * = 50, N T = 10 4 , M = 100 H * and µ = H * , we find that λ decoherence ∼ 10 −25 , 10 −20 and 10 −15 for n = 1 and m = 1, 2 and 3 respectively, λ decoherence ∼ 10 −27 , 10 −22 and 10 −17 for n = 2 and m = 1, 2 and 3 respectively, and λ decoherence ∼ 10 −28 , 10 −24 and 10 −19 for n = 3 and m = 1, 2 and 3 respectively. There are however cases where the critical value for the coupling constant λ above which decoherence occurs is not small. For instance, with n = 1 and m = 3, µ and N T are irrelevant and taking ∆N * = 50, one finds that λ decoherence > 1 as soon as M/H * > 7030. This is because if the environmental scalar field is too heavy, its condensate is too suppressed to yield efficient decoherence of the system. One concludes that parametrically small values of the interaction strength λ are in general enough to lead to successful decoherence of primordial cosmological perturbations, but only if the heavy test scalar field they couple to has a mass no more than a few orders of magnitude larger than the Hubble scale. Furthermore, if λ > λ P ζ , where the power spectrum is modified, but in a quasi scale-invariant (though model-dependent) way. With ∆N * 50, one has λ P ζ 10 32 λ decoherence , so there always is a wide range of value for λ such that decoherence occurs while leaving the power spectrum unchanged.

Conclusion
Let us now recap our main results and discuss possible extensions. In the early Universe, cosmological density perturbations are amplified from vacuum quantum fluctuations, and subsequently seed the formation of all structures in our Universe. This mechanism occurs in the presence of all degrees of freedom present in the standard model and beyond, to which cosmological fluctuations couple (at least gravitationally). They should therefore be described as an open quantum system (as opposed to an isolated one, as usually done), the evolution of which can be modelled with a Lindblad equation, under some conditions that we have clarified. This modified evolution leads to corrections to observable predictions such as the power spectrum of curvature fluctuations. Measurements of the CMB temperature and polarisation anisotropies [14] therefore constrain the properties of possible environments and place upper bounds on the interaction strengths, that we have derived. On the other hand, the Lindblad evolution also leads to decoherence of the system in the eigenbasis selected by the form of the interaction with the environment. Since quantum decoherence is thought to play a role in the quantum-to-classical transition of cosmological fluctuations, one may also require that decoherence has completed by the end of inflation, which places lower bounds on the interaction strength. We have then identified the viable scenarios where decoherence occurs without spoiling the quasi scale invariance of the power spectrum, see Figs. 6 and 8 .
In practice, we have considered local interactions of the formĤ int ∝ dxv n (x, t) ⊗ R(x, t), wherev is the Mukhanov-Sasaki variable that describes scalar cosmological fluctuations andR is the operator in the environment sector to whichv couples. In the case of linear interactions, n = 1, we have shown that the Lindblad equation can be solved completely, see Eqs. (3.7)-(3.10). In this case, the state remains Gaussian, and the power spectrum and the decoherence parameter can be calculated from the (modified) density matrix directly. For higher-order coupling, n ≥ 2, this is no longer the case, but we have shown how the power spectrum and the decoherence parameter can still be calculated exactly, see Eq. (5.12) and Eq. (5.17) respectively, in terms of a source functions that involves the correlator ofR in the environment sector and correlators of v in the (free limit of the) system, see Eq. (5.7). Since Eq. (5.12) and Eq. (5.17) linearly depend on the source function, this means that we have entirely solved the problem, for any interaction of the formĤ int ∝ dxf (v)(x, t) ⊗R(x, t) as long as f can be Taylor expanded.
As an illustration, we have discussed the situation where the environment is made of a heavy test scalar field ψ, see Appendix B, andR ∝ψ m . In that case, the time dependence of the effective interaction strength scales as a 7−2n , where a is the scale factor. We have shown that this precisely corresponds to the very peculiar configuration where the correction to the power spectrum is quasi scale invariant. The observational constraints on the parameters describing the environment (here the mass of the heavy test scalar field in Hubble units) and the interaction strength are therefore less straightforward, and, as was shown around Fig. 2 in the case of linear interactions n = 1, are in fact model dependent. Indeed, there are models for which the correction to the power spectrum coming from the environment is too small to be resolved by current CMB measurements and no constraint can be derived, there are models for which the correction improves the fit to the data and lower bounds on the interaction strength can be obtained, and there are models for which the correction deteriorates the fit to the data and upper bounds on the interaction strength can be derived. For higher-order interactions, n ≥ 2, such an analysis requires to calculate the source function beyond the de-Sitter limit, which again has to be done on a model-by-model basis, and which depends on the inflaton field dynamics over the entire inflationary period. In principle, this may allow one to probe the inflationary potential beyond the last ∼ 50 e-folds of inflation and might extend the range of scales one can access in the early Universe beyond the observable horizon [74,75]. We leave such an analysis for future work, and we now mention a few other possible prospects.
Let us first emphasise that our results assume that the power spectrum remains unaffected during preheating and reheating, as it is the case on large scales in the standard approach. In the presence of interactions with extra degrees of freedom, this should however be verified. Second, in the case of non-linear interactions, n ≥ 2, we have shown that the quantum state of cosmological fluctuations, modified by its interaction with the environment, is no longer Gaussian. Since tight upper bounds on the amount of primordial non-Gaussianities have been placed from recent CMB measurements [76], this constitutes another channel through which the environment properties and its interaction strength with cosmological fluctuations can be constrained [77]. Third, in the case where the effective interaction strength with the environment evolves in time as a p with p > 7 − 2n, the power spectrum is modified on small scales, with a blue tilt [comprised between n S = 1 and n S = 4 depending on the value of p, see Eq. (5.16)]. This implies that the amplitude of the power spectrum may reach sizeable values for scales that are smaller than the ones probed in the CMB but that are still of astrophysical interest. Observational bounds on the amount of primordial black holes therefore constitute yet another channel to constraint the environment and its interaction strength with the system. Conversely, this also suggests that quantum decoherence might be a promising candidate to produce such primordial black holes, if their role in the dark matter abundance of our Universe or in providing progenitors to the LIGO black-hole merging events is confirmed (see e.g. Ref. [78]).
Let us also mention that decoherence, per se, does not solve the quantum measurement problem [53,54], 5 but that there are alternatives to the standard formulation of quantum mechanics that do so, e.g. dynamical collapse models of the wavefunction such as the CSL proposal [79][80][81][82]. In this model, a non-linear and stochastic correction is added to the Schrödinger linear equation, that collapses to wavefunction towards one of the eigenstates of the operator appearing in the extra term. Interestingly, the averaged density matrix (where "averaged" here refers to the stochasticity of the theory) precisely satisfies a Lindblad equation. Therefore, the present work also allows one to compute observable corrections in the CSL theory, in a way that is complementary 6 to Ref. [26].
Finally, other effective approaches [83][84][85] have been proposed to incorporate heavy scalar fields and their effects on cosmological fluctuations in the early Universe, and it would interesting to compare them with the one used in this work. Let us however stress that our approach is not limited to the case where the environment consists of heavy test scalar fields but is fully generic, and does not assume anything about the environment apart from the validity conditions for the Lindblad equation. It therefore allows one to carry out a model-independent analysis of environmental influence and quantum decoherence in the early Universe.

A Deriving the Lindblad equation
In this section, we provide a detailed and generic derivation of the Lindblad equation. Following the usual text book discussions, see e.g. Refs. [46,[86][87][88][89], we pay special attention to the physical assumptions the Lindblad formalism relies on, and show how they concretely enter into the calculation.
Let a system "S" be in interaction with some environment "E". The Hilbert space E of the full system can be written as the tensorial product of the Hilbert space of the system, E S , with the Hilbert space of the environment, E E , namely E = E S ⊗ E E . Then, the corresponding Hamiltonian reads 7 Here, H S is the Hamiltonian of the system and acts in E S , while I E is the identity operator acting in E E . In the same manner, H E is the Hamiltonian of the environment and acts in E E , while I S is the identity operator acting in E S . They represent the free Hamiltonian H 0 acting in the full space E, while H int is an interaction term. It carries a (supposedly small) dimensionless coupling parameter g characterising the strength of the interactions between the system and the environment. The density matrix ρ of the full system (acting in the Hilbert space E) obeys the unitary Liouville-von Neumann equation In what follows, it will be convenient to factor out the time dependence of ρ due to the free Hamiltonian H 0 by going to the interaction picture. This is why we introducẽ As a consequence, the evolution equation ofρ(t) reads which can be formally integrated as This expression gives rise to a Born expansion in g of the solution of Eq. (A.5). Indeed, one can iteratively expand the integrand of Eq. (A.6) to obtaiñ This expression is an explicit solution of Eq. (A.5) at order g 2 . In the second term, the density matrix is evaluated at the time t, but evaluating it at any other time comprised between t and t + ∆t only gives a correction of order g 3 . This is why, for future convenience, we chose to evaluate it at the time t instead and to work with the expressioñ From now on, for display convenience, we will drop the O(g 3 ) and remember that the calculation is performed at order g 2 .
Let us now restrict the analysis to the reduced density matrix of the system,ρ S . It is obtained from the full density matrix by tracing out the environment degrees of freedom, i.
This relation defines the quantityρ correl , which describes the correlations between the system and the environment at time t. It satisfies 8 Tr E (ρ correl ) = 0 and Tr S (ρ correl ) = 0.
It is clear that, if we start from a situation where the density operator is factorised and ρ correl = 0, correlations can only appear if the interaction term is switched on. Hence, the termρ correl must carry some g-charge, which is explicitly displayed in Eq. (A.11) as g p , p being an unknown natural integer. Plugging Eq. (A.11) into Eq. (A.10), one then obtains four terms, (A.12) In order to determine which of these terms dominate, we now need to specify the interaction Hamiltonian H int . Let us first assume that it can be written as where A acts in E S and R acts in E E . More generic interaction Hamiltonians will be considered below, see Eq. (A.37). The evolution operator U can be factorised as U S ⊗ U E because it describes the time evolution in the case where H int = 0, that is to say when the system and the environment evolve independently. As a consequence, has (A. 18) In order to proceed further, we need to formulate a few approximations. Since the system is supposed to be "small" compared to the environment, let us first assume that its influence on the evolution of the environment is negligible. Under this condition, Here, the trace is taken in E E sinceρ E and R act in E E . Notice that, using the cyclic property of the trace and the fact that the density operatorρ E commutes with U E , this also means that 21) This implies that the right-hand side of Eq. (A.18), hence the first term of Eq. (A.12), vanishes.
This now allows us to determine the value of the positive integer p. Indeed, at leading order in g, the left-hand side of Eq. (A.12),ρ S (t + ∆t) −ρ S (t), is directly proportional to g pρ correl , since, in the absence of the interaction term,ρ S does not evolve in the interaction picture. The left-hand side of Eq. (A.12) is therefore of order p, while the right-hand side contains terms of order p + 1, 2 and p + 2. The only possibility that allows us to identify the dominant terms in both sides is that p = 2. As a consequence, the terms of order p + 1 and p + 2 in the right hand side of Eq. (A.12) are sub-dominant and This expression is valid at leading order in g, which is why we need to make a fourth assumption, namely that the interactions modify the dynamics of the system in the perturbative regime only.
Let us now evaluate the remaining term of Eq. (A.12), that is to say, the right hand side of Eq. (A.22). Plugging in Eq. (A.13) and expanding the double commutator, one has (A.23) where we have introduced the two-point correlation function Because the environment is in a stationary state, one can show that C R (t, t ) is in fact a function of τ ≡ t − t only. Indeed, one has where in the last step of the calculation, we have used the fact that the reduced density operator of the environment,ρ E , and the Hamiltonian for the environment, H E , commute. A more explicit form of the correlator C R (τ ) can be obtained by making use of Eq. (A.19). Indeed, one has In particular, one has C R (−τ ) = C * R (τ ). More specifically, one can see that C R (τ ) is a sum of exponentials oscillating at the Bohr frequencies of the environment. In the limit where the environment is large and contains an almost continuous set of energy levels, destructive interference occurs that quickly drives C R (τ ) to zero with a characteristic time t c , C R (τ ) C R (0)e −|τ |/tc .
Let us now further simplify the right-hand side of Eq. (A.22). Using Eq. (A.23), one has (A.27) The integration domain appearing in this expression is displayed in Fig. 11 as the hatched blue surface. When parametrised in terms of the variables t and τ = t − t , it is given by the integrand vanishes when |τ | t c , hence its support is given by the pale green stripe in Fig. 11. Let us now consider the extended integration domain where the upper bound on τ has been extended to infinity and the lower bound on t to t. It is displayed as the hatched red surface in Fig. 11. One can see that compared to the initial integration domain, two regions are added. The first one lies outside the integrand support and therefore negligibly contributes to the overall result. The second one is the small triangle that lies inside the integrand support. In the limit where t c ∆t, it corresponds to a small area compared to the integrand support comprised in the initial integration domain. As a consequence, its inclusion in the integration domain negligibly changes the result if one follows the evolution of the reduced density matrix for the system on time scales ∆t much larger than the typical correlation time of the environment, t c ∆t .
Under this fifth and last assumption, one then has (A.31) The time derivative ofρ S can then be obtained by dividing the left-hand side and the right-hand side by ∆t. If ∆t is much smaller than the typical time scale by which A varies, , the difference only giving rise to ∆t suppressed quantities. Note that because of that, the condition (A.30) actually means that the interaction operator A should vary on time scales much larger than the autocorrelation time of R in the environment. Moreover, the variation ofρ S between times t and t + ∆t in the above integrals is of order g 2 . Since the whole expression is already proportional to g 2 , in the integrals, one can simply writeρ S (t). As a consequence, the integral over t becomes trivial and ∆t can be factorised out of the right-hand side, so that This is a Markovian master equation, since the right hand side only depends onρ S at time t. Defining where the last equality is valid ifÃ is hermitian, it can be written in the more compact form The operators L 1 and L 2 can be further simplified under the condition t c ∆t. Indeed, because of the profile C R (τ ) = C R (0)e −|τ |/tc established above, the integrals L 1 and L 2 are dominated by their contribution on the interval τ ∈ [0, a few t c ]. In the limit where t c ∆t, the functionÃ(t − τ ) (which, as already said, varies with a typical time scale much larger than ∆t) does not evolve much in this interval and one then has (A. 34) and the same expression for L 2 . As a consequence, Eq. (A.34) reads Going back to the standard picture, one finally obtains This is the standard form of the Lindblad equation. It can be generalised to cases where the interaction term is not simply given by Eq. (A.13), namely by the product of an operator acting in E S and an operator acting in E E . Indeed, if one considers a generic interaction Hamiltonian one can define correlators in the environment as and all the previous steps can be repeated, leading to Here we have assumed that the assumptions discussed before still hold and that C R,ij = C R,ji . The characteristic times of the correlation functions C R,ij are denoted t c,ij and must all be much smaller than the typical time scale over which the system varies. If i and j are replaced by continuous indices x and y, the interaction Hamiltonian (A.37) is of the form (2.12), and Eq. (A.39) gives rise to the Lindblad equation (2.14). In this expression C R (x, y) denotes C R,ij (0), any potential dependence of t c on x and y is absorbed in the function C R (x, y), and γ = 2g 2 t c .
(A. 40) In conclusion, let us summarise the conditions under which the Lindblad equation has been obtained. This equation describes the evolution of the reduced density operator of the system, perturbatively coupled to the environment through an operator A ⊗ R where A acts on the system and R on the environment, on time scales much larger than the correlation time of R. It is valid only when, in the interaction picture,Ã varies with time scales much larger than this correlation time. The environment is assumed to be in a stationary state, which technically means that R (t) R(t ) must depend on t − t only, and sufficiently large not to be affected by its interaction with the system. Finally, the above equation is valid if the interaction Hamiltonian has vanishing mean value in the environment, R = 0, but this can always be achieved by properly redefining the free and interaction Hamiltonians.

B Concrete example: a massive scalar field as the environment
In this section, we present a concrete example where the Lindblad formalism can be applied. We consider the case where two scalar fields φ and ψ are present during inflation, φ being very light and ψ very heavy. The Fourier modes of φ play the role of the system, coupled to an environment made of the degrees of freedom contained in ψ. This allows us to show how the assumptions used in Appendix A in order to derive Eq. (2.14) can be realised in practice, how the parameters γ and C R (x, y) that appear in that equation can be related to microphysical quantities, and how the Lindblad operator A can be concretely identified.
Let us consider the following action is the potential of the field φ that we leave unspecified for the moment, M is the mass of the field ψ that we assume to be much larger than the Hubble scale, M H, λ is a (supposedly small) dimensionless coupling constant and µ is a mass scale parameter that appears in the power-law coupling between φ and ψ. For now we assume both φ and ψ to be test fields (i.e. they do not contribute much to the energy budget of the Universe) in an inflating background, but below we comment on how the result can be generalised to the case where φ is the inflaton field and the system corresponds to the observed cosmological perturbations.
A first important remark is that the action should be written in such a way that the quantum mean value of the interacting term vanishes in the stationary configuration of the environment, as required by Eq. (A.20). Following the procedure described in Appendix A, this can be easily done by adding and subtracting λµ 4−n−m ψ m st φ n , where ψ m st denotes the stationary quantum mean value of ψ m . This simply modifies the effective potential of the field φ according to On the other hand, it is now clear that the mean value of the interacting term taken in the environment sector, ψ m − ψ m st , vanishes. The action (B.2) can be decomposed as S = S φ + S ψ + S φψ , and S φψ gives rise to the interaction Hamiltonian Assuming that V eff (φ) = m 2 φ 2 /2, in Fourier space, one has where we have defined v(η, x) ≡ a(η)φ(η, x). If one lets m = 0 in Eq. (B.5), then one recovers the action for curvature perturbations if one ignores metric perturbations (and it is even exactly the same for the particular case of power-law inflation). This suggests that, by identifying the variable v(η, x) with the Mukhanov-Sasaki variable (which is the reason why we have used the same notation) in the uniform-curvature gauge, we can extend the present analysis to the case where φ is the inflaton field and the system corresponds to the curvature perturbations.
In terms of the Mukhanov-Sasaki variable, the interaction Hamiltonian is given by This is of the form (2.12) provided A = v n , R = ψ m − ψ m st , and the effective coupling constant reads g = λµ 4−n−m a 4−n which is, therefore, a time-dependent quantity. In Appendix A, we have showed that γ = 2g 2 τ c , see Eq. (A.40). In this expression, τ c is the correlation time of R. Given that, in the main text, the Lindblad equation is written in terms of conformal time, τ c in this context must be interpreted as a conformal correlation time and expressed as t c /a, t c being the cosmic correlation time. This implies that the ansatz (2.16) is satisfied, γ = γ * (a/a * ) p , if we make the identification and p = 7 − 2n. For a linear interaction, this leads to p = 5 and for a quadratic interaction, one has p = 3.
The explicit form of R also allows us to calculate the correlation function C R . Massive fields in de-Sitter space-times have long been studied [59][60][61]. In the smallseparation limit, 2 ≡ [(t 1 − t 2 ) 2 − a 2 (x 1 − x 2 ) 2 ]/4 min(1/H 2 , 1/M 2 ), point-splitting renormalisation yields the two-point correlation function given by Eq. (3.14) in Ref. [60]. In the regime where M H, expanding this formula in powers of H/M leads to at leading order in H/M , where Σ = ±1 depends on whether the separation between the two points (t 1 , x 1 ) and (t 2 , x 2 ) is timelike or spacelike. Assuming that ψ has Gaussian statistics (which is correct at leading order in perturbation theory), Wick theorem leads to In this expression, m/2 denotes the integer part of m/2, i.e. it is m/2 is m is even and (m − 1)/2 if m is odd, and a p,m are combinatory coefficients that can be calculated as follows. In the product of pairs written in the sum of Eq. (B.9), 2p is the number of ψ(t 1 , x 1 ) occurrences that go into auto-correlators and m − 2p is the number of ψ(t 1 , x 1 ) occurrences that go into cross-correlators. The number of ways to split the m occurrences of ψ(t 1 , x 1 ) between these two types of correlators is given by m 2p = m!/[(2p)!(m−2p)!], i.e. by the number of ways one can draw 2p elements out of m elements (or equivalently m − 2p elements out of m elements). The same applies for dispatching the occurrences of ψ(t 2 , x 2 ). Then, once the 2p occurrences of ψ(t 1 , x 1 ) that appear in auto-correlators are chosen, one has to arrange them into pairs, and the number of such arrangements is given by (2p − 1)!! = (2p)!/(2 p p!). The same applies for arranging into pairs the 2p occurrences of ψ(t 2 , x 2 ) that go into auto-correlators. Finally, the number of cross-correlators one can build from the m − 2p occurrences of ψ(t 1 , x 1 ) and ψ(t 2 , x 2 ) is given by (m − 2p)!.
Putting everything together, one obtains Let us also note that in Eq. (B.18),C R is proportional to H 6m . In a Universe that is inflating in the slow-roll regime, H is not strictly constant but scales as a − 1 , where 1 is the first slow-roll parameter. This slow time variation can be absorbed in p by replacing p → p − 6m 1 * , leading to p = 7 − 2n − 6m 1 * , (B.20) and replacing H with H * in Eq. (B.18). Let us note that this assumes that Eq. (B.8) is valid even when H varies with time. In that case however, there is no known analytical solution for the two-point correlation function of the field ψ. Nevertheless, since M H, the relaxation time of this correlation function, 1/M , is much smaller than the typical time scale over which H varies, 1/(H 1 ), and one can assume ψ to adiabatically follow Eq. (B.8).
In Appendix A, it was made clear that the Lindblad equation relies on the validity of a number of assumptions. Let us now derive the conditions under which those assumptions are verified in the simple model considered here.
First, one needs to check that ψ is a test field (strictly speaking this is not a required condition for the derivation of the Lindblad equation but is necessary for the consistency of the model presented here). In order for ψ to play a negligible role in the energy budget of the Universe, M 2 ψ 2 st has to be small compared to the total energy density of the Since M H, and since H/M Pl 10 −4 due to the current observational bound [68] on the tensor-to-scalar ratio, this number is necessarily very small and the condition is satisfied.
Second, if φ is taken to be a test field, a similar condition must be satisfied, V eff (φ) 3M 2 Pl H 2 , where V eff (φ) is given in Eq. (B.3). Assuming that the initial potential V (φ) where a prime here denotes derivation with respect to τ . In order to remove the contribution from the partial derivatives along d and r in this equation, let us now choose the characteristic line such that it satisfies η (τ ) = 1 , d (τ ) = −r(τ ) , r (τ ) = ω 2 (k) d(τ ), (C. 6) i.e. in such a way that Eq. (C.5) can be written as an ordinary differential equation of the variable τ only, namely Let us now specify the initial conditions. Requiring the initial state to be in the Bunch-Davies vacuum, 11 one has 12 (C.11) One still needs to solve for the functions η (τ ) and d(τ ) since these functions appear in the above expression. This can be done through the integration of Eqs. (C.6), which leads to η(τ ) = τ, (C.12) d (τ ) = u k (τ ) d in + w k (τ ) r in , (C.13) r (τ ) = −u k (τ ) d in − w k (τ ) r in . (C.14) Here we have chosen the first integration constant such that η in = τ in , and since d(τ ) and r(τ ) satisfy two coupled linear first-order differential equations, we have simply written that d is a linear combination of d in and r in , and calculated r = −d accordingly. The 11 Since the Lindblad correction to the standard dynamics ofρ s k vanishes when a/k < E, see Eq. (2.26), it does not spoil the ability to set the Bunch-Davies vacuum in the sub-Hubble limit, as long as one also is in the sub-E limit. 12 In the sub-Hubble limit, ω 2 (k) k 2 and the ground state is given by Eq. (A.30). In Appendix A, this condition is shown to be necessary in order for the Lindblad equation to be valid. As explained in Appendix B, in the present case the time scale over which the system evolves is the Hubble time, ∆t ∼ H −1 . Furthermore, if the environment correlation time and length are directly related, E ∼ t c , as is the case in the example discussed in Appendix B, this condition boils down to H E 1. Because of Eq. (C.30), it implies that −kη E 1. The second limit consists in evaluating the above expressions when the physical wavelength of the mode under consideration a/k is well outside the Hubble radius H −1 , which is the case at the end of inflation for all modes of astrophysical interest today. This amounts to taking −kη 1. Under these two conditions, the hypergeometric functions of Eqs. (C.31) and (C.32) can be expanded, in the large third argument limit for the first one and in the small third argument limit for the second one, and one obtains I 1 (ν) (−kη E ) α 1 α 1 π + 1 2 √ π Γ (−α 1 /2) Γ (1/2 + α 1 /2 + ν) Γ (1/2 − α 1 /2) Γ (1/2 − α 1 /2 + ν) − (−kη) 1+α 1 +2ν 4 ν (1 + α 1 + 2ν) Γ 2 (1 + ν) , (C.33) For the power spectrum, one has to evaluate the integral J k , see Eq. (3.13). Expanding Eq. (C.27) in the same limits as above, namely −kη E 1 and −kη 1, one has J k (η) π 2 k 2 γ 8k 3 sin 2 (πν) (−kη) 2−α 1 a a * p−3 One can then see which term dominates in the above expressions of I 1 and I 2 and, hence, in the power spectrum, depending on the value of α 1 (hence p). The critical values for α 1 are −2ν − 1, −2ν, −1, 0, 2ν − 1, 2ν so we have seven cases to distinguish a priori, but a more careful study reveals that, in fact, they can be casted into three cases only. If α 1 < −2ν − 1, the dominant terms come from the third ones of I 1 (ν), I 1 (−ν) and I 2 (ν). If −2ν − 1 < α 1 < 0, the dominant term comes from the second one of I 1 (ν). If α 1 > 0, the dominant term comes from the first one of I 1 (ν). These considerations lead to the expressions of ∆P i in the main text, see Eqs. (3.25), (3.26) and (3.27). For decoherence, as explained in the main text, one has to evaluate δ k = |v k | 2 J k + |v k | 2 I k − |v k | 2 K k , see Eq. (4.4). Using Eqs. (C.24), (C.26), (C.27) and (C.28), one can show that the connection and recurrence relations for Bessel functions [90,91] lead to several cancellations and this expression simplifies to δ k = I 1 (ν)+I 1 (−ν)−2 cos (πν) I 2 (ν), see Eq. (4.5). If α 1 < 0, the dominant contribution comes from the third term of I 1 (−ν), if 0 < α 1 < 2ν − 1, the dominant contributions come from the third term of I 1 (−ν) and the first terms of I 1 (ν), I 1 (−ν) and I 2 (ν), and if α 1 > 2ν −1, the dominant contributions come from the first terms of I 1 (ν), I 1 (−ν) and I 2 (ν). This leads to the expression of δ k in the main text, see Eq. (4.6).

D Power spectrum for quadratic interaction
In this appendix, we provide additional details for the calculation of the power spectrum in presence of quadratic interactions.

D.1 Calculation of the source
In Sec. 3.2.1, it is shown that the source term in the differential equation satisfied by the power spectrum, Eq. (3.51), involves the convolution product between the power spectrum itself and the Fourier transform of the environment correlator, where we have simply performed the change of integration variable k = p − k. Let us note that S 2 is nothing but 8γ(2π) −3/2 I. If θ denotes the angle between the vectors p and k, in spherical coordinates where θ is the polar angle, one has I = 2π ∞ 0 dp p 2 π 0 dθ sin θC R k 2 + p 2 − 2kp cos θ P vv (p) , where the integral over the azimuth angle has been performed. After changing the integration variable z = k 2 + p 2 − 2kp cos θ, one obtains I = π k ∞ 0 dp p P vv (p) This expression is useful since it allows one to perform one-dimensional integrals only.
In what follows, we perform an explicit calculation of this integral using the ansatz (2.17) for the environment correlator. The result can be generalised to any environment correlation function by using the argument presented below Eq. where Ei(x) ≡ − ∞ −x e −t dt/t is the exponential integral function. For the reason recalled above, we assume the correlation length of the environment to be much smaller than the Hubble radius, which amounts to − E /(aη) 1. This also implies that − E /(aη IR ) 1 since, by definition, η > η IR during inflation. As a consequence, the limit y → 0 is the relevant one, where (D.14) As mentioned above, the UV part is usually removed from the final result through adiabatic subtraction. In any case, one can see that, from Eq. (D.8), I UV is roughly constant when the mode under consideration is sub Hubble, from Eq. (D.9), I UV decreases when the mode is super Hubble, I UV ∝ 1/a, and from Eq. (D.14), I IR increases as a function of time when the mode has crossed out the environment correlation length, I IR ∝ a 2 . This implies that, at late time, I UV I IR and one can take I I IR anyway. This leads to the source function (D. 16) Then, using the explicit form of the Bunch-Davies normalised mode function given by Eq. (C. 24), one arrives at the following expression P vv (k) =|v k | 2 + π 3kC R sin 2 (πν) γ * with ν ≡ 3/2 (since, as explained in Sec. 3.2.3, we neglect all slow-roll corrections) and