1 Introduction↩︎

In general relativity, matter fields propagate on a dynamical background geometry, with the matter fields causing spacetime to curve via Einstein’s equation. Spacetime thus serves in a dual role; as a background for defining causal structures for quantum matter fields, as well as a dynamical variable influenced by those same fields. Whereas quantum field theory in curved space is well understood, it is not clear how to suitably modify general relativity so as to incorporate the back-reaction of matter fields which are now quantum in nature. Various proposals have been put forth on the strength of different underlying philosophical considerations. The dominant approach is quantum gravity, including string theory, loop quantum gravity and asymptotic safety. There are also a number of pragmatic approaches, such as mean field gravity approach (sometimes called “semi-classical gravity”), which involves sourcing the geometry from the expectation value of the energy momentum tensor. This approach has a domain of applicability, but cannot be considered as a fundamental theory, since it suffers from various pathologies including superluminal signaling [@gisin1989stochastic] and a breakdown in the statistical interpretation of the density matrix [@page1981indirect]. While the dominant view is that gravity is fundamentally quantum, the special status that spacetime holds within quantum field theory obliges us to consider the hypothesis that we should treat it as fundamentally classical.

While there has been considerable debate about whether there are consistent ways of coupling classical and quantum systems, a number of examples have been demonstrated since the mid 90’s, including via a master-equation [@blanchard1993interaction; @diosi1995quantum], and via a measurement and feedback approach [@diosi1998coupling; @kafri2014classical; @tilloy2016sourcing]. In [@oppenheim2023postquantum; @oppenheim2022two] the most general form of consistent classical-quantum dynamics was derived. The resulting dynamics preserves positivity and normalization of probabilities, and maps the state of the classical subsystem to a classical state and the state of the quantum subsystem to a quantum one. At longer distance scales, the back-reaction onto the classical degrees of freedom act as if the quantum systems is being weakly measured, resulting in decoherence, while at shorter distance scales, the dynamics can be highly non-Markovian with the classical degrees of freedom evolving stochastically. The two parameters (governing diffusion and decoherence) are thus related to each other via an inequality termed the ‘diffusion-decoherence’ tradeoff[@diosi1995quantum; @oppenheim2023gravitationally]. That any consistent classical-quantum coupling must be fundamentally irreversible has also been established using general probabilistic theories [@galley2023any].

The measurement and feedback approach has been used to describe Newtonian gravity[@kafri2014classical; @tilloy2016sourcing; @tilloy2017principle], and a theory coupling quantum fields to general relativity has been presented via a master-equation approach[@oppenheim2023postquantum], trajectories approach[@layton2024healthier] with Layton and Weller-Davies, and a manifestly covariant classical-quantum path-integral with Weller-Davies[@oppenheim2023covariant]. In [@grudka2024renormalisation] the formal renormalizability of the gravitational path integral was demonstrated, which stands in contrast to the non-renormalisability of perturbative quantum gravity. This further motivates probing the quantum vs classical nature of spacetime via experiment. The aforementioned decoherence-diffusion tradeoff allows one to tightly constrain the space of feasible experimental bounds [@oppenheim2022decoherence; @janse2024current], allowing one to possibly rule out classical gravity if the trade-off is experimentally found to be violated. Proposed tests which would witness the quantum nature of gravity include gravitationally-mediated entanglement experiments [@bose2017spin; @marletto2017gravitationally; @carney2019tabletop], non-entanglement witnesses of quantum gravity [@lami2024testing; @kryhin2025distinguishable].

One of the aims of the present work is to reformulate the theory in such a way as to constrain its parameters via observation. Constraints coming from tabletop interference experiments[@gerlich2011quantum] and precision gravity measurements such as LISA Pathfinder [@armano2024indepth; @janse2024current] have been derived via the Newtonian limit of the theory[@oppenheim2023gravitationally; @layton2023weak; @grudka2024renormalisation] via the decoherence-diffusion trade-off. Constraints coming from the production of stochastic gravitational waves have been estimated via dimensional analysis[@penington2024], and from using the stochastic Klein Gordon equation with Panella as a model[@oppenheim2025diffusion]. Stochastic violation of the Hamiltonian constraint has been shown to lead to the appearance of phantom cold dark matter with Panella and Pontzen in a Friedmann–Lemaître–Robertson–Walker model[@oppenheim2024emergence], which may give constraints from considering stochastic fluctuations as a candidate for dark matter.

Contemporaneous studies [@hirotani2026testing; @fabiano2026minimal] have also constrained its parameters for various models. More generally however, understanding relativistic extensions of hybrid dynamics, particularly in the case of gravity has been a challenge. Much work has been done on relativistic toy models [@oppenheim2023covariant; @oppenheim2024diffeomorphism; @diosi2022there; @grudka2024renormalisation; @carney2024classicalquantum] but the full theory, or even its linearized versions have additional subtleties.

A significant issue is one of interpretability. Historically, there have been two different approaches towards the study of gravity; the original, geometro-dynamical one, and the particle physicist one where, starting from a Lagrangian describing a massless spin-2 particle around a fixed spacetime background, one is able to recover general relativity. Thus, the fundamental object for particle physicists in the path integral formulation of quantum gravity (with the Einstein-Hilbert action) is the spin-2 particle, with its two degrees of freedom. However, this is not necessarily the case in alternative theories of gravity; for example all the degrees of freedom of the spin-2 particle are excited in massive gravity, whereas in quadratic gravity one encounters a spin-0 particle as well as an additional ghost particle. What then, are the degrees of freedom for classical-quantum gravity, and how do they relate to the Newtonian potential? The first objective of the present paper, is to understand what the degrees of freedom are.

Secondly, and not unrelatedly, we address the main obstacle to writing a classical gravitational path integral with the interpretation of a transition probability amplitude: the fact that the gravitational action we consider is unbounded from below. Were the unstable mode to represents a physical degree of freedom, it would leave us in the unenviable position of perturbing around a saddle-point, and would result in an indefinite two-point function. Since the two-point function has the interpretation of a covariance matrix in a classical stochastic theory, it must be positive semi-definite (PSD). This suggests a tension in the theory that must be resolved. The problem of unboundedness is not without precedent: in Euclidean Einstein-Hilbert gravity, the conformal mode has a wrong-sign kinetic term, rendering the action unbounded from below and the naive path integral divergent [@gibbons1978path]. Various contour and measure prescriptions have been proposed; in particular, Mazur and Mottola argued that the conformal mode is tied to the constraint structure of the theory and does not propagate, so that the physically relevant spin-2 part of the action is positive definite [@MAZUR]. Whether this fully resolves the conformal factor problem remains an active question [@marolf2022gravitational].

Here we show that if we perturb around flat spacetime, there exists a gauge invariant positive semi-definite sector that describes a massless, dynamical, classical spin-2 field undergoing diffusion. Additionally, we obtain a new dynamical, stochastic spin-0 field. We then address the issue of the unbounded sector, which corresponds to an off-shell vector mode, and show that it does not correspond to degrees of freedom, but to the momentum constraint. We then argue that it is normalized away.

Thirdly, we use our understanding of the different stochastic modes of the path integral, to place bounds on the parameters of the theory from observational constraints. We relate the spin-0 and the spin-2 fields to the Newtonian potential, and compute the two point function for the Newtonian field 16 29 as well as the power spectral density of the acceleration, Eqn 30 . This is then used to place an upper bound on the diffusion constants, using results from the Lisa Pathfinder[@Armano_2024]. This improves upon results from [@grudka2024renormalisation] through an analysis of how the different coupling constants and scalar fields contribute. Another upper bound, is obtained in Eq. 35 by computing the energy-density of stochastic gravitational waves. We then obtain a lower bound from the minimum rate of decoherence 5.4, and show that the theory is consistent with current observational constraints. We note that the decoherence rate 5.4 is non-Markovian due to coupling with the curvature scalar field \(\psi\).

The paper is structured as follows: in Section 2 we briefly review the classical quantum path integral. In Section 3 we give the two point function in terms of Barnes-Rivers projectors onto the spin-0 and spin-2 modes, and in Section 4 in terms of a Scalar-Vector-Tensor (SVT) decomposition of the metric. We show that the dynamical degrees of freedom have positive semi-definite two point functions. In particular, the action is PSD for the tensor mode corresponding to gravitational waves, and for the scalar sector. This is the condition we require in order to consistently interpret the action as giving transition probabilities for classical stochastic fields. We distinguish these dynamics modes of the theoryfrom those found in the Einstein Hilbert action and in Quadratic Gravity [@stelle1977renormalization],[@salvio2018quadratic],[@alvarez-gaume2016aspects],[@hell2023degreesfreedomr2gravity] in Subsection 4.1. In Subsection 4.2 we show that the only part of the action which is not PSD, corresponds to a constraint equation for the vector mode. As a result the path integral for the vector mode does not give transition probabilities, and we argue that this sector of the action is normalised away, resulting in a normalisable and well-behaved path integral.

In Section5 we constrain the parameter space of the theory using gravitational observations, as well as table-top decoherence experiments. Our bounds from gravitational waves is consistent with prior toy models as well as the study of Hirotani and Matsumura [@hirotani2026testing] which appeared while this draft was being prepared. We are additionally able to distinguish bounds which apply to the propagating modes, from those which place bounds on the constraints. Hirotani and Matsumara also introduce a noise model whose vacuum sector gives a theory related to that found in stochastic gravity[@hu2008stochastic; @moffat1997stochastic] and thus can describe a theory which is designed to be fundamentally quantum, but with a domain of validity which has a consistent and non-local classical limit.

Inspired by their work, we demonstrate that the Einstein tensor produced by our local path integral satisfies the Bianchi identity. We also resolve a discrepancy between the MSR and Onsager-Machlup formalisms of the theory, as well as in the stochastic differential equation (SDE) formalism. This is done in the Appendix. In Subsection7, we discuss the Martin Siggia Rose (MSR) form of the theory. In subsections 7.1, 7.2,7.3 we review results on the gravitational path integral, as well as on classical stochastic path integrals in general, while in subsections7.4,7.5 we review calculational details for Subsection5.2. Finally Subsection7.9 of the appendix defines the Barnes-Rivers projectors that we heavily use.

1.0.1 Conventions↩︎

We use the mainly plus \((-,+,+,+)\) metric, and natural units, \(c=\hbar=1\). Expressions of the form \(\frac{1}{\Box},\frac{1}{\triangle}\) are to be understood as the full and spatial Fourier transforms respectively, of these operators.

2 The Classical Quantum Path Integral↩︎

A classical-quantum system may be described by a quantum density matrix living in a Hilbert space, conditioned on a classical system being at a point in classical phase space. The state of the classical system may be represented by its probability distribution, \(p(z,\dot{z,}t)\), where we obtain \(p(z,\dot{z},t)\) by tracing over the Hilbert space of the classical quantum (CQ) state, \(\hat{\varrho}(z,\dot{z},t)\). The classical degrees of freedom can also be taken to live in phase space. The dynamics must by consistency preserve the phase space and map a CQ state to another, ensuring that classical degrees of freedom are mapped to classical degrees of freedom, and the quantum subsystem to a quantum subsystem. Since the CQ state is positive and normalized, the dynamics must be completely positive, and trace-preserving (CPTP). In [@oppenheim2022path] a path integral formulation for such dynamics was introduced, the key features of which we briefly summarize below

The classical-quantum (CQ) state of a system is given by \[\begin{align} \hat{\varrho}(z,\dot{z},t):=\int d\phi^+ d\phi^- \varrho(\phi^+,\phi^-,z,\dot{z},t)|\phi^+\rangle\langle\phi^-| \end{align}\] with the components, \(\varrho(\phi^+,\phi^-,z,\dot{z},t):=\langle\phi^+|\hat{\varrho}(z,\dot{z},t)|\phi^-\rangle\), where \(\phi^+,\phi^-\) represent the bra and ket fields respectively and can depend on \(z,\dot{z},t\). The evolution of these components is governed by a path integral: \[\begin{align} \varrho(\phi^{+},\phi^{-},z_f,\dot{z_f},t_f)= \frac{1}{\mathcal{N}} \int D\phi^{+}D\phi^{-}Dz e^{-\mathcal{I}_{CQ}}\varrho(\phi^{+},\phi^{-},z_i,\dot{z_i},t_i) \end{align}\] where \(\frac{1}{\mathcal{N}}\) is a normalization constant and \(\mathcal{I}_{CQ}\) is of the following form: \[\begin{align} \mathcal{I}_{CQ}=-i\mathcal{I}_Q(\phi^+,\phi^-)+\mathcal{I}_{FV}+\mathcal{I}_{Diff} \end{align}\] Here \(i\mathcal{I}_Q(\phi^+,\phi^-)\) governs the unitary dynamics of the bra and ket fields, \(\mathcal{I}_{FV}\) are Feynman Vernon, or Schwinger-Keldysh terms which cause decoherence, and \(\mathcal{I}_{Diff}\) is the diffusive part of the action which in general incorporates quantum back-reaction on a classical field. Actions containing the first two terms, \(i\mathcal{I}_Q(\phi^+,\phi^-)+\mathcal{I}_{FV}\) are used in a formalism known as the Schwinger-Keldysh or Feynman-Vernon  [@feynman1963theory], and are widely used in the study of open quantum systems. As a simple example, we may consider the following action [@grudka2024renormalisation], \[\begin{align} \mathcal{I}_{CQ}=\int^{t_f}_{t_i} dt-i((\dot{\phi}^+)^2-(\dot{\phi}^-)^2)-\frac{D_0}{8}((\phi^+)^2-(\phi^-)^2)^2-\frac{1}{2D_2}(\ddot{z}-\frac{D_1}{2}((\phi^+)^2+(\phi^-)^2))^2 \end{align}\] Here, \(\phi\) and \(z\) can also be taken to be fields in which case \(dt\rightarrow d^4x\). The imaginary part corresponds to \(i\mathcal{I}_Q(\phi^+,\phi^-)\) while the term \(\frac{D_0}{8}((\phi^+)^2-(\phi^-)^2)^2\) is the Feynman-Vernon term which acts to decohere the field with a strength governed by \(D_0\). The last term is a classical quantum version of the Onsager-Machlup (OM) functional, a probability density functional used to determine transition probabilities for stochastic processes [@onsager1953fluctuations], [@Chaichian2018PathI] and has the form of the square of an equation of motion. As \(D_2\) decreases, deviations from the equations of motion are suppressed. These deviations are stochastic kicks to the acceleration \(\ddot{z}\). If \(D_2\) is not positive (or PSD in the case of a matrix of equations), then the action is unbounded, and paths which have large deviations from the equations of motion dominate. Finally, \(D_1\) controls the strength of the back-reaction.

This generalizes to the field theoretic case, where the primary objects are the classical and quantum field  [@oppenheim2023covariant]. In the gravitational case, we have the following action: \[\label{eq:PQG-action} \begin{align} \mathcal{I}_{CQ}[ g_{\mu \nu}, T^{\mu\nu+}, T^{\mu\nu-}]& = \int d^4x \bigg[ - i\text{Det}[-g]\big(\mathcal{L}_{Q}[T^{\mu\nu+}] - \mathcal{L}_{Q}[T^{\mu\nu-}]\big) -\frac{\text{Det}[-g]}{8}( T^{\mu \nu + } - T^{\mu \nu - }) D_{0,\mu \nu \rho \sigma}(T^{\rho \sigma+ } - T^{\rho \sigma- }) \\ & - \frac{\text{Det}[-g]}{2}\bigg( G^{\mu \nu} +\Lambda g^{\mu\nu} - 8\pi G_N\frac{T^{ \rho \sigma +} + T^{ \rho \sigma -}}{2}\bigg)D_{2, \mu \nu \rho \sigma}^{-1}\bigg( G^{\rho \sigma} +\Lambda g^{\rho\sigma} - 8\pi G_N\frac{T^{ \rho \sigma +} + T^{ \rho \sigma -} }{2}\bigg)\bigg], \end{align}\tag{1}\]

This action governs the time evolution of the system through \[\begin{align} \varrho(T^{+},T^{-},g,t_f)= \frac{1}{\mathcal{N}}\int D\varphi^{+}D\varphi^{-}D'g e^{-\mathcal{I}_{CQ}}\varrho(T^{+},T^{-},g(z_i,\dot{z_i},t_i)) \end{align}\] where the stress-energy tensor \(T^+,T^-\) is a function of \(\varphi^+,\varphi^-\), and \(G\) is the Einstein tensor for spacetime metric \(g\), and the gravitational path integral is over continuous configurations, up to diffeomorphisms7.1. Let us break down the action in its constituent parts. The imaginary part of the action corresponds to the unitary action for the stress energy tensor, in component form, \(T^{\mu\nu+},T^{\mu\nu-}\), while the term \(( T^{\mu \nu + } - T^{\mu \nu - }) D_{0,\mu \nu \rho \sigma}(T^{\rho \sigma+ } - T^{\rho \sigma- })\) is a Feynman-Vernon term which acts to decohere the density matrix by causing off-diagonal terms to exponentially decay. The last term represents a classical-quantum version of the Onsager-Machlup (OM) functional. The classical OM functional is of the form, equations of motion (EOMs) squared, which in this case are Einstein’s equations. The classical-quantum version, has the dynamics of the classical system \(G^{\mu\nu}\) sourced by the average \(\bar T^{\mu\nu}:=\frac{1}{2}(T^{\mu\nu+}+T^{\mu\nu-})\) of the bra and ket fields. The classical field thus, diffuses around its equations of motion due to a stochastic force, with the covariance of the stochastic force given by the diffusion matrix \(\mathcal{D}_2\), with the components, \(D_{2}^{\rho\sigma,\mu\nu}\). One could instead write a more general, possibly non-local diffusion matrix, \(\mathcal{D}_2'(x,y)\). In this work however, we will stick to a local diffusion matrix of the form \(\mathcal{D}_2(x,x'):=\mathcal{D}_2\delta^4(x-x')\) for reasons of general covariance.

Additionally, we can characterize the strength of the classical-quantum interaction in the OM functional by the matrix, \(\mathcal{D}_1=8\pi G_N\). If we now set \(\mathcal{D}_2=\mathcal{D}_1\mathcal{D}_0^{-1}\mathcal{D}_1\), as we will do moving forward, no cross terms of the form \(T^{+}T^{-}\) survive, and quantum states retain their purity, conditioned on the classical trajectory. The bra states evolve independently of the ket fields. This is an example of saturating the inequality which we call the diffusion-decoherence trade-off, \(\mathcal{D}_0\geq 64\pi^2G_N^2\mathcal{D}_2^{-1}\) [@oppenheim2023covariant]. As long as we are to have a completely positive dynamics, this trade-off must be preserved[@oppenheim2023gravitationally].

Our primary object will be the inverse diffusion matrix, which will determine our action. Now we would like \(\mathcal{D}_{2}^{-1}(x,x')\) to have the following properties

1. Positivity: The matrix must be positive semi-definite (PSD) so that it is a bona-fide covariance matrix. Covariance matrices must be PSD i.e., have positive eigenvalues, which follows from the fact that the variance is the sum of squares, and must be positive. Another way to see this is that the matrix must be PSD in order to suppress deviations from the equations of motion in the OM action.

2. General covariance: This implies that the diffusion matrix does not introduce any structure additional to the background metric and is only delta-correlated in time and space, i.e., through \(\delta^4(x-x')\). This implies two conditions, (a), that the noise kernel is ultra-locality so that the noise contains no spatial or temporal derivatives in its covariance matrix. Secondly, (b) that it determines the geometry by determining the metric up to gauge transformations.

One of the proposals of [@oppenheim2023covariant] was to take it to be proportional to the Generalised DeWitt metric[@dewitt1967quantum] \[\begin{align} \label{eq:dewitt} \mathcal{D}_2^{-1 \rho\sigma,\mu\nu}=\frac{1}{\sqrt{-g^{}}}(\alpha\frac{g^{\rho\mu }g^{\sigma\nu}+g^{\rho\nu}g^{\sigma\mu}}{2}-\beta g^{\mu\nu}g^{\rho\sigma})\delta^{4}(x-x') \end{align}\tag{2}\] which is indefinite for non-zero \(\alpha\), and so on its surface, does not appear to satisfy condition (1). This may be contrasted with two proposals: (i), that \(\alpha\) is set to zero; this gives an infinite diffusion rate for the spin-2 sector, (ii), that \(C^2\), the square of the Weyl-tensor, is set to evolve deterministically; this then leads, for example, to an infinite decoherence rate for tensorial matter. Since both of the alternative proposals lead to infinite rates of diffusion or decoherence, we will continue with the DeWitt metric and show that it is indeed positive on all dynamical diffusive modes. We will further argue that the indefinite part is constrained to vanish in a gauge invariant way.

In Appendix 7, we will consider the stochastic differential equation form of Einstein’s equations, and will introduce another diffusion matrix for which we may relax the third condition, but which gives us the same action as below.

3 Spin-Projector Representation↩︎

Here, we explicitly compute that two point function around the Minkowski background for our choice of diffusion matrix and show that it consists of spin-0 and spin-2 parts. The gravitational two point function must be positive semi-definite if it is to indeed be a genuine covariance matrix, and here we find that the spin-0 part is PSD, but the spin-2 part is only PSD for non-tachyonic modes. The consequences of this indefiniteness will be discussed in 4.1 where we argue that it is benign.

The purely gravitational ‘action’ is given by \[\label{eq:LGR}\mathcal{S}^{(g)}=\frac{1}{2}\int d^4xd^4y\sqrt {-g} G_{\mu\nu}(x)(\mathcal{D}_2^{-1})^{\mu\nu,\rho\sigma}\delta^4(x-y)G_{\rho\sigma}(y)=\frac{1}{2}\int d^4x\sqrt -g (\alpha R_{\mu\nu}R^{\mu\nu}-\beta R^2)\tag{3}\] . This action is rather similar in form to that of quadratic gravity[@stelle1978classical] [@salvio2018quadratic] \[\begin{align} \label{eq:QuadGrav} \mathcal{S}_{QG}=\frac{1}{2}\int d^4x\sqrt{-g}( \alpha R^{\mu\nu}R_{\mu\nu}-\beta R^2+\gamma R^{\mu\nu\rho\sigma} R_{\mu\nu\rho\sigma}-\frac{M^2}{2}R ; \end{align}\tag{4}\] where \(R^{\mu\nu\rho\sigma}\) is the Riemann tensor. Since one combination of the terms above give us the topological Gauss-Bonnet term, which is a total derivative, 4 is equivalent to 3 up to the \(R\) term. However, the presence of that term, which corresponds to the Einstein-Hilbert term in quadratic gravity, and to the cosmological constant in classical-quantum gravity, spoils the PSD-ness of the spin-0 sector, and will thus be omitted.

In addition, we also note that \(C^2:=C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}\), Weyl-squared term, may be written as \(\frac{1}{2}C^2=\frac{1}{2}R^{\mu\nu\rho\sigma}R_{\mu\nu\rho\sigma}-R_{\mu\nu}R^{\mu\nu}+\frac{1}{6}R^2\), or up to the Gauss-Bonnet term as \(\frac{1}{2}C^2=R_{\mu\nu}R^{\mu\nu}-\frac{1}{3}R^2\). Note then that taking the limit, \(\alpha\rightarrow3\beta\) in 3 gives an action similar in form to that of pure Weyl-squared gravity.

Expanding 3 around Minkowski space to second order in \(h\), we obtain the following Lagrangian density: \[\begin{align} \label{eq:Mink} \mathcal{L}^{(g)}&= \frac{\alpha}{4}(\partial_{\lambda}\partial_\mu h^{\lambda}_{\nu}+\partial_{\lambda}\partial_\nu h^{\lambda}_{\mu}-\partial_\mu \partial_\nu h-\Box h_{\mu\nu})(\partial^{\sigma}\partial^\mu h_{\sigma}^{\nu}+\partial^{\sigma}\partial^\nu h_{\sigma}^{\mu}-\partial^\mu \partial^\nu h-\Box h^{\mu\nu})-\beta(\partial_{\mu}\partial_{\nu}h^{\mu\nu}-\Box h)^2\\&=h^{\mu\nu}[\frac{\alpha}{2}(\frac{1}{2}I_{\mu\nu,\rho\sigma}\Box^2+\frac{1}{2}\eta_{\mu\nu}\eta_{\rho\sigma}\Box^2-\eta_{\nu\beta}\Box\partial_{\rho}\partial_{\mu}-\frac{1}{2}\Box\eta_{\mu\nu}\partial_{\rho}\partial_{\sigma}-\frac{1}{2}\Box\eta_{\rho\sigma}\partial_{\mu}\partial_{\nu}+\partial_{\rho}\partial_{\sigma}\partial_{\mu}\partial_{\nu})\\&-\beta(\eta_{\rho\sigma}\eta_{\mu\nu}\Box^2-\eta_{\mu\nu}\Box\partial_\rho\partial_\sigma-\eta_{\rho\sigma}\Box\partial_{\mu}\partial_\nu+\partial_{\mu}\partial_{\nu}\partial_{\rho}\partial_\sigma)]h^{\rho\sigma}\\&=h^{\mu\nu}H_{\mu\nu,\rho\sigma}h^{\rho\sigma} \end{align}\tag{5}\] It is convenient to write \(\mathcal{H}\) in terms of the Barnes Rivers spin projectors 7.9: \[\begin{align} \label{eq:H}H_{\mu\nu,\rho\sigma}=\Box^2(\frac{\alpha}{4}(P^{(2)}+4P^{(0-s)})-3\beta P^{(0-s)})_{\mu\nu,\rho\sigma}=\Box^2(\frac{\alpha}{4}P^{(2)}+(\alpha-3\beta)P^{(0-s)})_{\mu\nu,\rho\sigma} \end{align}\tag{6}\] The three dimensional version of the spin-2 projector, \(P^{ij,kl}\) 22 is a familiar sight in classical gravity where it is used to project out the transverse-traceless part of the spatial part of the metric perturbation, \(h_{ij}\), to give us the two transverse,traceless modes corresponding to gravitational waves. In 4-dimensions, one obtains instead five transverse traceless modes, which one may write in terms of their polarization tensors [@Spin2]. After gauge-fixing and inverting 6 , one has, up to gauge-fixing terms, the following two-point function: \[\begin{align} \label{eq:Twopointfunction} \langle h_{\rho\sigma }h_{\mu\nu}\rangle=\frac{1}{\alpha\Box^2}(P^{(2)}+\frac{P^{(0-s)}}{4(1-3\frac{\beta}{\alpha})})_{\rho\sigma,\mu\nu} \end{align}\tag{7}\] One can easily check that the spin-0 part is positive definite for \(\alpha>3\beta\)¹. Explicitly, \[\frac{1}{(\alpha-3\beta)\Box^2}S^{\mu\nu} P^{(0s)}_{\mu\nu,\rho\sigma}S^{\rho\sigma} = \frac{1}{3(\alpha-3\beta)}\,(\frac{S^{\mu\nu}\theta_{\mu\nu}}{\Box})^2\]

For the spin-2 part to be PSD, we require \[\begin{align} S_{\rho\sigma}P^{\rho\sigma,\mu\nu(2)}S_{\mu\nu}\geq 0 \label{eq:32PSD} \end{align}\tag{8}\] for all \(S_{\rho\sigma}\) Acting on a tensor \(S^{\mu\nu}\) with the projection operator takes us to the (symmetric) transverse traceless subspace, i.e., \(\boldsymbol{P}S=S'\), \[\begin{align} S_{\rho\sigma}P^{\rho\sigma,\mu\nu}S_{\mu\nu}=S'_{\rho\sigma}S'^{\rho\sigma} \label{eq:32PSD2} \end{align}\tag{9}\] \[\begin{align} \label{eq:32transverse} p_\alpha S'^{\rho\sigma}=0, \quad \eta^{\rho\sigma}S'_{\rho\sigma}=0 \end{align}\tag{10}\] One can easily show for any 2 tensor that satisfies the above relation for all \[\begin{align} \label{eq:spacelike} p^2\leq0 \end{align}\tag{11}\] that the quadratic invariant, \[\begin{align} S_{\mu\nu}'S'^{\mu \nu} \geq 0 \end{align}\] To see this, first, note that for any tensor \(\boldsymbol{S}\) in Minkowski space, \[\begin{align} S'^{\mu\nu}S'_{\mu\nu}=(S'^{00})^2-2\sum_i(S'^{0i})^2+\sum_{ij}(S'^{ij})^2 \end{align}\]. Thus, it is the \(S^{'0i}\) term which could potentially spoil PSD-ness. By 10 and 11 there exists a coordinate system in which \(p=(1,0,0,0)\), and for that co-ordinate system, \(S'^{0\mu}=0\), and so \(S'^{\mu\nu}S'_{\mu\nu}\geq0\) in any coordinate system. For a null vector \(p\), the Barnes-Rivers projectors seem to be undefined. However, as long as the source tensor \(S^{\mu\nu}\) is conserved, i.e., \(p_\mu S^{\mu\nu}=0\) the projector is in fact well-defined. Again, we can work in a coordinate system where \(p=(1,0,0,\pm 1)\). Then, we see that \[\begin{align} \begin{aligned} S'^{0i}&=\mp S^{'3i}\\S^{00}\mp &=S^{30}\\ (S'^{00})^2+2\sum_i(-(S'^{0i})^2+(S'^{3i})^2)-(S^{30})^2+\sum_{i,j\neq3}(S'^{ij})^2&= \sum_{i,j\neq3}(S'^{ij})^2\ge0 . \end{aligned} \end{align}\]On the other hand, if we wish to remove this restriction from the source tensor, the projection operator will have to be suitably altered to ensure transversality: \[\theta^{\mu \nu}=\eta^{\mu\nu}-p^{\mu}n^{\nu}/p^{\alpha}n_\alpha-p^{\nu}n^{\mu}/p^{\alpha}n_\alpha\] where \(n^2,p^2=0\)

One can also check that this method does not work for spacelike p. Setting, \(p=(0,0,0,1)\) in some coordinate system, we have \(p_\alpha S'^{\rho\sigma}=S^{3\beta}=0\). Thus, now the best we can say is that \(S'^{30}=0\). More generally, one can only use a space-like vector to draw relations between the various components of \(S'^{0i}\), and so can say nothing about the relative size of the troublesome term, \(\sum_i(S'^{0i})^2\), which may be arbitrarily large. Thus, in this case, \(S'^{\mu\nu}S'_{\mu\nu}\) is indefinite.

Thus, for linearized gravity, the spin-2 sector has a positive semi-definite correlation function only for timelike, and null Fourier modes, unlike the spin-0 sector, which is positive semi-definite for all modes. However, in what follows, we will show using a scalar-vector-tensor decomposition that the actual degrees of freedom (which we will define in due course) have two point functions that are PSD.

4 The Scalar-Vector-Tensor Decomposition↩︎

In this section, our aim will be to identify the true degrees of freedom. To this object we will perform a standard scalar-vector-tensor decomposition of the metric, and consequently of the action, and show that the scalar and tensor sectors have PSD two-point functions. In addition, we will show that the indefinite vector sector, in our framework, does not represent a degree of freedom.

For performing the scalar-vector-tensor decomposition around Minkowski space, we write the perturbed metric as \[g_{\mu\nu}^{(1)}= -(1-2\phi)dt^2 + (B_{,i}+S_i)dtdx^{i}+((1+2\psi)\delta_{ij}+2E_{,ij}+F_{i,j}+F_{j,i}+h_{ij}^{TT})dx^{i}dx^{j}\] where the vector modes, \(S_i,F_{i}\) are transverse, while the tensor modes, \(h_{ij}^{TT}\) are both transverse and traceless. For infinitesimal coordinate changes, \(x^{\mu}\rightarrow x^{\mu}+\xi^{\mu}\), one has the following transformation of modes: \[\begin{align} \phi&\rightarrow \phi-\dot{\xi}_0 & B&\rightarrow B-\xi_0-\dot{\xi} & S_i&\rightarrow S_i-\dot{\xi}_i^{T} &h^{TT}_{ij}&\rightarrow h^{TT}_{ij}\\ \psi&\rightarrow\psi & E&\rightarrow E-\zeta & F_{i}&\rightarrow F_{i}-\xi_i^{T} \end{align}\] where \(\xi_\mu=(\xi_0,\xi_i^{T}+\zeta,_i)\). We can then combine the different modes to form the following gauge-invariant modes: \[\begin{align} \begin{aligned} \Phi&=\phi+\ddot{E}-\dot{B} \\V_{i}&=S_{i}-\dot{F_i} \end{aligned} \end{align}\] where \[\begin{align} S_{i,i}&=F_{i,i}=h^{TT}_{ii}=h^{TT}_{ij,i}=0 \end{align}\] The pure gravity action, to quadratic order, is then given by \[\begin{align} \mathcal{L}^{(g)}=\mathcal{L}_{tensor}+\mathcal{L}_{vector} +\mathcal{L}_{scalar} \end{align}\] with: \[\label{eq:tensor} \mathcal{L}_{tensor}=\frac{\alpha}{4} h^{ij{TT}}\Box^2h_{ij}^{TT}\tag{12}\] \[\label{eq:vector} \mathcal{L}_{vector}=-\frac{\alpha}{2}V_i\Box \Delta V^{i}\tag{13}\] \[\begin{align} \begin{aligned} \label{eq:scalar} \mathcal{L}_{scalar}&=(12\alpha-36\beta)(\partial^2_0\psi)^2+(2\alpha-4\beta)(\triangle\Phi)^2+(6\alpha-16\beta)(\triangle\psi)^2+ (-16\alpha+48\beta)(\triangle\psi\partial_0^2\psi)\\&+(8\alpha-24\beta)(\triangle\Phi\partial_0^2\psi)+(-4\alpha+16\beta)\triangle\Phi\triangle\psi \end{aligned} \end{align}\tag{14}\] up to integration by parts.Taking the Newtonian gauge amounts to setting \(\Phi=\phi, V_i=S_i\).

Clearly, the tensor sector is positive-semi definite, while the vector one is not. To check for the scalar sector, we pass to Fourier space (\(\partial_0^2\to -\omega^2\), \(\triangle\to -k^2\), \(p^2\equiv -\omega^2+k^2\)) and write 14 as a quadratic form in the column vector \((\psi,\Phi)^T\): \[\begin{align} \label{eq:scalarkernel} \mathcal{L}_{scalar} = \begin{pmatrix}\psi \\ \Phi\end{pmatrix}^T K \begin{pmatrix}\psi \\ \Phi\end{pmatrix}, \qquad K = \begin{bmatrix} K_{\psi\psi} & K_{\psi\Phi} \\ K_{\Phi\psi} & K_{\Phi\Phi}\end{bmatrix} \end{align}\tag{15}\] where, reading off the coefficients from 14 , \[\begin{align} \begin{aligned} K_{\psi\psi} &= (12\alpha-36\beta)\omega^4+(-16\alpha+48\beta)k^2\omega^2+(6\alpha-16\beta)k^4 \\ K_{\Phi\Phi} &= (2\alpha-4\beta)k^4 \\ K_{\psi\Phi} = K_{\Phi\psi} &= \tfrac{1}{2}\left[(8\alpha-24\beta)k^2\omega^2+(-4\alpha+16\beta)k^4\right]= (4\alpha-12\beta)k^2\omega^2+(-2\alpha+8\beta)k^4 \end{aligned} \end{align}\] The two-point function matrix \(\mathcal{M}=K^{-1}\) is then obtained by the standard \(2\times2\) inversion, with \(\det K = 8\alpha(\alpha-3\beta)k^4 p^4\): \[\begin{align} \label{eq:twopointscalar} \mathcal{M}= \frac{1}{8\alpha(\alpha-3\beta)k^4p^4} \begin{bmatrix} (2\alpha-4\beta)k^4 & -k^4(-2\alpha+8\beta)-k^2\omega^2(4\alpha-12\beta)\\ -k^4(-2\alpha+8\beta)-k^2\omega^2(4\alpha-12\beta)&k^4(6\alpha-16\beta)+k^2\omega^2(-16\alpha+48\beta)+\omega^4(12\alpha-36\beta) \end{bmatrix} \end{align}\tag{16}\] As long as \(\mathcal{M}\) and by extension its inverse is PSD for all Fourier modes, the scalar action is PSD. We can easily check if this matrix is PSD: if the determinant and a given diagonal element are positive for this matrix, then this matrix is PSD in Fourier space. The determinant of this matrix is given by \[\begin{align} \begin{aligned} Det[\mathcal{M}]&=\frac{1}{8k^4p^4\alpha(\alpha-3\beta)} \\ M_{00}&=(2\alpha-4\beta)k^4 \end{aligned} \end{align}\] Then, the condition for positive semi-definiteness is \[\begin{align} \begin{aligned} 3\beta &\leq \alpha \end{aligned} \end{align}\] where \(\alpha\) is a positive constant, and \(\beta\) may take negative values.

In order to obtain the effective scalar action, \(\mathcal{L}_{eff}\), which determines the dynamics of the scalar system, we vary the action with respect to \(\Phi\) to obtain the following constraint: \[\begin{align} \label{eq:constraint} 2(2\alpha-4\beta)\triangle^2\Phi+(8\alpha-24\beta)\triangle\partial_0^2\psi+(-4\alpha+16\beta)\triangle^2\psi=0 \end{align}\tag{17}\] which is a modified form of the Hamiltonian constraint [@arnowitt1959dynamical], and may be inserted back in the action to obtain \[\begin{align} \label{eq:MPP}\mathcal{L}_{eff}=\frac{4\alpha(\alpha-3\beta)\psi\Box^2\psi}{\alpha-2\beta}. \end{align}\tag{18}\] This is equivalent to integrating out \(\Phi\), and has the same form as the action for the tensor modes. We will discuss this procedure in further detail in Section 4.2.

We thus, see that the scalar and tensor sectors are PSD, while the vector sector is indefinite. We will argue that this is benign. In addition, we can now move towards identifying the true degrees of freedom which we will see correspond to \(\psi\) and \(h^{TT}_{ij}\).

4.1 The Degrees of Freedom↩︎

Before discussing gravitational degrees of freedom in our stochastic path integral, it is of benefit to recall the method of identifying degrees of freedom in quantum (Einstein-Hilbert and quadratic) gravity, so as to better illustrate procedural differences in the two cases and their physical nature. In the case of Einstein-Hilbert gravity, we have the action (up to constants) \[\begin{align} \label{eq:Einstein-Hilbert} \mathcal{S}=\int d^4x\sqrt{-g}R=\int d^4x \mathcal{L}^{{EH}}_{scalar}+\mathcal{L}^{{EH}}_{vector}+\mathcal{L}^{{EH}}_{tensor} \end{align}\tag{19}\] where in particular,\[\begin{align} \mathcal{L}^{{EH}}_{scalar}=(2\Phi\triangle\psi+3\ddot{\psi}\psi-\psi\triangle\psi) \end{align}\] \[\begin{align} \mathcal{L}^{{EH}}_{vector}=-\frac{1}{2}V_i\triangle V^i \end{align}\] \[\begin{align} \mathcal{L}^{{EH}}_{tensor}=\frac{1}{4}h_{ij}^{TT}\Box h^{ijTT} \end{align}\] We then vary the scalar action with respect to \(\Phi\) and substitute the constraint, \(\triangle\psi=0\), back in the action. From this we have that the entire scalar sector is constrained. A similar analysis of the vector and tensor sectors confirm that the only degrees of freedom are those of the tensorial sector. Thus, the procedure is to identify the constrained variable, and use the constraint to eliminate spurious modes. Similarly, in quadratic gravity, the constrained modes (which do not have time-derivatives acting on them, and so can be safely integrated out without adding non-local in time terms to the path integral) are integrated out to identify the actual degrees of freedom.

The procedure for our classical stochastic theory is slightly different. Here the ‘action’, 3 is not the actual Lagrangian of the gravitational sector; rather, it is the Onsager-Machlup (OM) action, representing diffusion around the gravitational equations of motion. Thus, identifying a degree of freedom is no longer equivalent to integrating out constrained modes and then counting the number of \(\partial_0^2\) operators in the effective OM action; instead from the effective OM action we obtain the stochastic differential equations, the squares of which constitute our effective OM action. It is those equations which determine the degrees of freedom. Thus, an equation of motion (EOM) obtained from the OM action, implies the diffusion of a degree of freedom around its EOM, while a constraint equation implies that the constrained mode diffuses around its constraint equation, and is not an independent degree of freedom.

Let us now apply the above analysis to our effective action. From 18 we see that we have a single degree of freedom in the scalar sector, which diffuses around the wave equation. This is the stochastic Klein Gordon equation[@oppenheim2025diffusion]
[@balan2010stochastic], which we can re-write as two first order SDEs, see 7.3, i.e., \[\begin{align} \label{SDE:scalar} \Box\psi=\xi \end{align}\tag{20}\] where \(\xi\) is a mean-zero, random noise process, with covariance given by \[\begin{align} \langle\xi\xi\rangle\propto\frac{\alpha-2\beta}{4\alpha(\alpha-3\beta)} \end{align}\] and similarly for the two tensor modes: \[\begin{align} \label{SDE:tensor} \Box h_{ij}^{TT}=\xi_{ij}^{'TT} \end{align}\tag{21}\] where again \(\xi_{ij}^{'TT}\) is another mean-zero, random noise process, with covariance given by \[\begin{align} \langle\xi^{ij'TT}\xi^{'kl TT}\rangle\propto\frac{1}{\alpha}P^{ij,kl} \end{align}\] where \[\begin{align} \begin{aligned} \label{eq:projector3D} P^{ij,kl}=&\frac{1}{2}(P^{ik}P^{jl}+P^{il}P^{jk})-\frac{1}{2}P^{ij}P^{kl}\\P^{ij}=&\delta^{ij}-\frac{k^{i}k^{j}}{k^2} \end{aligned} \end{align}\tag{22}\]

On the other hand, we will now see that the vectorial action does not correspond to a degree of freedom, but to the momentum constraint[@arnowitt2008republication]. The vectorial action 13 may be written as \[\begin{align} \mathcal{L}=-\frac{\alpha}{2}(\partial_\mu \partial_j V^{i})^2 \end{align}\] This is the sector which is not PSD. But crucially, it is also has fewer time derivatives, indicating that it corresponds to a constraint rather than a dynamical degree of freedom which must be summed over in the path integral.

One may be tempted to re-write as a diffusion equation of the form \[\begin{align} \label{eq:vectorsde} \partial_\mu \partial_j V^{i}&=\xi_{\mu ,j}^{i''} \end{align}\tag{23}\] or equivalently, \[\begin{align} \tag{24} \partial_{0}\partial_jV^{i}&=\xi_{0,j}^{i''}\\ \tag{25}\partial_{k}\partial_jV^{i}&=\xi_{k ,j}^{i''} \end{align}\] where \(\xi_{\mu ,i}^{i''}=0\) and \(\xi_{k ,j}^{i''}\) is symmetric in \(k \leftrightarrow j\), with a consistency requirement, \[\begin{align} \label{eq:consistency} \partial_k\xi_{0,j}^{i''}=\partial_0\xi_{k,j}^{i''} \end{align}\tag{26}\] with a non-definite ‘covariance’ matrix of the form \[\begin{align} \label{eq:vectorcovariance} \langle\xi^{''ij}\xi^{''kl }\rangle ^{\mu\nu}\propto-\frac{1}{\alpha}\eta^{\mu\nu}(\delta^{ik}\delta^{jl}-\frac{1}{3}\delta^{ij}\delta^{kl}) \end{align}\tag{27}\] One can thus see from whence the condition for positivity, \(k^2\leq0\) arises: the ‘covariance’ matrix above is indefinite on all such modes.

Additionally, one may note that@eq:eq:vectorsde2 is an equation of constraint in the presence of random noise, while 24 governs the evolution of that constraint through the consistency relations, 26 . Alternatively, working in phase space with \(\partial_0V^{i}:=\pi^{i}\) we note that the momentum, \(\pi_i\) diffuses around the noise process, \(\xi_{0,j}^{i''}\), but the noise process does not determine its evolution in time; thus, 24 is a constraint on the momentum. In the absence of a noise process, the constraint \(\partial_k\partial_jV^i=0\) is identically fulfilled, and 24 is redundant (assuming that \(V_i\) vanishes at the spatial boundary). The action 13 then, does not represent the diffusion of an actual degree of freedom as the actions in 14 , 12 do for the scalar and tensor modes. Despite appearances, the bona fide degrees of freedom are only the scalar and tensor ones. We will discuss this in further detail Section 4.2.

One may once more compare this with the degrees of freedom in classical Einstein-Hilbert gravity. There we only have two degrees of freedom, whereas we now have diffusion around an additional scalar degree of freedom, while the vector potential is fixed by continuity. One can see why instead of diffusion around only the tensorial degrees of freedom, we now have an additional, diffusive degree of freedom by looking at the tensorial two point function, \[\begin{align} \langle h^{TTij}h^{TTkl}\rangle\propto P^{ij,kl}/\Box^4=lim_{m\rightarrow0}\frac{P^{ij,kl}}{m^2}(\frac{1}{\Box^2-m^2}-\frac{1}{\Box^2}). \end{align}\] We can always write this or the corresponding action as that of two helicity-2 ‘modes’, one massless, and the other with vanishing mass. However, a massive helicity-2 mode must, by Lorentz covariance, have corresponding helicity-0, and helicity-1 modes, which do not disappear, even in the vanishing mass case. Therefore, one necessarily gains additional time-derivatives in the vector and scalar action for free when considering helicity-2 modes diffusing around their usual equation of motion.

4.2 Non-dynamical modes and positive semi-definiteness↩︎

The OM action determines the status of each SVT field through the number of its time derivatives: four (\(\psi\), \(h^{TT}_{ij}\)) correspond to dynamical stochastic differential equations, zero (\(\Phi\)) to an algebraic constraint, and two (\(V_i\)) to a first-order constraint, i.e. initial data on phase space. In this section, we explore the consequences of this by showing that the path integral for fields with two or zero time derivatives factorises over time slices, failing to enforce temporal continuity unless the field is fixed as initial data. As a result, only \(\psi\) and \(h^{TT}_{ij}\) are integrated over in the path integral; \(\Phi\) and \(V_i\) are specified by two constraints. The indefiniteness of the vector action is therefore benign.

While in quantum gravity, the unstable mode is associated with the Hamiltonian constraint, here, it is the vector mode whose action Eq 13 is not PSD. In particular it would appear we need to restrict the path integral for the field \(V^i\) to be non-tachyonic since the action does converge if \(\omega^2\geq k^2\). We will see however, that the action for the vector sector is only indefinite if we sum over both continuous and non-continuous geometries. If instead, the path integral consists of a sum over continuous geometries, the sector is no longer indefinite.

The first thing to note, is that the action for \(V_i\) has only second time-derivatives and corresponds to a pure constraint equation. While second time derivatives correspond to dynamical degrees of freedom for a quantum theory, it does not for an action of OM form, which is an equation of motion squared. Indeed, the action Eq. 18 corresponding to \(\psi\) is proportional to \(\int d^4x(\Box\psi)^2\) term in comparison to Eq 13 for the vector mode. The fact that \(\psi\) becomes dynamical is the result of two conditions, the first being that Hamiltonian constraint is no longer required to vanish, and the second being that the dynamical part of the action has second time derivatives of \(\psi\) as can be seen from Eq 14 . This is not the case for the action of \(V^i\). It would thus appear that the vector mode remains a constraint. That the vector sector represents equations of constraints@eq:eq:vectorsde , and not degrees of freedom then implies that the indefiniteness of this sector is relatively benign.

In order to see this, recollect that the Wiener path integral for simple Brownian motion, given by the SDE \(\dot{q}=\xi\) is over continuous, non-differentiable paths. Similarly, we argue that the gravitational path integral must be over continuous, non-differentiable paths. To be more precise, we note that although we have set \(R_{\mu\nu},R\) to be random variables, they describe the evolution of geometry. Now, although the spatio-temporal evolution of geometry may well be stochastic, we ask that the geometries evolved by \(R_{\mu\nu},R\) be spatio-temporally continuous. Since \(R_{\mu\nu},R\) are second order in the spatio-temporal derivatives of our gauge-invariant potentials, this requirement translates into continuity of potentials and those first order derivatives that are evolved by these two random variables.

However, one easily sees that for equations of constraint, the naive probability transition amplitude for a mode constrained in configuration space, \(q_c\), no longer enforces the continuity condition, \(\lim_{\delta t\rightarrow0}q_c(t+\delta t)-q_c(t)=0\) irrespective of whether the action is PSD or not. A simple example suffices to illustrate our argument. Starting with a linear equation of constraint \[\begin{align} \begin{aligned} q&=\xi \end{aligned} \end{align}\] we have the following time evolution for \(\rho(q)\), the classical probability density of the system: \[\begin{align} \begin{aligned} \label{eq:ToyPI} \rho(q_f)&=\int^{\infty}_{-\infty} \mathcal{D}q\exp{-\int dt\mathcal{L}_{c}(q)}\rho(q_i)\\&=\int^{\infty}_{-\infty}\Pi dq_k\exp{-\sum_k\epsilon \mathcal{L}_{c}(q_k)}\rho(q_i)\\&=\Pi\int^{\infty}_{-\infty}dq_k\exp{-\sum_k\epsilon \mathcal{L}_{c}(q_k)}\rho(q_i) \end{aligned} \end{align}\tag{28}\] where \(\mathcal{L}_c\) is the OM action for the constrained variable, \(q_i\) and \(q_f\) indicate its initial and final positions, and \(q_k\) being the value of \(q\) at each discretised time step \(k\). The path integral factorizes over \(q_k\) and does not propagate it. The exponential just serves as a multiplier at each time step, and is normalised away. One may argue that in the case of a force, whether the path integral factorizes or not depends on the discretization convention used. While this could be true, even if it does not factorize, we still have \(q_{k+1}\neq q_{k}\) and so the argument stands and in particular, for \(\epsilon\rightarrow0\), \(q_{k+1}\neq q_{k}\), regardless of the sign of \(\mathcal{L}_{c}(q_k)\).

Contrast this with the dynamical case, where the action, \(S_{dyn}\), is of the form \(\mathcal{S}_{dyn}=\frac{1}{2D}\sum_k\epsilon (\frac{q_{k+1}-q_k}{\epsilon}-F(q_k))^2\). Due to the time-derivative, the action does not factorize, and in the short time limit, \(\epsilon\rightarrow0\), we do indeed obtain \(q_{k+1}=q_k\). This indicates that a constrained particle follows continuous trajectories with a vanishing probability. To enforce continuity without significantly altering the dynamics (see 7.3), we add a delta function to the action, \(\delta(q_{k+1}-q_k)\) at every time slice. Continuous trajectories then evolve deterministically.

The argument above may also be easily generalized to constraints in phase space: if we are to have continuous trajectories on phase space, the action, \(\mathcal{L}(q,p)\) must either be a function of \(\dot{p}_k\) to ensure continuity (continuity in \(q\) is ensured by adding the delta function \(\delta(p_k-\frac{q_{k+1}-q_{k}}{\delta t})\) at every time slice), or continuity must be enforced through delta functions. Actions arising from constraints may thus, only give rise to continuous trajectories in phase space if continuity is actively enforced.

As a second example, we take a stochastic form of Newton’s equation in vacuum, such as those appearing in other stochastic theories of gravity[@hu2008stochastic; @moffat1997stochastic; @tilloy2016sourcing] \[\begin{align} \begin{aligned} \triangle\phi&=\xi\\\langle\xi\xi\rangle&=D_N \end{aligned} \end{align}\] with a PSD diffusion matrix, \(D_N\). The probability distribution at the spacial slice \(\Sigma_{t_f}\) at final time, given the distribution at the initial slice, \(\Sigma_{t_i}\) then reads \[\begin{align} \begin{aligned} \rho(\phi(\Sigma_{t_f}))=&\int D\phi\exp{-\frac{1}{2D_N}\int d^4x (\triangle\phi)^2} \rho(\phi(\Sigma_{t_i})). \end{aligned} \end{align}\] It is easy to see that the path integral only sums over spatially connected fields. However, the path integral factorizes in time, and so the field is discontinuous in time. If we want to make it continuous, we would need to add a delta function at each constant time slice: \(\delta(\phi(\Sigma_{t_{k+1}})-\phi(\Sigma_{t_{k}}))\). The final probability distribution, \(\rho(\phi(\Sigma_{t_f}))\) is then simply a function of \(\phi(\Sigma_{t_{i}})\), possibly up to a normalization.

In the case considered here, we can instead determine \(\Phi\) from its equation of constraint, 17 . This is because for a PSD action, \(\mathcal{S}(\theta)\), that is at most quadratic in the field \(\theta\), the most probable path, \(\mathcal{L}_{mpp} (\theta)\) obtained from setting \(\frac{\delta\mathcal{S}(\theta)}{\delta\theta}=0\) determines the transition amplitude [@Chaichian2018PathI]. One may thus, simply vary the PSD scalar action, 14 with respect to \(\Phi\) to obtain its equation of constraint, 17 . This equation is identically satisfied for the most probable path (which in turn is the only trajectory we consider), and thus, \(\Phi\) is completely determined by it at every constant time slice, and so is made to evolve continuously.

Let us now also consider the case for dynamical modes. There, we have the generalized action for a field \(\sigma\) \[\begin{align} \mathcal{S}_{dyn,mod}=\frac{1}{2D_{gen}}\int d^4x(\Box\sigma)^2 \end{align}\] This action is continuous in both phase space for a fixed \(\nabla^2\sigma\), as well as in \(\partial_j\sigma\) for a fixed \(\pi_{\sigma}\). If however, it had an overall minus sign, it would be discontinuous in \(\textit{both}\) space and time. To ensure continuity then, it would be kept fixed in both space and time, and would not diffuse at all. In a classical-quantum theory this would then imply instantaneous decoherence for the matter field coupled solely to \(\sigma\).

One can now finally turn to the case of the vector modes. It is clear from both 13 and 27 that the vector action does not enforce continuity in space for \(\partial_jV^i\) and by implication, for \(V^i\), whereas, continuity in time is maintained. Continuity in phase space, however, is not maintained. We must thus, run the same argument as above, with the difference that we are now concerned with time-like trajectories, \(\vec{x}_{t,\vec{q}}\) where \(t\) is the time index and \(\vec{q}\) denotes spatial position, instead of constant time slices. Although these trajectories evolve stochastically, at a constant \(t\), they must be continuous in space. Then by the fact that continuity in space is enforced through a delta function, \(\delta(V^i(\vec{x}_{\tau,{\vec{s+1}}})-V^i(\vec{x}_{\tau,\vec{s}}))\), a trajectory at position \(\vec{s}\) determines trajectories at all \(\vec{q}=\vec{s+i}\). After integrating out the delta functions, we thus, obtain the following action: \[\begin{align} \mathcal{S}_{vec,con}=-\frac{\alpha}{2}\int d^4x V^i_{\vec{s}}(t)\nabla^2\Box V^i_{\vec{s}}(t) \end{align}\] with the measure, \(\mathcal{D} V=\prod_{t} d V_{t,\vec{s}}\). On the other hand, if we demand continuity in phase space, \(V_i\) is completely determined in both space, and time.

5 Experimental Implications↩︎

Our purpose in this section will be to obtain upper and lower bounds on both diffusion constants \(\alpha\) and \(\beta\) in order to constrain the theory. We first compare the two point functions for the scalar and tensor sectors and then obtain bounds for the diffusion coefficients of the two sectors from LISA Pathfinder and LIGO. The two sectors have different diffusion constants, but the scalar one effectively upper bounds the tensorial one. Finally, for a point of comparison we reproduce the effective decoherence rate which we use to set a lower bound on the diffusion constant(s). In comparison to previous results, our analysis allows us to distinguish the role of both constants \(\alpha\) and \(\beta\) in determining observational constraints. Our results here should be considered as preliminary, since they rely on the linearised action, and we do not know at what scale, perturbation theory breaks down.

5.1 Comparison between the Scalar and the Tensor Sectors↩︎

How are the purely radiative scalar and tensor diffusion related, and can they be used to set a bound on each other, given that we have an additional parameter, \(\beta\) for scalar diffusion?

The correlation matrix for the scalar sector is given by \[\begin{align} \mathcal{M}=\frac{1}{8\alpha(\alpha-3\beta) k^4p^4}\begin{bmatrix} (2\alpha-4\beta)k^4 & -k^4(-2\alpha+8\beta)-k^2\omega^2(4\alpha-12\beta)\\ -k^4(-2\alpha+8\beta)-k^2\omega^2(4\alpha-12\beta)&k^4(6\alpha-16\beta)+k^2\omega^2(-16\alpha+48\beta)+\omega^4(12\alpha-36\beta) \end{bmatrix} \end{align}\] in the \(\psi,\Phi\) basis. One can decompose the \(\langle\Phi\Phi\rangle\) component into frequency-dependent, and frequency independent parts: \(\langle\Phi\Phi\rangle=\frac{3}{2\alpha k^4}+(1/\alpha p^4)(\frac{\omega^2}{k^2}+\frac{(-3\alpha+10\beta)}{4(\alpha-3\beta)})\). The frequency-independent part, \(\frac{3}{2\alpha k^4}\) is the two-point function that we would expect from diffusion around a constraint equation of the form \(\nabla^2\Phi=\xi\), \(\langle\xi(x)\xi(x')\rangle\propto \frac{2\alpha}{3}\delta^3(x-x')\). For \(\beta=(3/10)\alpha\) the purely radiative part (\(\frac{(-3\alpha+10\beta)}{4\alpha(\alpha-3\beta)p^4}\)) disappears, but a frequency dependent part still survives.

From 14 one notes that in the limit, \((\alpha-3\beta)\rightarrow0\),14 is singular. One can also see this directly from the two point function 16 by expanding \(\mathcal{M}\) in \(\alpha-3\beta\), \[\begin{align} \lim_{(\alpha-3\beta)\rightarrow0^+}\mathcal{M}=\frac{1}{8\alpha p^4}\begin{bmatrix} (2\beta)/(\alpha-3\beta) & (-2\beta)/(\alpha-3\beta)\\ (-2\beta)/(\alpha-3\beta)&(2\beta)/(\alpha-3\beta) \end{bmatrix}+O((\alpha-3\beta)^0) \end{align}\] The other limit is where \(\beta \rightarrow-\infty\): \[\begin{align} \lim_{\beta\rightarrow-\infty}\mathcal{M}=\frac{1}{3\alpha k^4p^4}\begin{bmatrix} 4k^4 & 8k^4-12k^2\omega^2\\ 8k^4-12k^2\omega^2&16k^4-48k^2\omega^2+36\omega^4 \end{bmatrix} \end{align}\] Clearly, for most of its allowed parameter space, the amount of radiative scalar diffusion is directly proportional to that of tensorial diffusion, before diverging as \((\alpha-3\beta)\rightarrow0\). Thus, radiative scalar diffusion may be used to set an upper bound on the amount of tensorial diffusion. Experiments based on the amount of tensorial diffusion, on the other hand, may only be used to set an upper bound on the value of the diffusion constant, \(\alpha^{-1}\), which is relevant only for those frequency branches in the scalar sector that are independent of \((\alpha-3\beta)\). In conclusion, it is to experiments based on tensorial diffusion that we must turn to obtain the most stringent upper bounds on noise.

Alternatively, one may consider the situation in terms of diffusion around the trace, and the traceless parts of the Einstein tensor. From \(\mathcal{M}\), one notes that as \(\alpha\rightarrow3\beta\) the dynamical part of \(\langle\Phi\Phi\rangle\) dominates, while for large, negative values of \(\beta\), the non-dynamical part is of increased importance, but the dynamical part does not disappear. Note from the action \(\int d^4x (\alpha G^{\mu\nu}-\beta G^2)\) that the latter is the regime in which the trace part of the Einstein equations is highly constrained. One can also compare our results with the regime where diffusion takes place around the trace, \(G\), and the action has the form, \(-\int d^4x\beta G^2\). From the positivity condition, only negative values of \(\beta\) are permissible. In this case, the correlation matrix does not exist, as the matrix of the quadratic action, \(\mathcal{M}^{-1}\) is singular, indicating a breakdown in the theory, as we would expect from taking the \(\alpha \rightarrow0\) limit in \(\mathcal{M}\). We can thus identify two regimes in which the theory is singular; when \(\alpha\rightarrow0\) and when \(\alpha\rightarrow3\beta^+\). In the former regime, the tensorial action will also be singular, and the theory as a whole breaks down. The latter limit corresponds to the action \(\int d^4x C^2\) where \(C\) is the Weyl-tensor; the scalar action now has an additional symmetry (conformal symmetry) which renders it degenerate.

5.2 Spectral Acceleration↩︎

We can now use the two-point function for the scalar sector to find the variance in acceleration for a non-relativistic particle. We may then bound the diffusion coefficient using precision gravity measurements from Lisa Pathfinder. The significant contribution to the variance in acceleration, \(\langle a_ia^i\rangle\) in the non-relativistic limit is then given by the spectral acceleration of \(\langle\Phi\Phi\rangle\), i.e., \(\langle\partial_i\Phi(x)\partial^i\Phi(y)\rangle\), where we note that in the absence of matter, \(\langle a_i \rangle=0\) for the linearised action. We first calculate the partial Fourier transform of the two point function (in space only), assuming that it is the convolution of two retarded Green’s functions (see 7.4 and [@oppenheim2024anomalous; @oppenheim2025diffusion]): \[\begin{align} \langle\Phi \Phi \rangle=\int ^{\infty}_{-\infty}d^3z\int^{t_f}_0dz^0G_{ret}(x-z)G_{ret}(y-z)=\int ^{\infty}_{-\infty}d^3z\int^{t_f}_0dz^0\int\frac{d^4p}{(2\pi)^4}\int\frac{d^4p'}{(2\pi)^4}G(p)G(p') \end{align}\] where \[\begin{align} G(p)G(p')=\frac{1}{8\alpha(\alpha-3\beta)} (\frac{8(\alpha-3\beta)}{k^2k'^2}+\frac{4\omega^2\omega'^2(\alpha-3\beta)}{k^2k'^2p^2p'^2}+\frac{(-2\alpha+8\beta)}{p^2p'^2}) \end{align}\]

Using \[\begin{align} \int^\infty_0e^{-iz^0(\omega+\omega')}dz^0=\pi\delta(\omega+\omega')-iP.V(1/(\omega+\omega')) \end{align}\] we obtain, \[\begin{align} \int ^{\infty}_{-\infty}d^3z\int^{t_f}_0dz^0\int\frac{dp}{(2\pi)^4}\int\frac{dp'}{(2\pi)^4}G(p)G(p')=\mathcal{G}_{stationary}+\mathcal{G}_{non-stationary} \end{align}\] where we have split the integral in a part that solely depends on \(x^0-y^0\), and a part that depends on both \(x^0,y^0\). Then, \[\begin{align} \mathcal{G}_{stationary}=\int\frac{d^4p}{(2\pi)^4}\frac{1}{8\alpha(\alpha-3\beta)} (8(\alpha-3\beta)/k^4+4\omega^4(\alpha-3\beta)/k^4p^4+(-2\alpha+8\beta)/p^4) \end{align}\]

Meanwhile, for \(\mathcal{G}_{non-stationary}\) we need only consider the non-stationary parts of terms of the form \(\frac{\omega^2\omega'^2}{k^2k'^2p^2p'^2},\frac{1}{p^2p'^2}\). For details, see Appendix 7.5.

We obtain the following time-dependent parts, in frequency-space \[\begin{align} \lim_{\epsilon\rightarrow0} \mathcal{FT}_{\omega}(1/p^4)|=-\frac{1}{8\pi^2|x-y|}(\frac{\sin{\omega|x-y|}}{\omega}(\frac{1}{\epsilon}-(x^0+y^0))+\frac{\pi|x-y|(x^0+y^0) }{|x^0+y^0-|x-y||+|x^0+y^0+|x-y||}2\pi\delta(\omega)) \end{align}\] where by \(\mathcal{FT}_{\omega}\) we refer to the full Fourier transform, followed by an inverse Fourier transform with respect to \(x^0-y^0\). This is equivalent to the spatial Fourier transform only for the stationary part, where we define \[\begin{align} \delta(\omega)=\frac{1}{2\pi}\int^\infty_{-\infty}d(x^0-y^0)e^{i\omega(x^0-y^0)} \end{align}\] and have only kept terms up to \(\epsilon^0\). We have a similar expression from the time-dependent term of the form \(\omega^4/k^4p^4\). For the stationary part, we have for \[\begin{align} \begin{aligned} \mathcal{FT}_{\omega}\frac{\omega^4}{k^4p^4}(x,\omega)&=\frac{\omega^4}{8\pi x(\omega^2+\epsilon^2)^4}((-x^2\epsilon^4+4\omega^2-x^2\omega^4-2\epsilon^2(2+x^2\omega^2))\\&+\frac{e^{-\epsilon x}}{\epsilon\omega}(\sin{\omega x}(\epsilon^4-6\epsilon^2\omega^2+\omega^4)+\epsilon\cos(\omega x)(4\epsilon^2\omega-4\omega^3))) \end{aligned} \end{align}\]

and \[\begin{align} \lim_{\epsilon\rightarrow0}\mathcal{FT}_{\omega}\frac{\omega^4}{k^4p^4}=&\frac{1}{8\pi x}(-x^2+\frac{4}{\omega^2}+(\frac{\sin{\omega x}}{\omega}(\frac{1}{\epsilon}-x)-\frac{4\cos(\omega x)}{\omega^2})) \end{align}\] \[\begin{align} \lim_{\epsilon\rightarrow0}\mathcal{FT}_{\omega}\frac{1}{p^4}=\frac{\sin{\omega|x|}}{8\pi\omega|x|}(1/\epsilon-|x|) \end{align}\]One thus obtains \[\begin{align} \label{eq:2pointNewton} \begin{aligned} \langle \Phi \Phi\rangle(|x-y|,x^0+y^0,\omega)&=\frac{1}{64(\alpha-3\beta)\alpha\pi}(-12\eta|x-y|+\frac{16(\alpha-3\beta)}{\omega^2|x-y|} +\frac{(2\alpha- 4\beta)\sin{\omega |x-y|}}{\omega|x-y|}(x^0+y^0-|x-y|)\\&-\frac{\cos{\omega |x-y|}}{\omega^2|x-y|}(16(\alpha-3\beta))-\frac{2\pi\delta(\omega)(14\alpha-40\beta)(x^0+y^0)}{(|x^0+y^0-|x-y||+|x^0+y^0+|x-y||)}) \end{aligned} \end{align}\tag{29}\] For the spectral density, we have \[\begin{align} \label{eq:2pointNewtonPSD} \begin{aligned} \langle\partial_i \Phi \partial^i\Phi\rangle&=-\nabla^2\langle \Phi \Phi\rangle \\&=\frac{1}{64(\alpha-3\beta)\alpha\pi}(\frac{24(\alpha-3\beta)}{|x-y|}+(2\alpha-4\beta)(\frac{\omega\sin{\omega|x-y|}}{|x-y|}(x^0+y^0-|x-y|))+\frac{(40\beta-12\alpha)\cos{\omega |x-y|}}{|x-y|}) \end{aligned} \end{align}\tag{30}\] where the second step follows by the linearity of the expectation value. One can integrate this over test masses with uniform density, \(\rho(x)=\frac{M}{V}\) to compute the covariance of the force exerted due to the Newtonian potential: \[\begin{align} \boldsymbol{E}(F_iF^i)(\omega)= \int_0^{R_2} d^3x\int_0^{R_1}d^3y\rho(\vec{x})\rho(\vec{y})\boldsymbol{E} (\partial_i \Phi(\vec{x}) \partial^i\Phi(\vec{y})) \end{align}\]

Matters may be simplified considerably by taking the low-frequency limit of \(-\nabla^2\langle \Phi\Phi\rangle\) before integrating it over a test mass: \[\begin{align} \lim_{\omega\rightarrow0}-\nabla^2\langle \Phi \Phi\rangle=\frac{1}{64(\alpha-3\beta)\alpha\pi}(\frac{4(3\alpha-8\beta)}{|x-y|}+ \omega^2(8(\alpha-3\beta)|x-y|+(2\alpha-4\beta)(x^0+y^0))+O(\omega^3)) \end{align}\] Then, we obtain, \[\begin{align} \boldsymbol{E}(F_iF^i)(\omega)=\frac{1}{64(\alpha-3\beta)\alpha\pi}(\frac{24(3\alpha-8\beta)M^2}{5R}+ \omega^2(8(\alpha-3\beta)(36M^2R/35)+(2\alpha-4\beta)(2T)M^2)+O(\omega^3)) \end{align}\] where we assume spherical geometry, and \(T=:\frac{x^0+y^0}{2}\).

One can compare this result with the noise results obtained during the LISA Pathfinder mission. The purpose of the mission was to show that a test-mass could be placed in near perfect free-fall, with the noise-levels below a certain threshold. For that purpose, the acceleration of a test mass, called TM1, was computed relative to another one a few centimeters away, and the power spectral density computed. Three branches were noted in the power spectral density (psd): a low-frequency branch with an approximate \(\omega^{-2}\) behavior, a frequency independent branch above 1mHz, and a rising branch above 10 mHz [@armano2024indepth]. We have in the LPF, \(\langle\Delta g\Delta g\rangle (\omega)=\langle g_1 g_1\rangle+\langle g_2 g_2\rangle-2\langle g_1 g_2\rangle\). Taking the radii of the spheres to be \(R\), and the distance between them to be \(D\), we obtain for the variance in acceleration, \[\begin{align} \langle\Delta g\Delta g\rangle (\omega)=\frac{1}{32(\alpha-3\beta)\alpha\pi}(4(3\alpha-8\beta)(6/5R-1/(2R+D))+ 8(\alpha-3\beta)\omega^2(-34R/35-2R^2/5(2R+D)-D)+O(\omega^3)) \end{align}\] where \(D+2R=D_{total}\), the total distance between the centres of the spheres, which for the sake of simplicity, we have treated as a constant. Notably, the average total time elapsed does not appear in this expression, irrespective of the geometry considered. Taking the leading \(\omega^{0}\) component, and reintroducing factors of c, we have \[\begin{align} \langle\Delta g\Delta g\rangle (\omega)=\frac{(3\alpha-8\beta)c^3}{8\pi\alpha(\alpha-3\beta)}(\frac{6}{5R}-\frac{1}{2R+D}) +O(\omega). \end{align}\].The upper bound for noise at \(\omega=10^{-4}\)Hz is \(10^{-30}\text{m}^2\text{s}^{-4}\text{Hz}^{-1}\). \(R=46\text{mm}, 2R+D=30\text{cm}\)[@armano2024indepth] giving us \[\begin{align} \label{eq:NewtonbOUND} \frac{(3\alpha-8\beta)}{\alpha(\alpha-3\beta)}:=D_{2, scalar}\lessapprox 10^{-56} \end{align}\tag{31}\]

5.3 Waves in FLRW↩︎

We now upper bound the tensorial diffusion constant using bounds from LIGO. One may recall that around Minkowski space the two point functions for the scalar and tensor functions are non-stationary (even in the wider sense as they do not depend solely on \((t_2-t_1)\), and indeed diverge if we integrate from \(t=-\infty\). However, the addition of a finite amount of friction is sufficient to make the two-point functions stationary. We may thus consider the MSR type path integral (see Appendix 7 ) for the tensorial sector around an FLRW background. Since the tensorial (or transverse traceless) components of the auxiliary field, \(\tilde{h}\) are gauge-invariant, we can identify the correct auxiliary variable without trouble. The tensorial action is then given by \[S=\int d^4x\sqrt{-g}(\tilde{h}_{ij}^{TT}\Box'{h}_{ij}^{TT}-\frac{\alpha}{4} (\tilde{h}_{ij}^{TT})^2)\] where \[\begin{align} \Box'=\partial_t^2+2\eta\partial_t-\triangle+2K \end{align}\] where \(\eta=aH_0\), \(K\) is the scalar, 3-dimensional curvature, and \(H_0\) is the Hubble constant. Assuming spatial-flatness, we set \(K=0\). The two point function, in Fourier space is then given by \[\begin{align} \label{HTT} \langle h^{ij,TT}h^{kl,TT}\rangle(p)=\frac{1}{(2\pi)^4}\frac{\alpha^{-1} P^{ij,kl}}{(k^2-\omega^2+2i\omega\eta)(k^2-\omega^2-2i\omega\eta)} \end{align}\tag{32}\] where \(P^{ij,kl}\) is given by 22 . It is simple enough to obtain the energy density per frequency: \[\begin{align} \begin{aligned} \rho_{GW}(\omega)&=(\partial_0h_{ij}^{TT})^2(\omega,0)+(\nabla h_{ij}^{TT})^2(\omega,0)= \lim_{(x-x')\rightarrow 0}(\partial_0^2+\nabla^2) \int \frac{d^3{k}}{(2\pi)^4}\langle h^{TT}_{ij}h^{ij,TT}\rangle(p))e^{ipx} \end{aligned} \end{align}\] The spatial Fourier transform is given by \[\begin{align} \langle h^{ij,TT}h^{TT}_{ij}\rangle(\omega,\vec{x}-\vec{x'}) \approx \frac{\alpha^{-1}e^{-\eta|x|}}{2(2\pi)^2\eta\omega x}\sin{x(\omega-\frac{\eta^2}{2\omega})}\approx \frac{e^{-\eta|x|}}{2(2\pi)^2\eta\omega x}\sin{x(\omega)} \end{align}\] where we have only kept terms up to order \(\eta^2\) in the argument of \(\sin\) and in the exponent, given that \(\eta \sim 10^{-18}\text{s}^{-1}\). The energy density per frequency in units of \(\text{kg}\text{m}^{-1}\text{s}^{-1}\) is given by \[\begin{align} \mathcal{\rho}_{GW}(x)(\omega)\approx\frac{\alpha^{-1}c^2e^{-\frac{\eta}{c}|x|}}{16\pi G_N(2\pi)^2\eta\frac{\omega}{c} x}((2\omega^2-\eta^2)\sin{x(\frac{\omega}{c})}+2\eta\omega\cos{x(\frac{\omega}{c}})) \end{align}\] where we have restored factors of \(8\pi G_N,c\). Taking the coincident limit, \(|x-x'|\rightarrow 0\) we have \[\begin{align} \label{eq:EnergyDensity} \mathcal{\rho}_{GW}(0)(\omega)\approx\frac{\alpha^{-1}c^2}{16\pi G(2\pi)^2}(\frac{((2\omega^2-\eta^2)}{\eta}+2c\delta(0)) \end{align}\tag{33}\] Keeping only the leading order terms in \(\eta\), we have the logarithmic energy density per logarithmic frequency given by \[\begin{align} \label{eq:tensorintegrated}d\text{ln}\mathcal{\rho}_{GW}(\omega)\approx\frac{\alpha^{-1}c^2}{8\pi G(2\pi)^2}(\frac{\omega^3}{\eta}+\omega c\delta(0))d\text{ln}\omega \end{align}\tag{34}\] If we regularize the delta function by assuming that the second term in 33 is of order, \(\frac{1}{\eta}\), then we obtain for the length scale in infra-red, \(l\sim 10^{26}\text{m}\), of the same order of magnitude as the radius of the visible universe. We can now use \(\Omega_{GW}(\omega)=\frac{d\ln\rho_{GW}}{\rho_c d \text{ln}\omega}=\frac{\omega\rho_{GW}}{\rho_c}\) , \(\rho_c=c^{2}\times 9.47\times 10^{-27}\text{kg}\text{m}^{-1}\text{s}^{-2}\) is the critical energy density required to close the universe, to put a bound on \(\alpha^{-1}\). From [@abac2025upper], we have \(\Omega_{GW}(25\text{Hz})\leq 10^{-9}\), from which we obtain \[\begin{align} \label{eq:tensorbound} \alpha^{-1}:=D_{2,tensor}\lessapprox10^{-66} \end{align}\tag{35}\] 10 orders of magnitude below \(D_{2,scalar}\) obtained from the noise threshold in LISA Pathfinder. A consistent estimate, \(D\lessapprox 10^{-68}\) (see also [@penington2024]), was obtained in [@oppenheim2025diffusion] using the stochastic Klein-Gordon equation as a toy model for the tensor sector, and a similar bound appears in [@hirotani2026testing].

Coherence experiments on molecular superpositions lower-bound the dimensionless diffusion coefficient at \(D_2\gtrapprox 10^{-71}\) [@oppenheim2025diffusion], via the decoherence-diffusion tradeoff. This leaves a viable window of roughly five orders of magnitude, \(10^{-71}\lesssim D_2 \lesssim 10^{-66}\) although further theoretical work on decoherence bounds could narrow this window.

If one instead assume the spectrum holds throughout, and use the integrated energy density \(\rho_{GW}\propto D_2\Lambda^3 T\) (as in [@oppenheim2025diffusion]), the resulting bound requires the UV cutoff to satisfy \(\ell\gtrsim 10^{2}\mathrm{m}\), which is far above the millimeter scale to which gravity has been probed. With scale-independent couplings, this would rule out the theory. However, asymptotic freedom [@grudka2024renormalisation] causes the coupling to run, and matter coupling plus non-linear effects and cosmological considerations may further suppress the integrated energy density; a definitive assessment requires a full treatment on an expanding background.

5.4 The effective decoherence rate↩︎

Given the dependence of the diffusion and decoherence terms on two parameters, \(\alpha,\beta\), does the diffusion-decoherence tradeoff hold for static potentials, if we integrate out \(\psi\)? More generally, in decoherence experiments, the measurement is of the decoherence rates, and not of the gravitational fields. Thus, to obtain the correct decoherence rate, we must integrate out all mediating gravitational fields. Such a calculation is beyond the scope of this work, but we include the analysis of [@oppenheim2026non] to better understand the role of the various coupling constants.

As an example, let us take 1 , with \(T^{\mu\nu}=\rho(x)U^{\mu}U^{\nu}\), i.e., pressure-less dust with energy density \(\rho(x)\). The trace, \(T\) is given by \(-\rho(x)\), and \(T_{\mu\nu}T^{\mu\nu}=\rho(x)^2\). Working in the rest-frame of this distribution, in units of \(8\pi G_N,\) we obtain the following Lagrangian density: \[\begin{align} \mathcal{S}_{dust}=S^{(g)}+\int d^4x\frac{(\alpha-\beta)}{8}(4\bar{\rho}^2+(\Delta\rho)^2)+{\frac{1}{2}}((-4\alpha+8\beta)\bar{\rho}\triangle\psi-4\beta\bar{\rho}\triangle\Phi)-i\Delta\rho(1-\phi+\nabla^2E+3\psi) \end{align}\] where \(\mathcal{S}^{(g)}\) given by 3 , is the pure gravitational action and \(\bar{\rho}=:(\rho^{+}+\rho^{-})/2\), and \(\Delta\rho=:\rho^{+}-\rho^{-}\). Since the vector and tensor fields do not interact with the dust distribution in its rest frame, we will ignore them forthwith, and only focus on the scalar sector. We then have the following Lagrangian density: \[\begin{align} \label{eq:DustAction} \begin{aligned} \mathcal{L}_{dust}&=\frac{1}{2}(\alpha((2\nabla^2\psi-\bar\rho)^2+2(\nabla^2(\psi-\Phi))^2+12(\partial_0\psi)^2-16\partial_0^2\psi\nabla^2\psi+8\partial_0^2\psi\nabla^2\Phi)-\beta(2\nabla^2(2\psi-\Phi)-6\partial_0^2\psi-\bar\rho)^2)\\&+\frac{\alpha-\beta}{8}(\Delta\rho)^2-i\Delta\rho(1-\phi+\nabla^2E+3\psi) \end{aligned} \end{align}\tag{36}\] In the case of slowly changing \(\psi,\Phi\), \(\partial_t \ll \partial_i\) the terms with time-derivatives drop out, allowing us to integrate out \(\psi\) without adding non-local in time terms to the action. We obtain the following, non-local in space effective action: \[\begin{align} \begin{aligned} \mathcal{S}_{dust,stat,eff}[\Phi]&=\int d^4x[(\frac{\alpha(\alpha-3\beta)}{2(3\alpha-8\beta)}(2\triangle\Phi-8\pi G_N\bar{\rho})^2+(\frac{64\pi^2G_N^2(\alpha-\beta)}{8})(\Delta{\rho})^2 \\&-i\Delta\rho[1-\phi+\nabla^2E-\frac{3}{(3\alpha-8\beta)}((-\alpha+4\beta)\Phi)])(x)\\&+\int d^3y(\frac{9}{4(3\alpha-8\beta)\nabla^4}(\Delta{\rho})^2 +i\Delta\rho\frac{3(8\pi G_N)}{(3\alpha-8\beta)}(\frac{(-\alpha+2\beta)\bar\rho}{\nabla^2}))(x,y)] \end{aligned} \end{align}\] where \(\frac{1}{\nabla^2}\) refers to the spatial Fourier transform of \(\frac{1}{k^2}\), and where we have explicitly written out the argument of the terms involved.

This has an obvious interpretation in the Newtonian gauge. The Newtonian potential \(\Phi\) is sourced by the average of the bra and ket dust fields but diffuses around the solution to Poisson’s equation, with an effective diffusion constant, \(D_{diff,eff}=\frac{(3\alpha-8\beta)}{\alpha(\alpha-3\beta)}\). This is the effective diffusion constant, \(D_{2,scalar}\) that was obtained in the \(\omega\rightarrow0\) limit in 31 , as one would indeed expect. Compare this with the approach of [@layton2023weak]. Therein, one assumed for the scalar sector that only the \(G_{00}\) component of Einstein’s equations were diffusive. That assumption would have allowed us to set \(\psi=\Phi\), and we would have consequently obtained an effective diffusion constant, \(D'_{diff,eff}=\frac{1}{(\alpha-\beta)}\).

In our case, the rate of decoherence, \(D_{dec,eff}=8\pi^2G_N^2(\alpha-\beta)+\frac{9}{4k^4(3\alpha-8\beta)}\), where the second term represents the secondary decoherence, does not set an upper bound on the rate of diffusion, \(D_{diff,eff}\) which diverges as \(\beta \rightarrow \alpha/3\), for a finite decoherence rate. Similarly, \(\lim_{\beta\rightarrow-\infty} D_{diff,eff}=\frac{8}{3\alpha}\), where clearly the decoherence rate, \(D_{dec,eff}\) diverges. The two rates are essentially independent of each other, although one should note that since \(D_{diff,eff}\) bounds \(\alpha^{-1}\) and thus the minimal decoherence rate, \(D_0>\frac{16\pi^2G_N^2\alpha}{3}\) obtained from decoherence experiments can only further constrain the parameters, \(\alpha,\beta\) of the theory.

By integrating out \(\psi\), which mediates interactions between the bra and the ket fields, we also obtain terms of the form \(\rho^{+}\rho^{-}\). A pure initial state will now be sent to a mixed one, even when conditioning on a particular classical trajectory. This is a direct consequence of the fact that the mass-field interaction term is of the form, \(\bar{T}^{\mu \nu}G_{\mu\nu}\), instead of \(\bar{T}^{\mu\nu}h_{\mu\nu}\). Physically, since we have integrated out a mediating field, we would expect an initially pure state to decohere. Additionally, the unitary part now has terms of the form \((\rho^{+})^2,(\rho^{-})^2\): however, note that there are no \(\rho^{+}\rho^{-}\) mixed terms in the unitary part.

Finally, setting \(\phi=\Phi\) and thus, \(\nabla^2E=0\), or in other words, taking the Newtonian gauge, and integrating out \(\Phi\), we obtain the following action: \[\begin{align} \begin{aligned} \mathcal{S}_{dust,static}[\bar{\rho},\Delta\rho]&=-\int d^4x[((\Delta\rho)^2\frac{(64\pi^2G_N^2)(\alpha-\beta)}{8}-i\Delta\rho)(x)\\&+\int d^3y(\frac{(3\alpha-\beta)}{4\alpha(\alpha-3\beta)\nabla^4})(\Delta\rho)^2-i\Delta\rho\frac{8\pi G_N\bar\rho}{\nabla^2})(x,y)] \end{aligned} \end{align}\] After integrating out all the interacting classical fields, the real part consists solely of a decoherence term, which by our initial positivity condition, \(\alpha\geq3\beta\) is always larger than the original decoherence term (which had the coefficient \(8\pi^2 G_N^2(\alpha-\beta)\)). We have thus made the theory only more constrained than before. One may also note that the limit, \(\alpha^+\rightarrow3\beta\) leads to instantaneous decoherence.

5.4.1 Non-Markovian Effects↩︎

If we drop the assumption, \(\partial_t \ll \partial_i\), and integrate out \(\psi\) and \(\Phi\), we obtain, in the Newtonian gauge, the following, non-local action: \[\begin{align} \begin{aligned} \mathcal{S}_{dust,effective}[\Delta\rho,\bar\rho]&=\int d^4x\frac{64\pi G_N^2(\alpha-\beta)}{8}((\Delta\rho)^2+4\bar{\rho}^2)-i\Delta\rho \\&+\int d^4y\Delta\rho\frac{((3\alpha-\beta)\nabla^4+2(\alpha-3\beta\nabla^2\partial_0^2+3(\alpha-3\beta)\partial_0^4))}{4\alpha(\alpha-3\beta)(\nabla^2-\partial_0^2)^2\nabla^4}\Delta\rho\\&-64\pi^2 G_N^2\bar\rho\frac{(\alpha-3\beta)(\alpha-\beta) \nabla^4+2(\alpha-3\beta)\nabla^2\partial_0^2+3\beta^2\partial_0^4)}{2(\alpha-3\beta)((\nabla^2-\partial_0^2)^2))}\bar\rho\\&-i\Delta\rho\frac{8\pi G_N((\alpha-3\beta)\nabla^2+(\alpha-\beta)\partial_0^2)}{(\alpha-3\beta)(\nabla^2-\partial_0^2)^2}\bar\rho \end{aligned} \end{align}\] In momentum-space, one thus, obtain the non-Markovian decoherence coefficient, \(\mathcal{D}_{\mathrm{decoh}}(\omega,k) = 8\pi^2 G_N^2(\alpha-\beta)\\+\frac{(3\alpha-\beta)k^4+2(\alpha-3\beta)k^2\omega^2+3(\alpha-3\beta)\omega^4}{4\alpha k^4(\alpha-3\beta)(\omega^2-k^2)^2}\), which, like the static decoherence rate, \(\mathcal{D}_{\mathrm{decoh,stat}}(k)=8\pi^2 G_N^2(\alpha-\beta)+\frac{(3\alpha-\beta)}{4\alpha(\alpha-3\beta)k^4}\) is larger than the initial minimum decoherence rate, \(\frac{\alpha-\beta}{8}\) and therefore leads to a tighter bound on the parameters, \(\alpha,\beta\).

6 sec:Discussion↩︎

The path integral formulation of postquantum gravity is locally Lorentz invariant, since the covariance matrix (the deWitt metric) only depends on the metric \(g\). This creates a tension with the action being indefinite, since the metric’s Lorentzian signature mean that in \(3+1\) dimensions the DeWitt has negative eigenvalues. In quantum gravity the action is also indefinite, but due to the conformal mode. Regardless, we have shown here that there do exist gauge-invariant, positive-semi-definite sectors that correspond to the gravitational degrees of freedom, and that the indefinite sector is non-dynamical: it corresponds to constraint equations rather than propagating modes. By performing a scalar-vector-tensor decomposition, we have been able to cleanly separate the Newtonian sector from the gravitational wave sector, and identify the physical degrees of freedom. In the scalar sector, the dynamical scalar \(\psi\) is promoted to a stochastic degree of freedom, while the Newtonian potential \(\Phi\) is not independently stochastic, but is fixed by the Hamiltonian constraint once \(\psi\) is specified. The same pattern holds in the vector sector: \(V_i\) satisfies a constraint equation, rather than an independent stochastic equation of motion. As a result, the indefiniteness of the gravitational action is benign: the sectors with positive-semi-definite two-point functions are precisely the dynamical ones, while the indefinite vector sector represents a constraint rather than a propagating mode.

This does not mean that \(V_i\) is determined once and for all. At second order, the vector-scalar mixing in the presence of matter sources means that \(V_i\) acquires stochastic contributions through the momentum constraint and has a dependence on \(\partial\psi\). Indeed, this must be the case if the momentum constraint is sourced by quantum fields in superposition and to preserve covariance. A fuller understanding of how the constraints are imposed in the path integral at higher order, where the vector-scalar coupling becomes dynamically relevant, remains an open question.

In addition, whereas classical, deterministic general relativity has only two dynamical modes, we now have an additional, scalar mode. This is to be expected, as the conservation of the energy-momentum tensor no longer suffices to make four of the remaining six modes vanish or cancel on-shell [@Spin2]. Further, we discussed possible experimental implications of a stochastic model of gravity, and whether such a model is ruled out by current experimental bounds on noise, and coherence times via the decoherence-diffusion trade-off. From the SVT decomposition, we are able to refine possible bounds from the LISA Pathfinder via the Newtonian potential, and from stochastic gravitational waves. We found that the theory is viable given current observational bounds, but care needs to be taken, given our poor understanding of a number of subtleties, such as the effects of nonlinearities, non-Markovianity, and scale-invariance. These calculations point the way to further theoretical improvements, which together with more precise experiments, can be used to test the quantum vs classical nature of gravity.

Indeed, this work highlights a subtlety with practical consequences. In the Newtonian limit with pressureless dust[@oppenheim2026non], our results differ significantly from those of [@layton2023weak], as well as [@tilloy2016sourcing]. The origin of this difference is illuminating. The work of [@layton2023weak] first reduces the gravitational phase space by solving the constraints classically, and then introduces stochasticity on the remaining degrees of freedom. Here, the stochastic dynamics acts on the full metric and the constraints emerge from the path integral. These two procedures do not commute, a tension that mirrors the long-standing debate in quantum gravity between Dirac quantisation and reduced phase space quantisation [@isham1992canonical; @kuchar1992time]. The present setting, being classical-stochastic rather than quantum, offers a considerably more tractable arena in which to study this inequivalence.

Acknowledgements: We would like to thank John Donoghue, Rhys Evans, Bob Holdom, Isaac Layton, Emanuele Panella, Geoffrey Penington, and Elizabeth Wilson for valuable discussions, as well as EPSRC (UKRI) for financial support. The calculations were performed by humans.

7 The Onsager Machlup and the Janssen–DeDominicis path integrals↩︎

In this section we motivate the stochastic gravitational path integral, and introduce an alternate formalism, the Janssen-de-Dominicis (JD) [@Janssen1976-hh],[@De_Dominicis1976-zy] or Martin-Siggia-Rose (MSR) formalism [@martin1973statistical]. We note that unless the diffusion matrix is conserved, the JD and the OM formalisms give us different two point functions. However, we show that even if we do take the ultra-local, generalized DeWitt metric Eq. 2 to be our diffusion matrix, despite the fact that it does not obey the Bianchi identity, we obtain the same OM action as that obtained from a conserved diffusion matrix. This is because saturating the inverse of the non-conserved DeWitt metric with the conserved Einstein tensor serves to cancel out the longitudinal parts of the DeWitt metric. This demonstrates the equivalence of the two formalisms. In addition, we discuss stochastic differential equation forms of the linearized Einstein equations, and present a local, conserved diffusion matrix which does satisfy the Bianchi identity. We thus see that we can have a local stochastic theory which nonetheless generates a valid spacetime.

A caveat before we start; many of the equations are formal, in the sense that they must be regularized before they can be used in the gravitational case. Our presentation mainly follows [@CriticalDynamics]. Consider a Langevin type equation of the form \[\begin{align} \partial_t^2\Phi(x,t)=\triangle\Phi(x,t)+\xi \end{align}\] where the covariance of \(\xi\), the Gaussian white noise is given by \[\begin{align} \langle \xi(x,t)\xi(x',t')\rangle=2D\delta(x-x') \end{align}\] where D could possibly be a function of \(t\) and \(x\). The probability distribution from which the correlations can be obtained is then given by \[\begin{align} \label{eq:IntroA1} W[\xi]=\mathcal{N}\int D\xi \exp({-\int d^4x \xi D_2^{-1}\xi/4}) \end{align}\tag{37}\] and so one obtains the OM path integral by a change of variables \[\begin{align} \label{eq:IntroA11} W[\Phi]\propto\int J D\Phi \exp{-\frac{1}{4}\int d^4x (\partial_t^2\Phi-\triangle\Phi) D_2^{-1}(\partial_t^2\Phi-\triangle\Phi)} \end{align}\tag{38}\] where the Jacobian \(J\) depends on the discretization used (see 7.3).

An entirely equivalent formalism is the Janssen–DeDominicis (JD), or the Martin-Siggia-Rose (MSR) response functional one, in which the covariance matrix, \(D_2\) does not appear in the denominator, and additionally, \(\Phi\) appears linearly. In this formalism, we have \[\begin{align} W[\Phi]\propto\int D[\Phi] D[i\tilde{\Phi}] \exp{-\int d^4x (\partial_t^2\Phi-\triangle\Phi)\tilde{\Phi}-\tilde{\Phi}D_2\tilde{\Phi}} \end{align}\] One can verify their equivalence by integrating over the Gaussian response functional.

We now replicate the steps above for the gravitational case twice: firstly for a non-conserved diffusion matrix, and secondly for a conserved one.

7.0.1 The Gravitational Case I: The deWitt Metric↩︎

This whole procedure is no longer as straightforward in the gravitational case. Suppose that we start from \[\begin{align} \label{eq:SD} G_{\mu \nu}^{(1)}&=\xi_{\mu \nu}+\xi'_{\mu\nu}\\ \langle\xi_{\mu \nu}\xi_{\rho\sigma}\rangle&=D_2(I_{\mu\nu,\rho\sigma}-\frac{\beta}{4\beta-D^{-1}_2}\eta_{\mu\nu}\eta_{\rho\sigma})\delta(x-x')^{4}\\\mathcal{D}_2&:= \langle\xi_{\mu \nu}\xi_{\rho\sigma}\rangle \end{align}\tag{39}\] where \(G_{\mu \nu}^{(1)}\) is the linearized Einstein tensor and \(I_{\mu\nu,\rho\sigma}:=\frac{1}{2}(\eta_{\rho\nu}\eta_{\sigma\mu}+\eta_{\rho\mu}\eta_{\sigma\nu})\). The first equation is a linearised version of that appearing in [@hu2008stochastic; @moffat1997stochastic] except here, the noise kernel is local. We here consider the vacuum case, since when matter is added, the dynamics takes on a different form in comparison to stochastic gravity.

There are two issues with this formulation. Firstly, \(\xi_{\mu\nu}\) is not a conserved tensor, as \(\partial^\mu\partial^\alpha\langle\xi_{\mu \nu}\xi_{\rho\sigma}\rangle\neq0\) unless the positivity constraint for the scalar sector, \(3\beta\leq \alpha\), is violated, and so we must add an additional non-conserved tensor for consistency such that 39 is consistent. Secondly, the deWitt metric can serve as a covariance matrix only in a formal sense given that it is not positive semi-definite; however, we have argued above 4.1 that after contracting the inverse diffusion matrix with the Einstein tensor, only the PSD sector of the action obtained is dynamical and so physically relevant.

Following the procedure of 37 ,38 we obtain the following, diffeomorphism invariant OM action \[\begin{align} \begin{aligned} \label{eq:OM}W[h^{\mu\nu}]&\propto\int J Dh^{\mu\nu} \exp{-\int d^4x \int d^4yG_{\mu\nu}(x)\mathcal{D}_2^{-1\mu\nu,\rho\sigma}(x,y)G_{\rho\sigma}(y)}\\&=\int d^4x (\alpha R^{\mu\nu}R_{\mu\nu}-\beta R^2) \end{aligned} \end{align}\tag{40}\]

Can we now claim that 40 implies the following JD action with the auxiliary, symmetric two-tensor \(\tilde{h}^{\rho\sigma}\), i.e.,, \[\begin{align} \label{eq:MSR}W[h^{\mu\nu}]&\propto^{?}\int J D[h] D[i\tilde{h}] \exp{-\int d^4x (G_{\rho\sigma})\tilde{h}^{\rho\sigma}-\frac{1}{4}\tilde{h}^{\mu\nu}\mathcal{D}_{2\mu\nu,\rho\sigma}\tilde{h}^{\rho\sigma}} \end{align}\tag{41}\] in the sense that both actions give us the same two point function?

Setting \(\alpha:=D_2^{-1}\), the JD action has the following Lagrangian density: \[\begin{align} \label{eq:ACTIONJD} \mathcal{L}_{JD}=G_{\rho\sigma}\tilde{h}^{\rho\sigma}-\frac{1}{4\alpha}(\tilde{h}^{\mu\nu}\tilde{h}_{\mu\nu}-\gamma \tilde{h}^2) \end{align}\tag{42}\] where \(\gamma=\frac{\beta}{4\beta-\alpha}\), which gives us the original action, \[\begin{align} \mathcal{L}_{OM}=\alpha R_{\mu\nu}R^{\mu\nu}-\beta R^2 \end{align}\]upon integrating out the auxiliary variable.

However, one does not recover the original propagator for any non-zero \(\alpha,\beta\). To see this, we can re-write the JD action in terms of spin-projection operators: \[\begin{align} \label{eq:JD} S_{JD}=\int d^4x h^{\mu\nu}\hat{O}_{\mu\nu,\rho\sigma}\tilde{h}^{\rho\sigma}-\frac{1}{4\alpha}\tilde{h}^{\mu\nu}(\frac{1}{2}(\eta_{\mu\rho}\eta_{\nu\sigma}+\eta_{\mu\sigma}\eta_{\nu\rho})-\gamma\eta_{\mu\nu}\eta_{\rho\sigma} )\tilde{h}^{\rho\sigma} \end{align}\tag{43}\] where \[\begin{align} \hat{O}_{\mu\nu,\rho\sigma}=\Box(-\frac{1}{2}P^{(2)}+P^{(0-s)})_{\mu\nu,\rho\sigma} \end{align}\] and the spin-projectors are given by 55 ,56 . Re-writing 43 in a matrix form, \[\begin{align} \begin{aligned} \label{eq:MJD} S_{JD}&=\begin{pmatrix}h \\ \tilde{h}\end{pmatrix}^T \mathcal{J} \begin{pmatrix}h \\ \tilde{h}\end{pmatrix}\\ \mathcal{J}&= \begin{bmatrix} \boldsymbol{0} &&\hat{\boldsymbol{O}}(i\omega,k)\\\hat{\boldsymbol{O}}(-i\omega,k) && -D_2\boldsymbol{I}(\gamma)' \end{bmatrix} \end{aligned} \end{align}\tag{44}\] with \[\begin{align} D_2{\boldsymbol{I}(\gamma)'}_{\rho\sigma,\mu\nu}= \frac{1}{4\alpha}(\frac{1}{2}(\eta_{\mu\rho}\eta_{\nu\sigma}+\eta_{\mu\sigma}\eta_{\nu\rho})-\gamma\eta_{\mu\nu}\eta_{\rho\sigma} ) \end{align}\] This matrix is non invertible, due to the non-invertibility of \(\hat{O}\). Adopting the gauge, \(\partial_{\mu}h^{\mu\nu}=0\), we have the altered operator, \[\begin{align} \hat{O}'_{\mu\nu,\rho\sigma}=\Box(-\frac{1}{2}P^{(2)}+P^{(0-s)}-\frac{1}{4\Delta}(P^{(1)}+2P^{(0-w)})_{\mu\nu,\rho\sigma} \end{align}\] where \(\frac{1}{\Delta}\) governs the strength of gauge-fixing. Inverting \(\mathcal{J}\) with \(\hat{O}'\) in place of \(\hat{O}\), we obtain for the propagator , \(\langle h_{\rho\sigma }h_{\mu\nu}\rangle\) \[\begin{align} \lim_{\Delta\rightarrow0} \langle h_{\rho\sigma }h_{\mu\nu}\rangle=\frac{1}{\alpha\Box^2}(P^{(2)}+\frac{(1-3\gamma)}{4}P^{(0)})_{\rho\sigma,\mu\nu} \end{align}\] One notes this is markedly different in its coefficients from the propagator obtained directly from the quadratic action, for any non-zero \(\alpha,\beta\) \[\begin{align} \label{eq:2ppp}\lim_{\Delta\rightarrow0} \langle h_{\rho\sigma }h_{\mu\nu}\rangle=\frac{1}{\alpha\Box^2}(P^{(2)}+\frac{P^{(0)}}{4(1-\frac{3\beta}{\alpha})})_{\rho\sigma,\mu\nu} \end{align}\tag{45}\] which one may verify by first integrating out \(\tilde{h}_{\mu\nu}\), \[\begin{align} \mathcal{L}(h^{\mu\nu})= (\Box(-\frac{1}{2}P^{(2)}+P^{(0-s)}))_{\alpha\beta,\rho\sigma}h^{\rho\sigma}\mathcal{D}_2^{-1\alpha\beta,\kappa\gamma}(\Box(-\frac{1}{2}P^{(2)}+P^{(0-s)}))_{\kappa\gamma,\mu\nu}h^{\mu\nu} \end{align}\] and then inverting the resulting action for the propagator.

One can also see this using gauge-invariant variables. Allowing for integration by parts, the scalar, vector and tensor sectors decouple. The action for the three sectors is of the generic form \[\begin{align} \mathcal{S}=\frac{1}{2}\int_{\omega,k}(h(-k,-\omega),\tilde{h}(-k,-\omega))\mathcal{J}'(h(k,\omega),\tilde{h}(k,\omega) \end{align}\] where \(\mathcal{J}'\) is given by \[\begin{align} \mathcal{J}'= \begin{bmatrix} \boldsymbol{0} &&\hat{\boldsymbol{G}}(i\omega,k)\\\hat{\boldsymbol{G}}(-i\omega,k) && -D_{2}\boldsymbol{I}'' \end{bmatrix} \end{align}\] and \((\hat{ G} h)_{\rho\sigma}\) is the scalar (respectively vector and tensor ) sector of the linearized Einstein action to which we have added an infinitesimal term of the form \(i\omega\) for regularization purposes and \(\mathbf{I''}\) is the coefficient matrix of the terms quadratic in scalar (vector and tensor) potentials.

However, it is no longer as simple to work with gauge-invariant potentials. The action in 40 is gauge invariant, whereas in 41 , only the first term is invariant under the following transformation of the auxiliary variables, \(\tilde{h}_{\mu\nu}\rightarrow\tilde{h}_{\mu\nu}+2\partial_{(\mu}\tilde{\xi}_{\nu)}\) (by the linearized Bianchi identity). Although the Einstein tensor can be written as a function of six variables, one has to contend with ten auxiliary variables corresponding to the ten potentials in a symmetric, rank-2 tensor in 4 dimensions. Further, even if one amends the second term in 41 by hand ² to an action of the form, for example for the scalar sector³, \[\begin{align} \label{eq:MSRSCALAR}S[\Phi,\psi,\tilde{\Phi},\tilde{\psi}]=\int d^4x G_{\rho\sigma}\tilde{h}^{\rho\sigma} (\Phi,\psi,\tilde{\Phi},\tilde{\psi})-(\tilde{\psi},\tilde{\Phi})^{T}\boldsymbol{A}(\alpha,\beta)(\tilde{\psi},\tilde{\Phi}) \end{align}\tag{46}\] one does not obtain 14 even with \(\beta\) set to zero on integrating out the auxiliary variables ,nor the correlation matrix 16 upon inverting 46 . Instead, the scalar JD action corresponding to 14 is given by an action of the form \[\begin{align} S[\Phi,\psi,\tilde{\Phi},\tilde{\psi}]_{JD}=-\int d^4x \Omega_1(a_1\partial_0^2\psi+\triangle(b_1\psi+c_1\Phi))\tilde{\psi}+\Omega_2(a_2\partial_0^2\psi+\triangle(b_2\psi+c_2\Phi))\tilde{\Phi}-(\alpha-\beta)(\tilde{\Phi}\tilde{\Phi}+\tilde{\psi}\tilde{\psi}) \end{align}\] where \(\Omega_i,a_i,b_i,c_i\) are functions of \(\alpha,\beta\) that may be obtained by diagonalizing the action in 14 in the \(\triangle\psi,\partial_0^2\psi,\Phi\) basis.⁴

The vector sector proves to be much more instructive. Based on the action 13 , \[\begin{align} \mathcal{L}_{vector}=V_i\triangle \tilde{V^i}-(-\tilde{S_i}\tilde{S^i}-\tilde{F_i}\triangle\tilde{F^i}) \end{align}\] One can easily verify that inverting the action does not give us the two-point function we would obtain from inverting 13 . Interestingly, if we ‘gauge-fix’ here by hand by setting \(\tilde{\dot{F_i}}=0\), the resulting action is non-dynamical. Let us recall why that is the case. The physically meaningful quantity is not the vectorial sector of \(h_{\mu\nu}\), but rather the special combination, \(V_i=S_i-\dot{F_i}\). The second term in 42 however, breaks diffeomorphism invariance of the action and thus makes meaningful all ten components of the auxiliary variable. This is of particular importance for the vectorial sector, for the original equation of motion (for \(V_i\)) is non-dynamical. For the action 42 , we then have two choices; either we preserve diffeomorphism invariance ⁵ and work with the invariant auxiliary variables, \(V_i\) or we break it, as in 42 . Neither choice recovers the action in 13 .

Finally, for the tensorial sector no complication arises, as \(\tilde{h}_{\mu\nu}^{TT}\) is diff-invariant.

7.0.2 The Gravitational Case II: A Conserved Diffusion Matrix↩︎

In this section, we perform the calculations of the above section using a conserved diffusion matrix instead. This approach thus, resolves one of the two issues with 39 , of 39 violating the Bianchi identity, while also resolving the inconsistency between the propagators obtained from the OM and the JD action.

Inspired by the work of Hirotani and Matsumura [@hirotani2026testing], we consider the alternative stochastic differential equation (SDE) formulation, \[\begin{align} \begin{aligned} \label{SD:BR} G_{\mu \nu}^{(1)}&=\xi_{\mu \nu}\\\langle\xi_{\mu \nu}\xi_{\rho \sigma}\rangle&=D_2((1-3\frac{\beta}{D_2^{-1}})^{-1}P^{(0-s)}+P^{(2)})_{\mu\nu,\rho\sigma}\delta^4(x-x')\\\langle\xi_{\mu \nu}\xi_{\rho \sigma}\rangle&=\mathcal{D}'_{2\rho\sigma,\mu\nu} \end{aligned} \end{align}\tag{47}\] This differs from the proposal of  [@hirotani2026testing] (c.f. [@hu2008stochastic]) in that it is a local kernel. Although this matrix is also indefinite, one can easily verify that the diffusion matrix, is indeed conserved, by the transversality of the projectors used here. This matrix is however, degenerate, and additional ‘gauge-fixing’ terms must be added to invert it to for example give us an inverse diffusion matrix of the form \[\begin{align} \label{eq:Diff2} (\langle\xi\xi\rangle)^{-1 \mu\nu,\rho\sigma}&=D_2^{-1}((1-3\frac{\beta}{D_2^{-1}})P^{(0-s)}+P^{(2)}+\frac{1}{\Delta}(P^{0-w}+P^{(1)}))^{\mu\nu,\rho\sigma}\delta^4(x-x') \end{align}\tag{48}\] or for example \[\begin{align} (\langle\xi\xi\rangle)^{-1 \mu\nu,\rho\sigma}&=D_2^{-1}((1-3\frac{\beta}{D_2^{-1}})P^{(0-s)}+P^{(2)}-(\frac{(4(1-\frac{3\beta}{D_2^{-1}})-3\Delta))P^{(0-w)}}{\Delta(1-3\frac{3\beta}{D_2^{-1}})}+\frac{2P^{(1)}}{\Delta}+\sqrt{3}\frac{P^{(sw)}+P^{(ws)}}{(1-3\frac{3\beta}{D_2^{-1}})}))^{\mu\nu,\rho\sigma}\delta^4(x-x') \end{align}\] where we can set the parameter \(\Delta\) to be as large as we please. The key here is that for all the choices above, the inverse diffusion matrix is saturated by a conserved tensors, and so the gauge-fixing terms will all drop out in the limit, \(\Delta\rightarrow\infty\). Thus, only the transverse part of the inverse diffusion matrix remains and so we obtain the same action as was obtained from the deWitt metric: \[\begin{align} S'_{G}= \int d^4x \int d^4y G_{\mu\nu}(x) \mathcal{D}_2^{'-1\mu\nu,\rho\sigma}G_{\rho\sigma}(y)=\int d^4x \alpha R^{\mu\nu}R_{\mu\nu}-\beta R^2. \end{align}\] Thus, in the OM formalism, both choices of diffusion matrix give us the same action, and consequently, the same two-point function.

Now, applying the JD formalism with \(\mathcal{D}_2'\), we have the JD action, \[\begin{align} \label{eq:JD2} S'_{JD}=\int d^4x h^{\mu\nu}\hat{O}_{\mu\nu,\rho\sigma}\tilde{h}^{\rho\sigma}-\frac{1}{4\alpha}\tilde{h}^{\mu\nu}((1-3\frac{\beta}{\alpha})^{-1}P^{(0-s)}+P^{(2)})_{\mu\nu,\rho\sigma}\tilde{h}^{\rho\sigma}\, , \end{align}\tag{49}\] giving us (upto gauge terms), the following two-point function, \[\begin{align} \langle h^{\mu\nu}h^{\rho\sigma}\rangle=\frac{1}{\alpha\Box^2}(P^{(2)}+\frac{1}{4(1-\frac{3\beta}{\alpha})}P^{(0-s)})^{\mu\nu,\rho\sigma} \end{align}\] This matches perfectly with 45 .

Note that in the JD approach, since we do not have to invert our diffusion matrix, there is no need to add ‘gauge-fixing’ terms to it.

7.0.3 The reason for the inconsistency↩︎

The inconsistency between the propagators obtained from the OM action, and the JD action when using the diffusion matrix 39 , arises from the fact that the action, 43 is not invariant under linearized diffeomorphisms,while 5 is. This follows from the fact that the diffusion matrix, 39 is not a conserved tensor. One can indeed write it in terms of projection operators: \[\begin{align} \begin{aligned} \label{eq:spindiff} D_2(I_{\mu\nu,\rho\sigma}-\frac{\beta}{4\beta-D^{-1}_0}\eta_{\mu\nu}\eta_{\rho\sigma})\delta(x-x')^{4}&=D_2(P^{(0-s)}+P^{(0-w)}+P^{(1)}+P^{(2)}-\frac{\beta}{4\beta-D^{-1}_2}(3P^{(0-s)}\\&+\sqrt{3}(P^{(ws)}+P^{(ws)})+P^{(0-w)})_{\mu\nu,\rho\sigma}\delta(x-x')^{4} \end{aligned} \end{align}\tag{50}\] The longitudinal terms thus render the entire action non-gauge invariant. The OM action on the other hand, remains gauge-invariant, due to the fact that it is saturated by the Einstein tensor, which is conversed. Saturating the inverse diffusion matrix with the Einstein tensor, and then inverting the resultant action is therefore not equivalent to saturating the diffusion matrix with the inverted Einstein-tensor for non-conserved matrices. In order to have a diffeomorphism invariant JD path integral, and thus consistent two point functions, one must use a conserved diffusion matrix. On the other hand, as we have noted above, this issue simply does not arise in the OM formalism.

7.1 The Gravitational Path Integral Measure↩︎

In the section, we review the path integral measure for the gravitational action, largely following the presentation in [@MAZUR]. The metric, \(g\) may be considered as a point in the Wheeler-DeWitt superspace, and the measure constructed on the space of the deformations at that given point. The space of metrics at a point (before quotienting by diffeomorphisms) is \(GL(D,R)/SO(D-1,1)\), which is not a vector space, however, we may instead consider the tangent space to \(g\) which is a vector space, and then use the exponential map to arrive at the full manifold of metrics (assuming that the space of metrics is contractible). Thus, for each topology, we would have \(Dg=D\delta g\), where \(\delta g\) spans the tangent space of \(g\). Now we would like to separate out the physical modes in \(\delta g\) from the gauge ones.The key point here is that the probability measure is not to be on the space of metrics, but on the space of metrics modulo diffeomorphisms. In fact, we must go further, and set it to be on the space of continuous geometries, i.e., the action must only sum over continuous metric configurations (modulo diffeomorphisms). We must thus factor out the volume of the diffeomorphism group, as in the quantum gravitational path integral. Further, we must be mindful of the fact that although the linearized gravitational action is of the form of a Wiener measure, it is unbounded from below and must be carefully regularized. As a reminder, the diffeomorphism \(x\rightarrow x+\xi\) results in the following metric transformation: \[\begin{align} \delta_\xi g_{\mu\nu}=2\nabla_{(\mu}\xi_{\nu})\; \end{align}\] which can be rewritten in the following form \[\begin{align} \delta_\xi g_{\mu\nu}=\frac{1}{2}\nabla_{\alpha}\xi^{\alpha}g_{\mu\nu}+(L\xi)_{\mu\nu} \end{align}\] where \((L\xi)_{\mu\nu}\) is the traceless part of \(\delta_\xi g_{\mu\nu}\) or equivalently, the operator \(L\) maps the vector \(\xi_{\mu}\) to the traceless subspace. It can further be shown that the kernel of the operator \(L^\dagger\) is spanned by transverse, traceless metric perturbations, i.e., \(h_{\alpha}^{\alpha TT}=\nabla^{\mu}h_{\mu\nu}^{TT}=0\). We can thus decompose the metric perturbation, as \[\begin{align} h_{\mu\nu}=\frac{1}{4}hg_{\mu\nu}+h_{\mu\nu}^{TT}+(L\xi)_{\mu\nu} \end{align}\] and check that is orthogonal with respect to our chosen metric on the Wheeler-deWitt superspace, \(G_{\rho\sigma\mu\nu}=\frac{1}{2}(g^{\rho\mu}g^{\sigma\nu}+g^{\rho\nu}g^{\sigma\mu}-\beta g^{\rho\sigma}g^{\mu\nu})\). Then one can define the functional measure with respect to this metric \[\begin{align} \int D\delta g e^{- i\langle \delta g,\delta g\rangle}=1 \end{align}\] where \(\langle \delta g,\delta g\rangle=\int_\mathcal{M} d^4x\delta gG\delta g\) Thus, starting from a path integral of the form \[\begin{align} \int D\delta g e^{-S[g]} \end{align}\] one has two equivalent choices. The first is the usual Faddeev-Poppov procedure; adding gauge fixing terms and factoring out the volume of diffeomorphism group from the measure, while in the second one can decompose \(\delta g\) in the form given above and factor out the part corresponding to range \(L\) . Thus, \[\begin{align} \int D\delta g e^{-S[\delta g]}=\int J_0[Dh]J_1[D\xi]J_2[Dh^{TT}] e^{-S[g]} \end{align}\] The only non-trivial Jacobian is \(J_1=Det(L^\dagger L)\) which in flat space is given by \(J_1={Det}(\Box)^{1/2}_S{Det}(\Box)^{1/2}_V\). Thus, the path integral is given by \[\begin{align} \int J_0[Dh]J_1[D\xi]J_2[Dh^{TT}] e^{-S[g]} =Vol(Diff(M))J_1\int DhDh^{TT} e^{-S[g]} \end{align}\] but this is evidently the same as \[\begin{align} Vol(Diff(M))J_1\int DhDh^{TT} e^{-S[g]}=K_a\int DVD ^3h^{TT}D \Phi D\psi e^{-S[g]} \end{align}\] where \(K_a\) is the product of the volume of the diffeomorphism group with the Jacobian of transformations from \(h_{\mu\nu}\) to the gauge invariant coordinates, \((\psi, \Phi, ^3h^{TT},V)\). What we have essentially done is to write the \(5+1\) modes of \((^4h^{TT},h)\) in a new, non-manifestly covariant basis. Thus, \(K_a\) is equivalent to \(J_1\) upto the group volume. Further, \({Det}(\Box)^{1/2}_V={Det}(\Box)^{3/2}\) is the determinant over the three transverse vectors in \(^4h^{TT}\). We can use this then to check if the non-manifestly covariant approach gives us the correct number of ‘dynamical’ modes (dynamical, in the usual, non-stochastic sense, where each box operator corresponds to a dynamical mode): each factor of \({Det}(\Box)^{1/2}\) transforms a mode into a constraint, and hence we have \(2*5+2-4=8\) modes or box operators. This is indeed the same number of box operators as we obtain in the non-manifestly covariant formulation.

7.2 Choice of discretization↩︎

Here we show that the choice of discretization does not matter, mainly based off the presentation in [@CriticalDynamics]. Consider the following Langevin equation \[\begin{align} \frac{\phi_{i+1}-\phi_{i}}{\epsilon}=\alpha F_{i+1}+(1-\alpha)F_{i}+\xi_i \end{align}\] where \(0\leq \alpha \leq1\) and \(F_i\) is the force at \(t_i\), and note that the choice of discretization is equivalent to setting the value of the Heaviside function \(\Theta=\alpha\) The Jacobian of the change of variable from \(\phi_{i+1}\rightarrow\xi_i\) is a lower triangular matrix with entries only on the main, and the lower sub-diagonal. Thus the functional determinant is given by \[\begin{align} \frac{1}{\epsilon^f} \prod_i(1-\alpha\epsilon\frac{\partial F_i}{\partial \phi_{i}})\approx\frac{1}{\epsilon^f}\exp(-\alpha\sum_i\epsilon\frac{\partial F_i}{\partial \phi_{i}}) \end{align}\] Taking the limit \(\epsilon\rightarrow0\) and absorbing the factor \(\frac{1}{\epsilon^f}\) in the functional measure, we obtain the term \(-\alpha\int d^3x\int dt\frac{\partial F}{\partial \phi}\). This is exactly canceled by the contribution from closed response loops contained in (the non-linear contribution to) terms of the form \(\int d^3x\int dt F\tilde{\phi}\), which necessarily involve a contraction of \(\tilde{\phi}\) with the internal field \(\phi\) in the JD formalism \[\begin{align} \int d^3x\int dt \phi\tilde{\phi}\frac{\partial F}{\partial \phi}=\int d^3x\int dt \Theta(0)\frac{\partial F}{\partial \phi} \end{align}\] as long as we set \(\alpha=\Theta(0)\).

7.3 Second Order SDEs and Heat Baths↩︎

Since the stochastic path integral in either of its iterations (OM or JD) is formulated using Stochastic differential equations that are first order in time, let us explicitly write down the second order equations as first order ones. It is convenient to do so, as working with second order SDEs requires one to compute a Jacobian that is non-lower triangular. Our exposition mainly follows [@Chaichian2018PathI]. Working in spatial Fourier space we have \[\begin{align} \begin{aligned} \partial_t\varphi&=\pi_{\varphi}\\ \partial_t\pi_{\varphi}-k^2\varphi&=\xi \end{aligned} \end{align}\] We can consider this as a coupled system of the variables, \(\varphi,\pi_{\varphi}\), undergoing stochastic diffusion in baths of different temperature, \(T_1,T_2\). \[\begin{align} \begin{aligned} \label{eq:SDE11}\partial_t\varphi-\pi_{\varphi}&=T_1\xi\\ \partial_t\pi_{\varphi}-k^2\varphi&=T_2\xi \end{aligned} \end{align}\tag{51}\] Compare this with the SDE for Brownian motion of two particles:\[\begin{align} \begin{aligned} \label{eq:SDE111}\partial_t\vartheta_1&=T_1\xi \\ \partial_t\vartheta_2&=T_2\xi\end{aligned} \end{align}\tag{52}\] We can write 51 , 52 in matrix form: \[\begin{align} \begin{aligned} \partial_t\boldsymbol{\vartheta}&=\boldsymbol{\xi}\\\partial_t\boldsymbol{\varphi}+\boldsymbol{D}\boldsymbol{\varphi}&=\boldsymbol{\xi} \end{aligned} \end{align}\] where \[\begin{align} \boldsymbol{D}=\begin{bmatrix} 0&-1\\-k^2&0 \end{bmatrix} \end{align}\] we obtain the Jacobian, \[\begin{align} J_N= \begin{bmatrix} \boldsymbol{1}+\boldsymbol{D}\frac{\epsilon}{2}\;&.&.&.&0\\\boldsymbol{D}\frac{\epsilon}{2}&\boldsymbol{1}+\boldsymbol{D}\frac{\epsilon}{2}&.&.&0\\.&.&.&.&.\\.&.&.&.&.\\\boldsymbol{D}\frac{\epsilon}{2} &.&.&.&\boldsymbol{1}+\boldsymbol{D}\frac{\epsilon}{2} \end{bmatrix}=\exp{\frac{1}{2}t\text{Tr}\boldsymbol{D}}=1 \end{align}\] where we have used the fact that the coefficient matrix is lower triangular. The Jacobian only gives us unity in this case. We thus have the following probability distribution, \[\begin{align} W(\varphi_f,\pi_{\varphi_f}|\varphi_i,\pi_{\varphi_i}) =\lim_{T_1\rightarrow 0, T_2\rightarrow1}\int D\varphi D\pi_{\varphi}\exp{-\int dt d^3k (\frac{1}{T_1}( \partial_t\varphi-\pi_{\varphi})^2+\frac{1}{T_2}( \partial_t\pi_{\varphi}-k^2\varphi)^2)} \end{align}\] This gives us a delta function at every time slice, which we may use to integrate out \(\pi_{\varphi}\) for intermediate times. In addition one must integrate over all states, \(\pi_{\varphi_i},\pi_{\varphi_f}\) if we want the probability distribution to be solely a function of \(\varphi_i,\varphi_f\).

7.3.1 Ensuring continuity↩︎

Based on the discussion above on the equivalence of diffusion constants and temperatures of the heat bath, we argue that the unique choice to ensure continuity in 4.2 is by adding a delta function at each time slice. Clearly, adding a general term of the form \(\mathcal{S}_{dyn}=\frac{1}{2D}\sum_k\epsilon (\frac{q_{k+1}-q_k}{\epsilon}-F(q_k))^2\) would suffice to ensure continuity, but with what diffusion constant or temperature, D, and which \(F(q_k)\)? Clearly the unique choice that ensures continuity is to set the temperature, D, to zero, and \(F(q_k)=0\), with the second equality also serving to preserve the mean. This then gives us a delta function, \(\delta(q_{k+1}-q_k)\) at every time step.

7.4 Pole Prescription↩︎

Here we discuss the pole-prescriptions for the two-point function, as justified in [@oppenheim2025diffusion]. The pole-prescription for the infinite time two-point function is given by \[\begin{align} \label{eq:pole-prescription} \frac{1}{p^4}=\frac{1}{(k^2-(\omega-i\epsilon)^2)(k^2-(\omega+i\epsilon)^2)} \end{align}\tag{53}\] We can obtain this directly from the action 43 , where the limits of integration are from \(\pm\infty\), and we can work in momentum space: \[\begin{align} \hat{O}_{\mu\nu,\rho\sigma}=\frac{1}{2}k^2(i\epsilon)(-P^{(2)}+P^{(0-s)})_{\mu\nu,\rho\sigma} \end{align}\] where we directly add to the (Fourier transform of the) box-operator the term in \(i\epsilon\) corresponding to the retarded prescription for its Green’s function. Then, we can integrate out the auxiliary variable tensor field, \(\tilde{h}^{\mu\nu}\) to obtain the following action, \[\begin{align} \mathcal{S}_{G}=\int d^4k\int d^4k'\frac{1}{4}k^2(i\epsilon)(-P^{(2)}+P^{(0-s)})_{\mu\nu,\kappa\gamma}h^{\kappa\gamma}(k) D_2^{-1}I^{\mu\nu,\alpha\beta}(\beta,(k,k'))k^2(i\epsilon)(-P^{(2)}+P^{(0-s)})_{\alpha\beta,\rho\sigma}h^{\rho\sigma}(k') \end{align}\]where \[\begin{align} I^{\mu\nu,\rho\sigma} (\gamma,(k,k'))=\delta^4(k+k')(\frac{1}{2}(\eta^{\mu\rho}\eta^{\nu\sigma}+\eta^{\mu\sigma}\eta^{\nu\rho})-\frac{\beta}{\alpha}\eta^{\mu\nu}\eta^{\rho\sigma}) \end{align}\]Integrating over \(k'\), we obtain \[\begin{align} \begin{aligned} \mathcal{S}_{G}&=\int d^4k\frac{1}{4}k^2(i\epsilon)(-P^{(2)}+P^{(0-s)})_{\mu\nu,\kappa\gamma}h^{\kappa\gamma}(k) D_2^{-1}\delta^{\mu\nu,\alpha\beta}(\beta,(k,k'))k^2(-i\epsilon)(-P^{(2)}+P^{(0-s)})_{\alpha\beta,\rho\sigma}h^{\rho\sigma}(k')\\&=\int d^4k (\alpha R^{\mu\nu}R_{\mu\nu}(k)-\beta R^2(k)) \end{aligned} \end{align}\]From the equivalence of the two point functions obtained from the scalar-vector-tensor decomposition, and the covariant two point function, the pole prescription for the Box-operator \(p^2\) must be identical, and given by 53 . That this is indeed the correct pole-prescription may also be directly verified via 44 , without integrating out auxiliary fields.

Similarly, for the finite time case, one may treat the spin-0 propagator as the convolution of two retarded Green’s functions, \(G(p)_{\rho\sigma}=\frac{\theta_{\rho\sigma}}{\sqrt{3(\alpha-3\beta)}p^2}\). The pole-prescription for the pole coming from the projection operator, \(\theta_{\rho\sigma}\) is also retarded: the projector does not give us an additional pole.

We can also go back to the scalar-vector-tensor decomposition, and consider the action is given by 20 ,21 . Then, we simply assume that the two-point function is the convolution of two-retarded Green’s function. If we do not wish to integrate out \(\Phi\), we can instead consider 16 , and consider it to be of the form \(\mathcal{G}(p)\mathcal{G}(p')\) where \(\mathcal{G}(p)\) is the matrix of retarded Green’s function for the stochastic differential equations giving rise to the action, 14 . The exact form of \(\mathcal{G}\) does not concern us apart from the fact that it is proportional to \(\frac{1}{p^2}\), where we must use the retarded prescription.

Additionally, given that the matrix \(\mathcal{G}(p)\mathcal{G}(p')\) must be symmetric in \(p,p'\), we can rewrite the \(\langle\Phi\Phi\rangle\) element in 16 as \[\begin{align} \langle\Phi(p)\Phi(p')\rangle=\frac{1}{8\alpha(\alpha-3\beta)} (\frac{8(\alpha-3\beta)}{k^2k'^2}+\frac{4\omega^2\omega'^2(\alpha-3\beta)}{k^2k'^2p^2p'^2}+\frac{(-2\alpha+8\beta)}{p^2p'^2}) \end{align}\]

7.5 Calculational Details for the Spectral Acceleration↩︎

For the two point function, we assume that it is the convolution of two retarded Green’s functions: \[\begin{align} \langle\Phi \Phi \rangle=\int ^{\infty}_{-\infty}d^3z\int^{t_f}_0dz^0G_{ret}(x-z)G_{ret}(y-z)=\int ^{\infty}_{-\infty}d^3z\int^{t_f}_0dz^0\int\frac{d^4p}{(2\pi)^4}\int\frac{d^4p'}{(2\pi)^4}G(p)G(p') \end{align}\] where \[\begin{align} G(p)G(p')=\frac{1}{8\alpha(\alpha-3\beta)} (\frac{8(\alpha-3\beta)}{k^2k'^2}+\frac{4\omega^2\omega'^2(\alpha-3\beta)}{k^2k'^2p^2p'^2}+\frac{(-2\alpha+8\beta)}{p^2p'^2}) \end{align}\] Using \[\begin{align} \label{eq:time} \int^\infty_0e^{-iz^0(\omega+\omega')}dz^0=\pi\delta(\omega+\omega')-iP.V(1/(\omega+\omega')) \end{align}\tag{54}\] we obtain, \[\begin{align} \int ^{\infty}_{-\infty}d^3z\int^{t_f}_0dz^0\int\frac{d^4p}{(2\pi)^4}\int\frac{d^4p'}{(2\pi)^4}G(p)G(p')=\mathcal{G}_{stationary}+\mathcal{G}_{non-stationary} \end{align}\] where \[\begin{align} \mathcal{G}_{stationary}=\int\frac{d^4p}{(2\pi)^4}\frac{1}{8\alpha(\alpha-3\beta)} (8(\alpha-3\beta)/k^4+4\omega^4(\alpha-3\beta)/k^4p^4+(-2\alpha+8\beta)/p^4) \end{align}\] where \(p^4=(k^2-(\omega-i\epsilon)^2)(k^2-(\omega+i\epsilon)^2)\) Meanwhile, for \(\mathcal{G}_{non-stationary}\) we need only consider the non-stationary parts of terms of the form \(\omega^2\omega'^2/k^2k'^2p^2p'^2,1/p^2p'^2\).Let us start with \(\omega^2\omega'^2/k^2k'^2p^2p'^2\): the fourier transform of the non-stationary part is given by \[\begin{align} \mathcal{FT}(\omega^2\omega'^2/k^2k'^2p^2p'^2)_{non-stationary}=-\int\frac{d^3k}{(2\pi)^3} (e^{-i\vec{k}.(\vec{x}-\vec{y})})e^{-\epsilon (x^0+y^0)}\frac{1}{4k^6}(4k^3\sin(k(x^0+y^0))+k^4\cos(k(x^0-y^0))/\epsilon) \end{align}\] where we have only kept terms up to \(e^0\) (for simplicity, we have kept the exponential is it is, but later on, when expanding, we will only keep terms up to an overall \(e^0\) ). Integrating over angles; \[\begin{align} \mathcal{FT}(\omega^2\omega'^2/k^2k'^2p^2p'^2)_{non-stationary}= - \int\frac{dke^{-\epsilon (x^0+y^0)}}{8\pi^2|x-y|}\frac{1}{k^2}(4\sin(k(x^0+y^0))+k\cos(k(x^0-y^0))/\epsilon)\sin{(k(|x-y|))} \end{align}\] Similarly, for \(1/p^4\), \[\begin{align} \mathcal{FT}(1/p^2p'^2)_{non-stationary}= - \int\frac{dke^{-\epsilon (x^0+y^0)}}{8\pi^2|x-y|}\frac{1}{k^2}(\sin(k(x^0+y^0))+k\cos(k(x^0-y^0))/\epsilon)\sin{(k(|x-y|))} \end{align}\] Use the following integrals: \[\begin{align} \begin{aligned} \int^\infty_0 \frac{1}{k^2}\sin{ka}\sin{kb}=ab\int^\infty_0 \text{sinc}ka \text{sinc}kb=\frac{ab\pi}{|a-b|+|a+b| }&\\ \int^\infty_0 \frac{1}{k}\sin{ka}\cos{kb}=a\int^\infty_0 \text{sinc} (ka)\cos{kb}==\frac{\pi}{4}(-\text{sgn}(a-b)+\text{sgn}(a+b)) &\\\int^{\infty}_{-\infty}\frac{\pi}{4}(-\text{sgn}(a-b)+\text{sgn}(a+b))e^{i\omega a}da=\pi\sin{\omega b}/(\omega) \end{aligned} \end{align}\]

Thus, we obtain the following time-dependent parts, in frequency-space where by \(\mathcal{FT}_{\omega}\) we refer to the full Fourier transform, followed by an inverse Fourier transform with respect to \((x^0-y^0)\) \[\begin{align} \lim_{\epsilon\rightarrow0} \mathcal{FT}_{\omega}(1/p^4)|=-\frac{1}{8\pi^2|x-y|}(\frac{\sin{\omega|x-y|}}{\omega}(\frac{1}{\epsilon}-(x^0+y^0))+\frac{\pi|x-y|(x^0+y^0) }{|x^0+y^0-|x-y||+|x^0+y^0+|x-y||}2\pi\delta(\omega)) \end{align}\] where we have \[\begin{align} \delta(\omega)=\frac{1}{2\pi}\int^\infty_{-\infty}d(x^0-y^0)e^{i\omega(x^0-y^0)} \end{align}\] and a similar expression from the time-dependent term of the form \(\omega^4/k^4p^4\). For the time-independent part, we have for \[\begin{align} \begin{aligned} \mathcal{FT}_{\omega}\frac{\omega^4}{k^4p^4}(x,\omega)=\frac{\omega^4}{8\pi x(\omega^2+\epsilon^2)^4}((-x^2\epsilon^4+4\omega^2-x^2\omega^4-2\epsilon^2(2+x^2\omega^2))+&\\\frac{e^{-\epsilon x}}{\epsilon\omega}(\sin{\omega x}(\epsilon^4-6\epsilon^2\omega^2+\omega^4)+\epsilon\cos(\omega x)(4\epsilon^2\omega-4\omega^3))) \end{aligned} \end{align}\] and \[\begin{align} \lim_{\epsilon\rightarrow0}\mathcal{FT}_{\omega}\frac{\omega^4}{k^4p^4}=&\frac{1}{8\pi x}(-x^2+\frac{4}{\omega^2}+(\frac{\sin{\omega x}}{\omega}(\frac{1}{\epsilon}-x)-\frac{4\cos(\omega x)}{\omega^2})) \end{align}\] \[\begin{align} \lim_{\epsilon\rightarrow0}\mathcal{FT}_{\omega}\frac{1}{p^4}=\frac{\sin{\omega|x|}}{8\pi\omega|x|}(1/\epsilon-|x|) \end{align}\] We thus, obtain for the (wide-sense) stationary part \[\begin{align} \int d^3k\frac{1}{(2\pi)^3(k^2-(\omega-i\epsilon)^2)(k^2-(\omega+i\epsilon)^2)}=\frac{1}{(2\pi)^2i|x|}\int ^{\infty}_{-\infty}dk(\frac{ke^{ik|x|}}{(k^2-(\omega-i\epsilon)^2)(k^2-(\omega+i\epsilon)^2)})&\\=\frac{\sin{\omega|x|e^{-\epsilon|x|}}}{8\pi\omega\epsilon|x|} \end{align}\]And similarly, \[\begin{align} \int d^3k\frac{k^2}{(2\pi)^3(k^2-(\omega-i\epsilon)^2)(k^2-(\omega+i\epsilon)^2)}=\frac{1}{(2\pi)^2i|x|}\int ^{\infty}_{-\infty}dk(\frac{k^3e^{ik|x|}}{(k^2-(\omega-i\epsilon)^2)(k^2-(\omega+i\epsilon)^2)})&\\=\frac{((\omega^2-\epsilon^2)\sin{\omega|x|+2\omega\epsilon\cos{\omega|x|})e^{-\epsilon|x|}}}{8\pi\omega\epsilon|x|} , \end{align}\]\[\begin{align} \int d^3k\frac{1}{(2\pi)^3k^2(k^2-(\omega-i\epsilon)^2)(k^2-(\omega+i\epsilon)^2)}=\frac{1}{(2\pi)^2i|x|}\int ^{\infty}_{-\infty}dk(\frac{e^{ik|x|}}{(k)(k^2-(\omega-i\epsilon)^2)(k^2-(\omega+i\epsilon)^2)})&\\=\frac{e^{-\epsilon|x|}}{8\pi\omega\epsilon|x|} (\frac{(\omega^2-\epsilon^2)\sin{\omega|x|}-2\omega\epsilon\cos{\omega|x|}+2e^{\epsilon|x|}\omega\epsilon}{(\omega^2+\epsilon^2)^2})) \end{align}\] where we have taken the principal value. Finally, \[\begin{align} \lim_ {\epsilon \rightarrow0} S(\vec{x},\omega)=\frac{1}{64\pi\alpha(\alpha-3\beta)|x-y|}((2\alpha-4\beta)\omega\sin{\omega|x-y|(\frac{1}{\epsilon}-|x-y|)}+2(20\beta-6\alpha)\cos{\omega|x-y|}+2(12\alpha-36\beta)) \end{align}\]One thus obtains \[\begin{align} \begin{aligned} \langle \Phi \Phi\rangle(|x-y|,x^0+y^0,\omega)&=\frac{1}{64\eta\alpha\pi}(-12\eta|x-y|+\frac{16\eta}{\omega^2|x-y|} +\frac{(2\alpha- 4\beta)\sin{\omega |x-y|}}{\omega|x-y|}(x^0+y^0-|x-y|)\\&-\frac{\cos{\omega |x-y|}}{\omega^2|x-y|}(16\eta)-\frac{2\pi\delta(\omega)(14\alpha-40\beta)(x^0+y^0)}{(|x^0+y^0-|x-y||+|x^0+y^0+|x-y||)}) \end{aligned} \end{align}\] For the spectral density, we have \[\begin{align} \begin{aligned} \langle\partial_i \Phi \partial^i\Phi\rangle&=-\nabla^2\langle \Phi \Phi\rangle \\&=\frac{1}{64\eta\alpha\pi}(\frac{24\eta}{|x-y|}+(2\alpha-4\beta)(\frac{\omega\sin{\omega|x-y|}}{|x-y|}(x^0+y^0-|x-y|))+\frac{(40\beta-12\alpha)\cos{\omega |x-y|}}{|x-y|}) \end{aligned} \end{align}\] where the second step follows by the linearity of the expectation value.

7.6 Spin Projectors in Linearized Gravity↩︎

The Barnes–Rivers spin projection operators provide a standard decomposition of symmetric rank-2 tensor fields in momentum space [@barnes1963thesis; @rivers1964lagrangian], and have been widely used in higher-derivative gravity theories [@van1973ghost; @julve1978quantum].

7.7 Transverse and Longitudinal Projectors↩︎

Given a momentum \(k^\mu\) with \(k^2 = k^\mu k_\mu\), define the transverse and longitudinal rank-2 projectors: \[\begin{align} \theta_{\mu\nu} &= \eta_{\mu\nu} - \frac{k_\mu k_\nu}{k^2}, \\ \omega_{\mu\nu} &= \frac{k_\mu k_\nu}{k^2}, \end{align}\] satisfying \(\theta_{\mu\nu} + \omega_{\mu\nu} = \eta_{\mu\nu}\), \(\theta_{\mu}{}^{\alpha}\theta_{\rho\nu} = \theta_{\mu\nu}\), \(\omega_{\mu}{}^{\alpha}\omega_{\rho\nu} = \omega_{\mu\nu}\), \(\theta_{\mu}{}^{\alpha}\omega_{\rho\nu} = 0\).

7.8 Space-time Representation↩︎

We may write the following in space-time as \[\begin{align} \theta_{\mu\nu} &= \eta_{\mu\nu} - \frac{\partial_{\mu}\partial_{\nu}}{\Box^2}, \\ \omega_{\mu\nu} &= \frac{\partial_{\mu}\partial_{\nu}}{\Box^2}, \end{align}\] where the expressions above are to be understood more properly as the Fourier transform of the momentum-space projectors.

7.9 Barnes–Rivers Spin Projection Operators↩︎

The identity operator on the space of symmetric rank-2 tensors is \[I_{\mu\nu,\rho\sigma} = \frac{1}{2}\bigl(\eta_{\mu\rho}\eta_{\nu\sigma} + \eta_{\mu\sigma}\eta_{\nu\rho}\bigr).\]

It decomposes into four orthogonal spin projectors:

7.9.0.1 Spin-2 (transverse traceless):

\[\label{eq:Spin-2} P^{(2)}_{\mu\nu,\rho\sigma} = \frac{1}{2} \bigl(\theta_{\mu\alpha}\theta_{\nu\beta} + \theta_{\mu\beta}\theta_{\nu\alpha}\bigr) - \frac{1}{3}\,\theta_{\mu\nu}\,\theta_{\rho\sigma}.\tag{55}\]

7.9.0.2 Spin-1 (transverse vector):

\[P^{(1)}_{\mu\nu,\rho\sigma} = \frac{1}{2} \bigl(\theta_{\mu\alpha}\omega_{\nu\beta} + \theta_{\mu\beta}\omega_{\nu\alpha} + \theta_{\nu\alpha}\omega_{\mu\beta} + \theta_{\nu\beta}\omega_{\mu\alpha}\bigr).\]

7.9.0.3 Spin-0, scalar (\(s\)):

\[\label{eq:spin-0} P^{(0s)}_{\mu\nu,\rho\sigma} = \frac{1}{3}\,\theta_{\mu\nu}\,\theta_{\rho\sigma}.\tag{56}\]

7.9.0.4 Spin-0, scalar (\(w\)):

\[\label{eq:spin0}P^{(0w)}_{\mu\nu|\rho\sigma} = \omega_{\mu\nu}\,\omega_{\rho\sigma}.\tag{57}\]

7.9.1 Completeness and Orthogonality↩︎

\[\begin{align} I_{\mu\nu|\rho\sigma} &= P^{(2)}_{\mu\nu|\rho\sigma} + P^{(1)}_{\mu\nu|\rho\sigma} + P^{(0s)}_{\mu\nu|\rho\sigma} + P^{(0w)}_{\mu\nu|\rho\sigma}, \\[4pt] P^{(i)}_{\mu\nu|}{}^{\rho\sigma}\,P^{(j)}_{\rho\sigma|\kappa\lambda} &= \delta^{ij}\,P^{(i)}_{\mu\nu|\kappa\lambda}. \end{align}\]

The degrees of freedom count (\(d=4\)) is \(5 + 3 + 1 + 1 = 10\), consistent with a symmetric rank-2 tensor.

7.9.2 Off-Diagonal Scalar Mixer↩︎

When diagonalising the scalar sector it is useful to define the transition operator \[P^{(0sw)}_{\mu\nu|\rho\sigma} = \frac{1}{\sqrt{3}}\,\theta_{\mu\nu}\,\omega_{\rho\sigma},\] satisfying \(P^{(0sw)} = \bigl(P^{(0ws)}\bigr)^\dagger\) and \(P^{(0s)}\,P^{(0sw)} = P^{(0sw)}\), \(P^{(0sw)}\,P^{(0w)} = P^{(0sw)}\).

---

1. \(\alpha\geq0\) throughout↩︎

2. This is equivalent to adding the term \(-\beta g^{\mu \nu}\) to the deWitt metric and setting all scalar potentials apart from \(\psi,\phi\) to zero↩︎

3. \(G_{\rho\sigma}\tilde{h}^{\rho\sigma}=4(\tilde{\Phi}\triangle\psi+\tilde{\psi}((\triangle-3\partial_0^2)\psi-\triangle\Phi)\)↩︎

4. The equations of motion also clearly differ from the Einstein equations for the scalar sector↩︎

5. without introducing additional fields as in the Stueckelberg trick↩︎