Nonperturbative semiclassical stability of de Sitter spacetime for small metric deviations

We consider the linearized semiclassical Einstein equations for small deviations around de Sitter spacetime including the vacuum polarization effects of conformal fields. Employing the method of order reduction, we find the exact solutions for general metric perturbations (of scalar, vector and tensor type). Our exact (nonperturbative) solutions show clearly that in this case de Sitter is stable with respect to small metric deviations and a late-time attractor. Furthermore, they also reveal a breakdown of perturbative solutions for a sufficiently long evolution inside the horizon. Our results are valid for any conformal theory, even self-interacting ones with arbitrarily strong coupling.


I. INTRODUCTION
De Sitter space can be understood as an exponentially expanding cosmological spacetime entirely driven by a positive cosmological constant. It plays a central role in most models of cosmological inflation, where the potential-dominated energy density and pressure of the inflaton field act approximately as a cosmological constant and lead to a quasi-exponential accelerated expansion. The inflationary scenario provides a natural mechanism, through parametric amplification of quantum vacuum fluctuations, for generating a nearly scale-invariant spectrum of adiabatic primordial inhomogeneities, which can successfully explain the observed CMB anisotropies and the large scale structure of the universe [1,2]. Furthermore, the physics of de Sitter could also be important for elucidating the final fate of the universe if its current accelerated expansion is entirely due to a small cosmological constant, a possibility compatible with observations so far.
The exponential expansion quickly redshifts away any initial perturbations. This offers a simple means of establishing natural initial conditions for subsequent evolution once such an accelerated expansion has already started over a region with a size larger than the Hubble radius and lasts for a sufficiently large number of e-foldings. Under these conditions the initial classical perturbations are effectively erased and the quantum state for modes with wavelengths much smaller than the Hubble radius is very close to the Bunch-Davies or Euclidean vacuum, which locally is essentially equivalent to the Minkowski vacuum * mfroeb@ffn.ub.edu † papadop@astro.auth.gr ‡ albert.roura@uni-ulm.de § enric.verdaguer@ub.edu since at those length-scales the spacetime appears almost flat. Such a scenario and the late-time attractor character of local de Sitter spacetime, often referred to as the "no-hair" property of de Sitter, is supported by a number of results and theorems in classical general relativity, both for linear perturbations [3,4] as well as for the full nonlinear case [5][6][7][8].
It is, therefore, of great interest to establish whether those classical no-hair results can be extended to the quantum mechanical case. Solid conclusions have recently been obtained within the framework of quantum field theory in curved spacetime [9,10]. Specifically, given a fairly general class of massive interacting theories with sufficiently weak coupling evolving on a fixed (nondynamical) de Sitter background, it has been shown to all orders in perturbation theory [11,12] that quantum correlators within a spacetime region of bounded physical size become at sufficiently late times arbitrarily close to those of the Euclidean vacuum (the generalization of the de Sitter-invariant Bunch-Davies vacuum to interacting theories).
Considering test fields evolving on a fixed background, however, offers an incomplete answer: addressing the full dynamical problem requires taking into account the backreaction of the quantum fields on the dynamics of the spacetime geometry. A number of studies have explored this question in the context of semiclassical gravity, where the metric is still treated classically, but its dynamics is governed by a generalization of the Einstein equation which includes the expectation value of the stress tensor operator of the quantum matter fields as a source [10,13]. Focusing on the backreaction of conformal matter fields, the dynamics of scalar-type metric perturbations around de Sitter has been analyzed in ref. [14], where the importance of considering also perturbations of the initial state of the matter fields has been emphasized. (The linear stability of de Sitter including initial classical stress tensor sources had earlier been studied in refs. [15,16].) In addition, the evolution of tensor metric perturbations for the same situation has been studied in ref. [17], where the semiclassical Einstein equation was solved perturbatively. (As we will show below, however, these perturbative solutions cease to be valid for a sufficiently long evolution inside the horizon.) The stability of tensor perturbations including the backreaction from conformal fields has also been considered [18] in investigations on the existence of unstable runaway solutions of the corresponding higher-order equations in the context of inflationary models driven by the trace anomaly [19]. Nevertheless, those analysis neglected the contribution of nonlocal terms which play a key role in the existence of runaway solution for perturbations around flat space [13,20].
Related studies have also been carried out for nonconformal fields in de Sitter spacetime. The case of massless minimally coupled free scalar fields has been considered in ref. [21] and a vanishing correction to the classical modes was found when solving perturbatively the linearized semiclassical equation for tensor perturbations. Both massless and massive nonconformal scalar fields were studied in refs. [22,23] and the stability of de Sitter spacetime with respect to spatially isotropic perturbations was established. A fairly general class of Gaussian initial states was considered and de Sitter was found to be a late-time attractor in all cases. Moreover, the importance of taking into account the contribution of nonlocal terms for light massive fields when analyzing the stability in the infrared regime was elucidated [23].
Here we consider the linearized semiclassical Einstein equation around a de Sitter background including the vacuum polarization effects of conformal matter fields, and solve it exactly for general metric perturbations (of scalar, vector and tensor type). In doing so, we make use of the method of order reduction [13,24], which eliminates the spurious solutions associated with higher-order derivatives while capturing the right dynamics in the infrared regime. Moreover, the method generates a backreaction equation which is equivalent to the original semiclassical Einstein equation up to the same order in inverse powers of the Planck mass at which the latter is valid within an effective field theory (EFT) approach to perturbative quantum gravity [25,26], but which can be significantly simpler to solve, a fact that we exploit in our calculation. Our exact solutions clearly show that de Sitter spacetime is also stable in this case and a latetime attractor as far as local geometrical properties are concerned. Furthermore, it reveals a breakdown of perturbation theory when solving the semiclassical equation for a long time evolution inside the horizon.
It should be stressed that our results are valid for any conformal field theory (CFT), even self-interacting ones with arbitrary strong coupling, as explained in sec. IX. In addition to metric perturbations, perturbed initial states of the matter fields have also been considered.
Although semiclassical gravity does take into account the backreaction of the quantum matter fields on the dynamics of the mean spacetime geometry, it does not provide a complete analysis because it does not include the quantum mechanical effects of the metric itself. Indeed, one needs to quantize the metric perturbations in order to account for certain relevant phenomena: doing so is necessary, for instance, for a proper description of the generation of primordial cosmological perturbations, and it has even been suggested that radiative corrections involving higher-order graviton loops could lead to a secular screening of the cosmological constant [27,28]. However, detailed calculations including graviton loops are technically complex and one is, in addition, confronted with the need to consider appropriate observables which are not only gauge-invariant beyond linear order but also infrared safe [29][30][31]. Because of such difficulties only partial progress has been made in this direction. It is, therefore, important to consider also somewhat less ambitious problems, but obtain solid results (and, if possible, exact) which can provide a robust foundation for further developments. One such example is the exact calculation of one-loop corrections from matter fields to the correlator of the Riemann tensor [32][33][34] for quantized metric perturbations around de Sitter. The results support the existence of quantum states for metric perturbations interacting with matter fields which exhibit (appropriately defined) de Sitter invariance, at least when graviton loops are neglected. In contrast, the results presented here only apply to the mean field geometry, but explore the effect of different (non-de Sitter-invariant) initial states on the dynamics and show not only that a self-consistent de Sitter-invariant solution exists, but also that it is a late-time attractor. Furthermore, the methods described below for obtaining nonperturbative solutions valid for long evolution times could prove helpful in order to extend the calculation of the Riemann correlator so that it correctly captures the details of its behavior for large separations (both spatial and temporal), which seems to require a nonperturbative treatment.
The rest of the paper is organized as follows. Semiclassical gravity and the method of order reduction are briefly reviewed in sec. II. The semiclassical Einstein equation for metric perturbations around a spatially flat FLRW spacetime including the quantum back-reaction of conformal fields is presented in sec. III. Given a cosmological constant and fields in the Bunch-Davies vacuum, there is a self-consistent semiclassical de Sitter background. In sec. IV the metric perturbations are decomposed into scalar, vector and tensor contributions. Next, we fix the gauge and write the (decoupled) semiclassical equations for the three types of perturbations, before and after employing the order reduction method. The exact (nonperturbative) solutions are obtained in sec. V and their implications for the stability of de Sitter spacetime as well as the breakdown of the perturbative solutions are analyzed in sec. VI. The effects of perturbing also the initial state of the matter fields are studied in sec. VII and we show that all our main conclusions remain unchanged.
Finally, in sec. IX we summarize and discuss our results as well as explaining their applicability to any CFT for the matter fields. Several useful formulae concerning the perturbative expansion of curvature tensors, their conformal transformations and some special functions are provided in the first three appendices. In addition, a method for generating a family of regular Gaussian initial states is described in appendix D, and in appendix E we compare our linearized semiclassical Einstein equation for tensor perturbations (before order reduction) with the one previously obtained by Starobinsky [35].
Throughout the paper we use natural units with c = = 1 and take κ 2 = 16πG N . We employ the "+++" sign convention of ref. [36] and Greek indices range over space and time, while Latin indices denote spatial components only.

A. Semiclassical Einstein equations
Semiclassical gravity can be regarded as a mean field approximation to a quantum theory of gravity where the mean gravitational field is treated classically and only matter is quantized. In contrast to quantum field theory in curved spacetime, it also includes the back-reaction of the matter fields on the mean geometry, given by the background metric g µν . To achieve this, the stress tensor on the right-hand side of the Einstein equation is replaced by the expectation value of an appropriate quantum stress tensor operator. This expectation value needs to be renormalized, and counterterms local in the gravitational field have to be included in the bare gravitational action. At all loop order for the matter fields (but no graviton loops) those counterterms are quadratic in the curvature and the renormalized semiclassical Einstein equation reads and B µν = (−g) −1/2 (δ/δg µν ) R 2 √ −g d 4 x are obtained by functionally differentiating the finite parts of the gravitational counterterms. The parameters a 1 (µ) and a 2 (µ) together with Λ(µ) and κ 2 (µ) are in general renormalized parameters which have to be determined by experiment, and µ is the renormalization scale. Note, nevertheless, that the backreaction equation (1) is renormalization group invariant and the dependence on µ of the different parameters appearing in the equation and the renormalized expectation value T µν (µ) ren cancel out.
The tensors A µν and B µν are explicitly given by where the two forms of A µν given above were obtained by using the identity (A5) and the definition of the Weyl tensor (A4) as well as the second Bianchi identity.

B. Order reduction
The semiclassical Einstein equation (1) contains terms with up to fourth-order derivatives of the metric, as seen from eq. (2). (The expectation value T µν (µ) ren also involves similar terms, as shown below.) Such kind of higher-order time derivatives are common in backreaction problems. A well known example is the Abraham-Lorentz-Dirac equation, which describes the effect of radiation reaction on the motion of a point-like charge in classical electrodynamics [37,38] (i.e. without considering the internal structure of the particle nor a finite size for the charge density distribution). In fact, they are a generic feature of effective field theories (EFTs), where the effects of the UV sector on the dynamics of the lowenergy degrees of freedom are encoded at the level of the action through an expansion of local terms with an increasing number of derivatives. The validity of the EFT expansion relies on the fact that for length-scales much larger than the inverse cut-off scale of the UV sector the higher-order terms in the expansion become increasingly smaller. In this regime their contribution amounts to a small correction to the equation of motion which results, when treated perturbatively, into locally small perturbations of the classical solutions. In contrast, solving the corresponding higher-order equations exactly gives rise to additional solutions exhibiting exponential instabilities with characteristic time-scales comparable to the inverse cutoff scale of the EFT (or sometimes fast oscillations with the same kind of characteristic timescale), often referred to as "runaway" solutions. These are spurious solutions which should not be taken seriously since they involve characteristic scales for which the EFT expansion breaks down and the contributions from the higher-order terms to the equation of motion no longer correspond to small corrections but to dominant terms.
The simplest way of avoiding such spurious solutions is by solving the corrected equations of motion perturbatively. However, perturbative solutions may not be valid for long times. This happens when quantities like the total time appear multiplying the perturbative parameter so that the expansion contains so-called secular terms which grow with time and lead to a breakdown for sufficiently long times of the truncated perturbative expansion. Those limitations can be overcome with the order reduction method, which consists in taking the equation of motion with corrections up to a finite order and writing an alternative equation which is equivalent up to that order but contains no higher derivative terms (this is achieved by taking successive derivatives of the original equation and substituting the higher-order derivatives in the correction terms to the appropriate order). The exact solutions of the equation obtained with this method agree locally with the perturbative solutions constructed around different times (each one with a finite domain of validity) and provides an interpolation between all of them valid for long times. This is particularly important when considering situations where the effects of the corrections are locally small, but can build up over long times and give rise to substantial accumulated effects. Two examples of such situations are an electric charge following a quasi-circular trajectory in a uniform magnetic field and emitting electromagnetic radiation for a sufficiently long time so that the radius of its orbit decreases, say, to half of its initial value due to radiation reaction, or an evaporating black hole emitting Hawking radiation for such a long time that its mass (or horizon size) decreases to a small fraction of its initial value.
Related alternative methods which have been employed in the literature for discarding the spurious solutions mentioned above involve finding the exact solutions of the original backreaction equation and then selecting the appropriate subset either by demanding analyticity of the solutions with respect to the perturbative parameter or checking explicitly which solutions exhibit unphysical characteristic scales and disregarding them. However, the latter method is less systematic and requires a case-by-case analysis, whereas the analyticity requirement may be too restrictive in some cases [13]. Furthermore, the order reduction method leads to equations of motion which are equivalent up to the order under consideration, but are often easier to solve, as will be the case for the problem analyzed in the remaining sections.
The order reduction method has been applied to electromagnetic [38] and gravitational [39] radiation reaction problems as well as higher derivative gravity [40]. It has also been employed in semiclassical gravity [13,24] and in this context it has been argued [41] that traceanomaly-driven inflationary models (with no cosmological constant and driven entirely by the vacuum polarization of large number of matter fields [19]) correspond to spurious solutions which lie beyond the EFT's domain of applicability and are automatically discarded when using order reduction.
It should be noted that order reduction cannot be always applied in a straightforward way. It may be ambiguous in integro-differential equations, or may lead to covariance breaking if the time derivatives and spatial derivatives are not simultaneously reduced; see [13] for a detailed discussion of these issues.
The order reduction method can be illustrated in a nutshell with the following simple example of a first order differential equation in time for a function f (η) with a perturbative correction of order κ 2 . Given where b is a constant and P is an arbitrary function, order reduction uses that f ′ = −bf + O κ 2 , and by deriving one more time Substituting those two equations into the right hand side, we get which is an equation of first order which is valid to the same order in κ 2 as the original equation (3), but does not have unphysical solutions. Rather than considering a truncated perturbative expansion, this equation can now be solved exactly. It is clear how the method works for equations of more derivatives or partial differential equations: one takes the lowest order equation and substitutes it in the higher order terms (in κ 2 ), taking additional derivatives if necessary.

III. CONFORMAL FIELDS IN A PERTURBED FLRW UNIVERSE
A. The model Our model consists of N massless free scalar fields conformally coupled to the spacetime curvature: whereR is the Ricci scalar associated with the metricg µν and ξ cc = (d − 2)/4(d − 1), which reduces to ξ cc = 1/6 in four dimensions. We will use dimensional regularization, but after the renormalization procedure has been carried out we will take d = 4. Furthermore, we will specialize the physical metricg µν to a slightly perturbed spatially flat FLRW spacetime, which is conformal to an almost flat metric g µν : where η denotes the conformal time.
The renormalized semiclassical Einstein equation (1) for this model with the fields in the conformal vacuum state was derived by Campos and Verdaguer in refs. [42,43] using the closed-time-path (CTP) effective action, and is given to linear order in the perturbation h µν by where α = N/(2880π 2 ). Here and throughout the rest of the paper ∇ µ denotes the covariant derivative with respect to the metric g µν and g = ∇ µ ∇ µ , whereas objects with a tilde are evaluated with the conformally related metricg µν . Furthermore, background quantities will be denoted by a superscript (0) such as g (0) µν = a 2 (η)η µν , and quantities linearized in the perturbation by a superscript (1) as ing For conformal fields the renormalized parameter a 2 (µ) does not depend on µ and we denote it by β. Moreover we have chosen a renormalization scaleμ such that a 1 (μ) = 0. Similarly, both Λ and κ 2 are also independent of µ in this case. The tensors A µν and B µν are defined in equation (2) and Finally, the kernel H(x − y; µ) depending on the renormalization scale µ is given in appendix E, but its exact form will not be needed in the bulk of the paper. Eq. (7) coincides with those derived by alternative methods [35,44,45].

B. The FLRW background
The semiclassical generalization of the Friedmann equation can be obtained by setting the perturbation h µν to zero in the 00 component of eq. (7), which gives and using the order reduction method it becomes Defining an effective cosmological constant Λ eff as Eq. (10) has the solution where the Hubble parameter H is given by 3H 2 = Λ eff , and −∞ < η ≤ 0. This solution is unique up to a shift of the origin of conformal time, η → η − η 0 , and its sign.
Hence, de Sitter spacetime, given here in spatially flat coordinates (the Poincaré patch), is a self-consistent solution of the semiclassical Friedmann equation (9), with the effective cosmological constant (11) having a small positive shift of quantum origin. The existence of such selfconsistent solutions follows straightforwardly from the fact that the renormalized expectation value of the stress tensor for the Bunch-Davies vacuum must be proportional to the metric, as implied by de Sitter invariance, and has been know for a long time [46,47]. When Λ = 0, eq. (9) still admits a de Sitter solution with H 2 = 2/(ακ 2 ) (closely connected to Starobinsky's original model of inflation [19,48]), but its characteristic scale lies beyond the domain of validity of semiclassical gravity when regarded as part of an EFT approach to quantum gravity, as briefly discussed in sec. II B. Such solutions are automatically discarded by the method of order reduction [24].

A. Gauge fixing
When considering metric perturbations around a given background, there is a gauge freedom (corresponding to local diffeomorphisms) associated with the mapping between the background and the perturbed geometry. Infinitessimal diffeomorphisms generated by an arbitrary vector fieldξ ν induce the following gauge transformation of the physical perturbation a 2 h µν introduced in Eq. (6): Rescaling the arbitrary vector fieldξ µ = a 2 ξ µ and using Eq. (12) for the scale factor a(η), Eq. (13) becomes Before proceeding any further, it is convenient to decompose the perturbation h µν exploiting the fact that the spatial sections of the background metric are maximally symmetric spaces [49][50][51]. The spatial part h ij decomposes as where δ ij ∂ i h TT jk = 0 = δ ij h TT ij and δ ij ∂ i w T j = 0. The temporal components can be similarly decomposed as where δ ij ∂ i v T j = 0. In total, we have four scalars φ, ψ, σ, τ , two transverse vectors w T i and v T i (with two independent components each) and a transverse traceless tensor h TT ij (with two independent components as well). Decomposing also the spatial part of the vector field ξ µ as where δ ij ∂ i ξ T j = 0, we can see the behavior of the various components under a gauge transformation: where primes denotes derivatives with respect to the conformal time η. Choosing ξ T i , ξ and ξ 0 appropriately, we can set w T i , σ and τ to zero. We will work in this gauge, where the perturbation of the spatial metric is entirely given by the tensorial component: This is the transverse traceless gauge, also known as spatially flat gauge when focusing on the scalar perturbations. This fixes completely the gauge if we restrict ourselves to metric perturbations that fall off at spatial infinity.
If we consider perturbations which do not necessarily fall off, there is still some residual gauge freedom which is not fixed by condition (19). On one hand, there are transformations which leave h ij invariant. One possibility are those generated by ξ T i which are functions only of the conformal time; this corresponds to translations on the spatial sections which can change for each surface of the foliation and lead to changes of v T i . A second possibility involves transformations generated by ξ i = b(η)δ ij x j , ξ 0 = −b(η)η; this corresponds to dilations on the spatial sections while changing at the same time the surface of the foliation so that the expansion of the FLRW background compensates for that, and leads to changes of φ and ψ while leaving h ij invariant. On the other hand, the transformations generated by ξ T i = E ij x j with E ij constant and traceless induce changes of h ij but leave it transverse and traceless. These residual gauge transformations will play a role in sec. V to show that certain solutions are pure gauge.

B. Semiclassical equations of motion
Using the metric decomposition and the gauge fixing introduced in the previous subsection, one can obtain from the semiclassical equation (7) the dynamical equations for the scalar, vector and tensor perturbations, which decouple form each other. From the 00 component of eq. (7) and the scalar part of its 0i component, one gets the following two equations for the scalar perturbations ψ and φ: Similarly, from the transverse part of the 0i component one gets the equation for the vector perturbation v T i : Finally, the equation for the tensor perturbations is obtained from the transverse and traceless part of the ij components: Employing order reduction as explained in sec. II B, the equations can be rewritten in the much simpler form where we have introduced the following parameter, which controls the expansion in powers of κ 2 : It is worth emphasizing that those equations are independent of the arbitrary parameter β and the renormalization scaleμ of the semiclassical theory. They involve only the semiclassical parameter α, which depends on the matter field content.

C. Nonlocal terms
As we have calculated explicitly, when we use order redution the nonlocal terms do not contribute to the semiclassical equations of motion (23) for the perturbation h µν . We now want to show that one can see this in general, without choosing a gauge or expanding the semiclassical equations explicitly in terms of the perturbation h µν . Basically this amounts to showing that A µν , given by eq. (2), is of order κ 2 when using order reduction. In order to do so, it is convenient to consider its definition as a functional derivative of the integral of the square of the Weyl tensor with where the last equality follows from the conformal invariance of the Weyl tensor with one raised index. This also implies that and we can equivalently show thatÃ µν is of order κ 2 .
Another consequence of the invariance of F under conformal transformations ofg µν is the vanishing trace of A µν : so thatg µνÃ µν = 0. Expressing now the Weyl tensor in terms of the Riemann tensor and its contractions according to its defining equation (A4), we can write where E 4 is the integrand of the four-dimensional Euler invariant and is given by The generalized Gauß-Bonnet theorem establishes that the integral of E 4 is a topological invariant, namely 32π 2 times the Euler characteristic. From eq. (29) and the fact that the variational derivative of a topological invariant vanishes, we see thatÃ µν can be expressed entirely in terms of the Ricci tensor, the Ricci scalar and covariant derivatives acting on them.
We are now ready to use order reduction. Taking into account thatR and substitutingR µν = Λg µν into the expression forÃ µν , the result can only be proportional tog µν up to order κ 2 . However, sinceg µνÃ µν = 0, we conclude thatÃ µν is of order κ 2 when order reduction is employed. One can alternatively check this fact by substituting eq. (31) into the explicit expression forÃ µν in terms of the Ricci tensor in eq. (2).
In conclusion, we see that we only have to consider the local terms in eq. (7) when using order reduction.

A. Scalar and vector perturbations
The solutions of eq. (23b) for the components of the vector perturbation v T i are arbitrary functions of time, which can be eliminated by a gauge transformation. Indeed, by using the residual gauge freedom described at the end of sec. IV A and choosing ξ T i as an appropriate function of time only, we can set v T i = 0. For φ, the solution of eq. (23c) is also an arbitrary function of time. The solution for ψ is then given by Since ψ enters into the perturbation h µν only through a spatial derivative according to eq. (16), the arbitrary function f (t) does not change the perturbation h µν and we can set it to zero. If we want to start with bounded initial perturbations, we must exclude solutions which are unbounded and have to take φ = 0.
On the other hand, if we had not excluded such unbounded solutions, we would first have to choose ξ 0 appropriately to make φ vanish, and would need to take a similar unbounded function ξ = −1/(2η)ξ 0 r 2 so that the combination ∂ i ∂ j σ + τ δ ij which enters into the decomposition of the perturbation (15) still vanishes. The solution for ψ is then an arbitrary function of time which we can set to zero as above.
Thus, we see that when order reduction is employed, both vector and scalar parts after solving the constraints are pure gauge and can be eliminated by a residual gauge transformation of the kind mentioned at the end of sec. IV A.
Note that Anderson et al. [14], who investigated scalar perturbations without using order reduction, also concluded that those perturbations (which they refer to as perturbations of the first kind) have to vanish because the corresponding solutions that they found lie outside the range of validity of the semiclassical theory. The gauge-invariant variables that they introduce, however, take a simple form in a gauge rather different from the one we employ, and so their intermediate expressions are not directly comparable to ours.

B. Tensor perturbations
To find the solutions for the tensor perturbations, we take the Fourier transform with respect to the spatial coordinates, where e ± ij (p) are a pair of transverse and traceless tensors corresponding to two different polarizations. Equation (23a) then becomes Setting ω 2 = (1 − 2ν)p 2 , s = −ωη and g ± = s 3 2 −ν f ± (s), this reduces to a Bessel equation for f ± , whose general solution is where C ± 1 and C ± 2 are integration constants. For the particular case p = 0, which corresponds to no spatial dependence in position space, eq. (23) can be solved directly and the general solution is given by where D ij and E ij are traceless tensors with respect to the induced background metric of the spatially flat sections, and independent of the spatial coordinates and the conformal time. 1 The first term on the right-hand side of eq. (36) corresponds to a Bianchi I anisotropic deformation of de Sitter, whereas the second one is pure gauge and can be eliminated by the residual gauge transformation generated by the transverse vector ξ T i = E ij x j /2.

VI. DE SITTER STABILITY AND SECULAR TERMS
A. Stability with respect to linear perturbations In this section we analyze the stability of the semiclassical de Sitter geometry with respect to small metric perturbations including the back-reaction due to quantum vacuum effects from conformal fields. We do so by focusing on the evolution of the Riemann tensor associated with the linearly perturbed metricg µν , which has a number of appealing properties. First of all, the linear perturbation around de Sitter of the Riemann tensorR αβ γδ with appropriately raised indices is a gauge-invariant object. This follows from the fact that for the unperturbed background it can be written asR (0) δ] , whose Lie derivative with respect to an arbitrary vector field vanishes. Furthermore, with this index structure the components remain unchanged when rescaling by the same constant the basis vectors of the tangent space at a given point. This implies that the components coincide with those in the physical basis of orthonormal vectors {a −1 ∂ 0 , a −1 ∂ i } of the background metric. Finally, the Riemann tensor provides a suitable characterization of the local geometry, in terms of which the stability and attractor nature of semiclassical de Sitter spacetime can be naturally formulated, as further discussed at the end of this subsection.
In terms of the metric perturbations the linearized Riemann tensor is given bỹ Using the gauge transformation (14), one can explicitly check that it is indeed gauge invariant at linear order. Fourier transforming with respect to the spatial coordinates and specializing to tensor perturbations, we get 1 These solutions can also be obtained by taking the limit ω → 0 of eq. (35) after rewriting C 1 =C 1 /ω 3−2ν so that a finite non-vanishing limit is obtained for the solutions associated with the two integration constants. In addition, the different ways of taking the limit p → 0 of e ± ij (p) and the possibility of considering arbitrary linear combinations gives rise to the general traceless tensors D ij and E ij . where Hence, we can see that all the Riemann components can be written in terms of g ± and g ′ ± . Since everything that will be said is entirely equivalent for both polarizations, for ease of notation we will omit in the remainder of this section the subindices ± labeling the two transverse polarizations associated with each momentum p.
Let us consider first the evolution of the Riemann tensor for modes well outside the horizon, i.e. with |ωη| ≪ 1. In this case one needs to evaluate the Bessel functions in eq. (35) using eqs. (C3) and (C5), which leads to Substituting into eq. (39) we see that the components of the Riemann perturbation in a physical basis decay like 1/a = −Hη or higher order at late times, i.e. in the limit η → 0. On the other hand, for modes inside the horizon, with |ωη| ≫ 1, one can use eqs. (C4) and (C6) to see that g and g ′ are of the form times an oscillatory factor corresponding to a linear combination of sin(ωη) and cos(ωη). Thus, from eq. (39) it follows that inside the horizon the components of the Riemann perturbation oscillate with an amplitude that decays like 1/a 1−ν .
Putting these results together we can conclude that de Sitter spacetime remains stable with respect to small metric perturbations of the semiclassical mean geometry when the quantum back-reaction of conformal matter fields is included. This is guaranteed by the fact that for any Fourier mode with comoving momentum p the perturbation of the Riemann tensor decays like 1/a 1−ν (times an oscillatory factor) when the corresponding physical wavelength 2πa/|p| is smaller than the de Sitter radius 1/H and like 1/a when it is larger, together with the regularity of the perturbation around horizon crossing (when the wavelength is comparable to 1/H). This extends the conclusions of the no-hair theorem for de Sitter spacetime, which is not only stable with respect to small metric perturbations but also a late-time attractor in classical general relativity, to the case where radiative corrections from loops of conformal fields are considered. In fact, the main effect of the radiative corrections compared to the classical case for pure gravity, which corresponds to taking ν = 0 in our results, is simply to alter slightly the exponent of the power-law decay for modes inside the horizon.
Our result can be used to illustrate in a simple way the fact that the stability and the character of late-time attractor of de Sitter spacetime applies to sufficiently localized observables characterizing the geometry within a region of fixed physical size (as opposed to comoving). The tensor perturbation h TT ij (η, x) and the amplitude g(η, p) associated with a given momentum and polarization are gauge-invariant objects with well-defined geometrical meaning. However, as their characteristic physical wavelength gets exponentially redshifted, at late times one would need to measure them over regions with a physical size that becomes arbitrarily large. Instead, the deviations of the geometric properties within a region of fixed physical size compared to those of de Sitter decay exponentially with the (proper) cosmological time. These features are adequately captured by the behavior of the Riemann tensor, which provides a good characterization of the local geometry.

B. Perturbative vs. nonperturbative solutions
In this subsection we will compare the exact nonperturbative solutions of the linearized semiclassical equation (34), given by eq. (35), to those that result from solving the equation perturbatively in ν. We will see that the perturbative solutions cease to be valid when the modes evolve inside the horizon for a sufficiently long time.
Here we concentrate on initial perturbations corresponding to Bunch-Davies positive frequency modes, but our conclusions can be easily generalized to arbitrary initial conditions. This choice amounts to setting C 2 = iC 1 , so that the linear combination within the square bracket on the right-hand side of eq. (35) becomes a Hankel function of the first kind, with a purely positive frequency oscillatory behavior at early times. Indeed, using eqs. (C4) and (C6) one can see that in this case the behavior for |ωη| ≫ 1 of the non-perturbative solution is given by On the other hand, if we first expand the solution (35) in powers of ν employing eqs. (C8)-(C9), we get At late times (for |ωη| ≪ 1) the exact solution is well approximated by the perturbative solution (43): one recovers the asymptotic behavior in eq. (40) but with the constant term and the coefficients of the higher-order ones given by their expansion in powers of ν truncated at linear order. That is, however, not the case at early times, with |ωη| ≫ 1. This can be seen by considering the early-time limit of the perturbative solution (43): which coincides with the result that one obtains by expanding in powers of ν the asymptotic expression (42) for the exact solution. It is clear that the perturbative solution deviates significantly from the exact one when ν|p||η| 1. The implications can be more easily understood if we normalize the modes so that they have a fixed amplitude at the initial time η 0 independently of the particular value of η 0 , which can be implemented by dividing the mode by its value at η 0 . Proceeding in this way with the exact solution, one finds for example that the amplitude of a mode which was initially well within the horizon has decreased by a factor 1/(−ωη 0 ) 1−ν by the time of horizon crossing. In contrast, repeating the procedure with the perturbative solution and treating ν perturbatively when normalizing by the value at η 0 , one obtains an amplitude at horizon crossing of order (−1/|p|η 0 ) (1 + ν ln(−|p|η 0 ) − iν|p|η 0 ). For modes which have spent a long time inside the horizon, so that ν|p||η 0 | 1, this amplitude at horizon crossing can be significantly larger.
The reason for the potentially large deviation of the perturbative solution with respect to the exact one is a breakdown of perturbation theory: the actual condition for the validity of the perturbative solution is ν|p|η 0 ≪ 1, which can be violated even for ν ≪ 1 when considering |p||η 0 | large enough. The term proportional to iν|p|η in eq. (44), responsible for the main deviations, is a secular term arising from the truncated perturbative expansion in powers of ν of the oscillatory factor e −iωη in eq. (42). Such a breakdown of perturbation theory for long times is very common when determining the evolution of a system by solving perturbatively the corresponding dynamical equations. This can be illustrated with the simple example of a harmonic oscillator with frequency Ω + δΩ. If one solves perturbatively in δΩ the corresponding equation of motion, x+(Ω+δΩ) 2 x = 0, to first order one finds a solution of the form x(t) ≈ A sin(Ωt + ϕ) [1 + δΩt], where the correction grows with time. In this case, however, the exact solution is obviously known: x(t) = A sin (Ω + δΩ)t + ϕ . It is, therefore, clear that the exact solution is qualitatively very similar to the unperturbed one but with a slightly corrected frequency. The growing terms in the perturbative solution, which are commonly know as secular terms, reflect the fact that the perturbative solution is not valid for arbitrarily long times but restricted instead to times such that δΩ t ≪ 1. This is a rather simple example, but the situation is completely analogous for our semiclassical solution.
We close this subsection by comparing the perturbative solution obtained by expanding the exact solution, which has been discussed above, with the one obtained by solving perturbatively eq. (34) as done in ref. [17]. Substituting the classical solution for g into the terms proportional to ν in eq. (34), treating them as a source, and using the retarded propagator for the unperturbed equation (with ν = 0), one gets the following result for the perturbative correction, already obtained in ref. [17] (see their eq. (48)): The integral can be done exactly with the help of appendix C. For large negative values of η 0 one gets where we have set η = 0 for the upper limit, whose ef-fect is unimportant for a late-time solution. (Note that the logarithmic term and the phase were not included in ref. [17], where only the dominant term was retained.) If we had normalized the classical solution at the initial time, which essentially amounts to dividing the unperturbed solution by 1/η 0 , the perturbative correction in eq. (46) would agree with the perturbative solution (44) when normalized at the initial time η 0 as discussed below that equation. There is actually a slight discrepancy: one does not get the term e −2i|p|η0 /2 appearing in eq. (46). This is due to a different choice of initial conditions. Whereas our exact solution (42) results from choosing C 2 = iC 1 exactly (at all orders in perturbation theory), eqs. (45)- (46) correspond to making the same choice for the zeroth-order solution but requiring the first-order correction and its derivative to vanish at the initial time. By imposing the same conditions when determining our exact and perturbative solutions, one obtains a new pair of constants C 1 and C 2 which differ at order ν and lead to a perturbative solution in full agreement with eq. (46).
It should be pointed out that the perturbative correction can be large only when the corresponding mode was trans-Planckian at the initial time, i.e. its physical wavelength was much smaller than the Planck length (otherwise one has −ν|p|η 0 ≪ 1 and the perturbative correction is always small), as recognized in ref. [17]. At such scales semiclassical gravity is not guaranteed to provide an accurate description. Furthermore, one would need a rather small amplitude of the initial perturbations so that nonlinear gravitational effects do not become important: otherwise for such short wavelengths the effective stress tensor quadratic in the metric perturbations [52,53] can generate a strong back-reaction on the background expansion, and even dominate over the cosmological constant. In any case, even if one carries out a linearized analysis without much concern for these issues, as done in ref. [17], the corrections to the classical solution will always be small as shown by our exact solutions and discussed above.

VII. PERTURBED INITIAL STATE
So far we have considered the dynamics of metric perturbations, but for a fixed initial state of the matter field, namely the Bunch-Davies vacuum. The state-dependent expectation value of the stress tensor is affected by the metric perturbations, but this is due to the their effect on the evolution of the scalar field operator in the Heisenberg picture or, alternatively, on the evolution of the state in the Schrödinger picture. However, if one wants to allow changes in the initial state, eq. (7) needs to be generalized. In fact, as discussed in appendix D, the Bunch-Davies vacuum of the matter fields is no longer a Hadamard state (free of excitations at arbitrarily short wavelengths) for nonvanishing metric perturbations at the initial time: Therefore, the initial state needs to be modified so that it is a Hadamard state in that case and, in particular, the renormalized expectation value of the stress tensor is finite at the initial time.
Starting with eq. (1) and considering not only linear perturbations around the de Sitter metric but also small perturbations of the initial state, one obtains where the the right-hand side corresponds to the perturbation of the expectation value of the stress tensor evaluated on the background metric due to the perturbation of the initial state. Here it has been assumed that such a term is of the same order as the remaining terms, which are linear in the metric perturbations. Hence, no further terms should be considered since those corresponding to the perturbation of the initial state in the stress tensor expectation value evaluated on the perturbed metric would be of higher order. Note that eq. (47) and the procedure employed below can also be used to consider the effect of any other additional stress tensor sources (even classical ones) which can be treated perturbatively and regarded of the same order as the terms linear in the metric perturbations.
For conformal fields in a FLRW background eq. (47) for the linear metric perturbations reduces again to eq. (7) plus the source term involving δT µν . Furthermore, one can introduce a decomposition of δT µν analogous to that for the metric perturbations in eqs. (15)- (16), which is complemented now by the conservation requirement for δT µν with respect to the background metric, i.e.
(0) ∇ µ δT µν = 0 with (0) ∇ µ being the covariant derivative associated with the background metric g (0) µν = a 2 η µν . More specifically, we can write and take as independent quantities the transverse and traceless tensor δT TT ij , the transverse vector δT T 0i and the two scalar functions δT 00 and δT . The temporal component of the conservation equation determines the scalar χ and, similarly, the scalar part of the spatial projection of the conservation equation, which can be written as the gradient of a function, determines the scalar ρ. On the other hand, the transverse part of this spatial projection determines the transverse vector W T j . Applying this decomposition to the generalization of eq. (7) including the source δT µν , one obtains eqs. (20)- (22) plus the corresponding sources. Finally, using order reduction eqs. (23) are generalized to which correspond, respectively, to the transverse and traceless part of the ij component of the order-reduced equations, the transverse part of the 0i component, the spatial divergence of the 0i component, and the 00 component. Moreover, the conservation equation has been used to express the spatial divergence of δT 0i in terms of δT 00 and δT on the right-hand side of eq. (49c).

A. Tensor perturbations
Let us start with eq. (49a) for the tensor perturbations. In addition to employing eq. (33) for h TT ij , one can use the analogous expression for the stress tensor perturbation and obtains the following equation for each one of the two polarizations: This equation is identical to eq. (34) but with an inhomogeneous source term. The general solution can be written as g ± = g ± h + g ± i , a sum of a homogeneous solution containing the information on the initial conditions at the time η 0 and a particular solution of the inhomogeneous equation with vanishing initial conditions, where G ret is the retarded propagator associated with eq. (51). Given two independent homogeneous solutions, u 1 and u 2 , the retarded propagator for such a linear second-order differential equation can be expressed as where is the Wronskian for this pair of solutions, which is nonzero for independent solutions. Therefore, the inhomogeneous solution can be written as with Provided that C ± j (η) are regular and have a finite limit when η → 0, as we will discuss next, the late-time behavior of the contribution form the inhomogeneous solution g ± i is the same as for the homogeneous solution, already analyzed in sec. VI A. We can consider the same pair of homogenous solutions as in that section and choose Taking into account that for this choice W (η ′ ) = (2ω/π)(−ωη ′ ) 2−2ν , we can immediately see that C ± j (η) will be regular and have a finite limit when η → 0 provided that J ± (η ′ ) is also regular and decays faster than (−η ′ ) 1−2ν as that limit is approached. This means that the components of the stress tensor δT TT ij should decay slightly faster than 1/a in conformal coordinates, or 1/a 3 in a physical basis. Such kind of behavior is fulfilled by classical radiation, which decays like 1/a 4 , and it seems plausible that small excitations of the vacuum state for the conformal fields considered here also exhibit the same decay; it is indeed the case for the class of regular states considered in appendix D, as shown there. We stress again that δT µν has been assumed to be small enough so that it is at most of the same order as the terms linear in the metric perturbation in eq. (47) and nonlinear gravitational effects are not important (otherwise one could consider for instance an excitation with sufficiently large energy overdensity to form a black hole by gravitational collapse, contradicting the conclusions above).
As seen above, if the quantities in eq. (55) have a regular limit as η → 0, the inhomogeneous solutions g ± i behave in the same way as the homogeneous ones and the same conclusions drawn in sec. VI apply to the inhomogeneous solution as well.

B. Scalar and vector perturbations
Unlike tensor perturbations, vector and scalar ones do not give any nontrivial contribution in the classical case, but they get a nonvanishing semiclassical correction. Nevertheless, since they are governed by dynamical constraints which become simple algebraic equations when working in Fourier space for the spatial coordinates, their long time behavior can be directly determined from the fall-off properties of the stress tensor perturbation δT µν .
Let us consider the vector perturbations first. From eq. (49b) it follows that in spatial Fourier spaceṽ T i = (κ 2 /p 2 ) δT T 0i and the metric perturbations decay in the same way as the stress tensor. From eq. (37) it follows that their contribution to the linear perturbation of the Riemann tensor is given bỹ and we can conclude that at late-times the curvature perturbations decay like 1/a δT T 0i and, hence, for any decaying (or even asymptotically constant) stress tensor perturbation the stability and attractor character of de Sitter spacetime are not altered.
The situation is similar for scalar perturbations but with some slight differences. From eqs. (49c)-(49d) one can see that φ decays at least like δT 00 or η µν δT µν , whereas ψ is only guaranteed to behave at late times like a δT 00 or a η µν δT µν . Given the contributions of scalar perturbations to the linearized Riemann tensor, one can see that they basically fall off like δT µν at late times. Therefore, if δT µν decays at least like 1/a the conclusions about the stability of de Sitter space remain unchanged. Moreover, when it falls off as 1/a 2 , like classical radiation or as shown to be the case for the class of initial state perturbations described in appendix D, the curvature perturbations due to scalar and vector perturbations fall off faster (by a factor 1/a) than those due to tensor perturbations, which become the dominant contribution at late times.

VIII. GENERAL CONFORMAL FIELD THEORIES
The results found for free conformal scalar fields in the previous sections can be straightforwardly generalized to any CFT for the matter sector (even strongly coupled ones). This is because when the background g (0) µν is Minkowski spacetime, the key ingredient for obtaining the linearized stress tensor expectation value , and hence the right-hand side of eq. (7), is entirely determined (up to a constant factor) by conformal as well as Poincaré invariance. Moreover, it transforms in a relatively simple way under conformal transformations and can be easily extended to the case of metric perturbations around a FLRW background, as described in some more detail below. The stress tensor expectation value can be obtained by functionally differentiating the so-called CTP effective action S eff [g, g ′ ] [54][55][56]: The effective action S eff [g, g ′ ] can be written as where the third term which includes nonlocal contributions results from functionally integrating out the matter fields. Performing those path integrals for the matter fields gives rise to UV divergences and one needs to introduce an appropriate regularization procedure; dimensional regularization is a good choice because it is compatible with general covariance (and in a number of cases with conformal invariance as well). Such divergences in Σ[g, g ′ ] can be absorbed by local counterterms in the bare gravitational action. For massless fields the cosmological constant and the Einstein-Hilbert term do not get renormalized in dimensional regularization and only counterterms quadratic in the curvature are necessary, which have been denoted by S div [g] in eq. (60). More specifically, the gravitational counterterms for a generic CFT on a curved background [9,57] are given in dimensional regularization by where d denotes the spacetime dimension and E 4 is the integrand of the Euler invariant defined in eq. (30). Although the (regulated) bare Σ[g, g ′ ] is invariant under conformal transformations of the metric 2 (see the ap-pendix in ref. [58]), the counterterms in S div [g] are not, which is the origin of the trace anomaly. In fact, the constant parameters b and b ′ in eq. (61), which take specific values for each CFT, are also the coefficients of the Weyl-squared and E 4 terms in the trace anomaly (the nonvanishing trace of the quantum stress tensor for conformal fields). So far the statements in this paragraph are valid for an arbitrary metric g. If one is, however, interested in the linearized stress tensor for metric perturbations around a background g µν one only needs to consider terms in S eff [g, g ′ ] quadratic in the metric perturbations or, equivalently, focus on its second functional derivative, which is related to the two-point function of the stress tensor in the background geometry.
Let us now show the universality of the linearized stress tensor expectation value in two steps. The first step is to establish it for metric perturbations around a Minkowski background. By requiring conservation of the linearized stress tensor together with Poincaré and conformal invariance, it has been shown [59] that the form of the properly renormalized two-point function of the stress tensor in a Minkowski background for any conformal theory is essentially unique (up to a constant factor), and so is the renormalized vacuum expectation value of the stress tensor in the linearly perturbed metric. Note that for perturbations around a Minkowski background the second term on the right-hand side of eq. (61) is at least cubic in the metric perturbations and does not contribute to the renormalized expectation value of the linearized stress tensor operator: it, therefore, depends only on the constant b.
The second step is to extend this result to the case of perturbations around a FLRW background via a conformal transformation. As already mentioned above, the counterterms in eq. (61) are not invariant under a conformal transformation of the metric and lead to the following change of the effective action: A conformal transformation leaves the two terms within the square brackets on the right-hand side of eq. (61) invariant in four dimensions and gives rise to terms of order (d − 4) otherwise, which amounts to a finite contribution in eq. (61). The difference between the gravitational counterterms in two conformally related geometries is, therefore, finite. The difference between the Weyl-squared terms gives a term proportional to C µνρσ C µνρσ ln a in the effective action, whose functional derivative corresponds to the first term inside the second square bracket on the right-hand side of eq. (7). On the other hand, since the E 4 term in eq. (61) vanishes up to quadratic order in the metric perturbations for perturbations around Minkowski space, only theẼ 4 term for the conformally transformed metricg contributes at this order. Its functional derivative corresponds to the terms inside the first square bracket on the right-hand side of eq. (7), which are conserved because they result from functionally differentiating a diffeomorphisminvariant integral. In addition to a finite contribution to the coefficient of the Weyl-squared term, which can be absorbed by redefining the renormalization scaleμ, the coefficient of the squared Ricci scalar term can take arbitrary finite values. Its functional derivative corresponds to the first term on the right-hand side of eq. (7) and it generates a term˜ R in the trace anomaly.
In summary, the linearized semiclassical equation (7) will have the same form for any CFT and only the numerical coefficients in front of the two square brackets will change depending on the values of the parameters b ′ and b, respectively.

IX. CONCLUSIONS
Employing the method of order reduction, we have solved nonperturbatively the semiclassical Einstein equation governing the dynamics of linear metric perturbations around de Sitter spacetime when the quantum backreaction of conformal scalar fields on the mean geometry is included. Our exact solutions establish the stability of de Sitter with respect to general linear metric perturbations (of scalar, vector and tensor type) and extend some of the existing "no-hair" results for de Sitter in classical general relativity to the case in which the effects of the quantum vacuum polarization of conformal fields on the semiclassical geometry are included. Indeed, we confirm the late-time attractor character of de Sitter space (for geometrical properties within a region of fixed physical size) by showing that the perturbations of the Riemann tensor, which characterizes entirely the local geometry, fall off with an inverse power of the scale factor.
Perturbative solutions of the semiclassical equation for tensor perturbations have recently been obtained [17] and it was found that the correction to the classical solution can grow arbitrarily large for modes spending a long time inside the horizon. In contrast, our exact nonperturbative solutions exhibit oscillations with decaying amplitude inside the horizon and reveal a breakdown of perturbation theory for long times inside the horizon due to secular terms that arise when expanding perturbatively oscillatory factors with a perturbatively corrected frequency. In addition, we have considered the effects due to perturbations of the initial state of the matter fields. In fact, in order for the state of the fields to continue being a Hadamard state (with no unphysical excitations of arbitrarily short-wavelength modes) when metric perturbations are present at the initial time, it is in general necessary to correct the states that one would consider in the absence of perturbations and consider instead properly "dressed" states which are adiabatic on the perturbed geometry at sufficiently high order.
As explained in sec. VIII, our results are applicable to general conformal field theories with arbitrarily strong self-interaction couplings. It would also be interesting to extend the results for free scalar fields to nonvanishing masses and arbitrary curvature couplings, particularly to massless (or sufficiently light) minimally coupled fields, which typically give rise to larger IR effects in de Sitter. We plan to return to this question in future investigations.
Besides semiclassical perturbations of the mean geometry, it is interesting to study the quantum fluctuations around it, which can be achieved by quantizing the metric perturbations around the mean geometry and dealing with them as an effective field theory. This was done in ref. [33], where the two-point quantum correlation function for tensor metric perturbations around de Sitter, including one-loop corrections from conformal fields, was calculated. A suitable method for selecting the adiabatic vacuum of the interacting theory was employed and twopoint functions compatible with de Sitter invariance were obtained. In particular, all secular terms cancelled out despite being a perturbative calculation. The reason for that can be understood with a simpler but qualitatively similar example. Let us consider a system whose dynamics is governed by a time-independent Hamiltonian with a small interaction term, and let us focus on the following two-point correlation function: where A and B are time-local operators in the Schrödinger picture. Evolving an arbitrary initial state |Ψ 0 for a long period (t 1 − t 0 ) will require in general a nonperturbative calculation. However, for eigenstates of the full Hamiltonian (including the small interaction term) one simply gets a phase factor,Û (t 1 , t 0 )|Ψ 0 = e −iE(t1−t0) |Ψ 0 , which cancels exactly with the complex conjugate counterpart arising from the evolution of Ψ 0 |. Such an exact cancellation also implies a cancellation at every order in perturbation theory and guarantees the cancellation of any possible secular terms associated with the period (t 1 − t 0 ) which may arise in a perturbative calculation at finite order provided that an energy eigenstate of the full Hamiltonian is properly selected as the initial state. Nevertheless, additional secular terms may arise due to the time evolutionÛ † (t 2 , t 1 ) in case of large time differences between the arguments of the correlation function. Thus, a perturbative calculation of the groundstate correlation function will be valid for small (t 2 − t 1 ) no matter how large (t 1 − t 0 ) is, but will break down for large (t 2 − t 1 ). The situation is analogous for metric perturbations around de Sitter. In general a nonperturbative calculation is needed to evolve arbitrary initial states other than the adiabatic vacuum for a sufficiently long time, but it is also required in order to calculate loop corrections to the two-point function with respect to the adiabatic vacuum which are valid for large invariant intervals (both for time-like or space-like separations). In this respect, the connection provided by stochastic grav-ity [60] between the solutions of the semiclassical Einstein equation and the two-point quantum correlation functions for metric perturbations with (resummed) matter loop corrections [61] suggests that the nonperturbative semiclassical solutions found here, and the methods employed to obtain them, could be exploited to compute matter loop corrections to the two-point function of the metric perturbations valid for arbitrarily large separations. Given the perturbed metric g µν and expanding through quadratic order in the metric perturbation we have where indices are raised and lowered with the unperturbed metric η µν , i.e. we regard h µν as a tensor field in flat space. For the Christoffel symbols we get The calculation of the curvature tensors can be done straightforwardly and we obtain The Weyl tensor is given by In four dimensions it obeys the identity which can be proved [62,63] by expanding the equality Alternatively, it can also be proved using the Gauß-Bonnet theorem. (C9) Appendix D: Regular initial states In deriving the linearized semiclassical Einstein equation (7) from the renormalized CTP effective action [42,43], a number of integration by parts with respect to the spacetime variables y µ were performed in order to write the nonlocal term as it appears in eq. (7). After renormalization, the kernel H(x − x ′ ;μ) in the effective action is a well-defined distribution provided that the metric perturbations fall-off sufficiently fast at spatial infinity and at the asymptotic initial time (there is no such requirement for the asymptotic future due to the causal nature of the kernel). However, when giving initial conditions at some finite time η 0 , the linearized stress tensor expectation value gets boundary contributions from that lower integration limit which diverge when the stress tensor is evaluated at η 0 .
Working in Fourier space for the spatial coordinates, the nonlocal term in eq. (7) can be written as where the Fourier transformed kernel is given bỹ see ref. [33] for its computation and for further details. The structure of the nonlocal term for each spatial Fourier mode is then analogous to that for the case of spatially isotropic and homogeneous metric perturbations studied in refs. [22,23]. There it was shown in detail that boundary terms at η 0 arise when writing the result in the same form as in eq. (D1) and that they give a contribution to T (1) µν (η) which diverges for η = η 0 . Furthermore, it was clarified that the reason for such divergences is the fact that although the Bunch-Davies vacuum is a Hadamard state [9] with regular UV behavior in de Sitter spacetime, in general that is no longer the case in a perturbed geometry: with respect to well-behaved adiabatic vacua associated with this geometry it exhibits excitations of modes with arbitrarily short wavelengths.
In refs. [22,23] a simple method was employed for constructing a family of properly "dressed" Gaussian initial states which are regular on the perturbed geometry. The states are prepared by evolving an asymptotic Bunch-Davies vacuum state from −∞ to η 0 in a given (nondynamical) perturbed geometry which is asymptotically de Sitter and matches the dynamical geometry at η 0 . The metric perturbations during this preparation period can be fairly arbitrary, which allows the generation of a wide family of Gaussian states, but need to fulfill a few requirements: they need to fall off sufficiently fast as η → −∞, so that the time integral in the nonlocal term converges, and they need to be small enough so that their contribution to T (1) µν (η) can be treated as a small perturbation. Moreover, the matching at η 0 between the nondynami-cal metric perturbations during the preparation period and the dynamical ones has to be smooth enough: up to fourth order, which is the maximum number of time derivatives that can appear in A µν . In fact, as shown in detail in refs. [22,23], requiring a smooth matching up to this order is exactly equivalent to demanding that the states generated in this way are of fourth adiabatic order, the standard requirement for regular states with a finite renormalized stress tensor expectation value [9]. When the preparation method described in the previous paragraph is employed, the nonlocal contribution in eq. (D1) can be naturally separated into two contributions which result from the following splitting of the time integral: The first term on the right-hand side will contain the dynamical metric perturbations in the integrand and it will vanish when using order reduction since A µν vanishes in that case, as described in sec. IV C. On the other hand, the second term will give a finite contribution in eq. (D1), even for η = η 0 , provided that A µν is regular and vanishes at η = η 0 , which follows from the condition of sufficiently smooth matching required above. This term gives a contribution to the right-hand side of eq. (7) which can be interpreted as a perturbation of the stress tensor expectation value associated with the modified initial state. In fact, the equation can then be rewritten as eq. (47) with δT µν (η, p) = 3α a 2 (η) η0 −∞ A µν (η ′ , p)H(η − η ′ , p;μ) dη ′ .
(D4) The integral is finite for η ≥ η 0 (remember that the metric perturbations are required to fall off sufficiently fast as η → −∞ so that the lower integration limit is convergent) and so is its limit η → 0. Therefore, at late times δT µν decays like 1/a 2 , clearly fulfilling the requirement of secs. VII A-VII B so that the inhomogeneous solutions of the semiclassical equation do not alter the conclusions about the semiclassical stability of de Sitter spacetime in this context.
In Starobinsky's inflationary model the expansion is driven by the trace anomaly of the quantized matter fields, so that the actual value of Starobinsky's Hubble parameter H is different from our H, where the expansion is driven by the cosmological constant Λ. Nevertheless, this has no effect on the form of eq. (E3).