The Cosmological Semiclassical Einstein Equation as an Infinite-Dimensional Dynamical System

We develop a comprehensive framework in which the existence of solutions to the semiclassical Einstein equation (SCE) in cosmological spacetimes is shown. Different from previous work on this subject, we do not restrict to the conformally coupled scalar field and we admit the full renormalization freedom. Based on a regularization procedure, which utilizes homogeneous distributions and is equivalent to Hadamard point-splitting, we obtain a reformulation of the evolution of the quantum state as an infinite-dimensional dynamical system with mathematical features that are distinct from the standard theory of infinite-dimensional dynamical systems (e.g., unbounded evolution operators). Nevertheless, applying methods closely related to Ovsyannikov's method, we show existence of maximal/global solutions to the SCE for vacuum-like states, and of local solutions for thermal-like states. Our equations do not show the instability of the Minkowski solution described by other authors.


Introduction
The semiclassical Einstein equation (SCE) is the equation where G µν denotes the Einstein tensor for the (classical) metric g µν , κ is a gravitational coupling constant constant (conventionally, κ = 8πG), and 〈T ren µν 〉 ω is the renormalized stress-energy tensor for a quantum field theory (QFT) in the state ω. That is, matter is described by a quantum field and gravity is described by a classical Lorentzian manifold. The SCE has been studied since the early 1960's by a number of authors, see [10,27,61] for an overview. It is typically introduced in an ad hoc manner as a minimal change of the classical Einstein equation by replacing the stress-energy tensor of a classical field by that of a quantum field to take into account the quantum nature of matter. In particular, it is not considered a fundamental equation but rather an approximation of a more fundamental theory within some domain of validity that is sufficiently remote from the Planck scale. Some possible derivations of the SCE from a quantum gravity are critically discussed in Sect. II.B of [21].
Many equations in QFT are plagued by ultraviolet divergences and the SCE is no different, because (naïvely) the expectation value of the stress-energy tensor involves the evaluation of singular quantum fields at a point. As already discussed in [59], a renormalization of the stress-energy tensor needs to be coordinate independent and thus has to follow the principle of general covariance.
A procedure which satisfies these requirements is the point-splitting formalism by Christensen [12,13]. Here one subtracts the singular part, given by the Hadamard parametrix (essentially a Lorentzian version of the heat kernel), from the two-point function of the state. For this reason one restricts the class of states to Hadamard states, viz., states with two-point functions that match the Hadamard parametrix up to smooth contributions. (Actually it is sufficient that the leading singularities of the two-point function match those of the Hadamard parametrix, as in the case of adiabatic states [33].) Due to the ambiguity in the regularization procedure (satisfying certain conditions or axioms [29][30][31]), a renormalization freedom arises. Some terms renormalize the gravitational constant or the cosmological constant by a finite amount. Thus it can be seen that the widespread belief that quantum matter automatically leads to a very large cosmological constant is not correct; instead the cosmological constant corresponds to a renormalization freedom. Other terms contain higher than second order derivatives in the metric. Such terms are likely to change the entire characteristic of the SCE, especially with respect to the classical Einstein equation, which is only of second order, and there seems to be to be no justifiable reason why these higher order terms should be discarded [60]. We would like to mention, however, the method of order reduction by which higher order derivatives in the SCE can be replaced by lower order derivatives in a systematic way, see [21] and references therein.
Moreover, it was noted early on, see e.g. [12,60], that the renormalized stress-energy tensor is not traceless, even if the classical action of the quantum field is conformally invariant. In fact, in this case one obtains the famous trace (or conformal) anomaly where the constants c • depend on spin of the (free) quantum field. Notably, the right-hand side does not depend on the quantum state at all, and contains higher than second order derivatives in the metric.
In the following we shall specialize our discussion to the special case of a free scalar field on flat Friedmann-Lemaître-Robertson-Walker (FLRW) spacetimes M = I × R 3 , I ⊂ R, with the metric (in conformal time τ and with signature convention −+++) where a(τ) > 0 is called the scale factor and dx 2 denotes the Euclidean metric on R 3 . This case already shows many features that distinguish it from QFT on the maximally symmetric Minkowski and de Sitter spacetime, or static spacetimes. Additionally, these spacetimes are of importance in cosmology as they represent the observed homogeneous and isotropic structure of our universe at the scale of several megaparsec, as well as the observed flatness [62].
Despite some mathematical problems in the pre-1990's literature on the SCE, by the beginning of the 1980's the approach has been developed to a stage where numerical solutions and cosmological applications of the SCE were in reach. This situation was exploited by Anderson in a series of four papers [1][2][3][4] starting with the conformally coupled and massless scalar field following up prior work by Starobinski [54]. Depending on the values of the aforementioned renormalization parameters, Anderson discovered a very rich behavior of the SCE ranging from big bang, big bounce to divergence of the scale factor to infinity in finite time. The non-conformally coupled case was investigated by analytical and numerical methods by Suen [57,58]. In his work he discusses an instability of Minkowski spacetime as a global solution to the SCE (in the sense of continuous dependence of the solution on its initial data). Further analytical and numerical work on the stability of the SCE has been performed by Hänsel and Verch [28].
As shown by Fulling, Sweeny and Wald in [23] (with some improved results by Radzikowski [49,50] and others), a state which is Hadamard on a time-slice around a Cauchy surface is Hadamard on the entire globally hyperbolic manifold containing this surface. From the perspective of the SCE, this result gives an important hint for the well-posedness of the former equation ensuring that the stress-energy tensor will remain well-defined on the entire spacetime manifold.
The proper description of the stress-energy tensor in the mathematically rigorous framework paved the way for new cosmological investigations on the SCE. Dappiaggi, Fredenhagen and Pinamonti [14] investigated the stability of the SCE on FLRW spacetimes using certain effective large-mass approximations of the quantum state. As a major milestone in the mathematical theory of the SCE, Pinamonti proved the existence of solutions for the trace equation for short times in the conformally coupled case for certain states defined on past null infinities [47]. This has been significantly extended by Pinamonti and one of us [48] to full solutions of the SCE on FLRW spacetimes, including the energy constraint with initial values specified on a Cauchy surface, but still limited to the conformally coupled scalar field, and a particular choice of renormalization parameters and quantum state. Other recent work related to the SCE include: [26], where numerical solutions to the equations presented here are analyzed; [40], where yet another approach to solving the SCE is developed; [36], where a simplified semiclassical problem is studied as an initial value problem; [51] and [35], where solutions in static spacetimes are analyzed.
We remark that in this paper we deal with the SCE without a classical Klein-Gordon field. However, as such a field can be important in inflationary cosmology (see e.g. [27,42,62], we give remarks on the minor changes that are needed to include it without changing the results of our paper.

Outline
After this introduction, in the second section we first outline the quantization of the (real, free) scalar field in curved spacetimes with emphasis on flat FLRW spacetimes. In particular, we present an initial value formulation for homogeneous isotropic states for the quantum scalar field in FLRW spacetimes [39]. Then, we develop a point-splitting regularization specially adapted to FLRW spacetimes and show its equivalence to the conventional Hadamard pointsplitting.
In the third section, we briefly discuss the renormalized stress-energy tensor for the quantized scalar field. Since the trace and the energy component of the stress-energy tensor play an important role for the (semiclassical) Einstein equation in FLRW spacetimes, we present the corresponding expressions. We also discuss the problem posed by higher derivatives present in the semiclassical theory, partially due to regularization and renormalization.
Next, in the fourth section, we show how to formulate and solve a dynamical system for the coincidence limit of the regularized two-point function and its derivatives. This system shows various interesting physical and mathematical features: First, it hides the higher derivatives present in the regularization, thereby circumventing one of the challenges faced when solving the SCE. Next, it distinguishes a class of 'generalized' vacuum states from other classes of states such as thermal states. Finally, its evolution is given by a (generally unbounded) evolution operator which can be understood as acting between differently weighted sequence spaces. As was pointed out to us after the completion of this work, the construction of the evolution operator can be seen as an application of Ovsyannikov's method [22,44] of solving Cauchy problems in scales of Banach spaces. 1 The fact that this evolution operator exists at all is intimately tied to the hyperbolicity of the Klein-Gordon operator, which expresses itself in our dynamical system through the nilpotence of a certain matrix. Depending on the choice of the weights in the sequence spaces, the evolution can be shown to exist for all time (geometrically growing weights) or for a finite amount of time (factorially growing weights).
Finally, in the fifth section, we use the just developed dynamical system for the quantum state to solve the SCE. In fact, we consider an abstract class of quasilinear equations which includes as a special case the SCE. For this abstract equation we show existence and uniqueness of solutions, as well as continuous dependence on the initial values and parameters of the equation. Existence is shown in finite time intervals with a priori bounds depending on the choice of initial values, in particular the initial values for the quantum state. We analyze this equation for various possible choices of parameters in the case of the SCE in FLRW spacetimes. In general, the resulting equation is of fourth order but in the special case of conformal coupling it can reduce to a second order equation. We also remark on the instability of Minkowski spacetime stating that arbitrarily small perturbations in the initial conditions can lead to finite effects in the solution of the energy equation [57,58] -no such instability appears in our approach based on the trace equation, which we show to depend continuously on the initial data. At the end, we explain a method of constructing physical initial data and give a short outlook on future research topics. Let us state our main results, summarizing the contents of Thms. 5.7 and 5.11: Theorem. There exists a non-empty set of self-consistent initial data for the SCE for flat FLRW spacetimes, and, given such self-consistent initial data, the SCE has locally a unique smooth solution which depends continuously on the initial data and parameters. Furthermore, the quantum state part of the solution is of Hadamard type. For a class of 'vacuum-like' initial data for the state, the solution is maximal or even global.
In the appendix we list several auxiliary results (concerning homogeneous distributions, weighted sequence spaces, combinatorial inequalities and Synge's world function) used throughout this work.

Klein-Gordon equation
Consider the (homogeneous) Klein-Gordon equation with mass m ≥ 0 and curvature coupling ξ (ξ = 0 is called minimal coupling and ξ = 1 6 conformal coupling). We define the d'Alembert operator as □ := −g µν ∇ µ ∇ ν , where ∇ is the covariant derivative and we employ Einstein summation convention. In conformal time τ, the flat FLRW metric takes the form (1.2). The d'Alembertian for this where we denote derivatives with respect to conformal time by primes and ∆ is the Laplace operator on R 3 . Introducing the conformally rescaled field ϕ and the potential V , the Klein-Gordon equation (2.1) can thus be written as (

Quantization
For the quantization of the (real, free) scalar field φ on a globally hyperbolic spacetime (M , g), we follow the algebraic approach, see e.g. [17,27,38]: We generate a (non-commutative, where f , g ∈ C ∞ c (M ) and α, β ∈ C. Moreover, 〈· | ·〉 is the canonical L 2 product on the spacetime, and G PJ denotes the uniquely defined Pauli-Jordan propagator [9,16], namely the difference of forward (retarded) and backward (advanced) propagator of the Klein-Gordon operator K.
A state on A is a linear functional ω : A → C, which is The two-point function ω 2 of ω is defined as Due to the positivity of the state we have ω 2 ( f , f ) ≥ 0. Meanwhile, the canonical commutation relations imply On Minkowski spacetime, and more generally on static spacetimes, one additionally imposes a spectrum condition (positive frequency condition) for the state, which distinguishes a vacuum state. On generic spacetimes, no 'natural' preferred state exists [20] and (for free fields) the spectrum condition is replaced by a condition on the singular structure of the two-point function. Typically one requires that the state is a Hadamard state, viz., its two-point function satisfies the microlocal spectrum condition -a condition on the smooth wave-front set [49]. In applications it is sometimes useful to relax the microlocal spectrum condition. For example, on FLRW spacetimes one often considers the class of adiabatic states which are obtained via a WKB-type approach [33,45]. These states satisfy a Sobolev version of the microlocal spectrum condition [33]. Remark 2.1. Here and in the following we assume that the state ω has a vanishing one-point function ω(φ( f )) = 0. Thus we do not distinguish between the two-point function and the connected two-point function. A non-vanishing one-point function φ bg (τ, x ) := ω(φ(τ, x )) is interpreted as a classical Klein-Gordon 'background' field. In the context of homogeneous and isotropic spacetimes, φ bg (τ) does not depend on x . Setting ϕ bg := aφ bg and π bg = ∂ τ ϕ bg , we see that the dynamics of the additional two degrees of freedom introduced by the classical Klein-Gordon field is given by (2.2), where the spatial Laplacian ∆ can be omitted.

Two-point functions
Consider the two-point function ω 2 of a homogeneous and isotropic state. That is, it holds that Since the antisymmetric part of the two-point function is fixed by the commutator 'function' (viz., the Pauli-Jordan propagator G PJ ), ω 2 is completely determined by its symmetric part Therefore, the Cauchy data at conformal time τ of a homogeneous and isotropic state can be given by (2.4) and its first time derivatives. We define which represents the Cauchy data of the two-point function at τ. Using Synge's rule (to pull derivatives into the coincidence limit for the time variables) and the Hamiltonian form of the Klein-Gordon equation (2.2), one can then rewrite the equations of motion for the two-point function as where ∆ r := r −2 ∂ r r 2 ∂ r is the (three dimensional) radial Laplacian. It is sometimes convenient to perform calculations in momentum space. We define the mode functions with k ∈ [0, ∞). The mode functions are simply the Fourier transform of G(τ, r = |x |) in x . Indeed, using the conventionf (k) := R 3 f (x )e −ik·x dx for the Fourier transform, we have for radial functions (or distributions) Conversely, the mode functions specify the two-point function by It follows immediately from (2.6) that the mode functions solve the dynamical equations A straightforward computation shows that J := G ϕϕ G ππ − G 2 (ϕπ) is a conserved quantity. This reduces the degrees of freedom in (2.5) (and (2.9)) to two. Moreover, it follows from the positivity of the state that J ≥ 1 4 with equality for pure states [39].

A point-splitting regularization
Many expressions of physical relevance in QFT involve expectation values of products of quantum fields at a point. Naïvely, such expressions are ill-defined because of the distributional nature of quantum fields. By restricting to a class of states which share a common singular structure, we can define a renormalization scheme that allows to make sense of these expressions. Below we develop a regularization scheme for two-point functions on FLRW spacetimes (or, in fact, for the kernels G(τ, r)) which carries features of both the Hadamard point-splitting method [12,13] and the WKB-type approach named adiabatic regularization [7,11,45]. More concretely, the aim of this subsection is to define kernels H n (τ, r) such that G(τ, r) − H n (τ, r) is sufficiently regular in the limit r → 0.
Further details about these functions can be found in the appendix (Sect. A.1), although without the constant length scale µ. Tacitly, the so defined distributions and the resulting regularization scheme depend on µ.
We make the Ansatz (equivalently either in position or momentum space with relations analogous to (2.7) and (2.8)) where we fix the lowest order parameters to We wish to find coefficients α • , β • and γ • such that the two constraints are satisfied. That is, H n satisfies the Klein-Gordon equation and has the properties of the two-point function of a pure state up to any desired order. Using the homogeneity property (A.4) and (2.11), we find Consequently, the coefficients must satisfy the equations Moreover, the left-hand side of (2.13) evaluates to Therefore, the coefficients must additionally satisfy the constraint It is remarkable that (with this constraint) the differential equations (2.14) can be solved recursively without solving any integrals: Proposition 2.2. The differential equations (2.14) with initial values (2.11) and constraint (2.15) are solved by the recurrence relations Proof. To find (2.16b), we eliminate γ ′ j from (2.14c) by adding the derivative of (2.14b) and using (2.14a). The other two equations obtained from (2.15) and −α j+1 + γ j = V α j + β ′ j , which follows from (2.14b), thus yield (2.16a) and (2.16c).

Comparison to Hadamard point-splitting
The regularization procedure proposed in the previous subsection is equivalent to the typically considered Hadamard point-splitting method, in which a truncation of the Hadamard parametrix is subtracted [12,13]. The truncated Hadamard parametrix (at order n and scale λ) for the Klein-Gordon equation (2.1) is defined as [19,41] H n (x, y) := lim where x, y are in a geodesically convex neighbourhood, is the regularized Synge world function (i.e., half the signed squared geodesic distance) for a time-function t, ∆(x, y) is the van Vleck-Morette determinant, v j (x, y) are smooth coefficient functions satisfying the recurrence relations [19,41] for j ∈ N. These relations are obtained by demanding that the Hadamard parametrix satisfies the Klein-Gordon equation up to an error of order σ n+1 . Note that the coefficient functions v j are entirely determined by the local geometry. For example, at lowest order v 1 is given by Since a FLRW spacetime is spatially isotropic and homogeneous, and the Hadamard parametrix is constructed from local geometric quantities, it inherits these properties so that analogously to (2.3). (Naturally, the same holds true for σ, ∆ 1 2 and v j .) As before for the two-point function, this implies that the Hadamard parametrix can be represented bỹ It follows from results in [18], that the most singular terms ofH n and H n agree, thus justifying the choice (2.11) of the lowest order coefficients. (The results of [18] are stated for cosmological time but a translation to conformal time is not difficult.) Equivalence ofH and H to arbitrary order then follows from the fact that both are approximate solutions of the Klein-Gordon equation with the same singular terms (of the form r − j and r j log r).
For n ≥ 2, one can compute the expansions To obtain this result, one uses a covariant expansion of the coefficient functions ∆ 1 2 and v j , see e.g. [12,13,15], together with an expansion of Synge's world function, see Sect. A.4.
The formulas above suggest the convention where γ is the Euler-Mascheroni constant and λ 0 > 0 is an arbitrary dimensionless constant.
With this convention, we find the differences (n ≥ 2) (2.20d) We will use these differences in the following section.

Semiclassical Einstein equation
On FLRW spacetimes, it is sufficient to look at the traced SCE where 〈T ren 〉 ω := g µν 〈T ren µν 〉 ω is the trace of the quantum stress-energy tensor, and at the energy-component (viz., the 00-or t t-component) The latter equation can be regarded as a constraint equation, which, if imposed at a fixed time, holds everywhere because of the covariant conservation ∇ µ 〈T ren µν 〉 ω = 0 of the stress-energy tensor [31,41], see also Sect. 3.3.
We will focus most of our attention on the traced equation (3.1), which completely determines the dynamics of the geometry, given by the scale factor. However, also the energy equation plays an important role; it is relevant when imposing consistent initial conditions.

Renormalized stress-energy tensor
The (classical) stress-energy tensor for a (real) free scalar field is To quantize this expression, we replace products of (derivatives of) the classical fields by their renormalized quantum counterparts. That is, using Hadamard point-splitting, we have for the expectation value of the stress-energy tensor in a state ω (cf. [31,41,53]): where, for fixed but arbitrary n ≥ 1, we set ω is the coincidence limit of the Hadamard coefficient v 1 ) with implicit parallel transport, c • are renormalization parameters, and I µν , J µν are the two fourth order conserved curvature tensors: These tensors can be obtained as the variations with respect to the metric of R 2 and R µν R µν [60]. For conformally flat spacetimes like FLRW, it holds I µν = 3J µν . The term involving the Hadamard coefficient v 1 in (3.4) ensures that the quantized stress-energy tensor satisfies ∇ µ 〈T ren µν 〉 ω = 0 [31,41]. The term c 1 m 4 g ab is a renormalization of the cosmological constant, and c 2 m 2 G ab corresponds to a renormalization of the gravitational coupling constant κ; the remaining two renormalization terms have no classical interpretation. A change of the scale λ in the Hadamard parametrix corresponds to a change of the renormalization parameters [27], and thus we are in principle at liberty to set λ to any value if we change the renormalization parameters accordingly.
We also remark that in case the state ω includes a non-zero one-point function (classical Klein-Gordon field), the corresponding contribution to the two point function is After performing the coincidence limit, this leads to an additional contribution given by the classical stress-energy tensor (3.3) with φ replaced by φ bg .

Trace of the stress-energy tensor
Taking the trace of (3.4) yields 2 where c • are the same constants as above, and we used that For homogeneous states and isotropic states (i.e., satisfying (2.3)) on FLRW spacetimes, we thus find where and (2.18) specializes to In the case of a non-vanishing one-point function yielding the classical field φ bg , the following expression needs to be added to 〈T ren 〉 ω :

Energy component of the stress-energy tensor
Next, let us state the energy component of the stress-energy tensor in a homogeneous and isotropic state on an FLRW spacetime: Consequently we also have that is, we can express the trace in terms of the energy component if a ′ is bounded away from zero. If (3.1) holds, then (3.9) and (3.10) imply that holds. In particular, this implies that the constraint equation (3.2) given by the energy component of the SCE is satisfied everywhere if it is satisfied at one instance of time. Note also that (3.11) implies that G 00 − κ〈T ren 00 〉 ω is asymptotically pushed to zero if a ′ > 0. If a non-vanishing one-point function is present, it contributes to 〈T ren 00 〉 ω the following expression:

Problem of higher derivatives
Note that the SCE contains, via the renormalized stress-energy tensor, up to fourth order derivatives of the metric and thus, in the case of a FLRW spacetime, fourth order derivatives of the scale factor. This is in striking contrast to the classical Einstein equation, which contains only second order derivatives of the metric. Therefore, it can be argued that the SCE has solutions which diverge significantly from solutions with similar initial data for the classical Einstein equation, see e.g. the discussion in [21]. In particular, it is generally expected that these higher derivatives cause additional unphysical runaway solutions to appear. A naïve strategy to solve the SCE for FLRW spacetimes, i.e., the equation −R = 〈T ren 〉 ω coupled with the dynamics of the quantum state, would be to move the highest order derivatives of the scale factor to the left-hand side and all remaining terms to the right-hand side. Unfortunately this is not possible (in the case of non-conformal coupling) because the regularization includes singular terms with fourth order derivatives of the scale factor, and thus it is not clear how to proceed, cf. [56] where the same problem is discussed. For this reason, earlier approaches to solving the SCE, e.g. [48], focused on the conformally coupled case, where this problem can be avoided.
Actually, if one wishes to solve the SCE such that the state is Hadamard, the problem is even more severe. In that case there is no obvious possibility to specify initial values for the state such that it is Hadamard as long as the spacetime is not known in a neighbourhood of the Cauchy surface. In the 'nicest' scenario, the analytic case, one needs to know all the derivatives of metric (here, the scale factor) at the Cauchy surface. A possible strategy to avoid this problem is to work with adiabatic states, as done e.g. in [48], but it comes at the cost that the solutions cannot be guaranteed to be smooth.
Below we suggest an alternative approach which relies on a dynamical system for the regularized two-point function and its derivatives in the coincidence limit and thus avoids the above mentioned problems.

Stability of solutions
There are various notions of stability for solutions of differential equations, the most basic but also the most important being the continuous dependence of the solution on its initial data and the parameters 3 in the differential equation. Together with existence and uniqueness of solutions, continuous dependence on the initial data is required for a well-posed problem in the sense of Hadamard.
A differential equation can typically be rewritten in many more or less equivalent forms. In the case of the SCE on FLRW spacetimes, we have seen in (3.10) that the trace can be expressed in terms of the energy component. In the derivation of that equation we divided by a ′ , and thus it is only valid if a ′ is bounded away from zero. Hence, solving SCE in terms of the energy equation leads to problems in the neighbourhood of Minkowski spacetime where a ′ = 0. This fact is related to the failure of well-posedness of the energy equation near the Minkowski solution described e.g. in [57,58] as we will also discuss in Sect. 5.3. Since this problem does not appear if the trace equation is taken as the governing dynamical equation, it should be preferred over the energy equation in the generic case.

Dynamical system for sequences of coincidence limits 4.1 Dynamics
For n ∈ N 0 , we define That this definition is indeed independent of l follows immediately from (A.4) and the definition of H n : Occasionally, we call M n moments. To understand this nomenclature, consider the momentum space representation of M n : That is, the sequence of M n is given by the moments ofĜ −Ĥ. Note, however, thatĜ −Ĥ cannot be expected to be positive.
To formulate a dynamics for M n , recall (2.6) and (2.12). It follows that where we defined Sometimes we write A(τ) to emphasize the dependence on τ for a fixed potential V . Henceforth we shall suppose that V is an arbitrary (but sufficiently regular) function and can forget its relation to the Klein-Gordon equation.
Considering M n as a sequence M = (M n ), we can also write the equation above as where L denotes the left-shift operator. Below we study this equation on the weighted sequence spaces where p ≥ 1 and w is a sequence of weights; we denote the norms by ∥ · ∥ p,w . See Sect. A.2 for an introduction and our conventions. Note in particular that in our convention the weights appear as inverses in the norm and thus they directly translate to the maximum growth rate of elements in the sequence space.
The careful reader will have noticed that above we have essentially performed a Taylor expansion of G − H in the radial variable r. Indeed, as we will see below, our approach can be applied in cases where that difference is analytic, with the best result in the entire case. Therefore our resulting reformulation of the initial value problem for the SCE can partially be solved with methods related to the Cauchy-Kovalevskaya theorem. However, unlike in the classical Cauchy-Kovalevskaya theorem, we do not require analyticity in the time variable (continuous differentiability of sufficient order is enough) and thus our methods more closely resemble the more abstract ones by Ovsyannikov [22,44]. For the same reason, the subset of moments that we consider cannot be compared with the set of analytic Hadamard states (i.e., Hadamard states which satisfy the analytic wavefront set condition) [55] except under additional assumptions including the analyticity of the scale factor. However, one can expect that requiring analyticity also in time would be too restrictive for most physically interesting solutions.
In the following two subsections, we calculate M in two relevant examples to motivate the choice of weights later on. Afterwards, the remainder of this section is concerned with solving (4.2). The solution is given by a mathematically rigorous definition of the time-ordered exponential of A ⊗ 1 + B ⊗ L.
Finally, we wish to remark that the set of possible sequences M contains many unphysical examples (in particular, positivity is not guaranteed), and, furthermore, not all possible Hadamard states can be represented as a sequence in (4.3) for the weights considered below.
Finally, we note that in the massless case m = 0 all moments vanish exactly: M = 0.

Properties of the matrices A and B
The matrix B is nilpotent: B 3 = 0. For this reason, many products of the matrices A and B vanish.
We are interested in the non-zero words in A and B with the largest number of B-factors. It is already an easy consequence of the nilpotency of B, that products of length 3n contain at most 2n B-factors (e.g., (B 2 A) n ). This implies that a non-zero word of length n contains at most ⌈ 2n 3 ⌉ B-factors. This estimate can be substantially improved using further relations between the matrices A and B.
In the following, we denote by dots the non-zero entries of a 3 × 3 matrix. For instance, ( · · · · ) denotes any 3 × 3 matrix (a ij ) whose only non-zero entries are a 12 , a 21 , a 23 , a 32 . The matrix A introduced in the previous subsection is an example of such a matrix and we write A = ( · · · · ). Another example is B = (· · ). Products and powers of 3 × 3 matrices will also be represented in this notation. For example, we write ( · · · · ) 3 for A(t 1 )A(t 2 )A(t 3 ). Note that (· · ) 3 = 0.
Proof. By induction, it follows easily from (4.11), that a non-zero word of odd length with a maximal number of (· · )-factors has the shape (· · ), (· ) or ( · ). This also shows that such a (· · )-maximal word contains at most one occurrence of two consecutive (· · )-factors. Hence we find that a non-zero word of length 2n + 1 contains at most n + 1 factors of the shape (· · ). It is then obvious, that a non-zero word of length 2n can also not contain more than n + 1 factor of the shape (· · ). This completes the proof. □ We conclude this subsection on the properties of the matrices A and B by noting that their matrix (max) norms are given by: (4.12)

Evolution operator
In this subsection, we define the evolution operator for as an operator on the weighted sequence spaces ℓ p (w) (see (4.3) and Sect. A.2 for the definition of these spaces) for certain weight sequences w and p ≥ 1.
Recall that evolution operator for S(τ) is the 'solution operator' associated to the solutions of the initial value problem for the dynamical equation That is, defining for n > 0 it is given by the Dyson series (−1) n U n (τ 0 , τ) for τ < τ 0 .
for τ, τ 0 , τ 1 in an appropriately chosen interval. Let us briefly remark on the relation of the spaces in (W1) and (W2) to spaces of analytic functions. Note that ℓ p (w) ⊂ ℓ ∞ (w) for any p ∈ [1, ∞], and hence we will restrict our attention to the case p = ∞. Hence, on the one hand, the moments define a power series of an even (i.e., all non-even powers are zero) analytic function in the neighbourhood of the origin. On the other hand, if G − H is analytic in a neighbourhood of the origin, then the corresponding moments will be contained in an appropriate weighted sequence space of type (W2). If the moments are in a space with geometrically growing weights w n = ω 2n (W1), then the moments grow so slow compared to 1/(2n)! that the radius of convergence of the corresponding power series is infinite, and thus they define an entire analytic function.

Evolution operator for geometrically growing weights
In this subsection we consider the case (W1), i.e., the evolution operator on ℓ p (w) for the geometrically growing weights w n = cω n with c, ω > 0. As described above, in this case G − H must be entire analytic in the radial variable. In this case, S(τ) is a bounded operator on ℓ p (w) with (4.16) The following theorem is easily shown, see e.g. Thms. 5.1 and 5.2 in [46] (the proofs of the analogous Thm. 4.10 in the next subsection can also be adapted to the case of geometrically growing weights discussed here): Theorem 4.4. Let c, ω > 0 and p ≥ 1. Suppose that I ⊊ R and V ∈ C(I). The evolution operator U(τ, τ 0 ), defined by (4.15), is the unique bounded operator on ℓ p (w) with weights w n = cω n , such that the properties (E1)-(E3) hold for τ, τ 0 , τ 1 ∈ I. Moreover, we have the bound 17) and the evolution operator is norm-continuous: 4 lim τ→τ 0 ∥U(τ, τ 0 ) − 1∥ p,w = 0.
Proof. Let us sketch a proof. Using the Banach fixed point theorem and the boundedness of the operators S(τ), it can be seen that the integral equation has a unique solution which is given the Dyson series (4.15). Property (E1) and the first part of (E3) are obvious, while (E2) follows from the uniqueness of the solution. For the second part of (E3) note that (E2) implies U(τ, τ 0 ) −1 = U(τ 0 , τ) and, together with the first part of (E3), we thus find ∂ τ 0 U(τ, τ 0 ) = −U(τ, τ 0 )S(τ 0 ). Next, the bound can be obtained via Grönwall's inequality from the integral equality together with the bound (4.16). Finally, norm-continuity can also easily be derived from the integral equation. □ In particular, this theorem implies the following: Remark 4.5. If we are given initial data of geometrically growing moments (e.g., the moments of the vacuum state on Minkowski spacetime), they will continue to be geometrically growing under time-evolution independent of the concrete potential V (t).
Remark 4.6. Given a continuous potential V , the solution exists in ℓ p (w) for arbitrarily large time intervals I.
The following perturbation result can be shown by a standard application of Grönwall's inequality, but can also be shown along similar lines as Thm. 4.11: Under more stringent assumptions on the potential V , one obtains the following regularity result: Proof. By property (E3), ∂ τ U(τ, τ 0 ) is bounded on ℓ p (w). Repeated differentiation of (E3), together with the boundedness of S(τ) and its derivatives, yields the desired result by iteration. □ This implies, in particular, that the solutions of (4.2) with initial conditions in ℓ p (w) are smooth if V is smooth.

Evolution operator for factorially growing weights
For weights growing faster than the geometric growth considered in the previous subsection, the left-shift operator and also S(τ) are unbounded (see also (A.5)). Moreover, the following shows that a construction based on the standard Hille-Yosida theory of C 0 -semigroups can not be used to resolve this: Proposition 4.9. If the weights grow faster than geometrically, the resolvent set of S(τ) is empty.
Proof. For any λ, ξ ∈ C, because it grows at most geometrically, and a short calculation shows However, even then it is sometimes possible to understand S(τ) and, somewhat surprisingly, also the evolution operator (4.15) as a bounded operator between two different weighted sequence spaces ℓ p (v) and ℓ p (w). In other words, both S(τ) and the evolution operator U(τ, τ 0 ) can be considered as unbounded operators on ℓ p (v) with domain ℓ p (w).
Let us consider the case (W2), i.e., the weights w n = (2n)!ω 2n and v n = (2n)!υ 2n with υ > ω ≥ 1. As explained earlier, this can essentially be understood as the case where G − H is analytic (with radius of convergence related to ω and υ) in the radial variable. We apply Prop. 4.3, (4.12) and (4.14) to obtain for n > 0 the bound where we set We compute r 0 = 1 and, for m > 0, odd n, n 2 + 1 even n. Thus, estimating the remaining factors in the sum by their supremum and recalling the definition (4.17) of C ω (τ, τ 0 ), we obtain for n > 0 Consequently, the time-ordered exponential (4.15) exists and can be bounded for sufficiently small time intervals by where we defined Note that these inequalities only guarantee the existence of the exponential for a bounded interval of time, even if ω = 1 and υ is arbitrarily large. The existence of such a bound is quite natural in a hyperbolic problem. Since the radius of convergence is initially bounded, it must shrink as the area of dependence grows with the propagation speed. Using the result above, we can now show: The evolution operator U(τ, τ 0 ), defined by (4.15), is the unique bounded operator from ℓ p (w) to ℓ p (v) with weights w n = (2n)!ω 2n and v n = (2n)!υ 2n , such that the properties (E1)-(E3) hold in the interval I (in the strong sense from ℓ p (w) to ℓ p (v)). Moreover, for M ∈ ℓ p (w), we have the bound (4.20) and the strong continuity property We have already shown that the evolution operator U(τ, τ 0 ) is well-defined as an operator from ℓ p (w) to ℓ p (v). In fact, the same estimates show that, for sufficiently small where v ϵ,n = (2n)!(υ − ϵ) 2n . Moreover, applying these estimates again, we find as τ → τ 0 because C 0 (τ, τ 0 ) → 0 in that limit. Property (E1) is obvious and (E2) is shown (as for groups generated by bounded operators) by multiplying the series (4.15) for U(τ, r) and U(r, τ 0 ). The final property (E3) is proven by differentiating (4.15) term by term using the formal relation ∂ τ U n (τ, Finally, we prove a perturbation theorem: Theorem 4.11. In addition to the assumptions of Thm. 4.10, suppose thatṼ ∈ C(I) such that 5 Denote by U(τ, τ 0 ) andŨ(τ, τ 0 ) the evolution operators for V andṼ . If M,M ∈ ℓ p (w), we have Proof. Setw Further, note that 2C 0 (τ, τ 0 )K(υ, ω) < 1 by assumption. Therefore, Now, note the elementary identity For the first summand we find 5 Note that this condition adds no additional restriction because the role of V andṼ can be exchanged.
For the second summand we find, using the fundamental theorem of calculus (for Banach space-valued integrals), and thus As for the case of geometrically growing weights, under more stringent assumptions on the potential V , one obtains the following: Proof. As in the proof of Thm. 4.10, U(τ, τ 0 ) is bounded from ℓ p (w) to ℓ p (v ϵ ). Moreover, it is not difficult to see that products of S(τ) and its derivatives are bounded from ℓ p (v ε ) to ℓ p (v). (Concrete estimates similar to those before Thm. 4.10 could be computed, but because no infinite sums are involved here this is not necessary.) The result follows now by iteratively differentiating property (E3). □ As before, this implies that solutions of (4.2) with initial conditions in ℓ p (w) are smooth if V is smooth.

Abstract semiclassical Einstein equation
The SCE on FLRW spacetimes contains only one geometric degree of freedom -the scale factor a = a(τ). Therefore, as described in Sect. 3, it turns into an ODE for the scale factor, coupled to a dynamical system describing the evolution of the state. In this section, we develop a scheme to solve such systems, encompassing also a large class of modifications of the SCE. Here we always solve the SCE forward in time but equivalent results hold for solutions backward in time.
For fixed k ∈ N and τ ∈ I ⊂ R, consider the initial value problem for the quasilinear system S(τ, a)M(τ), factor and a = (a, a ′ , . . . , a (k) ) its k-jet, We note that the results in this section generalize to states including classical fields, if one simply includes the additional degrees of freedom from Rem. 2.1 in the system (5.1) and adds (3.7) to f (τ, a, M). These modifications do not change the structure of the proofs below.
After developing this abstract formalism in the following two subsections, it will be applied to traced semiclassical Einstein equation in Sect. 5.3, by choosing k = 3 and a particular function f = f tr coming from the trace of the quantized stress-energy tensor.

Existence and uniqueness of solutions
Existence of solutions to (5.1) can be shown by a (partial) linearization and employing a fixed-point argument via the construction of a contraction map. With the preparatory results from the previous sections at hand, this theorem is a relatively straightforward adaption of standard results (see e.g. [37]) on quasilinear systems to the case of Eq. (5.1). The (standard) proofs are deferred to Sect. A.5.
where we use the Euclidean norm on R k+1 . For fixed R > 0 and weight sequence w, consider a init ∈ B R/2 and M init ∈ ℓ p (w). Suppose that there is I := [τ init , τ fin ] ⊂ R and a weight sequence v such that the following holds: (ii) There is L M > 0 such that, for all a,ã ∈ C(I; B R ), 3) uniformly for τ ∈ I.
Then there exists τ stop ∈ (τ init , τ fin ] such that (5.1) has a unique (local) solution In particular, we can apply the above theorem to the cases studied in the previous section. In the case of geometrically growing weights we can even obtain maximal or global solutions.
Further, suppose that τ → V (τ, b) is continuous for all b ∈ B R , and there is L V > 0 such that

extends to a maximal solution (the solution exists up to a singularity of V or f ) or to a global solution (the solution exists for arbitrarily large times).
Proof. This is shown by gluing local solutions of the initial value problem (using Thm. 5.1 and Prop. 5.2). □ Remark 5.4. In the case of factorially growing weights, the situation is considerably more complicated because, a priori, already the solution for the dynamics of the moments M exists only for a bounded time interval.
In the case of the SCE as described in Sect. 3, V and f (as we will see below) are rational functions in the derivatives of the scale factor a and log(a). Therefore, the assumptions of Thm. 5.1 imply that we can solve (5.1) with the k-jet a of a inside a ball B R and thus away from the poles of V and f . In the case of geometrically growing weights, we even have maximal (or global) solution. Moreover, if V and f are smooth (as in our application), also the solution is smooth: Proposition 5.5. Suppose that (a, M) ∈ C k+1 (I) × C 1 (I; ℓ p (w)) is a solution of (5.1) in the interval I ⊂ R with a ∈ C 1 (I; B R ) (for B R as in Thm. 5.1). If V ∈ C l (I × B R ) and f ∈ C l (I × B R × ℓ p (w)) for some l ∈ N, then a ∈ C k+l (I) and M ∈ C l (I; ℓ p (w)). In particular, if V and f are smooth, the solution is smooth.
Proof. Follows by iterated derivation and substitution of (5.1) together with Prop. 4.12. □ The regularity of the scale factor is especially relevant when discussing the regularity of the state. In fact, the state can only be Hadamard if the scale factor is smooth. In that case, if the state is initially Hadamard, it will be Hadamard for all time, because the Hadamard property is propagated on smooth spacetimes [50].
The proof can be found in Sect. A.5. The continuous dependence of the solution on the initial data also implies that the effect of any error in the initial data is bounded, with an error bound that increases in time. This is of particular importance in this approach, as the state represented via its moments essentially requires an infinite amount of initial data (even if this can in principle be encoded in a single function). The error caused by a truncation of the sequence of moments at some order can thus be controlled.

Application to the traced SCE
In this subsection we show how the results above may be applied to the SCE as discussed in Sect. 3.
Recall that ω reg 2 = ω 2 − H n , where n is sufficiently large. It follows by a straightforward calculation from the definitions of Sect. 2 that Thus, combining the results from Sects. 2 and 3, the traced SCE (3.1) can be expanded to the rather long equation We remark that the first line of this equation is due to terms proportional to □R.
In the general case, this equation can be rewritten as quasilinear fourth order equation of the form which can be solved using Thm. 5.1 (see also Props. 5.2 and 5.3). Note that the right-hand side has poles at a = 0 and for a such that 30(6ξ − 1) 2 log(aλ 0 ) = 11 + 5760π 2 (3c 3 + c 4 ) − 60ξ, (5.8) which must be taken into account for the choice of the ball B R in Thm. 5.1. We do not see the instability near the Minkowski solution described in [57,58], viz., the equation is continuous in a neighbourhood of the Minkowski solution. The reason is, clearly, that we based our analysis on the fourth order equation given by the traced SCE. Indeed, inspection of (3.10) immediately reveals, that the third order equation given by the energy component of the stress-energy tensor can only be used if a ′ is bounded away from zero and in that case it is sufficient to work with the third order equation, see also the next subsection.
However, in the non-conformally coupled case the singularity (5.8) appears. It can be seen that the position of this singularity depends on the choice of λ 0 , relating the length scales µ and λ in (2.19). While the choice of λ is related to the value of the renormalization parameters c 1 , c 2 , c 3 and c 4 , the value of µ (and thus λ 0 ) is related to our construction of M. Of course, changing λ 0 is not without effect and requires a corresponding change of M. It is important to note that such a modification of M will generally cause it to fall out of the sequence space.
There is only one special case in which this equation reduces to a lower than fourth order equation: if ξ = 1 6 and 3c 3 + c 4 = − 1 5760π 2 , (5.9) the fourth and third order terms drop out and the equation can be rewritten as and thus has a the correct form to be solved by the methods above. This equation is equivalent to that already considered in [48]. Note that the right-hand side has poles at a = 0 ('big bang/crunch') and a ′ 2 a 4 = 1440π 2 κ −1 − (1440π 2 c 2 − 5)m 2 , (5.10) viz., for a certain value of the square of the Hubble parameter (a −2 a ′ in conformal time), which must be taken into account for the choice of the ball B R in Thm. 5.1. Also the instability (5.10) can be considered irrelevant because c 2 should be set to zero as it corresponds to a renormalization of Newton's gravitational constant, which has already been measured, and thus, if the scalar field mass is much smaller than the Planck mass, the singularity occurs for a Hubble parameter close to the inverse Planck time (very many orders of magnitude larger than the currently observed value). We sum up the results of this subsection: Let us also briefly explain the philosophy behind the choice of the renormalization parameters c • . The approach described above yields a family of equations for each value of these parameters. The solutions for each parameter can then, in principle, be compared with experimental/observational data to fix these parameters. In this context it should be stressed that the same parameters must be chosen for all spacetimes and all experiments in order for the renormalization procedure to satisfy various standard assumptions, including local covariance [29][30][31].

Application to the energy component of the SCE
Analogous to the traced SCE, the energy component of the SCE can be expanded as In the special case (5.9), already considered above, this equation simplifies to i.e., degree lower than the corresponding equation for the trace. Just like the analogous equation for the traced SCE, this equation can be plugged into the abstract machinery developed in the first half of this section. However this comes at the price that f en has a pole at a ′ = 0 and thus it cannot be applied there.
Recall that the energy component of the SCE must be implemented as a constraint on the initial data. An approach to constructing initial data that satisfies the energy constraint and also the Hadamard condition is discussed in Sect. 5.6.

Minkowski solution
Before we continue to discuss the general case again and study the problem of constructing physically reasonable initial data, let us consider the simplest solution of the equations above: the Minkowski solution.
For proper physical initial data it should be required that the initial data for the state (resp. the moments) satisfy the energy (constraint) equation given by (3.8) or, equivalently, by (5.11). Therefore, on Minkowski spacetime, the following relation needs to hold: where for the last equality we assumed a time-translation invariant state (so that the left-hand side of (4.2) vanishes). A curious observation is that the energy constraint equation fixes the renormalization parameter c 1 .
For for the thermal state. That is, in these examples, c 1 is fixed once the mass m and the regularization scale λ, resp. the inverse temperature β are fixed. Note that the arbitrary scale λ 0 cancels out. It follows then from (3.9) and (3.10), for any time-translation invariant state on Minkowski spacetime, that G 00 = κ〈T ren 00 〉 ω ⇐⇒ −R = κ〈T ren 〉 ω . Therefore we have:

Construction of non-trivial physical initial data
It is not clear how to give initial values for a (Hadamard) state unless the scale function is known in a neighbourhood of the Cauchy surfaces. This is another reason for why it is not clear how to pose a satisfying initial value problem for the SCE. Our use of the moments M instead of a state does not completely solve this problem as it is not clear which sequences of moments belong to positive two-point functions of physical states.
To escape this problem, we propose the following approach: We consider a (fixed) FLRW spacetime on which we know how to construct a Hadamard state, e.g., a spacetime which is Minkowskian in some time interval (in the context of Lem. 5.9 this would be a neighbourhood of τ tow ), and then calculate the corresponding moments. Then, at some point in time (denoted τ init below), we start to switch on the SCE. If this is done smoothly, the resulting solution will be smooth and the propagated state will remain Hadamard. Hence, at any point after the SCE is fully switched on (denoted τ free below), we have a solution to the SCE with a Hadamard state.
To put this idea on a solid footing, we prove the following technical lemma, some aspects of which are also illustrated in Fig. 1: a init a χ Figure 1. Schematic illustration of the tow-in procedure described in Lem. 5.9, including the various relevant time intervals. In the past of τ init the SCE is 'turned off' and the solution is determined by a chosen scale factor a tow which evolves to the desired initial value a init . Then, in the interval (τ init , τ free ), the SCE is smoothly turned on using the function χ. Shrinking this time interval, the solution at τ free can be brought arbitrarily close to its value at τ init . Finally, after τ free the solution evolves freely via the SCE without the forcing caused by a tow .
Lemma 5.9. Let [τ tow , τ fin ] ⊂ R, τ init ∈ (τ tow , τ fin ) and a tow ∈ C k+1 [τ tow , τ fin ]. Suppose that M tow ∈ C 1 ([τ tow , τ fin ]; ℓ p (w)) is a solution of (5.1b) for the scale factor a tow and the weight sequence w. Further, given appropriately chosen b c , R, v, suppose that the assumptions of Thm. 5.1 hold for the initial data a init = a tow (τ init ), M init = M tow (τ init ). Then there exists τ stop ∈ (τ init , τ fin ] such that for any ϵ > 0 one can find τ free ∈ (τ init , τ stop ) and (a, M) satisfying which holds for τ ∈ I, and the right-hand side of (5.4) is unchanged. Hence Thm. 5.1 can be applied to find τ stop ∈ (τ init , τ fin ] and a solution Define (a, M) by gluing (a tow , M tow ) and (a χ , M χ ). Choosing χ such that χ(τ) = 0 for τ ≤ τ init and χ(τ) = 1 for τ ≥ τ free with τ free ∈ (τ init , τ stop ) arbitrary, properties (i)-(iii) are obviously satisfied. Since estimates in the previous paragraph are uniform in χ, (iv) follows by continuity of the solution, viz., if τ free is chosen sufficiently close to τ init , the estimates hold. Finally, in the smooth case, apply Prop. 5.5 to see that (v) holds. □ Applied to the traced SCE, where V and f are smooth for appropriately chosen B R such that a is bounded away from zero, this lemma together with the propagation of the Hadamard property on smooth spacetimes [50] shows: Written as g(a, M, τ) = a ′′′ (τ) − f en (a, a ′ , a ′′ , M) τ with the right-hand side as defined in Sect. 5.4, it is evident that the energy constraint is continuous in a for a > 0 and a ′ ̸ = 0. Since V and f en only depend on lower than third order derivatives of the scale factor, it then follows from Thm. 4.7 and (5.14) that for any γ > 0 there are c ± ∈ R and δ > 0 such that ±g(a c ± ,δ , M c ± ,δ , τ init ) ≥ γ It then follows by (iv) of Lem. 5.9, that one can choose τ free such that ±g(a c ± ,δ , M c ± ,δ , τ free ) > 0. Moreover, given sufficiently small δ > 0, Thm. 5.6 implies the continuous dependence of g(a c,δ , M c,δ , τ free ) on c. An application of the intermediate value theorem thus shows that there exists c 0 ∈ (c − , c + ) such that the associated solution (a, M) = (a c 0 δ , M c 0 ,δ ) satisfies the energy constraint g(a, M) = 0 at τ free .
As the energy constraint is propagated by solutions to the SCE, (a, M) is a proper solution to the SCE (i.e., both the energy and the trace equation) in [τ free , τ stop ]. As the Hadamard property is propagated on smooth spacetimes, M is associated to a Hadamard state, namely the state ω tow propagated on a. This shows (i) and (ii). Finally, (iii) follows from the triangle inequality using (5.14) and (iv) of Lem. 5.9. □ This theorem shows, in particular, that the set of solutions satisfying the Hadamard condition to the SCE is non-empty. In fact, for any choice of approximate initial conditions for the scale factor, a nearby solution of the SCE exists.

Reconstruction of the quantum state
There is no obvious way to directly relate a sequence of moments M to a quasi-free state. Note that this is not a classical moment problem, as the degree l of the counter terms in (4.1) is neither fixed nor unique. Here we can easily circumvent this problem: Suppose that in addition to the initial data for the moments M we are given initial data for the associated state, or rather the two-point function. (Concretely, this is possible, for instance, in the approach described in Sect. 5.6.) Then, once we have solved the quasilinear system (5.1) for the scale factor and the moments, we can evolve the initial data for the two-point function with the obtained scale factor. Necessarily, the evolved two-point function is compatible with the evolved moments. Alternatively, we can augment (5.1) by adding a third equation for the two-point function and directly co-evolve it with the moments.

Outlook
Here we restricted ourselves to the SCE with a free scalar field on flat cosmological spacetimes. An extension of our approach to other non-interacting types of matter (e.g., the Dirac field) and non-flat cosmological spacetime seems feasible. The former requires some modifications of the moment spaces, while the latter requires a careful treatment of homogeneous distributions on maximally symmetric spaces. Currently we are also investigating the specialization of some of our methods to (cosmological) de Sitter spacetime [25], where we expect to extend some of the results of [34].
The introduction of potential energy of the field (neglecting self-interaction) −λ〈:φ 4 :〉 ω would lead to quadratic terms in M ϕϕ and further correction terms in the trace equation that nevertheless still fall in the class of the abstract SCE (5.1). Such modifications could be of interest in the cosmology of the early universe on time scales well above the Planck time but below characteristic times of nuclear reactions, see [5,6] for related work.
A much more ambitious would be the passage to non-cosmological spacetimes. Spacedependent germs of distributional tensor structures that remain closed under successive applications of the Klein-Gordon operator would have to replace the expansion in radially symmetric homogeneous distributions. The point-splitting limit of such an approach would result in an infinite hierarchy of coupled PDEs for the germ coefficients. It is not obvious, if such a system can be set up or even be solved.
Another open question is if it possible to improve upon the techniques employed here to avoid what are essentially assumptions on the spatial analyticity of the state. Since the difference G − H does not satisfy a PDE away from the diagonal, it is currently not clear how that can be accomplished while at the same time still ensuring the correct regularization of the state and without adding arbitrarily high derivatives of the scale factor to system. Working in the opposite direction, one could attempt to develop a fully analytic theory (also in time). Such an analytic theory would inevitably involve analytic Hadamard states, it might offer other ideas to solve the SCE in non-cosmological spacetimes, using e.g. the construction of analytic Hadamard states developed in [24].
It is often argued that the higher derivatives of the metric appearing in the SCE can lead to runaway solutions which deviate significantly from classical solutions of the Einstein equation. Although recent numerical work in collaboration with N. Rothe suggests otherwise [26], it might still be worthwhile to apply our methods to the order-reduced SCE (cf. [21]), which is sometimes suggested as a better behaved approximation of the interaction between quantum matter and classical spacetime geometry.
Potentially, our new formulation of the SCE also has numerical consequences, as it avoids time integration of rapidly oscillating modes of the quantum state and integration in momentum space. Theorem 5.6 establishes certain bounds for the change of the solution under modification of the initial conditions. This can be used to truncate M at a certain order with a controlled error for the solution of the SCE. As the space of moments with zero entries after a prescribed order is invariant under the dynamics (4.2), such approximated initial conditions give rise to a finite dimensional ODE which can be numerically integrated in time using standard methods. However, convergence rates as a function of the time interval need to be carefully examined to judge numerical viability.
Defined in this way, k z + satisfies for z ∈ C the homogeneity property 〈k z + , k f 〉 = 〈k z+1 + , f 〉. Note that (A.2) is not the only possible extension of (A.1) to all of z ∈ C -different extensions differ by derivatives of the delta distribution at zero.
To calculate the (inverse) Fourier transform of k −n + , note first that Then we compute where we recall that ∆ r denotes the (three dimensional) radial Laplacian.

A.4 Expansion of Synge's world function
To compute the (truncated) Hadamard parametrix in a given spacetime, it is necessary (among other things) to find Synge's world function σ. The world function σ(x, y) is defined as half the square signed geodesic distance between the points x and y (if the two points lie in a geodesically convex neighbourhood) and satisfies the relation 2σ = (∇ µ ⊗ 1)σ (∇ µ ⊗ 1)σ . (A. 10) In fact, the relation (A.10) together with the coincidence limits uniquely defines Synge's world function. An expansion of Synge's world function σ(x, y) in terms of the coordinate distance δx between the points x and y can be obtained in the following way [53]: We make the Ansatz (in the sense of formal power series) σ(x, y) = n 1 n! ς µ 1 ···µ n (x)δx µ 1 δx µ n .