Linear stability of semiclassical theories of gravity

The linearization of semiclassical theories of gravity is investigated in a toy model, consisting of a quantum scalar field in interaction with a second classical scalar field which plays the role of a classical background. This toy model mimics also the evolution induced by semiclassical Einstein equations, such as the one which describes the early universe in the cosmological case. The equations governing the dynamics of linear perturbations around simple exact solutions of this toy model are analyzed by constructing the corresponding retarded fundamental solutions, and by discussing the corresponding initial value problem. It is shown that, if the quantum field which drives the back-reaction to the classical background is massive, then there are choices of the renormalization parameters for which the linear perturbations with compact spatial support decay polynomially in time for large times, thus indicating stability of the underlying semiclassical solution.


Introduction
In semiclassical theories of gravity, the back-reaction of quantum matter fields is studied by means of the so-called semiclassical Einstein equations, where the Einstein tensor of a four-dimensional spacetime (M, g) is equated to the expectation value of the stress-energy tensor of matter fields evaluated on a quantum state ω. Any couple formed by a spacetime and a quantum state satisfying these equations constitutes a solution in semiclassical gravity. The question about existence and uniqueness of solutions was analyzed in several physical scenarios and with different approaches, which take advantage of recent developments in the study of locally covariant quantum field theories on globally hyperbolic spacetimes [DFP08, Pin11, PS15, Hac16, JA19, JA21, San21, GS21, MPS21, GRS22a, GRS22b, JAM22] (see also [MPRZ21] for an application of the semiclassical Einstein equations in black hole evaporation). In this framework, the renormalization of quadratic fields like the Wick square or the stress-energy tensor T µν entering the semiclassical models is always guaranteed when sufficiently regular states are taken into account. This is the case of Hadamard states, which fulfil the microlocal spectrum condition and possess two-point functions which have an universal singular structure [Rad96,BFK96]. Hence, in Hadamard states Wick observables can be constructed by a covariant point-splitting regularisation which removes the universal singular part, and generalises the usual normal-ordering prescription. For further details about this topic we refer to [Wal77, KW91, BF00, HW01, HW02, HW05] (see also [HW15] for a recent review and some physical applications).
The study of the stability of solutions of semiclassical equations remains problematic even on this class of states and even at linear order, because of the presence of higher-order derivative terms in the expectation value of the quantum stress-energy tensor. In the past years, the issue about stability of solutions of the semiclassical Einstein equations was addressed in several works [HW78,Hor80,Kay81,Yam82,Jor87,Sue89,MW20]. It was argued that, because of the higher order derivative terms, linearized semiclassical Einstein equations around chosen backgrounds admit exponentially growing solutions. These exponentially growing linear perturbations are called runaway solutions and indicate that the chosen background is unstable. It is remarkable that runaway solutions are present even on flat spacetimes. A prescription of reduction of order was presented in [Sim91,PS93,FW96] to eliminate runaway unstable solutions, whereas a criterion for the validity of semiclassical gravity in the linear regime was proposed in [AMPM03]. More recently, the stability of semiclassical solutions has been treated in the framework of the so-called stochastic gravity (see [HV20] and references therein). The irregularity issues of the semiclassical Einstein equations were also deeply studied in [MPS21] in the case of cosmological spacetimes, for arbitrary values of the coupling-to-curvature parameter. In this case, higher order derivatives of the metric appear, and furthermore the expectation value of the traced stress-energy tensor contains a non local quantum contribution at the linear order in the perturbative potential. This term represents the source of instability of the model, and it forbids to solve the semiclassical equations in a direct way. However, it was proved that unique local solutions exist after applying to the semiclassical equation an inversion formula associated to that unbounded operator.
In this paper, we shall take inspiration from the ideas presented in [RDKK80] and [JAMS20], and we analyse the stability issue at linear order in a simple semiclassical toy model in flat spacetime, consisting of a classical background scalar field ψ coupled to a quantum free scalar field φ. The system of equations governing the dynamics is displayed in eq. (2), and the back-reaction of the quantum field φ on ψ is estimated by substituting the classical field φ 2 with φ 2 ω , the expectation value of the quantum Wick square in the quantum state ω. With this picture in mind, we are interested in analyzing the stability of the solutions of this semiclassical system against linear perturbations. To this end, we discuss the equation obtained by linearizing the semiclassical equations over full solutions, which are formed by a quantum state ω for the quantum field φ and by a classical background ψ 0 for the classical field ψ. For the sake of simplicity, in this toy model, the state which is chosen is the Minkowski vacuum, however, similar results hold with other choices for ω. The obtained equation (9) is a linear equation for the perturbations ψ 1 over the background classical field ψ 0 , and in some cases it gives origin to a well posed initial value problem for spatially compact solutions despite the presence of certain unavoidable non local contributions in it. These non local contributions are manifest in the linear equation for ψ 1 written in eq. (23), or in the corresponding non-homogeneous equation (24) obtained when a smooth compactly supported source f for ψ 1 is considered.
In a first step, we prove that this equation is of hyperbolic nature. We then explicitly construct the retarded fundamental solutions D R . Out of it we write the most general retarded solution of eq. (24), we discuss the well posedness of the corresponding initial value problem, and we study the decay of these solutions for large times.
Notably, there are several choices of the parameters governing the dynamics for which every linearized solution having smooth compactly supported initial values, or emerging from a compactly supported source, decays polynomially in time for large times, thus showing that perturbations tend to disperse, or better to disappear in time. This result of stability is ensured by both the non-vanishing mass of the quantum field, by the spatial support of the initial data and the source. On the contrary, in agreement with previous observations [HW78,Hor80,Jor87], if the quantum field φ is massless, then solutions of the linearized semiclassical equation which grow exponentially in time may always exist, even if the initial values and the source are of compact support.
The toy model studied in this paper formally mimics both the cosmological semiclassical model studied in [MPS21] and the weak field theory discussed in [Hor80], after interpreting the background solution as a degree of freedom of the metric, and the coupling constants as the renormalization freedoms appearing in a semiclassical theory of gravity. This analogy indicates that such semiclassical theories may have stable solutions at least when the quantum fields describing matter are massive. This paper is organized as follows, in the Section 2 we introduce the interacting toy model considered in this paper, and we discuss how to obtain a linearized semiclassical equation governing the back-reaction of the quantum field on the background. In the first part of the Section 3, we analyze the linearized equation, we construct its retarded fundamental solution and we show how to formulate a well-posed initial value problem out of it. In the last part of this section, we present a correspondence between the presented toy model and the semiclassical cosmological problem. Finally, the main results of the paper are summarized in Section 4. Appendix A contains some technical results about the decay at large time of certain functions.

Notations
In this paper, the units convention is G = c = 1, and the Lorentzian signature of the spacetime (M, g) is (−, +, +, +). Thus, the d'Alembert operator reads where ∆ x denotes the spatial Laplace-Beltrami operator. We shall employ the following conventions on the Fourier transforms: respectively. Hence, the convolution theorems read as F{f * g} = F{f }F{g} where * denotes the convolution operator such that f * g = f (x − y)g(y)dy.

The semiclassical model
In this work, we study the coupling between a quantum scalar field φ and another classical background scalar field ψ in the Minkowski spacetime (M, η). The equations of motion of the corresponding free theory are where g 1 , g 2 , λ, λ 1 , λ 2 denote the real coupling constants of the theory. With the choice of λ 2 = 0, the system is Lagrangian. On the other hand, if λ 2 = 0, then there is sufficient freedom in the definition of the coupling constants to fix λ 2 = 1. However, we do rather not impose any further constraint on the coupling constants, in order to describe as many semiclassical theories as possible. In particular, keeping λ 2 = 0, eq. (1b) mimics the form of equations arising in the study of semiclassical Einstein equations, as we shall discuss in subsection 3.4. The first equation (1a) is the equation of motion of the linear, Klein-Gordon like massive field φ, where m is the mass, and λψ plays the role of an external potential. The quantization of φ follows straightforwardly once the external potential field ψ is known. To this avail, we shall consider the unital * −algebra of field observables A generated by the identity I and the abstract smeared field φ(f ), with f ∈ C ∞ 0 (M) any compactly supported smooth function [BDFY15,Haa12]. The product in this * −algebra encodes the canonical commutation relations given in terms of the causal propagator ∆ . = ∆ R − ∆ A , i.e., the difference of the retarded and advanced propagators uniquely obtained as fundamental solutions of eq. (1a). The algebra of field observables can be enlarged to contain also Wick powers like φ 2 when a properly defined normal ordering procedure is taken into account [HW01,HW05,Mor03].
In light of the quantum nature of φ, eq. (1b) can be interpreted in the semiclassical approximation, namely by taking the expectation values of φ 2 in a suitable quantum state. More precisely, once a state ω on A is chosen for the quantum field theory, we have that = ω(φ 2 ) is the expectation value of the properly normal ordered quantum field φ 2 in the quantum state ω. Thus, the semiclassical system corresponding to equations (1a) (1b) turns out to be (2) To analyze the linearization of this system around a given simple solution of eq. (2) determined by (ψ 0 , ω), we decompose the classical field ψ in two parts, the background contribution ψ 0 plus a perturbation ψ 1 , so that ψ = ψ 0 + ψ 1 .
The state ω, or more specifically its two-point function ∆ +,ω appearing in the second equation of the system (2) through φ 2 ω , can also be decomposed in two parts i.e., the background contribution plus its perturbation. Notice that the linear order contribution in ∆ +,ω,1 is formed by two terms: one which takes into account the modified evolution equation induced by ψ 1 and satisfied by ∆ +,ω , and one which consists in w 1 , a symmetric bi-solution of the zeroth order equation of motion satisfied by φ. Since ∆ +,ω,0 is a bi-solution of the same equation, the effect due to w 1 can be reabsorbed in a redefinition of the background theory. For this reason, in the following, we shall take into account explicitly the effect due to the modified evolution only. To control the linear contribution in ∆ +,ω,1 due to the change of dynamics, we shall use perturbation theory for quantum fields. As we shall see later, if instead one takes into account the effect of w 1 explicitly, an inhomogenous (source) term has to be added to the linearized equation satisfied by ψ 1 . In fact, this extra term will not alter the discussion we are going to present. The quantization of φ is performed on the fixed background ψ 0 by considering the * −algebra of field observables A with a product that implements the canonical commutation relations emerging from to the zeroth-order equation Furthermore, both the quantum state ω and the background field ψ 0 satisfy the semiclassical equation, namely the second equation in the system given in eq. (2), which thus represents a constraint for the couple (ψ 0 , ω). We recall here that the expectation value of φ 2 in the state ω on this background theory is obtained by an ordinary point-splitting regularisation procedure where ∆ +,ω is the two-point function of the state ω, H(x, y) . = u(x, y)/(x−y) 2 +v(x, y) log (x − y) 2 /µ 2 is the universal divergent contribution present in any Hadamard two-point function, with µ a fixed length scale. Furthermore, c is a constant which encodes the regularisation freedom present in the construction of Wick powers like φ 2 [HW01,HW02].
Notice that the influence of ψ 1 on φ is governed by the first equation in the system (2); for a given ψ 1 , this equation descends from a Lagrangian for φ. In this Lagrangian, ψ 1 acts as a mass perturbation for φ, hence we may use Lagrangian methods to analyze the influence of ψ 1 to φ even if the full system formed by both equations in eq. (2) is not Lagrangian. Thanks to the principle of perturbative agreement (see e.g. [HW05,DTPP]) this approach is equivalent to directly analyze the effect of ψ 1 on the two-point function, as it was done, e.g., in [Pin11,PS15,MPS21], or to evaluate φ 2 ω in [EG11] in cosmological spacetimes. Thus, we pass to analyze the influence of the perturbation ψ 1 on φ by means of perturbation theory considering the interaction Lagrangian A perturbative construction of interacting quantum field theories on a generic curved spacetime (M, g) was rigorously formulated in a local and covariant way in the framework of perturbative algebraic quantum field theorycf. [BF00, DF01, HW01, HW02, FL14, FR16, GHP16] to which we refer for further details on the construction we are going to present. In particular, the expectation value of the interacting field φ 2 in the state ω is obtained by means of the Bogoliubov map where the perturbative potential is obtained by smearing the interaction Lagrangian given in eq. (5) with a smooth cut-off f ∈ C ∞ 0 (M) which is equal to 1 on the compact spacetime region where we want to test the semiclassical equation. This cut-off f is eventually removed by considering a suitable limit in which f tends to 1 on every point of M.
The Bogoliubov map R V is used to represent local field observables of the interacting theory as formal power series in λ, whose coefficients are well defined elements of A, the extended algebra of free fields. In particular, where T is the map which realizes the time ordering, while S(V ) is the time ordered exponential of the smeared interaction Lagrangian V , namely For the precise construction of the time ordering map T we refer to [BF00,HW01,HW02], where the old construction of Epstein and Glaser in [EG73] is generalized to a generic curved spacetime.
With the perturbative construction of φ 2 ω at disposal, the perturbative expansion of φ 2 ω in the interacting theory obtained by means of the Bogoliubov map reads at the linear order in V as where the contributions higher than the linear one are not displayed. Here, The factor i at the right hand side of the second equation in eq. (7) is present because S(V ) in eq. (6) is formally unitary. As expected by the principle of perturbative agreement, one can notice by direct computation that φ 2 (lin) ω matches the linear contribution obtained in [Pin11,PS15,GS21,MPS21] in the cosmological case.
The linearization studied in this paper consists of studying the semiclassical theory described by the system of equations given in eq. (2), where the expectation value of the Wick square is approximated by truncating at first order the formal power series in the interaction Lagrangian (5) occurring in the map R V defined in eq. (6). The state for the interacting quantum theory is constructed by means of the state on the linear quantum theory, i.e., on A, and it is fixed once and forever, no matter the form of the linear perturbation ψ 1 .
Hence, on the one side the background theory is described by (ψ 0 , ω), in which ψ 0 fulfils the semiclassical equation and the the quantum state ω constrained by eq. (8) is fixed once and for all in the linear algebra A of fields satisfying eq. (3). On the other side, the linear perturbation theory of the background is described by the classical field ψ 1 , which fulfils the following linearized semiclassical equation where φ 2 (lin) ω was constructed in eq. (7). Eq. (9) is the only equation we have at linear order, and it must be seen as a dynamical equation for the linear perturbation ψ 1 . If a perturbation of the state which modifies ω to ω + ω 1 is considered explicitly at the linear order, then an inhomogoenous contribution appears in eq. (9), which consists in adding . This extra contribution does not depend on ψ 1 , and hence it acts as a source term in eq. (9). We shall see in the next section that also solutions of the inhomogeneous version of eq. (9) have nice decay properties for several choices of the parameters a, λ, λ i , g i .

Stability of solutions of the linearized semiclassical equation
With the results presented in [RDKK80] in mind, our goal in this section is to show that perturbations ψ 1 over (ψ 0 , ω) which solves eq. (2) decay for large time at linear order in perturbation theory and for proper choices of the coupling constants of the model. To achieve this, it is crucial to assume compactly supported initial data on the fields, and to consider perturbations around the background solution which have spatial compact support, in order to avoid the class of exponentially growing solutions already seen, e.g., in [HW78, Hor80, Jor87].

Zeroth-order solution
Before discussing the perturbative construction of φ, we give here some details on the chosen background solution (ψ 0 , ω). For the sake of simplicity, we shall assume that the background ψ 0 ∈ R is constant, hence the state ω must be such that the Wick square has constant expectation value In view of the renormalization freedom present in the definition of the Wick power expressed by the constant c in eq. (4), it is always possible to fulfil the previous equation whenever ω is a translation invariant state, for any choice of ψ 0 ∈ R. A constant background external field ψ 0 corresponds to a mass renormalization, that is, As m 2 λ denotes the new renormalized mass of the background field φ, we need to impose the constraint that m 2 λ ≥ 0, otherwise the reference state for the system may not exist. This inequality always holds for sufficiently small λ 1 /g 1 , and it holds trivially in the case of vanishing expectation value of the Wick square. Besides, there is always the possibility of setting φ 2 (bac) ω = 0 by means of the choice of the renormalization of the field φ 2 (the constant c in eq. (4)).
To simply further the analysis, we shall select the Minkowski vacuum state ω 0 as reference state on A, whose corresponding two-point function is Θ is the Heaviside step function, and δ the Dirac delta function. However, other choices of quantum states, such that φ 2 ω is regular in M, do not alter significantly our analysis.

The linearized expectation value of the Wick square
Using the Bogoliubov map given in eq. (6), the expectation value of the Wick square can be evaluated at every perturbation order. The linearized contribution defined in eq. (7) in the adiabatic limit (f = 1) takes the following form where ∆ F,ω (y, x) = T (φ(y)φ(x)) ω and ∆ +,ω (y, x) = φ(y)φ(x) ω are the Feynman propagator and the two-point function associated to ω, respectively. Notice that the definition of the Feynman propagator ∆ F,ω . = ∆ +,ω +i∆ A , where ∆ A is the advanced propagator, employed here differs by a factor i from others constructions. Furthermore, the squares of ∆ F,ω and ∆ +,ω correspond to the pointwise multiplication of the integral kernels of the two distributions. We shall discuss how to construct these products below. A diagrammatic representation of the propagators in the integrand at the right hand side of eq. (12) is given in Figure 1.
The various propagators of the theory, and in particular those appearing in eq. (12), are in general not invariant under translation. However, they acquire translation invariance when they are referred to the Minkowski vacuum ω 0 . To keep this in mind, we shall denote the Feynman propagator and the two-point function referred to the Minkowski vacuum as ∆ F (y − x) = 0|T (φ(y)φ(x))|0 and ∆ + (y − x) = 0|φ(y)φ(x)|0 , respectively. Furthermore, we denote by W . = ∆ +,ω − ∆ + and we observe that W is a smooth function whenever ω is an Hadamard state, namely the singular part of two-point function ∆ +,ω is the same as the one in the Hadamard parametrix, or, equivalently, ω fulfils the microlocal spectrum condition [Rad96,BF00]. Hence, recalling that ∆ F = ∆ + + i∆ A , we get is represented by a non-oriented line because it is symmetric in the exchange of x ↔ y, while the two-point function Since W is smooth, we just need to construct ∆ 2 + and ∆ 2 F to give meaning to the right hand side of eq. (12). Contrary to ∆ 2 + which is well defined everywhere, ∆ 2 F is a well defined distributions only for test functions which are not supported on the origin. Therefore, a renormalization (extension) procedure is required to have a well defined ∆ 2 F also in this case. The extension of ∆ 2 F on test functions supported on the origin can be obtained keeping fixed the Steinmann scaling degree [Ste71,BF00]; however, this extension is not unique, and the remaining freedom amounts to an additional c 0 δ, where c 0 is a real parameter. This freedom is compatible with the ambiguity present in the definition of φ 2 [HW01]. To obtain explicit expressions of ∆ 2 F and ∆ 2 + , we make use of the known Källen-Lehmann spectral representations. In particular, using the convolution theorem and the definition of ∆ + (x) given in eq. (11), we obtain that where the spectral density ( is the Minkowski vacuum two-point function given in eq. (11) with mass M . On the other hand, the very same representation for ∆ 2 F (x) cannot be given because the integral over M 2 is divergent in that case. However, using the fact that for This expression coincides with ∆ 2 + (x) for x / ∈ J − (0) and with ∆ 2 + (−x) for x / ∈ J + (0). Moreover, it is well defined also on functions whose support contains 0, and thus it represents one of the possible extensions of ∆ 2 F ∈ D (M \ {0}) to D (M). Also, the difference of two ∆ 2 F constructed with a and a = a is a non vanishing distribution supported in the origin with scaling degree 4, and hence it must be proportional to the Dirac delta function. Finally, we can saturate the freedom in the construction of ∆ 2 F with various choices of a. In light of this, the constant a encodes the renormalization freedom present in the construction of ∆ 2 F (x). Therefore, recalling that where ∆ A (x, M 2 ) is the advanced propagator of the Klein-Gordon field of mass M . A detailed derivation of eq. (13) for the massless case can be found in Appendix C of [DF04]. Hence, eq. (12) can be written at the linear order outside x = 0 as where the operator K a maps compactly supported smooth function (or Schwartz function) to smooth functions. Moreover, its regularized integral kernel takes the form where the d'Alembert operator is taken in the distributional sense, and ∆ R (·, M 2 ) is the retarded propagator of the Klein Gordon equation with mass M . Its spatial Fourier kernel reads as where ω 0 = |p| 2 + M 2 , and Θ is the Heaviside step function. The operator W maps compactly supported smooth functions to smooth functions, and its integral kernel is the pointwise multiplication of the advanced propagator with the smooth part of the two-point function W , i.e., W . = 2∆ A W . With the choice of the Minkowski vacuum as reference state, both W and W vanish, and thus the linearized expectation value of the Wick square given in eq. (14) simplifies as For later purposes, we need to control the evolution of φ 2 (lin) 0 in time under the influence of ψ 1 : to this end, we study the kernel given in eq. (15) in the Fourier domain by means of the following proposition.
Proposition 3.1. Let ψ 1 ∈ S(M) be a Schwartz function on M, and letψ 1 be its Fourier transform. Then the Fourier transform of the linearized expectation value of the Wick square given in eq. (17) can be written as given for strictly positive mass m > 0 and for a > −4m 2 . The function F a (z) admits the following integral representation: and it has the following properties: a) F a (z) is analytic for z ∈ C \ (−∞, −4m 2 ] and continuous at z = −4m 2 ; b) the domain F a (z) has a branch cut on z ∈ (−∞, −4m 2 ) because there the imaginary part is discontinuous (the real part is continuous but not differentiable); it is strictly increasing for s ∈ −4m 2 , ∞ , and it diverges for large |s|; e) The imaginary part of F a admits the following integral representation: it is strictly positive for Im(z) > 0, and strictly negative for Im(z) < 0. Furthermore, it vanishes for z ∈ (−4m 2 , ∞), and it is discontinuous on z ∈ (−∞, −4m 2 ) (the absolute value is finite).
The second graph is the qualitative behavior of |ImFa(x)|.
Proof. Using the definition of K a in eq. (15), the Fourier transform of the distribution ∆ R given in eq.
Hence taking the Fourier transform on both sides of eq. (17) and applying the convolution theorem to K a (ψ 1 ) yields whereψ 1 (p 0 , p) ∈ S(M). Then, the first part of the thesis follows after recalling the form of 2 (M 2 ) = 1 16π 2 1 − 4m 2 M 2 and the definition of F a (z). The list of properties of F a (z) can be inferred directly from its integral representation. To check the validity of the representation (20) of F a , consider From the expression of F a given in eq. (20), F a (−p 2 0 + |p| 2 ) = (p 2 0 − |p| 2 + a)A a (−p 2 0 + |p| 2 ). We take the inverse Fourier transform in time of A a (−p 2 0 + |p| 2 ), that is, The integrand inÃ has two cuts for (p 0 − i ) 2 > (|p| 2 + 4m 2 ) located in the upper half complex plane and it is analytic outside the two cuts (it has no poles); furthermore, |A(−(w − i ) 2 + p 2 )| for w ∈ C vanishes in the limit |w| → ∞. Hence, that inverse Fourier transform can be obtained by standard results of complex analysis, including Jordan's lemma and Cauchy residue theorem. In particular, to evaluate the integral over the real line, for t < 0 we can close the contour in the lower half plane, and thusÃ = 0 because A(−(w − i ) 2 + p 2 ) is analytic in the lower half plane. On the other hand, if t > 0, then the contour is closed in the upper half plane, and thus only the two cuts matter in the evaluation of the integral over the real line which givesÃ. The contributions due to the two cuts for t > 0 can be combined to giveÃ Changing variable of integration to M 2 = w 2 − |p| 2 and computing the discontinuity of log along its cut, we obtain for t > 0 that and hence ( + a)A is equal to K a up to a constant factor given in eq. (15). Therefore, the expression of F a given in eq. (19) follows, thus proving that it coincides with eq. (20).
Remark 3.2. In the limit of vanishing mass m 2 = 0 and for −p 2 0 + |p| 2 > 0, the function F a (−p 2 0 + |p| 2 ) given in eq. (19) takes the form This logarithmic behaviour is similar to the case studied in [Hor80] for the linearized semiclassical Einstein equations, in the weak-field limit of gravity and after considering massless quantum fields (see also [HH81]). The function F a is also similar to the Fourier transform of the Green function associated to the first order equation analyzed in [FW96] in the context of semiclassical back-reaction.

Linearized solutions
The semiclassical equation (9) governing the dynamics of the perturbations at the linear order is a linear equation in the perturbation field ψ 1 . As we shall see below the properties of the solution space of that linear equation depends strictly on the parameters a, λ, λ i , g i , i = 1, 2. The non local state-dependent contribution φ 2 (lin) 0 (x) defined in eq. (17) was constructed in terms of the linear operator K a introduced in eq. (14), and thus it was expressed in terms of the function F a studied in Proposition 3.1. Therefore, eq. (9) evaluated in the Minkowski vacuum state reads To highlight its mathematical structure, eq. (22) can be rewritten in the following form: where P λ . = λ 2 − λ 1 , and P g . = g 2 − g 1 .
Because of the presence of a second d'Alembert operator in the expression of φ 2 (lin) 0 (x) given in eq. (14) through eq. (15), eq. (22) contains fourth-order derivatives in ψ 1 . It has thus a form similar to the semiclassical equations which usually appear in semiclassical theories of gravity, see, e.g., [Kay81] and the next Section. Thus, it manifests the same conceptual issues already known in semiclassical gravity as higher-order theory of gravity. In particular, there are cases where similar equations admit the so-called runaway solutions, which make the classical background field unstable in the semiclassical approach [HW78,Hor80,Jor87,PS93,Sim91,FW96].
Runaway solutions usually consist of solutions of the linearized system around some background which grow exponentially in time. This class of linearized solutions become dominant over the background at large times. Thus, the full solution of the system acquires, in principle, a very different form from the chosen background, and at the same time it is expected to be very sensitive to the chosen initial conditions. Therefore, the background solution cannot be assumed to be stable. On the contrary, if all the linearized solutions decay sufficiently fast to zero for large times, then the perturbations become negligible with respect to the background solution, thus indicating the stability of the background.
In the next part we analyze eq. (23) equipped with a compactly supported smooth source term, namely where K a is the linear operator introduced in eq. (14), and f ∈ C ∞ 0 (M) is a compactly supported source. The strategy is as follows. In a first step, we shall show that this equation manifests an hyperbolic nature, and we shall construct its retarded fundamental solutions as an operator D R : C ∞ 0 (M) → C ∞ (M). Afterwards, thanks to the regularity properties of D R , we shall prove that past compact solutions of the form ψ 1 = D R (f ) decay as 1/t 3/2 for large times t, hence getting the stability of the corresponding backgrounds against perturbation sourced by f .
In a second step, we shall study the smooth spatially compact solutions of the homogeneous equation (23) corresponding to eq. (9). To determine uniquely a solution in the future and in the past of t 0 , we shall equip eq. (23) with suitable initial conditions at t 0 = 0, i.e., smooth compactly supported initial data of the form ψ (j) We shall see that the number of initial conditions necessary to determine a spatially compact solution depends on the choice of the parameters: in some cases, four initial conditions have to be imposed, while, in other cases, only two initial conditions are sufficient to obtain a solution; finally, there are also cases where no solution exists. Thus, we shall prove that there are wide ranges of values of (a, g 1 , g 2 , λ, λ 1 , λ 2 ) for which solutions of eq. (23) with compactly supported initial data decay faster than 1/t 3/2 for large times. Therefore, the stability of the linearized back-reacted system is restored even in this case.
Our analysis starts from showing that eq. (22) written as eq. (23) manifests an hyperbolic nature.
Proof. As the solution ψ 1 is past compact by hypothesis, we can apply the retarded operator associated to P λ on both sides of eq. (24), thus obtaining λK a (ψ 1 ) + ∆ R,λ P g ψ 1 = ∆ R,λ f.
The constraints on the parameters λ 2 and g 2 given by hypothesis in Proposition 3.3 can be easily removed adapting the first part of the proof. For example, if g 2 = 0, there is no need to add P λ ∆ R,λ in the right hand side of eq. (25). On the other hand, if λ 2 = 0, then the proof starts with getting an analog of eq. (25) after applying the retarded operator ∆ R,g at the place of ∆ R,λ on both sides of eq. (24). Finally, if both λ 2 and g 2 vanish, then there is no need of preliminary applying any retarded operator to eq. (24). Proposition 3.3, and more precisely eq. (25), suggests that the form of a past compact solution ψ 1 of eq. (24) in x cannot be influenced by any modification of ψ 1 or f outside of J − (x). We shall see a posteriori that this indication is actually correct, because we shall prove in Theorem 3.5 that a retarded fundamental solution of eq. (24) exists.
We proceed in the analysis of the form of the solution of the linearized semiclassical equation by studying the associated retarded fundamental solution. We use Fourier techniques to analyze the fundamental solution of eq. (24). Hence, where f ∈ C ∞ 0 (M). Eq. (24) can equivalently be written in a more compact form as In order to obtain the retarded fundamental solutions associated to the kernel S(z), we need to study the set of points in the complex plane in which S(z) vanishes: we denote this set by S. Then, we shall prove that, if the parameters (λ, g i , λ i ) satisfy certain conditions, this set contains only real elements, and, furthermore, it includes only negative elements in some special cases. Among them, we shall impose as a constraint on the parameters the following inequality: The characterization of elements in S is studied in the following proposition.
Proposition 3.4. Let S ⊂ C be set of zeros of S(z) given in eq. (27). Fix the parameters in such a way that at least one of the two λ i , i ∈ {1, 2} is non vanishing and g 2 λ 1 − λ 2 g 1 ≥ 0, λ > 0, −4m 2 < a. Then we distinguish two cases: In particular, if λ 2 = 0, −λ 1 /λ 2 is in S only if g 2 λ 1 /λ 2 = −g 1 . Furthermore, S contains one, two or no elements depending on the parameters λ i , g i , and a.
To prove that the solution set is as stated in item a) and b), we proceed as follows. After multiplying both sides of the equation by (λ 1 + λ 2 z), taking the imaginary part yields the following equation: where the global sign of right-hand side depends only on Im(z), because by hypothesis g 2 λ 1 − λ 2 g 1 ≥ 0. In particular, for strictly positive Im(z) the right hand side is negative, while the left hand side is strictly positive thanks to item e) of Proposition 3.1; similarly, for strictly negative Im(z) the right hand side is positive, while the left hand side is strictly negative. Moreover, as Im(F a (x)) is not 0 also for x < −4m 2 (see Figure 2), we have that the only possible solutions of eq. (29) needs to be searched within (−4m 2 , ∞) ∪ {−λ 1 /λ 2 } when λ 2 = 0, or in (−4m 2 , ∞) otherwise. Furthermore, if both λ 2 = 0 and z = −λ 1 /λ 2 , then the left hand side of eq. (29) vanishes, whereas the right hand side vanishes only when g 2 λ 1 /λ 2 = −g 1 .
We are thus left with the analysis the following real equation: The properties of the function F a (s) for s ∈ (−4m 2 , ∞) ⊂ R are listed in Proposition 3.1, and the plot of this function is reported in Figure 2. Notice in particular that F a (s) is a concave function, because the second derivative of the integrand in eq. (19) is −2 (M 2 )(M 2 + s) −3 , which is a strictly negative integrable function, and thus F a (s) is strictly negative. Hence, for λ 2 = 0, λ 1 = 0 by hypothesis, (λ 1 + λ 2 s)F a (s) has a definite concavity. In this case, the maximum number of distinct solutions of eq. (30) is two. If λ 2 = 0, we rewrite eq. (30) as Notice that λ 16π 2 F a (s) + g2 λ2 is monotonically increasing. Using the hypothesis that λ 2 g 1 − λ 1 g 2 ≥ 0, we observe that I(s) is constant if λ 2 g 1 − λ 1 g 2 = 0 or monotonically decreasing if λ 2 g 1 − λ 1 g 2 > 0; in this latter case, it has also a discontinuity (a vertical asymptote) in s = −λ 1 /λ 2 . Hence, also in this case there are at most two solutions, thus concluding the proof. Taking into account all the previous statements, we are now ready to write the explicit form of the retarded fundamental solution of eq. (24), and to show that past compact solutions decay to zero for sufficiently large times, as expected for a perturbation over a stable background.
Fix as non-vanishing constants at least one of the two g i , and at least one of the λ i , assume that the inequality g 2 λ 1 − λ 2 g 1 ≥ 0 holds, and set −4m 2 < a < 0. Suppose that the set S defined in Proposition 3.4 contains only real negative elements, then the Fourier transform of the retarded fundamental solution D R of eq. (24) readsD where S(z) was defined in eq. (27). Hence where the elements of S are the zeros of S(s), with s ∈ (−4m 2 , ∞) ∪ {−λ 1 /λ 2 }. The retarded fundamental solution D R is a linear operator which maps smooth compactly supported functions to smooth functions, and with this operator at disposal the solution ψ 1 of eq. (24) with past compact support is For λ 2 = 0, ψ 1 (t, x) decays as 1/t 3/2 for large t.
Proof. We analyze the equation (24) in the Fourier domain. Using the results given in Proposition 3.1, it takes the form of where S was defined in eq. (27). Thus, the Fourier transform of the retarded operator D R yieldŝ and hence its Fourier inverse transform can be evaluated by means of standard methods of complex analysis. In view of the properties of the function F a (z) given in eq. (19), the function h(z) .
For |z| > R, it holds that |h(z)| < l(|z|), where the function l(r) vanishes in the limit of large positive r, because, in the worse case, 1/S(−z 2 + |p| 2 ) is dominated by c/F a (−z 2 + |p| 2 ) for large |z|, for some constant c, and |F a (−z 2 + |p| 2 )| grows as log |z| for large |z|. Therefore, from Jordan's lemma we may close the contour γ in the upper or lower plane, according to the sign of t, with a semicircle which does not contribute to the integral in the limit R → ∞.
The function 1/S(−(z − i ) 2 + |p| 2 ) has two poles for each element s ∈ S (the set of zeros of S(z)). Since s ∈ S is negative by hypothesis, we have that the poles are located on the line Im(z) = i , and correspond to the complex numbers z = i ± |p| − s.
Furthermore, the function 1/S(−(z − i ) 2 + |p| 2 ) has two branch cuts located at z = x + i , where x 2 ≥ |p| 2 + 4m 2 . Thus, 1/S is analytic in the lower half plane, and hence, we obtain thatD R = 0 for t < 0 by Jordan's lemma, because we close the contour in the lower half plane for t < 0.
On the other hand, if t > 0, then we close the contour γ in the upper half plane, and hence we need to take care of both the poles and the branch cuts. In this case, if we deform the previous contour γ to a newγ in such a way to avoid both the poles and the cuts, then the result of the contour integral over γ vanishes. Therefore, the only two non-vanishing contributions inD R , denoted byÕ andC, are due to the poles and the branch cuts, respectively.
The contribution due to the poles can be directly evaluated using the Cauchy residue theorem, which yieldsÕ where w s . = |p| 2 − s, and hence, in view of eq. (16), The contribution due to the cuts can be combined in the following form Thus, recalling eq. (16) again, it can be written in the position domain as The discontinuity in the two cuts is only due to the imaginary part of F a (z), hence thus getting eq. (31) by combining eqs. (34) and (35). Furthermore, the integral over M 2 present in C(x) in eq. (35) can always be taken, even when both λ 2 and g 2 vanish, because∆ R and 1/|F a (−M )| 2 decay as 1/M 2 and 1/| log(M )| 2 for large M , respectively.
The decay of ψ 1 = D R (f ) for large t descends straightforwardly from Lemma A.1 applied to D R (f ), which implies that ∆ R (f, −s) decays as 1/t 3/2 for large t, with s ≤ 0. Hence, the contribution O(f ) due the poles given in eq. (34) has the desired time decay property. The same holds for the contribution C(f ) due to the cuts given in eq. (35), because for λ 2 = 0 the function M 2 /|S(−M 2 )| is integrable in dM 2 , and, furthermore, f is smooth and compactly supported in time, so its time Fourier transform is a Schwartz function.
Remark 3.6. From the form of the kernel of D R obtained in eq. (31), a generic past compact solution ψ 1 = D R (f ) defined in eq. (32) can be decomposed into two parts, so that denotes the contribution due to the poles of 1/S, while ψ C 1 is the contribution due the cuts. We observe that, while there is a chance to determine ψ O 1 by means of a finite number of initial conditions given at some time t 1 in the future of suppf , we expect that it is not possible to determine ψ C 1 with a finite number of initial conditions, because the integration of M 2 is over uncountably many points.
In spite of this fact, we notice that the homogeneous equation (23) may still, in some cases, give origin to a well-posed initial value problem to uniquely determine spatially compact solutions. Actually, the contribution due the cuts cannot enter the construction of the solutions of the homogeneous equation on the whole space. The reason is that the kernel of the multiplicative operator T , which acts on S(R 4 ) and is defined as , contains only 0, with z = −(p 0 − i ) 2 + |p| 2 . Therefore, only the contributions due the poles can give origin to non trivial solutions of the homogeneous equation (23) written as S(z)ψ 1 = 0.
The decay rate of the smooth past compact solutions ψ 1 proved in Theorem 3.5 is the same which was obtained by means of Strichartz estimates for real, massive quantum scalar fields in four-dimensional Minkowski spacetime [Str77]. Actually, this behaviour is justified by the form of the fourth-order differential equation (24), which is composed of massive Klein-Gordon like operators on (M, η). Thus, the retarded fundamental solution is still a combination of Klein-Gordon like fundamental solutions, and hence the past compact solutions given in eq. (32) inherit the same late-time behaviour estimated for real Klein-Gordon fields.
Eventually, we can now to discuss the solutions of the linearized semiclassical equation without source given in eq. (24).
Theorem 3.7. Consider the semiclassical equation (22) written as Fix as non-vanishing constants at least one of the two g i , and at least one of the λ i , assume that the inequality g 2 λ 1 − λ 2 g 1 ≥ 0 holds, and set −4m 2 < a < 0.
Let S ⊂ (−4m 2 , ∞) ∪ {−λ 1 /λ 2 } ⊂ R be the set of zeros of S(z) given in eq. (27). As discussed in Proposition 3.4, S contains one, two or no elements depending on the parameters λ i , g i and a. If S = ∅, then eq. (22) admits no solutions. If S = ∅, let ψ 1 (t, x) be a smooth solution of eq. (22) with spatial compact support. Then its spatial Fourier transform is of the form Moreover, if S contains only negative elements, then each solution ψ 1 of eq. (22) is uniquely fixed by the initial values at t = 0 ψ (j) where |S| is the cardinality of S, and ϕ j ∈ C ∞ 0 (R 3 ). Furthermore, in this case, ψ 1 (t, x) decays for large time at least as 1/t 3/2 . Proof. We analyze eq. (22) in the Fourier domain. Using the results given in Proposition 3.1, it takes the form: Letψ(t, p) be the spatial Fourier transform of a generic solution ψ 1 (t, x) of the form given in eq. (32). This is a linear combination of e ip j 0 t with z = −p j 0 2 + |p| 2 ∈ S, where S is the set of points of the complex plane in which the function S given in eq. (27) vanishes, namely in which eq. (29) holds. According to Proposition 3.4, we have that S must be contained either in (−4m 2 , ∞) ⊂ C, or in (−4m 2 , ∞) ∪ {−λ 1 /λ 2 } ⊂ C for λ 2 = 0, g 2 λ 1 /λ 2 = −g 1 , and λ 1 /λ 2 ≤ 4m 2 .
According to the number of negative solutions of eq. (22), the explicit form of ψ 1 (t, x) reads as follows. If S contains only one negative solution x = −ñ, withñ ≥ 0, then any solution having smooth compactly supported initial data at t = 0 is of the form where wñ(p) = |p| 2 +ñ, and C ± are obtained from (φ 0 (p),φ 1 (p)) by solving Thus, thanks to Lemma A.1, we obtain the desired decay of ψ(t, x) for large t.
If S contains only two distinct negative elements, four initial data are needed to fix the solution. Denoting with s 1 = −n 1 and s 2 = −n 2 , n i ≥ 0, the two distinct elements of S, the linearized solution of the semiclassical equation (22) is a combination of two solutions of the Klein Gordon equation with different square masses n i . In this case, the solution with smooth compactly supported initial data (ϕ 0 (x), ϕ 1 (x), ϕ 2 (x), ϕ 3 (x)) at t = 0 is of the form where w i (p) = |p| 2 + n i , and C i ± (p) are obtained from (φ 0 (p),φ 1 (p),φ 2 (p),φ 3 (p)) by solving  The determinant of that matrix is equal to −4(n 1 − n 2 ) 2 w 1 w 2 , and hence C i ± (p) can be written as linear combinations of (φ 0 (p),φ 1 (p),φ 2 (p),φ 3 (p)) as Notice that these coefficients have either w 1 or w 2 in the denominator, in the worse case. Thus, the desired decay of ψ(t, x) for large t is obtained by applying Lemma A.1 as before.
(λ 1 + λ 2 x) F a (x) x a −(g 1 + g 2 x)/λ 1  (22) with compactly supported initial values decay at large times. In the following corollary, we identify certain sufficient (but not necessary) conditions on the parameters λ, λ i , g i , a which ensure such a behavior.
If these conditions hold, then the corresponding solutions of eq. (22) with compact spatial support decay at least as 1/t 3/2 at large times.

Proof.
a) Case λ 2 = 0, g2 λ1 ≥ 0. Eq. (30) takes the form The real part of that equation has now a positive solution if − g1 λ1 > λ 16π 2 F a (0). From the plot displayed in Figure 4 we infer that a single negative solution appears if − g1 λ1 ≤ λ 16π 2 F a (0) and λ 16π 2 F a (−4m 2 ) ≤ − g1 λ1 + 4m 2 g2 λ1 , while no solutions exist otherwise. b) Case a < 0, −λ 1 /λ 2 ≤ −g 1 /g 2 ≤ a < 0. From the plot displayed in Figure 3, it is found that all the possible solutions are negative, and in particular either one or two solutions exist. c) Case −λ 1 /λ 2 < a < 0 and 0 ≤ − g1 λ1 ≤ λ 16π 2 F a (0). This case is similar to the one displayed in Figure 3, but now the line −g 2 x + g 1 intercepts the vertical line between y = 0 and y = λ 16π 2 F a (0). Hence, all the possible solutions are negative, and there is always at least one solution.
Finally, the proof follows by applying the results of Theorem 3.7.
To summarize, we have thus proven that the desired decay as 1/t 3/2 for large time t of the solutions of the linearized Einstein equation (22) with compact spatial support holds if and only if the zeros of S defined in eq. (27) are all contained in the negative real axis. On the contrary, if some zeros were located in the positive real axis, then unstable runaway solutions would destabilize the background configuration. Finally, if no zeros are present in S, then eq. (22) does not admit solutions, but its counterpart with source given in eq. (24) still has non vanishing past compact solutions which decay in time, due to the contribution given by the source through the branch cuts.
Remark 3.9. According to the results presented in Proposition 3.4, every solution of eq. (30) is located in the positive real axis whenever the quantum field φ is massless, i.e., when m = 0, even if the inequality (28) holds (recalling that F a (z) is reduced to eq. (21) in the massless case). For this reason, we may expect that any stability result cannot be achieved for massless fields, at least when homogeneous equations are taken into account, even if compactly supported initial data are selected.
This observation is in accordance with the stability issue established in [Hor80], where it was shown that exponentially growing runaway solutions appear when the back-reaction of a quantum Maxwell field interacting with a weak gravitational field is taken into account in the framework of semiclassical gravity.

Applications in the cosmological model
The analysis employed in the toy model presented in eqs. (1a) and (1b) can be used to guess the behaviour of the linearized solutions of the semiclassical Einstein equations in cosmological spacetimes, where matter is modelled by a massive quantum scalar field φ.
In the cosmological scenario we are considering here, the spacetime geometry is described by the metric g µν of the flat Friedmann-Lemaître-Robertson-Walker spacetime. In conformal coordinates (τ, x 1 , x 2 , x 3 ), the metric is conformally flat and reads where the scale factor a(τ ) of the universe represents the unique degree of freedom of the spacetime. The dynamics of this universe is governed by the back-reaction of a linear quantum scalar field φ, whose equation of motion is where ξ denotes the coupling constant to the scalar curvature. Since there is an unique degree of freedom in this class of spacetimes, the dynamics of a is determined, up to a constraint, by the trace of the semiclassical Einstein equations [MPS21] − R + 4Λ = 8πG T ω .
The expectation value of the trace of the quantum stress-energy tensor associated to φ in a quantum state ω is where the renormalization constants α i are not fixed by the model, but describe the regularisation freedom present in the definition of the stress-energy tensor as normal-ordered Wick observable [HW01,HW05,Hac16]. The constants α 1 and α 2 correspond to renormalizations of the cosmological constant and the Newton constant, respectively, and thus they can be reabsorbed in a redefinition of Λ and G; on the contrary, α 3 is of pure quantum nature and cannot be reabsorbed in any corresponding classical parameter of the theory. Moreover, the coefficient [v 1 ] appearing in eq. (40) corresponds to the so-called quantum trace anomaly of the model, and reads up to contributions which can be reabsorbed by a redefinition of α i . Here, C abcd is the Weyl tensor, R ab the Ricci tensor and R the Ricci scalar [Wal78, Mor03,HW05]. We are interested in studying the linearized perturbation of this cosmological system around a spacetime which is a solution of the semiclassical Einstein equations written in the form of eq. (39). As a first step of this analysis, and for the sake of simplicity, we consider as background solution a spacetime with vanishing curvature, namely the Minkowski spacetime. Under this assumption, we can obtain a formal correspondence between the linearization of eq. (39) and the linearized semiclassical equation (22), viewing R as the perturbative external field ψ 1 over a vanishing background ψ 0 = 0. In view of this correspondence, the cosmological constant Λ is a zeroth-order contribution which can be assumed to vanish, whereas the trace anomaly given in eq. (41) is at least quadratic in the components of the Riemann curvature tensor R abcd , and thus it is negligible at linear order in R.
Taking into account all of this and eq. (40), eq. (39) takes the form of the linearized semiclassical equation (22) through the following correspondence between the cosmological parameters and the set of constants (g 1 , g 2 , λ 1 , λ 2 ): In the cosmological framework, g 1 turns to be a fixed negative parameter (in Planck's units, = 1 and (8πG) −1 = m 2 P /8π, where m P is the Planck mass), while, on the other hand, λ can be fixed to be strictly positive and different from 1/6 by assuming non-minimally and non-conformally coupled fields, i.e, ξ = 0, 1/6; hence, λ 2 = 0. On the contrary, both g 2 and λ 2 are free parameters of the semiclassical theory, whose signs can be chosen such that the inequality (28) holds, i.e., Under these assumptions, there are choices of the parameters (m 2 , ξ, a) for which the cosmological version of eq. (30) admits only negative solutions. For example, by choosing ξ > 1/6, α 3 > 0, a > −4m 2 , and sufficiently large m 2 we may apply Corollary 3.8, which ensures that only negative solutions exist. Namely, in these cases solutions of the cosmological linearized semiclassical Einstein equations written as in eq. (39) with spatial compact support decay to zero for large times, thus showing the stability of the chosen background.
On the other hand, one may expect that too large values of m 2 , even beyond the Planck scale m P , would be physically unacceptable for quantum fields describing elementary particles. In this viewpoint, the result is similar to the one obtained in [RD81] for massive quantum scalar fields in flat spacetime. Firstly, the conditions stated here are only sufficient, so other cases which provide stable solutions cannot be excluded a priori, for different choices of the parameters (m 2 , ξ, γ, a). Secondly, and most importantly, it is expected that the linearized perturbations in a more realistic cosmological model should be sourced by f ∈ C ∞ 0 (M) localized somewhere in the past. This source may have a quantum origin, for instance, related to some anisotropic or stochastic fluctuations at microscopic levels. This is the case which occurs for example in Stochastic Gravity, where a noise kernel bi-tensor modelling the stress-energy tensor fluctuations is added to the semiclassical Einstein equations, obtaining in this way the so-called Einstein-Langevin equations (see [HV20] and reference therein). In this picture, the stochastic source in the past drives the fluctuations of the gravitational field, and thus gives origin to the external perturbation which enters the cosmological linearized semiclassical Einstein equations as external source.
In this model with external sources, the cosmological counterpart of the linearized semiclassical equation (24) should be taken into account, with parameters λ 1 , λ 2 , g 1 , g 2 fixed as in eq. (42) and satisfying the inequality (43). Under these assumptions, and based on the results shown in Theorem 3.5 and Remark 3.6, the linearized curvature solution R depends on both the contributions due to the poles and the branch cuts of S(z). However, the contribution arising from poles are not present for several, apparently more physically acceptable values of the parameters (m 2 , ξ, a): for example, for sufficiently large ratio m 2 P /m 2 and 0 < ξ < 1/6, α 3 > 0, a > −4m 2 , the condition (43) holds. Furthermore, with this choice of parameters, the past compact linearized solution induced by a smooth compactly supported source f has no poles contribution, and hence it decays to zero for large times according to the results stated in Theorem 3.5.

Conclusions
In this paper, we have analyzed the stability problem of semiclassical theories in flat spacetime, using a toy model consisting of a quantum scalar field coupled to a second entirely classical scalar field. This toy model mimics other semiclassical theories of gravity described by the semiclassical Einstein equations, where a quantum matter field propagates over a classical curved background. It is known that higher order derivatives appearing in the semiclassical equations can destabilize the system, giving rise exponentially growing linearized solutions (see the references given in the Introduction). The main result stated in this paper consists of proving that, if the quantum field driving the back-reaction is massive, then the stability of background solutions can be restored at the linear order in the interaction, for spatially compact perturbations and for large values of the coupling constants, after assuming some sufficient (but not necessary) conditions on the parameters of the theory. On the other hand, removing the assumption of massive quantum fields seems to give rise to runaways solutions which may alter stability, in accordance with other results present in the literature about semiclassical theories. Namely, it is shown that unique solutions of this semiclassical initial value problem tend to disperse in time, namely at fixed position in space they decay polynomially in time.
As this toy model mimics the dynamics of the linearized semiclassical Einstein equations, our analysis indicates a possible mechanism to get stability in several linearized semiclassical theories of gravity, even for different conditions than the ones stated in this paper. For example, in the case of a semiclassical theory in cosmological spacetimes, which has been already investigated by the authors in [Pin11,PS15,MPS21].

A Large time decays of certain functions
This appendix contains a Lemma with the proof of the decay at large times of certain functions. The main idea of this proof is already known in the literature, see, e.g., [BB02], however since this Lemma is a key result for the stability discussed in the main text, we report its proof here for completeness.
Lemma A.1. Letf ∈ S(R 3 ) be a Schwartz function, and consider its spherical average f (|p|) = 1 4π f (p)dΩ, where Ω denotes the standard measure on the two-dimensional surface of the unit sphere. Consider e iw m 2 t w m 2 p 2 dp, w m 2 = p 2 + m 2 , then W (t) decays for large times as 1/t 3/2 if m > 0, and as 1/t 2 if m = 0.
Proof. Let us start discussing the massless case m = 0, In this case w 0 = p then p e ipt dp.
After integrating by parts twice, and using the rapid decay of f ∈ S(R 3 ), where ∂ 2 p (pf (p)) ∈ L 1 (R + ) because f (p), together with its derivatives, is of rapid decrease. Hence, by Riemann-Lebesgue lemma lim t→∞ ∞ 0 ∂ 2 p (pf (p))e ipt dp = 0, which implies that the second contribution in eq. (44) vanishes more rapidly than 1/t 2 for large times.
We pass now to analyze the case m > 0. After changing variable of integration p = (y/t + 2m)(y/t), W (t) = e imt t 3/2 ∞ 0 f y t y t + 2m e iy y t + 2m √ ydy.
To evaluate this integral, we insert an regulator, and we divide the integral in two parts, where g(y) = √ y + 2mf ( y(y + 2m)). Notice that the integral in the first contribution tends to a constant in the limit → 0, because lim →0 + ∞ 0 e iy− y √ ydy = √ π 2(−i) 3/2 .
After integrating by parts twice, and in view of the decay properties of g and its derivatives for large arguments, ) − g(0) e iy− y − 1 dy.
Hence, for 0 < δ < 1/2 we obtain that is a bounded function, because (g(x) − g(0)) vanishes as x for x near 0, and g decays faster to 0 for large x. Therefore, ∞ 0 e iy− y − 1 y 3/2−δ dy for a suitable constant C. Since C is uniform in , and 0 < δ < 1/2, we can apply dominated convergence theorem to take the limit as → 0 before computing the integral, and hence |A| ≤C t δ , which concludes the proof.