1 Introduction

The recent experimental discovery of gravitational waves (GWs) [1] has marked a new era for both observational and theoretical physics. With the new coming data from LIGO and from future experiments like LISA, it will become possible to test modified gravity theories, establishing for which range of parameters these theories agree with observations. Particularly, it may be even possible to test quantum gravity in its low-energy limit, even though a complete quantum theory for gravity remains one of the greatest problems in modern physics.

A natural observable to consider is the GW energy. As a non-linear phenomenon, gravity couples to itself and thus gravitates, which means that GWs—being a manifestation of gravity—produce a backreaction into the spacetime. Hence, one should be able to find a stress-energy tensor for the GWs that accounts for this phenomenon. In the case of classical general relativity (GR), such a stress-energy tensor is known:

$$\begin{aligned} t_{\mu \nu }^{\text {GR}} = \frac{1}{32\pi G} \langle \partial _\mu h_{\alpha \beta }\partial _\nu h^{\alpha \beta } \rangle , \end{aligned}$$
(1)

where the brackets denote an average over spacetime, which is responsible for taking only the long-wavelength modes; its precise definition will be explained later on. The GW stress-energy tensor has also been calculated for some other theories, including f(R), Chern–Simons and higher-derivative gravity [2,3,4,5]. In [6], it was indicated how the parameters of an analytic f(R) theory could be constrained by the measurement of the energy or momentum carried away by the GWs.

The phenomenology, however, is not the only motivation. An alternative for dark energy has been proposed based on the effective stress-energy tensor [4, 7,8,9]. Although this is not possible in GR because of the vanishing trace of \(t_{\mu \nu }^{\text {GR}}\), it was pointed out it could be possible in modified gravity theories. However, it was also found that in some models such as Starobinsky gravity, the effective stress-energy tensor could not be the only factor as it does not produce the right value for the cosmological constant [4]. We will show that the large contributions from the Standard Model cannot be canceled by the quantum gravitational effects, thus requiring the existence of another mechanism able to reconcile the discrepancy between theory and observation.

The purpose of this paper is, then, twofold: we will establish new phenomenological bounds and discuss the possibility of generating a contribution to the cosmological constant in this framework. Effective field theory techniques will be used to calculate quantum contributions to the GW backreaction and to the wave equation in an arbitrary background. The short-wave formalism will be employed, consisting of an averaging procedure that separates the low-frequency modes from the high ones, in order to calculate the GW stress-energy tensor in quantum GR. These theoretical findings will be useful to constrain some of the parameters of effective quantum gravity by the direct comparison with LIGO’s observations. Furthermore, on the theoretical side, they give us new insights of gravity at the quantum level since this approach is model independent and, as such, leads to genuine predictions of quantum gravity.

This paper is organized as follows. In Sect. 2, we will review the main results of the effective field theory approach applied to gravity. In Sect. 3, we use the short-wave formalism to calculate the leading order quantum corrections to the GW stress-energy tensor. The result allows us to constraint the amplitude of the massive mode present in effective quantum gravity. In Sect. 4, we discuss the quantum corrections to the propagation of GWs and we show that the equation describing the propagation in curved spacetime can be obtained by performing a minimal coupling prescription to the equation in Minkowski space. We draw the conclusions in Sect. 5.

2 Effective quantum gravity

The quantum effective action of gravity up to quadratic order in curvature is given by [10]

$$\begin{aligned} \Gamma&= \int \mathrm {d}^4x\sqrt{-g}\bigg (\frac{M_p^2}{2}R + b_1 R^2 + b_2 R_{\mu \nu }R^{\mu \nu } \nonumber \\&\quad + c_1 R\log \frac{\Box }{\mu ^2}R + c_2 R_{\mu \nu }\log \frac{\Box }{\mu ^2}R^{\mu \nu }\nonumber \\&\quad + c_3 R_{\mu \nu \rho \sigma }\log \frac{\Box }{\mu ^2}R^{\mu \nu \rho \sigma }\bigg ), \end{aligned}$$
(2)

where \(M_p = (8\pi G)^{-1/2}\) is the reduced Planck mass, G is the Newton’s constant, \(\mu \) is the renormalization scale and the kernel R denotes the Riemann tensor and its contractions (Ricci tensor and Ricci scalar) depending on the number of indices it carries. The signature \((-+++)\) will be adopted. We set the bare cosmological constant to zero as it is not important to our considerations. The coefficients \(b_i\) are free parameters and must be fixed by observations, while the coefficients \(c_i\) are predictions of the infra-red theory and depend on the field content under consideration (see Table 1 in [10] for their precise values). The log operators are known to lead to acausal effects that need to be removed by resolving the non-local operator as

$$\begin{aligned} \log \frac{\Box }{\mu ^2} = \int _0^\infty \mathrm {d}s\left( \frac{1}{\mu ^2+s} - G(x,x',\sqrt{s})\right) , \end{aligned}$$
(3)

where \(G(x,x';\sqrt{s})\) is a Green’s function for

$$\begin{aligned} (\Box + k^2)G(x,x';k) = \delta ^4(x-x'), \end{aligned}$$
(4)

and imposing proper boundary conditions on \(G(x,x';k)\) so that the result respects causality. Moreover, in the weak field limit, the log terms are not independent due to the following relation (see [3]):

$$\begin{aligned}&\delta \int \mathrm {d}^4x\sqrt{-g}\left( R_{\mu \nu \rho \sigma }\log \frac{\Box }{\mu ^2}R^{\mu \nu \rho \sigma }-4 R_{\mu \nu }\log \frac{\Box }{\mu ^2}R^{\mu \nu } \right. \nonumber \\&\qquad \left. +\, R\log \frac{\Box }{\mu ^2}R\right) {\mathop {=}\limits ^{\text {weak}}} 0. \end{aligned}$$
(5)

This can also be seen by linearizing the field equations [11]. The log operators in the above expression certainly break the topological invariance given by the Gauss–Bonnet theorem. Nonetheless, such expression still provides a useful relation that can be used to simplify calculations in the weak field limit. Therefore, since we will be interested only in the weak field scenario, the last term in (2) will be eliminated in favour of the other two log terms, which translates into a shift of their coefficients:

$$\begin{aligned} c_1&\rightarrow \alpha \equiv c_1-c_3, \end{aligned}$$
(6)
$$\begin{aligned} c_2&\rightarrow \beta \equiv c_2+4 c_3. \end{aligned}$$
(7)

Hence, from now on, \(\alpha \) will denote the coefficient of \(R\log \frac{\Box }{\mu ^2}R\) and \(\beta \) the coefficient of \(R_{\mu \nu }\log \frac{\Box }{\mu ^2}R^{\mu \nu }\). Note, however, that the last term in (2) will give independent contributions in the non-linear regime and, in particular, the background equations of motion [left-hand side of (20) below] will be changed, but none of this affects the right-hand side of (20).

The quantum action (2) yields the equations of motion (EOM)

$$\begin{aligned} G_{\mu \nu } + \Delta G_{\mu \nu }^L + \Delta G_{\mu \nu }^{NL} = 8\pi G T_{\mu \nu }, \end{aligned}$$
(8)

where \(\Delta G_{\mu \nu }^L\) denotes the local contribution to the modification of Einstein’s tensor and \(\Delta G_{\mu \nu }^{NL} = \Delta G_{\mu \nu }^\alpha + \Delta G_{\mu \nu }^\beta \) is the non-local one (due to the log operator), coming from the terms proportional to \(\alpha \) and \(\beta \), denoted by \(\Delta G_{\mu \nu }^\alpha \) and \(\Delta G_{\mu \nu }^\beta \), respectively. Here we will show only the calculation of the non-local part \(\Delta G_{\mu \nu }^{NL}\) as the local contribution can be straightforwardly obtained from it. However, our final results will be completely general, including both local and non-local physics. \(\Delta G_{\mu \nu }^\alpha \) has been calculated in the literature [12]:

$$\begin{aligned} -\xi \Delta G_{\mu \nu }^\alpha= & {} 2\left( R_{\mu \nu }-\frac{1}{4}g_{\mu \nu }R\right) \left( \log \frac{\square }{\mu ^{2}}R\right) \nonumber \\&-\,2 (\nabla _{\mu }\nabla _{\nu }-g_{\mu \nu }\square ) \left( \log \frac{\square }{\mu ^{2}}R\right) , \end{aligned}$$
(9)

where \(\xi = \frac{1}{16\pi G\alpha }\). Note that the integral term appearing in [12], which comes from the variation of the D’Alembert operator, is not present here. This is because in the weak field limit the variation of the D’Alembert operator leads to negligible contributions [13]. The other contribution to \(\Delta G_{\mu \nu }\) is given by

$$\begin{aligned} \zeta \Delta G_{\mu \nu }^\beta&= -\,\frac{1}{2} g_{\mu \nu }R_{\rho \sigma }\log \left( \frac{\Box }{\mu ^2}\right) R^{\rho \sigma } + \Box \log \left( \frac{\Box }{\mu ^2}\right) R_{\mu \nu }\nonumber \\&\quad + g_{\mu \nu }\nabla _\rho \nabla _\sigma \log \left( \frac{\Box }{\mu ^2}\right) R^{\rho \sigma }\nonumber \\&\quad + R_\mu ^\sigma \log \left( \frac{\Box }{\mu ^2}\right) R_{\nu \sigma } + R_\nu ^\sigma \log \left( \frac{\Box }{\mu ^2}\right) R_{\mu \sigma }\nonumber \\&\quad - \nabla _\rho \nabla _\mu \log \left( \frac{\Box }{\mu ^2}\right) R^{\rho }_{\nu } - \nabla _\rho \nabla _\nu \log \left( \frac{\Box }{\mu ^2}\right) R_{\mu }^{\rho } \end{aligned}$$
(10)

where \(\zeta = \frac{1}{16\pi G\beta }\).

3 Gravitational wave backreaction

The first step is to separate the fluctuations \(h_{\mu \nu }\) (GWs) from the background geometry \(\bar{g}_{\mu \nu }\), via \(g_{\mu \nu } = \bar{g}_{\mu \nu } + h_{\mu \nu }\). This separation is only meaningful in the limit where the GW wavelength \(\lambda \) is much smaller than the background radius L, i.e. \(\lambda \ll L\), so that a clear distinction between background and GW can be made. As a first approximation, the background metric \(\bar{g}_{\mu \nu }\) will be used to raise/lower indices as well as to build all the operators, e.g. \(\Box = \bar{g}^{\mu \nu }\nabla _\mu \nabla _\nu \). The connection is also assumed to be compatible with \(\bar{g}_{\mu \nu }\) instead of \(g_{\mu \nu }\).

The separation of gravity into background and fluctuations allows one to expand the Ricci tensor as

$$\begin{aligned} R_{\mu \nu } = \bar{R}_{\mu \nu } + R_{\mu \nu }^{(1)} + R_{\mu \nu }^{(2)} + O(h^3), \end{aligned}$$
(11)

where the bar quantities are calculated with respect to the background and the rest depends only on the fluctuation. The superscript (n) is used to indicate the order in h of the underlying term. Naively, one could think that the EOM could be calculated order by order in h, giving no backreaction into the background. The problem is that there are two small parameters in the game, namely the fluctuation amplitude h and \(\varepsilon \equiv \frac{\lambda }{L}\), so that one can compensate the other. Their relation is fixed by the EOMFootnote 1 and in the presence of external matter

$$\begin{aligned} h\ll \varepsilon \ll 1, \end{aligned}$$
(12)

as can be seen from Eq. (8).

To obtain the GW backreaction, one then needs to calculate the average of tensor fields over a region of length scale d, where \(\lambda \ll d\ll L\). This makes the high-frequency modes go away due to their rapid oscillation, but they leave the low modes intact. The subtle point is that there is no canonical way of summing tensors based on different points of a manifold. Here Isaacson’s definition [14, 15] of the average of a tensor will be used, which is based on the idea of parallel transporting the tensor field along geodesics from each spacetime position to a common point where its different values can be compared:

$$\begin{aligned} \langle A_{\mu \nu }(x) \rangle= & {} \int j_\mu ^{\alpha '}(x,x') j_\nu ^{\beta '}(x,x') A_{\alpha '\beta '}(x')\nonumber \\&f(x,x')\sqrt{-\bar{g}(x')}\mathrm{d}^4x', \end{aligned}$$
(13)

where \(j_\mu ^{\alpha '}\) is the bivector of geodesic parallel displacement and \(f(x,x')\) is a weight function that falls quickly and smoothly to zero when \(|x-x'|>d\), such that

$$\begin{aligned} \int _\text {all space} f(x,x')\sqrt{-\bar{g}(x')}\mathrm{d}^4x'=1. \end{aligned}$$
(14)

From this definition, the following rules can be proven [2]:

  • The average of an odd product of short-wavelength quantities vanishes.

  • The derivative of a short-wavelength tensor averages to zero, e.g., \( \langle \nabla _\mu T_{\alpha \beta }^\mu \rangle =0\).

  • As a corollary, integration by parts can be performed and one can flip derivatives, e.g., \( \langle R_\alpha ^\mu \nabla _\mu S_\beta \rangle = - \langle S_\beta \nabla _\mu R_\alpha ^\mu \rangle \).

Therefore, to obtain the backreaction one has to calculate

$$\begin{aligned} \left\langle G_{\mu \nu }\right\rangle + \left\langle \Delta G_{\mu \nu }^{NL}\right\rangle = 8\pi G \left\langle T_{\mu \nu }\right\rangle \end{aligned}$$
(15)

up to second order in h (higher orders are extremely small)Footnote 2. Taking the average of Eq. (9) gives

$$\begin{aligned} -\xi \left\langle \Delta G_{\mu \nu }^\alpha \right\rangle&= 2\left( \left\langle R_{\mu \nu }\log \left( \frac{\Box }{\mu ^2}\right) R\right\rangle \right. \nonumber \\&\quad \left. -\,\frac{1}{4} \bar{g}_{\mu \nu }\left\langle R\log \left( \frac{\Box }{\mu ^2}\right) R\right\rangle \right) \nonumber \\&\quad -2\left\langle (\nabla _\mu \nabla _\nu -g_{\mu \nu }\Box )\log \left( \frac{\Box }{\mu ^2}\right) R\right\rangle . \end{aligned}$$
(16)

It follows from the rules that

$$\begin{aligned} \left\langle R_{\mu \nu }\log \left( \frac{\Box }{\mu ^2}\right) R^{\mu \nu }\right\rangle= & {} \bar{R}_{\mu \nu }\log \left( \frac{\Box }{\mu ^2}\right) \bar{R}^{\mu \nu }\nonumber \\&+\,\left\langle R_{\mu \nu }^{(1)}\log \left( \frac{\Box }{\mu ^2}\right) R^{(1)\mu \nu }\right\rangle , \end{aligned}$$
(17)

since the average of linear terms in h vanishes. Cross terms (e.g. \(\bar{R} R^{(2)}\)) are absent as they are negligible due to the condition (12). In addition, the last line of Eq. (16) has a global derivative so that the high-frequency contribution averages to zero.

The combination of Eqs. (16) and (17) results in

$$\begin{aligned} -\,\xi \langle \Delta G_{\mu \nu }^\alpha \rangle&= 2\left( \bar{R}_{\mu \nu } -\frac{1}{4} \bar{g}_{\mu \nu }\bar{R}\right) \log \left( \frac{\Box }{\mu ^2}\right) \bar{R} \nonumber \\&\quad + 2\left( \left\langle R_{\mu \nu }^{(1)}\log \left( \frac{\Box }{\mu ^2}\right) R^{(1)}\right\rangle \right. \nonumber \\&\quad \left. -\, \frac{1}{4} \bar{g}_{\mu \nu } \left\langle R^{(1)}\log \left( \frac{\Box }{\mu ^2}\right) R^{(1)}\right\rangle \right) \nonumber \\&\quad -2(\nabla _\mu \nabla _\nu -\bar{g}_{\mu \nu }\Box )\log \left( \frac{\Box }{\mu ^2}\right) \bar{R}. \end{aligned}$$
(18)

Similarly, taking the average of Eq. (10) gives

$$\begin{aligned} \zeta \langle \Delta G_{\mu \nu }^\beta \rangle&= -\,\frac{1}{2} \bar{g}_{\mu \nu }\left( \bar{R}_{\rho \sigma }\log \left( \frac{\Box }{\mu ^2}\right) \bar{R}^{\rho \sigma } \right. \nonumber \\&\quad \left. +\, \left\langle R_{\rho \sigma }^{(1)}\log \left( \frac{\Box }{\mu ^2}\right) R^{(1)\rho \sigma }\right\rangle \right) \nonumber \\&\quad + \Box \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}_{\mu \nu } + \bar{g}_{\mu \nu }\nabla _\rho \nabla _\sigma \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}^{\rho \sigma }\nonumber \\&\quad + \bar{R}_\mu ^\sigma \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}_{\nu \sigma } + \bar{R}_\nu ^\sigma \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}_{\mu \sigma } \nonumber \\&\quad + 2\left\langle R_\mu ^{(1)\sigma }\log \left( \frac{\Box }{\mu ^2}\right) R_{\nu \sigma }^{(1)}\right\rangle \nonumber \\&\quad - \nabla _\rho \nabla _\mu \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}^{\rho }_\nu - \nabla _\rho \nabla _\nu \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}_\mu ^{\rho }. \end{aligned}$$
(19)

Combining Eqs. (15), (18) and (19) leads to the background EOM

$$\begin{aligned}&\bar{R}_{\mu \nu }-\frac{1}{2} \bar{g}_{\mu \nu }\bar{R} -\frac{2}{\xi }\left[ \left( \bar{R}_{\mu \nu }-\frac{1}{4} \bar{g}_{\mu \nu }\bar{R}\right) \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}\right. \nonumber \\&\qquad \left. -\, (\nabla _\mu \nabla _\nu -\bar{g}_{\mu \nu }\Box )\log \left( \frac{\Box }{\mu ^2}\right) \bar{R}\right] \nonumber \\&\qquad -\frac{1}{2\zeta } \bar{g}_{\mu \nu }\bar{R}_{\rho \sigma }\log \left( \frac{\Box }{\mu ^2}\right) \bar{R}^{\rho \sigma } + \frac{1}{\zeta }\bar{R}_\mu ^\sigma \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}_{\nu \sigma } \nonumber \\&\qquad + \frac{1}{\zeta }\bar{R}_\nu ^\sigma \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}_{\mu \sigma } + \frac{1}{\zeta }\Box \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}_{\mu \nu }\nonumber \\&\qquad + \frac{1}{\zeta }\bar{g}_{\mu \nu }\nabla _\rho \nabla _\sigma \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}^{\rho \sigma } - \frac{1}{\zeta }\nabla _\rho \nabla _\mu \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}^{\rho }_\nu \nonumber \\&\qquad - \frac{1}{\zeta }\nabla _\rho \nabla _\nu \log \left( \frac{\Box }{\mu ^2}\right) \bar{R}_\mu ^{\rho }\nonumber \\&\quad = 8\pi G( \langle T_{\mu \nu } \rangle + t_{\mu \nu }^{GR} + t_{\mu \nu }^{NL}), \end{aligned}$$
(20)

where \(t_{\mu \nu }^{GR}\) is the classical contribution to the GW stress-energy tensor:

$$\begin{aligned} t_{\mu \nu }^{GR} = -\,\frac{1}{8\pi G}\left( \langle R_{\mu \nu }^{(2)} \rangle -\frac{1}{2}\bar{g}_{\mu \nu }\left\langle R^{(2)}\right\rangle \right) , \end{aligned}$$
(21)

and \(t_{\mu \nu }^{NL}\) is the non-local one:

$$\begin{aligned} t_{\mu \nu }^{NL}&= -\,\frac{1}{8\pi G}\Bigg [-\frac{2}{\xi }\left( \left\langle R_{\mu \nu }^{(1)}\log \left( \frac{\Box }{\mu ^2}\right) R^{(1)}\right\rangle \right. \nonumber \\&\quad \left. -\,\frac{1}{4}\bar{g}_{\mu \nu }\left\langle R^{(1)}\log \left( \frac{\Box }{\mu ^2}\right) R^{(1)}\right\rangle \right) \nonumber \\&\quad +\frac{2}{\zeta }\left\langle R_\mu ^{(1)\sigma }\log \left( \frac{\Box }{\mu ^2}\right) R_{\nu \sigma }^{(1)}\right\rangle \nonumber \\&\quad -\frac{1}{2\zeta }\bar{g}_{\mu \nu }\left\langle R_{\rho \sigma }^{(1)}\log \left( \frac{\Box }{\mu ^2}\right) R^{(1)\rho \sigma }\right\rangle \Bigg ]. \end{aligned}$$
(22)

Similarly, the local contribution is given by

$$\begin{aligned} t_{\mu \nu }^L&= -\,\frac{1}{8\pi G}\Bigg [-32\pi G b_1\left( \langle R_{\mu \nu }^{(1)}R^{(1)} \rangle -\frac{1}{4}\bar{g}_{\mu \nu }\left\langle R^{(1)}R^{(1)}\right\rangle \right) \nonumber \\&\quad +32\pi G b_2\left\langle R_\mu ^{(1)\sigma }R_{\nu \sigma }^{(1)}\right\rangle -8\pi G b_2\bar{g}_{\mu \nu } \langle R_{\rho \sigma }^{(1)}R^{(1)\rho \sigma } \rangle \Bigg ]. \end{aligned}$$
(23)

Therefore, the total GW stress-energy tensor is \(t_{\mu \nu } = t_{\mu \nu }^{GR} + t_{\mu \nu }^{L} + t_{\mu \nu }^{NL}\).

At this point some comments are necessary. First of all, observe that the left-hand side of Eq. (20) corresponds solely to the background effect, which we interpret as pure gravity. In fact, the left-hand side is exactly the same as in Eq. (8) when \(g_{\mu \nu }\) is replaced by \(\bar{g}_{\mu \nu }\). The right-hand side represents the matter sector, as usual, but with the inclusion of the GW contribution. Such a contribution represents the most general stress-energy tensor to leading order, accounting for both classical and quantum effects. Note that, due to the diffeomorphism invariance of the theory, the total energy-momentum tensor is covariantly conserved,

$$\begin{aligned} \nabla ^\mu (T_{\mu \nu }+t_{\mu \nu })=0, \end{aligned}$$
(24)

which shows that energy and momentum are exchanged between matter sources and GWs. Far away from the source, this gives the conservation of the GW energy-momentum tensor

$$\begin{aligned} \partial ^\mu t_{\mu \nu } = 0. \end{aligned}$$
(25)

Up to this point, no gauge conditions have been applied and \(t_{\mu \nu }\) also accounts for spurious degrees of freedom. To eliminate them, we shall take the limit where the GW is far away from the source, so that the background is nearly flat and the linear EOM becomes [11]

$$\begin{aligned} \Box _\eta h_{\mu \nu } + 16\pi G\left[ b_2 + \beta \log \left( \frac{\Box _\eta }{\mu ^2}\right) \right] \Box _\eta ^2 h_{\mu \nu } = 0, \end{aligned}$$
(26)

where \(\Box _\eta =\eta ^{\mu \nu }\partial _\mu \partial _\nu \) is the flat D’Alembert operator. Note the absence of the parameter \(\alpha \) in Eq. (26). This happens because \(\alpha \) is proportional to terms depending on the trace h, which can be taken as zero far away from the source. Using the gauge conditions \(\partial ^\nu h_{\mu \nu }=0\) and \(h=0\) (only valid outside the source) together with Eq. (26) in the definition of \(t_{\mu \nu }\) gives

$$\begin{aligned} t_{\mu \nu }= & {} \frac{1}{8\pi G}\left[ \frac{1}{4} \langle \partial _\mu h_{\alpha \beta }\partial _\nu h^{\alpha \beta } \rangle + \frac{1}{2} \langle h_\mu ^\sigma \Box _\eta h_{\nu \sigma } \rangle \right. \nonumber \\&\left. -\, \frac{1}{8} \eta _{\mu \nu } \langle h_{\rho \sigma }\Box _\eta h^{\rho \sigma } \rangle \right] . \end{aligned}$$
(27)

Comparing this to Eq. (1), it is clearly seen that the first term in \(t_{\mu \nu }\) corresponds to GR, while the other two come from quantum corrections. Observe that the log operators do not appear explicitly in Eq. (27) as the gravitational field is on shell. This means that their contribution will only show up in the field solutions. For the same reason, the procedure (3) of imposing causality need not be pursued at this stage as the non-local effects are only reflected in the solutions for \(h_{\mu \nu }\). The parameters \(b_2\) and \(\beta \) now only appear in the mass m of \(h_{\mu \nu }\).

The GW energy density is then

$$\begin{aligned} \rho \equiv t_{00}= & {} \frac{1}{8\pi G}\left[ \frac{1}{4} \langle \dot{h}_{\alpha \beta }\dot{h}^{\alpha \beta } \rangle + \frac{1}{2} \langle h_0^\alpha \Box _\eta h_{0\alpha } \rangle \right. \nonumber \\&\left. +\, \frac{1}{8} \langle h_{\rho \sigma }\Box _\eta h^{\rho \sigma } \rangle \right] . \end{aligned}$$
(28)

As a concrete example, take a plane wave solution propagating in the z direction,

$$\begin{aligned} h_{\mu \nu } = \epsilon _{\mu \nu }\cos (\omega t - k z). \end{aligned}$$
(29)

Plugging this into Eq. (28) gives

$$\begin{aligned} \rho = \frac{1}{16\pi G}\left[ \frac{\epsilon ^2\omega ^2}{4} + \frac{1}{2}\left( \epsilon _0^\alpha \epsilon _{0\alpha } + \frac{\epsilon ^2}{4}\right) (\omega ^2-k^2)\right] . \end{aligned}$$
(30)

Therefore, modifications in the dispersion relations lead to measurable differences into the GW energy. In the case of the classical wave, i.e. \(\omega ^2=k^2\), the second term vanishes identically, resulting in the classical energy as expected. In the most general case, there could be complex frequencies leading to damping as was shown in [17,18,19,20]. In such case, Eq. (30) can be straightforwardly generalized. Note that the second term in (30) is proportional to the particle’s mass m and, therefore, is constant as any change in the frequency would be compensated by a change in the momentum. Dividing the constant term by the critical density \(\rho _c=\frac{3H_0^2}{8\pi G}\), where \(H_0\) is today’s Hubble constant, leads to the frequency-independent gravitational wave density parameter \(\Omega _0\), which was constrained to be smaller than \(1.7\times 10^{-7}\) by LIGO [21]:

$$\begin{aligned} \Omega _0 = \frac{1}{12}\left( \epsilon _0^\alpha \epsilon _{0\alpha } + \frac{\epsilon ^2}{4}\right) \frac{m^2}{H_0^2} < 1.7\times 10^{-7}. \end{aligned}$$
(31)

We remind the reader that the initial parameters \(b_2\) and \(\beta \) appear only in terms of the mass m as the field \(h_{\mu \nu }\) is on shell. Figure 1 shows the allowed region of the parameter space \((m,\epsilon )\).

Fig. 1
figure 1

The blue area in the graph represents the allowed region of the parameter space \((m,\epsilon )\)

Full size image

Using the lower bound on the mass of the complex poleFootnote 3 found in [17], i.e. \(m>5\times 10^{-13}\) GeV, we can translate the above constraint to

$$\begin{aligned} \epsilon < 1.4\times 10^{-33}. \end{aligned}$$
(32)

To the best of our knowledge, this is the first bound ever found on the amplitude of the massive mode. It is 12 orders of magnitude smaller than the strain sensibility of LIGO’s interferometer, which can probe amplitudes up to \(\sim 10^{-22}\) in the frequency range from 10 Hz to 10 kHz. Although it seems hopeless to reach such small distances, the Chongqing University detector (currently under development) will be able to probe amplitudes as small as \(10^{-32}\) [22] in the high-frequency range 0.1–10 GHz, which is not far from the bound just found. Observe, however, that we have found an upper bound on \(\epsilon \) and not a lower one, thus \(\epsilon \) could be arbitrarily small and not be detectable by the Chongqing detector. Should the existence of these extra modes be confirmed in future experiments, this would be the first evidence for a massive mode.

As it was stressed before, the effective energy-momentum tensor may lead to a contribution to the accelerated expansion of today’s universe if its trace is different from zero. The trace of the GW energy-momentum tensor (27) is non-vanishing:

$$\begin{aligned} t = -\,\frac{1}{32\pi G} \langle h_{\alpha \beta }\Box _\eta h^{\alpha \beta } \rangle \ne 0, \end{aligned}$$
(33)

as the gravitational field now satisfies the modified EOM (26). Therefore, the energy-momentum tensor \(t_{\mu \nu }\) can be split into a traceful and a traceless component,

$$\begin{aligned} t_{\mu \nu } = t_{\mu \nu }-\frac{1}{4} \eta _{\mu \nu }t + \frac{1}{4} \eta _{\mu \nu }t, \end{aligned}$$
(34)

and the cosmological constant can be identified as

$$\begin{aligned} \Lambda \equiv \frac{1}{16} \langle h_{\alpha \beta }\Box _\eta h^{\alpha \beta } \rangle = \frac{1}{16}\epsilon ^2 m^2, \end{aligned}$$
(35)

where in the second equality the plane wave solution (29) was used. After taking the average, \(\Lambda \) depends very mildly on space and time. In fact, it is precisely constant across any region of length d and its variation only becomes appreciable in a region containing several lengths of size d. Therefore, for our purposes, we can safely neglect the spacetime dependence of the emergent cosmological constant \(\Lambda \) and consider it a constant indeed. Remember that, initially, the cosmological constant was set to zero. A non-zero initial or bare cosmological constant \(\Lambda _b\) would just be shifted by the \(\Lambda \) found above and the physical cosmological constant would be \(\Lambda _\mathrm{eff}\equiv \Lambda _b + \Lambda \). The important proposition here is that quantum gravitational waves give a non-zero contribution to the cosmological constant \(\Lambda _\mathrm{eff}\). In this scenario, the new gravitational interactions and degrees of freedom appearing in high energies are represented by non-local effects in the low-energy limit. The latter, combined with the local interactions, yields a gravitational energy-momentum tensor whose trace is non-vanishing and which contributes to the total cosmological constant.

4 Gravitational wave propagation

Up to now, only the physics of the low-frequency waves has been considered. For completeness, we shall turn our attention to the high-frequency ones, which will lead to the equation describing the GW propagation in curved spacetime. This is easily achieved by subtracting the background equation (20) from the total EOM (8)

$$\begin{aligned} G_{\mu \nu } + \Delta G_{\mu \nu } - \langle G_{\mu \nu } + \Delta G_{\mu \nu } \rangle = 8\pi G(T_{\mu \nu } - \left\langle T_{\mu \nu }\right\rangle ), \end{aligned}$$
(36)

where \(\Delta G_{\mu \nu } = \Delta G_{\mu \nu }^L + \Delta G_{\mu \nu }^{NL}\). Ignoring the local part for a moment and keeping only the terms up to linear order in h and \(\lambda /L\) gives

$$\begin{aligned}&R_{\mu \nu }^{(1)} - \frac{1}{2} \bar{g}_{\mu \nu }R^{(1)} + \frac{2}{\xi }(\nabla _\mu \nabla _\nu -\bar{g}_{\mu \nu }\Box )\log \left( \frac{\Box }{\mu ^2}\right) R^{(1)} \nonumber \\&\quad + \frac{1}{\zeta }\Bigg [ \Box \log \left( \frac{\Box }{\mu ^2}\right) R^{(1)}_{\mu \nu }\nonumber \\&\quad + \bar{g}_{\mu \nu }\nabla _\rho \nabla _\sigma \log \left( \frac{\Box }{\mu ^2}\right) R^{(1)\rho \sigma } - \nabla _\rho \nabla _\mu \log \left( \frac{\Box }{\mu ^2}\right) R^{(1)\rho }_\nu \nonumber \\&\quad - \nabla _\rho \nabla _\nu \log \left( \frac{\Box }{\mu ^2}\right) R^{(1)\rho }_\mu \Bigg ]=0. \end{aligned}$$
(37)

Outside the matter source, we can use the gauge \(\nabla ^\nu h_{\mu \nu }=0\) together with \(h=0\), leading to

$$\begin{aligned} \Box h_{\mu \nu } + 16\pi G\beta \log \left( \frac{\Box }{\mu ^2}\right) \Box ^2 h_{\mu \nu } = 0. \end{aligned}$$
(38)

Analogously, including the local terms gives

$$\begin{aligned} \Box h_{\mu \nu } + 16\pi G\left[ b_2 + \beta \log \left( \frac{\Box }{\mu ^2}\right) \right] \Box ^2 h_{\mu \nu } = 0. \end{aligned}$$
(39)

When deriving Eq. (39), we made use of the commutation relation of covariant derivatives and we discarded terms proportional to the background curvature as they only contribute to higher orders in \(\lambda /L\). Equation (39) describes the propagation of GWs in an arbitrary curved background in the absence of external matter, when the only source for a non-vanishing Ricci tensor is the GW energy-momentum tensor. The curvature terms do not appear as they provide no contribution to leading order. Therefore, the case where curvature is present can be obtained by applying a simple “minimal coupling” prescription to Eq. (26) where spacetime is flat, that is, by performing the following substitution:

$$\begin{aligned} \eta _{\mu \nu }\rightarrow \bar{g}_{\mu \nu },\end{aligned}$$
(40)
$$\begin{aligned} \partial _\mu \rightarrow \nabla _\mu . \end{aligned}$$
(41)

Equations (20) and (39) together describe the entire classical and quantum process (to leading order) of the GW self-gravitation: small perturbations around spacetime change the curvature, which in turn modify the GW’s trajectory and vice versa.

5 Conclusions

We showed in this paper how to calculate the quantum corrections to the GW stress-energy tensor. The result shows that quantum effects promote the traceless tensor \(t_{\mu \nu }^{\text {GR}}\) to a traceful quantity that contributes to the current accelerated expansion of the universe. In addition, the energy density component acquires a dependence on modifications to the dispersion relation, indicating a useful observable to probe when looking for quantum gravitational effects. In fact, by using the latest LIGO’s data, it was obtained a new upper bound on the amplitude of the massive mode. We also showed that the high-frequency mode equation led to a generalization of the wave equation (26) to arbitrary curved spacetimes (39). Such a generalization is important to the study of quantum GW solutions in cosmology and in the early universe where the spacetime was highly curved. Lastly, it must be stressed once again that these quantum contributions are model independent (since they are derived from an effective field theory) and constitute actual predictions of quantum gravity, shedding new light on Quantum Gravity as a whole and giving us some hints of how a complete theory, if such a theory exists, should behave below the Planck scale.