# Benchmarking the cosmological master equations

Thomas Colas, Julien Grain, Vincent Vennin

► To cite this version:

Thomas Colas, Julien Grain, Vincent Vennin. Benchmarking the cosmological master equations. Eur.Phys.J.C, 2022, 82 (12), pp.1085. 10.1140/epjc/s10052-022-11047-9. hal-03782906

HAL Id: hal-03782906

https://hal.science/hal-03782906

Submitted on 22 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Benchmarking the cosmological master equations

Thomas Colas<sup>1,2,3,a</sup>, Julien Grain<sup>1,2,b</sup>, Vincent Vennin<sup>2,3,c</sup>

¹ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France

² Laboratoire Astroparticule et Cosmologie, CNRS Université Paris Cité, 10 rue Alice Domon et Léonie Duquet, 75013 Paris, France

³ Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 75005 Paris, France

Received: 7 October 2022 / Accepted: 16 November 2022 / Published online: 1 December 2022
© The Author(s) 2022

**Abstract** Master equations are commonly employed in cosmology to model the effect of additional degrees of freedom, treated as an “environment”, onto a given “system”. However, they rely on assumptions that are not necessarily satisfied in cosmology, where the environment may be out of equilibrium and the background is dynamical. In this work, we apply the master-equation program to a model that is exactly solvable, and which consists of two linearly coupled scalar fields evolving on a cosmological background. The light field plays the role of the system and the heavy field is the environment. By comparing the exact solution to the output of the master equation, we can critically assess its performance. We find that the master equation exhibits a set of “spurious” terms that explicitly depend on the initial conditions, and which arise as a consequence of working on a dynamical background. Although they cancel out in the perturbative limit of the theory (i.e. at leading orders in the interaction strength), they spoil resummation. However, when those terms are removed, the master equation performs impressively well to reproduce the power spectra and the amount of the decoherence of the light field, even in the strongly decohered regime. We conclude that master equations are able to perform late-time resummation, even though the system is far from the Markovian limit, provided spurious contributions are suppressed.

## Contents

1 Introduction .................................... 2

2 The master-equation bestiary .................... 3

2.1 An exact master equation: the Nakajima-Zwanzig equation .................................... 3

2.2 Born approximation: the time-convolutionless cumulant expansion ................................ 4

2.3 Markovian approximation: the Lindblad equation ........................................ 5

2.4 Link with perturbative methods ................................................................. 6

3 Curved-space Caldeira-Leggett model ................................................................. 6

3.1 Exact description ........................................................................................................ 7

3.2 Effective description: the TCL<sub>2</sub> master equation .................................................... 9
Interpretating the master equation ................................................................. 9

3.3 Transport equations .................................................................................................... 11
TCL<sub>2</sub> transport equation .................................................................................... 11
Exact transport equation .................................................................................... 11
SPT transport equation .................................................................................... 11

3.4 Spurious terms ........................................................................................................ 12

4 Non-perturbative resummation ................................................................. 13

4.1 Power spectra ........................................................................................................ 13
Growth rate ........................................................................................................ 14
Relative deviation to the exact result ................................................................. 14

4.2 Decoherence ........................................................................................................ 15

5 Conclusion .................................................................................................................... 16

Appendix A: Microphysical derivation of the TCL<sub>2</sub> master equation .................................... 17

Appendix B: Phase-space representation of the TCL<sub>2</sub> master equation .................................... 19

Appendix C: Coefficients of the transport equation for TCL<sub>2</sub> .................................... 19

C.1 Exact results ........................................................................................................ 20

D<sub>11</sub> coefficient ........................................................................................................ 20

D<sub>12</sub> coefficient ........................................................................................................ 21

Δ<sub>11</sub> coefficient ........................................................................................................ 21

Δ<sub>12</sub> coefficient ........................................................................................................ 22

C.2 Sub-Hubble limit ........................................................................................................ 22

C.3 Super-Hubble limit ........................................................................................................ 22

D<sub>11</sub> coefficient ........................................................................................................ 23

<footnote><sup>a</sup> e-mail: thomas.colas@universite-paris-saclay.fr (corresponding author)</footnote>

<footnote><sup>b</sup> e-mail: julien.grain@universite-paris-saclay.fr</footnote>

<footnote><sup>c</sup> e-mail: vincent.vennin@ens.fr</footnote>
<table><tr><td>D<sub>12</sub> coefficient</td><td>24</td></tr><tr><td>Δ<sub>11</sub> coefficient</td><td>24</td></tr><tr><td>Δ<sub>12</sub> coefficient</td><td>24</td></tr><tr><td>Appendix D: Comparison between TCL and perturbation theory in the curved-space Caldeira-Leggett model</td><td>25</td></tr><tr><td>D.1 Perturbation theory</td><td>25</td></tr><tr><td>D.2 Perturbative solution of TCL</td><td>27</td></tr><tr><td>Appendix E: Comparison with other late-time resummation techniques</td><td>28</td></tr><tr><td>References</td><td>30</td></tr></table>

# 1 Introduction

According to the standard model of cosmology, all structures in our universe emerge from the gravitational amplification of vacuum quantum fluctuations at early times. This idea is supported by the data, e.g. the measurements of the cosmic microwave background anisotropies [1], which reveal that primordial fluctuations are almost scale invariant, quasi Gaussian and adiabatic. Those observations are consistent with a phase of primordial inflation, driven by a single scalar field along a smooth potential.

However, most physical setups that have been proposed to embed inflation contain a large number of additional degrees of freedom [2]. Even if they provide negligible contributions to the dynamics of the universe expansion, they may affect the emergence of cosmic structures in various ways. For instance, they could lead to entropic fluctuations, or to deviations from Gaussian statistics, that future cosmological surveys might be able to detect [3,4]. They may also contribute to processes occurring after inflation (such as the production of curvature perturbations [5], dark matter [6], or dark energy [7,8]) but that crucially depend on the way those extra fields are excited during inflation. At the more fundamental level, additional degrees of freedom may also alter the quantum state in which primordial density fluctuations are placed, in particular through the mechanism of decoherence [9–18]. Decoherence [19–21] is usually associated with the erasure of genuine quantum signatures so this may affect our ability to prove or disprove that cosmic structures are of quantum-mechanical origin [22,23].

For those reasons, it has become of increasing importance to design reliable tools to model the presence of additional degrees of freedom in the early universe [24–40]. One such approach is the so-called master equation program (see for instance Refs. [41,42]), where an effective equation of motion is obtained for the reduced density matrix of a “system” of interest, once the degrees of freedom contained in the “environment” have been traced out. One of its appealing advantages is its ability to resum late-time secular effects [43–49], hence to go beyond standard perturbation theory

and implement non-perturbative resummations in cosmology.

However, master equations were primarily developed in the context of quantum optics, so they rely on assumptions (e.g. that the environment comprises a large reservoir in thermal equilibrium) that are not necessarily satisfied in cosmology. There, since the background is dynamical, the Hamiltonian is time-dependent [50] and the environment is generally out-of-equilibrium [51]. This is why, in this work, we want to understand under which conditions the master-equation program can be employed in cosmology, and what physical insight one shall expect to get out of it.

We address this issue by considering a toy model that is exactly solvable, such that the output of master equations can be compared to the exact result and examined in a critical way. This allows us to benchmark master equations. In practice, we consider two linearly coupled scalar field evolving on a homogeneous and isotropic universe. The model has been solved exactly in Refs. [52,53], where it has been shown that each Fourier sector is placed in a four-mode squeezed state, which is a Gaussian state. By tracing over the heaviest field, one obtains the reduced state of the lightest field, which follows a non-unitary evolution, and which can be compared with the predictions of different approaches, such as master equations or standard perturbative techniques. In this model, the environment does not reach thermal equilibrium, and as we will show the Markovian limit [54] is not attained either. This is why it is a priori challenging for conventional master-equation approaches to properly describe its dynamics.

The rest of this article is organised as follows. In Sect. 2, we introduce the master-equation formalism, and clarify the levels at which the different approximations enter the calculation. In Sect. 3, we introduce the cosmological model mentioned above, and show how it can be solved exactly. We then apply the master-equation program to this setting, and find that it exhibits a set of terms that we dub “spurious”. These terms do not exist in the perturbative limit of the theory, and they prevent resummation due to their dependence on the initial conditions. In Sect. 4 we then analyse the ability of the master equation to reproduce the power spectra of the model, as well as to predict the amount of quantum decoherence, when spurious terms are removed “by hand”. We find that master equations are impressively efficient in that case, even in the strongly decohered regime, and that they perform much better than standard perturbative methods (such as e.g. the in-in formalism). This also leads us to draw a few conclusions as to whether a heavy scalar field can efficiently decohere cosmological fluctuations. In Sect. 5, we summarise our main findings and further discuss the status of the spurious terms. The paper ends by a few technical appendices, to which we defer the derivation of some of the results given in the main text.
## 2 The master-equation bestiary

The master-equation program proposes to describe the quantum state of a system when it weakly couples to an environment. In practice, one considers a Hamiltonian of the form

$$ \hat{\mathcal{H}} = \hat{\mathcal{H}}_S + \hat{\mathcal{H}}_E + g\hat{\mathcal{H}}_{\text{int}}, \qquad (2.1) $$

where $\hat{\mathcal{H}}_S$ and $\hat{\mathcal{H}}_E$ respectively denote the Hamiltonians for the system and the environment in the absence of interactions, and $g\hat{\mathcal{H}}_{\text{int}}$ is the interaction term, controlled by the coupling constant $g$. The system alone is described by the reduced density matrix, which is obtained from the full density matrix by tracing over the environmental degrees of freedom,

$$ \hat{\rho}_{\text{red}} = \text{Tr}_{\text{E}}(\hat{\rho}). \qquad (2.2) $$

An evolution equation for $\hat{\rho}_{\text{red}}$ can be derived with different levels of approximation, corresponding to as many different master equations. In this section, we review the most common master equations, see Ref. [41] for further details (readers already familiar with the master-equation basic tools may want to skip this section and jump to Sect. 3).

### 2.1 An exact master equation: the Nakajima–Zwanzig equation

Our first step is to derive an exact, formal master equation, before applying an approximation scheme. Hereafter we work in the interaction picture, where quantum states evolve with the interaction Hamiltonian $g\hat{\mathcal{H}}_{\text{int}}$ and operators evolve with the free Hamiltonian, i.e. the Hamiltonian in the absence of interactions $\hat{\mathcal{H}}_0 \equiv \hat{\mathcal{H}}_S + \hat{\mathcal{H}}_E$. Operators in the interaction picture are denoted with an overall tilde, in order to make the distinction with the Schrödinger and Heisenberg pictures where they carry an overall hat. The link between the Schrödinger and the interaction picture is given by

$$ \begin{aligned} \tilde{\rho}(\eta) &= \hat{U}_0^\dagger(\eta)\hat{\rho}(\eta)\hat{U}_0(\eta) \quad \text{and} \\ \tilde{\mathcal{H}}_{\text{int}}(\eta) &= \hat{U}_0^\dagger(\eta)\hat{\mathcal{H}}_{\text{int}}(\eta)\hat{U}_0(\eta), \end{aligned} \qquad (2.3) $$

where $\eta$ denotes time and where we have introduced the free evolution operator

$$ \begin{aligned} \hat{U}_0(\eta) &= \mathcal{T} \exp \left[ -i \int_{\eta_0}^{\eta} \hat{\mathcal{H}}_0(\eta') d\eta' \right] \\ &= \mathcal{T} \exp \left[ -i \int_{\eta_0}^{\eta} \hat{\mathcal{H}}_S(\eta') d\eta' \right] \\ &\quad \otimes \mathcal{T} \exp \left[ -i \int_{\eta_0}^{\eta} \hat{\mathcal{H}}_E(\eta') d\eta' \right], \end{aligned} \qquad (2.4) $$

with $\mathcal{T}$ indicating time ordering (time arguments increase from right to left). In this work we employ natural units where $\hbar = c = 1$. As mentioned above, in the interaction picture the

total density matrix evolves with the interaction Hamiltonian,

$$ \frac{d\tilde{\rho}}{d\eta} = -ig[\tilde{\mathcal{H}}_{\text{int}}(\eta), \tilde{\rho}(\eta)] \equiv g\mathcal{L}(\eta)\tilde{\rho}(\eta), \qquad (2.5) $$

which defines the Liouville–Von-Neumann super-operator<fnref n="1" /> $\mathcal{L}$.

Let us now introduce the projection super-operator $\mathcal{P}$, defined as

$$ \mathcal{P}\tilde{\rho} = \text{Tr}_{\text{E}}(\tilde{\rho}) \otimes \tilde{\rho}_{\text{E}}, \qquad (2.6) $$

where $\tilde{\rho}_{\text{E}}$ is a fixed reference state in the environment. In practice, it is taken as the state of the environment in the absence of interactions with the system, which is indeed constant in the interaction picture. One can check that $\mathcal{P}$ is a projector, i.e. $\mathcal{P}^2 = \mathcal{P}$, and that $\mathcal{P}\tilde{\rho}$ contains the relevant information to reconstruct the reduced state (2.2) of the system. Upon applying the super-projector $\mathcal{P}$ and its complementary projector $\mathcal{Q} = \text{Id} - \mathcal{P}$ to Eq. (2.5), one obtains

$$ \frac{\partial}{\partial \eta} \mathcal{P} \tilde{\rho}(\eta) = g \mathcal{P} \mathcal{L}(\eta) \tilde{\rho}(\eta), \qquad (2.7) $$

$$ \frac{\partial}{\partial \eta} \mathcal{Q} \tilde{\rho}(\eta) = g \mathcal{Q} \mathcal{L}(\eta) \tilde{\rho}(\eta). \qquad (2.8) $$

Here we have used that since the reference state $\tilde{\rho}_{\text{E}}$ is independent of time, $\mathcal{P}$ and $\mathcal{Q}$ commute with $\partial/\partial\eta$. Inserting the identity $\text{Id} = \mathcal{P} + \mathcal{Q}$ between the Liouville operator and the density matrix, one obtains

$$ \frac{\partial}{\partial \eta} \mathcal{P} \tilde{\rho}(\eta) = g \mathcal{P} \mathcal{L}(\eta) \mathcal{P} \tilde{\rho}(\eta) + g \mathcal{P} \mathcal{L}(\eta) \mathcal{Q} \tilde{\rho}(\eta), \qquad (2.9) $$

$$ \frac{\partial}{\partial \eta} \mathcal{Q} \tilde{\rho}(\eta) = g \mathcal{Q} \mathcal{L}(\eta) \mathcal{P} \tilde{\rho}(\eta) + g \mathcal{Q} \mathcal{L}(\eta) \mathcal{Q} \tilde{\rho}(\eta). \qquad (2.10) $$

A formal solution of Eq. (2.10) is given by

$$ \begin{aligned} \mathcal{Q} \tilde{\rho}(\eta) = & \mathcal{G}_{\mathcal{Q}}(\eta, \eta_0) \mathcal{Q} \tilde{\rho}(\eta_0) \\ & + g \int_{\eta_0}^{\eta} d\eta' \mathcal{G}_{\mathcal{Q}}(\eta, \eta') \mathcal{Q} \mathcal{L}(\eta') \mathcal{P} \tilde{\rho}(\eta'), \end{aligned} \qquad (2.11) $$

where $\eta_0$ is some initial time and $\mathcal{G}_{\mathcal{Q}}(\eta, \eta')$ is the propagator defined as

$$ \mathcal{G}_{\mathcal{Q}}(\eta, \eta') \equiv \mathcal{T} \exp \left[ g \int_{\eta'}^{\eta} d\eta'' \mathcal{Q} \mathcal{L}(\eta'') \right]. \qquad (2.12) $$

Plugging Eq. (2.11) into Eq. (2.9), one then obtains a closed equation for the time evolution of the projected density matrix $\mathcal{P}\tilde{\rho}$, namely

<footnote n="1">¹ In this work, following Ref. [41], “super-operator” denotes an operation which maps positive operators to positive operators.</footnote>
$$ \frac{\partial}{\partial \eta} \mathcal{P} \tilde{\rho}(\eta) = g \mathcal{P} \mathcal{L}(\eta) \mathcal{G}_Q(\eta, \eta_0) \mathcal{Q} \tilde{\rho}(\eta_0) + g \mathcal{P} \mathcal{L}(\eta) \mathcal{P} \tilde{\rho}(\eta) \\ + g^2 \int_{\eta_0}^{\eta} d\eta' \mathcal{P} \mathcal{L}(\eta) \mathcal{G}_Q(\eta, \eta') \mathcal{Q} \mathcal{L}(\eta') \mathcal{P} \tilde{\rho}(\eta'). \qquad (2.13) $$

This is the Nakajima–Zwanzig equation. Although formal, it provides an exact master equation for the reduced state of the system. It can be further simplified by assuming that the initial state does not contain correlations between the system and the environment, i.e. $\tilde{\rho}(\eta_0) = \text{Tr}_E(\tilde{\rho}) \otimes \text{Tr}_S(\tilde{\rho}) = \text{Tr}_E(\tilde{\rho}) \otimes \tilde{\rho}_E$, hence $\mathcal{Q} \tilde{\rho}(\eta_0) = 0$. Moreover, without loss of generality one can assume that the expectation value of the interaction Hamiltonian vanishes in the reference state, i.e. $\text{Tr}_E(\tilde{\mathcal{H}}_{\text{int}}\tilde{\rho}_E) = 0$ [if this is not satisfied, one simply redefines $\tilde{\mathcal{H}}_S$ by adding $g\text{Tr}_E(\tilde{\mathcal{H}}_{\text{int}}\tilde{\rho}_E) \otimes \text{Id}_E$ to it]. This leads to $\mathcal{P} \mathcal{L}(\eta) \mathcal{P} = 0$, so the Nakajima–Zwanzig equation reduces to

$$ \frac{\partial}{\partial \eta} \mathcal{P} \tilde{\rho}(\eta) = g^2 \int_{\eta_0}^{\eta} d\eta' \mathcal{K}(\eta, \eta') \mathcal{P} \tilde{\rho}(\eta'), \qquad (2.15) $$

where we have introduced the memory kernel $\mathcal{K}(\eta, \eta')$ defined as

$$ \mathcal{K}(\eta, \eta') = \mathcal{P} \mathcal{L}(\eta) \mathcal{G}_Q(\eta, \eta') \mathcal{Q} \mathcal{L}(\eta') \mathcal{P}. \qquad (2.16) $$

In this form, the master equation is as difficult to solve as the Liouville equation (2.5) of the full setup. However, it allows efficient approximation schemes to be designed, as we shall now see. The first approximation relies on the assumption of weak coupling between the system and the environment and is discussed in Sect. 2.2, the second approximation concerns properties of the environment itself and is developed in Sect. 2.3.

## 2.2 Born approximation: the time-convolutionless cumulant expansion

An effective description of the system alone is in general possible only when it weakly couples to its environment. This naturally provides a small parameter, namely the interaction strength, in which to perform an expansion. This is the so-called Born approximation, which also addresses one

of the difficulties inherent to the Nakajima–Zwanzig equation (2.15), namely the fact that it is non-local in time, i.e. the time derivative of $\mathcal{P} \tilde{\rho}(\eta)$ depends on its past history $\mathcal{P} \tilde{\rho}(\eta')$ for $\eta' < \eta$. The Time-ConvolutionLess projection operator method (TCL in the following) consists in expanding the dynamics of the system in powers of the coupling constant $g$, rendering the equation local in time (while preserving its non-Markovian nature<fnref n="3" />). One thus obtains an equation of the form

$$ \frac{\partial}{\partial \eta} \mathcal{P} \tilde{\rho}(\eta) = \sum_{n=2}^{\infty} g^n \mathcal{K}_n(\eta) \mathcal{P} \tilde{\rho}(\eta), \qquad (2.17) $$

where the $\mathcal{K}_n$ operators are called the TCL$_n$ operators and can be computed iteratively. This can be done by expanding Eq. (2.12) in $g$, and by using Eq. (2.17) to express $\mathcal{P} \tilde{\rho}(\eta')$ in terms of $\mathcal{P} \tilde{\rho}(\eta)$ in the right-hand side of Eq. (2.15), at the required order. For instance, at leading order in $g$, $\mathcal{G}_Q(\eta, \eta') = \text{Id}$, see Eq. (2.12), so Eq. (2.16) leads to $\mathcal{K}(\eta, \eta') = \mathcal{P} \mathcal{L}(\eta) \mathcal{Q} \mathcal{L}(\eta') \mathcal{P} = \mathcal{P} \mathcal{L}(\eta) \mathcal{L}(\eta') \mathcal{P}$, where we have used that $\mathcal{Q} = 1 - \mathcal{P}$ and that $\mathcal{P} \mathcal{L} \mathcal{P} = 0$, see footnote 2. At that order, Eq. (2.15) also indicates that $\mathcal{P} \tilde{\rho}$ is constant hence

$$ \mathcal{K}_2(\eta) = \int_{\eta_0}^{\eta} d\eta' \mathcal{P} \mathcal{L}(\eta) \mathcal{L}(\eta') \mathcal{P}, \qquad (2.18) $$

and truncating Eq. (2.17) at order $n = 2$ leads to the TCL$_2$ master equation

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -g^2 \int_{\eta_0}^{\eta} d\eta' \text{Tr}_E [\tilde{\mathcal{H}}_{\text{int}}(\eta), [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}_{\text{red}}(\eta) \otimes \tilde{\rho}_E]]. \qquad (2.19) $$

This expansion can be carried on. At order $n = 3$, one needs to expand the memory kernel $\mathcal{K}(\eta, \eta')$ at order $g$ and keep $\mathcal{P} \tilde{\rho}(\eta') \simeq \mathcal{P} \tilde{\rho}(\eta)$ in the right-hand side of Eq. (2.15), given that $\mathcal{P} \tilde{\rho}(\eta') - \mathcal{P} \tilde{\rho}(\eta) = \mathcal{O}(g^2)$ as shown above. One obtains $\mathcal{K}_3(\eta) = \int_{\eta_0}^{\eta} d\eta' \int_{\eta'}^{\eta} d\eta'' \mathcal{P} \mathcal{L}(\eta) \mathcal{Q} \mathcal{L}(\eta'') \mathcal{Q} \mathcal{L}(\eta') \mathcal{P}$, so

$$ \mathcal{K}_3(\eta) = \int_{\eta_0}^{\eta} d\eta' \int_{\eta'}^{\eta} d\eta'' \mathcal{P} \mathcal{L}(\eta) \mathcal{L}(\eta'') \mathcal{L}(\eta') \mathcal{P} \qquad (2.20) $$

where we have used again that $\mathcal{Q} = 1 - \mathcal{P}$ and that $\mathcal{P} \mathcal{L} \mathcal{P} = 0$. Note that, if the odd moments of the interaction Hamiltonian vanish in the environment (as will be the case for the model studied in the rest of this work),

$$ \begin{aligned} \mathcal{P} \mathcal{L} \mathcal{P} \tilde{\rho} &= -i \mathcal{P} [\tilde{\mathcal{H}}_{\text{int}}, \mathcal{P} \tilde{\rho}] \\ &= -i \mathcal{P} [\tilde{\mathcal{H}}_{\text{int}}, \text{Tr}_E(\tilde{\rho}) \otimes \tilde{\rho}_E] \\ &= -i [\text{Tr}_E(\tilde{\mathcal{H}}_{\text{int}}\tilde{\rho}_E), \text{Tr}_E(\tilde{\rho})] \otimes \tilde{\rho}_E \\ &= 0. \end{aligned} \qquad (2.14) $$

<footnote n="2">² This can be shown by computing</footnote>

<footnote n="3">³ In this work, following Ref. [41], the dynamical map $\tilde{\rho}(\eta) \rightarrow \tilde{\rho}(\eta') = \mathcal{M}_{\eta \rightarrow \eta'} \tilde{\rho}(\eta)$ is said to be Markovian if its generators form a semi-group, i.e. $\mathcal{M}_{\eta \rightarrow \eta'} = \mathcal{M}_{\eta'' \rightarrow \eta'} \mathcal{M}_{\eta \rightarrow \eta''}$. Note that a Markovian master equation is necessarily local in time, but the reverse is not necessarily true.</footnote>
i.e. Tr_E[ ℋ_{int}(η₁) ... ℋ_{int}(η₂p+1) ̃ρ_E ] = 0, a similar calculation as the one performed in footnote 2 for p = 1 then shows that P L(η₁) ... L(ηp+1) P = 0. This implies that K₃ vanishes, as well as all odd TCLₙ generators.

In that case, the leading correction comes from TCL₄, which receives two contributions. The first one comes from the term of order g² in the memory kernel K(η, η') while keeping P̃ρ(η') ≃ P̃ρ(η) in the right-hand side of Eq. (2.15). The second contribution comes from keeping the memory kernel at leading order but expand P̃ρ(η') at order g². The latter can be formally obtained from the TCL₂ equation, the solution of which reads P̃ρ(η') = P̃ρ(η₀) + g² ∫η' η₀ K₂(η'')P̃ρ(η'')dη''. Together with Eq. (2.18), this leads to

$$
\begin{equation}
\begin{split}
K_4(\eta) = & \int_{\eta_0}^{\eta} d\eta_1 \int_{\eta_0}^{\eta_1} d\eta_2 \int_{\eta_0}^{\eta_2} d\eta_3 \\
& \times [\mathcal{P}\mathcal{L}(\eta)\mathcal{L}(\eta_1)\mathcal{L}(\eta_2)\mathcal{L}(\eta_3)\mathcal{P} \\
& \quad - \mathcal{P}\mathcal{L}(\eta)\mathcal{L}(\eta_1)\mathcal{P}\mathcal{L}(\eta_2)\mathcal{L}(\eta_3)\mathcal{P} \\
& \quad - \mathcal{P}\mathcal{L}(\eta)\mathcal{L}(\eta_2)\mathcal{P}\mathcal{L}(\eta_1)\mathcal{L}(\eta_3)\mathcal{P} \\
& \quad - \mathcal{P}\mathcal{L}(\eta)\mathcal{L}(\eta_3)\mathcal{P}\mathcal{L}(\eta_1)\mathcal{L}(\eta_2)\mathcal{P}].
\end{split}
\tag{2.21}
\end{equation}
$$

This expansion can be carried on to the required level of
accuracy, which allows one to work out Eq. (2.17) when
truncated at the corresponding order TCL<sub>n</sub>. Note that, even
if the TCL<sub>2</sub> order may be sufficient for practical purposes, the
derivation of the fourth-order generator is useful to control
the validity of the cumulant expansion, by evaluating the error
estimate g<sup>2</sup>||K<sub>4</sub>||/||K<sub>2</sub>|| and checking that it is indeed small.

2.3 Markovian approximation: the Lindblad equation

The TCL₂ master equation (2.19) is in general not Markovian
in the sense given in footnote 3, since it involves a convolu-
tion over the past history through the integral over η'. How-
ever, a further approximation can be performed that renders
the dynamics Markovian. This leads to the so-called Gorini–
Kossakowski–Sudarshan–Lindblad equation, in short Lind-
blad equation in what follows. It can be obtained by first
decomposing the interaction Hamiltonian as

$$
\hat{\mathcal{H}}_{\text{int}}(\eta) = \sum_{i} \hat{\mathbf{o}}_{i}^{(S)}(\eta) \otimes \hat{\mathbf{o}}_{i}^{(E)}(\eta), \qquad (2.22)
$$

where $\hat{\mathbf{o}}_i^{(S)}$ and $\hat{\mathbf{o}}_i^{(E)}$ form a basis of operators acting on
the system and the environment respectively. Plugging this
decomposition into Eq. (2.19), the TCL₂ master equation
reads

$$
\begin{equation}
\begin{split}
\frac{d\tilde{\rho}_{\text{red}}}{d\eta} = & -\sum_{i,j} g^2 \int_{\eta_0}^{\eta} d\eta' \{ \Re \left[ \mathcal{K}_{ij}^{>}(\eta, \eta') \right] \\
& \qquad \times \left[ \tilde{\mathbf{o}}_i^{(S)}(\eta), \left[ \tilde{\mathbf{o}}_j^{(S)\dagger}(\eta'), \tilde{\rho}_{\text{red}}(\eta) \right] \right] \\
& + i \Im \left[ \mathcal{K}_{ij}^{>}(\eta, \eta') \right] \\
& \qquad \times \left[ \tilde{\mathbf{o}}_i^{(S)}(\eta), \left\{ \tilde{\mathbf{o}}_j^{(S)\dagger}(\eta'), \tilde{\rho}_{\text{red}}(\eta) \right\} \right] \},
\end{split}
\tag{2.23}
\end{equation}
$$

where {A, B} ≡ AB + BA denotes the anticommutator and
the memory kernel K̂_{ij}(η, η') is defined as

$$
\mathcal{K}_{ij}^{>} (\eta, \eta') = \mathrm{Tr}_E \left[ \hat{\mathbf{o}}_i^{(E)} (\eta) \hat{\mathbf{o}}_j^{(E)\dagger} (\eta') \hat{\rho}_E \right]. \quad (2.24)
$$

This expression is given in the Heisenberg picture. It involves
the two-point correlation functions of the $\hat{\mathbf{o}}_i^{(E)}$ operators in
the environment, and thus depends on the environment prop-
erties.

Typical environments contain a large number of degrees of freedom, hence they behave as reservoirs in which these correlation functions quickly decay with |η − η'|. More precisely, if the relaxation time of the environment is small compared to the typical time scales over which the system evolves, one may coarse-grain the evolution of the system on scales larger than the environment relaxation time. The memory kernel is then sharply peaked, such that the integral over η' only receives contributions close to its upper bound η. In this limit, the past history (η' < η) is not involved in the dynamics anymore, which therefore becomes Markovian.

Formally, if $\mathcal{K}_{ij}^{>} (\eta, \eta') \propto \delta(\eta - \eta')$, in the Schrödinger picture Eq. (2.23) takes the form

$$
\begin{equation}
\begin{split}
\frac{d\hat{\rho}_{\text{red}}}{d\eta} ={}& -i [\hat{\mathcal{H}}_S(\eta), \hat{\rho}_{\text{red}}(\eta)] \\
& + \sum_{i,j} \mathcal{D}_{ij} \left[ \hat{\mathbf{o}}_i^{(S)} \hat{\rho}_{\text{red}}(\eta) \hat{\mathbf{o}}_j^{(S)} \right. \\
& \qquad \left. - \frac{1}{2} \{ \hat{\mathbf{o}}_j^{(S)} \hat{\mathbf{o}}_i^{(S)}, \hat{\rho}_{\text{red}}(\eta) \} \right],
\end{split}
\tag{2.25}
\end{equation}
$$

where the dissipator matrix Dij is a positive semi-definite
matrix. This entails that it can be diagonalised by a unitary
transformation (due to the hermiticity implied by the positive
semi-definiteness), and in this basis Eq. (2.25) becomes<fnref n="4" />

<footnote n="4">⁴ Another approximation known as the rotating-wave approximation is sometimes performed to obtain the Lindblad equation. Since the evolution of the system is coarse-grained over time scales larger than those describing the dynamics of the environment, this approximation consists in removing the quickly oscillating terms appearing in the master equation, for consistency. The implementation of this approach is however challenging in cosmology, where the dynamical background prevents the existence of a natural frequency basis [55].</footnote>
$$ \frac{d\hat{\rho}_{\text{red}}}{d\eta} = -i [\hat{\mathcal{H}}_S(\eta), \hat{\rho}_{\text{red}}(\eta)] + \sum_k \gamma_k \left[ \hat{\mathcal{L}}_k \hat{\rho}_{\text{red}}(\eta) \hat{\mathcal{L}}_k^\dagger - \frac{1}{2} \{ \hat{\mathcal{L}}_k^\dagger \hat{\mathcal{L}}_k, \hat{\rho}_{\text{red}}(\eta) \} \right] \quad (2.26) $$

where $\hat{\mathcal{L}}_k$ are the so-called jump operators and $\gamma_k$ are the positive eigenvalues of the dissipator matrix. This is called a Lindblad equation and is the most generic form of a Markovian dynamical equation that preserves trace, Hermiticity and positivity of the density matrix [54]. This is why Lindblad equations play a key role when studying environmental effects. However, they rely on strong hypotheses regarding the decay rate of the memory kernel in the environment, which may or may not be always satisfied. Indeed, in the cosmological context, fields evolve on a dynamical background, which implies that the environment does not necessarily reach a stationary state in which fluctuations swiftly decay. One of the goals of this article is to check the reliability of the master-equation approach for cosmological systems.

## 2.4 Link with perturbative methods

Later on in this work, we will investigate the extent to which TCL master equations go beyond perturbative effects and enable some non-perturbative resummation. At this stage however, it is important to stress that, when solved perturbatively, they reduce to standard perturbative results. This is because, when deriving the TCL$_n$ equation, no contribution of order lower than $g^n$ has been dropped.

More explicitly, the Liouville–Von-Neumann equation (2.5) can be formally solved as

$$ \tilde{\rho}(\eta) = |\phi\rangle\langle\phi| - ig \int_{-\infty}^{\eta} d\eta' [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}(\eta')], \quad (2.27) $$

where $|\phi\rangle$ denotes the initial state of the combined system-environment setup. By recursively evaluating $\tilde{\rho}$ in the right-hand side with Eq. (2.27) itself, one obtains

$$ \tilde{\rho}(\eta) = \sum_{n=0}^{\infty} (-ig)^n \int_{-\infty}^{\eta} d\eta_1 \int_{-\infty}^{\eta_1} d\eta_2 \cdots \int_{-\infty}^{\eta_{n-1}} d\eta_n \times [\tilde{\mathcal{H}}_{\text{int}}(\eta_1), [\tilde{\mathcal{H}}_{\text{int}}(\eta_2), \cdots [\tilde{\mathcal{H}}_{\text{int}}(\eta_n), |\phi\rangle\langle\phi| \cdots]], \quad (2.28) $$

which displays all contributions to the quantum state order-by-order in $g$. In turn, this allows one to compute corrections to the observables at all orders, as in the in-in formalism.<fnref n="5" />

$$ \langle \tilde{O} \rangle (\eta) = |\phi| \bar{T} \left[ e^{ig \int_{-\infty}^{\eta} d\eta' \tilde{\mathcal{H}}_{\text{int}}(\eta')} \right] \tilde{O}(\eta) \quad (2.29) $$

Let us see how this compares with a perturbative solution of TCL$_n$. For TCL$_2$, since the right-hand side of Eq. (2.19) is proportional to $g^2$, one has $\tilde{\rho}_{\text{red}}(\eta) \otimes \tilde{\rho}_{\text{E}} = \tilde{\rho}_{\text{red}}(\eta_0) \otimes \tilde{\rho}_{\text{E}} + \mathcal{O}(g^2) = |\phi\rangle\langle\phi| + \mathcal{O}(g^2)$, and Eq. (2.19) leads to

$$ \tilde{\rho}_{\text{red}}(\eta) = |\phi\rangle\langle\phi| - g^2 \int_{\eta_0}^{\eta} d\eta' \text{Tr}_E \times [\tilde{\mathcal{H}}_{\text{int}}(\eta), [\tilde{\mathcal{H}}_{\text{int}}(\eta'), |\phi\rangle\langle\phi|] + \mathcal{O}(g^4). \quad (2.31) $$

Assuming that $\text{Tr}_E(\tilde{\mathcal{H}}_{\text{int}}\tilde{\rho}_{\text{E}}) = 0$ as done above Eq. (2.15), this reduces to Eq. (2.28) when traced over the environmental degrees of freedom and truncated at order $g^2$. This shows that solving TCL$_2$ at order $g^2$ is equivalent to Standard Perturbation Theory (SPT hereafter) at that same order. Likewise, one can show that solving TCL$_n$ perturbatively at order $g^n$ is equivalent to SPT$_n$. Therefore, TCL$_n$ contains all terms of order $g^n$, and some terms of order $g^{m>n}$.<fnref n="6" />

This is why TCL is at least as good as SPT, and one of our goals is to determine how much better it is when employed in a cosmological context. In other words, when master equations are used as bona fide dynamical maps (i.e. when they are taken per se and solved without further perturbative expansion), we want to investigate their ability to resum late-time secular effects in situations of cosmological interest [44,46,59].

## 3 Curved-space Caldeira–Leggett model

Let us now apply the master-equation program to two massive test fields $\phi$ and $\chi$ in a Friedmann–Lemaître–Robertson–Walker geometry, described by the metric

$$ ds^2 = a^2(\eta) (-d\eta^2 + d\mathbf{x}^2), \quad (3.1) $$

where $a$ is the scale factor and $\eta$ is conformal time. For convenience we restrict the analysis to a de-Sitter background

$$ \times \mathcal{T} \left[ e^{-ig \int_{-\infty}^{\eta} d\eta'' \tilde{\mathcal{H}}_{\text{int}}(\eta'')} \right] |\phi\rangle, $$

where $\bar{T}$ denotes anti time-ordering. By Taylor expanding the exponential functions, one obtains

$$ \langle \tilde{O} \rangle (\eta) = \sum_{n=0}^{\infty} (ig)^n \int_{-\infty}^{\eta} d\eta_1 \int_{-\infty}^{\eta_1} d\eta_2 \cdots \int_{-\infty}^{\eta_{n-1}} d\eta_n |\phi| \times [\tilde{\mathcal{H}}_{\text{int}}(\eta_n), [\tilde{\mathcal{H}}_{\text{int}}(\eta_{n-1}), \cdots [\tilde{\mathcal{H}}_{\text{int}}(\eta_1), \tilde{O}(\eta)] \cdots] |\phi\rangle. \quad (2.30) $$

Using that $\langle \tilde{O} \rangle (\eta) = \text{Tr}[\tilde{O}(\eta)\tilde{\rho}(\eta)]$, together with $\text{Tr}[\tilde{O}(\eta)[\tilde{\mathcal{H}}_{\text{int}}(\eta_i), |\phi\rangle\langle\phi|]] = -\langle\phi|[\tilde{\mathcal{H}}_{\text{int}}(\eta_i), \tilde{O}(\eta)]|\phi\rangle$, this is indeed consistent with Eq. (2.28).

<footnote n="5">⁵ This can also be shown in the in-in formalism, where the expectation value of an operator $\tilde{O}$ at time $\eta$ reads</footnote>

<footnote n="6">⁶ Let us stress that since the TCL expansion is organised differently from the one of SPT, it does not admit a straightforward diagrammatic representation. In this sense it is more comparable to the Dynamical Renormalisation Group (DRG) resummation [56–58] where diagrams are partially resummed.</footnote>
for which $a(\eta) \equiv -1/(H\eta)$, where $H$ is the constant Hubble parameter and $\eta$ varies between $-\infty$ and 0. We consider the case where the fields are minimally coupled to gravity and where their self-interaction is quadratic, so the action is of the form

$$S = - \int d^4 x \sqrt{-\det g} \left[ \left( \frac{1}{2} g^{\mu\nu} \partial_{\mu} \varphi \partial_{\nu} \varphi + \frac{1}{2} m^2 \varphi^2 \right) + \left( \frac{1}{2} g^{\mu\nu} \partial_{\mu} \chi \partial_{\nu} \chi + \frac{1}{2} M^2 \chi^2 \right) + \lambda^2 \varphi \chi \right]. \quad (3.2)$$

In this expression, $m$ and $M$ are the masses of the two fields and we assume that they satisfy $m < 3H/2 < M$. So $\varphi$ and $\chi$ can be respectively considered as light and heavy, in the cosmological sense. Having in mind possible applications to cosmological perturbations, where the adiabatic degree of freedom is light, in what follows they will respectively play the role of the system and of the environment. The parameter $\lambda$, which also has dimension of a mass, controls their interaction. If those fields were to describe cosmological perturbations, higher-order interaction terms would be parametrically suppressed, and this setting would correspond to the leading order in cosmological perturbation theory. This model, referred to as the curved-space Caldeira-Leggett model [53,60-62], is therefore of physical interest, and as we shall now see it has the advantage to be exactly solvable.

The quantum state of the fields $\varphi$ and $\chi$ was studied in details in Refs. [52,63], where it was shown that each Fourier sector is placed in a four-mode squeezed state. On super-Hubble scales, the dynamical background leads to the creation of pairs of particles with opposite wave-momenta in each field, and the interaction then entangles these particles, leading to correlations between the two fields. Four-mode squeezed states are Gaussian states, and since the action (3.2) is quadratic Gaussianity is indeed preserved throughout the evolution. Such states are fully described by their covariance matrix (i.e. their quantum two-point expectation values). This is why our goal is now to compute the covariance matrix of the system.

### 3.1 Exact description

The action (3.2) being quadratic, different Fourier modes decouple on a homogeneous background, which makes it useful to introduce

$$v_{\varphi}(\eta, \mathbf{k}) \equiv a(\eta) \int_{\mathbb{R}^3} \frac{d^3\mathbf{x}}{(2\pi)^{3/2}} \varphi(\mathbf{x})e^{-ik \cdot \mathbf{x}} \quad \text{and} \\ v_{\chi}(\eta, \mathbf{k}) \equiv a(\eta) \int_{\mathbb{R}^3} \frac{d^3\mathbf{x}}{(2\pi)^{3/2}} \chi(\mathbf{x})e^{-ik \cdot \mathbf{x}}. \quad (3.3)$$

An additional prefactor $a$ is introduced in these expressions for later convenience. The conjugate momenta can be obtained from Eq. (3.2) and read

$$p_{\varphi} = v'_{\varphi} - \frac{a'}{a} v_{\varphi} \quad \text{and} \quad p_{\chi} = v'_{\chi} - \frac{a'}{a} v_{\chi}, \quad (3.4)$$

where hereafter a prime denotes derivation with respect to the conformal time $\eta$. A Legendre transform gives the Hamiltonian

$$H = \int_{\mathbb{R}^{3+}} d^3 \mathbf{k} \, \mathbf{z}^{\dagger} \mathbf{H}(\eta) \mathbf{z}, \quad (3.5)$$

where the phase-space variables have been arranged into the vector $\mathbf{z} \equiv (v_{\varphi}, p_{\varphi}, v_{\chi}, p_{\chi})^{\mathrm{T}}$, and where $\mathbf{H}$ is a four-by-four matrix given by

$$\mathbf{H}(\eta) = \begin{pmatrix} \mathbf{H}^{(\varphi)} & \mathbf{V} \\ \mathbf{V} & \mathbf{H}^{(\chi)} \end{pmatrix}, \quad (3.6)$$

with

$$\begin{align}
\mathbf{H}^{(\varphi)}(\eta) &= \begin{pmatrix} k^2 + m^2 a^2 \frac{a'}{a} & 1 \\ \frac{a'}{a} & 1 \end{pmatrix}, \notag \\
\mathbf{H}^{(\chi)}(\eta) &= \begin{pmatrix} k^2 + M^2 a^2 \frac{a'}{a} & 1 \\ \frac{a'}{a} & 1 \end{pmatrix}, \notag \\
\mathbf{V}(\eta) &= \begin{pmatrix} \lambda^2 a^2 & 0 \\ 0 & 0 \end{pmatrix}. \tag{3.7}
\end{align}$$

Note that, since $\varphi$ and $\chi$ are real fields, one has $z^*(\eta, \mathbf{k}) = z(\eta, -\mathbf{k})$. This explains why, in order to avoid double counting, the integral in Eq. (3.5) is performed over $\mathbb{R}^{3+} \equiv \mathbb{R}^2 \times \mathbb{R}^+$.

Following the canonical quantisation prescription, field variables are promoted to quantum operators. In order to work with hermitian operators, we split the fields into real and imaginary components, that is

$$\hat{\mathbf{z}} = \frac{1}{\sqrt{2}} (\hat{\mathbf{z}}^{\mathrm{R}} + i \hat{\mathbf{z}}^{\mathrm{I}}), \quad (3.8)$$

such that $\hat{\mathbf{z}}^s$ is Hermitian for $s = \mathrm{R}, \mathrm{I}$. These variables are canonical since $[\hat{v}_i^s(\mathbf{k}), \hat{p}_j^{s'}(\mathbf{q})] = i\delta^3(\mathbf{k} - \mathbf{q})\delta_{i,j}\delta_{s,s'}$ where $i, j = \varphi, \chi$. In this basis, the Hamiltonian takes the same form as in Eq. (3.5), i.e.

$$\hat{H} = \frac{1}{2} \sum_{s=\mathrm{R},\mathrm{I}} \int_{\mathbb{R}^{3+}} d^3 \mathbf{k} \, (\hat{\mathbf{z}}^s)^{\mathrm{T}} \mathbf{H}(\eta) \hat{\mathbf{z}}^s. \quad (3.9)$$

Being separable, there is no mode coupling nor interactions between the $\mathrm{R}$ and $\mathrm{I}$ sectors and the state is factorisable in this decomposition. Hence, from now on, we focus on a given wavenumber $\mathbf{k}$ and a given $s$-sector, and to make notations lighter we leave the $\mathbf{k}$ and $s$ dependence implicit.
A further factorisation can be performed under the field-space rotation

$$ \hat{z} = \underbrace{\begin{pmatrix} \cos \theta & 0 & -\sin \theta & 0 \\ 0 & \cos \theta & 0 & -\sin \theta \\ \sin \theta & 0 & \cos \theta & 0 \\ 0 & \sin \theta & 0 & \cos \theta \end{pmatrix}}_{\mathcal{P}} \hat{z}_{\ell-\text{h}}, \quad \text{where} \quad (3.10) $$

$$ \theta = \frac{1}{2} \arctan \left( \frac{2\lambda^2}{m^2 - M^2} \right), $$

where $\ell$ and $\text{h}$ stand for “light” and “heavy” respectively. In this basis the two fields decouple, and their masses are given by

$$ m_{\ell}^{2} = \frac{1}{2} \left[ m^{2} + M^{2} - (M^{2} - m^{2}) \sqrt{1 + \left( \frac{2\lambda^{2}}{M^{2} - m^{2}} \right)^{2}} \right], \quad (3.11) $$

$$ m_{\text{h}}^{2} = \frac{1}{2} \left[ m^{2} + M^{2} + (M^{2} - m^{2}) \sqrt{1 + \left( \frac{2\lambda^{2}}{M^{2} - m^{2}} \right)^{2}} \right]. \quad (3.12) $$

These expressions imply that $m_{\ell}^2 < m^2 < M^2 < m_{\text{h}}^2$ so after the field-space rotation it remains true that $m_{\ell}^2 < 9H^2/4 < m_{\text{h}}^2$, hence the notation.

In this basis, the problem reduces to the dynamics of two uncoupled free fields evolving in a de-Sitter background. In the Heisenberg picture, this can be cast in terms of the mode-function decomposition

$$ \hat{v}_i(\eta) = v_i(\eta)\hat{a}_i + v_i^*(\eta)\hat{a}_i^\dagger \quad (3.13) $$

for $i = \ell, \text{h}$ and where $\hat{a}_i$ and $\hat{a}_i^\dagger$ are the creation and annihilation operators of the uncoupled fields. Heisenberg’s equation yields the classical equation of motion for the mode functions, i.e.

$$ \begin{aligned} v''_{\ell} + \left( k^2 - \frac{v_{\ell}^2 - \frac{1}{4}}{\eta^2} \right) v_{\ell} &= 0 && \text{and} \\ v''_{\text{h}} + \left( k^2 - \frac{v_{\text{h}}^2 - \frac{1}{4}}{\eta^2} \right) v_{\text{h}} &= 0. \end{aligned} \quad (3.14) $$

In these expressions, $v_\ell = \frac{3}{2}\sqrt{1 - \left(\frac{2m_\ell}{3H}\right)^2}$ and $v_\text{h} = \frac{3}{2}\sqrt{1 - \left(\frac{2m_\text{h}}{3H}\right)^2} \equiv i\mu_\text{h}$. By normalising the mode functions to the Bunch-Davies vacuum [64] in the asymptotic, sub-Hubble past, one obtains<fnref n="7" />

$$ \begin{aligned} v_{\ell}(\eta) &= \frac{1}{2} \sqrt{\frac{\pi z}{k}} e^{i \frac{\pi}{2} (v_{\ell} + \frac{1}{2})} H_{v_{\ell}}^{(1)}(z) \quad \text{and} \\ v_{\text{h}}(\eta) &= \frac{1}{2} \sqrt{\frac{\pi z}{k}} e^{-\frac{\pi}{2} \mu_{\text{h}} + i \frac{\pi}{4}} H_{i\mu_{\text{h}}}^{(1)}(z). \end{aligned} \quad (3.15) $$

In these expressions, $z \equiv -k\eta$ and $H_\nu^{(1)}$ is the Hankel function of the first kind and of order $\nu$. The mode functions of the momenta operators can be obtained by using Eq. (3.4), which still applies in the $\ell - \text{h}$ basis, and one finds

$$ p_{\ell}(\eta) = -\frac{1}{2} \sqrt{\frac{k\pi}{z}} e^{i\frac{\pi}{2}(v_{\ell} + \frac{1}{2})} \times \left[ \left(v_{\ell} + \frac{3}{2}\right) H_{v_{\ell}}^{(1)}(z) - zH_{v_{\ell}+1}^{(1)}(z) \right], \quad (3.16) $$

$$ p_{\text{h}}(\eta) = -\frac{1}{2} \sqrt{\frac{k\pi}{z}} e^{-\frac{\pi}{2}\mu_{\text{h}} + i\frac{\pi}{4}} \times \left[ \left(i\mu_{\text{h}} + \frac{3}{2}\right) H_{i\mu_{\text{h}}}^{(1)}(z) - zH_{i\mu_{\text{h}}+1}^{(1)}(z) \right]. \quad (3.17) $$

As mentioned above, the state being Gaussian, it is fully characterised by the covariance matrix

$$ \Sigma(\eta) = \frac{1}{2} \text{Tr} \left[ \{ \hat{z}(\eta), \hat{z}^{\text{T}}(\eta) \} \hat{\rho}_0 \right], \quad (3.18) $$

where $\hat{\rho}_0$ is the Schrödinger state at initial time, $\hat{\rho}_0 = \hat{\rho}(\eta_0)$. In the uncoupled basis, this leads to a block-diagonal covariance matrix of the form

$$ \begin{aligned} \Sigma_{\ell-\text{h}}(\eta) &= \begin{pmatrix} \Sigma_{\ell}(\eta) & \mathbf{0} \\ \mathbf{0} & \Sigma_{\text{h}}(\eta) \end{pmatrix} && \text{where} \\ \Sigma_{i}(\eta) &= \begin{pmatrix} |v_i(\eta)|^2 & \Re e \left[ v_i(\eta) p_i^*(\eta) \right] \\ \Re e \left[ v_i(\eta) p_i^*(\eta) \right] & |p_i(\eta)|^2 \end{pmatrix} \end{aligned} \quad (3.19) $$

for $i = \ell, \text{h}$. In the $\varphi - \chi$ basis, the covariance matrix can be readily obtained by performing the rotation

$$ \Sigma(\eta) = \mathbf{P} \cdot \Sigma_{\ell-\text{h}}(\eta) \cdot \mathbf{P}^{\text{T}} \equiv \begin{pmatrix} \Sigma_{\varphi\varphi}(\eta) & \Sigma_{\varphi\chi}(\eta) \\ \Sigma_{\varphi\chi}(\eta) & \Sigma_{\chi\chi}(\eta) \end{pmatrix}, \quad (3.20) $$

with

$$ \Sigma_{\varphi\varphi}(\eta) = \cos^2(\theta)\Sigma_{\ell}(\eta) + \sin^2(\theta)\Sigma_{\text{h}}(\eta), \quad (3.21) $$

$$ \Sigma_{\varphi\chi}(\eta) = \cos(\theta) \sin(\theta) [\Sigma_{\ell}(\eta) - \Sigma_{\text{h}}(\eta)], \quad (3.22) $$

$$ \Sigma_{\chi\chi}(\eta) = \cos^2(\theta)\Sigma_{\text{h}}(\eta) + \sin^2(\theta)\Sigma_{\ell}(\eta). \quad (3.23) $$

Finally, the reduced state of the system $\varphi$ is obtained by tracing out the $\chi$ field, see Eq. (2.2). It is still a Gaussian

<footnote n="7">⁷ Note that, since all mass parameters (including $\lambda$) are negligible compared to $k/a$ in the asymptotic past, the Bunch-Davies vacuum can be set both in the $\varphi - \chi$ and in the $\ell - \text{h}$ basis [65]. The vacuum state being invariant under rotations (see Ref. [52]), those two prescriptions are identical.</footnote>
state, with covariance matrix given by $\Sigma_{\varphi\varphi}$ [52]. We have thus found an exact solution to the problem at hand (namely compute the reduced state of the system), to which we will now compare effective methods in order to test their robustness.

As explained in Sect. 1, one of the main physical effects driven by the interaction with an environment is decoherence, namely the transition from a pure quantum state into a statistical mixture. The loss of quantum coherence can be measured with the so-called purity parameter $\gamma(\eta) \equiv \text{Tr}(\hat{\rho}_{\text{red}}^2)$ which measures the amount of quantum entanglement between the system and the environment. Pure states correspond to $\gamma = 1$, and mixed states have $\gamma < 1$ (with $\gamma = 0$ corresponding to a maximally mixed state). The amount by which $\chi$ decoheres $\varphi$ is given by

$$ \gamma(\eta) = \frac{1}{4} \det [\Sigma_{\varphi\varphi}(\eta)]^{-1}, \qquad (3.24) $$

the expression being valid for any Gaussian state [52]. In the absence of interactions between the system and the environment, $\det \Sigma_{\varphi\varphi} = 1/4$ so $\gamma = 1$. Otherwise, the system is said to have decohered when $\gamma \ll 1$.

## 3.2 Effective description: the TCL$_2$ master equation

We now turn our attention to the TCL$_2$ master equation (2.23). We remind that it is formulated in the interaction picture, where the interaction Hamiltonian reads $\tilde{H}_{\text{int}}(\eta) = a^2(\eta)\tilde{v}_{\varphi}(\eta)\tilde{v}_{\chi}(\eta)$. The TCL$_2$ master equation (2.19) thus takes the form

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \left\{ \begin{aligned} & \Re \left[ \mathcal{K}^{>}(\eta, \eta') \right] \\ & \times \left[ \tilde{v}_{\varphi}(\eta), \left[ \tilde{v}_{\varphi}(\eta'), \tilde{\rho}_{\text{red}}(\eta) \right] \right] \\ & \times + i \Im \left[ \mathcal{K}^{>}(\eta, \eta') \right] \left[ \tilde{v}_{\varphi}(\eta), \left\{ \tilde{v}_{\varphi}(\eta'), \tilde{\rho}_{\text{red}}(\eta) \right\} \right] \end{aligned} \right\}, \qquad (3.25) $$

where the memory kernel is given by

$$ \mathcal{K}^{>}(\eta, \eta') \equiv \text{Tr}_E [\hat{v}_{\chi}(\eta) \hat{v}_{\chi}(\eta') \hat{\rho}_E] \qquad (3.26) $$

and we recall that $\hat{\rho}_E$ corresponds to the state of the environment in the absence of interactions with the system [a derivation of Eq. (3.25) following microphysical considerations is also presented in Appendix A]. Since $\hat{v}_{\chi}(\eta)\hat{v}_{\chi}(\eta')$ is not hermitian for $\eta \neq \eta'$, the kernel $\mathcal{K}^{>}(\eta, \eta')$ is complex and can be evaluated as follows. In the interaction picture, operators evolve with the free Hamiltonian, so one can use the results obtained in Sect. 3.1 in the uncoupled basis. More precisely, a similar mode-function decomposition as in Eq. (3.13) can be introduced,

$$ \tilde{v}_i(\eta) = v_i(\eta)\hat{a}_i + v_i^*(\eta)\hat{a}_i^\dagger \qquad (3.27) $$

where $i = \varphi, \chi$, and an analogous expression for $\tilde{p}_i(\eta)$. The mode functions are still given by Eqs. (3.15)–(3.17), where $m_\ell$ and $m_h$ are simply replaced with $m$ and $M$. This leads to

$$ \mathcal{K}^{>}(\eta, \eta') = v_{\chi}(\eta)v_{\chi}^{*}(\eta'). \qquad (3.28) $$

### Interpreting the master equation

While the above form (3.25) of the cosmological master equation is compact, it makes the connection with quantum Brownian motion [62,66–69] less apparent. A form that is easier to interpret can be obtained by expressing all operators at the same time. This can be achieved by inverting the mode-function expansion to yield $\hat{a}_\varphi$ and $\hat{a}_\varphi^\dagger$ in terms of $\tilde{v}_\varphi(\eta)$ and $\tilde{p}_\varphi(\eta)$. Inserting those expressions in Eq. (3.27) evaluated at time $\eta'$ leads to

$$ \tilde{v}_{\varphi}(\eta') = -2\Im \left[ p_{\varphi}(\eta) v_{\varphi}^{*}(\eta') \right] \tilde{v}_{\varphi}(\eta) + 2\Im \left[ v_{\varphi}(\eta) v_{\varphi}^{*}(\eta') \right] \tilde{p}_{\varphi}(\eta). \qquad (3.29) $$

Here we have used that $\Im \left[ v_{\varphi}(\eta) p_{\varphi}^{*}(\eta) \right] = -1/2$, which comes from the canonical commutation relation $[\tilde{v}_{\varphi}(\eta), \tilde{p}_{\varphi}(\eta)] = 1$. Plugging Eq. (3.29) into Eq. (3.25), one finds

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -i \left[ \begin{array}{c} \overbrace{\tilde{H}^{(\text{LS})}(\eta)}^{\tilde{H}} \\ -\frac{1}{2} \sum_{i,j} \mathbf{D}_{ij}(\eta) [\tilde{z}_i(\eta), [\tilde{z}_j(\eta), \tilde{\rho}_{\text{red}}(\eta)]] \\ -\frac{i}{2} \Delta_{12}(\eta) \sum_{i,j} \omega_{ij} [\tilde{z}_i(\eta), [\tilde{z}_j(\eta), \tilde{\rho}_{\text{red}}(\eta)]] \end{array} \right] \qquad (3.30) $$

which defines the “Lamb-shift” Hamiltonian $\tilde{H}^{(\text{LS})}$ (see below), where $\omega = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$, and where we have used the canonical commutation relation again. In this expression, $\tilde{z}(\eta) \equiv (\tilde{v}_{\varphi}(\eta), \tilde{p}_{\varphi}(\eta))^{\text{T}}$ and the two-by-two matrices $\mathbf{D}$ and $\mathbf{\Delta}$ are given by

$$ \mathbf{D}_{11}(\eta) = -4\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \times \Im \left[ p_{\varphi}(\eta) v_{\varphi}^*(\eta') \right] \Re \left[ v_{\chi}(\eta) v_{\chi}^*(\eta') \right], \qquad (3.31) $$

$$ \begin{aligned} \mathbf{D}_{12}(\eta) &= \mathbf{D}_{21}(\eta) \\ &= 2\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \\ & \times \Im \left[ v_{\varphi}(\eta) v_{\varphi}^*(\eta') \right] \Re \left[ v_{\chi}(\eta) v_{\chi}^*(\eta') \right], \end{aligned} \qquad (3.32) $$

$$ \mathbf{D}_{22}(\eta) = 0, \qquad (3.33) $$
and

$$ \Delta_{11}(\eta) = -4\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \\ \times \mathfrak{Im} [p_{\varphi}(\eta)v_{\varphi}^*(\eta')] \mathfrak{Im} [v_{\chi}(\eta)v_{\chi}^*(\eta')], \quad (3.34) $$

$$ \begin{aligned} \Delta_{12}(\eta) &= \Delta_{21}(\eta) \\ &= 2\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \\ & \quad \times \mathfrak{Im} [v_{\varphi}(\eta)v_{\varphi}^*(\eta')] \mathfrak{Im} [v_{\chi}(\eta)v_{\chi}^*(\eta')], \quad (3.35) \\ \Delta_{22}(\eta) &= 0. \end{aligned} \quad (3.36) $$

The corresponding equation in the Schrödinger picture can be obtained using the fact that operators are mapped between the two pictures with the free Hamiltonian of the system, see Eq. (2.3), and one finds<fnref n="8" />

$$ \frac{d\hat{\rho}_{\text{red}}}{d\eta} = -i \left[ \hat{H}^{(\varphi)}(\eta) + \hat{H}^{(\text{LS})}(\eta), \hat{\rho}_{\text{red}}(\eta) \right] \\ + \sum_{i,j} \mathcal{D}_{ij}(\eta) \left[ \hat{z}_i \hat{\rho}_{\text{red}}(\eta) \hat{z}_j - \frac{1}{2} \{ \hat{z}_j \hat{z}_i, \hat{\rho}_{\text{red}}(\eta) \} \right]. \quad (3.38) $$

In this expression, the dissipator matrix is defined as

$$ \begin{aligned} \mathcal{D}(\eta) &\equiv \mathcal{D}(\eta) - i \Delta_{12}(\eta) \omega \\ &= \begin{pmatrix} \mathcal{D}_{11}(\eta) & \mathcal{D}_{12}(\eta) - i \Delta_{12}(\eta) \\ \mathcal{D}_{12}(\eta) + i \Delta_{12}(\eta) & 0 \end{pmatrix}. \end{aligned} \quad (3.39) $$

One can see that Eq. (3.38) has the same form as the Lindblad equation (2.25), with the crucial difference that the dissipator matrix $\mathcal{D}(\eta)$ is not positive semi-definite in the present case.<fnref n="9" /> It is also worth stressing that Eq. (3.30) has the same form as the master equation obtained by Hu, Paz and Zhang in their seminal paper [67] and that describes quantum Brownian

motion. The first term in the right-hand side of Eq. (3.30) provides a unitary contribution, which renormalises the energy levels of the system due to the interaction with the environment [41,79,80]. This is why it is often referred to as the Lamb-shift Hamiltonian. In our case, it reads

$$ \begin{aligned} \hat{H}^{(\varphi)}(\eta) + \hat{H}^{(\text{LS})}(\eta) \\ = \frac{1}{2} \left[ \hat{p}_{\varphi} \hat{p}_{\varphi} + (k^2 + m^2 a^2 + \Delta_{11}) \hat{v}_{\varphi} \hat{v}_{\varphi} \right. \\ \left. + \left( \frac{a'}{a} + \Delta_{12} \right) \{ \hat{v}_{\varphi}, \hat{p}_{\varphi} \} \right]. \end{aligned} \quad (3.40) $$

One can thus see that $\Delta_{11}$ renormalises the mass of the field $\varphi$, while $\Delta_{12}$ renormalises the comoving Hubble parameter. Note that, in the context of effective-field theoretic calculations, these contributions are usually re-absorbed in an effective speed of sound $c_S^2$ [27,81,82]. The second and the third terms in Eq. (3.30) are of a different nature, since they capture the non-unitary evolution of the system and thus cannot be described by an effectively local Lagrangian. This is due to dissipation and decoherence, which respectively correspond to the imaginary and the real part of the dissipator matrix in Eq. (3.38).<fnref n="10" />

Finally, in phase space, the TCL₂ master equation takes the form of a Fokker-Planck equation for the reduced Wigner function $W_{\text{red}}$. The latter is defined by the Wigner-Weyl transform of the reduced density matrix [96], and provides a quantum analogue of a phase-space quasi probability distribution. In Appendix B, we derive general results on the phase-space representation of the TCL₂ master equation. In particular, we find that performing the Wigner-Weyl transform of Eq. (3.30) leads to

$$ \frac{dW_{\text{red}}}{d\eta} = \left\{ \tilde{H}^{(\varphi)} + \tilde{H}^{(\text{LS})}, W_{\text{red}} \right\} + \Delta_{12} \sum_{i} \frac{\partial}{\partial z_{i}} (z_{i} W_{\text{red}}) \\ - \frac{1}{2} \sum_{i,j} [\omega \mathcal{D} \omega]_{ij} \frac{\partial^2 W_{\text{red}}}{\partial z_{i} \partial z_{j}}, \quad (3.41) $$

where brackets correspond to Poisson brackets (not to be confused with the anti-commutator). Only the term involving $\tilde{H}^{(\varphi)} + \tilde{H}^{(\text{LS})}$ is unitary, as mentioned above. The second term, proportional to $\Delta_{12}$, is dissipative: it is a drift (or friction) term that accounts for the energy transfer from the system into the environment [22]. Finally, the term proportional

$$ \begin{aligned} \mathcal{D}_{ij} [\hat{z}_i, [\hat{z}_j, \hat{\rho}_{\text{red}}]] &= \mathcal{D}_{ij} (-2\hat{z}_i \hat{\rho}_{\text{red}} \hat{z}_j + \{\hat{z}_j \hat{z}_i, \hat{\rho}_{\text{red}}\}), \quad (3.37) \\ \omega_{ij} [\hat{z}_i, \{\hat{z}_j, \hat{\rho}_{\text{red}}\}] &= \omega_{ij} (2\hat{z}_i \hat{\rho}_{\text{red}} \hat{z}_j - \{\hat{z}_j \hat{z}_i, \hat{\rho}_{\text{red}}\}). \end{aligned} $$

<footnote n="8">⁸ Here we use that since $\mathcal{D}$ is symmetric, $\omega$ is anti-symmetric given the canonical commutation relations between phase-space variables $[\hat{z}_i, \hat{z}_j] = w_{ij}$, one has</footnote>

<footnote n="9">⁹ If the dynamical map generated by Eq. (3.38) were Markovian in the sense introduced in footnote 3, i.e. if it described a semi-group evolution, then according to Lindblad theorem [54] the fact that its dissipator is not semi-definite positive would imply that it is not Completely Positive and Trace Preserving (CPTP). However, Eq. (3.38) belongs to the class of so-called “Gaussian master equations”, which were shown to be CPTP in Refs. [70,71] (and thus map a quantum state to another proper quantum state). The contrapositive of Lindblad’s theorem thus imposes that our master equation is non-Markovian [32,41,72–78].</footnote>

<footnote n="10">¹⁰ The fact that the real and the imaginary part of the memory kernel lead to distinct physical effects is also encountered in the influence-functional approach [45,66,83–92], of which the master equation is the dynamical generator [93,94]. Indeed, in the influence functional description, $\mathfrak{Im}[\mathcal{K}^>(\eta, \eta')]$ is related to the retarded and advanced Green’s function of the environment and can be interpreted as a dissipation kernel, while $\mathfrak{Re}[\mathcal{K}^>(\eta, \eta')]$ is related to the Keldysh-Green’s function [42,95] and can be interpreted as a noise kernel [51].</footnote>
to $\omega D\omega$ corresponds to diffusion and leads to decoherence. One can show that this equation admits Gaussian solutions, hence the reduced state of the system is still Gaussian in TCL.

### 3.3 Transport equations

As mentioned above, the state being Gaussian, it is fully characterised by its covariance matrix. Since the initial covariance matrix is the same in all approaches (TCL2, exact, SPT) a first strategy to benchmark the cosmological master equation consists in comparing the equation of motion for the covariance of the system, usually refereed to as the transport equations.

#### TCL2 transport equation

In the TCL approach, the transport equations can be obtained by differentiating Eq. (3.18) with respect to time in the Schrödinger picture, and using Eq. (3.38) to evaluate $d\hat{\rho}_{\text{red}}/d\eta$. This gives

$$ \frac{d\Sigma_{\text{TCL}}}{d\eta} = \omega (H^{(\varphi)} + \Delta) \Sigma_{\text{TCL}} - \Sigma_{\text{TCL}} (H^{(\varphi)} + \Delta) \omega \\ - \omega D\omega - 2\Delta_{12}\Sigma_{\text{TCL}}, \qquad (3.42) $$

where $D$ and $\Delta$ were introduced in Eqs. (3.31)–(3.36). The first two terms correspond to the unitary evolution, which as stressed above receives an additional contribution from the Lamb-shift Hamiltonian. The last two terms respectively correspond to the diffusion (a source term proportional to $D$) and the dissipation (a damping term proportional to $\Delta_{12}$).

#### Exact transport equation

In the exact approach presented in Sect. 3.1, the transport equations for the full system-plus-environment setup can be obtained by differentiating Eq. (3.18) with respect to time in the Heisenberg picture, and using the Heisenberg equations to evaluate $d\hat{z}/d\eta$. The Hamiltonian (3.9) being quadratic, one finds

$$ \frac{d\Sigma}{d\eta} = \Omega H \Sigma - \Sigma H \Omega, \qquad (3.43) $$

where $H$ was defined in Eq. (3.6) and $\Omega$ is a four-by-four block-diagonal matrix where each 2×2 block on the diagonal is $\omega$.

Using blockwise multiplication we can split the above into a set of coupled differential equations for the covariance of the system ($\Sigma_{\varphi\varphi}$), of the environment ($\Sigma_{\chi\chi}$), and for their cross-covariance ($\Sigma_{\varphi\chi}$). Using Eq. (3.7), it reads

$$ \frac{d\Sigma_{\varphi\varphi}}{d\eta} = \omega H^{(\varphi)} \Sigma_{\varphi\varphi} - \Sigma_{\varphi\varphi} H^{(\varphi)} \omega \\ + \omega V \Sigma_{\varphi\chi}^{\text{T}} - \Sigma_{\varphi\chi} V \omega, \qquad (3.44) $$

$$ \frac{d\Sigma_{\chi\chi}}{d\eta} = \omega H^{(\chi)} \Sigma_{\chi\chi} - \Sigma_{\chi\chi} H^{(\chi)} \omega \\ + \omega V^{\text{T}} \Sigma_{\varphi\chi} - \Sigma_{\varphi\chi}^{\text{T}} V \omega, \qquad (3.45) $$

$$ \frac{d\Sigma_{\varphi\chi}}{d\eta} = \omega H^{(\varphi)} \Sigma_{\varphi\chi} - \Sigma_{\varphi\chi} H^{(\varphi)} \omega \\ + \omega V \Sigma_{\chi\chi} - \Sigma_{\varphi\varphi} V \omega. \qquad (3.46) $$

Note that these transport equations can also be obtained in the phase-space representation (i.e. using Wigner functions), as explained in Appendix B. In the present case, a first integral of the above system can be easily constructed, since we know that, in spite of having three covariance matrices ($\Sigma_{\varphi\varphi}$, $\Sigma_{\chi\chi}$ and $\Sigma_{\varphi\chi}$), only two combinations are independent (namely $\Sigma_{\ell}$ and $\Sigma_{h}$). More precisely, from Eq. (3.21) one can show that $\Sigma_{\varphi\chi} = \Sigma_{\varphi\chi}^{\text{T}} = \tan(2\theta)(\Sigma_{\varphi\varphi} - \Sigma_{\chi\chi})/2$. Focusing on the dynamics of the reduced system, Eq. (3.44) can thus be written as

$$ \frac{d\Sigma_{\varphi\varphi}}{d\eta} = \omega (H^{(\varphi)} + \Delta_{\text{ex}}) \Sigma_{\varphi\varphi} \\ - \Sigma_{\varphi\varphi} (H^{(\varphi)} + \Delta_{\text{ex}}) \omega - \omega D_{\text{ex}} \omega, \qquad (3.47) $$

where

$$ \Delta_{\text{ex}} \equiv -\frac{\lambda^2}{M^2 - m^2} V \quad \text{and} \\ D_{\text{ex}} \equiv -\frac{\lambda^2}{M^2 - m^2} (\omega \Sigma_{\chi\chi} V - V \Sigma_{\chi\chi} \omega). \qquad (3.48) $$

The reason why we write the exact transport equation in this form is to allow for an easy comparison with its TCL2 counterpart (3.42). This suggests to interpret $\Delta_{\text{ex}}$ as a Lamb-shift contribution in the exact approach, and $D_{\text{ex}}$ as a diffusion matrix. From Eq. (3.48), the only non-vanishing entries of those matrices are given by

$$ \Delta_{\text{ex},11} = -\frac{\lambda^4 a^2}{M^2 - m^2}, \qquad (3.49) $$

$$ D_{\text{ex},11} = -2 \frac{\lambda^4 a^2}{M^2 - m^2} \Sigma_{\chi\chi,12}, $$

$$ D_{\text{ex},12} = \frac{\lambda^4 a^2}{M^2 - m^2} \Sigma_{\chi\chi,11}. \qquad (3.50) $$

Note that, in the asymptotic past, when $a \to 0$, the above coefficients vanish, which confirms that the two fields become effectively uncoupled and that Bunch-Davies initial conditions can be safely set, see footnote 7.

#### SPT transport equation

In the perturbative approach introduced in Sect. 2.4, at leading order, the transport equation is simply given by the exact
transport equation, Eq. (3.47), where the right-hand side is truncated at order $\lambda^4$:

$$ \frac{d\Sigma_{\text{SPT}}}{d\eta} = \omega H^{(\varphi)} \Sigma_{\text{SPT}} - \Sigma_{\text{SPT}} H^{(\varphi)} \omega + \omega \Delta_{\text{ex}} \Sigma_{\varphi\varphi}^{\text{free}} - \Sigma_{\varphi\varphi}^{\text{free}} \Delta_{\text{ex}} \omega - \omega D_{\text{SPT}} \omega. \quad (3.51) $$

Here, $\Sigma_{\varphi\varphi}^{\text{free}}$ corresponds to $\Sigma_{\varphi\varphi}$ evaluated in the free theory and is given by the second part of Eq. (3.19) with the mode functions $v_\varphi$ and $p_\varphi$. Similarly, $D_{\text{SPT}}$ is given by Eq. (3.50) where $\Sigma_{\chi\chi}$ is replaced with $\Sigma_{\chi\chi}^{\text{free}}$, which is given by the second part of Eq. (3.19) with the mode functions $v_\chi$ and $p_\chi$. Note that $\Delta_{\text{ex}}$ does not need to be expanded since it is already of order $\lambda^4$, see Eq. (3.49).

Even though the covariance matrix in SPT can be obtained by integrating the above transport equation, in the present situation an exact solution to the full theory is known, so it can also be obtained by expanding Eq. (3.21) in $\lambda$. Here, not only $\theta^2 = \lambda^4/(m^2 - M^2)^2 + O(\lambda^8)$ needs to be expanded, see Eq. (3.10), but also $m_\ell^2 = m^2 - \lambda^4/(M^2 - m^2) + O(\lambda^8)$ and $m_h^2 = M^2 + \lambda^4/(M^2 - m^2) + O(\lambda^8)$ in $\Sigma_h$ and $\Sigma_\ell$, see Eqs. (3.11) and (3.12). On the numerical results presented below, we have checked that these two approaches coincide.

## 3.4 Spurious terms

The TCL$_2$ coefficients are expressed as integrals between $\eta_0$ and $\eta$, see Eqs. (3.31)–(3.36), where $\eta_0 \to -\infty$ if Bunch-Davies initial conditions are chosen. Formally, they can be written as

$$ D_{11} = F_{D_{11}}(\eta, \eta) - F_{D_{11}}(\eta, \eta_0), \quad (3.52) $$

where $F_{D_{11}}(\eta, \cdot)$ is the primitive of the integrand appearing in Eq. (3.31), which itself depends on $\eta$, and with similar notations for the other TCL$_2$ coefficients. The $F$ functions are derived explicitly in Appendix C, where it is shown that the integrals (3.31)–(3.36) can be performed analytically and involve products of four Hankel functions. The second term in Eq. (3.52), the one of the form $F(\eta, \eta_0)$, features several properties that we now describe and that will lead us to dub it “spurious”.

First, the spurious terms involve the initial time $\eta_0$, which implies that they carry explicit dependence on the initial conditions. If the environment memory kernel (3.26) is sufficiently peaked around $\eta' = \eta$, that is if the integrands in Eqs. (3.31)–(3.36) are much smaller around $\eta' = \eta_0$ than around $\eta' = \eta$, then this contribution should be suppressed compared to the non-spurious one. This is similar to the Lindbladian limit discussed in Sect. 2.3. Whether or not this is the case can be verified explicitly in the super-Hubble regime (i.e. at late time, $-k\eta \ll 1$) where the $F$ functions take

simple forms. The expansion in the limit $-k\eta \ll 1$ is performed in Appendix C.3, where it is shown that the spurious terms dominate for all coefficients. More precisely, $F_{D_{11}}(\eta, \eta) \propto (-k\eta)^{-2}$ while $F_{D_{11}}(\eta, \eta_0) \propto (-k\eta)^{-7/2}$, $F_{D_{12}}(\eta, \eta) \propto (-k\eta)^{-1}$ while $F_{D_{11}}(\eta, \eta_0) \propto (-k\eta)^{-5/2}$, $F_{\Delta_{11}}(\eta, \eta) \propto (-k\eta)^{-2}$ while $F_{\Delta_{11}}(\eta, \eta_0) \propto (-k\eta)^{-7/2}$, and $F_{\Delta_{12}}(\eta, \eta)$ vanishes while $F_{\Delta_{12}}(\eta, \eta_0) \propto (-k\eta)^{-5/2}$. Let us stress that the late-time domination of the spurious terms is strongly related to having a dynamical background. This is the first indication we encounter that applying the master-equation program to cosmology may not be as straightforward as in other situations.

Second, in Appendix C.1, we notice that, using various identities satisfied by the Hankel functions, the expressions for the non-spurious contributions can be vastly simplified. More precisely, after a lengthy though straightforward calculation we find that

$$ \begin{aligned} F_{D_{11}}(\eta, \eta) &= D_{\text{SPT},11}(\eta), \\ F_{D_{12}}(\eta, \eta) &= D_{\text{SPT},12}(\eta), \\ F_{\Delta_{11}}(\eta, \eta) &= \Delta_{\text{ex},11}(\eta), \\ F_{\Delta_{12}}(\eta, \eta) &= 0. \end{aligned} \quad (3.53) $$

Let us now recall the result obtained in Sect. 2.4, namely the fact that the perturbative version of TCL is strictly equivalent to SPT. This implies that $\Sigma_{\text{TCL}} = \Sigma_{\text{SPT}} + O(\lambda^8)$, where $\Sigma_{\text{SPT}}$ only contains terms of order $\lambda^0$ (namely $\Sigma_{\varphi\varphi}^{\text{free}}$) and $\lambda^4$. As a consequence, the right-hand sides of Eqs. (3.42) and (3.51) coincide at order $\lambda^4$. The terms of order $\lambda^0$ are trivially identical, and for the terms of order $\lambda^4$ one obtains (recalling that both $D$ and $\Delta$ are of order $\lambda^4$)

$$ \begin{aligned} & \omega \Delta \Sigma_{\varphi\varphi}^{\text{free}} - \Sigma_{\varphi\varphi}^{\text{free}} \Delta \omega - \omega D \omega - 2 \Delta_{12} \Sigma_{\text{SPT}} \\ &= \omega \Delta_{\text{ex}} \Sigma_{\varphi\varphi}^{\text{free}} - \Sigma_{\varphi\varphi}^{\text{free}} \Delta_{\text{ex}} \omega - \omega D_{\text{SPT}} \omega. \end{aligned} \quad (3.54) $$

Each term in the left-hand side can be decomposed into a non-spurious part and a spurious part, see Eq. (3.52). An important remark is that, thanks to Eq. (3.53), the non-spurious part exactly coincides with the right-hand side, hence the spurious contributions cancel out. We have therefore proven that the spurious terms are absent from the perturbative limit of TCL and only arise at higher order. This is obviously consistent with the fact that, at leading order, TCL coincides with the exact theory, which is not plagued by any spurious contribution.

Third, we have checked that if one includes the spurious terms when solving the TCL transport equation (3.42), then the result quickly blows up. This is due to the late-time divergences of the spurious contributions mentioned above. On the contrary, as we will see below, if one removes them, then the result is remarkably well-behaved.

To summarise, spurious terms cancel out at leading order in the interaction strength, and at higher order, the fact that they carry an explicit dependence on the initial time, com-
bined with their late-time divergent behaviour, indicates that they cannot be resummed. This leads us to conclude that, for the simple model we have considered here, resummation cannot be efficiently performed with the standard master-equation program.

However, this may be due to the over-simplicity of that particular model, which contains a single degree of freedom in the environment. As we further argue in Sect. 5, if “larger” environment are considered, initial-time dependent terms may be parametrically suppressed. In order to gain insight on such situations, in what follows we analyse the consequences of removing the spurious terms “by hand”.

If spurious contributions are removed, Eq. (3.53) indicates that $D = D_{SPT}$ and that $\Delta = \Delta_{ex}$. The $\Delta$ matrix is perfectly captured by TCL since $\Delta_{ex}$ only contains contributions proportional to $\lambda^4$, see Eq. (3.48). In particular, there is no damping term, i.e. $\Delta_{12} = 0$ (note however that this is due to the specifics of the interaction we consider, which is such that $V_{12} = 0$). The diffusive part, i.e. the one driven by $D$, is however only partly contained in TCL, where $D = D_{SPT}$, whereas $D_{ex}$ contains terms of higher-order in $\lambda$. We therefore expect spurious-free TCL to lie somewhere between SPT and the exact theory, which we now further investigate.

## 4 Non-perturbative resummation

In Sect. 2.4, we have shown that the TCL master equation reduces to standard perturbation theory when solved at leading order in the interaction strength. In Appendix D this equivalence is shown explicitly for the toy model introduced in Sect. 3. However, the TCL master equation can also be treated as a bona fide dynamical map for the quantum state of the system, and solved as it is. In that case, its ability to resum secular effects has been investigated in various contexts [22, 43, 44, 46, 59], and we now want to study how late-time resummation proceeds in the (spurious-free) cosmological Caldeira-Leggett model.

### 4.1 Power spectra

As mentioned above, both in the exact and TCL descriptions, the state of the system remains Gaussian, hence it is fully characterised by its covariance matrix, i.e. by its power spectra. This is why we first compare these setups at the level of their power spectra. If the cosmological Caldeira-Leggett model were to describe cosmological perturbations, note that the configuration-configuration power spectrum would be directly related to cosmological observables, such as the CMB temperature anisotropies.

The power spectra in the exact theory are given by Eq. (3.21), and as explained above, by expanding these formulas at first order in $\lambda^4$ one obtains their SPT counterpart. In

the TCL setup, the power spectra can be obtained by solving the transport equation (3.42). In the model under consideration, there is no damping term, $\Delta_{12} = 0$, but in general it can be absorbed by introducing

$$ \sigma_{\text{TCL}} \equiv e^{\Gamma(\eta, \eta_0)} \Sigma_{\text{TCL}} \quad \text{with} \\ \Gamma(\eta, \eta_0) \equiv 2 \int_{\eta_0}^{\eta} d\eta' \Delta_{12}(\eta'), \qquad (4.1) $$

which is solution of a damping-free transport equation, namely

$$ \frac{d\sigma_{\text{TCL}}}{d\eta} = \omega (H^{(\varphi)} + \Delta) \sigma_{\text{TCL}} - \sigma_{\text{TCL}} (H^{(\varphi)} + \Delta) \omega - e^{\Gamma(\eta, \eta_0)} \omega D\omega. \qquad (4.2) $$

This equation can be seen as a homogeneous part, describing unitary evolution, and a source term, describing diffusion. The homogeneous part is generated by the Hamiltonian $H^{(\varphi)} + H^{(\text{LS})}$, and by denoting $g_{\text{LS}}(\eta, \eta_0)$ the associated Green's matrix, the solution of Eq. (4.2) reads

$$ \sigma_{\text{TCL}}(\eta) = g_{\text{LS}}(\eta, \eta_0) \sigma_{\text{TCL}}(\eta_0) g_{\text{LS}}^{\text{T}}(\eta, \eta_0) \\ - \int_{\eta_0}^{\eta} d\eta' e^{\Gamma(\eta', \eta_0)} g_{\text{LS}}(\eta, \eta') \\ \times [\omega D(\eta')\omega] g_{\text{LS}}^{\text{T}}(\eta, \eta'). \qquad (4.3) $$

Note that $g_{\text{LS}}$ is obtained from the Lamb-shift corrected mode functions

$$ g_{\text{LS}}(\eta, \eta') = 2 \left( \begin{array}{c} \Im m \left[ v_{\text{LS}}(\eta) p_{\text{LS}}^{*}(\eta') \right] - \Im m \left[ v_{\text{LS}}(\eta) v_{\text{LS}}^{*}(\eta') \right] \\ \Im m \left[ p_{\text{LS}}(\eta) p_{\text{LS}}^{*}(\eta') \right] - \Im m \left[ p_{\text{LS}}(\eta) v_{\text{LS}}^{*}(\eta') \right] \end{array} \right), \qquad (4.4) $$

where $v_{\text{LS}}$ is the solution of $v_{\text{LS}}'' + \omega_{\text{LS}}^2 v_{\text{LS}} = 0$ where $\omega_{\text{LS}}^2 = k^2 + m^2 a^2 + \Delta_{11} - \Delta_{12} + \Delta_{12}^2 - 2 \Delta_{12} a'/a$, see Eq. (3.40), with Bunch-Davies initial conditions, and $p_{\text{LS}} = v_{\text{LS}}' - (a'/a) v_{\text{LS}}$ as in Eq. (3.4).$^{11}$ This leads to

$$ \Sigma_{\text{TCL}}(\eta) = e^{-\Gamma(\eta, \eta_0)} g_{\text{LS}}(\eta, \eta_0) \Sigma_{\text{TCL}}(\eta_0) g_{\text{LS}}^{\text{T}}(\eta, \eta_0) \\ - \int_{\eta_0}^{\eta} d\eta' e^{-\Gamma(\eta, \eta')} g_{\text{LS}}(\eta, \eta') \\ \times [\omega D(\eta')\omega] g_{\text{LS}}^{\text{T}}(\eta, \eta'). \qquad (4.5) $$

<footnote>$^{11}$ In the present case, since $\Delta_{12} = 0$ and $\Delta_{11} = \Delta_{\text{ex},11}$, where $\Delta_{\text{ex},11}$ is given in Eq. (3.49), one has $\omega_{\text{LS}}^2 = k^2 + [m^2 - \lambda^4/(M^2 - m^2)]a^2$. This implies that $v_{\text{LS}}$ and $p_{\text{LS}}$ can be expressed in terms of Hankel functions as in Eqs. (3.15) and (3.16), with $v$ replaced by $v_{\text{LS}} = \frac{3}{2} \sqrt{1 - \left(\frac{2m_{\text{LS}}}{3H}\right)^2}$ where $m_{\text{LS}}^2 = m^2 - \lambda^4/(M^2 - m^2)$. This is consistent with effective-field theoretic approaches where the masses of light scalar fields are renormalised by heavy fields with contributions $O(\lambda^4/M^2)$ [97].</footnote>
In practice, this integral is computed numerically from a large negative value of $\eta_0$ (sufficiently large that we check the result does not depend on $\eta_0$).

## Growth rate

First we compare in Fig. 1 the growth rate of the configuration-configuration power spectrum, $d \ln \Sigma_{11}/d \ln a$. The result is given in the free theory (i.e. setting $\lambda = 0$, grey line), in the exact theory (black line), in TCL (blue line) and in SPT (orange line). The difference between these different setups becomes more pronounced at late time, on which the inset zooms in. One can see that TCL provides an excellent approximation, better than SPT, which itself is closer to the exact result than the free theory.

The behaviour of the TCL covariance matrix can be further understood by investigating the super-Hubble (i.e. late time, $-k\eta \ll 1$) limit of the transport equation (3.42). In this regime, an expansion of the coefficients can be found in Appendix C.3. By inserting power-law ansatz for the entries of the covariance matrix, one finds that the diffusion term becomes negligible at large scales, and that

$$ \Sigma_{\text{TCL},11} \propto a^{2v_{\text{LS}}-1}, \quad \Sigma_{\text{TCL},12} \propto a^{2v_{\text{LS}}}, \\ \Sigma_{\text{TCL},22} \propto a^{2v_{\text{LS}}+1}, \qquad (4.6) $$

where $v_{\text{LS}}$ was introduced in footnote 11. The corresponding growth rate, $2v_{\text{LS}} - 1$, is displayed in Fig. 1 with the dotted blue line, and one can check that it asymptotes the TCL result at late time indeed.

In the exact theory, the term involving $\Sigma_\ell$ dominates over the one involving $\Sigma_h$ in Eq. (3.21), so the growth rate is

<figure><img src="image_3.png" /><figcaption>Fig. 1 Growth rate of the configuration-configuration power spectrum as a function of time, labeled with the scale factor ($a_* = k/H$ corresponds to the time of Hubble exit, i.e. when $\eta_* = -1/k$). The result is displayed in the free (grey), exact (black), TCL$_2$ (blue) and SPT (orange) theories. The blue dotted line corresponds to the super-Hubble expansion for TCL, see Eq. (4.6), which leads to the growth rate $d \ln \Sigma_{\text{TCL},11}/d \ln a = 2v_{\text{LS}} - 1$. The parameters are set to $\lambda = H$, $m = H/10$ and $M = \sqrt{10}H$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2 Relative error in the configuration-configuration power spectrum in TCL$_2$ ($|\Sigma_{\text{TCL},11} - \Sigma_{\varphi\varphi,11}|/\Sigma_{\varphi\varphi,11}$, blue lines) and SPT ($|\Sigma_{\text{SPT},11} - \Sigma_{\varphi\varphi,11}|/\Sigma_{\varphi\varphi,11}$, orange lines). The result is displayed as a function of time, labeled with the scale factor, and for $M^2 = 10H^2$ (solid lines) and $M^2 = 100H^2$ (dashed lines). The dotted lines correspond to the super-Hubble formula (4.6), which indeed provide a good fit at late time. The parameters are taken as $m^2 = 10^{-4}H^2$ and $\lambda^2 = H^2$. The grey-shaded area is where the error is larger than 100%</figcaption></figure>

given by $2v_\ell - 1$, where $v_\ell$ is given below Eq. (3.14). It is worth stressing that by expanding $v_\ell$ at leading order in $\lambda^4$, one recovers $v_{\text{LS}}$ [namely $m_\ell^2 = m_{\text{LS}}^2 + O(\lambda^8)$]. As a consequence, TCL correctly reproduces the growth rate at first order in $\lambda^4$.

Although this may seem as a perturbative result, let us stress that the resummed non-perturbative feature lies in Eq. (4.6). Indeed, in SPT, expanding $\Sigma_{\varphi\varphi}$ at leading order in $\lambda^4$ leads to

$$ \Sigma_{\text{SPT},11} \propto a^{2v_{\varphi}-1} \left[ 1 + \frac{\lambda^4}{H^2 v_{\varphi} (M^2 - m^2)} \ln a \right] \qquad (4.7) $$

at late time, where $v_\varphi = \frac{3}{2} \sqrt{1 - \left(\frac{2m}{3H}\right)^2}$. This matches Eq. (4.6) at leading order in $\lambda^4$, but Eq. (4.6) contains all higher-order terms in $\lambda^4$ that allow the logs to be resumed. In particular, Eq. (4.7) implies that at late time, the growth rate in SPT approaches the one of the free theory, $2v_\varphi - 1$, while as stated above the growth rate of TCL incorporates the first correction in $\lambda^4$.

## Relative deviation to the exact result

The performance reached by TCL or SPT is given by the relative deviation of their covariance matrices to the exact result. This is displayed in Fig. 2 for $m^2 = 10^{-4}H^2$ and $\lambda^2 = H^2$, which purposely corresponds to a large coupling. One can check that TCL is always more accurate than SPT, and that the difference in accuracy becomes more pronounced at larger $M$. This can be understood as follows. In the super-Hubble regime, TCL behaves according to Eq. (4.6), which is super-imposed in Fig. 2 and indeed provides a good fit. It
<figure><img src="image_5.png" /><figcaption>Fig. 3 Relative error in all entries of the covariance matrix of the system. The blue lines correspond to the TCL₂ result |ΣTCL,ij − Σφφ,ij|/|Σφφ,ij|, while the orange lines correspond to standard perturbation theory (SPT) |Σ(2)φφ,ij − Σφφ,ij|/|Σφφ,ij|. Different line styles correspond to different entries of the covariance matrix, and the parameters of the model are chosen as m² = 10⁻⁴H², M² = 10²H² and a/a* = e⁵. The grey-shaded area is where the error is larger than 100%. The peaky features correspond to where the exact power spectrum Σφφ,12 vanishes and Σφφ,22 goes through a local minimum</figcaption></figure>

leads to

$$
\begin{align}
\frac{|\Sigma_{\text{TCL},11} - \Sigma_{\varphi\varphi,11}|}{\Sigma_{\varphi\varphi,11}} &\approx a^2(v_{\text{LS}} - v_{\ell}) - 1 \nonumber \\
&= \frac{\lambda^8 \ln(a)}{v_{\varphi} H^2 (M^2 - m^2)^3} + \mathcal{O}(\lambda^{12}). \tag{4.8}
\end{align}
$$

The last result is expanded at leading order in λ (hence in ln a), which provides a good approximation as long as the relative error is much smaller than one, as in Fig. 2. In SPT, Eq. (4.7) gives rise to

$$
\frac{|\Sigma_{\text{SPT},11} - \Sigma_{\varphi\varphi,11}|}{\Sigma_{\varphi\varphi,11}} \approx \frac{\lambda^8 \ln^2(a)}{2(M^2 - m^2)^2 H^4 v_{\varphi}^2} + \mathcal{O}(\lambda^{12}) \quad (4.9)
$$

at late time. There are two main differences between Eqs. (4.8) and (4.9). First, when the environment is heavy, M ≫ H, the relative error in TCL decays as λ⁸/M⁶ while it is suppressed by λ⁸/M⁴ in SPT. This explains why, when going from M² = 10H² to M² = 100H² in Fig. 2, the relative error decreases by a factor 10³ in TCL and by a factor 10² in SPT. This indicates that, although both results become more accurate as the environment is heavier, the gain in accuracy is much stronger for TCL. Second, the relative error in SPT increases as ln²(a) at late time, while it only increases as ln(a) in TCL. This is why in Fig. 2, the difference in accuracy between these two approaches becomes even larger as time proceeds.

Finally, in Fig. 3 we display the relative error for all power spectra (i.e. all entries of the covariance matrix), as a function of the interaction strength λ. When λ is small, the relative

<figure><img src="image_6.png" /><figcaption>Fig. 4 Purity parameter (3.24) in the free (grey), exact (black), TCL₂ (blue), and SPT (orange) approaches as a function of time. The dotted line corresponds to the super-Hubble expansion of the TCL result, see Eq. (4.11). The lower panel shows the error of TCL₂ and SPT relative to the exact result. The parameters are set to m² = 10⁻⁴H² and M² = 10²H² with λ² = 10⁻¹H². The perturbative result rapidly diverges, while TCL₂ accurately predicts the amount of decoherence even in the fully decohered regime</figcaption></figure>

error scales as λ⁸ for both SPT and TCL, in agreement with the fact that both methods match the exact result at order λ⁴ [see Sect. 2.4, see also Eqs. (4.8) and (4.9)]. One can also see that both in TCL and in SPT, the reconstruction of the configuration–configuration power spectrum is better than for the configuration–momentum power spectrum, which is itself better than the momentum–momentum power spectrum. In TCL, all power spectra are accurately computed up to large values of λ. For instance, even when λ/H = 1, the relative error is smaller than 10⁻⁴ for all power spectra. In SPT however, the momentum–momentum power spectrum is already out of control for such values of λ. Indeed, the correlators involving the momentum are given with less precision in SPT, and the perturbative expansion breaks down for the momentum–momentum power spectrum much sooner than for the configuration–configuration power spectrum. This will be of prime importance below, since those correlators play an essential role in the process of decoherence.

4.2 Decoherence

We turn our attention to decoherence that we measure using the purity whose expression for Gaussian states is given by Eq. (3.24). The result is displayed in Fig. 4. As time proceeds, the system entangles with its environment, decoherence occurs (i.e. γ decreases away from 1), and the system becomes maximally mixed soon after Hubble-crossing for the parameters used in the figure.

The lower panel displays the error relative to the exact result. One can see that, when time proceeds, the SPT result
quickly diverges. So perturbation theory is only able to describe quasi pure states, for which $1 - \gamma \ll 1$, and breaks down when decoherence proceeds. The reason for the weak performance of SPT is that the purity parameter is driven by the so-called cosmological decaying mode, which is encoded in the power spectra involving the momentum. Around Fig. 3 we saw that those are precisely the correlators that SPT predicts with the least accuracy. On the contrary, TCL₂ remarkably describes the full decoherence process, and is able to approximate the full quantum state even in the strongly decohered regime. The relative error freezes to a tiny value at large scales (here of the order $10^{-6}$), which is a manifestation of the resummation occurring in TCL.

The simplest way to access the late-time behaviour of the purity is to derive an equation of motion for $\det(\Sigma_{\text{TCL}})$ from the transport equation (3.42), namely

$$ \frac{d \det(\Sigma_{\text{TCL}})}{d\eta} = D_{11} \Sigma_{\text{TCL},11} + 2 D_{12} \Sigma_{\text{TCL},12} - 4 \Delta_{12} \det(\Sigma_{\text{TCL}}). \quad (4.10) $$

All unitary contributions (i.e. those involving $H^{(\varphi)}$ and $\Delta_{11}$) have cancelled out (indeed, only non-unitary contributions can change the purity). This implies that diffusion, controlled by $D$, is crucial in the process of decoherence (since $\Delta_{12} = 0$ in the present case). It contrasts with Sect. 4.1 where we had found that $D$ gives negligible corrections to the power spectra on large scales – those negligible corrections are precisely the ones driving decoherence. The results of Appendix C.3 together with Eq. (4.6) indicate that the two first terms of Eq. (4.10) are of the same order $a^{2v_{\text{LS}}-1}$ at late time. In this limit Eq. (4.10) can be integrated, and one obtains

$$ \det(\Sigma_{\text{TCL}}) \approx \frac{1}{4} + \frac{2^{2v_{\text{LS}}-3}}{\pi} \Gamma^2(v_{\text{LS}}) \left(\frac{\lambda}{H}\right)^4 \left(\frac{H}{M}\right)^3 \left(\frac{a}{a_*}\right)^{2v_{\text{LS}}}, \quad (4.11) $$

where the prefactors in Eq. (4.6) have been set by neglecting diffusion (alternatively, they can be set by asymptotic matching at Hubble crossing and this gives a very similar result) and we have neglected contributions exponentially suppressed by $M/H$ to reach a concise expression (they be easily kept but do not bring any particular insight). The purity $\gamma_{\text{TCL}} = 1/(4 \det \Sigma_{\text{TCL}})$ obtained from this expression is displayed in Fig. 4 with the dotted line. One can check that it provides an excellent approximation to the full TCL result, hence to the exact result too.

The above formula (4.11) also allows us to study under which conditions decoherence occurs for the model at hand. It is non perturbative in $\lambda$ since one should recall that $v_{\text{LS}}$ depends on $\lambda$, see footnote 11, although the rate of decoherence is mostly proportional to $(\lambda/H)^4$. Similarly, although $v_{\text{LS}}$ depends on $M$, decoherence occurs at a rate mostly proportional to $(H/M)^3$, so it is slower for heavier environments.

Finally, it is very efficient on super-Hubble scales, since it scales as $a^{2v_{\text{LS}}-3}$, so roughly as the spatial volume, as often encountered [14, 17]. For instance, for scales of astrophysical interest today that are such that $a/a_* \sim e^{50}$ at the end of inflation, for $M/H = 100$ and $m/H = 10^{-2}$, one finds that decoherence proceeds during inflation as soon as $\lambda/H > 10^{-15}$, a very small value indeed.

## 5 Conclusion

Let us now summarise our main results and open up a few prospects. In this work, we have investigated how the master-equation program can be implemented in cosmology. To this end, we have used a toy model where two scalar fields are linearly coupled and evolve on a de-Sitter background. It has the advantage of being exactly solvable, an “integrable system”, in which the performance of effective methods can be assessed and compared with more traditional, perturbative techniques.

We have derived the second-order Time-ConvolutionLess (TCL) equation in this setup, which is a master equation for the reduced density matrix of the system (here the lighter field), and which features the memory kernel of the environment (here the heavier field). It possesses three contributions: a unitary “Lamb-shift” term (renormalisation of the bare Hamiltonian), a dissipation term (energy exchange with the environment) and a diffusive term (driving the quantum decoherence process). They can all be expressed in terms of integrals ranging from the initial time to the time at which the master equation is written.

Usually, the memory kernel is sufficiently peaked around the coincident configuration that these integrals are dominated by their upper bound, hence they carry negligible dependence on the initial time. This is the case if the relaxation time of the environment around its stationary configuration is small compared to the time scale over which the evolution of the system is tracked. This is the so-called Markovian, or Lindbladian limit. In the present case however, due to the presence of a dynamical background, there is no such thing as a stationary configuration for the environment, which strongly departs from being a thermal bath. In practice we find that these integrals carry a non-negligible dependence on the initial time, through a set of terms that we have dubbed “spurious”.

We have then shown that these spurious terms cancel out when the TCL equation is solved perturbatively in the coupling constant, i.e. they are absent from the perturbative version of the theory. This is consistent with the fact that the perturbative solution to the TCL equation is strictly equivalent to standard perturbation theory (such as the in-in formalism for instance). When solving the TCL equation non-
perturbatively however, they lead to unphysical diverging behaviours, which clearly signals their problematic nature.

However, if one removes them “by hand” (which does not necessarily makes the dynamics Markovian, see footnote 9), one finds that the TCL equation provides an excellent approximation to the full theory: it successfully reproduces all power spectra up to large values of the interaction strength, and it tracks the amount of decoherence very accurately, including at late time when the system is in a strongly mixed state. This is due to an explicit resummation of logarithmic terms (i.e. of powers of $\ln a$, where $a$ is the scale factor of the universe), and in Appendix E we show that this resummation is more efficient than the late-time resummation technique proposed in Ref. [44]. The incorporation of these late-time secular effects makes TCL vastly superior to perturbative methods. Although we have found that it does not require particularly heavy environment, the advantage of TCL compared to perturbative methods is even more pronounced when the mass $M$ of the environmental field is larger than the Hubble scale $H$, since the relative error of the former scales as $(H/M)^6$ while it scales as $(H/M)^4$ for the latter.

To summarise, we have found that the master-equation program can be successfully applied in cosmological backgrounds, provided spurious terms are suppressed.

The presence of the spurious terms may be related to the simplicity of our toy model, where only one field is contained in the environment, which can therefore not be considered as a proper reservoir. If multiple fields were present indeed, all with different masses, thus oscillating at different frequencies, the memory kernel would be suppressed away from the coincident limit through the accumulation of random phases [41,73] (technically, the memory kernel would involve some Fourier transform of the mass distribution of the environmental fields, which may be peaked if that distribution is sufficiently broad). This mechanism was studied e.g. in the context of black-hole physics in Ref. [98]. Another possibility would be to consider non-linear interactions between the two fields, which would imply that one Fourier mode in the system couples to all Fourier modes in the environment, hence making the number of environmental degrees of freedom to which the system couples infinite. One could also consider situations in which non-linearities only arise within the environmental sector,<fnref n="12" /> as in quasi-single field models [28,99]. The same mechanism of random phase addition would presumably occur in those cases, which would also lead to a suppression of the spurious terms. Whether or not that suppression is enough should be the subject of further investigations. Another, maybe more adventurous question, is whether or not one can design an improved master equation,

where the removal of spurious contributions is automatically taken care of. Indeed, our results show that master equations free from spurious terms are extremely powerful at deriving reliable predictions for cosmology, and perform much better than perturbative methods. We plan to address these issues in future works.

**Acknowledgements** It is a pleasure to thank Sudhasattwa Brahma, Cliff Burgess, Denis Comelli, Jaime Calderón-Figueroa, Bei-Lok Hu, Christian Käding, Greg Kaplanek and Amaury Micheli for interesting discussions.

**Data Availability Statement** This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data sharing not applicable to this article as no datasets were generated or analysed during the current study.]

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Funded by SCOAP³. SCOAP³ supports the goals of the International Year of Basic Sciences for Sustainable Development.

## Appendix A: Microphysical derivation of the TCL₂ master equation

In this appendix, we present an alternative derivation of the TCL₂ master equation (3.25) in the curved-space Caldeira-Leggett model, which does not rely on the cumulant expansion of the Nakajima-Zwanzig equation. We start from the Liouville-Von-Neumann equation in the interaction picture (2.5), namely

$$ \frac{d\tilde{\rho}}{d\eta} = -i\lambda^2 [\tilde{\mathcal{H}}_{\text{int}}(\eta), \tilde{\rho}(\eta)]. \quad (A.1) $$

As noted in Eq. (2.27), it can be solved formally as

$$ \tilde{\rho}(\eta) = \tilde{\rho}(\eta_0) - i\lambda^2 \int_{\eta_0}^{\eta} d\eta' [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}(\eta')]. \quad (A.2) $$

Inserting this expression into Eq. (A.1), one obtains

$$ \begin{aligned} \frac{d\tilde{\rho}}{d\eta} = & -i\lambda^2 [\tilde{\mathcal{H}}_{\text{int}}(\eta), \tilde{\rho}(\eta_0)] \\ & - \lambda^4 \int_{\eta_0}^{\eta} d\eta' [\tilde{\mathcal{H}}_{\text{int}}(\eta), [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}(\eta')]] + \mathcal{O}(\lambda^6). \end{aligned} \quad (A.3) $$

<footnote n="12">¹² Let us note that the presence of non-linearities, even if confined to the environmental sector, would leave an imprint on the non-Gaussian statistics of the system [28,30].</footnote>
This procedure could be iterated to obtain higher-order nested commutators, controlled by higher powers of the interaction strength. If the coupling constant $\lambda$ is small (Born approximation), one may stop at order $O(\lambda^4)$ where the first non-unitary effects appear.

Our next task is to turn Eq. (A.3) into an ordinary differential equation that is local in time for the reduced density matrix $\tilde{\rho}_{\text{red}}(\eta)$. By tracing Eq. (A.3) over the environmental degrees of freedom, one finds

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} \approx -i\lambda^2 \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), \tilde{\rho}(\eta_0)] \\ - \lambda^4 \int_{\eta_0}^{\eta} d\eta' \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}(\eta')]]. \qquad \text{(A.4)} $$

In the interaction picture, the deviation of $\tilde{\rho}$ from its initial configuration is necessarily controlled by some positive power $p$ of the interaction strength,

$$ \tilde{\rho}(\eta) = \tilde{\rho}_{\text{red}}(\eta_0) \otimes \tilde{\rho}_{\text{E}}(\eta_0) + \lambda^p \tilde{\rho}_{\text{correl}}(\eta) \qquad \text{(A.5)} $$

where $\text{Tr}_{\text{E}}(\tilde{\rho}_{\text{correl}}) = \text{Tr}_{\text{S}}(\tilde{\rho}_{\text{correl}}) = 0$. Consequently,

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -i\lambda^2 \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), \tilde{\rho}_{\text{red}}(\eta_0) \otimes \tilde{\rho}_{\text{E}}(\eta_0)] \\ - i\lambda^{p+2} \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), \tilde{\rho}_{\text{correl}}(\eta_0)] \\ - \lambda^4 \int_{\eta_0}^{\eta} d\eta' \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), \\ [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}_{\text{red}}(\eta') \otimes \tilde{\rho}_{\text{E}}(\eta')]] \\ - \lambda^{p+4} \int_{\eta_0}^{\eta} d\eta' \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), \\ [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}_{\text{correl}}(\eta')]]. \qquad \text{(A.6)} $$

Which term dominates depends on the value of $p$, which can be determined as follows. Let us first recall that $\tilde{\mathcal{H}}_{\text{int}}(\eta)$ was given above Eq. (3.25) and reads

$$ \tilde{\mathcal{H}}_{\text{int}}(\eta) = a^2(\eta) \tilde{v}_{\varphi}(\eta) \tilde{v}_{\chi}(\eta), \qquad \text{(A.7)} $$

which leads to

$$ \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), \tilde{\rho}_{\text{red}}(\eta_0) \otimes \tilde{\rho}_{\text{E}}(\eta_0)] \\ = a^2(\eta) [\tilde{v}_{\varphi}(\eta), \tilde{\rho}_{\text{red}}(\eta_0)] \text{Tr}_{\text{E}}[\tilde{v}_{\chi}(\eta) \tilde{\rho}_{\text{E}}(\eta_0)]. \qquad \text{(A.8)} $$

Note that

$$ \text{Tr}_{\text{E}}[\tilde{v}_{\chi}(\eta)\tilde{\rho}_{\text{E}}(\eta_0)] = \langle \tilde{v}_{\chi}(\eta - \eta_0) \rangle \qquad \text{(A.9)} $$

which is the mean value of the environment field operator. Such a mean value can always be absorbed in a redefinition of the field, such that the first term in the right-hand side of Eq. (A.6) vanishes. From Eq. (A.5), $d\tilde{\rho}_{\text{red}}/d\eta$ necessarily

contains a term of order $O(\lambda^p)$, so the only possibility is that $p = 4$. As a consequence, the terms of order $O(\lambda^{p+2})$ and $O(\lambda^{p+4})$ can be neglected, and one finds

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -\lambda^4 \int_{\eta_0}^{\eta} d\eta' \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), \\ \times [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}_{\text{red}}(\eta') \otimes \tilde{\rho}_{\text{E}}(\eta')]]. \qquad \text{(A.10)} $$

At leading order in $\lambda$, one can safely replace $\tilde{\rho}_{\text{red}}(\eta')$ by $\tilde{\rho}_{\text{red}}(\eta)$ and $\tilde{\rho}_{\text{E}}(\eta')$ by $\hat{\rho}_{\text{E}}$ in Eq. (A.10). This leads to a manifestly time-local equation, namely

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -\lambda^4 \int_{\eta_0}^{\eta} d\eta' \text{Tr}_{\text{E}}[\tilde{\mathcal{H}}_{\text{int}}(\eta), [\tilde{\mathcal{H}}_{\text{int}}(\eta'), \tilde{\rho}_{\text{red}}(\eta) \otimes \hat{\rho}_{\text{E}}]] \qquad \text{(A.11)} $$

which is consistent with Eq. (2.19). Replacing $\mathcal{H}_{\text{int}}$ by Eq. (A.7) and expanding the commutators yields the result

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \\ \times \left\{ \begin{aligned} & [\tilde{v}_{\varphi}(\eta) \tilde{v}_{\varphi}(\eta') \tilde{\rho}_{\text{red}}(\eta) - \tilde{v}_{\varphi}(\eta') \tilde{\rho}_{\text{red}}(\eta) \tilde{v}_{\varphi}(\eta)] \mathcal{K}^{>}(\eta, \eta') \\ & - [\tilde{v}_{\varphi}(\eta) \tilde{\rho}_{\text{red}}(\eta) \tilde{v}_{\varphi}(\eta') - \tilde{\rho}_{\text{red}}(\eta) \tilde{v}_{\varphi}(\eta') \tilde{v}_{\varphi}(\eta)] \mathcal{K}^{<}(\eta, \eta') \end{aligned} \right\}. \qquad \text{(A.12)} $$

where the memory kernels are defined as in Eq. (3.26), namely

$$ \mathcal{K}^{>}(\eta, \eta') \equiv \text{Tr}_{\text{E}}[\hat{v}_{\chi}(\eta) \hat{v}_{\chi}(\eta') \hat{\rho}_{\text{E}}], \qquad \text{(A.13)} $$

$$ \mathcal{K}^{<}(\eta, \eta') \equiv \text{Tr}_{\text{E}}[\hat{v}_{\chi}(\eta') \hat{v}_{\chi}(\eta) \hat{\rho}_{\text{E}}]. \qquad \text{(A.14)} $$

As in Eq. (3.28), they can be expressed in terms of the mode functions

$$ \mathcal{K}^{>}(\eta, \eta') = v_{\chi}(\eta) v_{\chi}^{*}(\eta') \qquad \text{(A.15)} $$

$$ \mathcal{K}^{<}(\eta, \eta') = v_{\chi}^{*}(\eta) v_{\chi}(\eta') = \mathcal{K}^{>*}(\eta, \eta'), \qquad \text{(A.16)} $$

and the master equation reads

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \\ \times \left\{ \begin{aligned} & [\tilde{v}_{\varphi}(\eta) \tilde{v}_{\varphi}(\eta') \tilde{\rho}_{\text{red}}(\eta) - \tilde{v}_{\varphi}(\eta') \tilde{\rho}_{\text{red}}(\eta) \tilde{v}_{\varphi}(\eta)] \mathcal{K}^{>}(\eta, \eta') \\ & - [\tilde{v}_{\varphi}(\eta) \tilde{\rho}_{\text{red}}(\eta) \tilde{v}_{\varphi}(\eta') - \tilde{\rho}_{\text{red}}(\eta) \tilde{v}_{\varphi}(\eta') \tilde{v}_{\varphi}(\eta)] \mathcal{K}^{>*}(\eta, \eta') \end{aligned} \right\}. \qquad \text{(A.17)} $$

Expanding $\mathcal{K}^{>}$ into its real and imaginary part, one recovers Eq. (3.25).
## Appendix B: Phase-space representation of the TCL₂ master equation

An alternative representation of the quantum state is given in the phase-space by the Wigner function (see Ref. [96] for a brief introduction). For Gaussian states, the Wigner function takes the simple form of a multivariate Gaussian [100], which makes it particularly convenient to work with.

The Wigner function is defined as the inverse Weyl transform of the density matrix. For a generic quantum operator $\hat{O}$, the inverse Wigner–Weyl transform reads

$$W_{\hat{O}}(v_{\varphi}, p_{\varphi}) = 2 \int_{-\infty}^{\infty} dy e^{-2ip_{\varphi}y} \langle v_{\varphi} + y | \hat{O} | v_{\varphi} - y \rangle \quad (B.1)$$

and is a function of the phase-space variables $v_{\varphi}$ and $p_{\varphi}$. The above formula is written in the configuration representation, it can also be written in the momentum representation,

$$W_{\hat{O}}(v_{\varphi}, p_{\varphi}) = 2 \int_{-\infty}^{\infty} dk e^{2ikv_{\varphi}} \langle p_{\varphi} + k | \hat{O} | p_{\varphi} - k \rangle. \quad (B.2)$$

In this way, commutators of quantum operators are mapped to the Poisson brackets of their phase-space representations. Indeed, using the above formulas, one finds

$$W_{[\hat{v}_{\varphi}, \hat{O}]} = i \frac{\partial}{\partial p_{\varphi}} W_O \quad \text{and} \\ W_{[\hat{p}_{\varphi}, \hat{O}]} = -i \frac{\partial}{\partial v_{\varphi}} W_O, \quad (B.3)$$

$$W_{\{\hat{v}_{\varphi}, \hat{O}\}} = 2v_{\varphi} W_O \quad \text{and} \\ W_{\{\hat{p}_{\varphi}, \hat{O}\}} = 2p_{\varphi} W_O. \quad (B.4)$$

This leads to

$$i \omega_{ij} W_{[\hat{z}_j, \hat{O}]} = \frac{\partial W_O}{\partial z_i}, \quad (B.5)$$

$$\frac{1}{2} W_{\{\hat{z}_i, \hat{O}\}} = z_i W_O, \quad (B.6)$$

where we have introduced the phase-space vector $z = (v_{\varphi}, p_{\varphi})^T$.

These relations can be used to compute the inverse Weyl transform of the TCL₂ master equation (3.30). Using that $\omega$ is antisymmetric, one finds

$$\frac{d W_{\text{red}}}{d \eta} = \{ \tilde{H}_0 + \tilde{H}^{(\text{LS})}, W_{\text{red}} \} \\ + \Delta_{12} \sum_i \frac{\partial}{\partial z_i} (z_i W_{\text{red}}) \\ - \frac{1}{2} \sum_{i,j} [\omega D \omega]_{ij} \frac{\partial^2 W_{\text{red}}}{\partial z_i \partial z_j}, \quad (B.7)$$

where $W_{\text{red}} = W_{\hat{\rho}_{\text{red}}}$ is the reduced Wigner function, i.e. the inverse Wigner–Weyl transform of the reduced density matrix $\hat{\rho}_{\text{red}}$. The curly brackets now represent Poisson’s brackets, not to be confused with the anticommutators for quantum operators. This coincides with Eq. (3.41).

The first term in Eq. (B.7) corresponds to the free evolution dressed by the Lamb-shift Hamiltonian. This part of the equation only captures unitary/time-reversible evolution. The second term is a damping term reading as a total derivative and the last term is a diffusion term. These last two terms can be combined into a single second-order differential operator involving the dissipator matrix defined in Eq. (3.39), and they induce a non-unitary evolution.

Let us finally mention that the TCL₂ transport equation can be simply obtained from Eq. (B.7) using the Gaussianity of the state. Indeed, the state being Gaussian, the reduced Wigner function is given by

$$W_{\text{red}} = \sqrt{\frac{1}{4\pi^2 \det \boldsymbol{\Sigma}_{\text{TCL}}} \times \exp \left[ -\frac{1}{2} \sum_{i,j} z_i (\boldsymbol{\Sigma}_{\text{TCL}})_{ij}^{-1} z_j \right]}, \quad (B.8)$$

where $\boldsymbol{\Sigma}_{\text{TCL}}$ is the covariance of the reduced system. Upon inserting Eq. (B.8) into Eq. (B.7), one obtains

$$\frac{d\boldsymbol{\Sigma}_{\text{TCL}}}{d\eta} = \boldsymbol{\omega} (\boldsymbol{H}^{(\varphi)} + \boldsymbol{\Delta}) \boldsymbol{\Sigma}_{\text{TCL}} \\ - \boldsymbol{\Sigma}_{\text{TCL}} (\boldsymbol{H}^{(\varphi)} + \boldsymbol{\Delta}) \boldsymbol{\omega} - \boldsymbol{\omega} \boldsymbol{D} \boldsymbol{\omega} \\ - 2\boldsymbol{\Delta}_{12} \boldsymbol{\Sigma}_{\text{TCL}}, \quad (B.9)$$

which indeed coincides with Eq. (3.42).

## Appendix C: Coefficients of the transport equation for TCL₂

In this appendix, we work out the coefficients of the transport equation for TCL₂, defined in Eqs. (3.31), (3.32), (3.34) and (3.35). They involve the scale factor, which in a de-Sitter universe is given by $a = k/(Hz)$, as well as the mode functions

$$v_{\varphi}(\eta) = \frac{1}{2} \sqrt{\frac{\pi z}{k}} e^{i \frac{\pi}{2} (v_{\varphi} + \frac{1}{2})} H_{v_{\varphi}}^{(1)}(z), \quad (C.1)$$

$$p_{\varphi}(\eta) = -\frac{1}{2} \sqrt{\frac{k\pi}{z}} e^{i\frac{\pi}{2}(v_{\varphi} + \frac{1}{2})} \times \left[ \left(v_{\varphi} + \frac{3}{2}\right) H_{v_{\varphi}}^{(1)}(z) - z H_{v_{\varphi}+1}^{(1)}(z) \right], \quad (C.2)$$

$$v_{\chi}(\eta) = \frac{1}{2} \sqrt{\frac{\pi z}{k}} e^{-\frac{\pi}{2}\mu_{\chi} + i\frac{\pi}{4}} H_{i\mu_{\chi}}^{(1)}(z), \quad (C.3)$$
$$p_{\chi}(\eta) = -\frac{1}{2}\sqrt{\frac{k\pi}{z}}e^{-\frac{\pi}{2}\mu_{\chi}+i\frac{\pi}{4}} \times \left[ \left(i\mu_{\chi} + \frac{3}{2}\right) H_{i\mu_{\chi}}^{(1)}(z) - zH_{i\mu_{\chi}+1}^{(1)}(z) \right]. \quad (C.4)$$

Here, we recall that $z = -k\eta$, $H_v^{(1)}$ is the Hankel function of the first kind and of order $v$ and

$$v_{\varphi} = \frac{3}{2}\sqrt{1 - \left(\frac{2m}{3H}\right)^2} \quad \text{and} \quad \mu_{\chi} = \frac{3}{2}\sqrt{\left(\frac{2M}{3H}\right)^2 - 1}. \qquad (C.5)$$

### C.1 Exact results

In order to perform the integrals involved in Eqs. (3.31)–(3.35), we will make use of the formula

$$\int C_{v_1}(Az)D_{v_2}(Az) \frac{dz}{z} = \frac{C_{v_1}(Az)D_{v_2}(Az)}{v_1 + v_2} + \frac{Az}{v_1^2 - v_2^2} [C_{v_1}(Az)D_{v_2+1}(Az) - C_{v_1+1}(Az)D_{v_2}(Az)], \quad (C.6)$$

see Eq. (10.22.6) of Ref. [101], where $C_{v_1}$ and $D_{v_2}$ are any of the Bessel functions, and $A$ is a fixed arbitrary parameter. Anticipating the computation, let us finally define

$$F_{v,\mu}(z) = \frac{H_v^{(2)}(z)H_{i\mu}^{(1)}(z)}{\nu + i\mu} + \frac{z}{\nu^2 + \mu^2} \times \left[ H_v^{(2)}(z)H_{i\mu+1}^{(1)}(z) - H_{v+1}^{(2)}(z)H_{i\mu}^{(1)}(z) \right], \quad (C.7)$$

$$G_{v,\mu}(z) = \frac{H_v^{(2)}(z)H_{-i\mu}^{(2)}(z)}{\nu - i\mu} + \frac{z}{\nu^2 + \mu^2} \times \left[ H_v^{(2)}(z)H_{-i\mu+1}^{(2)}(z) - H_{v+1}^{(2)}(z)H_{-i\mu}^{(2)}(z) \right], \quad (C.8)$$

in terms of which it will be convenient to express our results.

#### $D_{11}$ coefficient

We start with $D_{11}$ defined in Eq. (3.31), namely

$$D_{11}(\eta) = -4\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \times \mathfrak{Im} \left[ p_{\varphi}(\eta) v_{\varphi}^*(\eta') \right] \Re \left[ v_{\chi}(\eta) v_{\chi}^*(\eta') \right]. \quad (C.9)$$

Expanding the real and the imaginary parts, and replacing $a = k/(Hz)$, it is given by

$$D_{11}(z) = i \frac{k^3 \lambda^4}{z^2 H^4} \int_z^{z_0} \frac{dz'}{(z')^2} \left[ p_{\varphi}(z) v_{\varphi}^*(z') - p_{\varphi}^*(z) v_{\varphi}(z') \right] \times \left[ v_{\chi}(z) v_{\chi}^*(z') + v_{\chi}^*(z) v_{\chi}(z') \right]. \quad (C.10)$$

We thus have four terms,

$$D_{11}(z) = i \frac{k^3 \lambda^4}{z^2 H^4} p_{\varphi}(z) v_{\chi}(z) \int_z^{z_0} \frac{dz'}{(z')^2} v_{\varphi}^*(z') v_{\chi}^*(z') \\ + i \frac{k^3 \lambda^4}{z^2 H^4} p_{\varphi}(z) v_{\chi}^*(z) \int_z^{z_0} \frac{dz'}{(z')^2} v_{\varphi}^*(z') v_{\chi}(z') \\ - i \frac{k^3 \lambda^4}{z^2 H^4} p_{\varphi}^*(z) v_{\chi}(z) \int_z^{z_0} \frac{dz'}{(z')^2} v_{\varphi}(z') v_{\chi}^*(z') \\ - i \frac{k^3 \lambda^4}{z^2 H^4} p_{\varphi}^*(z) v_{\chi}^*(z) \int_z^{z_0} \frac{dz'}{(z')^2} v_{\varphi}(z') v_{\chi}(z'), \quad (C.11)$$

which can be re-organised as

$$D_{11}(z) = 2 \frac{k^3 \lambda^4}{z^2 H^4} \mathfrak{Im} \left[ p_{\varphi}(z) v_{\chi}(z) \int_{z_0}^{z} \frac{dz'}{(z')^2} v_{\varphi}^*(z') v_{\chi}^*(z') \right. \\ \left. + p_{\varphi}(z) v_{\chi}^*(z) \int_{z_0}^{z} \frac{dz'}{(z')^2} v_{\varphi}^*(z') v_{\chi}(z') \right]. \quad (C.12)$$

Therefore, we have two integrals to compute. Making use of Eq. (C.6), they are given by

$$\int_{z_0}^{z} \frac{dz'}{(z')^2} v_{\varphi}^*(z') v_{\chi}^*(z') = -i \frac{\pi}{4k} e^{-\frac{\pi}{2} \mu_{\chi} - i \frac{\pi}{2} v_{\varphi}} \times [G_{v_{\varphi}, \mu_{\chi}}(z) - G_{v_{\varphi}, \mu_{\chi}}(z_0)], \\ \int_{z_0}^{z} \frac{dz'}{(z')^2} v_{\varphi}^*(z') v_{\chi}(z') = \frac{\pi}{4k} e^{-\frac{\pi}{2} \mu_{\chi} - i \frac{\pi}{2} v_{\varphi}} \times [F_{v_{\varphi}, \mu_{\chi}}(z) - F_{v_{\varphi}, \mu_{\chi}}(z_0)]. \quad (C.13)$$

We conclude that $D_{11}$ can be written as

$$D_{11}(z) = F_{D_{11}}(z, z) - F_{D_{11}}(z, z_0), \quad (C.14)$$

where

$$F_{D_{11}}(z_1, z_2) = \frac{\pi k^2 \lambda^4}{2 z^2 H^4} e^{-\frac{\pi}{2} \mu_{\chi}} \times \mathfrak{Im} \left[ -ip_{\varphi}(z_1) v_{\chi}(z_1) G_{v_{\varphi}, \mu_{\chi}}(z_2) e^{-i \frac{\pi}{2} v_{\varphi}} \right. \\ \left. + p_{\varphi}(z_1) v_{\chi}^*(z_1) F_{v_{\varphi}, \mu_{\chi}}(z_2) e^{-i \frac{\pi}{2} v_{\varphi}} \right]. \quad (C.15)$$

It is worth noting that in the case where $z_1 = z_2$, this function can be further simplified by making repeated use of the
Wronskian identity (see Eq. (10.5.5) of Ref. [101]), namely

$$H_{\nu+1}^{(1)}(z)H_{\nu}^{(2)}(z) - H_{\nu}^{(1)}(z)H_{\nu+1}^{(2)}(z) = -\frac{4i}{\pi z}. \quad (C.16)$$

Recalling that $[H_ν^{(1)}(z)]^* = H_ν^{(2)}(z)$ and $[H_ν^{(2)}(z)]^* = H_ν^{*(1)}(z)$, after a tedious but straightforward calculation it leads to

$$F_{D_{11}}(z, z) = -\frac{2}{v_{\varphi}^2 + \mu_{\chi}^2} \left(\frac{k}{z}\right) \times \left(\frac{\lambda}{H}\right)^4 \mathrm{Re} \left[v_{\chi}(z) p_{\chi}^{*}(z)\right]. \quad (C.17)$$

Given that Eq. (C.5) leads to $v_φ^2 + μ_χ^2 = (M^2 - m^2)/H^2$, and since $\Sigma_{\chi\chi,12}^{\text{free}}(z) = \text{Re}[v_χ(z)p_χ^*(z)]$, this can be rewritten as

$$F_{D_{11}}(z, z) = -2 \frac{\lambda^4 a^2}{M^2 - m^2} \Sigma_{\chi\chi,12}^{\text{free}}(z), \quad (C.18)$$

where we have also used that $a = -k/(Hz)$ in a de-Sitter universe. This corresponds to $D_{11}$ in the exact theory when evaluated at leading order in the interaction strength, see Eq. (3.50). In other words, we have shown that

$$F_{D_{11}}(z, z) = D_{SPT,11}(z). \quad (C.19)$$

### $D_{12}$ coefficient

The other coefficients can be computed similarly. For $D_{12}$ defined in Eqs. (3.32), one has

$$\begin{aligned} D_{12}(\eta) &= 2\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \Im [v_{\phi}(\eta) v_{\phi}^*(\eta')] \\ & \quad \times \Re [v_{\chi}(\eta) v_{\chi}^*(\eta')] \\ &= \frac{i}{2} \frac{k^3}{z^2} \frac{\lambda^4}{H^4} \int_{z_0}^{z} \frac{dz'}{(z')^2} [v_{\phi}(z) v_{\phi}^*(z') - v_{\phi}^*(z) v_{\phi}(z')] \\ & \quad \times [v_{\chi}(z) v_{\chi}^*(z') + v_{\chi}^*(z) v_{\chi}(z')]. \end{aligned} \quad (C.20)$$

This leads to

$$D_{12}(z) = F_{D_{12}}(z, z) - F_{D_{12}}(z, z_0), \quad (C.21)$$

where

$$\begin{aligned} F_{D_{12}}(z_1, z_2) = & -\frac{\pi k^2 \lambda^4}{4 z^2 H^4} e^{-\frac{\pi}{2} \mu_\chi} \\ & \times \Im \left[ -i v_\phi(z_1) v_\chi(z_1) G_{v_\phi, \mu_\chi}(z_2) e^{-i \frac{\pi}{2} v_\phi} \right. \\ & \left. \qquad + v_\phi(z_1) v_\chi^*(z_1) F_{v_\phi, \mu_\chi}(z_2) e^{-i \frac{\pi}{2} v_\phi} \right]. \end{aligned} \quad (C.22)$$

As for $F_{D_{11}}$, this expression can be simplified in the coincident configuration $z_1 = z_2$ by repeatedly using the Wronskian identity (C.16), and one finds

$$F_{D_{12}}(z, z) = \left(\frac{k}{z}\right)^2 \left(\frac{\lambda}{H}\right)^4 \frac{|v_{\chi}(z)|^2}{\mu_{\chi}^2 + v_{\phi}^2}. \quad (C.23)$$

Using again that $v_φ^2 + μ_χ^2 = (M^2 - m^2)/H^2$, this can be written as

$$F_{D_{12}}(z, z) = \frac{\lambda^4 a^2}{M^2 - m^2} \Sigma_{\chi\chi,11}^{\text{free}}(z) = D_{SPT,12}(z), \quad (C.24)$$

where we recognise the leading-order contribution in $D_{ex,12}$, see Eq. (3.50).

### $\Delta_{11}$ coefficient

For $\Delta_{11}$ defined in Eqs. (3.34), one has

$$\begin{aligned} \Delta_{11}(\eta) &= -4\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \\ & \quad \times \Im [p_{\phi}(\eta) v_{\phi}^*(\eta')] \Im [v_{\chi}(\eta) v_{\chi}^*(\eta')] \\ &= -\frac{k^3}{z^2} \frac{\lambda^4}{H^4} \int_{z_0}^{z} \frac{dz'}{(z')^2} [p_{\phi}(z) v_{\phi}^*(z') - p_{\phi}^*(z) v_{\phi}(z')] \\ & \quad \times [v_{\chi}(z) v_{\chi}^*(z') - v_{\chi}^*(z) v_{\chi}(z')]. \end{aligned} \quad (C.25)$$

This leads to

$$\Delta_{11}(z) = F_{\Delta_{11}}(z, z) - F_{\Delta_{11}}(z, z_0), \quad (C.26)$$

where

$$\begin{aligned} F_{\Delta_{11}}(z_1, z_2) = & -\frac{\pi k^2 \lambda^4}{2 z^2 H^4} e^{-\frac{\pi}{2} \mu_\chi} \\ & \times \Re \left[ -i p_\phi(z_1) v_\chi(z_1) G_{v_\phi, \mu_\chi}(z_2) e^{-i \frac{\pi}{2} v_\phi} \right. \\ & \left. \qquad - p_\phi(z_1) v_\chi^*(z_1) F_{v_\phi, \mu_\chi}(z_2) e^{-i \frac{\pi}{2} v_\phi} \right]. \end{aligned} \quad (C.27)$$

This expression can be simplified when $z_1 = z_2$ using the Wronskian identity (C.16), if one further uses two additional properties of the Hankel functions. The first one is the recurrence relation (see Eq. (10.6.1) of Ref. [101])

$$H_{\nu-1}^{(2)}(z) + H_{\nu+1}^{(2)}(z) = \frac{2\nu}{z} H_{\nu}^{(2)}, \quad (C.28)$$

and the second one is the inversion formula

$$H_{-\nu}^{(1)}(z) = e^{i\pi\nu} H_{\nu}^{(1)}(z), \quad (C.29)$$

$$H_{-\nu}^{(2)}(z) = e^{-i\pi\nu} H_{\nu}^{(2)}(z).$$
After a tedious but straightforward calculation, this leads to

$$F_{\Delta_{11}}(z, z) = - \left(\frac{k}{z}\right)^2 \left(\frac{\lambda}{H}\right)^4 \frac{1}{\nu_{\phi}^2 + \mu_{\chi}^2} \qquad (\text{C.30}) \\ = \Delta_{\text{ex,11}},$$

where we have recognised $\Delta_{\text{ex,11}}$, see Eq. (3.49), using again that $\nu_{\phi}^2 + \mu_{\chi}^2 = (M^2 - m^2)/H^2$. Note that, here, the agreement between $F_{\Delta_{11}}(z, z)$ and $\Delta_{\text{ex,11}}$ is valid at all orders, given that $\Delta_{\text{ex,11}}$ only contains terms of order $\lambda^4$.

### $\Delta_{12}$ coefficient

Finally, for $\Delta_{12}$ defined in Eqs. (3.35), one has

$$\begin{aligned} \Delta_{12}(\eta) = & 2\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \Im[v_{\phi}(\eta)v_{\phi}^*(\eta')] \\ & \times \Im[v_{\chi}(\eta)v_{\chi}^*(\eta')] \\ = & \frac{1}{2} \frac{k^3}{z^2} \frac{\lambda^4}{H^4} \int_{z_0}^{z} \frac{dz'}{(z')^2} [v_{\phi}(z)v_{\phi}^*(z') - v_{\phi}^*(z)v_{\phi}(z')] \\ & \times [v_{\chi}(z)v_{\chi}^*(z') - v_{\chi}^*(z)v_{\chi}(z')], \end{aligned} \qquad (\text{C.31})$$

and this leads to

$$\Delta_{12}(z) = F_{\Delta_{12}}(z, z) - F_{\Delta_{12}}(z, z_0), \qquad (\text{C.32})$$

where

$$F_{\Delta_{12}}(z_1, z_2) = \frac{\pi k^2}{4 z^2} \frac{\lambda^4}{H^4} e^{-\frac{\pi}{2} \mu_{\chi}} \times \Re \left[ -i v_{\phi}(z_1) v_{\chi}(z_1) G_{v_{\phi}, \mu_{\chi}}(z_2) e^{-i \frac{\pi}{2} v_{\phi}} - v_{\phi}(z_1) v_{\chi}^*(z_1) F_{v_{\phi}, \mu_{\chi}}(z_2) e^{-i \frac{\pi}{2} v_{\phi}} \right]. \qquad (\text{C.33})$$

This expression can be simplified when $z_1 = z_2$ using the Wronskian identity (C.16), and we find that it vanishes,

$$F_{\Delta_{12}}(z, z) = 0. \qquad (\text{C.34})$$

In particular, it implies that $F_{\Delta_{12}}(z, z) = \Delta_{\text{ex,12}}(z)$ and that, as for $F_{\Delta_{11}}(z, z)$, this is valid at all orders in $\lambda^4$ given that both quantities identically vanish.

### C.2 Sub-Hubble limit

In order to gain analytic insight, let us expand the coefficients derived above in the sub-Hubble ($z \gg 1$) and super-Hubble ($z \ll 1$) limits. In the sub-Hubble limit, one can use the asymptotic expansion

$$H_{\nu}^{(1)}(z) = \sqrt{\frac{2}{\pi z}} e^{-iz - i \frac{\pi}{2} \nu - i \frac{\pi}{4}} \sum_{k=0}^{\infty} a_k(\nu) \left(\frac{i}{z}\right)^k, \qquad (\text{C.35})$$

$$H_{\nu}^{(2)}(z) = \sqrt{\frac{2}{\pi z}} e^{i z + i \frac{\pi}{2} \nu + i \frac{\pi}{4}} \sum_{k=0}^{\infty} a_k(\nu) \left(\frac{i}{z}\right)^k, \qquad (\text{C.36})$$

see Eq. (10.17.5) of Ref. [101], with

$$a_k(\nu) = \frac{\left(\frac{1}{2} - \nu\right)_k \left(\frac{1}{2} + \nu\right)_k}{(-2)^k k!} \qquad (\text{C.37})$$

where the parenthesis with lower index indicate the Pochhammer's symbol, i.e. $(x)_k = \Gamma(x + k)/\Gamma(x)$. Inserting these formulas into Eqs. (C.7) and (C.8), one obtains

$$F_{v,\mu}(z) = - \frac{e^{\frac{\pi}{2}(\mu+iv)}}{\pi} \left( \frac{4i}{\nu^2 + \mu^2} + \frac{2}{z} - i \frac{\nu^2 + \mu^2}{2z^2} \right) + \mathcal{O}(z^{-3}), \qquad (\text{C.38})$$

$$G_{v,\mu}(z) = - \frac{e^{\frac{\pi}{2}(\mu+iv)}}{\pi} \frac{e^{-2iz}}{z^2} + \mathcal{O}(z^{-3}). \qquad (\text{C.39})$$

Note that $F_{v,\mu}(z)$ is non vanishing in the sub-Hubble regime. Let us now expand Eqs. (C.15), (C.22), (C.27) and (C.33) in the limit $z_1, z_2 \gg 1$. At leading order, one obtains

$$F_{D_{11}}(z_1, z_2) \approx \frac{k^2 \lambda^4}{2H^4 z_1^3} \left( \frac{2}{\nu_{\phi}^2 + \mu_{\chi}^2} - 1 + \frac{z_1}{z_2} \right), \qquad (\text{C.40})$$

$$F_{D_{12}}(z_1, z_2) \approx \frac{k \lambda^4}{2H^4 \left(\nu_{\phi}^2 + \mu_{\chi}^2\right) z_1^2}, \qquad (\text{C.41})$$

$$F_{\Delta_{11}}(z_1, z_2) \approx - \frac{k^2 \lambda^4}{H^4 \left(\nu_{\phi}^2 + \mu_{\chi}^2\right) z_1^2}, \qquad (\text{C.42})$$

$$F_{\Delta_{12}}(z_1, z_2) \approx \frac{k \lambda^4 (z_1 - z_2)}{4H^4 z_1^3 z_2}, \qquad (\text{C.43})$$

### C.3 Super-Hubble limit

To organise the super-Hubble expansion, we introduce the quantities

$$\alpha_{\nu}(z) \equiv \frac{1 + i \cot \pi \nu}{\Gamma(1 + \nu)} \left(\frac{z}{2}\right)^{\nu - \frac{3}{2}}, \\ \beta_{\nu}(z) \equiv \frac{-i}{\sin \pi \nu} \frac{1}{\Gamma(1 - \nu)} \left(\frac{z}{2}\right)^{\frac{3}{2} - \nu}, \qquad (\text{C.44})$$

$$\gamma_{\mu}(z) \equiv \frac{1 + \coth \pi \mu}{\Gamma(1 + i \mu)} \left(\frac{z}{2}\right)^{i \mu}, \\ \delta_{\mu}(z) \equiv \frac{-1}{\sinh \pi \mu} \frac{1}{\Gamma(1 - i \mu)} \left(\frac{z}{2}\right)^{-i \mu}, \qquad (\text{C.45})$$

together with the function

$$f_{\chi}(z) \equiv \sum_{k=0}^{\infty} \frac{(-1)^k \left(\frac{z}{2}\right)^{2k}}{k! (x+1)_k} \qquad (\text{C.46})$$
$$ = 1 - \frac{\left(\frac{z}{2}\right)^2}{x+1} + \frac{\left(\frac{z}{2}\right)^4}{2(x+1)(x+2)} + \mathcal{O}(z^6) \quad (C.47) $$

such that

$$ H_{\nu}^{(1)}(z) = \alpha_{\nu}(z) f_{\nu}(z) \left(\frac{z}{2}\right)^{\frac{3}{2}} + \beta_{\nu}(z) f_{-\nu}(z) \left(\frac{z}{2}\right)^{-\frac{3}{2}}, \quad (C.48) $$

$$ H_{i\mu}^{(1)}(z) = \gamma_{\mu}(z) f_{i\mu}(z) + \delta_{\mu}(z) f_{-i\mu}(z), \quad (C.49) $$

and

$$ H_{\nu+1}^{(1)}(z) = \frac{\alpha_{\nu}(z)}{\nu+1} f_{\nu+1}(z) \left(\frac{z}{2}\right)^{\frac{5}{2}} + \nu \beta_{\nu}(z) f_{-\nu-1}(z) \left(\frac{z}{2}\right)^{-\frac{5}{2}}, \quad (C.50) $$

$$ H_{i\mu+1}^{(1)}(z) = \frac{\gamma_{\mu}(z)}{i\mu+1} f_{i\mu+1}(z) \frac{z}{2} + i\mu\delta_{\mu}(z) f_{-i\mu-1}(z) \left(\frac{z}{2}\right)^{-1}, \quad (C.51) $$

see Eqs. (10.2.2), (10.4.7) and (10.4.8) of Ref. [101]. This allows one to expand Eqs. (C.7) and (C.8), and one finds

$$ F_{\nu,\mu}(z) = -2\sqrt{2}z^{-3/2} \beta_{\nu}^{*} \left[ \frac{(\nu + i\mu)\gamma_{\mu} + (\nu - i\mu)\delta_{\mu}}{\nu^2 + \mu^2} \right] \\ + \frac{z^{1/2}}{\sqrt{2}} \beta_{\nu}^{*} \left[ \frac{(1 - i\mu)\gamma_{\mu} + (1 + i\mu)\delta_{\mu}}{(1 + \mu^2)(\nu - 1)} \right] \\ + \frac{z^{3/2}}{2\sqrt{2}} \alpha_{\nu}^{*} \left[ \frac{(\nu - i\mu)\gamma_{\mu} + (\nu + i\mu)\delta_{\mu}}{\nu^2 + \mu^2} \right] + \mathcal{O}(z^{5/2}). \quad (C.52) $$

and

$$ G_{\nu,\mu}(z) = -2\sqrt{2}z^{-3/2} \beta_{\nu}^{*} \left[ \frac{(\nu - i\mu)\gamma_{\mu}^{*} + (\nu + i\mu)\delta_{\mu}^{*}}{\nu^2 + \mu^2} \right] \\ + \frac{z^{1/2}}{\sqrt{2}} \beta_{\nu}^{*} \left[ \frac{(1 + i\mu)\gamma_{\mu}^{*} + (1 - i\mu)\delta_{\mu}^{*}}{(1 + \mu^2)(\nu - 1)} \right] \\ + \frac{z^{3/2}}{2\sqrt{2}} \alpha_{\nu}^{*} \left[ \frac{(\nu + i\mu)\gamma_{\mu}^{*} + (\nu - i\mu)\delta_{\mu}^{*}}{\nu^2 + \mu^2} \right] + \mathcal{O}(z^{5/2}). \quad (C.53) $$

Hereafter, to lighten the notation, we have dropped the explicit z-dependence of $\alpha_\nu$, $\beta_\nu$, $\gamma_\mu$ and $\delta_\mu$. This is because,

since $v_\phi$ is close to 3/2 in practice, see Eq. (C.5), this does not affect the power counting in z, see Eqs. (C.44) and (C.45).

It is worth noting that the terms of orders $z^{-3/2}$ and $z^{1/2}$ cancel out in $F_{\nu,\mu}^*(z) + G_{\nu,\mu}(z)$ since $\beta_\nu$ is pure imaginary, see Eq. (C.45). One indeed has

$$ F_{\nu,\mu}^*(z) + G_{\nu,\mu}(z) = \frac{1}{\nu^2 + \mu^2} \frac{z^{3/2}}{\sqrt{2}} \operatorname{Re} (\alpha_\nu) \times [(\nu + i\mu)\gamma_\mu^* + (\nu - i\mu)\delta_\mu^*] + \mathcal{O}(z^{5/2}). \quad (C.54) $$

Let us now expand the coefficients of the transport equation in the super-Hubble limit, i.e. when $z \ll 1$ (but keeping $z_0$ arbitrary).

### D11 coefficient

For $D_{11}$, one finds

$$ \begin{aligned} D_{11} = z^{-7/2} S_{D_{11}}^{(-7/2)}(z, z_0) &+ \frac{\pi}{2} \frac{e^{-\pi \mu_{\chi}}}{\nu_{\phi}^2 + \mu_{\chi}^2} \frac{\lambda^4}{H^4} \\ & \times \left[ \frac{3}{2} |\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}|^2 + 2\mu_{\chi} \Im \left( \gamma_{\mu_{\chi}}^* \delta_{\mu_{\chi}} \right) \right] \frac{k^2}{z^2} \\ & \times + z^{-3/2} S_{D_{11}}^{(-3/2)}(z, z_0), \end{aligned} \quad (C.55) $$

where

$$ S_{D_{11}}^{(-7/2)}(z, z_0) = -\frac{\pi^2 k^2}{2\sqrt{2}} \left(v_\phi - \frac{3}{2}\right) \frac{\lambda^4}{H^4} \\ \times \Im \left\{ \beta_{v_\phi} \left[ F_{v_\phi, \mu_\chi}^*(z_0) + G_{v_\phi, \mu_\chi}(z_0) \right] \right. \\ \left. \times (\gamma_{\mu_\chi} + \delta_{\mu_\chi}) \right\} e^{-\pi \mu_\chi} \quad (C.56) $$

and

$$ S_{D_{11}}^{(-3/2)}(z, z_0) = -\frac{\pi^2}{16\sqrt{2}} \frac{1}{(1 + \mu_{\chi}^2)(v_{\phi} - 1)} \frac{\lambda^4}{H^4} k^2 \\ \times \Im \left[ \beta_{v_{\phi}} \left[ F_{v_{\phi}, \mu_{\chi}}^{*}(z_0) + G_{v_{\phi}, \mu_{\chi}}(z_0) \right] \right. \\ \times \left( \gamma_{\mu_{\chi}}(\mu_{\chi} + i) \{ \mu_{\chi}(-7 + 2v_{\phi}) \right. \\ \left. + i[10 + v_{\phi}(-7 + 2v_{\phi})] \} \right. \\ \left. + \delta_{\mu_{\chi}}(\mu_{\chi} - i) \{ \mu_{\chi}(-7 + 2v_{\phi}) \right. \\ \left. - i[10 + v_{\phi}(-7 + 2v_{\phi})] \} \right) \right] e^{-\pi \mu_{\chi}} $$
$$ + \frac{\pi^2}{16\sqrt{2}} \left(\frac{3}{2} + \nu_{\varphi}\right) \frac{\lambda^4}{H^4} \frac{k^2}{z^{1/2}} \\ \times \mathfrak{Im} \left\{ \left[ -\alpha_{v_{\varphi}}^* F_{v_{\varphi}, \mu_{\chi}}^* (z_0) + \alpha_{v_{\varphi}} G_{v_{\varphi}, \mu_{\chi}} (z_0) \right] \right. \\ \left. (\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}) \right\} e^{-\pi \mu_{\chi}} + \mathcal{O}(z^{1/2}) \quad (C.57) $$

are spurious contributions, i.e. they arise from the term $F_{D_{11}}(z, z_0)$ in Eq. (C.14).

## D$_{12}$ coefficient

One finds

$$ \mathbf{D}_{12} = z^{-5/2} S_{\mathbf{D}_{12}}^{(-5/2)}(z, z_0) \\ + \frac{\pi}{4} \frac{|\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}|^2}{\nu_{\varphi}^2 + \mu_{\chi}^2} \\ \times \frac{\lambda^4}{H^4} e^{-\pi \mu_{\chi}} \frac{k}{z} + z^{-1/2} S_{\mathbf{D}_{12}}^{(-1/2)}(z, z_0), \quad (C.58) $$

where

$$ S_{\mathbf{D}_{12}}^{(-5/2)}(z, z_0) = \frac{\pi^2}{4\sqrt{2}} \frac{\lambda^4}{H^4} k \mathfrak{Im} \left\{ \beta_{v_{\varphi}} \left[ F_{v_{\varphi}, \mu_{\chi}}^{*}(z_0) + G_{v_{\varphi}, \mu_{\chi}}(z_0) \right] (\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}) \right\} e^{-\pi \mu_{\chi}} \quad (C.59) $$

and

$$ S_{\mathbf{D}_{12}}^{(-1/2)}(z, z_0) = \frac{\pi^2 k}{16\sqrt{2}} \frac{e^{-\pi \mu_{\chi}}}{(1 + \mu_{\chi}^2)(\nu_{\varphi} - 1)} \frac{\lambda^4}{H^4} \\ \times \mathfrak{Im} \left\{ \beta_{v_{\varphi}} \left[ F_{v_{\varphi}, \mu_{\chi}}^{*}(z_0) + G_{v_{\varphi}, \mu_{\chi}}(z_0) \right] \right. \\ \left. \times \left[ (1 + \mu_{\chi}^2)(\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}) + (1 - i\mu_{\chi}) \right. \right. \\ \left. (1 - \nu_{\varphi})\gamma_{\mu_{\chi}} + (1 + i\mu_{\chi})(1 - \nu_{\varphi})\delta_{\mu_{\chi}} \right] \right\} \quad (C.60) $$

are again spurious contributions.

## Δ$_{11}$ coefficient

For Δ$_{11}$, we have

$$ \begin{aligned} \mathbf{\Delta}_{11} = & z^{-7/2} S_{\mathbf{\Delta}_{11}}^{(-7/2)}(z, z_0) \\ & - \frac{\pi}{2} \frac{\mu_{\chi} e^{-\pi \mu_{\chi}}}{\nu_{\varphi}^2 + \mu_{\chi}^2} \frac{\lambda^4}{H^4} \left( |\gamma_{\mu_{\chi}}|^2 - |\delta_{\mu_{\chi}}|^2 \right) \frac{k^2}{z^2} \\ & + z^{-3/2} S_{\mathbf{\Delta}_{11}}^{(-3/2)}(z, z_0), \end{aligned} \quad (C.61) $$

where

$$ \begin{aligned} S_{\mathbf{\Delta}_{11}}^{(-7/2)}(z, z_0) = & \frac{\pi^2 k^2}{2\sqrt{2}} \left( \nu_{\varphi} - \frac{3}{2} \right) \frac{\lambda^4}{H^4} \\ & \times \mathfrak{Im} \left\{ \beta_{v_{\varphi}} \left[ F_{v_{\varphi}, \mu_{\chi}}^{*}(z_0) + G_{v_{\varphi}, \mu_{\chi}}(z_0) \right] \right. \\ & \left. (\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}) \right\} e^{-\pi \mu_{\chi}} \end{aligned} \quad (C.62) $$

and

$$ \begin{aligned} S_{\mathbf{\Delta}_{11}}^{(-3/2)}(z, z_0) = & \frac{\pi^2 k^2}{16\sqrt{2}} \frac{1}{(1 + \mu_{\chi}^2)(\nu_{\varphi} - 1)} \frac{\lambda^4}{H^4} \\ & \times \mathfrak{Im} \left[ \beta_{v_{\varphi}} \left[ F_{v_{\varphi}, \mu_{\chi}}^{*}(z_0) + G_{v_{\varphi}, \mu_{\chi}}(z_0) \right] \right. \\ & \times \left( \gamma_{\mu_{\chi}}(\mu_{\chi} + i) \{ \mu_{\chi}(-7 + 2\nu_{\varphi}) \right. \\ & \left. + i[10 + \nu_{\varphi}(-7 + 2\nu_{\varphi})] \} \right. \\ & \left. + \delta_{\mu_{\chi}}(\mu_{\chi} - i) \{ \mu_{\chi}(-7 + 2\nu_{\varphi}) \right. \\ & \left. - i[10 + \nu_{\varphi}(-7 + 2\nu_{\varphi})] \} \right] e^{-\pi \mu_{\chi}} \\ & - \frac{\pi^2}{16\sqrt{2}} \left( \frac{3}{2} + \nu_{\varphi} \right) \frac{\lambda^4}{H^4} \frac{k^2}{z^{1/2}} \\ & \times \mathfrak{Im} \left\{ \left[ -\alpha_{v_{\varphi}}^* F_{v_{\varphi}, \mu_{\chi}}^{*}(z_0) + \alpha_{v_{\varphi}} G_{v_{\varphi}, \mu_{\chi}}(z_0) \right] \right. \\ & \left. \times (\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}) \right\} e^{-\pi \mu_{\chi}} \end{aligned} \quad (C.63) $$

are spurious contributions. It is also worth noting that, in Eq. (C.61), one can simplify

$$ |\gamma_{\mu_{\chi}}|^2 - |\delta_{\mu_{\chi}}|^2 = 2 \frac{e^{\pi \mu}}{\pi \mu}. \quad (C.64) $$

## Δ$_{12}$ coefficient

Finally, for Δ$_{12}$, one obtains

$$ \begin{aligned} \mathbf{\Delta}_{12}(z) = & z^{-5/2} S_{\mathbf{\Delta}_{12}}^{(-5/2)}(z, z_0) \\ & + z^{-1/2} S_{\mathbf{\Delta}_{12}}^{(-1/2)}(z, z_0) \end{aligned} \quad (C.65) $$

which only contains spurious terms as shown in Eq. (C.34), given by

$$ S_{\mathbf{\Delta}_{12}}^{(-5/2)}(z, z_0) = -\frac{\pi^2 k}{4\sqrt{2}} \frac{\lambda^4}{H^4} \mathfrak{Im} \left\{ \beta_{v_{\varphi}} \left[ F_{v_{\varphi}, \mu_{\chi}}^{*}(z_0) + G_{v_{\varphi}, \mu_{\chi}}(z_0) \right] (\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}) \right\} e^{-\pi \mu_{\chi}} \quad (C.66) $$
and

$$
\begin{align}
S_{\Delta_{12}}^{(-1/2)} (z, z_0) & = - \frac{\pi^2 k}{16\sqrt{2}} \frac{e^{-\pi \mu_{\chi}}}{(1 + \mu_{\chi}^2)(v_{\varphi} - 1)} \frac{\lambda^4}{H^4} \nonumber \\
& \quad \times \Re \left\{ \beta_{v_{\varphi}} \left[ F_{v_{\varphi}, \mu_{\chi}}^*(z_0) + G_{v_{\varphi}, \mu_{\chi}}(z_0) \right] \right. \tag{C.67} \\
& \qquad \left. \times \left[ (1 + \mu_{\chi}^2)(\gamma_{\mu_{\chi}} + \delta_{\mu_{\chi}}) + (1 - i\mu_{\chi})(1 - v_{\varphi})\gamma_{\mu_{\chi}} \right. \right. \nonumber \\
& \qquad \left. \left. + (1 + i\mu_{\chi})(1 - v_{\varphi})\delta_{\mu_{\chi}} \right] \right\}. \nonumber
\end{align}
$$

Appendix D: Comparison between TCL and perturbation theory in the curved-space Caldeira–Leggett model

In this appendix, we compare Standard Perturbation Theory (SPT) to the perturbative solutions of the TCL master equation, in the context of the curved-space Caldeira–Leggett model introduced in Sect. 3. This will allow us to exhibit a concrete manifestation of the generic statement proven in Sect. 2.4, that TCL<sub>n</sub> solved perturbatively at order n coincides with SPT<sub>n</sub>.

D.1 Perturbation theory

The two-field system detailed in Sect. 3.1 being linear, the field operators admit a decomposition of the form

$$
\begin{align}
\hat{v}_{\varphi}(\eta) = v_{\varphi\varphi}(\eta)\hat{a}_{\varphi} + v_{\varphi\varphi}^{*}(\eta)\hat{a}_{\varphi}^{\dagger} & \nonumber \\
& + v_{\varphi\chi}(\eta)\hat{a}_{\chi} + v_{\varphi\chi}^{*}(\eta)\hat{a}_{\chi}^{\dagger}, \tag{D.1}
\end{align}
$$

$$
\begin{align}
\hat{v}_{\chi}(\eta) = v_{\chi\varphi}(\eta)\hat{a}_{\varphi} + v_{\chi\varphi}^{*}(\eta)\hat{a}_{\varphi}^{\dagger} & \nonumber \\
& + v_{\chi\chi}(\eta)\hat{a}_{\chi} + v_{\chi\chi}^{*}(\eta)\hat{a}_{\chi}^{\dagger}, \tag{D.2}
\end{align}
$$

where (â<sub>φ</sub>; â<sub>φ</sub>†) and (â<sub>χ</sub>; â<sub>χ</sub>†) are the creation and annihilation
operators of the φ and χ quanta respectively. This generalises
the decomposition (3.13) to the case where fields interact and
exchange quanta. A similar decomposition can be introduced
for the momenta operators p̂<sub>φ</sub> and p̂<sub>χ</sub>, where the Hamilto-
nian (3.6) and (3.7) gives the mode functions

$$
p_{ij}(\eta) = v'_{ij}(\eta) - \frac{a'}{a} v_{ij}(\eta) \tag{D.3}
$$

for i, j ∈ {φ, χ}. Using Heisenberg’s equations, one finds
that the mode functions evolve according to

$$
v_{ij}'' + \omega_i^2(\eta)v_{ij} = -\lambda^2 a^2(\eta)v_{ij}, \quad (D.4)
$$

where we have introduced ω₂φ(η) ≡ k² + m²a²(η) − a''/a
and ω₂χ(η) ≡ k² + M²a²(η) − a''/a, and where ī = χ
when i = φ and ī = φ when i = χ. This constitutes a
set of coupled differential equations, where the coupling is
mediated by λ². It can thus be solved perturbatively in λ.

* **Zeroth order:** The right-hand side of Eq. (D.4) vanishes, hence the uncoupled dynamics is recovered, namely $v_{ii}^{(0)}(\eta) = v_i(\eta)$ and $v_{ii}^{(0)}(\eta) = 0$, where $v_\varphi$ and $v_\chi$ are the free-field mode functions [i.e. they are given by Eq. (3.15) if one replaces $v_\ell$ by $v_\varphi$ and $\mu_h$ by $\mu_\chi$]. One also has $p_{ii}^{(0)}(\eta) = p_i(\eta)$ and $p_{ii}^{(0)}(\eta) = 0$.

* **First order:** At first order, the right-hand side of Eq. (D.4) needs to be replaced with the zeroth-order solution. This does not change the diagonal mode functions $v_{ii}^{(1)}(\eta) = v_{ii}^{(0)}(\eta)$ and $p_{ii}^{(1)}(\eta) = p_{ii}^{(0)}(\eta)$, while the cross mode functions now obey $v_{ii}^{(1)\prime\prime} + \omega_i^2 v_{ii}^{(1)} = \lambda^2 a^2 v_i$. Using the Green's functions of the homogeneous (hence uncoupled) system of differential equation, $g_i(\eta, \eta') = 2\Re[m[v_i(\eta)v_i^*(\eta')]]$, this gives rise to

$$
\begin{equation}
\begin{split}
v_{ii}^{(1)}(\eta) = {}& -2\lambda^2 \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \\
& \times \Re[m[v_i(\eta)v_i^*(\eta_1)]v_i^*(\eta_1)].
\end{split}
\tag{D.5}
\end{equation}
$$

Using Eq. (D.3), this leads to

$$
\begin{equation}
\begin{split}
p_{ii}^{(1)}(\eta) = {}& -2\lambda^2 \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \\
& \times \Re[m[p_i(\eta)v_i^*(\eta_1)]v_i^*(\eta_1)].
\end{split}
\tag{D.6}
\end{equation}
$$

* **Second order:** At second order, Eq. (D.4) is sourced by the first-order solution, so the diagonal mode functions obey $v_{ii}^{(2)\prime\prime} + \omega_i^2 v_{ii}^{(2)} = -\lambda^2 a^2 v_{ii}^{(1)}$. Using again the homogeneous Green functions, together with Eq. (D.5), this gives rise to

$$
\begin{equation}
\begin{split}
v_{ii}^{(2)}(\eta) ={}& v_i(\eta) + 4\lambda^4 \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \\
& \times \int_{\eta_0}^{\eta_1} d\eta_2 a^2(\eta_2) \Re[m[v_i(\eta)v_i^*(\eta_1)] \\
& \times \Re[m[v_i^*(\eta_1)v_i^*(\eta_2)] v_i(\eta_2)].
\end{split}
\tag{D.7}
\end{equation}
$$

Using Eq. (D.3), this leads to

$$
\begin{equation}
\begin{split}
p_{ii}^{(2)}(\eta) ={}& p_i(\eta) + 4\lambda^4 \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \\
& \times \int_{\eta_0}^{\eta_1} d\eta_2 a^2(\eta_2) \Re[m[p_i(\eta)v_i^*(\eta_1)] \\
& \times \Re[m[v_i^*(\eta_1)v_i^*(\eta_2)] v_i(\eta_2)].
\end{split}
\tag{D.8}
\end{equation}
$$

One may also compute the cross mode functions, and
carry on the expansion, but that would lead to subdomi-
nant corrections to the power spectra.

The covariance matrix can be computed using Eq. (3.18), and
onefinds
$$ \Sigma_{\varphi\varphi}(\eta) = \begin{pmatrix} |v_{\varphi\varphi}(\eta)|^2 + |v_{\varphi\chi}(\eta)|^2 & \Re e [v_{\varphi\varphi}(\eta) p_{\varphi\varphi}^*(\eta)] + \Re e [v_{\varphi\chi}(\eta) p_{\varphi\chi}^*(\eta)] \\ \Re e [v_{\varphi\varphi}(\eta) p_{\varphi\varphi}^*(\eta)] + \Re e [v_{\varphi\chi}(\eta) p_{\varphi\chi}^*(\eta)] & |p_{\varphi\varphi}(\eta)|^2 + |p_{\varphi\chi}(\eta)|^2 \end{pmatrix}, \quad (D.9) $$

$$ \Sigma_{\chi\chi}(\eta) = \begin{pmatrix} |v_{\chi\chi}(\eta)|^2 + |v_{\chi\varphi}(\eta)|^2 & \Re e [v_{\chi\chi}(\eta) p_{\chi\chi}^*(\eta)] + \Re e [v_{\chi\varphi}(\eta) p_{\chi\varphi}^*(\eta)] \\ \Re e [v_{\chi\chi}(\eta) p_{\chi\chi}^*(\eta)] + \Re e [v_{\chi\varphi}(\eta) p_{\chi\varphi}^*(\eta)] & |p_{\chi\chi}(\eta)|^2 + |p_{\chi\varphi}(\eta)|^2 \end{pmatrix}, \quad (D.10) $$

$$ \Sigma_{\varphi\chi}(\eta) = \begin{pmatrix} \Re e [v_{\varphi\varphi}(\eta) v_{\chi\varphi}^*(\eta)] + \Re e [v_{\chi\chi}(\eta) v_{\varphi\chi}^*(\eta)] & \Re e [v_{\varphi\varphi}(\eta) p_{\chi\varphi}^*(\eta)] + \Re e [p_{\chi\chi}(\eta) v_{\varphi\chi}^*(\eta)] \\ \Re e [p_{\varphi\varphi}(\eta) v_{\chi\varphi}^*(\eta)] + \Re e [v_{\chi\chi}(\eta) p_{\varphi\chi}^*(\eta)] & \Re e [p_{\varphi\varphi}(\eta) p_{\chi\varphi}^*(\eta)] + \Re e [p_{\chi\chi}(\eta) p_{\varphi\chi}^*(\eta)] \end{pmatrix}. \quad (D.11) $$

By inserting the mode functions obtained above into these expressions, one obtains the first perturbative corrections to the power spectra. For the configuration-configuration power spectrum of the $\varphi$ field, one finds

$$ \Sigma_{\varphi\varphi,11}^{(2)}(\eta) = |v_{\varphi\varphi}^{(0)}(\eta)|^2 + |v_{\varphi\chi}^{(1)}(\eta)|^2 + 2\Re e \left[ v_{\varphi\varphi}^{(2-0)}(\eta) v_{\varphi\varphi}^{(0)*}(\eta) \right], \quad (D.12) $$

where we have introduced the short-hand notation $v_{\varphi\varphi}^{(2-0)}(\eta) = v_{\varphi\varphi}^{(2)}(\eta) - v_{\varphi\varphi}^{(0)}(\eta)$, which selects the terms of order $\lambda^2$ in $v_{\varphi\varphi}^{(2)}(\eta)$. This gives rise to

$$ \begin{aligned} \Sigma_{\varphi\varphi,11}^{(2)}(\eta) &= |v_{\varphi}(\eta)|^2 \\ &+ 4\lambda^4 \left| \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \Im [v_{\varphi}(\eta) v_{\varphi}^*(\eta_1)] v_{\chi}(\eta_1) \right|^2 \\ &+ 8\lambda^4 \Re e \left\{ v_{\varphi}(\eta) \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \right. \\ &\qquad \left. \times \int_{\eta_0}^{\eta_1} d\eta_2 a^2(\eta_2) \Im [v_{\varphi}(\eta) v_{\varphi}^*(\eta_1)] \right. \\ &\qquad \left. \times \Im [v_{\chi}(\eta_1) v_{\chi}^*(\eta_2)] v_{\varphi}^*(\eta_2) \right\}. \end{aligned} \quad (D.13) $$

For the configuration-momentum power spectrum, one obtains

$$ \Sigma_{\varphi\varphi,12}^{(2)}(\eta) = \begin{aligned}[t] & \Re e \left[ v_{\varphi\varphi}^{(0)}(\eta) p_{\varphi\varphi}^{(0)*}(\eta) \right] \\ & + \Re e \left[ v_{\varphi\chi}^{(1)}(\eta) p_{\varphi\chi}^{(1)*}(\eta) \right] \\ & + \Re e \left[ v_{\varphi\varphi}^{(0)}(\eta) p_{\varphi\varphi}^{(2-0)*}(\eta) \right. \\ & \qquad \left. + v_{\varphi\varphi}^{(2-0)}(\eta) p_{\varphi\varphi}^{(0)*}(\eta) \right], \end{aligned} \quad (D.14) $$

namely

$$ \begin{aligned} \Sigma_{\varphi\varphi,12}^{(2)}(\eta) = & \Re e [v_{\varphi}(\eta) p_{\varphi}^*(\eta)] \\ & + 4\lambda^4 \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \Im [v_{\varphi}(\eta) v_{\varphi}^*(\eta')] v_{\chi}(\eta') \\ & \times \int_{\eta_0}^{\eta} d\eta'' a^2(\eta'') \Im [p_{\varphi}(\eta) v_{\varphi}^*(\eta'')] v_{\chi}(\eta'') \\ & + 4\lambda^4 \Re e \left\{ v_{\varphi}(\eta) \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \int_{\eta_0}^{\eta_1} d\eta_2 a^2(\eta_2) \right. \\ & \qquad \left. \times \Im [p_{\varphi}(\eta) v_{\varphi}^*(\eta_1)] \Im [v_{\chi}(\eta_1) v_{\chi}^*(\eta_2)] v_{\varphi}^*(\eta_2) \right\} \\ & + 4\lambda^4 \Re e \left\{ p_{\varphi}(\eta) \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \int_{\eta_0}^{\eta_1} d\eta_2 a^2(\eta_2) \right. \\ & \qquad \left. \times \Im [v_{\varphi}(\eta) v_{\varphi}^*(\eta_1)] \Im [v_{\chi}(\eta_1) v_{\chi}^*(\eta_2)] v_{\varphi}^*(\eta_2) \right\}. \end{aligned} \quad (D.15) $$

Finally, for the momentum-momentum power spectrum, one has

$$ \Sigma_{\varphi\varphi,22}^{(2)}(\eta) = |p_{\varphi\varphi}^{(0)}(\eta)|^2 + |p_{\varphi\chi}^{(1)}(\eta)|^2 + 2\Re e \left[ p_{\varphi\varphi}^{(2-0)}(\eta) p_{\varphi\varphi}^{(0)*}(\eta) \right], \quad (D.16) $$

which leads to

$$ \begin{aligned} \Sigma_{\varphi\varphi,22}^{(2)}(\eta) = & |p_{\varphi}(\eta)|^2 + 4\lambda^4 \left| \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \Im [p_{\varphi}(\eta) v_{\varphi}^*(\eta_1)] v_{\chi}(\eta_1) \right|^2 \\ & + 8\lambda^4 \Re e \left\{ p_{\varphi}(\eta) \int_{\eta_0}^{\eta} d\eta_1 a^2(\eta_1) \int_{\eta_0}^{\eta_1} d\eta_2 a^2(\eta_2) \right. \\ & \qquad \left. \times \Im [p_{\varphi}(\eta) v_{\varphi}^*(\eta_1)] \Im [v_{\chi}(\eta_1) v_{\chi}^*(\eta_2)] v_{\varphi}^*(\eta_2) \right\}. \end{aligned} \quad (D.17) $$
## D.2 Perturbative solution of TCL

Let us start with the TCL$_2$ master equation written in the form

$$ \frac{d\tilde{\rho}_{\text{red}}}{d\eta} = -\lambda^4 a^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \times \left\{ \begin{aligned} & [\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta')\tilde{\rho}_{\text{red}}(\eta) - \tilde{v}_{\varphi}(\eta')\tilde{\rho}_{\text{red}}(\eta)\tilde{v}_{\varphi}(\eta)] K^{>}(\eta, \eta') \\ & - [\tilde{v}_{\varphi}(\eta)\tilde{\rho}_{\text{red}}(\eta)\tilde{v}_{\varphi}(\eta') - \tilde{\rho}_{\text{red}}(\eta)\tilde{v}_{\varphi}(\eta')\tilde{v}_{\varphi}(\eta)] K^{>*(}\eta, \eta') \end{aligned} \right\}. \quad (\text{D.18}) $$

This equation was obtained in Eq. (A.17) from microphysical considerations and is just a convenient rewriting of Eq. (3.25). We want to solve it at order $\lambda^4$, i.e. drop all contributions of higher order. Since the right-hand side is already proportional to $\lambda^4$, this implies that it can be evaluated in the free theory, where $\tilde{\rho}_{\text{red}}(\eta) \simeq \tilde{\rho}_{\text{red}}(\eta_0)$. One can thus integrate Eq. (D.18), which leads to

$$ \tilde{\rho}_{\text{red}}^{(2)}(\eta) = \tilde{\rho}_{\text{red}}(\eta_0) - \lambda^4 \int_{\eta_0}^{\eta} d\eta' a^2(\eta') \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') \times \left\{ \begin{aligned} & [\tilde{v}_{\varphi}(\eta')\tilde{v}_{\varphi}(\eta'')\tilde{\rho}_{\text{red}}(\eta_0) - \tilde{v}_{\varphi}(\eta'')\tilde{\rho}_{\text{red}}(\eta_0)\tilde{v}_{\varphi}(\eta')] v_{\chi}(\eta')v_{\chi'}^{*}(\eta'') \\ & - [\tilde{v}_{\varphi}(\eta')\tilde{\rho}_{\text{red}}(\eta_0)\tilde{v}_{\varphi}(\eta'')] - \tilde{\rho}_{\text{red}}(\eta_0)\tilde{v}_{\varphi}(\eta'')\tilde{v}_{\varphi}(\eta')] v_{\chi}^{*}(\eta')v_{\chi'}(\eta'') \end{aligned} \right\}, \quad (\text{D.19}) $$

where we have used that the memory kernels are related to the free mode functions via Eq. (A.15).

Let us now compute the entries of the covariance matrix using this expression for $\tilde{\rho}_{\text{red}}^{(2)}$. The configuration-configuration power spectrum reads

$$ \Sigma_{\text{TCL},11}^{(2)}(\eta) = \text{Tr} \left[ \tilde{v}_{\varphi}(\eta) \tilde{v}_{\varphi}(\eta) \tilde{\rho}_{\text{red}}^{(2)}(\eta) \right], \quad (\text{D.20}) $$

that is

$$ \begin{aligned} \Sigma_{\text{TCL},11}^{(2)}(\eta) ={}& \text{Tr} [\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta)\tilde{\rho}_{\text{red}}(\eta_0)] \\ & - \lambda^4 \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{\chi}(\eta') \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') v_{\chi}^{*}(\eta'') \\ & \times \{ \text{Tr} [\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta')\tilde{v}_{\varphi}(\eta'')\tilde{\rho}_{\text{red}}(\eta_0)] \\ & - \text{Tr} [\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta'')\tilde{\rho}_{\text{red}}(\eta_0)\tilde{v}_{\varphi}(\eta')] \\ & + \lambda^4 \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{\chi}^{*}(\eta') \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') v_{\chi}(\eta'') \\ & \times \{ \text{Tr} [\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta')\tilde{\rho}_{\text{red}}(\eta_0)\tilde{v}_{\varphi}(\eta'')] \\ & - \text{Tr} [\tilde{v}_{\varphi}(\eta)\tilde{v}_{\varphi}(\eta)\tilde{\rho}_{\text{red}}(\eta_0)\tilde{v}_{\varphi}(\eta'')\tilde{v}_{\varphi}(\eta')] \end{aligned} $$

(D.21)

Since the initial state is the Bunch-Davies vacuum, $\tilde{\rho}_{\text{red}}(\eta_0) = |\phi\rangle\langle\phi|$, using the mode-function decomposition (3.27) one obtains

$$ \Sigma_{\text{TCL},11}^{(2)}(\eta) = |v_{\varphi}(\eta)|^2 - 4\lambda^4 \Re \left[ v_{\varphi}^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{\varphi}^*(\eta') v_{\chi}(\eta') \right. \\ \times \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') v_{\varphi}^*(\eta'') v_{\chi'}^*(\eta'') \\ \left. - |v_{\varphi}(\eta)|^2 \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{\varphi}(\eta') v_{\chi}(\eta') \times \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') v_{\varphi}^*(\eta'') v_{\chi'}^*(\eta'') \right]. \quad (\text{D.22}) $$

This expression matches Eq. (D.13), as can be shown by expanding the real and imaginary parts and relabeling the integration domain. Following the same method, one finds

$$ \begin{aligned} \Sigma_{\text{TCL},12}^{(2)}(\eta) ={}& \Re \left[ v_{\varphi}(\eta) p_{\varphi}^{*}(\eta) \right] \\ & - 4\lambda^4 \Re \left\{ v_{\varphi}(\eta) p_{\varphi}(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{\varphi}^*(\eta') v_{\chi}(\eta') \right. \\ & \times \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') v_{\varphi}^*(\eta'') v_{\chi'}^*(\eta'') \\ & \left. - \Re \left[ v_{\varphi}(\eta) p_{\varphi}^*(\eta) \right] \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{\varphi}(\eta') v_{\chi}(\eta') \right. \\ & \left. \times \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') v_{\varphi}^*(\eta'') v_{\chi'}^*(\eta'') \right\}, \end{aligned} \quad (\text{D.23}) $$

which can be shown to match Eq. (D.15), and

$$ \begin{aligned} \Sigma_{\text{TCL},22}^{(2)}(\eta) ={}& |p_{\varphi}(\eta)|^2 \\ & - 4\lambda^4 \Re \left[ p_{\varphi}^2(\eta) \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{\varphi}^*(\eta') v_{\chi}(\eta') \right. \\ & \times \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') v_{\varphi}^*(\eta'') v_{\chi'}^*(\eta'') \\ & \left. - |p_{\varphi}(\eta)|^2 \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{\varphi}(\eta') v_{\chi}(\eta') \right. \\ & \left. \times \int_{\eta_0}^{\eta'} d\eta'' a^2(\eta'') v_{\varphi}^*(\eta'') v_{\chi'}^*(\eta'') \right], \end{aligned} \quad (\text{D.24}) $$

which can be shown to match Eq. (D.17).
## Appendix E: Comparison with other late-time resummation techniques

In this section, we compare TCL with the late-time resummation technique proposed in Ref. [44] and also studied in Ref. [59]. The idea is to keep track of the growing mode only, in order to simplify the analysis in the late-time limit. As we will make clear, the method also implicitly performs an additional layer of approximation compared to TCL, which makes it less efficient.

The starting point is to rewrite the free mode function

$$v_{\varphi}(z) = \frac{1}{2} \sqrt{\frac{\pi z}{k}} e^{i \frac{\pi}{2} \left(v_{\varphi} + \frac{1}{2}\right)} H_{v_{\varphi}}^{(1)}(z), \quad (E.1)$$

where we recall that $z = -k\eta$, as (see Eq. (10.4.3) of Ref. [101])

$$v_{\varphi}(z) = e^{i \frac{\pi}{2} \left(v_{\varphi} + \frac{1}{2}\right)} \frac{v_{-}(z) + i v_{+}(z)}{\sqrt{2}} \quad (E.2)$$

where

$$v_{+}(z) = \sqrt{\frac{\pi z}{2k}} Y_{v_{\varphi}}(z) \quad \text{and} \\ v_{-}(z) = \sqrt{\frac{\pi z}{2k}} J_{v_{\varphi}}(z) \quad (E.3)$$

are real functions. Here, $J_ν$ and $Y_ν$ are the Bessel functions of the first and second kind respectively, and of order $v$. The reason why this decomposition is convenient is that $v_-$ corresponds to the cosmological “decaying mode” [i.e. $v_-(\eta)$ decreases on super-Hubble scales], while $v_+$ stands for the growing mode. Let us recall that the heavy-field mode function cannot be divided into a growing mode and a decaying mode, since both modes oscillate with similar amplitude on super-Hubble scales.

In the interaction picture, where operators evolve as in the free theory, the mode-function expansion (3.27) of the field operators can then be written as

$$\tilde{v}_{\varphi}(\eta) = v_{\varphi}(\eta) \hat{a}_{\varphi} + v_{\varphi}^{*}(\eta) \hat{a}_{\varphi}^{\dagger} \quad (E.4)$$

$$= v_{-}(\eta) \hat{P}_{\varphi} + v_{+}(\eta) \hat{Q}_{\varphi}, \quad (E.5)$$

where

$$\hat{P}_{\varphi} = \frac{1}{\sqrt{2}} \left[ e^{i \frac{\pi}{2} \left(v_{\varphi} + \frac{1}{2}\right)} \hat{a}_{\varphi} + e^{-i \frac{\pi}{2} \left(v_{\varphi} + \frac{1}{2}\right)} \hat{a}_{\varphi}^{\dagger} \right], \quad (E.6)$$

$$\hat{Q}_{\varphi} = \frac{i}{\sqrt{2}} \left[ e^{i \frac{\pi}{2} \left(v_{\varphi} + \frac{1}{2}\right)} \hat{a}_{\varphi} - e^{-i \frac{\pi}{2} \left(v_{\varphi} + \frac{1}{2}\right)} \hat{a}_{\varphi}^{\dagger} \right]. \quad (E.7)$$

One can check that they constitute a set of canonical variables since $[\hat{Q}_{\varphi}, \hat{P}_{\varphi}] = i$.

The idea proposed in Refs. [44,59] is to insert the decomposition (E.5) into the TCL$_2$ master equation (3.25) in order

to identify the leading late-time contribution. One finds

$$\frac{d\hat{\rho}_{\text{red}}^{\text{IR}}}{d\eta} = -\lambda^4 a^2(\eta) \left\{ \begin{aligned} & v_{-}(\eta) X_{-}^{*}(\eta) v_{\chi}(\eta) \\ & \times [\hat{P}_{\varphi}^2 \tilde{\rho}_{\text{red}}(\eta) - \hat{P}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) \hat{P}_{\varphi}] \\ & + v_{-}(\eta) X_{-}(\eta) v_{\chi}^{*}(\eta) \\ & \times [\tilde{\rho}_{\text{red}}(\eta) \hat{P}_{\varphi}^2 - \hat{P}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) \hat{P}_{\varphi}] \\ & + v_{+}(\eta) X_{-}^{*}(\eta) v_{\chi}(\eta) \\ & \times [\hat{Q}_{\varphi} \hat{P}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) - \hat{P}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) \hat{Q}_{\varphi}] \\ & + v_{+}(\eta) X_{-}(\eta) v_{\chi}^{*}(\eta) \\ & \times [\tilde{\rho}_{\text{red}}(\eta) \hat{P}_{\varphi} \hat{Q}_{\varphi} - \hat{Q}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) \hat{P}_{\varphi}] \\ & + v_{-}(\eta) X_{+}^{*}(\eta) v_{\chi}(\eta) \\ & \times [\hat{P}_{\varphi} \hat{Q}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) - \hat{Q}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) \hat{P}_{\varphi}] \\ & + v_{-}(\eta) X_{+}(\eta) v_{\chi}^{*}(\eta) \\ & \times [\tilde{\rho}_{\text{red}}(\eta) \hat{Q}_{\varphi} \hat{P}_{\varphi} - \hat{P}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) \hat{Q}_{\varphi}] \\ & + v_{+}(\eta) X_{+}^{*}(\eta) v_{\chi}(\eta) \\ & \times [\hat{Q}_{\varphi}^2 \tilde{\rho}_{\text{red}}(\eta) - \hat{Q}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) \hat{Q}_{\varphi}] \\ & + v_{+}(\eta) X_{+}(\eta) v_{\chi}^{*}(\eta) \\ & \times [\tilde{\rho}_{\text{red}}(\eta) \hat{Q}_{\varphi}^2 - \hat{Q}_{\varphi} \tilde{\rho}_{\text{red}}(\eta) \hat{Q}_{\varphi}] \end{aligned} \right\}, \quad (E.8)$$

where

$$X_{+}(\eta) = \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{+}(\eta') v_{\chi}(\eta'), \quad (E.9)$$

$$X_{-}(\eta) = \int_{\eta_0}^{\eta} d\eta' a^2(\eta') v_{-}(\eta') v_{\chi}(\eta'). \quad (E.10)$$

The authors of Refs. [44,59] argue that dropping the decaying-mode contributions constitutes a valid approximation in the infrared (IR) limit, and for this reason hereafter we label the quantities computed in this scheme with the superscript “IR”.

In the interaction picture, the configuration–configuration power spectrum reads

$$\langle \tilde{v}_{\varphi}(\eta) \tilde{v}_{\varphi}(\eta) \rangle = v_{-}(\eta) v_{-}(\eta) \langle \hat{P}_{\varphi}^2 \rangle \\ + v_{-}(\eta) v_{+}(\eta) \langle \hat{Q}_{\varphi} \hat{P}_{\varphi} + \hat{P}_{\varphi} \hat{Q}_{\varphi} \rangle \\ + v_{+}(\eta) v_{+}(\eta) \langle \hat{Q}_{\varphi}^2 \rangle \\ \approx v_{+}(\eta) v_{+}(\eta) \langle \hat{Q}_{\varphi}^2 \rangle, \quad (E.11)$$

where in the second line we have neglected the decaying mode contribution. The next step is to compute $\langle \hat{Q}_{\varphi}^2 \rangle (\eta) = \text{Tr} [\hat{Q}_{\varphi}^2 \tilde{\rho}_{\text{red}}(\eta)]$ with the IR master equation (E.8). Upon dif-
ferentiating this expression with respect to time, one obtains

$$ \frac{d\langle\hat{Q}_{\phi}^2\rangle}{d\eta} = \Gamma(\eta)\langle\hat{Q}_{\phi}^2\rangle \qquad (E.12) $$

where

$$ \Gamma(\eta) = 4\lambda^4 a^2(\eta) v_-(\eta) \Im \left[ v_\chi(\eta) X_+^*(\eta) \right], \qquad (E.13) $$

which gives rise to

$$ \langle \hat{Q}_{\phi}^2 \rangle (\eta) = e^{\int_{\eta_*}^{\eta} d\eta' \Gamma(\eta')} \langle \hat{Q}_{\phi}^2 \rangle (\eta_*). \qquad (E.14) $$

Since we are interested in the late-time behaviour of the power spectra, we can assume $-k\eta \ll 1$ and let $\eta_*$ denote the Hubble crossing time, $\eta_* \equiv -1/k$, if the above integral is dominated by its upper bound (hence does not depend much on the choice of the lower bound). If the effect of the interaction with the environment is small in sub-Hubble scales, as argued in Ref. [44] one can evaluate $\langle \hat{Q}_{\phi}^2 \rangle (\eta_*)$ in the free theory, which simply yields

$$ \langle \hat{Q}_{\phi}^2 \rangle (\eta_*) \approx 1. \qquad (E.15) $$

In the super-Hubble limit, using the results derived in Appendix C.3, one can also approximate

$$ X_+^*(z) \approx \frac{\pi}{H^2} \frac{(-1)^{3/4}}{\nu_{\phi}^2 + \mu_{\chi}^2} \times \frac{z^{-\nu_{\phi}}}{\sin(\pi \nu_{\phi}) \Gamma(1 - \nu_{\phi})} \times \left[ \nu_{\phi} \left( \gamma_{\mu_{\chi}}^* + \delta_{\mu_{\chi}}^* \right) - i \mu_{\chi} \left( \gamma_{\mu_{\chi}}^* - \delta_{\mu_{\chi}}^* \right) \right] e^{-\pi \mu_{\chi}}, \qquad (E.16) $$

where $\gamma_\mu$ and $\delta_\mu$ were defined in Eq. (C.45). This leads to

$$ \Gamma(z) = \frac{\pi}{2\nu_{\phi}} \frac{1}{\nu_{\phi}^2 + \mu_{\chi}^2} \frac{\lambda^4 k}{H^4 z} \mu_{\chi} \times \left( |\gamma_{\mu_{\chi}}|^2 - |\delta_{\mu_{\chi}}|^2 \right) e^{-\pi \mu_{\chi}}, \qquad (E.17) $$

where $|\gamma_{\mu_{\chi}}|^2 - |\delta_{\mu_{\chi}}|^2 = 2e^{\pi \mu_{\chi}} / (\pi \mu_{\chi})$. One thus has

$$ \int_{\eta_*}^{\eta} d\eta' \Gamma(\eta') \approx -\frac{1}{\nu_{\phi}} \frac{1}{\nu_{\phi}^2 + \mu_{\chi}^2} \frac{\lambda^4}{H^4} \ln(-k\eta), \qquad (E.18) $$

which one can check does not depend on the detailed choice of $\eta_*$ as announced above. Combining the above results, one obtains

$$ \Sigma_{\text{IR},11}(\eta) = e^{-\frac{1}{\nu_{\phi}} \frac{1}{\nu_{\phi}^2 + \mu_{\chi}^2} \frac{\lambda^4}{H^4} \ln(-k\eta)} $$

<figure><img src="image_7.png" /><figcaption>Fig. 5 Relative error in the three power spectra for TCL_2 (blue curves) and the IR resummation method presented in Appendix E (green curves). The parameters are taken as $m^2 = 10^{-4}H^2$, $M^2 = 10^3 H^2$ and $\lambda^2 = 10^{-3}H^2$</figcaption></figure>

$$ \times |v_{\phi}(\eta)|^2. \qquad (E.19) $$

The configuration–momentum and momentum–momentum power spectra can be computed along similar lines. Starting from

$$ \tilde{p}_{\phi}(\eta) = p_{+}(\eta) \hat{Q}_{\phi} + p_{-}(\eta) \hat{P}_{\phi} \qquad (E.20) $$

and using the fact that $p_+ = v'_+ - (a'/a)v_+$ and $p_- = v'_- - (a'/a)v_-$ are still growing and decaying respectively, one has

$$ \langle \tilde{v}_{\phi}(\eta) \tilde{p}_{\phi}(\eta) \rangle \approx v_{+}(\eta) p_{+}(\eta) \langle \hat{Q}_{\phi}^2 \rangle, \qquad (E.21) $$

$$ \langle \tilde{p}_{\phi}(\eta) \tilde{p}_{\phi}(\eta) \rangle \approx p_{+}(\eta) p_{+}(\eta) \langle \hat{Q}_{\phi}^2 \rangle. \qquad (E.22) $$

This implies that the same correction is obtained for all power spectra, i.e.

$$ \Sigma_{\text{IR},11}(\eta) = e^{-\frac{1}{\nu_{\phi}} \frac{1}{\nu_{\phi}^2 + \mu_{\chi}^2} \frac{\lambda^4}{H^4} \ln(-k\eta)} |v_{\phi}(\eta)|^2, \qquad (E.23) $$

$$ \Sigma_{\text{IR},12}(\eta) = e^{-\frac{1}{\nu_{\phi}} \frac{1}{\nu_{\phi}^2 + \mu_{\chi}^2} \frac{\lambda^4}{H^4} \ln(-k\eta)} \times \Re \left[ v_{\phi}(\eta) p_{\phi}^*(\eta) \right], \qquad (E.24) $$

$$ \Sigma_{\text{IR},22}(\eta) = e^{-\frac{1}{\nu_{\phi}} \frac{1}{\nu_{\phi}^2 + \mu_{\chi}^2} \frac{\lambda^4}{H^4} \ln(-k\eta)} |p_{\phi}(\eta)|^2. \qquad (E.25) $$

These expressions feature manifest resummations over powers of $\ln(a)$, which we now compare with the resummation performed by the TCL_2 master equation. The relative difference between the three power spectra and their exact counterpart is displayed in Fig. 5, both for TCL_2 (blue curves)<fnref n="13" /> and IR (green curves).

<footnote n="13">¹³ Let us note that at late time, the relative error in TCL asymptotes a constant in Fig. 5, hence it is not described by Eq. (4.8). The reason</footnote>
Let us first note that the growth rate of the power spectra is correctly captured in the IR approach, even at strong coupling where the perturbative result usually breaks down. This can be further understood by noting that Eqs. (E.23)–(E.25) take the same form as Eq. (4.6) with

$$v_{\text{IR}} = \frac{1}{2v_{\varphi}} \frac{1}{v_{\varphi}^2 + \mu_{\chi}^2} \left(\frac{\lambda}{H}\right)^4 + v_{\varphi}, \qquad (\text{E.26})$$

while according to footnote 11, in TCL₂ one has

$$v_{\text{LS}} = \frac{3}{2} \sqrt{1 - \left(\frac{2m_{\text{LS}}}{3H}\right)^2} \quad \text{where} \qquad (\text{E.27})$$

$$m_{\text{LS}}^2 = m^2 - \frac{\lambda^4}{M^2 - m^2}$$

and we recall that in the exact theory

$$\begin{gathered} v_{\ell} = \frac{3}{2} \sqrt{1 - \left(\frac{2m_{\ell}}{3H}\right)^2} \quad \text{where} \\ m_{\ell}^2 = \frac{1}{2} \left[ m^2 + M^2 - (M^2 - m^2) \times \sqrt{1 + \left(\frac{2\lambda^2}{M^2 - m^2}\right)^2} \right]. \end{gathered} \qquad (\text{E.28})$$

Since Eqs. (E.26), (E.27) and (E.28) coincide when expanded at first order in $\lambda^4$, one concludes that, at the level of the growth rate, the Lamb-shift renormalisation of the mass is correctly accounted for in the IR approach [44] as for TCL, at least at leading order in the coupling constant. This is similar to the dynamical renormalisation group (DRG) treatment of late-time secular divergences in de Sitter performed in Refs. [56–58], as pointed out in Refs. [44,46,59].

The IR approach however fails to reproduce the overall amplitude of the power spectra beyond the perturbative level, which explains why it does not perform as well as TCL. Let us also note that another disadvantage of the IR method is that it does not allow one to track decoherence, which as explained in Sect. 4.2 is not driven by the growing modes.

## References

1. Planck collaboration, Planck 2018 results. X. Constraints on inflation. Astron. Astrophys. 641, A10 (2020). https://doi.org/10.1051/0004-6361/201833887. arXiv:1807.06211

2. J. Martin, C. Ringeval, V. Vennin, Encyclopædia Inflationaris. Phys. Dark Univ. 5–6, 75 (2014). https://doi.org/10.1016/j.dark.2014.01.003. arXiv:1303.3787

3. L. Amendola et al., Cosmology and fundamental physics with the Euclid satellite. Living Rev. Relativ. 21, 2 (2018). https://doi.org/10.1007/s41114-017-0010-3. arXiv:1606.00180

4. SKA Cosmology SWG collaboration, Overview of Cosmology with the SKA. PoS AASKA14, 016 (2015). https://doi.org/10.22323/1.215.0016. arXiv:1501.04076

5. D.H. Lyth, D. Wands, Generating the curvature perturbation without an inflaton. Phys. Lett. B 524, 5 (2002). https://doi.org/10.1016/S0370-2693(01)01366-1. arXiv:hep-ph/0110002

6. K. Enqvist, R.J. Hardwick, T. Tenkanen, V. Vennin, D. Wands, A novel way to determine the scale of inflation. JCAP 02, 006 (2018). https://doi.org/10.1088/1475-7516/2018/02/006. arXiv:1711.07344

7. C. Ringeval, T. Suyama, T. Takahashi, M. Yamaguchi, S. Yokoyama, Dark energy from primordial inflationary quantum fluctuations. Phys. Rev. Lett. 105, 121301 (2010). https://doi.org/10.1103/PhysRevLett.105.121301. arXiv:1006.0368

8. C. Kiefer, F. Queisser, A.A. Starobinsky, Cosmological constant from decoherence. Class. Quantum Gravity 28, 125022 (2011). https://doi.org/10.1088/0264-9381/28/12/125022. arXiv:1010.5331

9. R.H. Brandenberger, V.F. Mukhanov, T. Prokopec, Entropy of a classical stochastic field and cosmological perturbations. Phys. Rev. Lett. 69, 3606 (1992). https://doi.org/10.1103/PhysRevLett.69.3606. arXiv:astro-ph/9206005

10. A.O. Barvinsky, A.Y. Kamenshchik, C. Kiefer, I.V. Mishakov, Decoherence in quantum cosmology at the onset of inflation. Nucl. Phys. B 551, 374 (1999). https://doi.org/10.1016/S0550-3213(99)00208-4. arXiv:gr-qc/9812043

11. F.C. Lombardo, D. Lopez Nacir, Decoherence during inflation: the generation of classical inhomogeneities. Phys. Rev. D 72, 063506 (2005). https://doi.org/10.1103/PhysRevD.72.063506. arXiv:gr-qc/0506051

12. C. Kiefer, I. Lohmar, D. Polarski, A.A. Starobinsky, Pointer states for primordial fluctuations in inflationary cosmology. Class. Quantum Gravity 24, 1699 (2007). https://doi.org/10.1088/0264-9381/24/7/002. arXiv:astro-ph/0610700

13. P. Martineau, On the decoherence of primordial fluctuations during inflation. Class. Quantum Gravity 24, 5817 (2007). https://doi.org/10.1088/0264-9381/24/23/006. arXiv:astro-ph/0601134

14. C.P. Burgess, R. Holman, D. Hoover, Decoherence of inflationary primordial fluctuations. Phys. Rev. D 77, 063534 (2008). https://doi.org/10.1103/PhysRevD.77.063534. arXiv:astro-ph/0601646

15. T. Prokopec, G.I. Rigopoulos, Decoherence from isocurvature perturbations in inflation. JCAP 11, 029 (2007). https://doi.org/10.1088/1475-7516/2007/11/029. arXiv:astro-ph/0612067

16. E. Nelson, Quantum decoherence during inflation from gravitational nonlinearities. JCAP 03, 022 (2016). https://doi.org/10.1088/1475-7516/2016/03/022. arXiv:1601.03734

17. J. Martin, V. Vennin, Observational constraints on quantum decoherence during inflation. JCAP 05, 063 (2018). https://doi.org/10.1088/1475-7516/2018/05/063. arXiv:1801.09949

18. J. Martin, V. Vennin, Non Gaussianities from quantum decoherence during inflation. JCAP 06, 037 (2018). https://doi.org/10.1088/1475-7516/2018/06/037. arXiv:1805.05609

19. W.H. Zurek, Pointer basis of quantum apparatus: into what mixture does the wave packet collapse? Phys. Rev. D 24, 1516 (1981). https://doi.org/10.1103/PhysRevD.24.1516

20. W.H. Zurek, Environment induced superselection rules. Phys. Rev. D 26, 1862 (1982). https://doi.org/10.1103/PhysRevD.26.1862

21. E. Joos, H. Zeh, The emergence of classical properties through interaction with the environment. Z. Phys. B 59, 223 (1985). https://doi.org/10.1007/BF01725541

22. T.J. Hollowood, J.I. McDonald, Decoherence, discord and the quantum master equation for cosmological perturbations. Phys.

<footnote>is that Eq. (4.8) captures the error in the growth rate, while for the parameters displayed in Fig. 5 the error in the overall amplitude provides the dominant contribution.</footnote>
Rev. D 95, 103521 (2017). https://doi.org/10.1103/PhysRevD.95.103521. arXiv:1701.02235

23. J. Martin, A. Micheli, V. Vennin, Discord and decoherence. JCAP 04, 051 (2022). https://doi.org/10.1088/1475-7516/2022/04/051. arXiv:2112.05037

24. D. Koks, A. Matacz, B. Hu, Entropy and uncertainty of squeezed quantum open systems. Phys. Rev. D 55, 5917 (1997). https://doi.org/10.1103/PhysRevD.55.5917. arXiv:quant-ph/9612016

25. C. Anastopoulos, B.L. Hu, A master equation for gravitational decoherence: probing the textures of spacetime. Class. Quantum Gravity 30, 165007 (2013). https://doi.org/10.1088/0264-9381/30/16/165007. arXiv:1305.5231

26. M. Fukuma, Y. Sakatani, S. Sugishita, Master equation for the Unruh-DeWitt detector and the universal relaxation time in de Sitter space. Phys. Rev. D 89, 064024 (2014). https://doi.org/10.1103/PhysRevD.89.064024. arXiv:1305.0256

27. C. Cheung, P. Creminelli, A.L. Fitzpatrick, J. Kaplan, L. Senatore, The effective field theory of inflation. JHEP 03, 014 (2008). https://doi.org/10.1088/1126-6708/2008/03/014. arXiv:0709.0293

28. X. Chen, Y. Wang, Quasi-single field inflation and non-Gaussianities. JCAP 04, 027 (2010). https://doi.org/10.1088/1475-7516/2010/04/027. arXiv:0911.3380

29. L. Senatore, M. Zaldarriaga, The effective field theory of multifield inflation. JHEP 04, 024 (2012). https://doi.org/10.1007/JHEP04(2012)024. arXiv:1009.2093

30. V. Assassini, D. Baumann, D. Green, L. McAllister, Planck-suppressed operators. JCAP 01, 033 (2014). https://doi.org/10.1088/1475-7516/2014/01/033. arXiv:1304.5226

31. N. Arkani-Hamed, J. Maldacena, Cosmological collider physics. arXiv:1503.08043

32. S. Shandera, N. Agarwal, A. Kamal, Open quantum cosmological system. Phys. Rev. D 98, 083535 (2018). https://doi.org/10.1103/PhysRevD.98.083535. arXiv:1708.00493

33. S. Akhtar, S. Choudhury, S. Chowdhury, D. Goswami, S. Panda, A. Swain, Open quantum entanglement: a study of two atomic system in static patch of de Sitter space. Eur. Phys. J. C 80, 748 (2020). https://doi.org/10.1140/epjc/s10052-020-8302-2. arXiv:1908.09929

34. C.P. Burgess, *Introduction to Effective Field Theory* (Cambridge University Press, Cambridge, 2020). https://doi.org/10.1017/9781139048040

35. S. Banerjee, S. Choudhury, S. Chowdhury, R.N. Das, N. Gupta, S. Panda et al., Indirect detection of cosmological constant from interacting open quantum system. Ann. Phys. 443, 168941 (2022). https://doi.org/10.1016/j.aop.2022.168941. arXiv:2004.13058

36. L. Pinol, S. Aoki, S. Renaux-Petel, M. Yamaguchi, Inflationary flavor oscillations and the cosmic spectroscopy. arXiv:2112.05710

37. J. Oppenheim, C. Sparaciari, B. Šoda, Z. Weller-Davies, Gravitationally induced decoherence vs space-time diffusion: testing the quantum nature of gravity. arXiv:2203.01982

38. G.L. Pimentel, D.-G. Wang, Boostless cosmological collider bootstrap. arXiv:2205.00013

39. S. Jazayeri, S. Renaux-Petel, Cosmological bootstrap in slow motion. arXiv:2205.10340

40. S. Brahma, A. Berera, J. Calderón-Figueroa, Quantum corrections to the primordial tensor spectrum: open EFTs & Markovian decoupling of UV modes. arXiv:2206.05797

41. H.P. Breuer, F. Petruccione, *The Theory of Open Quantum Systems* (Oxford University Press, Oxford, 2002). https://doi.org/10.1093/acprof:oso/9780199213900.001.0001

42. E.A. Calzetta, B.-L.B. Hu, *Nonequilibrium Quantum Field Theory*, Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 2008). https://doi.org/10.1017/CBO9780511535123

43. C.P. Burgess, R. Holman, G. Tasinato, M. Williams, EFT beyond the horizon: stochastic inflation and how primordial quantum fluctuations go classical. JHEP 03, 090 (2015). https://doi.org/10.1007/JHEP03(2015)090. [arXiv:1408.5002]

44. D. Boyanovsky, Effective field theory during inflation: reduced density matrix and its quantum master equation. Phys. Rev. D 92, 023527 (2015). https://doi.org/10.1103/PhysRevD.92.023527. arXiv:1506.07395

45. D. Boyanovsky, Effective field theory during inflation. II. Stochastic dynamics and power spectrum suppression. Phys. Rev. D 93, 043501 (2016). https://doi.org/10.1103/PhysRevD.93.043501. arXiv:1511.06649

46. C. Burgess, R. Holman, G. Tasinato, Open EFTs, IR effects & late-time resummations: systematic corrections in stochastic inflation. JHEP 01, 153 (2016). https://doi.org/10.1007/JHEP01(2016)153. arXiv:1512.00169

47. G. Kaplanek, C.P. Burgess, R. Holman, Qubit heating near a hotspot. JHEP 08, 132 (2021). https://doi.org/10.1007/JHEP08(2021)132. arXiv:2106.10803

48. S. Chaykov, N. Agarwal, S. Bahrami, R. Holman, Loop corrections in Minkowski spacetime away from equilibrium 1: Late-time resummations. arXiv:2206.11288

49. S. Chaykov, N. Agarwal, S. Bahrami, R. Holman, Loop corrections in Minkowski spacetime away from equilibrium 2: finite-time results. arXiv:2206.11289

50. J. Martin, Inflationary perturbations: the cosmological Schwinger effect. Lect. Notes Phys. 738, 193 (2008). https://doi.org/10.1007/978-3-540-74353-8_6. arXiv:0704.3540

51. J.-T. Hsiang, B.-L. Hu, Fluctuation-Dissipation relation for a quantum Brownian oscillator in a parametrically squeezed thermal field. arXiv:2107.13343

52. T. Colas, J. Grain, V. Vennin, Four-mode squeezed states: two-field quantum systems and the symplectic group Sp(4, mathbb{R}). Eur. Phys. J. C 82, 6 (2022). https://doi.org/10.1140/epjc/s10052-021-09922-y. arXiv:2104.14942

53. S. Banerjee, S. Choudhury, S. Chowdhury, J. Knaute, S. Panda, K. Shirish, Thermalization phenomena in quenched quantum brownian motion in De Sitter space. arXiv:2104.10692

54. G. Lindblad, On the generators of quantum dynamical semigroups. Commun. Math. Phys. 48, 119 (1976). https://doi.org/10.1007/BF01608499

55. G. Kaplanek, E. Tjoa, Mapping Markov: on effective master equations for two accelerated qubits. arXiv:2207.13750

56. D. Boyanovsky, H.J. de Vega, R. Holman, M. Simonato, Dynamical renormalization group resummation of finite temperature infrared divergences. Phys. Rev. D 60, 065003 (1999). https://doi.org/10.1103/PhysRevD.60.065003. arXiv:hep-ph/9809346

57. C.P. Burgess, L. Leblond, R. Holman, S. Shandera, Super-hubble de Sitter fluctuations and the dynamical RG. JCAP 03, 033 (2010). https://doi.org/10.1088/1475-7516/2010/03/033. arXiv:0912.1608

58. D. Green, A. Premkumar, Dynamical RG and critical phenomena in de Sitter space. JHEP 04, 064 (2020). https://doi.org/10.1007/JHEP04(2020)064. arXiv:2001.05974

59. S. Brahma, A. Berera, J. Calderón-Figueroa, Universal signature of quantum entanglement across cosmological distances. arXiv:2107.06910

60. A.O. Caldeira, A.J. Leggett, Influence of dissipation on quantum tunneling in macroscopic systems. Phys. Rev. Lett. 46, 211 (1981). https://doi.org/10.1103/PhysRevLett.46.211

61. A.O. Caldeira, A.J. Leggett, Quantum tunneling in a dissipative system. Ann. Phys. 149, 374 (1983). https://doi.org/10.1016/0003-4916(83)90202-6

62. A.O. Caldeira, A.J. Leggett, Path integral approach to quantum Brownian motion. Phys. A 121, 587 (1983). https://doi.org/10.1016/0378-4371(83)90013-4
63. S. Choudhury, S. Panda, N. Pandey, A. Roy, Four-mode squeezed states in de Sitter space: a study with two field interacting quantum system. arXiv:2203.15815

64. T.S. Bunch, P.C.W. Davies, Quantum field theory in de sitter space: renormalization by point splitting. Proc. R. Soc. Lond. A360, 117 (1978). https://doi.org/10.1098/rspa.1978.0060

65. J. Grain, V. Vennin, Squeezing formalism and canonical transformations in cosmology. JCAP 2002, 022 (2020). https://doi.org/10.1088/1475-7516/2020/02/022. arXiv:1910.01916

66. B.L. Hu, A. Matacz, Quantum Brownian motion in a bath of parametric oscillators: a model for system—field interactions. Phys. Rev. D 49, 6612 (1994). https://doi.org/10.1103/PhysRevD.49.6612. arXiv:gr-qc/9312035

67. B.L. Hu, J.P. Paz, Y. Zhang, Quantum brownian motion in a general environment: exact master equation with nonlocal dissipation and colored noise. Phys. Rev. D 45, 2843 (1992). https://doi.org/10.1103/PhysRevD.45.2843

68. J.J. Halliwell, T. Yu, Alternative derivation of the hu-paz-zhang master equation of quantum Brownian motion. Phys. Rev. D 53, 12 (1996). https://doi.org/10.1103/PhysRevD.53.12

69. Y.-W. Huang, W.-M. Zhang, Exact master equation for quantum Brownian motion with generalization to momentum-dependent system-environment couplings. arXiv:2204.09965

70. L. Ferialdi, Exact closed master equation for gaussian non-Markovian dynamics. Phys. Rev. Lett. (2016). https://doi.org/10.1103/physrevlett.116.120402

71. L. Diósi, L. Ferialdi, General non-Markovian structure of gaussian master and stochastic schrödinger equations. Phys. Rev. Lett. (2014). https://doi.org/10.1103/physrevlett.113.200403

72. R.S. Whitney, Staying positive: going beyond Lindblad with perturbative master equations. J. Phys. A Math. Gen. 41, 175304 (2008). https://doi.org/10.1088/1751-8113/41/17/175304. [arXiv:0711.0074]

73. H.-P. Breuer, E.-M. Laine, J. Piilo, B. Vacchini, Colloquium: non-markovian dynamics in open quantum systems. Rev. Mod. Phys. 88, 021002 (2016). https://doi.org/10.1103/RevModPhys.88.021002

74. D. Moustos, C. Anastopoulos, Non-Markovian time evolution of an accelerated qubit. Phys. Rev. D 95, 025020 (2017). https://doi.org/10.1103/PhysRevD.95.025020. arXiv:1611.02477

75. F. Nicacio, R.N.P. Maia, Gauge quantum thermodynamics of time-local non-Markovian evolutions. arXiv:2204.02966

76. S. Prudhoe, S. Shandera, Classifying the non-Markovian, non-time-local, and entangling dynamics of an open quantum system. arXiv:2201.07080

77. G. Spaventa, P. Verrucchi, Nature and origin of the operators entering the master equation of an open quantum system. arXiv:2209.14209

78. D. Chruściński, Dynamical maps beyond Markovian regime. arXiv:2209.14902

79. C.A. Brasil, F.F. Fanchini, R.d.J. Napolitano, A simple derivation of the lindblad equation. Revista Brasileira de Ensino de Física 35, 01-09 (2013). https://doi.org/10.1590/s1806-11172013000100003

80. D. Manzano, A short introduction to the lindblad master equation. AIP Adv. 10, 025106 (2020). https://doi.org/10.1063/1.5115323

81. D. Baumann, D. Green, Equilateral non-Gaussianity and new physics on the horizon. JCAP 09, 014 (2011). https://doi.org/10.1088/1475-7516/2011/09/014. arXiv:1102.5343

82. S. Garcia-Saenz, S. Renaux-Petel, Flattened non-Gaussianities from the effective field theory of inflation with imaginary speed of sound. JCAP 11, 005 (2018). https://doi.org/10.1088/1475-7516/2018/11/005. arXiv:1805.12563

83. F.C. Lombardo, Influence functional approach to decoherence during inflation. Braz. J. Phys. 35, 391 (2005). https://doi.org/10.1590/S0103-97332005000300005. arXiv:gr-qc/0412069

84. M.G. Jackson, K. Schalm, Model independent signatures of new physics in the inflationary power spectrum. Phys. Rev. Lett. 108, 111301 (2012). https://doi.org/10.1103/PhysRevLett.108.111301. arXiv:1007.0185

85. M.G. Jackson, Integrating out heavy fields in inflation. arXiv:1203.3895

86. D. Boyanovsky, Information loss in effective field theory: entanglement and thermal entropies. Phys. Rev. D 97, 065008 (2018). https://doi.org/10.1103/PhysRevD.97.065008. arXiv:1801.06840

87. D. Boyanovsky, Imprint of entanglement entropy in the power spectrum of inflationary fluctuations. Phys. Rev. D 98, 023515 (2018). https://doi.org/10.1103/PhysRevD.98.023515. arXiv:1804.07967

88. C. Burrage, C. Käding, P. Millington, J. Minář, Open quantum dynamics induced by light scalar fields. Phys. Rev. D 100, 076003 (2019). https://doi.org/10.1103/PhysRevD.100.076003. arXiv:1812.08760

89. C. Burrage, C. Käding, P. Millington, J. Minář, Influence functionals, decoherence and conformally coupled scalars. J. Phys. Conf. Ser. 1275, 012041 (2019). https://doi.org/10.1088/1742-6596/1275/1/012041. arXiv:1902.09607

90. L. Pinol, S. Renaux-Petel, Y. Tada, A manifestly covariant theory of multifield stochastic inflation in phase space: solving the discretisation ambiguity in stochastic inflation. JCAP 04, 048 (2021). https://doi.org/10.1088/1475-7516/2021/04/048. arXiv:2008.07497

91. S. Choudhury, S. Dey, R.M. Gharat, S. Mandal, N. Pandey, Schwinger-Keldysh path integral formalism for a quenched quantum inverted oscillator. arXiv:2210.01134

92. C. Käding, M. Pitschmann, A new method for directly computing reduced density matrices. arXiv:2204.08829

93. H.-P. Breuer, A. Ma, F. Petruccione, Time-local master equations: influence functional and cumulant expansion. arXiv e-prints (2002) arXiv:quant-ph/0209153

94. D. Boyanovsky, Effective field theory out of equilibrium: Brownian quantum fields. New J. Phys. 17, 063017 (2015). https://doi.org/10.1088/1367-2630/17/6/063017. [arXiv:1503.00156]

95. A. Kamenev, A. Levchenko, Keldysh technique and non-linear sigma model: basic principles and applications. Adv. Phys. 58, 197-319 (2009). https://doi.org/10.1080/00018730902850504

96. W.B. Case, Wigner functions and weyl transforms for pedestrians. Am. J. Phys. 76, 937 (2008). https://doi.org/10.1119/1.2957889

97. C.P. Burgess, Introduction to effective field theory. Ann. Rev. Nucl. Part. Sci. 57, 329 (2007). https://doi.org/10.1146/annurev-nucl.56.080805.140508. arXiv:hep-th/0701053

98. C.P. Burgess, R. Holman, G. Kaplanek, Quantum hotspots: mean fields, open EFTs, nonlocality and decoherence near black holes. Fortsch. Phys. 2022, 2200019 (2021). https://doi.org/10.1002/prop.202200019. arXiv:2106.10804

99. S. Pi, M. Sasaki, Curvature perturbation spectrum in two-field inflation with a turning trajectory. JCAP 10, 051 (2012). https://doi.org/10.1088/1475-7516/2012/10/051. arXiv:1205.0161

100. R. Simon, E.C.G. Sudarshan, N. Mukunda, Gaussian-wigner distributions in quantum mechanics and optics. Phys. Rev. A 36, 3868 (1987). https://doi.org/10.1103/PhysRevA.36.3868

101. NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.1.6 of 2022-06-30