The Λ cold dark matter ( Λ CDM) concordance model (Cervantes-Cota and Smoot, 2011; Deruelle and Uzan, 2018; Dodelson and Schmidt, 2021), based on the pillars of the Standard Model of particle physics and general relativity, has shaped our current view of the Universe, and has been the driving force behind many of the breakthroughs of modern cosmology, for example the prediction and the discovery of the cosmic microwave background radiation (Aghanim et al., 2019, 2020; Alpher and Herman, 1948a,b; Gamov 1948a,b; Penzias and Wilson 1965). Modeled by only six parameters (Spergel 2015; Aghanim et al., 2020), it features an impressive simplicity while correctly predicting and fitting most of the cosmological data (Aghanim et al., 2019, 2020). One of the most important assumptions within the Λ CDM paradigm is that the Universe is almost spatially homogeneous and isotropic, especially during its earliest phases, but even today when considered on its largest scales. The resulting simplification of Einstein’s equations is remarkable as it reduces the ten coupled non-linear partial differential equations in four variables to two ordinary equations in one variable, with solutions known as the Friedmann–Lemaître–Robertson–Walker (FLRW) solutions (Friedman 1922, 1924; Lemaître 1931; Robertson 1933; Walker 1937).

Obviously, a look at the night sky reveals that the Universe is not homogeneous and isotropic, but is characterized by clusters of galaxies and stars, and large voids inbetween (Blumenthal et al., 1984; Cole et al., 2005; Colless et al., 2001; Ross et al., 2020; Zel'dovich et al., 1982). For explanation, the concordance model assumes that smallest quantum fluctuations of the primordial matter and geometry have been stretched to the present time, thereby generating the observable large scale structure. Importantly, these inhomogeneities on any scale smaller than the observable Universe are presumed to evolve following the underlying FLRW background structure, but conversely their evolution does not affect the global FLRW evolution. More precisely, it is assumed that effects from the small scale inhomogeneities onto the largest scales can be neglegted, i.e., there is no substantial backreaction.

Doubts regarding the simplistic nature and the question of backreaction have gained momentum in recent years. In fact, the Λ CDM model, as appealing it may be, leads to the conclusion that approximately 69% of the energy budget of our Universe consists of a yet unknown fluid, dubbed “dark energy,” (Aghanim et al., 2020), and which drives the very recent accelerated expansion of the Universe (Riess et al., 1998; Perlmutter et al., 1999; Peebles and Ratra 2003). Most of the remaining 31% of the energy budget is credited to another yet unknown form of cold “dark” matter (Peebles 1982; Blumenthal et al., 1984; Aghanim et al., 2020), which provides an explanation for the characteristic rotation and motion of the remaining 6% of ordinary matter in the Universe. In summary, we are faced with the problem that we are literally in the dark about 94% of the energy and matter content of the observable Universe.

In recent years, these conceptual problems have been accompanied by important tensions in the estimates of certain cosmological parameters as made by different collaborations (Di Valentino et al., 2020a,b; Pesce et al., 2020). The evaluation of the Hubble constant H 0 as performed by the Planck collaboration (explicitely assuming a Λ CDM model) gives a value of H 0 = ( 67.27 ± 0.60 ) km / ( s ⋅ Mpc ) (Aghanim et al., 2020), while the SH0ES collaboration finds H 0 = ( 74.03 ± 1.42 ) km / ( s ⋅ Mpc ) (Riess et al., 2019), which in turn is based on the measurements of the Hubble Space Telescope. This leads to a tension at the 4.4 σ level (Di Valentino et al., 2020a). Furthermore, we point to the (albeit weaker) tensions regarding the measurement of the parameter S 8 , a measure for the matter energy density Ω m and the amplitude of structure growth σ 8 (Aghanim et al., 2020; Di Valentino et al., 2020b).

On the other hand, the theoretical modeling of the early and very early Universe turns out to be a difficult undertaking, in particular the faithful consideration of all interactions within coupled quantum cosmological–matter systems. Since classical cosmological perturbation theory and its various applications to the physics of our Universe (Durrer, 2004; Mukhanov, 2005), represents a successful formalism to model (most of) the cosmological data today, one of the most promising approaches to make progress in the field is to consider an inhomogeneous, but perturbative, quantum cosmology, i.e., to establish a quantization of the well-known (possibly) gauge-invariant cosmological perturbation theories (Brandenberger et al., 1993; Brandenberger, 2004; Elizaga Navascués et al., 2016). In fact, there has been tremendous progress in developing such quantum cosmological perturbation theories, for example, in quantum geometrodynamics (Kiefer, 2007; Brizuela and Krämer, 2018), in string cosmology (Erdmenger, 2009), as well as in loop quantum cosmology (LQC) and spinfoam cosmology (Bianchi et al., 2010; Vidotto 2011; Cailleteau et al., 2012; Agullo et al., 2013; Elizaga Navascués et al., 2016), to mention but a few. Unfortunately, the majority of these approaches neglect backreaction effects from the inhomogeneous quantum fields on the homogeneous, dynamical degrees of freedom, or incorporate a series of assumptions which are hard to control, similar to the situation in classical cosmological perturbation theory.

It seems hence very timely to scrutinize and question the various assumptions of the concordance model of cosmology, and to develop suitable formalisms which are able to take interactions in coupled (quantum) cosmological models more realistically and unambigiously into account. In this review, we start by assessing the question of backreaction, i.e., whether cosmological inhomogeneities have an effect on the large scale evolution of the Universe, especially in view of the occurent inconsistencies within the standard model. We consider different aspects of backreaction, in particular we discuss backreaction in classical, semiclassical and quantum mechanical models. Our main focus is on the purely quantum mechanical backreaction and we discuss one recent approach to including backreaction in quantum cosmology in more detail (Schander and Thiemann, 2019a). The structure of the paper is then as follows.

In Section 2, we provide an overview of the results in the field of classical backreaction, which is particularly relevant for late time cosmological models. In Section 3, we consider semiclassical backreaction which occurs when considering quantum fields on classical curved space times. Section 4 gives an overview of quantum backreaction, i.e., backreaction that occurs in purely quantum theoretical models. In Section 5, we focus on one particular approach to quantum backreaction which uses mathematical tools inspired by the Born–Oppenheimer approximation. Section 6 provides a final discussion and an outlook.

Standard perturbative approaches to cosmological perturbation theory implicitely conjecture that backreaction, i.e., the effects of cosmological inhomogeneities on the global or macroscopic evolution of the Universe can be ignored. For purely classical models of the Universe that are particularly relevant for its late time evolution, this conjecture has generated an intense debate over the last decades. And still, there is no consensus on the question of backreaction in the classical regime , see for example the reviews by Clarkson et al. (2011), Ellis (2011), and Bolejko and Korzynski (2017).

The question of backreaction is closely related to the fitting problem, (Ellis and Stoeger, 1987), and the problem of averaging, (Clarkson et al., 2011). In fact, an intuitive way to access the effects of small inhomogeneities on the macroscopic scales is to construe an averaging procedure that defines new homogeneous variables by integrating the inhomogeneous fields over a certain space time domain, and to compare their properties and dynamics to the assumed FLRW Universe, (Ellis, 2011). However, it is inadmissible to conclude from the validity of Einstein’s equations for the inhomogeneous fields on the smallest scales (where they have been excellently checked), that the averaged fields satisfy the Einstein equations (Paranjape, 2012). This is because evaluating the Einstein tensor and taking a space (time) average does not commute in general. Hence, the averaging procedure can lead to additional contributions to Einstein’s equations that might be considered as effective source terms for the geometry, see for example (Buchert, 2000, 2001).

As it turns out, the results regarding the form and strength of backreaction depend heavily on the averaging procedure and the matter model being chosen, as well as on the choice of space time volumes to be integrated over. In the non-perturbative regime, the two most discussed averaging procedures are the scalar averaging scheme by Buchert (2008) and the Macroscopic Gravity approach by Zalaletdinov (1997, 2008). Buchert’s scheme focuses on building spatial averages of scalar fields and derives effective (scalar) equations of motion for the averaged quantities, for example an improved Raychaudhuri equation for the averaged scale factor that includes a kinematical backreaction term. While being technically easy to implement, the Buchert scheme relies on a system of scalar equations that is not closed (Clarkson et al., 2011), and consequently requires additional information to fix the solution. Besides, the averaging demands to fix suitable spatial domains and hence, a hypersurface slicing. In contrast, the Macroscopic Gravity approach is manifestly covariant but requires to define an auxiliary so–called bi–local transport operator, (Ellis, 2011). Physical applications of these schemes yield a range of different results, ranging from explaining the recent accelerated expansion of the Universe or the H 0 -tension (Buchert and Räsänen, 2012; Heinesen and Buchert, 2020), to negligible backreaction effects (Paranjape and Singh, 2007, 2008).

Many of the afore-mentioned approaches (among many others) assume the matter content to be modeled by a fluid, which is likely to be a poor approximation to the true lumpy late time Universe. Models with more realistic matter distributions are for example the Timescape Cosmology by Wiltshire (2009), who separates the Universe into underdense expanding regions bounded by overdense virialized structures, the Swiss Cheese Model (Kantowski, 1969; Tomita, 2000; Biswas and Notari, 2008), or modifications of FLRW Universes that cut spherically symmetric Lemaître–Tolman–Bondi or Szekeres dust space time regions (Marra et al., 2008; Bolejko and Celerier, 2010), to mention but a few. By construction, many of these models follow the evolution of an appropriately fitted FLRW model since they assume a background structure from the beginning. Consequently, they do not attack the backreaction problem outlined before. In contrast, the model by Lindquist and Wheeler (1957) assembles static Schwarzschild regions without relying on any background, and which has been further investigated by Clifton (2011) and Clifton and Ferreira (2009). In both cases, the models provide insights into backreaction effects on light propagation (Krasinski and Bolejko, 2011; Sussman, 2011), which points to another important topic.

In fact, cosmological observations such as the distance–redshift relation or the angular diameter distance rely on measurements of light, traveling along our past lightcone in a very inhomogeneous Universe. The seminal work by Kristian and Sachs (1966) laid out the basis for analyzing backreaction on light propagation. Flanagan (2005), for example, used these ideas to compute the deceleration parameter as measured by comoving observers. Gasperini et al. (2011) define a covariant light–cone average for the backreaction problem, see also (Fanizza et al., 2020) for a more recent generalized proposal. Räsänen (2009) and Räsänen (2010) derives a relationship of the redshift and the angular diameter distance to the average expansion rate for statistically homogeneous and isotropic universes, based on Buchert's approach, and Barausse et al. (2005) and Bonvin et al. (2006) evaluate the distance–redshift relation and the luminosity distance in a perturbative framework. Most recently, Heinesen (2021a), Heinesen (2021b) and Koksbang (2019), Koksbang (2020), Koksbang (2021) investigated the effects of inhomogeneities and averaging on a possible redshift drift.

Many approaches that attempt to make direct contact with cosmological observations restrict their analysis to cosmological perturbation theory in an FLRW Universe, as opposed to the above–mentioned non-perturbative approaches. Most of them consider flat Λ CDM models with Gaussian scalar perturbations as initital conditions. To evaluate backreaction, they compute the deviations to the Hubble expansion rate or similar variables that are caused by backreaction (Brandenberger et al., 2018; Clarkson et al., 2009; Kolb et al., 2010; Kolb et al., 2005; Li and Schwarz, 2008), or give effective Friedmann equations with additional contributions (Paranjape and Singh, 2007; Behrend et al., 2008; Brown et al., 2009; Peebles, 2010; Baumann et al., 2012). The idea is to perform appropriate spatial averages of the perturbed quantities and to use the given statistical information of the perturbation fields in guise of their power spectra. It turns out that due to the smallness of the gravitational potential and the power suppression of modes on large scales, backreaction for the expansion rate is always small. However, the backreaction to the deceleration parameter q and the variance of the Hubble rate depend on an auxiliary UV–cutoff that might lead to large backreaction even if it is set by scales larger than the non-linearity scale (Clarkson et al., 2011). Other approaches to backreaction in the linear regime are (Baumann et al., 2012) who propose a reformulation of perturbation theory that leads to small backreaction on the largest scales but affects the baryon accoustic oscillations, and (Green and Wald, 2011, 2012, 2013, 2014) who claim, using a point limit process, that backreaction can never mimic dark energy and put strong constraints on its strength.

We also point to the quite recent advent of numerical tools that allow to simulate increasingly realistic models of the Universe, including relativistic effects (Löffler et al., 2012; Mertens et al., 2016) and N-body simulations (Adamek et al., 2016; Barrera-Hinojosa and Li, 2020). Using the N-body relativistic code “gevolution,” Adamek et al. (2019) find that backreaction on the expansion rate in a Λ CDM and an Einstein–de Sitter Universe remains small if one chooses averaging volumes related to the Poisson gauge, while when choosing comoving gauge backreaction is of the order of 15 % . Other works in this respect were done by Macpherson et al. (2019), who also claim that backreaction effects are small, however based on a fluid approximation which breaks down as soon as it comes to shell crossing.

Finally, let us also point to the consideration of backreaction from long wavelength modes of the early Universe. In this respect, early contributions were notably made by Tsamis and Woodard (1993, 1996), as well as by Abramo et al. (1997) and Mukhanov et al. (1997). The latter works pursue a gauge—invariant formulation of the backreaction problem associated with an effective long wavelength energy momentum tensor, and within a slow–roll inflationary scenario. Unruh (1998) subsequently examined the question of whether this effect is indeed locally measurable, and it was found that such backreaction effects (in single field inflationary theories) can be absorbed by a gauge transformation, (Abramo and Woodard, 2002; Geshnizjani and Brandenberger, 2002). However, backreaction of such fluctuations becomes locally measurable after introducing an additional (subdominant) clock feld, (Geshnizjani and Brandenberger, 2005). This approach was extended by Marozzi et al. (2013) based on the gauge—invariant formalism by Finelli et al. (2011), and secondly by Brandenberger et al. (2018) beyond perturbation theory. Further contributions were made by Losic and Unruh (2005), Losic and Unruh (2008) who support the idea that backreaction represents a real and measurable effect.

For considerations of backreaction in models of the very early Universe, the standard model of cosmology suggests that (at least) the matter fields should be studied in a quantum mechanical framework. The implementation of such ideas can be realized via different paths, and we consider here the approaches of semiclassical gravity (Ford 2005; Wald 1977, 1978) and stochastic (semiclassical) gravity (Calzetta and Hu 1987; Hu and Verdaguer 2008; Jordan 1986, 1987). Both approaches rely on the framework of quantum field theory on curved space times (QFT on CST) (Birrell and Davies, 1984; Fulling 1989), which itself takes the effects of the classical curved space times on the quantum matter fields into account but in general not the backreaction effects of those quantum fields on the classical background.

The backreaction problem in semiclassical gravity was first brought in by Wald (1977) who considered the backreaction from particle creaction on a gravitational field. The idea of the semiclassical program is to consistently define an improved set of Einstein field equations in which the expectation value of the quantum stress–energy tensor T μ ν of the matter fields with respect to an appropriate quantum state of the matter fields ω appears as a source term, R μ ν + 1 2 g μ ν R = 8 π G ω ( : T μ ν : ) , (1) where the quantization with respect to the matter fields is expressed using bold letters and the dots indicate the normal ordering of the stress–energy tensor. The state ω should be considered as a positive linear functional in the sense of the algebraic approach to quantum field theory (QFT) (Haag 1992; Araki 1999).

The first goal of semiclassical gravity is to define a procedure that leads to a meaningful expression for the expectation value of the stress–energy tensor. In fact, the latter depends on products of operator-valued distributions, even for the simple case of a real-valued Klein–Gordon field, and its expectation value is in general a divergent expression. Wald (1977) gave a set of axioms that are required to hold for a suitable renormalization scheme. Possible proposals are the Hadamard point-splitting method (Brunetti and Fredenhagen, 2000; Hollands and Wald, 2001), and the adiabatic regularization procedure (Fulling and Parker 1974; Fulling et al., 1974; Parker and Fulling, 1974). In either scheme, the result of the regularization procedure is a set of modified “semiclassical” Einstein equations. These equations are substantially harder to solve than the original Einstein equations and many studies restrict to cases of conformally coupled matter to avoid problems regarding the well-posedness and the stability of the solutions (Ford, 2005). Caution is also required regarding the question of self-consistency of the backreaction effects, as has been discussed by Flanagan and Wald (1996).

Many applications of the semiclassical gravity approach to early Universe cosmology have been considered. For example, Fischetti et al. (1979) analyzed the backreaction effects from a conformally invariant matter field in an FLRW Universe with classical radiation, and found that the trace anomaly can soften the cosmological singularity, but not avoid it. Other works in this direction were done by Anderson (1983, 1984, 1985), who also considered the trace effects on the particle horizon. A well-known example of trace anomaly effects from semiclassical gravity is the Starobinsky (1987) cosmological model, which has however not survived the observational scrutiny of the Planck data (Ade et al., 2016).

Another application of semiclassical gravity is the study of backreaction of particle creation on the dynamics of the early Universe, as already conceived by Wald (1977). Grishchuk (1977) as well as Hu and Parker (1977) considered the effect of gravitons around the Planck time in an FLRW Universe with a classical, isotropic fluid. The model leads to a timely non-local (i.e., history-dependent) backreaction effect, (Hu and Verdaguer, 2020). Similar studies were performed for anisotropic FLRW Universes and it was shown that particle production due to the shear anisotropy isotropizes space time (Hu and Parker, 1978; Zel'dovich and Starobinsky, 1972). Regarding the effects of particle creation in a spatially inhomogeneous but isotropic Universe, we refer to the work by Campos and Verdaguer (1994).

We also point to the more recent works by Finelli et al. (2002), Finelli et al. (2004) who specifically consider a slow—roll (almost de Sitter) phase of the very early Universe and compute a(n adiabatically) renormalized energy momentum tensor of quantum inflaton, respectively cosmological scalar fluctuations. In case of the cosmological scalar perturbations, they found that the energy momentum tensor is characterized by a negative energy density which grows during inflation, and also that backreaction is not a mere gauge artifact.

Further contributions to the topic of semiclassical gravity for cosmological situations were notably made by Dappiaggi et al. (2008); Dappiaggi et al. (2010); Eltzner and Gottschalk (2011); Gottschalk and Siemssen (2018); Hack (2013); Matsui and Watamura (2020); Parker and Raval (1999); Pinamonti (2011), to mention but a few. Most recently, Meda et al. (2020) and Pinamonti and Siemssen (2015) have made progress on the definition of the semiclassical theory for general couplings by proving existence and uniqueness of solutions in flat cosmological space times with a massive quantum scalar field. The idea of relating backreaction effects to the decay of a cosmological constant has for example been promoted by Dymnikova and Khlopov (2001). We also point to the work by Matsui and Watamura (2020) (and references therein) who claim that the approach of semiclassical gravity is not appropriate to describe the early Universe.

The second approach to evaluating backreaction in semiclassical cosmology that we present here, denoted stochastic gravity (Hu and Verdaguer 2020), creates a link to open system concepts and statistical features such as dissipation, fluctuations, noise and decoherence. It employs a so-called closed time path coarse grained effective action (CTP CGEA) (Calzetta and Hu 1987; Jordan 1986, 1987), in order to derive a set of modified semiclassical Einstein equations, denoted as Einstein–Langevin equations. It includes the semiclassical approach but extends it by a stochastic noise term (Hu and Matacz, 1995).

Some of the first applications of the CTP CGEA formalism to the backreaction problem in cosmology were made by Calzetta and Hu (1987, 1989, 1994). Hu and Matacz (1995) derived the Einstein–Langevin equations for the case of a free massive scalar field in a flat FLRW background, as well as for a Bianchi Type-I Universe. The case of a massless conformally coupled field was discussed in (Campos and Verdaguer, 1994). The scope of works includes topics such as stochastic inflation, where quantum fluctuations present in the noise term backreact on the inflaton field (Calzetta and Hu, 1995; Lombardo and Mazzitelli, 1996), as well as studies of the reheating phase in inflationary cosmology (Boyanovsky et al., 1995; Ramsey and Hu, 1997). The formalism was also used by Sinha and Hu (1991) to check the validity of the minisuperspace approximation in quantum cosmology. Further applications can be found in the paper and textbook by Hu and Verdaguer (2008), Hu and Verdaguer (2020).

We also point to one of the most prominent applications of stochastic methods to early Universe cosmology by Starobinsky (1982) and Starobinsky (1988). His stochastic inflationary model evaluates backreaction of small scalar field quantum perturbations on the corresponding long wavelength modes (which are assumed to behave classically) by additional stochastic terms in the long wavelength equations of motion. A slow-roll behavior of the background is assumed. Interestingly, it has been shown that the stochastic and the quantum feld theoretic approaches to perturbations in the early Universe yield the same results, (Starobinsky and Yokoyama, 1994; Finelli et al., 2009; Tsamis and Woodard, 2005). For recent considerations of stochastic inflation beyond the (strict) slow-roll conditions, we refer to the work by Pattison et al. (2019) and references therein.

Both approaches, semiclassical as well as stochastic gravity regard the gravitational field as a classical entity from the start while the matter fields are considered to be of quantum nature. While this represents a seminal progress to incorporating quantum effects of the matter fields in the early Universe, it can and should be questioned whether this somehow incompatible approach (classical and quantum fields treated at the same level) survives the test of future observations, and whether it should be replaced by a more consistent approach - quantum gravity - at least for the earliest moments of the cosmic history.

The question of backreaction in quantum gravity and quantum cosmology encompasses a variety of different approaches and definitions of backreaction. In quantum cosmology, backreaction is usually identified as the effects from the inhomogeneous quantum perturbation fields on the (quantum) homogeneous and isotropic degrees of freedom, which is also the notion of backreaction used in the next section (Schander and Thiemann, 2019c). This approach is tightly related to a perturbative expansion with respect to the inverse Planck mass m Pl − 1 = ( G / ℏ ) 1 / 2 , where G is Netwon’s constant and we set ℏ ≡ 1 . More precisely, it employs a Born–Oppenheimer type scheme (Born and Oppenheimer, 1927), with respect to m Pl − 1 . We will thus focus on implementations of the Born–Oppenheimer method to quantum gravity and quantum cosmology.

In fact, the idea that quantum gravity can be considered as a perturbative theory with respect to m Pl was already introduced by Brout (1987). The first investigations of backreaction in quantum gravity that rely on this expansion were performed in the framework of quantum geometrodynamics (Wheeler, 1957; Kiefer, 2007). The idea is to expand the Wheeler–DeWitt equation in terms of the ratio of the Planck mass and the matter field mass (Kiefer and Singh, 1991). A different idea, conceptually similar to the schemes considered here, is to use a Born–Oppenheimer type approach, relying on the same perturbation parameter. Different considerations of the problem (giving rise to similar results) can be found in the works by Bertoni et al. (1996), Brout and Venturi (1989), and Kiefer (1994) (for a summary, see Kiefer, 2007). A review of the ideas of Born and Oppenheimer will be given in the next section, but to understand its use in the given context we present the key ideas.

For simplicity, let Q denote the gravitational and q the matter degrees of freedom. The Born–Oppenheimer scheme employs an ansatz solution for the quantum Hamiltonian and momentum constraint of the form (Kiefer, 2007), Ψ ( q , Q ) = ∑ n χ n ( Q ) ψ n ( q , Q ) , (2) where { ψ n ( q , Q ) } n is supposed to be a known orthonormal basis of the matter Hilbert space that solves the matter part of the constraint and Q is to be considered as an external parameter for this eigenvalue problem. Then, one applies the constraints to Ψ and applies some ψ k ( q , Q ) from the left (i.e., one considers the inner product of the matter states). This gives rise to constraint equations for the geometric factors χ n ( Q ) , which can be seen as an effective quantum problem for the geometric part, including the backreaction effects of the quantum matter system. In this scenario, the Born–Oppenheimer approximation consists in neglecting the contributions that enter with higher orders in m Pl − 1 .

In order to extract physical results from the formalism, one can additionally employ a semiclassical approximation (which is however independent of the Born–Oppenheimer approach). This should yield the semiclassical limit of quantum gravity, i.e., a matter QFT on CST. It is common to employ a WKB ansatz for the geometrical states χ n ( Q ) of the form, χ n ( Q ) = C n ( Q ) e i m Pl 2 S [ Q ] , (3) where S [ Q ] stands for the geometric action in the geometrodynamical approach. The perturbative scheme in m Pl − 1 eventually yields the semiclassical Einstein equations. In this sense, these approaches evaluate the backreaction of the quantum matter fields on the quantum or classical geometry.

One can apply the Born–Oppenheimer and WKB approximations in a different manner. Instead of taking the expectation value with respect to the quantum matter system, one applies the Wheeler–DeWitt constraints on the total Born–Oppenheimer ansatz function and uses the WKB approximation for the geometrical part. Restricting again to the lowest order with respect to the inverse Planck mass, this yields a quantum constraint for the matter wave function which depends on the classical action (through the WKB ansatz), and derivatives with respect to the spatial metric thereof. The idea of the above–cited works (and also of Briggs and Rost, 2000) is to introduce an external time parameter that depends on this derivative, hence giving rise to a Schrödinger equation for the matter system that includes the backreaction of the geometry through the geometry–dependent time derivative. In fact, this gives rise to a notion of time in a formerly background independent framework. Such ideas go back to DeWitt (1967) and have been applied to a variety of cosmological situations (see (Kiefer 2007) and references therein). It is however a different notion of backreaction than the one considered in the next section. Besides, the present approach uses a Born–Oppenheimer approach plus a semiclassical WKB approximation while the next section uses the purely quantum mechanical space adiabatic perturbation extension of the Born–Oppenheimer scheme. Applications of the former works to the inflationary paradigm with perturbations and a discussion of the question of unitary evolution of the perturbations can be found in the work by Chataignier and Krämer (2021) and references therein. They also consider cosmological perturbations that include gravitational contributions (i.e., the Mukhanov–Sasaki variables). Similar approaches that do not split the system into geometric and matter parts but include (perturbative) parts of the gravitational degrees of freedom in the fast subsystem and (homogeneous) matter parts in the slow sector were already presented by Halliwell and Hawking (1985) and Vilenkin (1989). This choice is also used in Schander and Thiemann (2019c).

The Born–Oppenheimer approximation was also considered within approaches to quantum gravity with other variable choices. Giesel et al. (2009) aimed at an application of the Born–Oppenheimer methods to loop quantum gravity (LQG) (Thiemann, 2008; Rovelli, 2010), using holonomy–flux variables or connection–flux variables. As it turns out, this choice of variables prevents the use of the Born–Oppenheimer methods since the flux operators are mutually non-commuting (which is a prerequisite for the Born–Oppenheimer scheme). Instead, they use commuting co-triad variables for the gravity sector and a scalar field for the matter sector to derive the semiclassical Einstein equations that take the backreaction of the quantum matter fields via an expectation value into account. Giesel et al. (2009) consider their model on a discrete lattice (as it is common practice for approaches to LQG), and thus formally obtain a lattice QFT on a discrete curved space time. They also point to the possibility of pursuing the formal Born–Oppenheimer scheme and computing quantum solutions to the gravity sector with the effective backreaction of the quantum matter fields. Besides, they introduce a hybrid approach (similar to the models we consider here) where the gravitational sector is restricted to FLRW solutions and the fast part of the system is given by the matter quantum fields. They also propose to introduce coherent states for the gravitational subsystem in order to make progress in finding solutions. Due to the complexity of the gravity-matter systems, the focus of this work lies on spelling out the conceptual ideas rather than technically carrying out the program in detail.

More recently, Stottmeister and Thiemann (2016a,b,c) considered similar questions in the context of LQG but employed the more general scheme of space adiabatic perturbation theory (SAPT) (Panati et al., 2003). Since in the latter approach, the variables of the slow, gravitational sector are not required to commute, it is in principle possible to apply the Born–Oppenheimer ideas also to LQG and related theories. The concrete implementation turns however out to be difficult due to the particular structure of the LQG phase space (which relies on a cotangent bundle of a compact Lie group rather than on a vector space) and its quantum represention. Other open issues of their attempts are related to the underlying graph structure of LQG models and the projective limits of finite dimensional truncations of the gravitational phase space that are needed in order to construct a continuum theory (Stottmeister and Thiemann, 2016c). They also point out that a major obstruction to the derivation of a QFT on CST from LQG lies in the inequivalent representations of quantum fields for different gravitational configurations (Stottmeister and Thiemann, 2016c). This problem is a generic feature of background dependent quantum field theories. In this work, we present a (perturbative) solution to this problem which makes the application of space adiabatic methods to quantum cosmology possible (Fernandez-Mendez et al., 2012; Castelló Gomar et al., 2015; Castelló Gomar et al., 2016).

For completeness, we also mention the application of Born–Oppenheimer methods within the spinfoam approach to LQC (Rovelli and Vidotto, 2008) (see Ashtekar et al., 2003; Bojowald, 2008 for the LQC approach), and that Castelló Gomar et al. (2016) consider a conceptually different kind of Born–Oppenheimer approximation in the hybrid approach to LQC. These seminal works make important progress by first considering the problem of backreaction, but must either remain on a rather formal level or include various assumptions which are hard to control. We therefore advocate to employ the SAPT scheme presented in the next section which serves as an unambigious, rigorous and perturbative approach, in principle applicable to any quantum cosmological framework and realizable up to any perturbative order, to taking quantum cosmological backreactions thoroughly into account. The approach applies to a much wider variety of quantum systems in comparison to the Born–Oppenheimer approach (in particular to quantum cosmological perturbation theory), does not rely on the introduction of semiclassical ansatz states and iteratively provides quantum constraints or equations of motion whose solutions approximate the true solutions up to, in principle, indefinitely small errors.

Computing and including backreaction in (perturbative) quantum cosmology requires an approximation scheme that ideally takes the physical characteristics of the system into account. SAPT as proposed by Panati et al. (2003) and extensions thereof are ideally suited to achieve this goal and to integrate backreaction effects into quantum cosmology (Schander and Thiemann, 2019a; Schander and Thiemann 2019b; Schander and Thiemann 2019c; Neuser et al., 2019).

SAPT is a generalization of the well-known Born–Oppenheimer approximation for non-relativistic molecular systems. Both approaches exploit the small ratio of two internal parameters such as the mass ratio of electrons and nuclei in a molecule to define a perturbation parameter, ε 2 : = m e m n ≈ 5.46 × 10 − 4 ≪ 1 , (4) with the electron mass m e ≈ 9.11 × 10 − 31 kg and the nuclei mass m n ≈ 1.67 × 10 − 27 kg. In the simplest atom with one nucleus and one electron (although exact solutions are known for this case), the Hamilton function has the form, H ( q , P ; x , y ) = ε 2 P 2 2 m e + y 2 2 m e + V ( q ; x ) , (5) where ( q , P ) and ( x , y ) are the canonically conjugate pairs of the nucleus and the electron respectively. V ( q ; x ) is a smooth potential, typically a Coulomb potential depending on the distance between nucleus and electron. In this molecular set up, the equipartition theorem states that the kinetic energies of nuclei and electrons are of the same order, and hence, on average, the nuclei move much slower than the electrons with correspondent statistically–averaged velocities, 〈 v n 〉 ≈ ε 〈 v e 〉 . Born and Oppenheimer used this fact to define suitable ansatz solutions for the quantum mechanical problem: On the typical electronic time scale, the nuclei are at rest and the non-trivial electronic contributions of the Hamilton operator can be considered at fixed nuclei positions q ∈ ℝ , H e ( q ) : = y 2 2 m e + V ( q , x ) , (6) where the bold letters y , x denote momentum and position operators of the electron defined on their respectively dense domains in L 2 ( ℝ , d x ) . Ideally, the operator function H e ( q ) admits a solvable q–dependent eigenvalue problem, H e ( q )   ξ n ( q ) = E n ( q )   ξ n ( q ) , (7) with a discrete q–dependent eigenbasis { ξ n ( q ) } n ∈ ℕ in the fast Hilbert space ℋ f : = L 2 ( ℝ , d x ) , for which the so–called electronic energy bands E n ( q ) are gapped functions, i.e., E n ( q ) − E m ( q ) ≠ 0 pointwise for m ≠ n . One can use this eigenbasis as an ansatz solution for the full Hamilton operator H ^ (the Weyl quantization with respect to the nuclei sector is labeled by hats), and ask whether it provides an approximate solution to the entire problem. Equivalently, we can define for every electronic eigensolution the direct integral operator, π ^ n : = ∫ ℝ ⊕ d q   ξ n ( q ) 〈 ξ n ( q ) , ⋅ 〉 e = ∫ ℝ ⊕ d q   π n ( q ) , (8) on the total Hilbert space ℋ = L 2 ( ℝ , d q ) ⊗ ℋ f and ask whether it commutes with H ^ . Of course, the answer is in the negative, but it turns out that the commutator scales like ε, [ H ^ , π ^ n ] ∼ ε . (9) This is because of the adiabatic relation between the electrons and the nuclei. By construction, the Weyl quantization H ^ e ( q ^ ) of the electronic Hamiltonian and π ^ n commute. However, the remaining contribution to the Hamilton operator H ^ , in particular the kinetic energy of the nucleus, scales like ε 2 and leads hence to the estimate in Eq. 9. The Born–Oppenheimer approximation builds on this result and proposes to use the ansatz functions, Ψ ( q ; x ) = ∑ n ψ n ( q ) ξ n ( q ; x ) , (10) to solve the full quantum problem. In its simplest version, the scheme neglects any of the contributions that arise from applying the kinetic energy operator of the nucleus to the electronic ansatz functions (as they enter with small ε-factors), and thus results in an effective eigenvalue problem for the nucleus only, ( P ^ 2 2 m n + E n ( q ^ ) ) ψ n ( q ) ≡ E ψ n ( q ) . (11) If this nucleonic eigenproblem can be solved, the scheme leads in fact to viable results for the stationary energy spectra of molecules, which are given by the energy solutions E, and which include the backreaction of the electrons via the potential energy E n ( q ^ ) .

Unfortunately, the Born–Oppenheimer approach comes with some limitations which preclude its application to more complicated systems. Firstly, the scheme explicitely uses that the electronic eigenfunctions depend only on the configuration variable of the nucleus. Was the coupling between electrons and nuclei provided by non-commuting slow operators, for example by q ^ and p ^ , the scheme would fail since the direct integral construction in Eq. 8 builds on the commutativity (i.e., the existence of a common spectrum) of the coupling operators. Also one could not define the ansatz functions in Eq. 10. Secondly, the scheme does not provide a simple extension to better error estimates. This becomes problematic if one is interested in the dynamical evolution of the system. The interesting dynamics of the nuclei happens on time scales t n ∼ t e / ε or larger, but considering the evolution generated by H ^ with respect to the above ansatz functions, the scheme cannot lead to trustworthy results due to the commutator relation Eq. 9.

It should be possible to do better. In fact, the adiabatic theorem (Teufel, 2003) states that under certain conditions (to be discussed in the sequel), there exists an orthogonal projection operator Π ^ ∈ ℬ ( ℋ ) in the bounded operators on the total Hilbert space ℋ such that, [ H ^ , Π ^ ] = O 0 ( ε ∞ ) , (12) where the right hand side means that for all m ∈ ℕ , there exists a constant C m < ∞ such that ‖ [ H ^ , Π ^ ] ‖ ℬ ( ℋ ) ≤ C m ε m , in the norm of bounded operators on ℋ . Most importantly, Π ^ can be constructed by a “semiclassical symbol” function, i.e., an operator-valued ansatz function like the almost–projector function π n ∈ C ∞ ( ℝ , ℬ ( ℋ f ) ) from the simple example above. This symbol function appears as an asymptotic series in the perturbation parameter ε, and–to anticipate the result–the equivalent of π n ( q ) will serve as the base clause to an iterative ε-scheme to compute better and better approximations to Π ^ . This is the idea of SAPT (Panati et al., 2003).

SAPT uses an ε-scaled phase space (or deformation) quantization scheme (Blaszak and Domanski, 2012) for the slow subsector of the system, while retaining a standard Hilbert space representation for the fast sector. Phase space quantum mechanics is a formulation of quantum mechanics that employs an algebra of phase space functions A Q instead of using the standard operator algebra in the Hilbert space representation of quantum mechanics. Quantum mechanical observables are thus represented by real-valued phase space functions. The pullback of the operator product to the phase space algebra gives rise to a non-commutative star product ⋆ . Since the star product reduces to the commutative multiplication of phase space functions in the limit ℏ → 0 , this formulation of quantum mechanics is also known as deformation quantization. The standard textbooks by Dimassi and Sjöstrand (1999), Folland (1989), and Hörmander (1985a,b) give thorough introductions to the usual scalar–valued phase space quantization scheme and pseudodifferential calculus, but the situation here is more subtle. SAPT requires to consider operator–valued symbol functions, in particular functions on the slow phase space with values in the operators on the fast Hilbert space. We will thus deal with symbol functions of the form A ( q , p ) ∈ C ∞ ( Γ s , ℬ ( ℋ f ) ) where Γ s denotes the slow phase space (Teufel, 2003; Appendix A). It is straightforward to map the symbol functions (operator–valued or not) to their operator representatives in the standard Hilbert space approach. The concrete prescription depends of course on the operator ordering that one chooses. In case of the symmetric Weyl quantization prescription, this relation is provided by the Weyl correspondence (Dubin et al., 1980), and a symbol function appears as the kernel of an integral operator that acts on an element of the Hilbert space (Teufel, 2003).

Symbol functions which give rise to admissible operators in the quantum theory can be classified by their asymptotic behavior on the slow phase space. One important class of symbols relevant for SAPT are the semiclassical symbols S ρ m , with m ∈ ℝ and 0 ≥ ρ ≥ 1 . An operator–valued function A ∈ C ∞ ( ℝ 2 , ℬ ( ℋ f ) ) is in the symbol class S ρ m ( ℬ ( ℋ f ) ) if for every α, β ∈ ℕ , there exists a positive constant C α , β such that, sup q ∈ ℝ ‖ ( ∂ q α ∂ p β A ) ( q , p ) ‖ ℬ ( ℋ f ) ≤ C α , β 〈 p 〉 m − ρ | β | , (13) for every p ∈ ℝ , and 〈 p 〉 = ( 1 + | p | 2 ) 1 / 2 (Teufel, 2003). For such symbols the Weyl ordering prescription for quantum theory gives rise to a specific star product, and we can finally make sense of the space adiabatic perturbation idea. For two such operator–valued symbols on the slow phase space A ∈ S ρ m 1 ( ℬ ( ℋ f ) ) , B ∈ S ρ m 2 ( ℬ ( ℋ f ) ) , their star product is given by, ( A ⋆ B ) ( q , P ) = exp ( i ℏ 2 ( ∂ x ∂ P − ∂ q ∂ D ) ) A ( x , D )   B ( q , P ) | x = q , D = P ∈ S ρ m 1 + m 2 ( ℬ ( ℋ f ) ) . (14) We note that the exponential has a series expansion which could be considered as a series with respect to ℏ in the given context. Alluding to the adiabaticity of the system, it is reasonable to define a rescaled momentum operator p : = ε P . Replacing P and D by their ε-scaled versions in the star product formula, the new expansion parameter is ε ℏ . Any star product of symbol functions can thus be written in a series expansion in ε. Comparing terms of the same polynomial order in ε, this defines a perturbation theory for quantum mechanical equations which will iteratively solve the eigenvalue problems of interest. The first two orders of the rescaled star product are given by, ( A ⋆ B ) ( q , p ) = ( A 0 ⋅ B 0 ) ( q , p ) + i ε ℏ 2 { A 0 ( q , p ) , B 0 ( q , p ) } + ε ( A 0 ⋅ B 1 + A 1 ⋅ B 0 ) ( q , p ) + O ( ε 2 ) , (15) where the symbol functions have been expanded with respect to ε according to A = ∑ k ε k A k and B = ∑ k ε k B k . Note that the ε-rescaling changes the whole symplectic structure (we now have, { q , p } s = ε ), as well as the canonical commutation relations since we obtain, [ q ^ , p ^ ] s = i ε ℏ   1 ^ s . (16) From now on, we will set ℏ ≡ 1 . Note that in the original Born–Oppenheimer approximation, the perturbative parameter occurs (after an appropriate rescaling) only in the Hamiltonian and the perturbation theory consists in splitting the Hamiltonian (and its spectrum) accordingly, while here the quantum algebra is redefined, giving rise to a (n in principle infinite) perturbation series in ε.