Quantum corrections to the primordial tensor spectrum: Open EFTs&Markovian decoupling of UV modes

Perturbative quantum corrections to primordial power spectra are important for testing the robustness and the regime of validity of inflation as an effective field theory. Although this has been done extensively for the density power spectrum (and, to some extent, for the tensor spectrum) using loop corrections, we do so in an open quantum system approach to the problem. Specifically, we calculate the first order corrections to the primordial gravitational wave spectrum due to (cubic) tensor interactions alone. We show that our results match expectations from standard loop corrections only in the strict Markovian limit, and therefore, establish a systematic way to relax this approximation in the future, as is generally necessary for gravitational systems.


Introduction
As with any successful theory of nature, it is important to question the realm of validity of inflation as an effective field theory (EFT) on curved spacetimes [1,2]. Inflation is generally considered as the relevant theory to describe long-wavelength degrees of freedom (dofs) as originating from quantum fluctuations in the very early universe. However, a natural question often posed is how far back can we trust inflation as an EFT of quantum fields on a quasi-de Sitter (dS) background. This typically goes under the name of the 'trans-Planckian problem' of inflationary cosmology [3,4]. Simply put, classical inhomogeneities observed at late times, if traced back long enough, could have arisen from trans-Planckian modes, a regime in which the predictions of inflation cannot be trusted anymore [5]. Instead of going into the debate of understanding the pros and cons of whether this is indeed a physical problem (see, e.g. [6-9] for recent perspectives), a safe statement is that going beyond linear interactions between the cosmological fluctuations should be a pointer towards understanding the depth of this problem. In other words, if this was to be a physical problem, then quantum corrections to the linear two-point function (i.e. the primordial power spectrum) would lead to large O(1) corrections to the free theory calculation since the interactions would become strongly coupled at the Planck scale (M Pl ) [10].
The reader might pause at this point and say that we already know that loop corrections to the inflationary power spectrum are small [11,12], and therefore, this unambiguously shows that there is no trans-Planckian problem of inflationary cosmology. However, the main point of this article will be to show that the standard arguments regarding the nonlinear quantum effects to the two-point function is correct to the leading order, and yet there are many subtleties which have to do with the limits of validity of the approximations which go into standard loop corrections. To identify these caveats, let us begin by first spelling out that we are interested in the EFT of the super-Hubble modes during inflation which later re-enter the Hubble horizon to become late-time inhomogeneities [13][14][15][16][17]. In other words, we are choosing to divide our full Hilbert space into system and environment dofs (as is always done in physics) -either in the form of 'heavy' vs. 'light' modes, 'slow' vs. 'fast' variables, 'high' vs. 'low' energies and so on. And even then, what is remarkable about our choice is that this immediately tells us that the usual formulation of organizing an EFT in terms of energies -the standard Wilsonian effective action -will not be suitable for us.
Let us unpack the last assertion in a bit more detail. The physical wavelength of primordial perturbations are stretched due to gravitational red-shifting over cosmological expansion. In the case of inflation, if we demarcate super-Hubble modes as the 'system' dofs 1 , then the Hilbert space corresponding to these keeps growing with time. In these types of scenarios, the system can gain or lose energy to the bath (or environment), and thus, the system undergoes non-unitary evolution [18][19][20]. In other words, there is no 'energy conservation' law implying that the heavy modes which are integrated out can reappear later on. This forbids us to write down an effective Wilsonian action for the low-energy dofs in terms of operators defined only in the low-energy Hilbert space. All of this points towards adopting a more general EFT formulation for cosmological perturbations.
In fact, it is often the case that in theories with horizons, it is better to use an open quantum systems approach [21,22]. In open EFTs, one only has access to the system modes while tracing out over the environment modes. The crucial difference with standard Wilsonian integrating out is that the unobserved sector is not demarcated on the basis of energy of the modes. Rather, we will consider the physical wavelength to be the discrepancy parameter: The super-Hubble modes during inflation will be our system and we will trace out over all sub-Hubble dofs, allowing for the possibility that some of the environment dofs will become part of the system Hilbert space due to gravitational red-shifting. In this way, we use the natural Hubble horizon provided to us by the dynamics of the system as the demarcating scale -one that is of extreme importance since it has been long argued that modes classicalize only after crossing the Hubble horizon [23][24][25][26].
When adopting such a new approach to inflation as an open EFT, several fresh perspectives are bound to come up. The one we will seek to address in this work is to what approximation are we allowed to ignore the quantum entanglement between different scales [19], and the associated dissipation, when calculating correlation functions? Functionally, in such an open EFT, one calculates the reduced density matrix of the system modes and then uses this density matrix to calculate correlation functions between operators living solely in the system Hilbert space. If indeed we were dealing with a standard QFT in flat space where our well-trusted energy conservation law was valid, then it can be shown that the low-energy density matrix, related to the vacuum state of the field theory, is equivalent to using the Wilsonian effective action [27]. However, this is precisely the point -such an effective action is not available in the case of inflation! To further stress this point, recall that gravitational interactions has this unique feature that they are universal and cannot be 'turned off' at will. Thus, unlike the quintessential example of collider experiments in high-energy physics, the interactions between the short and long-wavelength modes always remain active and thus one cannot assume the short-wavelength dofs to remain in their ground state. As shown in [28], this invalidates the standard Wilsonian argument and necessitates the generalization to Feynman-Vernon influence functionals.
In this paper, however, we will not use the path integral formulation but rather stick to the canonical version of using master equations which determine the evolution of the reduced density matrix. What we aim to do is calculate the first order quantum corrections to the power spectrum due to cubic interactions. so how does this differ from standard loop corrections? Since we have argued that the usual Wilsonian EFT formulation does not work anymore, we cannot then use the same tools to calculate radiative corrections in the standard manner. More physically, the open EFT we are introducing to the game takes into account the dissipation happening between the system and environment degrees of freedom. Of course, we still need to make approximations in order to do explicit calculation. Although abandoning energy hierarchies, one often uses the time scales associated with the system and environment. If the quantity of interest happens over a long enough time-scale such that the system can relax to equilibrium before being perturbed by the environment, then one can systematically ignore the entanglement between the system modes and the environment (even though the interaction is never turned off). We will show that when solving the master equation assuming such a 'Markovian approximation', the results, to leading order, indeed match with what one gets from the one-loop corrections in standard QFT in a quasi-dS background. This is the main finding of our work -although the leading-order quantum corrections due to one-loop effects (using the well-known 'in-in' formalism) match with those when using an open EFT approach, it does so only under the Markovian approximation. The validity of this assumption is not clearly understood for gravitational systems, and this gives us a systematic way to find corrections to this approximation in an open EFT setup. If using the standard 'in-in' formulation to calculate loop corrections to primordial correlation functions, one would be blind to such effects which arise due to the quantum entanglement between long and short dofs. And, of course, if trans-Planckian modes are ever to play a role in the quantum corrections, it will precisely be due to this fact that they cannot be safely 'disentangled' or 'renormalized' from the observable modes of interest. In this paper, we do not answer whether this is indeed the case, or not, for inflation; rather, we open a new direction for doing such calculations and show that the well-known results match only when imposing the 'Markovian' approximation.
In the next subsection, we outline the technical details of the problem at hand and give a quick overview of the existing literature on the topic. In section 2, we introduce the interaction Hamiltonian and set up our master equation in section 3. Section 4 elaborates the solution for the primordial tensor power spectrum and first order quantum corrections to it. We conclude in section 5 by discussing the significance of our result and ways to generalize it in future. Most of the technical details have been relegated to the various appendices.
1.1 Quantum corrections to the power spectrum: Review of existing literature Specifically, we are interested in calculating the first order quantum correction to the primordial tensor power spectrum due to cubic interactions between the tensor modes alone. Comparing with the standard EFT language, our result should be compared with the one-loop correction to the tensor power spectrum due to cubic tensor interactions. The reason for studying this specific type of corrections is that this arises from the purely gravitational part of the action and does not involve the slow-roll parameters for the dominant terms of interest. The cubic interaction is purely due to the non-linearity of gravity. In fact, this is the system of interest for most models of inflation which do not introduce any new rank-two massless tensor. To be clear, there can obviously be other one-loop effects to the tensor power spectrum which are bigger than the one we are considering; however, for any model of inflation, the interaction between the primordial gravity waves is going to be of this form (to leading order) as long as one has general relativity as the effective theory in mind. Contrast this with the case of scalar (or density) perturbations, whose leading order (cubic) interactions can easily be modified from the single-clock slow roll case by introducing new features in the potential or by adding higher derivative terms in the action.
In order to do a like-for-like comparison, we need contrast our result vis a vis the one loop correction to the gravitational wave power spectrum due to (cubic) pure tensor mode interactions. However, to the best of our knowledge, this calculation has not been done in the past. Nevertheless, there are similar calculations which have been done which can help guide us in our quest. Once again, we are looking for the result of such one loop corrections to the tensor power spectrum only to compare with our result and our calculation for computing the first-order quantum corrections due to the sub-Hubble (environment) tensor modes do not rely on this. Our sole motivation for this is to check how far standard loop corrections are able to capture the correct physics of inflation (when there is no longer any time-translation symmetry) when compared with the full machinery of open EFTs. In the process, we fill a gap in our current understanding of how quantum corrections from purely tensor perturbations affect the primordial tensor power spectrum.
To this end, we note that the case of loop corrections to primordial power spectra has been studied for more than a decade now. The initial motivation was indeed to test the limits of validity of inflation as an EFT and if such quantum corrections, howsoever tiny, can open new directions in understanding the physics of inflation (or, more generally, that of dS space) [1]. First, the case of loop corrections to the scalar power spectrum were studied in detail in [11,12]. It was found that loop corrections do not spoil the scale invariance of the spectrum [29]. Rather, the loop corrections to the inflationary power spectrum has a logarithmic dependence on the renormalization scale. Since then, further studies have only strengthened this result when considering different types of matter fields when coupled to the inflaton. In a similar way, the one loop correction to the tensor power spectrum has been calculated due to coupling to free fermions (both massive and massless) as well as due to coupling to other types of matter components such as isocurvature fields [30][31][32]. For our purposes, the case closest to us was studied in [32] where the author considered the one-loop correction of the tensor power spectrum due to a massless minimally-coupled scalar field. We shall return to a detailed comparison of their findings with our result in the discussion section, but let us summarize the main conclusions of loop corrections people have reached over the years with regards to inflationary perturbations.
Although it was initially thought that the scale invariance of the tree-level power spectrum will be spoiled due to a logarithmic dependence on the comoving momentum coming from one-loop effects [1,30,[33][34][35], the common consensus seems to be that such terms must cancel amongst themselves, and we would find no such enhancement to the power spectrum due to one-loop corrections [12,31,32,36]. We do not provide any new or alternative viewpoint on this and, in any case, we are not interested in computing a loop correction in this work. Rather, our goal is to provide a fresh perspective on what should be the physical way to compute quantum corrections for an open EFT such as inflation. In the beginning, as will be demonstrated later on, the calculations in the open EFT approach does seem to suggest that the perturbative quantum corrections would end up involving the comoving momentum of the subhorizon modes. However, at the very end, we do find that the final result is actually scale-invariant, as has been long speculated starting from the seminal work on [11]. What we add to this discussion, due to the nature of the new techniques employed in this work, is that how only under certain simplifying assumptions does one reach this familiar result, and how this opens up the possibility of going beyond such approximations to calculate next-to-leading order corrections for perturbative quantum effects during inflation. We will also show that this assumption of (time-local) Markovian behaviour is a good one, but not an exact one, and one needs to relax this for exploring further.

Action and mode expansions for tensor perturbations
We systematically introduce the quadratic action and the cubic interaction terms involving tensor modes alone and, in the process, set up our notation.

The quadratic action
We consider that the cosmological background is described by the flat, FLRW metric where t and τ are cosmic and conformal times respectively. For inflation, we have the usual relation: a(τ ) ∼ 1/Hτ , where H is the (approximately constant) Hubble scale during inflation. To consider the action for cosmological perturbations, we can go to the Hamiltonian (ADM) formalism: Since we are only interested in interactions purely between primordial gravitational waves, it is sufficient to consider the metric on the spatial hypersurface q ij = a 2 (δ ij + h ij ), where we have denoted tensor modes by h ij and have not considered any scalar (or vector) perturbations.
As is well known, the above prescription leads to a quadratic action for tensor perturbations which is given by [37] (primes denote derivatives with respect to conformal time τ ) where the tensor field, and its corresponding polarization tensors, are given by It is customary to work with the + and × polarizations, and throughout this work, we will mostly follow the conventions of [24]. In general, these tensors can be built from two vectors perpendicular to the momentum k in the following way, where e 1 × e 2 =k. Eqs. (83)-(85) in Appendix C show the explicit form of these tensors for three particular choices ofk. The linear equations of motion simplify significantly if we introduce the canonical variable (and also conveniently brings this problem closer to the case of scalar perturbations) such that the quadratic action becomes It reproduces the well-known fact that, at linear level, the gravitational wave action looks like two copies of that for a minimally-coupled scalar field. Consequently, the equation of motion takes the where, assuming the Bunch-Davies vacuum, the solution for the mode functions are In this way we can write the tensor field in the Heisenberg picture aŝ This expansion can be similarly found through the rotation where U 0 (τ ; τ 0 ) is the evolution operator corresponding to the free Hamiltonian which can be schematically expressed as where the quadratic Hamiltonian can easily be derived from the action (7) and T denotes timeordering. Note that the quadratic Hamiltonian is not time-independent due to a time-dependent mass (squeezing) term as discussed below.ĥ ij (x) is the Schrödinger operator given bŷ The quadratic Hamiltonian, corresponding to the action (7), has two pieces -the standard Hamiltonian for free scalar fields in flat space and a term which couples k with −k, known as the squeezing operator. These squeezing interaction can be thought of as a time-dependent mass term which result in a gravitational pump sourcing zero-momentum correlated pairs. The unitary evolution operators depend on the so-called squeezing parameters and are defined in the same fashion as for scalar perturbations [38]. Although all of these things are well-known, we repeat them to emphasize that the evolution generated by the quadratic Hamiltonian do not lead to any mode-coupling, and therefore, they do not contribute to the entanglement between modes of different wavelengths. Only cubic interactions, which we introduce in the next subsection, will be responsible for such couplings.

Cubic interaction Hamiltonian
The cubic action which describes all the self-interactions of primordial tensor modes is given by [24,39] where we have introduced the slow-roll parameter ǫ := −Ḣ/H 2 , and H := a ′ /a = aH, where H denotes the (approximately) constant Hubble parameter. From the above action, we can derive the interaction Hamiltonian as One thing which is important to point out here is that there are terms in the above cubic action which are independent of any 'so-called' slow-roll parameters. This is so because there are tensor mode perturbations even for pure dS whereas scalar perturbations can only appear for a quasi-dS background. 2 This fact implies some simplifications regarding gauge issues when dealing with pertubative quantum corrections to the tensor power spectrum, when considering tensor interactions alone. Next, we Fourier expand the operator above and go to the interaction picture, as it will simplify much of the computations. The standard way to pass over from the Schrödinger picture to the interaction one is given by is the unitary evolution operator corresponding to the quadratic Hamiltonian introduced in Eqn. (12). Following these steps, we obtain where we have defined and The tensor fields are expanded in terms of the creation and annihilation operators introduced in Eqn. (10), so we have all the ingredients ready for setting up the master equation for the super-Hubble modes. Before moving on, note the following two properties. Firstly, the interaction Hamiltonian, as written in Eqn. (16), is manifestly hermitian. Furthermore, we only consider the third term on the r.h.s. of Eqn. (15) since the others are subdominant. The reason is twofold. First, the chosen term is the only one with no time or spatial derivatives. Thus, requirements like slow-roll or background homogeneity pose no constraints on it. Moreover, momenta factors coming from the spatial derivatives become more relevant in the UV, and so it is expected that their contribution when involving IR modes (S) will be suppressed in comparison to terms with no such derivatives.

Master Equation and Lindblad Operators
Before writing down the master equation governing the evolution of the reduced density matrix corresponding to the super-Hubble modes, let us introduce some necessary notation. First, unless stated otherwise (through a subscript), every operator is written in the Schrödinger picture. Likewise, the subscripts E and S indicate that operators are acting on the environment or system Hilbert space, respectively. If no subscript is pointed out, then the operator acts on the entire Hilbert space. Notice that for any operator which spans the entire Hilbert space, we have the decomposition where the sum is over any possible combination of system and environment modes associated to the Fourier expansion of the operator. Naturally, this expansion is going to be very messy, but the resulting expressions can be reduced by relabeling momenta and using other symmetries (related to permutations of the indices), similar to what was done in Eqn. (16). A full derivation of the master equation is shown in Appendix A, where we find that where ρ where V := V Schr denotes an interaction term in the Schrödinger picture and . The sum over i is meant to schematically represent the tracing out of all the environment degrees of freedom (and would, more explicitly, correspond to integrating out of sub-Hubble momentum modes). L 1 and L 2 are the so-called Lindblad operators. The vectors |E 0 denote the initial environment state, which we take to coincide with the Bunch-Davies vacuum |0 . Likewise, |S 0 stands for the initial state of the system modes, which are also Bunch-Davies states for modes that start inside the horizon but that eventually cross it to form part of the system (at a given time τ ). On the other hand, |E i and |S(τ ) denote the environment and system states at time τ , respectively.
There are a couple of important things we should point out at this stage. Most crucially, note that the form of our equation automatically highlights the main assumption underlying our work -Markovian approximation. In general, a master equation for a reduced density matrix should take the following form [41]: The above expression is very similar to our master equation Eqn. (21) except that now we have allowed for arbitrary dissipation coefficients γ(t), which represent the exchange of energy between the system and environment and is the quintessential marker for the non-unitary evolution of the reduced density matrix. What is crucial here is that under the Markovian approximation, one finds numerical values for γ implying that they are constant and γ > 0. This is indeed the case for us and shows how the Markovian approximation is built into our formalism. In fact, we are considering the simplest case where although we are allowing for non-unitary dynamics of the super-Hubble modes (as is typical for open EFTs), we require time-independence of the dissipation coefficients [42]. Even for a time-local, Markovian system, one might allow for the γ's to be time-dependent, and this would imply γ(τ ) > 0 at all times. What we have shown in Appendix C is that our assumption of the Markovian evolution is, at least, self-consistent. There, we have evaluated the dissipation kernel corresponding to our system and shown that it is indeed sharply-peaked, thereby justifying our time-local approximation [43,44].
Of course, the more interesting and general case would be to consider a time nonlocal master equation but more on this will be discussed in Section 5.
Another conceptual point we wish to clarify is the choice of the basis states |E i , over which we trace out at some time τ . As has been emphasized in [24], these are not the same as the state one gets from evolving the initial environment state under the free evolution operator (corresponding to the sub-Hubble modes), i.e.
Rather, at a given time τ , we identify a complete set of basis states |E i , over which we take the trace of the environment modes. Therefore, these are not time-dependent (unlike the right hand side of the above equation). On the other hand, such a choice of a basis for the sub-Hubble modes at a given time τ indicates that there will be considerable overlap between |E 0 and |E i . The subtleties of this choice for tracing out the environment dofs have been described in detail in [24,45].

Perturbative solution
In order to derive the perturbative solution of the reduced density matrix, we refer to Eqn. (56) in Appendix A, and then trace over the environment degrees of freedom, in the same fashion as in Eqn. (58). Since Eqn. (56) is given in the interaction picture, we must also multiply each side of it by the free theory unitary operator and its conjugate, as in ρ(τ ) = U 0 (τ ; τ 0 )ρ I (τ )U † 0 (τ ; τ 0 ). First, let us work at the zeroth-order approximation in some detail, which will help us learn some of the technicalities needed when we go to second order. Using the fact that, at this order where we have used U 0 (τ ; τ 0 ) = U E (τ ; τ 0 ) U S (τ ; τ 0 ), and the fact that the states |E i form a complete basis. Further, notice that we have recovered the equation ρ (0) red (τ ) = ρ S (τ ), which was used in the previous section to find the master equation.
Before going to the second-order solution, let us warn the reader we will relax our notation a bit from hereon. As mentioned before, any given operator can be expanded in terms of its sub-horizon and super-horizon Fourier modes asÂ = iÂ i E ⊗Â i S . Henceforth, we shall omit the sum and represent every operator asÂ =Â E ⊗Â S . Finally, we introduce the functions G I (τ ) = dτ ′ V I (τ ′ ), where the integrands are either h 0 (τ ) or h 1 (τ ), depending on the combination of creation and annihilation operators (as shown in Eqn. (16)) in action. With these considerations, schematically the correction to the reduced density matrix at second order is given by

Power Spectrum
In this section we will compute the corrections to the tensor power spectrum due to the evolution of the initial pure state into a mixed one, signaling the generation of entanglement entropy [19].

Linear power spectrum
As before, we start with the zeroth-order approximation, showing that the present formalism is obviously consistent with the expected outcome. First, the tensor power spectrum is defined to be where the two-point function at zeroth order is computed as follows In this way, obtain as expected at late times τ → 0. We show this simple derivation in detail to explain to the reader what we operationally mean by taking the trace over the environment dofs, as it will be done later on for the first order quantum correction in a more complicated context.

Corrections to the power spectrum: Seeting up the problem
Firstly, for the sake of brevity we abbreviate the product of tensor modes aŝ which takes the following form in the interaction picture: Even though we are working in the Schrödinger picture, the resulting combination of unitary operators will takeÔ q intoÔ q (τ ). Furthermore, notice that the creation and annihilation operators act on the initial (Bunch-Davies) vacuum. Considering this, the correction to the tensor spectrum is of the form In principle, we now must deal with every possible combination of system and environment modes in both V I and G I . Nevertheless, we must have at least one system mode, otherwise we are missing the effect of the environment on the system. Then, we would have to work with any possible combination of (EES) and (ESS) modes. Notice that due to momentum conservation (imposed by the delta functions), the sum of the three momenta must be 0. This constrains the contribution from the (ESS) combination (which can be thought of being of the 'folded' shape), and is thus subdominant in comparison to the 'squeezed limit' contribution (EES).

Inner products
The resulting scalar products in Eqn. (33) determine which combinations of creation and annihilation operators contribute to the corrections to the spectrum. Both V I and G I carry every possible combination of the product of three of such operators (ĉ andĉ † ), as shown in Eqn. (16) for V I and its time integration for G I . However, since we have to deal with scalar products on different Hilbert spaces (H E for the environment and H S for the system), the creation and annihilation operators must act on each space depending on the magnitude of the momentum associated with them. For notational convenience, whenĉ (ĉ † ) acts on H E we shall denote it byb (b † ), and byâ (â † ) when it acts on H S . In this way, the one-particle states in the system and environment at the time τ are respectively given byâ α † Looking at the scalar product on H E in Eqn. (33), we can see that the combinations that contribute to the product must be V I (τ ′ ) ∼bbâ or ∼bbâ † , and G I ∼b †b †â or ∼b †b †â † . 3 From the expression forÔ q (τ ) and from the trace, we can see that the kind of terms we need to compute are of the form ∆k ∆p where we have fixed k 1 and p 1 for the system modes and the remaining momenta for the environment modes. In the end, we will consider every case by multiplying the result by an appropriate multiplicity factor. Next, let us tackle each term above separately, where the first one is simply given by The second term does not present a 'contraction' between q and k 1 , and thus we've got The final term is the most complicated one, so we break it in parts as follows Notice that, from the last equality, we can combine one factor from the first line, together with the second, third and last lines to cancel out with the expression from the second term above, Eqn. (35), and thus we end up with where 18 is the multiplicity factor (3 ways of choosing the system mode in h 0 , 3 factors from h 1 , and 2 ways of contracting the inner product over environment modes).

Contribution of different polarizations to the integral
Each product above contains complicated sums over polarizations, for which we show the results below. Notice that due to the state (scalar) products, the first momenta on the left has the same polarization as the first on the right, and so on. In consequence, the sums of interest are, upon integration over φ, 4q 4 + 4k 4 2 + 11q 2 k 2 2 + 8qk 2 (q 2 + k 2 2 ) cos θ + q 2 k 2 2 cos(2θ) sin 4 θ (q 2 + k 2 2 + 2qk 2 cos θ) 2 where without loss of generality we have considered q = q(0, 0, 1) = ±k 1 , k 2 = k 2 (sin θ cos φ, sin θ sin φ, cos θ) , and due to momentum conservation k 3 = −(k 1 + k 2 ). Notice that the momentum we use can be ±k 3 depending on whether we used a creation or an annihilation operator. Then, we need to solve integrals of the type (say, for the third term on Eqn.(36)),

Integration over time and momenta
The correction to the spectrum follows from the rather complicated integral shown above and other terms which are similarly very involved. In principle, one can try and play with the order of integrations to facilitate the task ahead, although special care has to be put into the integration limits to be consistent with which modes belong to the environment and at which times. So, in principle, we could perform first the integration over the environment mode k 2 , which allows to define a dissipation kernel in an approach similar to [15,46]. We try and follow this path in Appendix C, which leads to a dead-end due to the complicated form of the remaining integrals (however, the kernel can provide useful information about the physics of the system and, in particular, show its time locality). To proceed further in that vein would require numerical solutions. Here, we shall follow a more fruitful path, where we start with the integration over τ ′′ and be able to evaluate our result analytically. In order to deal with the annoying pre-factors of the mode functions, let us define β k (τ ) := √ 2k v k (τ ), such that the integral takes the form where Ξ is the function resulting from the integration over τ ′′ , and is given by with y := cos θ. Next, we introduce the dimensionless variables ω ′ := −qτ ′ and κ i := −k i τ ′ , such that the integral Eqn. (42) becomes Then, we Taylor expand the entire integrand around κ 2 → ∞, which renders This is a good approximation for large values of κ 2 , and a reasonable one for κ 2 ∼ 1. Then, we can integrate over y, κ 2 and ω ′ (in that order), to find Pl where we have also included the complex conjugate term. The terms on the first line come from the upper limit of the κ 2 momentum integration, where aΛ represents a comoving cutoff for k 2 which translates into Λ/H for κ 2 (= −k 2 /(a(τ ′ )H), whereas the terms on the second line come from κ 2 → 1 (for k 2 , this lower limit would have been aH). Notice that the UV-divergences (going as ln Λ) have been canceled out due to adding in the complex conjugate (they appear as a purely imaginary contribution). However, that is not the case for the divergences for all the terms, for instance, the logarithmic ones coming from both terms on the second line of Eqn. (36). We have more to say on this later.

Quantum corrections
We have now all the ingredients to compute the final result describing the quantum corrections to the tensor power spectrum. Let us state the final result now and then discuss some of the intermediate steps as well as the consequences in the next subsection.
Following the same prescription as before, we can compute the contribution to the spectrum from each term in Eqn. (36), including their complex conjugates. In doing so, we find that the expressions coming from the limit κ 2 = 1 cancel out up to order O(ω) (at least). In other words, our final result depends on positive powers of ω which, in the late time limit τ → 0, goes to zero. This is quite remarkable since, looking at Eqn. (46), one would have been skeptical that the terms with inverse powers of ω might have survived and that would have been disastrous. However, they all cancel out neatly in the end.
Then, the only surviving terms are those associated with the UV-divergences, and a rather innocuous IR one as well. We will discuss these terms next, but for now our final result for the leading order correction to the spectrum is given by: where Ci is the cosine integral function, and we have not explicitly written some O(1) numerical factors above (involving the Euler-Mascheroni constant γ E and other such terms) which depend on the renormalization scheme.

UV and IR divergences
Admittedly, we have pulled our final result for the first order quantum correction to the (dimensionless) tensor power spectrum, in Eqn. (47), rather out of the hat. The gruesome details of each of the integrals in Eqn. (36) have been shown in Appendix B. Let us elaborate on some of the intermediate steps now. Firstly, note that we have multiplied by a factor of q 3 to get the dimensionless power spectrum. A factor of q 2 then combines with the factor of a 2 from the prefactor to give us a multiplicative factor of ω 2 H 2 while one factor of the leftover q had been used to convert the measure from dq to dω. After evaluating the integrals, the results can have two types of divergences as usual -UV and IR ones. To understand them fully, look at our final expressions for the four terms we get from performing these integrals, given in Eqns. (68) -(71) in Appendix B.
For the last two terms Eqns. (70) and (71), the beautiful cancellations happen as shown in Eqn. (46) above. The UV divergent terms appear as i ln Λ, and therefore, fall out of the final expression due to adding in its complex conjugate. 4 This is not strange at all -all this suggests is that although we need to put in a cutoff to regulate momentum integrals in the intermediate steps, the physical result is independent of such cutoffs simply due to how the UV divergences appears in them.
On the other hand, there are a plethora of inverse ω terms remaining in these two expressions which would lead to serious IR divergences when taking the late time limit of τ → 0. As mentioned in passing above for Eqn. (46), all of these terms cancel and we have a resulting expression which is O(ω 2 ) or higher for these terms. Note that the fact that these terms drop out only in the qτ → 0 limit, and not otherwise, does not tell us anything deep at all. Of course, even the linear power spectrum is only scale-invariant in this limit, which is the relevant one for observations. At the risk of being pedantic, we repeat this trivial argument to convince the reader that nothing strange in going on for the quantum correction in this aspect of it. In conclusion, these two terms (70) and (71) do not have any UV or IR divergences at all.
Let us now come to the more interesting terms given in Eqns. (68) and (69). These terms can be represented as giving a correction of the following generic form (we have omitted a constant prefactor to simplify the expression): where we have only written down the terms which survive after several cancellations amongst themselves from Eqns. (68), (69), (70) and (71) in Appendix -B. Let us first talk about the UV sensitivity of the above term and deal with the ln(2ω) term later. The UV divergence is logarithmic and of the form ln (Λ/H). 5 This is what is commonly expected in loop corrections for inflationary perturbations.
The common way to deal with them is separate out a divergent piece from the above term and keep the finite contribution, for some 'renormalization scale' at which the observations are made [11,29,32]. A more rigorous way to achieve the same would be to instead use a dimensional regularization scheme instead of the comoving cut-off used here. However, there is a simple way to formalize the equivalence between the dimensional and cutoff regularizations [47], This equivalence, proven in the context of Gauge theories, shows that diffeomorphism invariance is respected when a comoving cutoff has been introduced in the calculation. The divergent part is isolated ε → 0 and can be dispensed with. Up to finite factors that can be absorbed into the divergent part, and up to O(1) constants depending on the renormalization scheme, the correction to tensor spectrum now reads The introduction of the renormalization scale is quite well-known, and has been discussed in great detail in [11,29,32,36]. Instead to rehashing those arguments, let us point out that the crucial observation regarding the UV behaviour of our calculation goes as follows. Had we found a term which goes as ln(Λ/q) [2,30,33,34], that would have been a truly problematic term. As was first pointed out in [11], and has been repeated many times since [29,48], these type of terms are simply ruled out by the diffeomorphism symmetry of the problem. However, choosing the correct covariant regularization scheme led us to find a term of the form ln(Λ/H) and this is indeed what is expected from loop corrections.
However, there is another type of large logarithms appearing in Eqn. (48) coming from the terms Eqns. (68) and (69) in Appendix B. This is the term proportional to ln(|qτ |), which diverges in the IR, when taking the late time limit τ → 0. As mentioned, this is not a UV divergent term, but rather an IR one. 6 These type of divergences first appeared in [49][50][51] but have been since also shown to be harmless in the sense that they can never affect local observations [29,52,53]. The essence of the argument is that and the long wavelength can be reabsorbed, through the remaining gauge freedom, by change of coordinates (see [36] for a nice review of these arguments). In [12], a similar logic was used in the context of scalar perturbations, where the time dependence induced by a cubic interaction (with spatial derivatives) is suppressed through the effects of a quartic interaction necessary to maintain diffeomorphism invariance at late times. Finally, this IR growth appearing in the quantum correction to the tensor power spectrum can be related to the issue of the IR divergence of the Bunch-Davies vacuum for the free graviton mode [54] and is certainly not relevant for any local observable.
In conclusion, what we find is that taking the limit qτ → 0 does not lead to any divergence at all. This is all the more surprising since, naively, it did look like the perturbative quantum corrections (e.g. see Eqn. (46)) would give large corrections in this limit. However, the fact that the various terms in Eqn. (47) nicely cancel out is not an accident at all and, we believe, is sensitively dependent on our Markovian approximation. We will explain this is more detail in the next section.

Interpretation of our result
The most obvious conclusion to be drawn from our computation is that the leading order quantum corrections to the tensor power spectrum matches exactly with what is expected from loop corrections, when considering purely cubic tensor interactions. Let us elaborate on this statement a bit. Firstly, just from power counting one expects that the one-loop effect will be suppressed by a factor of O(H 2 /M 2 Pl ). Furthermore, more recent computations of such loop corrections to the tensor power spectrum (not from tensor loops but rather from other forms of matter loops) indicate that we should expect a correction of the form [31,32] where µ is some renormalization scale. Note that earlier works predicted much larger corrections where the logarithm would contain a term of the form log (k/µ) [30]. These type of large logarithmic corrections have since been ruled out on the basis of symmetry, and also correctly including different terms which should contribute to the same order [11,36].
Thus, our result is consistent with this broad consensus that the quantum corrections calculated in an open EFT approach does not contain any large logarithmic factors (coming from the UV) which also spoils the scale invariance of the spectrum. In fact, in the late time limit, which is what one should naturally consider for evaluating cosmological correlation functions, our result shows that the leading order quantum corrections exactly cancel. There is no extra time-dependence introduced in the tensor power spectrum coming from quantum corrections due to logarithmic terms depending on the comoving momentum, which has been the main finding of loop corrections [29].
For the benefit of the reader, let us rewrite our final result in the form: with a negative constant C and some O(1) constant C 2 . Similar logarithmic terms have also appeared for loop corrections to the tensor power spectrum when considering massless isocurvature fields coupled to tensor modes, resulting in loops arising from Dirac fermions, minimally coupled scalars or Gauge fields [32], or for loop corrections to the tensor power spectrum due to interactions with a conformal scalar field [31]. Moreover, as seen in [31], there can be additive constants appearing in the one-loop correction, in addition to the ln(H/µ) term, which depends on the renormalziation scheme.
This is exactly what we find in our case as well. Another important similarity is that the prefactor C of the logarithmic term is negative for all the different loop corrections considered in [31,32], and so is the case for us. Thus, it shows that our results are indeed along the expected lines of those coming from one loop corrections to the tensor power spectrum.
Having said that, keep in mind that the interaction terms for the aforementioned systems are determined by the coupling of the graviton to these matter fields, and not by pure GR alone (see e.g. [55] for a discussion on quantum corrections in GR as an EFT). This is the one of the novelties of our work. The cubic interaction of the tensor modes, which leads to the quantum corrections to the power spectrum, is completely specified by the nonlinearity of GR and is free from choices of the coupling. In fact, this is why our leading terms in the interaction Hamiltonian can be free of both slow-roll parameters as well as any derivative couplings. Since we were unable to find the corresponding calculation for loop corrections to the tensor power spectrum, arising from cubic tensor interactions alone, it is impossible to do an exact comparison of our computation with standard loop corrections evaluated in the in-in formalism. Nevertheless, this highlights a salient feature of our work: Even without going into the merits of abandoning standard loop corrections in favor of using an open EFT approach to compute quantum corrections to inflationary n-point functions, which is the main purpose of this work, we have managed to fill the lacuna in the existing literature by calculating quantum corrections to the gravitational wave spectrum due to tensor interactions alone.

Looking ahead: Relaxing the Markovian approximation
As is often done when exploring new techniques, we have verified that our result matches with standard loop corrections to the tensor power spectrum and have the same logarithmic dependence on the Hubble parameter. Although this is a nice sanity check to make sure that our computations are correct to the leading order, the whole point of undertaking this exercise is to be able to go beyond loop corrections which necessarily have Wilsonian EFT underpinning its technical implementation. As we have emphasized several times, the late-time limit of the first-order quantum correction calculated here is under the strict Markovian approximation which leads to our master equation having the standard Lindblad form. Nevertheless, this is an assumption baked into our analyses -one which we find to be self-consistent (see Appendix C).
However, this does not imply that this assumption must necessarily be true for cosmological dynamics of the form considered here. Counterexamples of cosmological setups displaying non-Markovian behaviour have been recently demonstrated in [18,56]. Therefore, our goal is to next consider a more general master equation for the superhorizon modes that is time-nonlocal (of the Nakajima-Zwanzig form) [43,57,58]. In fact, this can be seen from our calculation as well. One finds that the dissipation kernel is sharply peaked for our model (Appendix B) insofar that it goes as 1/(τ − τ ′ ) 6 (see Eqn. (89)). Nevertheless, exact Markovian behavious would demand a delta-function peakedness which obviously is not the case for us. In other words, this shows how a systematic study can now be undertaken to explore the deviation from Markovianity for our system by allowing for more general master equations.
Therefore, the general trend seems to suggest that the relevant question to ask would be the following: What is the regime in which we are interested in calculating quantum correction to the power spectrum due to cubic (and other higher-point) interactions? This is so because in certain regimes, the Nakajima-Zwanzig master equation does indeed reduce to the standard Markovian-Lindblad form, and if it so happens for the observational questions of interest, then we can conclude that our results are robust and our Markovian approximation justified. Indeed, this would imply that standard loop corrections implemented to calculate radiative corrections to the power spectrum of inflationary perturbations do capture the essential physics of these systems. However, this is yet to be proven and a systematic study of these systems have to be be undertaken in order to either prove the above statement or falsify it.
But the story does not end here. Apart from observational consequences, there are other conceptual lessons to be learnt from applying open EFT techniques to calculating such quantum corrections to the inflationary power spectra. Let us return to the entanglement of UV-modes to the IR ones. It might indeed be possible that the trans-Planckian modes are decoupled from observational dofs, even allowing for the non-Wilsonian character of the system. 7 However, if it is found that non-Markovian effects are small for observable modes and can yet become large for another subset of modes of the system, these would lead to deep ramifications for inflation (e.g. for eternal inflation) [59] and even for QFT in dS spacetimes (e.g. for the lifetime of metastable dS) [60][61][62]. The intriguing aspects of non-Markovian nature of gravitational interactions are yet to explored at all and we have simply set up the mechanism to systematically investigating this question in more detail in the future.

Acknowledgements:
The authors thank Thomas Colas for feedback on an earlier version of this draft. SB is supported in part by the Higgs Fellowship. AB is partially supported by STFC. JCF is supported by the Secretary of Higher Education, Science, Technology and Innovation of Ecuador (SENESCYT).

A.1 Notation
Before deriving the master equation, let us introduce some notation to avoid confusion. First, unless otherwise stated (through a subscript), every operator is written in the Schrödinger picture. Likewise, the subscripts E and S indicate that operators act on the environment or system Hilbert space, respectively. If no subscript is pointed out, the operator acts on the entire Hilbert space. Notice that for any operator like that we haveÂ = where the sum is over any possible combination of system and environment modes associated with the Fourier expansion of the operator. Naturally, this expansion is quite lengthy, but the resulting expressions can be reduced by relabeling momenta and using other symmetries, similarly to the process followed to arrive to Eqn. (16).

A.2 Building the Master Equation
In this section, we sketch the derivation of the master equation in terms of the so-called Lindblad operators. For this, we will work in the Schrödinger picture, although it is convenient to start from the von-Neumann equation in the interaction picture, where ρ I is the density matrix in the interaction picture. Next, we will re-write the equation above in a way more suitable for the upcoming approximations. For this, notice that the solution is given by The equation above is the full solution to the von-Neumann equation. We intend to go to second-order approximation, for which we have Now, it is turn to go to the Schrödinger picture. To do so, we use ρ(τ ) = U 0 (τ ; τ 0 )ρ I (τ )U † 0 (τ ; τ 0 ), so that where we have taken the time-derivative of (56). Next, in order to find the master equation governing the evolution of the reduced density matrix, we trace over the environment dofs as follows where |E i denote the environment states, in our case the sub-Hubble modes at time τ . Performing this operation on both sides of (57) we obtain where and V I (τ ′ − τ ) = U † 0 (τ ′ ; τ )V U 0 (τ ′ ; τ ) . L 1 and L 2 are the so-called Lindblad operators. In order to see how Eqn. (59) comes about, notice that every term containing an H 0 factor in (57) contributes to the commutator of the quadratic Hamiltonian in (59). Regarding the effective Hamiltonians, the first one comes directly from the third term on the first line, on the r.h.s. of (57). The second effective Hamiltonian comes (partially) from the last term on (57). To see this, notice It can be seen that in order to construct V eff2 we summed and subtracted 1 2 i L † 2 L 1 , ρ S (τ ) to the expression, which renders the formula above. In this way, the last term encompasses the non-unitary part of the evolution, as the remaining terms can be written as a commutator with a Hamiltonian, which can be seen as unitary evolution.
Then, upon the introduction of the dimensionless variables ω ′ and κ i , the integrals we need to solve are of the type where i is an index referring to the two functions f 1 (θ) and f 2 (θ) -Eqs. (37) and (40)-, which we found by summing over the polarisation tensors present in Eqn. (36). Next, as explained in the main text, we expand the integrand around κ 2 → ∞, which is the region that dominates the integral. We have checked that the same approximation gives a good order-of-magnitude approximation for the region near κ 2 = 1. However, and rather conveniently, the contributions from that limit cancel out, leaving us with the UV-divergent terms.

C Dissipation Kernel
In principle, one can also find the master equation in terms of a dissipation kernel, similar to the process followed by [15,46]. Such an approach proved to be less advantageous than the one used in this work when trying to find the corrections to the spectrum, although one can still find a kernel characterizing the dynamics of the environment in its interaction with the system. To do so, consider the interaction of the form where we have chosen the same combination of environment and system degrees of freedom as before.
Next, we re-write the interaction Hamiltonian in terms of the Bunch-Davies mode functions, such thatV where we have introduced the following operator for future convenience Finally, from here onwards, we denote the coupling as λ(τ ) := 4 √ 2H/(τ M Pl ).

C.1 Master equation
First, we shall work in the interaction picture in contrast to the process followed in the main text, so every operator should be understood to belong to that picture. With this consideration in mind, the density matrix is given by where Notice that environment and system degrees of freedom are defined exactly as before, i.e., they are delimited by the Hubble length. Furthermore, since the quadratic Hamiltonian governs the evolution of states in the interaction picture, the factorisation of the density matrix holds at later times. Moreover, the weak coupling between system and environment keeps the latter unperturbed, such that ρ E (τ ) ≈ ρ E (τ 0 ). Then, tracing over the environment degrees of freedom on Eqn. (56), we get e + (k 2 ) =    cos 2 θ cos 2 φ − sin 2 φ (1 + cos 2 θ) sin φ cos φ − sin θ cos θ cos φ (1 + cos 2 θ) sin φ cos φ cos 2 θ sin 2 φ − cos 2 θ − sin θ cos θ sin φ − sin θ cos θ cos φ − sin θ cos θ sin φ sin 2 θ