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

Abstract

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.

1 Introduction

The semiclassical Einstein equation (SCE) is the equation

$$\begin{aligned} G_{\mu \nu } = \kappa {\langle T^\mathrm {ren}_{\mu \nu } \rangle }_\omega , \end{aligned}$$

(1.1)

where $G_{\mu \nu }$ denotes the Einstein tensor for the (classical) metric $g_{\mu \nu }$, $\kappa$ is a gravitational coupling constant constant (conventionally, $\kappa = 8\uppi \mathrm {G}$), and $\langle T^\mathrm {ren}_{\mu \nu } \rangle _\omega$ is the renormalized stress-energy tensor for a quantum field theory (QFT) in the state $\omega$. 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 wide-spread 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

$$\begin{aligned} \langle T^\mathrm {ren} \rangle= & {} g^{\mu \nu } \langle T_{\mu \nu }^\mathrm {ren} \rangle = c_1 R^2 + c_2 R_{\mu \nu } R^{\mu \nu } + c_3 R_{\mu \nu \lambda \sigma } R^{\mu \nu \lambda \sigma } \\&+ c_4 \Box R + \text {renormalization freedom}, \end{aligned}$$

where the constants $c_\bullet$ 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 \times {{\mathbb {R}}}^3$, $I \subset {{\mathbb {R}}}$, with the metric (in conformal time $\tau$ and with signature convention ${-}{+}{+}{+}$)

$$\begin{aligned} g = a(\tau )^2 \bigl (-\mathop {}\!\mathrm {d}\tau ^2 + \mathop {}\!\mathrm {d}\varvec{x}^2\bigr ), \end{aligned}$$

(1.2)

where $a(\tau ) > 0$ is called the scale factor and $\mathop {}\!\mathrm {d}\varvec{x}^2$ denotes the Euclidean metric on ${{\mathbb {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 includes: [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.

1.1 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 point splitting.

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.¹ 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 Theorems 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 appendix, we list several auxiliary results (concerning homogeneous distributions, weighted sequence spaces, combinatorial inequalities and Synge’s world function) used throughout this work.

2 Scalar Field on Flat FLRW

2.1 Klein–Gordon Equation

Consider the (homogeneous) Klein–Gordon equation

$$\begin{aligned} K \phi :=(\Box + \xi R + m^2) \phi = 0 \end{aligned}$$

(2.1)

with mass $m \ge 0$ and curvature coupling $\xi$ ($\xi =0$ is called minimal coupling and $\xi =\frac{1}{6}$ conformal coupling). We define the d’Alembert operator as $\Box :=-g^{\mu \nu } \nabla _{\!\mu } \nabla _{\!\nu }$, where $\nabla$ is the covariant derivative and we employ Einstein summation convention.

In conformal time $\tau$, the flat FLRW metric takes the form (1.2). The d’Alembertian for this metric is $\Box = a^{-3} (\partial _\tau ^2 - a^{-1} a'' - \Delta ) a$, where we denote derivatives with respect to conformal time by primes and $\Delta$ is the Laplace operator on ${{\mathbb {R}}}^3$. Introducing the conformally rescaled field $\varphi$ and the potential V,

$$\begin{aligned} \varphi :=a \phi , \quad V :=(6\xi -1) \frac{a''}{a} + a^2 m^2, \end{aligned}$$

the Klein–Gordon equation (2.1) can thus be written as $(\partial _\tau ^2 - \Delta + V) \varphi = 0$, where we used $R = 6 a^{-3} a''$. Further rewritten in Hamiltonian form, this equation becomes

$$\begin{aligned} \partial _\tau \begin{pmatrix} \varphi \\ \pi \end{pmatrix} = \begin{pmatrix} 0 &{}\quad 1 \\ \Delta - V &{}\quad 0 \end{pmatrix} \begin{pmatrix} \varphi \\ \pi \end{pmatrix}, \quad \pi :=\varphi '. \end{aligned}$$

(2.2)

2.2 Quantization

For the quantization of the (real, free) scalar field $\phi$ on a globally hyperbolic spacetime (M, g), we follow the algebraic approach, see, e.g., [17, 27, 38]: We generate a (non-commutative, unital) $*$-algebra $\mathcal {A}$ by ‘smeared’ quantum fields $\hat{\phi }(f)$ for $f \in C^\infty _{\mathrm c}(M)$ satisfying

• (i) linearity: $\hat{\phi }(\alpha f + \beta g) = \alpha \hat{\phi }(f) + \beta \hat{\phi }(g)$,

• (ii) hermiticity: $\hat{\phi }(f)^* = \hat{\phi }(\overline{f})$,

• (iii) Klein–Gordon equation: $\hat{\phi }(Kf) = 0$,

• (iv) canonical commutation relations (CCR): $[\hat{\phi }(f),\hat{\phi }(g)] = -\mathrm {i}\langle f \,|\, G^\mathrm {PJ} g \rangle \mathbb {1}$,

where $f,g \in C^\infty _{\mathrm c}(M)$ and $\alpha ,\beta \in {{\mathbb {C}}}$. Moreover, $\langle \cdot \,|\,\cdot \rangle$ is the canonical $L^2$ product on the spacetime, and $G^\mathrm {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 $\mathcal {A}$ is a linear functional $\omega : \mathcal {A} \rightarrow {{\mathbb {C}}}$, which is

• (i) normalized [$\omega (\mathbb {1}) = 1$] and

• (ii) positive [$\omega (a^* a) \ge 0$ for all $a \in \mathcal {A}$],

The two-point function $\omega _2$ of $\omega$ is defined as $\omega _2(f, g) :=\omega (\hat{\phi }(f)\hat{\phi }(g))$. Due to the positivity of the state we have $\omega _2(\overline{f}, f) \ge 0$. Meanwhile, the canonical commutation relations imply

$$\begin{aligned} \omega _2(f,g)-\omega _2(g,f) = -\mathrm {i}\langle f \,|\, G^\mathrm {PJ} g \rangle . \end{aligned}$$

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 $\omega$ has a vanishing one-point function $\omega (\hat{\phi }(f))=0$. Thus, we do not distinguish between the two-point function and the connected two-point function. A non-vanishing one-point function $\phi ^\mathrm {bg}(\tau ,\varvec{x}) :=\omega (\hat{\phi }(\tau ,\varvec{x}))$ is interpreted as a classical Klein–Gordon ‘background’ field. In the context of homogeneous and isotropic spacetimes, $\phi ^\mathrm {bg}(\tau )$ does not depend on $\varvec{x}$. Setting $\varphi ^\mathrm {bg} :=a\phi ^\mathrm {bg}$ and $\pi ^\mathrm {bg} = \partial _\tau \varphi ^\mathrm {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 $\Delta$ can be omitted.

2.3 Two-Point Functions

Consider the two-point function $\omega _2$ of a homogeneous and isotropic state. That is, it holds that

$$\begin{aligned} \omega _2\bigl ((\tau ,\varvec{x}),(\eta ,\varvec{y})\bigr ) = \omega _2\bigl (\tau , \eta , r = |\varvec{x}-\varvec{y}|\bigr ). \end{aligned}$$

(2.3)

Since the antisymmetric part of the two-point function is fixed by the commutator ‘function’ (viz. the Pauli–Jordan propagator $G^\mathrm {PJ}$), $\omega _2$ is completely determined by its symmetric part

$$\begin{aligned} \frac{1}{2} \bigl ( \omega _2(\tau , \eta , r) + \omega _2(\eta , \tau , r) \bigr ). \end{aligned}$$

(2.4)

Therefore, the Cauchy data at conformal time $\tau$ of a homogeneous and isotropic state can be given by (2.4) and its first time derivatives.

We define

$$\begin{aligned} \mathcal {G}(\tau ,r) :=\begin{pmatrix} \mathcal {G}_{\varphi \varphi }(\tau ,r) \\ \mathcal {G}_{(\varphi \pi )}(\tau ,r) \\ \mathcal {G}_{\pi \pi }(\tau ,r) \end{pmatrix} :=\lim _{\eta \rightarrow \tau } \begin{pmatrix} \mathbb {1}\\ \frac{1}{2} (\partial _\tau + \partial _\eta ) \\ \partial _\tau \partial _{\eta } \end{pmatrix} a(\tau ) a(\eta ) \omega _2(\tau ,\eta ,r), \end{aligned}$$

(2.5)

which represents the Cauchy data of the two-point function at $\tau$. 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

$$\begin{aligned} \partial _\tau \mathcal {G} = \begin{pmatrix} 0 &{}\quad 2 &{}\quad 0 \\ \Delta _r - V &{}\quad 0 &{}\quad 1 \\ 0 &{}\quad 2(\Delta _r - V) &{}\quad 0 \end{pmatrix} \mathcal {G}, \end{aligned}$$

(2.6)

where $\Delta _r :=r^{-2} \partial _r r^2 \partial _r$ is the (three dimensional) radial Laplacian.

It is sometimes convenient to perform calculations in momentum space. We define the mode functions

$$\begin{aligned} \widehat{\mathcal {G}}(\tau ,k) :=\begin{pmatrix} \widehat{\mathcal {G}}_{\varphi \varphi }(\tau ,k) \\ \widehat{\mathcal {G}}_{(\varphi \pi )}(\tau ,k) \\ \widehat{\mathcal {G}}_{\pi \pi }(\tau ,k) \end{pmatrix} :=4\uppi \int _0^\infty \mathcal {G}(\tau ,r) \frac{\sin (k r)}{k r} r^2 \mathop {}\!\mathrm {d}r \end{aligned}$$

(2.7)

with $k \in [0,\infty )$. The mode functions are simply the Fourier transform of $\mathcal {G}(\tau ,r=|\varvec{x}|)$ in $\varvec{x}$. Indeed, using the convention $\hat{f}(\varvec{k}) :=\int _{{{\mathbb {R}}}^3} f(\varvec{x}) \mathrm {e}^{-\mathrm {i}\varvec{k}\cdot \varvec{x}} \mathop {}\!\mathrm {d}\varvec{x}$ for the Fourier transform, we have for radial functions (or distributions) $f(\varvec{x}) = f(r = |\varvec{x}|)$

$$\begin{aligned} \hat{f}(\varvec{k}) = \hat{f}(k = |\varvec{k}|) = 4\uppi \int _0^\infty f(r) \frac{\sin (kr)}{kr} r^2 \mathop {}\!\mathrm {d}r. \end{aligned}$$

Conversely, the mode functions specify the two-point function by

$$\begin{aligned} \mathcal {G}(\tau ,r) = \frac{1}{2\uppi ^2} \int _0^\infty \widehat{\mathcal {G}}(\tau ,k) \frac{\sin (k r)}{k r} k^2 \mathop {}\!\mathrm {d}k. \end{aligned}$$

(2.8)

It follows immediately from (2.6) that the mode functions solve the dynamical equations

$$\begin{aligned} \partial _\tau \widehat{\mathcal {G}} = \begin{pmatrix} 0 &{}\quad 2 &{}\quad 0 \\ -(k^2 + V) &{}\quad 0 &{}\quad 1 \\ 0 &{}\quad -2(k^2 + V) &{}\quad 0 \end{pmatrix} \widehat{\mathcal {G}}. \end{aligned}$$

(2.9)

A straightforward computation shows that $\widehat{\mathcal {J}} :=\widehat{\mathcal {G}}_{\varphi \varphi } \widehat{\mathcal {G}}_{\pi \pi } - \widehat{\mathcal {G}}_{(\varphi \pi )}^{\;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 $\widehat{\mathcal {J}} \ge \frac{1}{4}$ with equality for pure states [39].

2.4 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 $\mathcal {G}(\tau ,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 $\mathcal {H}_n(\tau ,r)$ such that $\mathcal {G}(\tau ,r) - \mathcal {H}_n(\tau ,r)$ is sufficiently regular in the limit $r \rightarrow 0$.

Let $\mu$ be an arbitrary (but fixed) length scale. On ${{\mathbb {R}}}$, define the piecewise function

$$\begin{aligned} k^z_+ :={\left\{ \begin{array}{ll} k^z &{} \hbox {if } k > 0, \\ 0 &{} \hbox {if } k \le 0. \end{array}\right. } \end{aligned}$$

For $z \in {{\mathbb {C}}}\setminus \{-1,-2,\dotsc \}$ this defines a distribution (by analytic continuation, cf. Chap. III.2 of [32]). It can be extended to all $z \in {{\mathbb {C}}}$ by defining for $n \in {{\mathbb {N}}}$,

$$\begin{aligned} \langle k_+^{-n}, f \rangle :=\frac{1}{(n-1)!} \biggl ( -\int _0^\infty \log (\mu k) f^{(n)}(k) \mathop {}\!\mathrm {d}k + f^{(n-1)}(0) \sum _{j=1}^{n-1} \frac{1}{j} \biggr ), \quad f \in C^\infty _{\mathrm c}({{\mathbb {R}}}). \end{aligned}$$

We also define for $r \ge 0$ and $z \in {{\mathbb {C}}}$ the distributions

$$\begin{aligned} h_z(r) :=\frac{\mathrm {e}^{\mathrm {i}z \uppi /2}}{2\uppi ^2} \frac{r^{z-2}}{\Gamma (z)} \biggl ( \log \Bigl (\frac{r}{\mu }\Bigr ) - \psi (z) \biggr ), \end{aligned}$$

where $\psi$ denotes the Digamma function. Note that $h_{-2} = -{(\uppi ^2r^4)}^{-1}$ and $h_{0} = 1/{(2\uppi ^2r^2)}^{-1}$. Further details about these functions are found in appendix (Appendix A.1), although without the constant length scale $\mu$. Tacitly, the so-defined distributions and the resulting regularization scheme depend on $\mu$.

We make the Ansatz (equivalently either in position or momentum space with relations analogous to (2.7) and (2.8))

$$\begin{aligned} \mathcal {H}_n(\tau ,r)&:=\left( \begin{array}{c} \mathcal {H}_{\varphi \varphi ,n}(\tau ,r) \\ \mathcal {H}_{(\varphi \pi ),n}(\tau ,r) \\ \mathcal {H}_{\pi \pi ,n}(\tau ,r) \end{array}\right) = \left( \begin{array}{c} 0 \\ 0 \\ \gamma _{-1}(\tau ) \end{array}\right) h_{-2}(r) + \sum _{j=0}^n \left( \begin{array}{c} \alpha _j(\tau ) \\ \beta _j(\tau ) \\ \gamma _j(\tau ) \end{array}\right) h_{2j}(r), \end{aligned}$$

(2.10a)

$$\begin{aligned} \widehat{\mathcal {H}}_n(\tau ,k)&:=\left( \begin{array}{c} \widehat{\mathcal {H}}_{\varphi \varphi ,n}(\tau ,k) \\ \widehat{\mathcal {H}}_{(\varphi \pi ),n}(\tau ,k) \\ \widehat{\mathcal {H}}_{\pi \pi ,n}(\tau ,k) \end{array}\right) = \left( \begin{array}{c} 0 \\ 0 \\ \gamma _{-1}(\tau ) \end{array}\right) k_+^1 + \sum _{j=0}^n \left( \begin{array}{c} \alpha _j(\tau ) \\ \beta _j(\tau ) \\ \gamma _j(\tau ) \end{array}\right) k_+^{-(2j+1)}, \end{aligned}$$

(2.10b)

where we fix the lowest order parameters to

$$\begin{aligned} \alpha _0 = \frac{1}{2}, \quad \beta _0 = 0, \quad \gamma _{-1} = \frac{1}{2}. \end{aligned}$$

(2.11)

We wish to find coefficients $\alpha _\bullet$, $\beta _\bullet$ and $\gamma _\bullet$ such that the two constraints

$$\begin{aligned} \partial _\tau \mathcal {H}_n - \begin{pmatrix} 0 &{}\quad 2 &{}\quad 0 \\ \Delta _r - V &{}\quad 0 &{}\quad 1 \\ 0 &{}\quad 2(\Delta _r - V) &{}\quad 0 \end{pmatrix} \mathcal {H}_n&= \mathcal {O}(r^{2(n-1)}) \end{aligned}$$

(2.12)

$$\begin{aligned} \widehat{\mathcal {H}}_{\varphi \varphi ,n} \widehat{\mathcal {H}}_{\pi \pi ,n} - \widehat{\mathcal {H}}_{(\varphi \pi ),n}^{\,2}&= \frac{1}{4} + \mathcal {O}(k^{-2(n+1)}) \end{aligned}$$

(2.13)

are satisfied. That is, $\mathcal {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

$$\begin{aligned}&\partial _\tau \mathcal {H}_n - \begin{pmatrix} 0 &{}\quad 2 &{}\quad 0 \\ \Delta _r - V &{}\quad 0 &{}\quad 1 \\ 0 &{}\quad 2(\Delta _r - V) &{}\quad 0 \end{pmatrix} \mathcal {H}_n \\&\quad = \sum _{j=0}^{n-1} \begin{pmatrix} \alpha _j' - 2 \beta _j \\ \beta _j' + \alpha _{j+1} + V \alpha _j - \gamma _j \\ \gamma _j' + 2 \beta _{j+1} + 2 V \beta _j \end{pmatrix} h_{2j} + \begin{pmatrix} \alpha _n' - 2 \beta _n, \\ \beta _n' + V \alpha _n - \gamma _n, \\ \gamma _n' + 2 V \beta _n \end{pmatrix} h_{2n}. \end{aligned}$$

Consequently, the coefficients must satisfy the equations

$$\begin{aligned} \alpha _j'&= 2 \beta _j, \end{aligned}$$

(2.14a)

$$\begin{aligned} \beta _j'&= -\alpha _{j+1} - V \alpha _j + \gamma _j, \end{aligned}$$

(2.14b)

$$\begin{aligned} \gamma _j'&= -2 \beta _{j+1} - 2 V \beta _j. \end{aligned}$$

(2.14c)

Moreover, the left-hand side of (2.13) evaluates to

$$\begin{aligned} \widehat{\mathcal {H}}_{\varphi \varphi ,n} \widehat{\mathcal {H}}_{\pi \pi ,n} - \widehat{\mathcal {H}}_{(\varphi \pi ),n}^{\,2}&= \alpha _0 \gamma _{-1} + \sum _{j=1}^n \Bigl ( \alpha _0 \gamma _j + \alpha _{j+1} \gamma _{-1} - \beta _0 \beta _j \\&\quad + \sum _{i=1}^j (\alpha _i \gamma _{j-i} - \beta _i \beta _{j-i}) \Bigr ) k_+^{-2j} + \mathcal {O}\bigl (k_+^{-2(n+1)}\bigr ) \\&= \frac{1}{4} + \frac{1}{2} \sum _{l=1}^n \Bigl ( \alpha _{j+1} + \gamma _j + 2 \sum _{i=1}^j (\alpha _i \gamma _{j-i} - \beta _i \beta _{j-i}) \Bigr ) k_+^{-2j}\\ {}&\quad + \mathcal {O}\bigl (k_+^{-2(n+1)}\bigr ). \end{aligned}$$

Therefore, the coefficients must additionally satisfy the constraint

$$\begin{aligned} \alpha _{j+1} + \gamma _j = -2 \sum _{i=1}^j (\alpha _i \gamma _{j-i} - \beta _i \beta _{j-i}). \end{aligned}$$

(2.15)

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

$$\begin{aligned} \alpha _{j+1}&= -\frac{1}{2} (V \alpha _j + \beta _j') - \sum _{i=1}^j (\alpha _i \gamma _{j-i} - \beta _i \beta _{j-i}), \end{aligned}$$

(2.16a)

$$\begin{aligned} \beta _{j+1}&= -\frac{1}{4} V' \alpha _j - \frac{1}{4} \beta _j'' - V \beta _j, \end{aligned}$$

(2.16b)

$$\begin{aligned} \gamma _j&= \frac{1}{2} (V \alpha _j + \beta _j') - \sum _{i=1}^j (\alpha _i \gamma _{j-i} - \beta _i \beta _{j-i}). \end{aligned}$$

(2.16c)

Proof

To find (2.16b), we eliminate $\gamma _j'$ from (2.14c) by adding the derivative of (2.14b) and using (2.14a). The other two equations obtained from (2.15) and $-\alpha _{j+1} + \gamma _j = V \alpha _j + \beta _j'$, which follows from (2.14b), thus yield (2.16a) and (2.16c).

To see that (2.15) is consistent with (2.14), we first subtract (2.14a) and (2.14c) to obtain $\alpha _{j+1}' + \gamma _j' = -2 V \beta _j$. Then, we calculate

$$\begin{aligned} \partial _\tau \sum _{i=1}^j (\alpha _i \gamma _{j-i} - \beta _i \beta _{j-i})&= \sum _{i=1}^j \bigl ( \alpha _i' \gamma _{j-i} + \alpha _i \gamma _{j-i}' - 2 \beta _i' \beta _{j-i} \bigr ) \\&= 2 \sum _{i=1}^j \bigl ( \beta _i \gamma _{j-i} - \alpha _i (\beta _{j-i+1} + V \beta _{j-i})\\&\quad - (\gamma _i - \alpha _{i+1} - V \alpha _i) \beta _{j-i} \bigr ) \\&= 2 \sum _{i=1}^j (\alpha _i \beta _{j-i+1} - \alpha _{i+1} \beta _{j-i}) + 2 \gamma _0 \beta _j \\&= 2 (\gamma _0 - \alpha _1) \beta _j = V \beta _j, \end{aligned}$$

where in the last step, we used that (2.14b) for $j=0$ implies $\gamma _0 - \alpha _1 = \frac{1}{2} V$.$\square$

2.5 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 $\lambda$) for the Klein–Gordon equation (2.1) is defined as [19, 41]

$$\begin{aligned} H_n(x,y) :=\lim _{\varepsilon \rightarrow 0^+} \frac{1}{8\uppi ^2} \biggl ( \frac{\Delta (x,y)^{\frac{1}{2}}}{\sigma _\varepsilon (x,y)} + \sum _{j=0}^{n-1} v_j(x,y) \sigma (x,y)^j \log \frac{\sigma _\varepsilon (x,y)}{\lambda ^2} \biggr ), \end{aligned}$$

where x, y are in a geodesically convex neighborhood,

$$\begin{aligned} \sigma _{\varepsilon }(x,y) :=\sigma (x,y) + \mathrm {i}\varepsilon \, \bigl (t(x) - t(y)\bigr ) + \frac{1}{2} \varepsilon ^2, \quad \varepsilon > 0, \end{aligned}$$

is the regularized Synge world function (i.e., half the signed squared geodesic distance) for a time function t, $\Delta (x,y)$ is the van Vleck–Morette determinant, $v_j(x,y)$ are smooth coefficient functions satisfying the recurrence relations [19, 41]

$$\begin{aligned} (K \otimes \mathbb {1}) \Delta ^\frac{1}{2}&= \bigl ( -(\Box \otimes \mathbb {1}) \sigma + \sigma (\Box \otimes \mathbb {1}) - 2 \bigr ) v_0, \end{aligned}$$

(2.17a)

$$\begin{aligned} (K \otimes \mathbb {1}) v_{j-1}&= \bigl ( -(\Box \otimes \mathbb {1}) \sigma + \sigma (\Box \otimes \mathbb {1}) + 2j - 2 \bigr ) j v_j \end{aligned}$$

(2.17b)

for $j \in {{\mathbb {N}}}$. These relations are obtained by demanding that the Hadamard parametrix satisfies the Klein–Gordon equation up to an error of order $\sigma ^{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

$$\begin{aligned}{}[v_1]&= \frac{1}{8} m^4 + \frac{6\xi -1}{24} m^2 R + \frac{(6\xi -1)^2}{288} R^2 \nonumber \\&\quad - \frac{1}{720} R_{\mu \nu } R^{\mu \nu } + \frac{1}{720} R_{\mu \nu \rho \sigma } R^{\mu \nu \rho \sigma } + \frac{5\xi -1}{120} \Box R. \end{aligned}$$

(2.18)

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

$$\begin{aligned} H_n\bigl ((\tau ,\varvec{x}), (\eta ,\varvec{y})\bigr ) = H_n\bigl (\tau , \eta , r=|\varvec{x}-\varvec{y}|\bigr ), \end{aligned}$$

analogously to (2.3). (Naturally, the same holds true for $\sigma$, $\Delta ^\frac{1}{2}$ and $v_j$.) As before for the two-point function, this implies that the Hadamard parametrix can be represented by

$$\begin{aligned} \tilde{\mathcal {H}}_n(\tau ,r) :=\begin{pmatrix} \tilde{\mathcal {H}}_{\varphi \varphi ,n}(\tau ,r) \\ \tilde{\mathcal {H}}_{(\varphi \pi ),n}(\tau ,r) \\ \tilde{\mathcal {H}}_{\pi \pi ,n}(\tau ,r) \end{pmatrix} :=\lim _{\eta \rightarrow \tau } \begin{pmatrix} \mathbb {1}\\ \frac{1}{2} (\partial _\tau + \partial _\eta ) \\ \partial _\tau \partial _\eta \end{pmatrix} a(\tau ) a(\eta ) H_n(\tau ,\eta ,r). \end{aligned}$$

It follows from results in [18] that the most singular terms of $\tilde{\mathcal {H}}_n$ and $\mathcal {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 of $\tilde{\mathcal {H}}$ and $\mathcal {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 \ge 2$, one can compute the expansions

$$\begin{aligned} \tilde{\mathcal {H}}_{\varphi \varphi ,n}(\tau ,r)= & {} \frac{1}{4\uppi ^2 r^2} + \frac{1}{48\uppi ^2} \frac{a''(\tau )}{a(\tau )} + \left( \frac{V(\tau )}{16\uppi ^2} + \frac{3V(\tau )^2+V''(\tau )}{384\uppi ^2} r^2 \right) \\&\quad \times \log \left( \frac{r^2 a(\tau )^2}{2\lambda ^2}\right) + \frac{r^2}{5760\uppi ^2} \left( 30 V(\tau ) \frac{a'(\tau )}{a(\tau )} + 11 \Bigl (\partial _\tau \frac{a'(\tau )}{a(\tau )}\Bigr )^2\right. \\&\quad \left. - 2 \frac{a''(\tau )^2}{a(\tau )^2} + 12 \frac{a'(\tau )}{a(\tau )} \partial _\tau \frac{a''(\tau )}{a(\tau )} \right) + \mathcal {O}(r^4),\\ \tilde{\mathcal {H}}_{(\varphi \pi ),n}(\tau ,r)= & {} \frac{V'(\tau )}{32\uppi ^2} \log \left( \frac{r^2 a(\tau )^2}{2\lambda ^2}\right) + \frac{1}{96\uppi ^2} \\&\quad \times \left( 6 V(\tau ) \frac{a'(\tau )}{a(\tau )} + \partial _\tau \frac{a''(\tau )}{a(\tau )} \right) + \mathcal {O}(r^2),\\ \tilde{\mathcal {H}}_{\pi \pi ,n}(\tau ,r)= & {} -\frac{1}{2\uppi ^2 r^4} + \frac{V(\tau )}{8\uppi ^2 r^2} + \frac{V(\tau )^2+V''(\tau )}{64\uppi ^2} \log \left( \frac{r^2 a(\tau )^2}{2\lambda ^2}\right) \\&\quad + \frac{1}{192\uppi ^2} \left( 3 V(\tau )^2 - 6 V(\tau ) \frac{a'(\tau )^2}{a(\tau )^2} + 4 V(\tau ) \frac{a''(\tau )}{a(\tau )}\right. \\&\quad \left. + 12 V'(\tau ) \frac{a'(\tau )}{a(\tau )} + V''(\tau ) \right) + \frac{1}{960\uppi ^2} \\&\quad \times \left( \Bigl (\partial _\tau \frac{a'(\tau )}{a(\tau )}\Bigr )^2 + 2 \frac{a''(\tau )^2}{a(\tau )^2} + 4 \partial _\tau ^2 \frac{a''(\tau )}{a(\tau )} \right) + \mathcal {O}(r^2). \end{aligned}$$

To obtain this result, one uses a covariant expansion of the coefficient functions $\Delta ^\frac{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

$$\begin{aligned} \mu ^2 = 2\mathrm {e}^{2\gamma -2} \lambda ^2 \lambda _0^2, \end{aligned}$$

(2.19)

where $\gamma$ is the Euler–Mascheroni constant and $\lambda _0 > 0$ is an arbitrary dimensionless constant. With this convention, we find the differences ($n \ge 2$)

$$\begin{aligned} \tilde{\mathcal {H}}_{\varphi \varphi ,n} - \mathcal {H}_{\varphi \varphi ,n} \Bigl |_{r=0}&= \frac{V}{8\uppi ^2} \log (a \lambda _0) + \frac{1}{48\uppi ^2} \frac{a''}{a}, \end{aligned}$$

(2.20a)

$$\begin{aligned} \tilde{\mathcal {H}}_{(\varphi \pi ),n} - \mathcal {H}_{(\varphi \pi ),n} \Bigl |_{r=0}&= \frac{V'}{16\uppi ^2} \log (a \lambda _0) + \frac{1}{96\uppi ^2} \left( 6 V \frac{a'}{a} + \partial _\tau \frac{a''}{a} \right) , \end{aligned}$$

(2.20b)

$$\begin{aligned} \tilde{\mathcal {H}}_{\pi \pi ,n} - \mathcal {H}_{\pi \pi ,n} \Bigl |_{r=0}&= \frac{V^2+V''}{32\uppi ^2} \log (a \lambda _0) + \frac{1}{960\uppi ^2} \left( \Bigl (\partial _\tau \frac{a'}{a}\Bigr )^2 + 2 \frac{a^{\prime \prime \,2}}{a^2} + 4 \partial _\tau ^2 \frac{a''}{a} \right) \nonumber \\&\quad + \frac{1}{192\uppi ^2} \left( 3 V^2 - 6 V \frac{a^{\prime \,2}}{a^2} + 4 V \frac{a''}{a} + 12 V' \frac{a'}{a} + V'' \right) , \end{aligned}$$

(2.20c)

$$\begin{aligned} \Delta _r (\tilde{\mathcal {H}}_{\varphi \varphi ,n} - \mathcal {H}_{\varphi \varphi ,n}) \Bigl |_{r=0}&= \frac{3V^2+V''}{32\uppi ^2} \left( \frac{5}{6} + \log (a \lambda _0) \right) + \frac{V}{32\uppi ^2} \frac{a^{\prime \,2}}{a^2} \nonumber \\&\quad + \frac{1}{960\uppi ^2} \left( 11 \Bigl (\partial _\tau \frac{a'}{a}\Bigr )^2 - 2 \frac{a^{\prime \prime \,2}}{a^2} + 12 \frac{a'}{a} \partial _\tau \frac{a''}{a} \right) . \end{aligned}$$

(2.20d)

We will use these differences in the following section.

3 Semiclassical Einstein equation

On FLRW spacetimes, it is sufficient to look at the traced SCE

$$\begin{aligned} -R = \kappa {\langle T^\mathrm {ren} \rangle }_\omega , \end{aligned}$$

(3.1)

where $\langle T^\mathrm {ren} \rangle _\omega :=g^{\mu \nu } \langle T_{\mu \nu }^\mathrm {ren} \rangle _\omega$ is the trace of the quantum stress-energy tensor and at the energy component (viz. the 00- or tt-component)

$$\begin{aligned} G_{00} = \kappa \langle T_{00}^\mathrm {ren} \rangle _\omega \end{aligned}$$

(3.2)

The latter equation can be regarded as a constraint equation, which, if imposed at a fixed time, holds everywhere because of the covariant conservation $\nabla ^\mu \langle T_{\mu \nu }^\mathrm {ren} \rangle _\omega = 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.

3.1 Renormalized Stress-Energy Tensor

The (classical) stress-energy tensor for a (real) free scalar field is

$$\begin{aligned} \begin{aligned} T_{\mu \nu }&= (1-2\xi ) (\nabla _{\!\mu } \phi ) (\nabla _{\!\nu } \phi ) - \frac{1}{2} (1-4\xi ) g_{\mu \nu } (\nabla ^{\sigma } \phi ) (\nabla _{\!\sigma } \phi ) - \frac{1}{2} g_{\mu \nu } m^2 \phi ^2 \\ {}&\quad + \xi \bigr (G_{\mu \nu } \phi ^2 - 2 \phi \nabla _{\!\mu } \nabla _{\!\nu } \phi - 2 g_{\mu \nu } \phi \Box \phi \bigr ). \end{aligned}\end{aligned}$$

(3.3)

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 $\omega$ (cf. [31, 41, 53]):

$$\begin{aligned} \begin{aligned} {\langle T_{\mu \nu }^\mathrm {ren} \rangle }_\omega&= (1-2\xi ) [(\nabla _{\!\mu } \otimes \nabla _{\!\nu }) \omega _2^\mathrm {reg}] - \frac{1}{2} (1-4\xi ) g_{\mu \nu } [(\nabla ^\sigma \otimes \nabla _{\!\sigma }) \omega _2^\mathrm {reg}] - \frac{1}{2} g_{\mu \nu } m^2 [\omega _2^\mathrm {reg}] \\ {}&\quad + \xi \bigl ( G_{\mu \nu } [\omega _2^\mathrm {reg}] - 2 [(\mathbb {1}\otimes \nabla _{\!\mu } \nabla _{\!\nu }) \omega _2^\mathrm {reg}] - 2 g_{\mu \nu } [(\mathbb {1}\otimes \Box ) \omega _2^\mathrm {reg}] \bigr ) \\ {}&\quad + \frac{1}{4\uppi ^2} g_{\mu \nu } [v_1] + c_1 m^4 g_{\mu \nu } + c_2 m^2 G_{\mu \nu } + c_3 I_{\mu \nu } + c_4 J_{\mu \nu }, \end{aligned}\end{aligned}$$

(3.4)

where, for fixed but arbitrary $n \ge 1$, we set $\omega _2^\mathrm {reg} :=\omega _2 - H_n$, $[\,\cdot \,]$ denotes the coincidence limit (e.g., $[v_1](x) = v_1(x,x)$ is the coincidence limit of the Hadamard coefficient $v_1$) with implicit parallel transport, $c_\bullet$ are renormalization parameters, and $I_{\mu \nu }$, $J_{\mu \nu }$ are the two fourth-order conserved curvature tensors:

$$\begin{aligned} I_{\mu \nu }&:=2 R R_{\mu \nu } - 2 \nabla _{\!\mu } \nabla _{\!\nu } R - \frac{1}{2} g_{\mu \nu } (R^2 + 4 \Box R), \\ J_{\mu \nu }&:=2 R^{\rho \sigma } R_{\rho \mu \sigma \nu } - \nabla _{\!\mu } \nabla _{\!\nu } R - \Box R_{\mu \nu } - \frac{1}{2} g_{\mu \nu } (R_{\rho \sigma } R^{\rho \sigma } + \Box R). \end{aligned}$$

These tensors can be obtained as the variations with respect to the metric of $R^2$ and $R_{\mu \nu } R^{\mu \nu }$ [60]. For conformally flat spacetimes like FLRW, it holds $I_{\mu \nu } = 3 J_{\mu \nu }$.

The term involving the Hadamard coefficient $v_1$ in (3.4) ensures that the quantized stress-energy tensor satisfies $\nabla ^\mu \langle T_{\mu \nu }^\mathrm {ren} \rangle _\omega = 0$ [31, 41]. The term $c_1 m^4 g_{a b}$ is a renormalization of the cosmological constant, and $c_2 m^2 G_{a b}$ corresponds to a renormalization of the gravitational coupling constant $\kappa$; the remaining two renormalization terms have no classical interpretation. A change of the scale $\lambda$ in the Hadamard parametrix corresponds to a change of the renormalization parameters [27], and thus, we are in principle at liberty to set $\lambda$ to any value if we change the renormalization parameters accordingly.

We also remark that in case the state $\omega$ includes a nonzero one-point function (classical Klein–Gordon field), the corresponding contribution to the two point function is $\phi ^\mathrm {bg}(x)\phi ^\mathrm {bg}(y)$, that is, in (3.4) we have to replace $\omega _2^\mathrm {reg}(x,y)$ by $\omega _2^\mathrm {reg}(x,y) + \phi ^\mathrm {bg}(x)\phi ^\mathrm {bg}(y)$. After performing the coincidence limit, this leads to an additional contribution given by the classical stress-energy tensor (3.3) with $\phi$ replaced by $\phi ^\mathrm {bg}$.

3.2 Trace of the Stress-Energy Tensor

Taking the trace of (3.4) yields²

$$\begin{aligned} \begin{aligned} {\langle T^\mathrm {ren} \rangle }_\omega&= \bigl ( (6\xi -1) (\xi R + m^2) - m^2 \bigr ) [\omega _2^\mathrm {reg}] + (6\xi -1) g^{\mu \nu } [(\nabla _{\!\mu } \otimes \nabla _{\!\nu }) \omega _2^\mathrm {reg}] \\ {}&\quad - \frac{9\xi -2}{2\uppi ^2} [v_1] + 4 c_1 m^4 - c_2 m^2 R - (6 c_3 + 2 c_4) \Box R. \end{aligned}\end{aligned}$$

(3.5)

where $c_\bullet$ are the same constants as above, and we used that $4 \uppi ^2 [(\mathbb {1}\otimes K) \omega _2^\mathrm {reg}] = 3 [v_1]$.

For homogeneous states and isotropic states (i.e., satisfying (2.3)) on FLRW spacetimes, we thus find

$$\begin{aligned} \begin{aligned} {\langle T^\mathrm {ren} \rangle }_\omega&= \bigl ( (6\xi -1) (\xi R + m^2) - m^2 \bigr ) [\omega _2^\mathrm {reg}] - \frac{6\xi -1}{a^2} \bigl ( [(\Delta \otimes \mathbb {1}) \omega _2^\mathrm {reg}] \\ {}&\quad + [(\partial _\tau \otimes \partial _\tau ) \omega _2^\mathrm {reg}] \bigr ) - \frac{9\xi -2}{2\uppi ^2} [v_1] + 4 c_1 m^4 - c_2 m^2 R - (6 c_3 + 2 c_4) \Box R, \end{aligned}\end{aligned}$$

(3.6)

where

$$\begin{aligned} R = 6 \frac{a''}{a^3}, \quad \Box R = 36 \frac{a'' a^{\prime \,2}}{a^7} - 18 \frac{a^{\prime \prime \,2}}{a^6} - 24 \frac{a^{(3)} a'}{a^6} + 6 \frac{a^{(4)}}{a^5} \end{aligned}$$

and (2.18) specializes to

$$\begin{aligned}\begin{aligned}{}[v_1]&= \frac{m^4}{8} + \frac{1}{60} \left( \frac{a^{\prime \,4}}{a^8} - \frac{a'' a^{\prime \,2}}{a^7} \right) + \frac{(6\xi -1) m^2}{4} \frac{a''}{a^3} + \frac{(6\xi -1)^2}{8} \frac{a^{\prime \prime \,2}}{a^6} \\ {}&\quad + \frac{5\xi -1}{20} \left( 6 \frac{a'' a^{\prime \,2}}{a^7} - 3 \frac{a^{\prime \prime \,2}}{a^6} - 4 \frac{a^{(3)} a'}{a^6} + \frac{a^{(4)}}{a^5} \right) . \end{aligned}\end{aligned}$$

In the case of a non-vanishing one-point function yielding the classical field $\phi ^\mathrm {bg}$, the following expression needs to be added to $\langle T^\mathrm {ren}\rangle _\omega$:

$$\begin{aligned} \bigl ( (6\xi -1) (\xi R + m^2) - m^2 \bigr ) \bigl (\phi ^\mathrm {bg}\bigr )^2 - \frac{6\xi -1}{a^2} \bigl (\dot{\phi }^\mathrm {bg}\bigr )^2. \end{aligned}$$

(3.7)

3.3 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:

$$\begin{aligned} \begin{aligned} {\langle T_{00}^\mathrm {ren} \rangle }_\omega&= \frac{1}{2} [(\partial _\tau \otimes \partial _\tau ) \omega _2^\mathrm {reg}] - \frac{1}{2} [(\mathbb {1}\otimes \Delta ) \omega _2^\mathrm {reg}] + \frac{1}{2} a^2 m^2 [\omega _2^\mathrm {reg}] \\ {}&\quad + \xi \biggl ( G_{00} [\omega _2^\mathrm {reg}] + 6 \frac{a'}{a} [(\mathbb {1}\otimes \partial _\tau ) \omega _2^\mathrm {reg}] \biggr ) \\ {}&\quad - \frac{a^2}{4\uppi ^2} [v_1] - c_1 a^2 m^4 + c_2 m^2 G_{00} + (3c_3 + c_4) J_{00}, \end{aligned}\end{aligned}$$

(3.8)

where

$$\begin{aligned} G_{00} = 3 \frac{a^{\prime \,2}}{a^2}, \quad J_{00} = -24 \frac{a'' a^{\prime \,2}}{a^5} - 6 \frac{a^{\prime \prime \,2}}{a^4} + 12 \frac{a^{(3)} a'}{a^4}. \end{aligned}$$

A straightforward calculation shows

$$\begin{aligned} R = \frac{1}{a^2} \left( \frac{a}{a'} \partial _\tau G_{00} + 2 G_{00} \right) . \end{aligned}$$

(3.9)

Consequently, we also have

$$\begin{aligned} {\langle T^\mathrm {ren} \rangle }_\omega = -\frac{1}{a^2} \biggl ( \frac{a}{a'} \partial _\tau {\langle T_{00}^\mathrm {ren} \rangle }_\omega + 2 {\langle T_{00}^\mathrm {ren} \rangle }_\omega \biggr ), \end{aligned}$$

(3.10)

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

$$\begin{aligned} \partial _\tau \bigl ( G_{00} - \kappa {\langle T_{00}^\mathrm {ren} \rangle }_\omega \bigr ) = -2 \frac{a'}{a} \bigl ( G_{00} - \kappa {\langle T_{00}^\mathrm {ren} \rangle }_\omega \bigr ) \end{aligned}$$

(3.11)

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} - \kappa {\langle T_{00}^\mathrm {ren} \rangle }_\omega$ is asymptotically pushed to zero if $a' > 0$.

If a non-vanishing one-point function is present, it contributes to $\langle T_{00}^\mathrm {ren} \rangle _\omega$ the following expression:

$$\begin{aligned} \frac{1}{2} \bigl (\partial _\tau \phi ^{\mathrm {bg}}\bigr )^2 + \frac{1}{2} a^2 m^2 \bigl (\phi ^\mathrm {bg}\bigr )^2 + \xi \biggl ( G_{00} \bigl (\phi ^\mathrm {bg}\bigr )^2 + 3 \frac{a'}{a} \partial _\tau \bigl (\phi ^\mathrm {bg}\bigr )^2 \biggr ). \end{aligned}$$

3.4 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 = \langle T^\mathrm {ren} \rangle _\omega$ 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 neighborhood 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.

3.5 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³ 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 neighborhood 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 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.

4 Dynamical System for Sequences of Coincidence Limits

4.1 Dynamics

For $n \in {{\mathbb {N}}}_0$, we define

$$\begin{aligned} \mathcal {M}_n :=\Delta _r^n (\mathcal {G} - \mathcal {H}_l) \Big |_{r=0}, \quad l \ge n+1. \end{aligned}$$

(4.1)

That this definition is indeed independent of l follows immediately from (A.4) and the definition of $\mathcal {H}_n$:

Proposition 4.1

For $j \ge l \ge n+1$, we have

$$\begin{aligned} \Delta _r^n (\mathcal {H}_j - \mathcal {H}_l) \Big |_{r=0} = 0. \end{aligned}$$

Occasionally, we call $\mathcal {M}_n$ moments. To understand this nomenclature, consider the momentum space representation of $\mathcal {M}_n$:

$$\begin{aligned} \mathcal {M}_n = (-1)^n \int _0^\infty k^{2(n+1)} (\hat{\mathcal {G}} - \hat{\mathcal {H}}_l) \mathop {}\!\mathrm {d}k. \end{aligned}$$

That is, the sequence of $\mathcal {M}_n$ is given by the moments of $\hat{\mathcal {G}} - \hat{\mathcal {H}}$. Note, however, that $\hat{\mathcal {G}} - \hat{\mathcal {H}}$ cannot be expected to be positive.

To formulate a dynamics for $\mathcal {M}_n$, recall (2.6) and (2.12). It follows that

$$\begin{aligned} \partial _\tau \mathcal {M}_n = \partial _\tau \Delta _r^n (\mathcal {G} - \mathcal {H}_l) \Big |_{r=0} = A \mathcal {M}_n + B \mathcal {M}_{n+1}, \end{aligned}$$

where we defined

Sometimes we write $A(\tau )$ to emphasize the dependence on $\tau$ 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 $\mathcal {M}_n$ as a sequence $\mathcal {M} = (\mathcal {M}_n)$, we can also write the equation above as

$$\begin{aligned} \partial _\tau \mathcal {M}(\tau ) = \bigl (A(\tau ) \otimes \mathbb {1}+ B \otimes L\bigr ) \mathcal {M}(\tau ), \end{aligned}$$

(4.2)

where L denotes the left-shift operator. Below we study this equation on the weighted sequence spaces

$$\begin{aligned} \varvec{\ell }^p(w) :={{\mathbb {R}}}^3 \otimes \ell ^p(w), \end{aligned}$$

(4.3)

where $p \ge 1$ and w is a sequence of weights; we denote the norms by ${\Vert \,\cdot \,\Vert }_{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 $\mathcal {G} - \mathcal {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 $\mathcal {M}$ in two relevant examples to motivate the choice of weights later on. Afterward, 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 \otimes \mathbb {1}+ B \otimes L$.

Finally, we wish to remark that the set of possible sequences $\mathcal {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.

4.2 $\mathcal {M}$ for the Massive Vacuum State on Minkowski Spacetime

First, let us recall the expansions of the modified Bessel functions, see, e.g., Chap. 10 of [43]:

$$\begin{aligned} I_\nu (z)&= \bigl (\tfrac{1}{2} z\bigr )^\nu \sum _{n=0}^\infty \frac{\bigl (\tfrac{1}{4} z^2\bigr )^j}{j!\Gamma (\nu +j+1)}, \\ K_n(z)&= \tfrac{1}{2} \bigl (\tfrac{1}{2} z\bigr )^{-n} \sum _{j=0}^{n-1} \frac{(n-j-1)!}{j!} \bigl (-\tfrac{1}{4} z^2\bigr )^j + (-1)^{n+1} \log \bigl (\tfrac{1}{2} z\bigr ) I_n(z) \\ {}&\quad + (-1)^n \tfrac{1}{2} \bigl (\tfrac{1}{2} z\bigr )^n \sum _{j=0}^\infty \bigl (\psi (j+1) + \psi (n+j+1)\bigr ) \frac{\bigl (\tfrac{1}{4} z^2\bigr )^j}{j!(n+j)!}, \end{aligned}$$

for $\nu \in {{\mathbb {C}}}$ and $n = 0,1,2,\dotsc$, where $\psi$ denotes the Digamma function.

The two-point function for the massive vacuum state on Minkowski spacetime has the initial data:

$$\begin{aligned} \mathcal {G}_{\varphi \varphi }(r)&= \frac{m}{4\uppi ^2r} K_1(mr) \nonumber \\&= \frac{1}{4\uppi ^2r^2} + \frac{m^2}{8\uppi ^2} \sum _{j=0}^\infty \Bigl ( \log (\tfrac{1}{2} mr) - \tfrac{1}{2} \bigl ( \psi (j+1) + \psi (j+2) \bigr )\Bigr ) \frac{\bigl (\tfrac{1}{2} mr\bigr )^{2j}}{j!(j+1)!}, \end{aligned}$$

(4.4a)

$$\begin{aligned} \mathcal {G}_{(\varphi \pi )}(r)&= 0, \end{aligned}$$

(4.4b)

$$\begin{aligned} \mathcal {G}_{\pi \pi }(r)&= -\frac{m^2}{4\uppi ^2r^2} K_2(mr) \nonumber \\&= -\frac{1}{2\uppi ^2r^4} + \frac{m^2}{8\uppi ^2r^2}\nonumber \\&\quad + \frac{m^4}{16\uppi ^2} \sum _{j=0}^\infty \Bigl ( \log (\tfrac{1}{2} mr) - \tfrac{1}{2} \bigl ( \psi (j+1) + \psi (j+3) \bigr )\Bigr ) \frac{\bigl (\tfrac{1}{2} mr\bigr )^{2j}}{j!(j+2)!}, \end{aligned}$$

(4.4c)

Now we determine $\mathcal {H}_n$ on Minkowski spacetime. For this purpose, we need to solve the recurrence relations (2.2) for constant $V = m^2$. Hence, also the coefficients $\alpha _\bullet , \beta _\bullet$ and $\gamma _\bullet$ are constant, and (2.2) imply $\beta _\bullet = 0$,

$$\begin{aligned} \alpha _{j+1} - \gamma _j = -m^2 \alpha _j, \quad \alpha _{j+1} + \gamma _j = -2\sum _{i=1}^j \alpha _i \gamma _{j-i}. \end{aligned}$$

(4.5)

This recurrence relation can be solved in closed form:

Proposition 4.2

If $V = m^2$ is constant,

$$\begin{aligned} \alpha _j = \frac{1}{2} \left( {\begin{array}{c}-\frac{1}{2}\\ j\end{array}}\right) m^{2j}, \quad \beta _j = 0, \quad \gamma _{j-1} = \frac{1}{2} \left( {\begin{array}{c}\frac{1}{2}\\ j\end{array}}\right) m^{2j}. \end{aligned}$$

(4.6)

Proof

Inserting (4.6) on the left-hand side of (4.5) yields

$$\begin{aligned} \alpha _{j+1} - \gamma _j&= \frac{1}{2} \biggl ( \left( {\begin{array}{c}-\frac{1}{2}\\ j+1\end{array}}\right) - \left( {\begin{array}{c}\frac{1}{2}\\ j+1\end{array}}\right) \biggr ) m^{2(j+1)} \\&= -\frac{1}{2} \left( {\begin{array}{c}-\frac{1}{2}\\ j\end{array}}\right) m^{2(j+1)} = -m^2 \alpha _j, \\ \alpha _{j+1} + \gamma _j&= \frac{1}{2} \biggl ( \left( {\begin{array}{c}-\frac{1}{2}\\ j+1\end{array}}\right) + \left( {\begin{array}{c}\frac{1}{2}\\ j+1\end{array}}\right) \biggr ) m^{2(j+1)} \\&= \frac{1}{2} \biggl ( -\sum _{i=0}^{j+1} \left( {\begin{array}{c}-\frac{1}{2}\\ i\end{array}}\right) \left( {\begin{array}{c}\frac{1}{2}\\ j-i+1\end{array}}\right) + \left( {\begin{array}{c}\frac{1}{2}\\ j+1\end{array}}\right) \biggr ) m^{2(j+1)} \\ {}&= -\frac{1}{2} m^{2(j+1)} \sum _{i=1}^j \left( {\begin{array}{c}-\frac{1}{2}\\ i\end{array}}\right) \left( {\begin{array}{c}\frac{1}{2}\\ j-i+1\end{array}}\right) = -2 \sum _{i=1}^j \alpha _i \gamma _{j-i}, \end{aligned}$$

where in the second step of the second equation we used the Chu–Vandermonde identity. $\square$

Then, with the coefficients (4.6), we find

$$\begin{aligned} \mathcal {H}_{\varphi \varphi ,n}(r)&= \frac{1}{4\uppi ^2r^2} + \frac{m^2}{8\uppi ^2} \sum _{j=0}^{n-1} \biggl ( \log \Bigl (\frac{r}{\mu }\Bigr ) - \psi (2j+2) \biggr ) \frac{\bigl (\tfrac{1}{2} mr\bigr )^{2j}}{j!(j+1)!}, \end{aligned}$$

(4.7a)

$$\begin{aligned} \mathcal {H}_{(\varphi \pi ),n}(r)&= 0, \end{aligned}$$

(4.7b)

$$\begin{aligned} \mathcal {H}_{\pi \pi ,n}(r)&= -\frac{1}{2\uppi ^2r^4} + \frac{m^2}{8\uppi ^2r^2} + \frac{m^4}{16\uppi ^2} \sum _{j=0}^{n-1} \biggl ( \log \Bigl (\frac{r}{\mu }\Bigr ) - \psi (2j+2) \biggr ) \frac{\bigl (\tfrac{1}{2} mr\bigr )^{2j}}{j!(j+2)!}. \end{aligned}$$

(4.7c)

Therefore, subtracting (4.7) (of sufficiently high order) from (4.4) and differentiating n times with $\Delta$, we obtain in the coinciding point limit $r \rightarrow 0$:

$$\begin{aligned} \mathcal {M}_{\varphi \varphi ,n}&= \frac{1}{2\uppi ^2} \bigl (\tfrac{1}{2} m\bigr )^{2n+2} \Bigl ( \log (\tfrac{1}{2} m \mu ) + \psi (2n+2) \nonumber \\&\quad - \tfrac{1}{2} \bigl ( \psi (n+1) + \psi (n+2) \bigr )\Bigr ) \left( {\begin{array}{c}2n+1\\ n+1\end{array}}\right) , \end{aligned}$$

(4.8a)

$$\begin{aligned} \mathcal {M}_{(\varphi \pi ),n}&= 0, \end{aligned}$$

(4.8b)

$$\begin{aligned} \mathcal {M}_{\pi \pi ,n}&= \frac{1}{\uppi ^2} \bigl (\tfrac{1}{2} m\bigr )^{2n+4} \Bigl ( \log (\tfrac{1}{2} m \mu ) + \psi (2n+2) \nonumber \\&\quad - \tfrac{1}{2} \bigl ( \psi (n+1) + \psi (n+3) \bigr )\Bigr ) \frac{(2n+1)!}{n!(n+2)!}. \end{aligned}$$

(4.8c)

Let $r \ge \frac{3}{2}$. It follows from the properties of the Digamma function (e.g., (5.4.14) and (5.4.15) of [43]) that

$$\begin{aligned} \psi (2n+2) - \frac{1}{2} \bigl ( \psi (n+1) + \psi (n+r) \bigr ) \le \log (2) \end{aligned}$$

and the bound is approached as $n \rightarrow \infty$. Moreover, by the duplication formula for the Gamma function,

$$\begin{aligned} \frac{(2n+1)}{n!(n+r)!} \le \frac{2^{2n+1}}{\sqrt{\uppi }}. \end{aligned}$$

Consequently, we have

$$\begin{aligned} |\mathcal {M}_{\varphi \varphi ,n}| \le \frac{|\log (m \mu )|}{4\uppi ^{5/2}} m^{2n+2}, \quad |\mathcal {M}_{\pi \pi ,n}| \le \frac{|\log (m \mu )|}{8\uppi ^{5/2}} m^{2n+4}, \end{aligned}$$

whence $\mathcal {M} \in \varvec{\ell }^p(w)$ for $p \ge 1$ and $w_n = \omega ^n$ with $\omega > m^2$.

As a further consistency check, we show that $(A \otimes \mathbb {1}+ B \otimes L) \mathcal {M} = 0$ for $V = m^2$. Indeed, inserting (4.8) into the left-hand side, a straightforward calculation shows that $-m^2 \mathcal {M}_{\varphi \varphi ,n} + \mathcal {M}_{\pi \pi ,n} + \mathcal {M}_{\varphi \varphi ,n+1} = 0$.

Finally, we note that in the massless case $m=0$ all moments vanish exactly: $\mathcal {M}=0$.

4.3 $\mathcal {M}$ for the Massless Thermal State on Minkowski Spacetime

The two-point function for the massless KMS-state on Minkowski spacetime at inverse temperature $\beta$ has the initial data

$$\begin{aligned} \mathcal {G}_{\varphi \varphi }(r)&= \frac{1}{4\uppi r\beta } \coth \Bigl (\frac{\uppi r}{\beta }\Bigr ) = \frac{1}{4\uppi ^2r^2} + \frac{1}{2\uppi ^2\beta ^2} \sum _{j=0}^\infty (-1)^j \zeta (2j+2) \Bigl (\frac{r}{\beta }\Bigr )^{2j}, \end{aligned}$$

(4.9a)

$$\begin{aligned} \mathcal {G}_{(\varphi \pi )}(r)&= 0, \end{aligned}$$

(4.9b)

$$\begin{aligned} \mathcal {G}_{\pi \pi }(r)&= -\frac{\uppi }{2r\beta ^3} \coth \Bigl (\frac{\uppi r}{\beta }\Bigr ) {{\,\mathrm{csch}\,}}\Bigl (\frac{\uppi r}{\beta }\Bigr )^2 \nonumber \\ {}&= -\frac{1}{2\uppi ^2r^4} + \frac{1}{2\uppi ^2\beta ^4} \sum _{j=0}^\infty (-1)^j (2j+2) (2j+3) \zeta (2j+4) \Bigl (\frac{r}{\beta }\Bigr )^{2j}, \end{aligned}$$

(4.9c)

where $\zeta$ denotes the (Riemann) Zeta function. These data can be obtained after an appropriate Bogoliubov transformation from the two-point function of the vacuum state. The coefficients for $\mathcal {H}_n$ with $V=0$ are all zero except for $\alpha _0$ and $\gamma _{-1}$. Therefore,

$$\begin{aligned} \mathcal {H}_{\varphi \varphi ,n}(r) = \frac{1}{4\uppi ^2r^2}, \quad \mathcal {H}_{(\varphi \pi ),n}(r) = 0, \quad \mathcal {H}_{\pi \pi ,n}(r) = -\frac{1}{2\uppi r^4}. \end{aligned}$$

(4.10a)

Subtracting (4.10) from (4.9) and differentiating n times with $\Delta$, we obtain in the coinciding point limit $r \rightarrow 0$:

$$\begin{aligned} \mathcal {M}_{\varphi \varphi ,n}&= (-1)^n \beta ^{-2n-2} \zeta (2n+2) (2n+1)!, \\ \mathcal {M}_{(\varphi \pi ),n}&= 0, \\ \mathcal {M}_{\pi \pi ,n}&= (-1)^n \beta ^{-2n-4} \zeta (2n+4) (2n+3)!. \end{aligned}$$

Note that $\zeta (2) = \frac{1}{6}\uppi ^2$, $\zeta (4) = \frac{1}{90}\uppi ^4$ and the Zeta function $\zeta (s)$ is monotonically decreasing for $s > 1$ with limit 1. Hence, $\mathcal {M} \in \varvec{\ell }^p(w)$ for $p \ge 1$ and $w_n = (2n)! \omega ^{2n}$ with $\omega >\beta ^{-1}$.

As a further consistency check, we show that $(A \otimes \mathbb {1}+ B \otimes L) \mathcal {M} = 0$ for $V=0$. Indeed, it is easily checked that $\mathcal {M}_{\pi \pi ,n} + \mathcal {M}_{\varphi \varphi ,n+1} = 0$.

4.4 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 nonzero 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^2A)^n$). This implies that a nonzero word of length n contains at most $\lceil \frac{2n}{3}\rceil$ 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 nonzero entries of a $3\times 3$ matrix. For instance, denotes any $3\times 3$ matrix $(a_{i\!j})$ whose only nonzero 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 . Another example is . Products and powers of $3\times 3$ matrices will also be represented in this notation. For example, we write for $A(t_1) A(t_2) A(t_3)$. Note that .

A quick calculation shows the following additional relations pertaining to products of the matrices A and B:

(4.11)

Using these relations, we can show that nonzero words in A and B of length n contain at most B-factors. More generally, we have:

Proposition 4.3

A nonzero word in and of length n contains at most -factors.

Proof

By induction, it follows easily from (4.11), that a nonzero 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 nonzero word of length $2n+1$ contains at most $n+1$ factors of the shape . It is then obvious that a nonzero word of length 2n can also not contain more than $n+1$ factor of the shape . This completes the proof. $\square$

We conclude this subsection on the properties of the matrices A and B by noting that their matrix (max) norms are given by:

$$\begin{aligned} \Vert A\Vert = 2 \sqrt{1 + V^2}, \quad \Vert B\Vert = 2. \end{aligned}$$

(4.12)

4.5 Evolution Operator

In this subsection, we define the evolution operator for

$$\begin{aligned} S(\tau ) :=A(\tau ) \otimes \mathbb {1}+ B \otimes L \end{aligned}$$

(4.13)

as an operator on the weighted sequence spaces $\varvec{\ell }^p(w)$ (see (4.3) and Sect. A.2 for the definition of these spaces) for certain weight sequences w and $p \ge 1$.

Recall that evolution operator for $S(\tau )$ is the ‘solution operator’

$$\begin{aligned} U(\tau , \tau _0) \mathcal {M}_0 = \mathcal {M}(\tau ) \end{aligned}$$

associated with the solutions of the initial value problem for the dynamical equation (4.2) with initial data $\mathcal {M}_0$ at $\tau _0$:

$$\begin{aligned} {\left\{ \begin{array}{ll} \partial _\tau \mathcal {M}(\tau ) = S(\tau ) \mathcal {M}(\tau ), \\ \mathcal {M}(\tau _0) = \mathcal {M}_0.\\ \end{array}\right. } \end{aligned}$$

It can formally be obtained via the time-ordered exponential

$$\begin{aligned} U(\tau ,\tau _0) :={\text {Texp}}\left( \int _{\tau _0}^\tau S(\eta ) \mathop {}\!\mathrm {d}\eta \right) . \end{aligned}$$

That is, defining for $n>0$

$$\begin{aligned} U_n(\tau ,\tau _0) :=\int _{\tau _0}^\tau \int _{\tau _0}^{\tau _1} \!\cdots \! \int _{\tau _0}^{\tau _{n-1}} S(\tau _1) \,\cdots \, S(\tau _n) \mathop {}\!\mathrm {d}\tau _n \cdots \mathop {}\!\mathrm {d}\tau _1, \end{aligned}$$

(4.14)

it is given by the Dyson series

$$\begin{aligned} U(\tau ,\tau _0) = {\left\{ \begin{array}{ll} \mathbb {1}+ \sum \nolimits _{n=1}^\infty U_n(\tau ,\tau _0) &{} \hbox {for } \tau \ge \tau _0 \\ \mathbb {1}+ \sum \nolimits _{n=1}^\infty (-1)^n U_n(\tau _0,\tau ) &{} \hbox {for } \tau < \tau _0. \end{array}\right. } \end{aligned}$$

(4.15)

Note that this is simply the solution of the Picard iteration method.

Motivated by the examples in Sects. 4.2 and 4.3, in the following two subsections we show that the series (4.15) converges as an operator

• (W1) on $\varvec{\ell }^{p(w)}$ for $w_n = c \omega ^n$ with $c,\omega > 0$,

• (W2) from $\varvec{\ell }^{p(w)}$ to $\varvec{\ell }^{p(v)}$ for $w_n = (2n)! \omega ^{2n}$ and $v_n = (2n)! \upsilon ^{2n}$ with $\upsilon > \omega \ge 1$ in bounded time intervals,

and that it has the properties of an evolution operator, i.e.,

• (E1) $U(\tau ,\tau ) = \mathbb {1}$,

• (E2) $U(\tau ,\tau _0) = U(\tau ,\tau _1) U(\tau _1,\tau _0)$,

• (E3) $\partial _\tau U(\tau ,\tau _0) = S(\tau ) U(\tau ,\tau _0)$ and $\partial _{\tau _0} U(\tau ,\tau _0) = -U(\tau ,\tau _0) S(\tau _0)$

for $\tau , \tau _0, \tau _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 $\ell ^p(w) \subset \ell ^\infty (w)$ for any $p \in [1,\infty ]$, and hence, we will restrict our attention to the case $p = \infty$. A sequence of moments $(\mathcal {M}_n)$ is in the sequence space $\ell ^\infty (w)$ of factorially growing weights $w_n = \omega ^{2n} (2n)!$ (W2), if there exists a constant $M > 0$ such that

$$\begin{aligned} |\mathcal {M}_n| \le M \omega ^{2n} (2n)! \end{aligned}$$

for all $n \in {{\mathbb {N}}}_0$. Compare this to the following well-known definition: a smooth function f is analytic in $K \Subset {{\mathbb {R}}}$ if and only if there exist constants $C > 0$ and $M > 0$ such that

$$\begin{aligned} \sup _{x \in K}\ |f^{(n)}(x)| \le M C^n n!. \end{aligned}$$

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 neighborhood of the origin. On the other hand, if $\mathcal {G} - \mathcal {H}$ is analytic in a neighborhood 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 = \omega ^{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.

4.6 Evolution Operator for Geometrically Growing Weights

In this subsection, we consider the case (W1), i.e., the evolution operator on $\varvec{\ell }^p(w)$ for the geometrically growing weights $w_n = c \omega ^n$ with $c,\omega > 0$. As described above, in this case $\mathcal {G}-\mathcal {H}$ must be entire analytic in the radial variable. In this case, $S(\tau )$ is a bounded operator on $\varvec{\ell }^p(w)$ with

$$\begin{aligned} {\Vert S(\tau ) \mathcal {M}\Vert }_{p,w} \le (\Vert A(\tau )\Vert + \Vert B\Vert \omega ) {\Vert \mathcal M\Vert }_{p,w} \le 2\bigl (\sqrt{1+V(\tau )^2} + \omega \bigr ) {\Vert \mathcal M\Vert }_{p,w}. \end{aligned}$$

(4.16)

The following theorem is easily shown, see, e.g., Thms. 5.1 and 5.2 in [46] (the proofs of the analogous Theorem 4.10 in the next subsection can also be adapted to the case of geometrically growing weights discussed here):

Theorem 4.4

Let $c,\omega > 0$ and $p \ge 1$. Suppose that $I \subsetneq {{\mathbb {R}}}$ and $V \in C(I)$. The evolution operator $U(\tau ,\tau _0)$, defined by (4.15), is the unique bounded operator on $\varvec{\ell }^p(w)$ with weights $w_n = c \omega ^n$, such that the properties (E1)–(E3) hold for $\tau ,\tau _0,\tau _1 \in I$. Moreover, we have the bound

$$\begin{aligned} {\Vert U(\tau ,\tau _0)\Vert }_{p,w} \le \mathrm {e}^{C_\omega (\tau ,\tau _0)}, \end{aligned}$$

where

$$\begin{aligned} C_\omega (\tau ,\tau _0) :=2 \omega |\tau - \tau _0| + 2 |*|{\int _{\tau _0}^\tau \sqrt{1+V(\eta )^2} \mathop {}\!\mathrm {d}\eta }, \end{aligned}$$

(4.17)

and the evolution operator is norm-continuous:⁴$\lim _{\tau \rightarrow \tau _0} {\Vert U(\tau ,\tau _0) - \mathbb {1}\Vert }_{p,w} = 0$.

Proof

Let us sketch a proof. Using the Banach fixed point theorem and the boundedness of the operators $S(\tau )$, it can be seen that the integral equation

$$\begin{aligned} U(\tau , \tau _0) = \mathbb {1}+ \int _{\tau _0}^\tau S(\eta ) U(\eta , \tau _0) \mathop {}\!\mathrm {d}\eta \end{aligned}$$

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(\tau ,\tau _0)^{-1} = U(\tau _0,\tau )$ and, together with the first part of (E3), we thus find $\partial _{\tau _0} U(\tau ,\tau _0) = -U(\tau ,\tau _0) S(\tau _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. $\square$

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 $\varvec{\ell }^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 Theorem 4.11:

Theorem 4.7

In addition to the assumptions of Theorem 4.4, suppose that $\tilde{V} \in C(I)$. Denote by $U(\tau ,\tau _0)$ and $\tilde{U}(\tau ,\tau _0)$ the evolution operators for V and $\tilde{V}$. If $\mathcal {M}, \tilde{\mathcal M} \in \varvec{\ell }^p(w)$, we have

$$\begin{aligned}&{\Vert U(\tau ,\tau _0) \mathcal {M} - \tilde{U}(\tau ,\tau _0) \tilde{\mathcal M}\Vert }_{p,w} \\&\quad \le \mathrm {e}^{C_\omega (\tau ,\tau _0)} {\Vert \mathcal {M} - \tilde{\mathcal M}\Vert }_{p,w} + 2 \mathrm {e}^{C_\omega (\tau ,\tau _0)+\tilde{C}_\omega (\tau ,\tau _0)} {\Vert \tilde{\mathcal M}\Vert }_{p,w} \int _{\tau _0}^\tau |V(\eta )-\tilde{V}(\eta )| \mathop {}\!\mathrm {d}\eta . \end{aligned}$$

Under more stringent assumptions on the potential V, one obtains the following regularity result:

Proposition 4.8

If, in addition to the assumptions of Theorem 4.4, $V \in C^k(I)$, then $\partial _\tau ^{k+1} U(\tau ,\tau _0)$ is bounded on $\varvec{\ell }^p(w)$.

Proof

By property (E3), $\partial _\tau U(\tau ,\tau _0)$ is bounded on $\varvec{\ell }^p(w)$. Repeated differentiation of (E3), together with the boundedness of $S(\tau )$ and its derivatives, yields the desired result by iteration. $\square$

This implies, in particular, that the solutions of (4.2) with initial conditions in $\varvec{\ell }^p(w)$ are smooth if V is smooth.

4.7 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(\tau )$ are unbounded (see also (A.5)). Moreover, the following shows that a construction based on the standard Hille–Yosida theory of $C_0$-semigroups cannot be used to resolve this:

Proposition 4.9

If the weights grow faster than geometrically, the resolvent set of $S(\tau )$ is empty.

Proof

For any $\lambda , \xi \in {{\mathbb {C}}}$,

$$\begin{aligned} \mathcal {N}_n = \frac{1}{4^n} \bigl (4V(\tau )+\lambda ^2\bigr )^n \begin{pmatrix} 4\xi \\ 2\lambda \xi \\ \lambda ^2\xi \end{pmatrix} \end{aligned}$$

is in $\varvec{\ell }^p(w)$ because it grows at most geometrically, and a short calculation shows

$$\begin{aligned} (S(\tau ) - \lambda ) \mathcal {N} = \bigl ((A(\tau ) - \lambda ) \otimes \mathbb {1}+ B \otimes L \bigr )\mathcal {N} = 0. \end{aligned}$$

$\square$

However, even then it is sometimes possible to understand $S(\tau )$ and, somewhat surprisingly, also the evolution operator (4.15) as a bounded operator between two different weighted sequence spaces $\varvec{\ell }^p(v)$ and $\varvec{\ell }^p(w)$. In other words, both $S(\tau )$ and the evolution operator $U(\tau ,\tau _0)$ can be considered as unbounded operators on $\varvec{\ell }^p(v)$ with domain $\varvec{\ell }^p(w)$.

Let us consider the case (W2), i.e., the weights $w_n = (2n)! \omega ^{2n}$ and $v_n = (2n)! \upsilon ^{2n}$ with $\upsilon > \omega \ge 1$. As explained earlier, this can essentially be understood as the case where $\mathcal {G}-\mathcal {H}$ is analytic (with radius of convergence related to $\omega$ and $\upsilon$) in the radial variable. We apply Proposition 4.3, (4.12) and (4.14) to obtain for $n > 0$ the bound

$$\begin{aligned} {\Vert U_n(\tau ,\tau _0)\mathcal {M}\Vert }_{p,v}&\le \frac{1}{n!} \sum _{m=0}^{\lceil \frac{n+1}{2}\rceil } \left( {\begin{array}{c}n\\ m\end{array}}\right) |*|{\int _{\tau _0}^\tau \Vert A(\eta )\Vert \mathop {}\!\mathrm {d}\eta }^{n-m} \Vert B\Vert ^m {\Vert (\mathbb {1}\otimes L^m) \mathcal {M}\Vert }_{p,v} \\&\le \frac{2^n}{n!} \sum _{m=0}^{\lceil \frac{n+1}{2}\rceil } \left( {\begin{array}{c}n\\ m\end{array}}\right) r_m |*|{\int _{\tau _0}^\tau \sqrt{1+V(\eta )^2} \mathop {}\!\mathrm {d}\eta }^{n-m} {\Vert \mathcal M\Vert }_{p,w}, \end{aligned}$$

where we set

$$\begin{aligned} r_m = \sup _n \frac{w_{n+m}}{v_n} = \sup _n \frac{(2m+2n)!}{(2n)!} \frac{\omega ^{2m+2n}}{\upsilon ^{2n}}. \end{aligned}$$

We compute $r_0 = 1$ and, for $m > 0$,

(4.18)

where we applied Proposition A.2, the duplication formula for the Gamma function, and Gautschi’s inequality (A.7).

In Proposition A.3, we show that

$$\begin{aligned} \sum _{m=1}^{\lceil \frac{n+1}{2}\rceil } \frac{(m-1)!}{(n-m)!} \le {\left\{ \begin{array}{ll} \frac{3}{2} &{} \hbox {odd } n, \\ \frac{n}{2} + 1 &{} \hbox {even } n. \end{array}\right. } \end{aligned}$$

Thus, estimating the remaining factors in the sum by their supremum and recalling the definition (4.17) of $C_\omega (\tau ,\tau _0)$, we obtain for $n > 0$

$$\begin{aligned} {\Vert U_n(\tau ,\tau _0)\mathcal {M}\Vert }_{p,v} < C_0(\tau ,\tau _0)^n \biggl (\frac{1}{n!} + \frac{1}{\uppi } \left( \frac{n}{2} + 1\right) \left( \frac{\upsilon }{\omega }\right) ^\frac{1}{2}\! \left( \frac{2\upsilon \omega }{\upsilon -\omega }\right) ^{n+2}\biggr ) {\Vert \mathcal M\Vert }_{p,w}. \end{aligned}$$

(4.19)

Consequently, the time-ordered exponential (4.15) exists and can be bounded for sufficiently small time intervals by

$$\begin{aligned} \begin{aligned} {\Vert U(\tau ,\tau _0) \mathcal {M}\Vert }_{p,v}&\le \mathrm {e}^{C_0(\tau ,\tau _0)} {\Vert \mathcal M\Vert }_{p,w} \\ {}&\quad + \frac{C_0(\tau ,\tau _0) K(\upsilon ,\omega )^3}{2\uppi } \left( \frac{\upsilon }{\omega }\right) ^\frac{1}{2} \frac{3 - 2C_0(\tau ,\tau _0) K(\upsilon ,\omega )}{\bigl (1 - C_0(\tau ,\tau _0) K(\upsilon ,\omega )\bigr )^2} {\Vert \mathcal M\Vert }_{p,w}, \end{aligned}\end{aligned}$$

(4.20)

where we defined

$$\begin{aligned} K(\upsilon ,\omega ) :=\frac{2\upsilon \omega }{\upsilon -\omega }. \end{aligned}$$

Note that these inequalities only guarantee the existence of the exponential for a bounded interval of time, even if $\omega =1$ and $\upsilon$ 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:

Theorem 4.10

Let $\upsilon > \omega \ge 1$ and $p \ge 1$. Suppose that $I = [\tau _0,\tau _1]$ and $V \in C(I)$ such that

$$\begin{aligned} C_0(\tau _0,\tau _1) K(\upsilon ,\omega ) < 1. \end{aligned}$$

The evolution operator $U(\tau ,\tau _0)$, defined by (4.15), is the unique bounded operator from $\varvec{\ell }^p(w)$ to $\varvec{\ell }^p(v)$ with weights $w_n = (2n)! \omega ^{2n}$ and $v_n = (2n)! \upsilon ^{2n}$, such that the properties (E1)–(E3) hold in the interval I (in the strong sense from $\varvec{\ell }^p(w)$ to $\varvec{\ell }^p(v)$). Moreover, for $\mathcal {M} \in \varvec{\ell }^p(w)$, we have the bound (4.20) and the strong continuity property

$$\begin{aligned} \lim _{\tau \rightarrow \tau _0} {\Vert U(\tau ,\tau _0) \mathcal {M} - \mathcal {M}\Vert }_{p,v} = 0. \end{aligned}$$

Proof

We have already shown that the evolution operator $U(\tau ,\tau _0)$ is well defined as an operator from $\varvec{\ell }^p(w)$ to $\varvec{\ell }^p(v)$. In fact, the same estimates show that, for sufficiently small $\varepsilon > 0$, $U(\tau ,\tau _0)$ is bounded from $\varvec{\ell }^p(w)$ to $\varvec{\ell }^p(v_\varepsilon )$, where $v_{\varepsilon ,n} = (2n)! (\upsilon -\varepsilon )^{2n}$. Moreover, applying these estimates again, we find

$$\begin{aligned}\begin{aligned}&{\Vert U(\tau ,\tau _0) \mathcal {M} - \mathcal {M}\Vert }_{p,v} \\&\quad \le \left( \mathrm {e}^{C_0(\tau ,\tau _0)}- 1 + \frac{C_0(\tau ,\tau _0) K(\upsilon ,\omega )^3}{2\uppi } \left( \frac{\upsilon }{\omega }\right) ^\frac{1}{2} \frac{3 - 2 C_0(\tau ,\tau _0) K(\upsilon ,\omega )}{\bigl (1 - C_0(\tau ,\tau _0) K(\upsilon ,\omega )\bigr )^2} \right) {\Vert \mathcal M\Vert }_{p,w} \rightarrow 0 \end{aligned}\end{aligned}$$

as $\tau \rightarrow \tau _0$ because $C_0(\tau ,\tau _0) \rightarrow 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(\tau ,r)$ and $U(r,\tau _0)$. The final property (E3) is proved by differentiating (4.15) term by term using the formal relation $\partial _\tau U_n(\tau ,\tau _0) = S(\tau ) U_{n-1}(\tau ,\tau _0)$. The resulting expression is well defined because $U(\tau ,\tau _0)$ is bounded from $\varvec{\ell }^p(w)$ to $\varvec{\ell }^p(v_\varepsilon )$ and $S(\tau )$ is bounded from $\varvec{\ell }^p(v_\varepsilon )$ to $\varvec{\ell }^p(v)$. $\square$

Finally, we prove a perturbation theorem:

Theorem 4.11

In addition to the assumptions of Theorem 4.10, suppose that $\tilde{V} \in C(I)$ such that⁵

$$\begin{aligned} C_0(\tau ,\tau _0) \ge 2 |*|{\int _{\tau _0}^\tau \sqrt{1+\tilde{V}(\eta )^2} \mathop {}\!\mathrm {d}\eta } \quad \text {and}\quad 2C_0(\tau ,\tau _0) K(\upsilon ,\omega ) < 1. \end{aligned}$$

Denote by $U(\tau ,\tau _0)$ and $\tilde{U}(\tau ,\tau _0)$ the evolution operators for V and $\tilde{V}$. If $\mathcal {M}, \tilde{\mathcal M} \in \varvec{\ell }^p(w)$, we have

$$\begin{aligned}&{\Vert U(\tau ,\tau _0) \mathcal {M} - \tilde{U}(\tau ,\tau _0) \tilde{\mathcal M}\Vert }_{p,v} \\&\quad \le c(\tau ) {\Vert \mathcal {M} - \tilde{\mathcal M}\Vert }_{p,w} + 2 c(\tau )^2 {\Vert \tilde{\mathcal M}\Vert }_{p,w} \int _{\tau _0}^\tau |V(\eta )-\tilde{V}(\eta )| \mathop {}\!\mathrm {d}\eta , \end{aligned}$$

where

Proof

Set

$$\begin{aligned} \tilde{w}_n = (2n)! \tilde{\omega }^{2n}, \quad \tilde{\omega }= \frac{2\upsilon \omega }{\upsilon +\omega }. \end{aligned}$$

Observe that

$$\begin{aligned} \frac{\tilde{\omega }\omega }{\tilde{\omega }-\omega } = \frac{2\upsilon \omega }{\upsilon -\omega } = \frac{\upsilon \tilde{\omega }}{\upsilon -\tilde{\omega }} \end{aligned}$$

and

$$\begin{aligned} \frac{\tilde{\omega }}{\omega } = \frac{2\upsilon }{\upsilon +\omega }< \frac{\upsilon }{\omega }, \quad \frac{\upsilon }{\tilde{\omega }} = \frac{\upsilon +\omega }{2\omega } < \frac{\upsilon }{\omega }. \end{aligned}$$

Further, note that $2C_0(\tau ,\tau _0)K(\upsilon ,\omega ) < 1$ by assumption. Therefore,

$$\begin{aligned}&{\Vert U(\tau ,\tau _0) \mathcal {M}\Vert }_{p,v} \le c(\tau ) {\Vert \mathcal M\Vert }_{p,\tilde{w}}, \quad {\Vert \tilde{U}(\tau ,\tau _0) \mathcal {M}\Vert }_{p,v} \le c(\tau ) {\Vert \mathcal M\Vert }_{p,\tilde{w}}, \nonumber \\&{\Vert U(\tau ,\tau _0) \mathcal {M}\Vert }_{p,\tilde{w}} \le c(\tau ) {\Vert \mathcal M\Vert }_{p,w}, \quad {\Vert \tilde{U}(\tau ,\tau _0) \mathcal {M}\Vert }_{p,\tilde{w}} \le c(\tau ) {\Vert \mathcal M\Vert }_{p,w}.\nonumber \\ \end{aligned}$$

(4.21)

Now, note the elementary identity

$$\begin{aligned} U(\tau ,\tau _0) \mathcal {M} - \tilde{U}(\tau ,\tau _0) \tilde{\mathcal M} = U(\tau ,\tau _0) (\mathcal {M} - \tilde{\mathcal M}) + \bigl (U(\tau ,\tau _0) - \tilde{U}(\tau ,\tau _0)\bigr ) \tilde{\mathcal M}. \end{aligned}$$

For the first summand, we find

$$\begin{aligned} {\Vert U(\tau ,\tau _0) (\mathcal {M} - \tilde{\mathcal M})\Vert }_{p,v} \le {\Vert U(\tau ,\tau _0) (\mathcal {M} - \tilde{\mathcal M})\Vert }_{p,\tilde{w}} \le c(\tau ) {\Vert \mathcal {M} - \tilde{\mathcal M}\Vert }_{p,w}. \end{aligned}$$

For the second summand, we find, using the fundamental theorem of calculus (for Banach space-valued integrals),

$$\begin{aligned} \bigl (U(\tau ,\tau _0) - \tilde{U}(\tau ,\tau _0)\bigr ) \tilde{\mathcal M} = \int _{\tau _0}^\tau U(\tau ,\eta ) \bigl ( \tilde{S}(\eta ) - S(\eta ) \bigr ) \tilde{U}(\eta ,\tau _0) \tilde{\mathcal M} \mathop {}\!\mathrm {d}\eta , \end{aligned}$$

and thus,

$$\begin{aligned} {\Vert U(\tau ,\tau _0) \tilde{\mathcal M} - \tilde{U}(\tau ,\tau _0) \tilde{\mathcal M}\Vert }_{p,v}&\le \int _{\tau _0}^\tau \Vert [\Vert \big ]{U(\tau ,\eta ) \bigl ( \tilde{S}(\eta ) - S(\eta ) \bigr ) \tilde{U}(\eta ,\tau _0) \tilde{\mathcal M}}_{p,v} \mathop {}\!\mathrm {d}\eta \\&\le 2 c(\tau ) \int _{\tau _0}^\tau |V(\eta ) - \tilde{V}(\eta )| {\Vert \tilde{U}(\eta ,\tau _0) \tilde{\mathcal M}\Vert }_{p,\tilde{w}} \mathop {}\!\mathrm {d}\eta \\&\le 2 c(\tau )^2 {\Vert \tilde{\mathcal M}\Vert }_{p,w} \int _{\tau _0}^\tau |V(\eta )-\tilde{V}(\eta )| \mathop {}\!\mathrm {d}\eta . \end{aligned}$$

$\square$

As for the case of geometrically growing weights, under more stringent assumptions on the potential V, one obtains the following:

Proposition 4.12

If, in addition to the assumptions of Theorem 4.10, $V \in C^k(I)$, then $\partial _\tau ^{k+1} U(\tau ,\tau _0)$ is bounded from $\varvec{\ell }^p(w)$ to $\varvec{\ell }^p(v)$.

Proof

As in the proof of Theorem 4.10, $U(\tau ,\tau _0)$ is bounded from $\varvec{\ell }^p(w)$ to $\varvec{\ell }^p(v_\varepsilon )$. Moreover, it is not difficult to see that products of $S(\tau )$ and its derivatives are bounded from $\varvec{\ell }^p(v_\epsilon )$ to $\varvec{\ell }^p(v)$. (Concrete estimates similar to those before Theorem 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). $\square$

As before, this implies that solutions of (4.2) with initial conditions in $\varvec{\ell }^p(w)$ are smooth if V is smooth.

5 Abstract Semiclassical Einstein Equation

The SCE on FLRW spacetimes contains only one geometric degree of freedom—the scale factor $a = a(\tau )$. 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 \in {{\mathbb {N}}}$ and $\tau \in I \subset {{\mathbb {R}}}$, consider the initial value problem for the quasilinear system

$$\begin{aligned} \partial _\tau ^{k+1} a(\tau )&= f(\tau , \varvec{a}, \mathcal {M}), \end{aligned}$$

(5.1a)

$$\begin{aligned} \partial _\tau \mathcal {M}(\tau )&= S(\tau , \varvec{a}) \mathcal {M}(\tau ), \end{aligned}$$

(5.1b)

where

• $a=a(\tau ) \in {{\mathbb {R}}}$ is the scale factor and $\varvec{a} = (a, a', \dotsc , a^{(k)})$ its k-jet,

• $\mathcal {M}=\mathcal {M}(\tau ) \in \varvec{\ell }^p(v)$ is the sequence of coincidence limits, as defined in Sect. 4, for some $p \ge 1$ and sequence of weights v,

• $S(\tau , \varvec{a}) = S\bigl (\tau , a(\tau ), a'(\tau ), \dotsc , a^{(k)}(\tau ) \bigr )$ is the generator of the dynamics of $\mathcal {M}$, as defined in (4.13), for a potential $V(\tau , \varvec{a}) = V\bigl (\tau , a(\tau ), a'(\tau ), \dotsc , a^{(k)}(\tau ) \bigr )$ which depends on the k-jet of the scale factor and may also have an explicit time-dependence, and

• $f(\tau , \varvec{a}, \mathcal {M}) = f\bigl (\tau , a(\tau ), a'(\tau ), \dotsc , a^{(k)}(\tau ), \mathcal {M}(\tau )\bigr )$ specifies the dynamics of the scale factor including a possible explicit time dependence and a back reaction by the quantum field via $\mathcal {M}$.

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(\tau , \varvec{a},\mathcal {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_\mathrm {tr}$ coming from the trace of the quantized stress-energy tensor.

5.1 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.

Theorem 5.1

Let $\varvec{b}_c \in {{\mathbb {R}}}^{k+1}$ and set

$$\begin{aligned} \mathcal {B}_r :=\bigl \{ \varvec{b} \in {{\mathbb {R}}}^{k+1} \;\big |\; \Vert \varvec{b}-\varvec{b}_c\Vert \le r \bigr \}, \quad r > 0, \end{aligned}$$

where we use the Euclidean norm on ${{\mathbb {R}}}^{k+1}$. For fixed $R > 0$ and weight sequence w, consider $\varvec{a}_\mathrm {init}\in \mathcal {B}_{R/2}$ and $\mathcal {M}_\mathrm {init}\in \varvec{\ell }^p(w)$. Suppose that there is $I :=[\tau _\mathrm {init},\tau _\mathrm {fin}] \subset {{\mathbb {R}}}$ and a weight sequence v such that the following holds:

• (i) For each $\varvec{a} \in C(I; \mathcal {B}_R)$, there exists a unique solution $\mathcal {M}[\varvec{a}] \in C^1(I; \varvec{\ell }^p(v))$ to (5.1b) with initial value $\mathcal {M}[\varvec{a}](\tau _\mathrm {init}) = \mathcal {M}_\mathrm {init}$. Set

  $$\begin{aligned} \mu _\mathcal {M} = \sup _{\tau \in I, \varvec{a} \in C(I;\mathcal {B}_R)} {\Vert \mathcal {M}[\varvec{a}](\tau )\Vert }_{p,v}. \end{aligned}$$

• (ii) There is $L_{\mathcal M}>0$ such that, for all $\varvec{a},\tilde{\varvec{a}} \in C(I; \mathcal {B}_R)$,

  $$\begin{aligned} {\Vert \mathcal {M}[\varvec{a}](\tau )-\mathcal {M}[\tilde{\varvec{a}}](\tau )\Vert }_{p,v} \le L_{\mathcal M} \int _{\tau _\mathrm {init}}^\tau \Vert \varvec{a}(\eta ) - \tilde{\varvec{a}}(\eta )\Vert \mathop {}\!\mathrm {d}\eta . \end{aligned}$$

  (5.2)

• (iii) For each $\varvec{b},\tilde{\varvec{b}} \in \mathcal {B}_R$ and $\mathcal {M},\tilde{\mathcal M} \in \varvec{\ell }^p(v)$ bounded by $\mu _\mathcal {M}$, the map $\tau \mapsto f(\tau ,\varvec{b},\mathcal {M})$ is continuous for $\tau \in I$, and there are $\mu _f > 0$, $L_f > 0$ such that

  $$\begin{aligned} |f(\tau ,\varvec{b},\mathcal {M})|&\le \mu _f, \end{aligned}$$

  (5.3)
  $$\begin{aligned} |f(\tau ,\varvec{b},\mathcal {M})-f(\tau ,\tilde{\varvec{b}},\tilde{\mathcal M})|&\le L_f \bigl (\Vert \varvec{b}-\tilde{\varvec{b}}\Vert + {\Vert \mathcal {M}-\tilde{\mathcal M}\Vert }_{p,v}\bigr ) \end{aligned}$$

  (5.4)

  uniformly for $\tau \in I$.

Then, there exists $\tau _\mathrm {stop} \in (\tau _\mathrm {init},\tau _\mathrm {fin}]$ such that (5.1) has a unique (local) solution

$$\begin{aligned} (a, \mathcal {M}) \in C^{k+1}(J) \times C^1(J; \varvec{\ell }^p(v)) \end{aligned}$$

in $J = [\tau _\mathrm {init},\tau _\mathrm {stop}]$, satisfying the initial conditions $\varvec{a}(\tau _\mathrm {init}) = \varvec{a}_\mathrm {init}, \mathcal {M}(\tau _\mathrm {init}) = \mathcal {M}_\mathrm {init}$ and the bounds $\varvec{a}(\tau ) \in \mathcal {B}_R$, $\Vert \mathcal {M}(\tau )\Vert _{p,v} \le \mu _\mathcal {M}$ for $\tau \in J$.

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.

Proposition 5.2

Consider Theorem 5.1 with the weights

• (a) $v_n = w_n = \omega ^n$ with $\omega > 0$ (geometrically growing weights), or

• (b) $v_n = (2n!) \upsilon ^{2n}$ and $w_n = (2n!) \omega ^{2n}$ with $\upsilon ,\omega > 0$ (factorially growing weights).

Further, suppose that $\tau \mapsto V(\tau ,\varvec{b})$ is continuous for all $\varvec{b} \in \mathcal {B}_R$, and there is $L_V > 0$ such that

$$\begin{aligned} |V(\tau ,\varvec{b})-V(\tau ,\tilde{\varvec{b}})| \le L_V \Vert \varvec{b}-\tilde{\varvec{b}}\Vert \end{aligned}$$

for each $\varvec{b}, \tilde{\varvec{b}} \in \mathcal {B}_R$, uniformly for $\tau$ in some interval $I'$. Then, there exists $I = [\tau _\mathrm {init},\tau _\mathrm {fin}] \subset I'$ such that the assumptions (i)–(iii) of Theorem 5.1 are satisfied.

Proof

The assertion follows from Theorems 4.4 and 4.7 for (a), and from Theorems 4.10 and 4.11 for (b). $\square$

Proposition 5.3

If the assumptions of Proposition 5.2 hold with geometrically growing weights (a), then the unique local solution $(a,\mathcal {M})$ of Theorem 5.1 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 Theorem 5.1 and Proposition 5.2). $\square$

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 $\mathcal {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 Theorem 5.1 imply that we can solve (5.1) with the k-jet $\varvec{a}$ of a inside a ball $\mathcal {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,\mathcal {M}) \in C^{k+1}(I) \times C^1(I; \varvec{\ell }^p(w))$ is a solution of (5.1) in the interval $I \subset {{\mathbb {R}}}$ with $\varvec{a} \in C^1(I; \mathcal {B}_R)$ (for $\mathcal {B}_R$ as in Theorem 5.1). If $V \in C^l(I \times \mathcal {B}_R)$ and $f \in C^l(I \times \mathcal {B}_R \times \varvec{\ell }^p(w))$ for some $l \in {{\mathbb {N}}}$, then $a \in C^{k+l}(I)$ and $\mathcal {M} \in C^l(I; \ell ^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 Proposition 4.12. $\square$

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].

5.2 Continuous Dependence on Initial Data and Parameters

To study the continuous dependence of the abstract SCE (5.1) on its initial data, the potential V and the ‘back-reaction’ function f, we consider another quasilinear equation of the form (5.1):

$$\begin{aligned} \partial _\tau ^{k+1} \tilde{a}(\tau )&= \tilde{f}(\tau , \tilde{\varvec{a}}, \tilde{\mathcal M}), \end{aligned}$$

(5.5a)

$$\begin{aligned} \partial _\tau \tilde{\mathcal M}(\tau )&= \tilde{S}(\tau , \tilde{\varvec{a}}) \tilde{\mathcal M}(\tau ), \end{aligned}$$

(5.5b)

with initial conditions $\tilde{\varvec{a}}_\mathrm {init}\in {{\mathbb {R}}}^{k+1}$ and $\tilde{\mathcal M}_\mathrm {init}\in \varvec{\ell }^p(w)$. Above, $\tilde{S}(\tau , \tilde{\varvec{a}})$ is given by (4.13) for a potential $\tilde{V}(\tau , \tilde{\varvec{a}})$.

Theorem 5.6

Suppose that the coefficients of (5.1) and (5.5) satisfy the assumptions (i), (iii) of Theorem 5.1 for the same constants. Furthermore, (with the obvious notation) assume that, instead of assumption (ii) of Theorem 5.1, we have:

• (ii’) There is $L_{\mathcal M}' > 0$ such that, for all $\varvec{a}, \tilde{\varvec{a}} \in C(I; \mathcal {B}_R)$,

  $$\begin{aligned} {\Vert \mathcal {M}[\varvec{a}](\tau ) {-} \tilde{\mathcal M}[\tilde{\varvec{a}}](\tau )\Vert }_{p,v} \le L_{\mathcal M}' \left( {\Vert \mathcal {M}_\mathrm {init}{-} \tilde{\mathcal M}_\mathrm {init}\Vert }_{p,w} + \int _{\tau _\mathrm {init}}^\tau \Vert \varvec{a}(\eta ) - \tilde{\varvec{a}}(\eta )\Vert \mathop {}\!\mathrm {d}\eta \right) . \end{aligned}$$

Then, the two solutions $(a, \mathcal {M})$ and $(\tilde{a}, \tilde{\mathcal M})$ on $J = [\tau _\mathrm {init},\tau _\mathrm {stop}] \subset I$ satisfy

$$\begin{aligned} \Vert \varvec{a}(\tau ) - \tilde{\varvec{a}}(\tau )\Vert\le & {} \Vert \varvec{a}_\mathrm {init}- \tilde{\varvec{a}}_\mathrm {init}\Vert \mathrm {e}^{K(\tau ,\tau _\mathrm {init})} \nonumber \\&+ (\mu _f' + L_f L_\mathcal {M}' \Vert \mathcal {M}_\mathrm {init}- \tilde{\mathcal M}_\mathrm {init}\Vert _{p,w}) \int _{\tau _\mathrm {init}}^\tau \mathrm {e}^{K(\tau ,\eta )} \mathop {}\!\mathrm {d}\eta \end{aligned}$$

(5.6)

for $\tau \in J$, where we set $K(\tau ,\eta ) = (1+L_f) (\tau - \eta ) + \frac{1}{2} L_f L_\mathcal {M}' (\tau - \eta )^2$ and

$$\begin{aligned} \mu _f' = \sup {} |f(\tau , \varvec{b}, \mathcal {M}) - \tilde{f}(\tau , \varvec{b}, \mathcal {M})| \end{aligned}$$

with the supremum being taken over over $\tau \in J$, $\varvec{b} \in \mathcal {B}_R$ and $\Vert \mathcal M\Vert _{p,v} \le \mu _\mathcal {M}$.

The proof is 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.

5.3 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 $\omega _2^\mathrm {reg} = \omega _2 - H_n$, where n is sufficiently large. It follows by a straightforward calculation from the definitions of Sect. 2 that

$$\begin{aligned}{}[\omega _2^\mathrm {reg}]&= \frac{1}{a^2} \bigl ( \mathcal {M}_{\varphi \varphi ,0} + \mathcal {H}_{\varphi \varphi ,1} - \tilde{\mathcal {H}}_{\varphi \varphi ,1} \bigr ) \Bigl |_{r=0}, \\ [(\mathbb {1}\otimes \Delta ) \omega _2^\mathrm {reg}]&= \frac{1}{a^2} \bigl ( \mathcal {M}_{\varphi \varphi ,1} + \Delta _r (\mathcal {H}_{\varphi \varphi ,2} - \tilde{\mathcal {H}}_{\varphi \varphi ,2}) \bigr ) \Bigl |_{r=0}, \\ [(\partial _\tau \otimes \partial _\tau ) \omega _2^\mathrm {reg}]&= \frac{1}{a^2} \bigl ( \mathcal {M}_{\pi \pi ,0} {+} \mathcal {H}_{\pi \pi ,2} - \tilde{\mathcal {H}}_{\pi \pi ,2} \bigr ) {+} \frac{a^{\prime \,2}}{a^4} \bigl ( \mathcal {M}_{\varphi \varphi ,0} + \mathcal {H}_{\varphi \varphi ,1} - \tilde{\mathcal {H}}_{\varphi \varphi ,1} \bigr ) \\ {}&\quad - 2\frac{a'}{a^3} \bigl ( \mathcal {M}_{(\varphi \pi ),0} + \mathcal {H}_{(\varphi \pi ),2} - \tilde{\mathcal {H}}_{(\varphi \pi ),2} \bigr ) \Bigl |_{r=0}. \end{aligned}$$

Thus, combining the results from Sects. 2 and 3, the traced SCE (3.1) can be expanded to the rather long equation

$$\begin{aligned}\begin{aligned} 0&= \left( -12 (3 c_3 + c_4) -\frac{1}{480\uppi ^2} + \frac{6\xi -1}{48\uppi ^2} + \frac{(6\xi -1)^2}{16\uppi ^2} \log (a \lambda _0) \right) \\ {}&\quad \times \left( \frac{a^{(4)}}{a^5} - 4 \frac{a^{(3)} a'}{a^6} - 3 \frac{a^{\prime \prime \,2}}{a^6} + 6 \frac{a'' a^{\prime \,2}}{a^7} \right) \\ {}&\quad + \frac{(6\xi -1)^2}{32\uppi ^2} \left( 4 \frac{a^{(3)} a'}{a^6} + 3 \frac{a^{\prime \prime \,2}}{a^6} - 10 \frac{a'' a^{\prime \,2}}{a^7} \right) + \frac{1}{240\uppi ^2} \left( -\frac{a'' a^{\prime \,2}}{a^7} + \frac{a^{\prime \,4}}{a^8} \right) \\ {}&\quad + \left( \frac{6}{\kappa } + m^2 \Bigl ( -6 c_2 + \frac{1}{48\uppi ^2} + \frac{6\xi -1}{8\uppi ^2} \bigl (1 + \log (a \lambda _0) \bigr ) \Bigr ) \right) \frac{a''}{a^3} \\ {}&\quad + \frac{(6\xi -1) m^2}{16\uppi ^2} \frac{a^{\prime \,2}}{a^4} + m^4 \left( 4 c_1 + \frac{1}{32\uppi ^2} + \frac{1}{8\uppi ^2} \log (a \lambda _0) \right) - \frac{m^2}{a^2} \mathcal {M}_{\varphi \varphi ,0} \\ {}&\quad + (6\xi -1) \left( \Bigl ( 6\xi \frac{a''}{a^5} - \frac{a^{\prime \,2}}{a^6} + \frac{m^2}{a^2} \Bigr ) \mathcal {M}_{\varphi \varphi ,0} + 2 \frac{a'}{a^5} \mathcal {M}_{(\varphi \pi ),0} - \frac{1}{a^4} \bigl ( \mathcal {M}_{\pi \pi ,0} + \mathcal {M}_{\varphi \varphi ,1} \bigr ) \right) . \end{aligned}\end{aligned}$$

We remark that the first line of this equation is due to terms proportional to $\Box R$.

In the general case, this equation can be rewritten as quasilinear fourth order equation of the form

$$\begin{aligned} \partial _\tau a^{(3)} = f_\mathrm {tr}(a, a', a'', a^{(3)}, \mathcal {M}_0, \mathcal {M}_1), \end{aligned}$$

(5.7)

which can be solved using Theorem 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

$$\begin{aligned} 30 (6\xi -1)^2 \log (a \lambda _0) = 11 + 5760\uppi ^2 (3c_3+c_4) - 60\xi , \end{aligned}$$

(5.8)

which must be taken into account for the choice of the ball $\mathcal {B}_R$ in Theorem 5.1.

We do not see the instability near the Minkowski solution described in [57, 58], viz. the equation is continuous in a neighborhood 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 $\lambda _0$, relating the length scales $\mu$ and $\lambda$ in (2.19). While the choice of $\lambda$ is related to the value of the renormalization parameters $c_1, c_2, c_3$ and $c_4$, the value of $\mu$ (and thus $\lambda _0$) is related to our construction of $\mathcal {M}$. Of course, changing $\lambda _0$ is not without effect and requires a corresponding change of $\mathcal {M}$. It is important to note that such a modification of $\mathcal {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

$$\begin{aligned} \xi = \frac{1}{6} \quad \text {and}\quad 3c_3+c_4 = -\frac{1}{5760\uppi ^2}, \end{aligned}$$

(5.9)

the fourth- and third-order terms drop out and the equation can be rewritten as

$$\begin{aligned}\begin{aligned} a''&= \left( \frac{a^{\prime \,2}}{a^4} - 1440\uppi ^2 \kappa ^{-1} + (1440\uppi ^2 c_2 - 5) m^2 \right) ^{-1} \\ {}&\quad \times \left( \frac{a^{\prime \,4}}{a^5} + \frac{1}{2} m^4 a^3 \bigl ( 1920\uppi ^2 c_1 + 15 + 60 \log (a \lambda _0) \bigr ) - 240\uppi ^2 m^2 a \mathcal {M}_{\varphi \varphi ,0} \right) \end{aligned}\end{aligned}$$

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

$$\begin{aligned} \frac{a^{\prime \,2}}{a^4} = 1440\uppi ^2 \kappa ^{-1} - (1440\uppi ^2 c_2 - 5) m^2, \end{aligned}$$

(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 $\mathcal {B}_R$ in Theorem 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:

Theorem 5.7

The traced SCE is of the form (5.1) and can be locally solved (using Theorem 5.1 and Proposition 5.2, 5.3, 5.5), yielding a unique smooth solution which depends continuously on the initial data and parameters. In the case of geometrically growing weights for the moments $\mathcal {M}$ (i.e., vacuum-like states), a maximal or global solution exists. In the generic fourth-order case, the maximal solution exists up to a big bang/crunch $a = 0$, the logarithmic singularity (5.8), or a blow-up $a \rightarrow \infty$ (big rip). In the second-order case (i.e., (5.9) are satisfied), the maximal solutions exist up to a big bang/crunch, the singularity (5.10), or a big rip.

Let us also briefly explain the philosophy behind the choice of the renormalization parameters $c_\bullet$. 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].

5.4 Application to the Energy Component of the SCE

Analogous to the traced SCE, the energy component of the SCE can be expanded as

$$\begin{aligned} \begin{aligned} 0&= \left( 6 (3 c_3 + c_4) + \frac{1}{960\uppi ^2} - \frac{6\xi -1}{96\uppi ^2} - \frac{(6\xi -1)^2}{32\uppi ^2} \log (a \lambda _0) \right) \\ {}&\quad \left( 2 \frac{a^{(3)} a'}{a^4} - \frac{a^{\prime \prime \,2}}{a^4} - 4 \frac{a'' a^{\prime \,2}}{a^5} \right) - \frac{(6\xi -1)^2}{16\uppi ^2} \frac{a'' a^{\prime \,2}}{a^5} + \frac{1}{960\uppi ^2} \frac{a^{\prime \,4}}{a^6} \\&\quad + \left( -\frac{3}{\kappa } + m^2 \Bigl ( 3 c_2 - \frac{1}{96\uppi ^2} - \frac{6\xi -1}{16\uppi ^2} \bigl (1 + \log (a \lambda _0) \bigr ) \Bigr ) \right) \frac{a^{\prime \,2}}{a^2} \\ {}&\quad - m^4 \left( c_1 + \frac{1}{32\uppi ^2} \log (a \lambda _0) \right) a^2 + \frac{m^2}{2} \mathcal {M}_{\varphi \varphi ,0} \\&\quad + (6\xi -1) \left( -\frac{a^{\prime \,2}}{2a^4} \mathcal {M}_{\varphi \varphi ,0} + \frac{a'}{a^3} \mathcal {M}_{(\varphi \pi ),0} \right) + \frac{1}{2a^2} \bigl ( \mathcal {M}_{\pi \pi ,0} - \mathcal {M}_{\varphi \varphi ,1} \bigr ), \end{aligned}\end{aligned}$$

(5.11)

In the special case (5.9), already considered above, this equation simplifies to

$$\begin{aligned}\begin{aligned} 0&= \frac{1}{960\uppi ^2} \frac{a^{\prime \,4}}{a^6} + \left( -\frac{3}{\kappa } + m^2 \Bigl ( 3 c_2 - \frac{1}{96\uppi ^2} \Bigr ) \right) \frac{a^{\prime \,2}}{a^2}\\&\quad - m^4 \left( c_1 + \frac{1}{32\uppi ^2} \log (a \lambda _0) \right) a^2 \\&\quad + \frac{m^2}{2} \mathcal {M}_{\varphi \varphi ,0} + \frac{1}{2a^2} \bigl ( \mathcal {M}_{\pi \pi ,0} - \mathcal {M}_{\varphi \varphi ,1} \bigr ). \end{aligned}\end{aligned}$$

The general equation can be rewritten as a quasilinear third-order equation of the form

$$\begin{aligned} \partial _\tau a'' = f_\mathrm {en}(a, a', a'', \mathcal {M}_0, \mathcal {M}_1), \end{aligned}$$

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_\mathrm {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.

5.5 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:

$$\begin{aligned} 0&= \frac{1}{2} \mathcal {M}_{\pi \pi ,0} - \frac{1}{2} \mathcal {M}_{\varphi \varphi ,1} + \frac{m^2}{2} \mathcal {M}_{\varphi \varphi ,0} - m^4 \left( c_1 + \frac{1}{32\uppi ^2} \log (\lambda _0) \right) \\&= \mathcal {M}_{\pi \pi ,0} - m^4 \left( c_1 + \frac{1}{32\uppi ^2} \log \lambda _0 \right) , \end{aligned}$$

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 example, using the results from Sects. 4.2 and 4.3, we find

$$\begin{aligned} \mathcal {M}_{\pi \pi ,0} = \frac{m^4}{32\uppi ^2} \left( \frac{1}{4} + \log \Bigl ( \frac{1}{2} m \mu \Bigr ) \right) = \frac{m^4}{32\uppi ^2} \left( -\frac{3}{4} + \gamma + \frac{1}{2} \log \Bigl ( \frac{1}{2} m^2 \lambda ^2 \Bigr ) + \log \lambda _0 \right) , \end{aligned}$$

for the Minkowski vacuum and

$$\begin{aligned} \mathcal {M}_{\pi \pi ,0} = \frac{\uppi ^4}{15 \beta ^4} \end{aligned}$$

for the thermal state. That is, in these examples, $c_1$ is fixed once the mass m and the regularization scale $\lambda$, resp. the inverse temperature $\beta$ are fixed. Note that the arbitrary scale $\lambda _0$ cancels out.

It follows then from (3.9) and (3.10), for any time translation-invariant state on Minkowski spacetime, that

$$\begin{aligned} G_{00} = \kappa {\langle T_{00}^\mathrm {ren} \rangle }_\omega \quad \Longleftrightarrow \quad -R = \kappa {\langle T^\mathrm {ren} \rangle }_\omega . \end{aligned}$$

Therefore, we have:

Theorem 5.8

Any time translation-invariant state on Minkowski spacetime satisfies the SCE on Minkowski spacetime, provided the renormalization parameter $c_1$ is chosen consistently with the energy constraint.

5.6 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 neighborhood 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 $\mathcal {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 Lemma 5.9 this would be a neighborhood of $\tau _\mathrm {tow}$), and then calculate the corresponding moments. Then, at some point in time (denoted $\tau _\mathrm {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 $\tau _\mathrm {free}$ below), we have a solution to the SCE with a Hadamard state.

Fig. 1

Schematic illustration of the tow-in procedure described in Lemma 5.9, including various relevant time intervals. In the past of $\tau _\mathrm {init}$ the SCE is ‘turned off’ and the solution is determined by a chosen scale factor $a_\mathrm {tow}$ which evolves to the desired initial value $a_\mathrm {init}$. Then, in the interval $(\tau _\mathrm {init}, \tau _\mathrm {free})$, the SCE is smoothly turned on using the function $\chi$. Shrinking this time interval, the solution at $\tau _\mathrm {free}$ can be brought arbitrarily close to its value at $\tau _\mathrm {init}$. Finally, after $\tau _\mathrm {free}$ the solution evolves freely via the SCE without the forcing caused by $a_\mathrm {tow}$

To put this idea on a solid footing, we prove the following technical lemma, some aspects of which are also illustrated in Fig. 1:

Lemma 5.9

Let $[\tau _\mathrm {tow}, \tau _\mathrm {fin}] \subset {{\mathbb {R}}}$, $\tau _\mathrm {init}\in (\tau _\mathrm {tow}, \tau _\mathrm {fin})$ and $a_\mathrm {tow}\in C^{k+1}[\tau _\mathrm {tow}, \tau _\mathrm {fin}]$. Suppose that $\mathcal {M}_\mathrm {tow}\in C^1([\tau _\mathrm {tow}, \tau _\mathrm {fin}]; \varvec{\ell }^p(w))$ is a solution of (5.1b) for the scale factor $a_\mathrm {tow}$ and the weight sequence w. Further, given appropriately chosen $\varvec{b}_c$, R, v, suppose that the assumptions of Theorem 5.1 hold for the initial data $\varvec{a}_\mathrm {init}= \varvec{a}_\mathrm {tow}(\tau _\mathrm {init})$, $\mathcal {M}_\mathrm {init}= \mathcal {M}_\mathrm {tow}(\tau _\mathrm {init})$. Then, there exists $\tau _\mathrm {stop} \in (\tau _\mathrm {init}, \tau _\mathrm {fin}]$ such that for any $\varepsilon > 0$ one can find $\tau _\mathrm {free}\in (\tau _\mathrm {init}, \tau _\mathrm {stop})$ and $(a, \mathcal {M})$ satisfying

• (i) $(a, \mathcal {M}) = (a_\mathrm {tow}, \mathcal {M}_\mathrm {tow})$ in $[\tau _\mathrm {tow},\tau _\mathrm {init}]$,

• (ii) $(a, \mathcal {M})$ solves (5.1) in $[\tau _\mathrm {free}, \tau _\mathrm {stop}]$,

• (iii) $(a, \mathcal {M}) \in C^{k+1}(J') \times C^1(J';\varvec{\ell }^p(v))$ with $J' = [\tau _\mathrm {tow}, \tau _\mathrm {stop}]$,

• (iv) $\Vert \varvec{a}(\tau _\mathrm {free}) - \varvec{a}(\tau _\mathrm {init})\Vert \le \varepsilon$ and $\Vert \mathcal {M}(\tau _\mathrm {free}) - \mathcal {M}(\tau _\mathrm {init})\Vert _{p,v} \le \varepsilon$,

• (v) $(a, \mathcal {M})$ is smooth if $a_\mathrm {tow}, f$ and V are smooth.

Proof

For any $\chi \in C^\infty ([\tau _\mathrm {tow}, \tau _\mathrm {fin}]; [0,1])$, define

$$\begin{aligned} f_\chi (\tau , \varvec{a}, \mathcal {M}) = \bigl ( 1-\chi (\tau ) \bigr ) a_\mathrm {tow}^{(k+1)}(\tau ) + \chi (\tau ) f(\tau , \varvec{a}, \mathcal {M}). \end{aligned}$$

We want to apply Theorem 5.1 with $f_\chi$ instead of f. Assumptions (i) and (ii) of Theorem 5.1 continue to hold with the same constants. Estimating $\chi$ by 1, (5.3) in assumption (iii) can be replaced by

$$\begin{aligned} |f_\chi (\tau , \varvec{a}, \mathcal {M})| \le \sup _{\eta \in I}{} |a_\mathrm {tow}^{(k+1)}(\eta )| + \mu _f, \end{aligned}$$

(5.12)

which holds for $\tau \in I$ and the right-hand side of (5.4) is unchanged. Hence, Theorem 5.1 can be applied to find $\tau _\mathrm {stop} \in (\tau _\mathrm {init}, \tau _\mathrm {fin}]$ and a solution

$$\begin{aligned} (a_\chi , \mathcal {M}_\chi ) \in C^{k+1}(J) \times C^1(J; \varvec{\ell }^p(v)) \end{aligned}$$

in $J = [\tau _\mathrm {init}, \tau _\mathrm {stop}]$ for (5.1) with f replaced by $f_\chi$.

Define $(a, \mathcal {M})$ by gluing $(a_\mathrm {tow}, \mathcal {M}_\mathrm {tow})$ and $(a_\chi , \mathcal {M}_\chi )$. Choosing $\chi$ such that $\chi (\tau ) = 0$ for $\tau \le \tau _\mathrm {init}$ and $\chi (\tau ) = 1$ for $\tau \ge \tau _\mathrm {free}$ with $\tau _\mathrm {free}\in (\tau _\mathrm {init}, \tau _\mathrm {stop})$ arbitrary, properties (i)–(iii) are obviously satisfied. Since estimates in the previous paragraph are uniform in $\chi$, (iv) follows by continuity of the solution, viz., if $\tau _\mathrm {free}$ is chosen sufficiently close to $\tau _\mathrm {init}$, the estimates hold. Finally, in the smooth case, apply Proposition 5.5 to see that (v) holds. $\square$

Applied to the traced SCE, where V and f are smooth for appropriately chosen $\mathcal {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:

Proposition 5.10

If $\mathcal {M}_\mathrm {tow}$ is given by a Hadamard state on the FLRW spacetime with smooth scale function $a_\mathrm {tow}$, the solution $(a, \mathcal {M})$ in Lemma 5.9 has moments $\mathcal {M}$ for a Hadamard state with respect to a.

While the ‘tow in’ construction described above can be used to prove the existence of Hadamard states fulfilling the trace equation of the SCE, the status of the energy constraint remains unsatisfactory: Although the SCE preserves the energy constraint, this is not case for the deformed version used in the proof of Lemma 5.9. In other words, although the solutions constructed by that lemma satisfy the traced SCE, they do not satisfy the energy equation and hence not the full SCE. To fix this issue, we need to amend this construction through an additional step. Namely, we give up control over the initial datum of the third derivative of the scale factor and change it in such a way that the energy constraint is satisfied again from $\tau _\mathrm {free}$ onwards. That this can indeed be done consistently is shown in the following theorem.

Theorem 5.11

Given $b_0, b_1, b_2 \in {{\mathbb {R}}}$ such that

$$\begin{aligned} 6 (3 c_3 + c_4) + \frac{1}{960\uppi ^2} - \frac{6\xi -1}{96\uppi ^2} - \frac{(6\xi -1)^2}{32\uppi ^2} \log (b_0 \lambda _0) \ne 0, \quad b_0 > 0, \quad b_1 \ne 0, \end{aligned}$$

(5.13)

and $\varepsilon > 0$, there exists a time interval $J = [\tau _\mathrm {free}, \tau _\mathrm {stop}] \subset {{\mathbb {R}}}$, a weight sequence w and $(a, \mathcal {M}) \in C^\infty (J) \times C^\infty (J; \varvec{\ell }^p(w))$ such that

• (i) there is a Hadamard state $\omega$ with associated moments $\mathcal {M}$ for the scale factor a,

• (ii) $(a, \omega )$ satisfies the SCE, viz., both trace equation (3.1) and the energy constraint (3.2),

• (iii) $\Vert [\Vert \big ]{a^{(j)}(\tau _\mathrm {free}) - b_j} \le \varepsilon$ for $j = 0,1,2$.

Proof

Let $\chi \in C^\infty ({{\mathbb {R}}}; [0,1])$ with ${{\,\mathrm{supp}\,}}\chi = [-1, 1]$ and $\chi (0) = 1$. For $\delta > 0$, define $\chi _\delta (x) = \chi (\frac{x}{\delta })$.

Fix $[\tau _\mathrm {tow}, \tau _\mathrm {fin}] \subset {{\mathbb {R}}}$ and $\tau _\mathrm {init}\in (\tau _\mathrm {tow}, \tau _\mathrm {fin})$. Let $a_\mathrm {tow}\in C^\infty [\tau _\mathrm {tow}, \tau _\mathrm {fin}]$ with $a_\mathrm {tow}> 0$ be initially Minkowskian (i.e., $a_\mathrm {tow}= 1$ in a neighborhood of $\tau _\mathrm {tow}$) and such that

$$\begin{aligned} a_\mathrm {tow}^{(j)}(\tau _\mathrm {init}) = b_j \quad \text {for j=0,1,2}. \end{aligned}$$

Then, the functions

$$\begin{aligned} a_{\mathrm {tow},c,\delta }(\tau ) = a_\mathrm {tow}(\tau ) + \frac{c}{2} \int _{\tau _\mathrm {tow}}^\tau (\tau - \eta )^2 \chi _\delta (\eta -\tau _\mathrm {init}) \mathop {}\!\mathrm {d}\eta , \end{aligned}$$

defined for $c \in {{\mathbb {R}}}$ and $\delta > 0$, satisfy

$$\begin{aligned} |[|\big ]{a_{\mathrm {tow},c,\delta }^{(j)}(\tau _\mathrm {init}) - b_j}&\le c \delta ^{3-j} \quad \text {for j=0,1,2}, \nonumber \\ a_{\mathrm {tow},c,\delta }'''(\tau _\mathrm {init})&= a_\mathrm {tow}'''(\tau _\mathrm {init}) + c. \end{aligned}$$

(5.14)

In particular, given any $c \in {{\mathbb {R}}}$, if $\delta > 0$ is sufficiently small, the conditions (5.13) hold with $b_0, b_1$ replaced by $a_{\mathrm {tow},c,\delta }(\tau _\mathrm {init})$ and $a_{\mathrm {tow},c,\delta }'(\tau _\mathrm {init})$.

Let $\omega _\mathrm {tow}$ be a translation-invariant Hadamard state on Minkowski spacetime such that the associated moments are (at most) geometrically growing (e.g., choose the vacuum state). Denote $(a_{c,\delta }, \mathcal {M}_{c,\delta })$ the solution obtained from Lemma 5.9 with $a_{\mathrm {tow},c,\delta }$ and initial moments $\mathcal {M}_\mathrm {tow}$ given by $\omega _\mathrm {tow}$. Observe that $\tau _\mathrm {free}$ and $\tau _\mathrm {stop}$ can be chosen uniformly for c from any fixed bounded interval.

Written as

$$\begin{aligned} g(\varvec{a}, \mathcal {M}, \tau ) = a'''(\tau ) - f_\mathrm {en}(a, a', a'', \mathcal {M}) \bigr |_\tau \end{aligned}$$

with the right-hand side as defined in Sect. 5.4, it is evident that the energy constraint is continuous in $\varvec{a}$ for $a > 0$ and $a' \ne 0$. Since V and $f_\mathrm {en}$ only depend on lower than third-order derivatives of the scale factor, it then follows from Theorem 4.7 and (5.14) that for any $\gamma > 0$ there are $c_\pm \in {{\mathbb {R}}}$ and $\delta > 0$ such that

$$\begin{aligned} \pm g(\varvec{a}_{c_\pm ,\delta }, \mathcal {M}_{c_\pm ,\delta }, \tau _\mathrm {init}) \ge \gamma \end{aligned}$$

It then follows by (iv) of Lemma 5.9 that one can choose $\tau _\mathrm {free}$ such that $\pm g(\varvec{a}_{c_\pm ,\delta }, \mathcal {M}_{c_\pm ,\delta }, \tau _\mathrm {free}) > 0$. Moreover, given sufficiently small $\delta > 0$, Theorem 5.6 implies the continuous dependence of $g(\varvec{a}_{c,\delta }, \mathcal {M}_{c,\delta }, \tau _\mathrm {free})$ on c. An application of the intermediate value theorem thus shows that there exists $c_0 \in (c_-,c_+)$ such that the associated solution $(a, \mathcal {M}) = (a_{c_0\delta }, \mathcal {M}_{c_0,\delta })$ satisfies the energy constraint $g(\varvec{a}, \mathcal {M}) = 0$ at $\tau _\mathrm {free}$.

As the energy constraint is propagated by solutions to the SCE, $(a, \mathcal {M})$ is a proper solution to the SCE (i.e., both the energy and the trace equation) in $[\tau _\mathrm {free}, \tau _\mathrm {stop}]$. As the Hadamard property is propagated on smooth spacetimes, $\mathcal {M}$ is associated with a Hadamard state, namely the state $\omega _\mathrm {tow}$ propagated on a. This shows (i) and (ii). Finally, (iii) follows from the triangle inequality using (5.14) and (iv) of Lemma 5.9. $\square$

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.

5.7 Reconstruction of the Quantum State

There is no obvious way to directly relate a sequence of moments $\mathcal {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 $\mathcal {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.

6 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) $-\lambda \langle {:}\phi ^4{:}\rangle _\omega$ would lead to quadratic terms in $\mathcal {M}_{\varphi \varphi }$ 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. Space-dependent 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 $\mathcal {G}-\mathcal {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 $\mathcal {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.

A Appendix

1.1 A.1 Homogeneous Distributions

Roughly following Chap. 3.2 of [32], we define for $z \in {{\mathbb {C}}}$ with ${\text {Re}}z > -1$ the function on ${{\mathbb {R}}}$

$$\begin{aligned} k_+^z :={\left\{ \begin{array}{ll} k^z &{} \hbox {if } k > 0, \\ 0 &{} \hbox {if } k \le 0. \end{array}\right. } \end{aligned}$$

(A.1)

Since this function is locally integrable, it defines a distribution. It can be extended to $z \in {{\mathbb {C}}}\setminus \{-1,-2,\dotsc \}$ by analytic continuation. To further extend $k_+^z$ to all $z \in {{\mathbb {C}}}$, we define for $n \in {{\mathbb {N}}}$ and any test function $f \in C^\infty _{\mathrm c}({{\mathbb {R}}})$

$$\begin{aligned} \langle k_+^{-n}, f \rangle :=\frac{1}{(n-1)!} \biggl ( -\int _0^\infty \log (k) f^{(n)}(k) \mathop {}\!\mathrm {d}k + f^{(n-1)}(0) \sum _{j=1}^{n-1} \frac{1}{j} \biggr ). \end{aligned}$$

(A.2)

Defined in this way, $k_+^z$ satisfies for $z \in {{\mathbb {C}}}$ the homogeneity property $\langle k_+^z, k f \rangle = \langle k_+^{z+1}, f \rangle$. Note that (A.2) is not the only possible extension of (A.1) to all of $z \in {{\mathbb {C}}}$—different extensions differ by derivatives of the delta distribution at zero.

To calculate the (inverse) Fourier transform of $k_+^{-n}$, note first that

$$\begin{aligned} \langle k_+^{-1}, \mathrm {e}^{\mathrm {i}k x} \rangle = \lim _{\nu \rightarrow 0} \int _0^\infty k^{\nu -1} \mathrm {e}^{\mathrm {i}k x} \mathop {}\!\mathrm {d}k - \frac{1}{\nu }. \end{aligned}$$

Then, we compute

$$\begin{aligned} \int _0^\infty k^{\nu -1} \mathrm {e}^{\mathrm {i}k x} \mathop {}\!\mathrm {d}k&= \Gamma (\nu ) (-\mathrm {i}x)^{-\nu } = \bigl ( \nu ^{-1} - \gamma + \mathcal {O}(\nu ) \bigr ) \bigl ( 1 - \nu \log (-\mathrm {i}x + 0) + \mathcal {O}(\nu ^2) \bigr ) \\&= \nu ^{-1} - \gamma - \log (-\mathrm {i}x + 0) + \mathcal {O}(\nu ) \end{aligned}$$

Therefore, we find

$$\begin{aligned} \langle k_+^{-1}, \mathrm {e}^{\mathrm {i}k x} \rangle = -\gamma - \log (-\mathrm {i}x + 0) = -\gamma - \log |x| + \frac{\mathrm {i}\uppi }{2} {{\,\mathrm{sgn}\,}}(x), \end{aligned}$$

which, by (A.2), immediately implies that, for $n \in {{\mathbb {N}}}$,

$$\begin{aligned} \langle k_+^{-n}, \mathrm {e}^{\mathrm {i}k x} \rangle = \frac{(\mathrm {i}x)^{n-1}}{\Gamma (n)} \left( \psi (n) - \log |x| + \frac{\mathrm {i}\uppi }{2} {{\,\mathrm{sgn}\,}}(x) \right) , \end{aligned}$$

(A.3)

where we used the relation between the harmonic numbers and the Digamma function

$$\begin{aligned} \psi (n) = \sum _{j=1}^{n-1} \frac{1}{j} - \gamma , \quad n \in {{\mathbb {N}}}. \end{aligned}$$

For completeness, let us also recall the definition of the Digamma function: $\psi (z) :=\Gamma '(z) / \Gamma (z)$ for $z \in {{\mathbb {C}}}\setminus \{0, -1, -2, \dotsc \}$.

This motivates the definition of the distributions (for $r \ge 0$)

$$\begin{aligned} h_z(r) :=\frac{\mathrm {e}^{\mathrm {i}z \uppi /2}}{2\uppi ^2} \frac{r^{z-2}}{\Gamma (z)} \bigl ( \log (r) - \psi (z) \bigr ), \end{aligned}$$

which extends analytically to $z \in {{\mathbb {C}}}$. In particular, we have for $n \in {{\mathbb {N}}}_0$

$$\begin{aligned} h_{-n}(r) = \frac{\mathrm {i}^n n!}{2\uppi ^2 r^{n+2}} \end{aligned}$$

because both $1/\Gamma (z)$ and $\psi (z)/\Gamma (z)$ are entire functions with

$$\begin{aligned} \frac{1}{\Gamma (-n)} = 0 \quad \text {and}\quad \frac{\psi (-n)}{\Gamma (-n)} = (-1)^{n+1} \Gamma (n+1). \end{aligned}$$

Taking the imaginary part of (A.3), we find

$$\begin{aligned} h_n(r) = \frac{1}{2\uppi ^2r} {\text {Im}}\langle k_+^{-n}, k \mathrm {e}^{\mathrm {i}k r} \rangle = \frac{1}{2\uppi ^2r} {\text {Im}}\langle k_+^{-n+1}, \sin (k r) \rangle \end{aligned}$$

for even n. Moreover, $h_z$ satisfies the important homogeneity property

$$\begin{aligned} -\Delta _r h_{z+2}(r) = h_z(r), \end{aligned}$$

(A.4)

where we recall that $\Delta _r$ denotes the (three dimensional) radial Laplacian.

1.2 A.2 Weighted Sequence Spaces

Let $w = (w_n)$ be a sequence of strictly positive numbers, called the weights. By $\ell ^p(w)$, $p \ge 1$, we denote the space of complex sequences $x = (x_n)$ with convergent norm

$$\begin{aligned} \Vert x\Vert _{p,w} :=\Bigl ( \sum _n {}|w_n^{-1} x_n|^p \Bigr )^{1/p}. \end{aligned}$$

If $p = \infty$, we denote by $\ell ^\infty (w)$ the space of complex sequences with convergent norm

$$\begin{aligned} \Vert x\Vert _{\infty ,w} :=\sup _n w_n^{-1} |x_n|. \end{aligned}$$

These are the weighted $\ell ^p$ spaces. If $w_n=1$, we omit the weight and denote by $\ell ^p$ with norm $\Vert \,\cdot \,\Vert _p$ the ordinary $\ell ^p$ spaces. Note that $\ell ^p(w)$ is reflexive for $1< p < \infty$.

Consider two weight sequences v, w. Then, it is easily seen that

$$\begin{aligned} \Vert x\Vert _{p,v} \le \sup _n \frac{w_n}{v_n} \Vert x\Vert _{p,w}. \end{aligned}$$

One of the most important operators on sequence spaces is the left-shift operator L, formally defined by $L (x_0, x_1, x_2, \dotsc ) = (x_1, x_2, \dotsc )$. If w are weights such that the sequence $(w_{n+1}/w_n)$ of ratios of consecutive weights is bounded, the left-shift operator is bounded on $\ell ^p(w)$. Indeed,

$$\begin{aligned} \Vert Lx\Vert _{p,w} \le \sup _n \frac{w_{n+1}}{w_n} \Vert x\Vert _{p,w}. \end{aligned}$$

(A.5)

However, if $(w_{n+1}/w_n)$ is unbounded, also the left-shift operator is unbounded.

More generally, for $m \in {{\mathbb {N}}}$ and two sequences v, w weights, we have

$$\begin{aligned} \Vert L^m x\Vert _{p,v} \le \sup _n \frac{w_{n+m}}{v_n} \Vert x\Vert _{p,w}. \end{aligned}$$

Finally, note that, if L is unbounded on $\ell ^p(w)$, its resolvent set is empty and its point spectrum fills the entire complex plane.

1.3 A.3 Some Inequalities

Lemma A.1

For $0 \le m \le n \in {{\mathbb {N}}}_0$ and $p \in (0,1)$,

$$\begin{aligned} \left( {\begin{array}{c}n\\ m\end{array}}\right) (1-p)^n \le \left( {\begin{array}{c}\lfloor m/p\rfloor \\ m\end{array}}\right) (1-p)^{\lfloor m/p\rfloor }, \end{aligned}$$

where $\lfloor \,\cdot \,\rfloor$ denotes the floor function.

Proof

The case $m = 0$ is obvious. For $m > 0$, we calculate the ratio of successive terms on the left-hand side:

$$\begin{aligned} \frac{\left( {\begin{array}{c}n+1\\ m\end{array}}\right) (1-p)^{n+1}}{\left( {\begin{array}{c}n\\ m\end{array}}\right) (1-p)^n} = \frac{n+1}{n+1-m} (1-p). \end{aligned}$$

Then, we observe that

$$\begin{aligned} \frac{n+1}{n+1-m} (1-p) \ge 1 \end{aligned}$$

if and only if $n+1 \le m/p$ to find the maximum at $n = \lfloor m/p\rfloor$. $\square$

Note that this lemma can also be stated in the language of probability theory: $n = \lfloor m/p\rfloor$ maximizes the probability of getting exactly m failures for a random variable following the binomial distribution with parameters n (number of trials) and p (probability of success).

An important inequality, accurately describing the asymptotics of the Gamma function, is Stirling’s inequality

$$\begin{aligned} \sqrt{2\uppi }\, x^{x+\frac{1}{2}} \mathrm {e}^{-x}< \Gamma (x+1) < x^{x+\frac{1}{2}} \mathrm {e}^{-x+1}, \quad x > 0. \end{aligned}$$

This double inequality can be improved in various ways, e.g., for $x \ge 1$,

$$\begin{aligned} \sqrt{2\uppi }\, x^{x+\frac{1}{2}} \mathrm {e}^{-x+(12x+1)^{-1}}< \Gamma (x+1) < \sqrt{2\uppi }\, x^{x+\frac{1}{2}} \mathrm {e}^{-x+(12x)^{-1}}, \end{aligned}$$

(A.6)

which can be obtained from [52], where also sharper bounds are presented.

Another useful inequality for the Gamma function is Gautschi’s inequality. For $x > 0$ and $s \in (0,1)$, we have

$$\begin{aligned} x^{1-s}< \frac{\Gamma (x+1)}{\Gamma (x+s)} < (x+1)^{1-s}, \end{aligned}$$

(A.7)

see, e.g., (5.6.4) of [43], which follows from the strict log-convexity of the Gamma function.

Combining the inequalities above, we find

Proposition A.2

For $m \in {{\mathbb {N}}}$ and $p \in (0,1)$,

$$\begin{aligned} \max _n \left( {\begin{array}{c}n+m\\ m\end{array}}\right) p^m (1-p)^n \le \min \biggl \{1,\frac{2}{\sqrt{2\uppi m}} (1-p)^{-\frac{1}{2}}\biggr \} \end{aligned}$$

(A.8)

and, as $m \rightarrow \infty$,

$$\begin{aligned} \max _n \left( {\begin{array}{c}n+m\\ m\end{array}}\right) p^m (1-p)^n \sim \frac{1}{\sqrt{2\uppi m}} (1-p)^{-\frac{1}{2}}. \end{aligned}$$

(A.9)

Proof

Denote by $\{x\} = x - \lfloor x\rfloor$ the fractional part of a real number x. We calculate

$$\begin{aligned} \max _n \left( {\begin{array}{c}n+m\\ m\end{array}}\right) p^m (1-p)^n&= \left( {\begin{array}{c}\lfloor m/p\rfloor \\ m\end{array}}\right) p^m (1-p)^{\lfloor m/p\rfloor -m} \\&\le \left( {\begin{array}{c}m/p\\ m\end{array}}\right) p^m (1-p)^{m/p-m} \biggl (1 + \frac{p}{m (1-p)}\biggr )^{\{m/p\}} \\&< \frac{1}{\sqrt{2\uppi m}} (1-p)^{-\frac{1}{2}} \biggl (1 + \frac{p}{m (1-p)}\biggr )^{\{m/p\}} \\&\le \frac{1}{\sqrt{2\uppi m}} (1-p)^{-\frac{1}{2}} \biggl (1 + \biggl \{\frac{m}{p}\biggr \} \frac{p}{m(1-p)} \biggr ), \end{aligned}$$

where we applied (in this order) Lemma A.1, Gautschi’s inequality (A.7), Stirling’s inequality (A.6) and Bernoulli’s inequality. Finally, an application of the inequality

$$\begin{aligned} \biggl \{\frac{m}{p}\biggr \} \frac{p}{m(1-p)} = \biggl \{\frac{m (1-p)}{p}\biggr \} \frac{p}{m(1-p)} \le \min \biggl \{1, \frac{p}{m(1-p)}\biggr \} \end{aligned}$$

yields (A.8), and from

$$\begin{aligned} \lim _{m \rightarrow \infty } \biggl \{\frac{m}{p}\biggr \} \frac{p}{m(1-p)} = 0 \end{aligned}$$

we obtain (A.9). $\square$

Proposition A.3

Let $n \in {{\mathbb {N}}}$. We have

$$\begin{aligned} \sum _{m=1}^{\lceil \frac{n+1}{2}\rceil } \frac{(m-1)!}{(n-m)!} \le {\left\{ \begin{array}{ll} \frac{3}{2} &{} \hbox {odd } n, \\ \frac{n}{2} + 1 &{} \hbox {even } n. \end{array}\right. } \end{aligned}$$

Proof

We consider only the odd case, the proof for the even case proceeds analogously. The smallest summand in the sum for n is $1/(n-1)!$, which is larger than the second smallest summand 1/n! in the sum for $n+2$. Therefore, we can bound all sums uniformly in n by summing up the smallest summands for each n. If we proceed like that for $n \ge 5$, we obtain

$$\begin{aligned} 1 + \frac{1}{6} + \sum _{n=2}^\infty \frac{1}{(2n)!} = \cosh 1 - \frac{1}{3} < \frac{3}{2}. \end{aligned}$$

For $n=1$ and $n=3$, the sums yield 1 and $\frac{3}{2}$, respectively. $\square$

1.4 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 $\sigma$. The world function $\sigma (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 neighborhood) and satisfies the relation

$$\begin{aligned} 2 \sigma = \bigl ((\nabla ^\mu \otimes \mathbb {1}) \sigma \bigr ) \bigl ((\nabla _\mu \otimes \mathbb {1}) \sigma \bigl ). \end{aligned}$$

(A.10)

In fact, the relation (A.10) together with the coincidence limits

$$\begin{aligned}{}[\sigma ] = 0, \quad [(\nabla _\mu \otimes \mathbb {1}) \sigma ] = 0, \quad [(\nabla _\mu \nabla _\nu \otimes \mathbb {1}) \sigma ] = g_{\mu \nu } \end{aligned}$$

(A.11)

uniquely defines Synge’s world function.

An expansion of Synge’s world function $\sigma (x,y)$ in terms of the coordinate distance $\delta 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)

$$\begin{aligned} \sigma (x,y) = \sum _n \frac{1}{n!} \varsigma _{\mu _1 \cdots \mu _n}(x) \delta x^{\mu _1} \delta x^{\mu _n}. \end{aligned}$$

As a consequence of (A.10) and (A.11), we find the recurrence relation

$$\begin{aligned} 2(1-n) \varsigma _{\mu _1 \cdots \mu _n}&= \sum _{j=2}^{n-2} \left( {\begin{array}{c}n\\ j\end{array}}\right) \, g^{\nu \rho } \bigl (\partial _\nu \varsigma _{(\mu _1 \cdots \mu _j|} - \varsigma _{(\mu _1 \cdots \mu _j| \nu }\bigr ) \\ {}&\quad \times \bigl (\partial _\rho \varsigma _{|\mu _{j+1} \cdots \mu _n)} - \varsigma _{|\mu _{j+1} \cdots \mu _n) \rho }\bigr )- 2n \partial _{(\mu _1} \varsigma _{\mu _2 \cdots \mu _n)} \end{aligned}$$

together with ‘initial’ coefficients $\varsigma = 0, \varsigma _\mu = 0, \varsigma _{\mu \nu } = g_{\mu \nu }$.

1.5 A.5 Some Technical Proofs

Proof of Theorem 5.1

For $\varvec{a} = (a_0, a_1, \dotsc , a_k) \in C(I; \mathcal {B}_R)$, define

$$\begin{aligned} \phi [\varvec{a}](\tau )&:=\bigl ( a_1(\tau ), a_2(\tau ), \dotsc , a_k(\tau ), f(\tau , \varvec{a}, \mathcal {M}[\varvec{a}]) \bigr ), \end{aligned}$$

(A.12a)

$$\begin{aligned} \Phi [\varvec{a}](\tau )&:=\varvec{a}_\mathrm {init}+ \int _{\tau _\mathrm {init}}^\tau \phi [\varvec{a}](\eta ) \mathop {}\!\mathrm {d}\eta . \end{aligned}$$

(A.12b)

By assumptions (ii) and (iii), we have

$$\begin{aligned} \Vert \Phi [\varvec{a}](\tau ) - \varvec{b}_c\Vert \le \Vert \varvec{a}_\mathrm {init}- \varvec{b}_c\Vert + \int _{\tau _\mathrm {init}}^\tau \Vert \phi [\varvec{a}](\eta )\Vert \mathop {}\!\mathrm {d}\eta \le \frac{R}{2} + (\tau - \tau _\mathrm {init}) (\Vert \varvec{b}_c\Vert + R + \mu _f), \end{aligned}$$

(A.13)

which proves that $\Phi$ closes on $C(J; \mathcal {B}_R)$, where $J = [\tau _\mathrm {init},\tau _\mathrm {stop}]$ for sufficiently small $\tau _\mathrm {stop} \le \tau _\mathrm {fin}$.

If $\varvec{a}, \tilde{\varvec{a}} \in C(J; \mathcal {B}_R)$, we calculate (using assumptions (iii) and (iv)) for $\tau \in J$

$$\begin{aligned} \Vert \Phi [\varvec{a}](\tau ) - \Phi [\tilde{\varvec{a}}](\tau )\Vert&\le \int _{\tau _\mathrm {init}}^\tau \Vert \phi [\varvec{a}](\eta )-\phi [\tilde{\varvec{a}}](\eta )\Vert \mathop {}\!\mathrm {d}\eta \\&\le \int _{\tau _\mathrm {init}}^\tau \bigl ( (1+L_f) \Vert \varvec{a}(\eta )-\tilde{\varvec{a}}(\eta )\Vert \\ {}&\quad + L_f {\Vert \mathcal {M}[\varvec{a}](\eta )-\mathcal {M}[\tilde{\varvec{a}}](\eta )\Vert }_{p,v} \bigr ) \mathop {}\!\mathrm {d}\eta \\&\le (1+L_f) (\tau _\mathrm {stop}-\tau _\mathrm {init}) {\Vert \varvec{a}-\tilde{\varvec{a}}\Vert }_{C(J; {{\mathbb {R}}}^{k+1})} \\&\quad + L_f L_{\mathcal M} \int _{\tau _\mathrm {init}}^\tau \int _{\tau _\mathrm {init}}^{\eta } \Vert \varvec{a}(\zeta )-\tilde{\varvec{a}}(\zeta )\Vert \mathop {}\!\mathrm {d}\zeta \mathop {}\!\mathrm {d}\eta \\&\le \bigl ( (1+L_f) (\tau _\mathrm {stop}-\tau _\mathrm {init})\\&\quad + \tfrac{1}{2} L_f L_{\mathcal M} (\tau _\mathrm {stop}-\tau _\mathrm {init})^2 \bigr ) {\Vert \varvec{a}-\tilde{\varvec{a}}\Vert }_{C(J; {{\mathbb {R}}}^{k+1})}. \end{aligned}$$

Therefore,

$$\begin{aligned} {\Vert \Phi [\varvec{a}] - \Phi [\tilde{\varvec{a}}]\Vert }_{C(J; {{\mathbb {R}}}^{k+1})}&\le \bigl ( (1+L_f) (\tau _\mathrm {stop}-\tau _\mathrm {init}) \\&\quad + \tfrac{1}{2} L_f L_{\mathcal M} (\tau _\mathrm {stop}-\tau _\mathrm {init})^2 \bigr ) {\Vert \varvec{a}-\tilde{\varvec{a}}\Vert }_{C(J; {{\mathbb {R}}}^{k+1})}. \end{aligned}$$

which shows that $\Phi$ is a contraction map if $\tau _\mathrm {stop}$ is sufficiently small.

It now follows by the Banach fixed-point theorem that $\Phi$ has a unique fixed point $\varvec{a} \in C^1(J; \mathcal {B}_R)$. By the form of (A.12), $\varvec{a}$ is such that $\varvec{a} = (a, a', \dotsc , a^{(k)})$ for some $a \in C^{k+1}(J)$. This yields the unique solution to (5.1) with the properties stated in the theorem. $\square$

Proof of Theorem 5.6

The two solutions $a(\tau )$ and $\tilde{a}(\tau )$ satisfy

$$\begin{aligned} \varvec{a}(\tau )&= \varvec{a}_\mathrm {init}+ \int _{\tau _\mathrm {init}}^\tau \bigl ( a'(\eta ), \dotsc , a^{(k)}(\eta ), f(\eta , \varvec{a}, \mathcal M[\varvec{a}]) \bigr ) \mathop {}\!\mathrm {d}\eta , \\ \tilde{\varvec{a}}(\tau )&= \tilde{\varvec{a}}_\mathrm {init}+ \int _{\tau _\mathrm {init}}^\tau \bigl ( \tilde{a}'(\eta ), \dotsc , \tilde{a}^{(k)}(\eta ), \tilde{f}(\eta , \tilde{\varvec{a}}, \tilde{\mathcal M}[\tilde{\varvec{a}}]) \bigr ) \mathop {}\!\mathrm {d}\eta \end{aligned}$$

and thus

$$\begin{aligned} \Vert \varvec{a}(\tau ) - \tilde{\varvec{a}}(\tau )\Vert \le \Vert \varvec{a}_\mathrm {init}- \tilde{\varvec{a}}_\mathrm {init}\Vert + \int _{\tau _\mathrm {init}}^\tau \bigl ( \Vert \varvec{a}(\eta ) - \tilde{\varvec{a}}(\eta )\Vert + |f(\eta , \varvec{a}, \mathcal {M}[\varvec{a}]) - \tilde{f}(\eta , \tilde{\varvec{a}}, \tilde{\mathcal M}[\tilde{\varvec{a}}])| \bigr ) \mathop {}\!\mathrm {d}\eta . \end{aligned}$$

Since

$$\begin{aligned}&|f(\tau , \varvec{a}, \mathcal {M}[\varvec{a}]) - \tilde{f}(\tau , \tilde{\varvec{a}}, \tilde{\mathcal M}[\tilde{\varvec{a}}])| \\&\quad \le |f(\tau , \varvec{a}, \mathcal {M}[\varvec{a}]) - \tilde{f}(\tau , \varvec{a}, \mathcal {M}[\varvec{a}])| + |\tilde{f}(\tau , \varvec{a}, \mathcal {M}[\varvec{a}]) - \tilde{f}(\tau , \tilde{\varvec{a}}, \tilde{\mathcal M}[\tilde{\varvec{a}}])| \\&\quad \le \mu _f' + L_f \bigl ( \Vert \varvec{a}(\tau ) - \tilde{\varvec{a}}(\tau )\Vert + \Vert \mathcal {M}[\varvec{a}](\tau ) - \tilde{\mathcal M}[\tilde{\varvec{a}}](\tau )\Vert _{p,v} \bigr ) \\&\quad \le \mu _f' + L_f \Vert \varvec{a}(\tau ) - \tilde{\varvec{a}}(\tau )\Vert + L_f L_\mathcal {M}' \left( \Vert \mathcal {M}_\mathrm {init}- \tilde{\mathcal M}_\mathrm {init}\Vert _{p,w} + \int _{\tau _\mathrm {init}}^\tau \Vert \varvec{a}(\eta ) - \tilde{\varvec{a}}(\eta )\Vert \mathop {}\!\mathrm {d}\eta \right) \end{aligned}$$

with the constants as defined in the theorem, we find

$$\begin{aligned}&\Vert \varvec{a}(\tau ) - \tilde{\varvec{a}}(\tau )\Vert - \Vert \varvec{a}_\mathrm {init}- \tilde{\varvec{a}}_\mathrm {init}\Vert \\&\quad \le \int _{\tau _\mathrm {init}}^\tau \biggl ( \mu _f' + L_f L_\mathcal {M}' \Vert \mathcal {M}_\mathrm {init}- \tilde{\mathcal M}_\mathrm {init}\Vert _{p,w} \\&\qquad + \bigl ((1 + L_f) + L_f L_\mathcal {M}' (\tau - \eta ) \bigr ) \Vert \varvec{a}(\eta ) - \tilde{\varvec{a}}(\eta )\Vert \biggr ) \mathop {}\!\mathrm {d}\eta . \end{aligned}$$

This implies inequality (5.6) by a generalization of Grönwall’s inequality (e.g., Theorem 1.4 of [8]). $\square$