FeynGrav and Recent Progress in Computational Perturbative Quantum Gravity

Abstract

This article reviews recent progress in computational quantum gravity caused by the framework that efficiently computes Feynman’s rules. The framework is implemented in the FeynGrav package, which extends the functionality of the widely used FeynCalc package. FeynGrav provides all the tools to study quantum gravitational effects within the standard model. We review the framework, provide the theoretical background for the efficient computation of Feynman rules, and present the proof of its completeness. We review the derivation of Feynman rules for general relativity, Horndeski gravity, Dirac fermions, Proca field, electromagnetic field, and $SU\left(N\right)$ Yang–Mills model. We conclude with a discussion of the current state of the FeynGrav package and discuss its further development.

1. Introduction

Perturbative quantum gravity is one of several approaches to quantum gravity. It exists in the effective field theory paradigm, and describes gravitational phenomena at energies below the Planck scale [1,2,3]. To a certain extent, the theory associates gravity with small metric perturbations propagating about the flat spacetime. The Planck scale arises naturally as the characteristic scale of small perturbations. Consequently, the theory can perturbatively treat perturbations with small amplitudes. However, the structure of the theory is more sophisticated and extends beyond this simple setup, as we discuss below.

Perturbation theory operates with propagating degrees of freedom, polarisation operators, and interaction operators describing the interaction between different degrees of freedom. On the practical ground, defining a theory means explicitly providing these components. The graviton propagator and the polarisation operators were derived and studied in classical papers [4,5,6]. The derivation of the interaction rules was the most challenging part of the theory. Because of the effective nature of the theory, it generates an infinite set of interaction operators suppressed by different powers of the same gravitational coupling. For a given level of the perturbation theory, only a finite number of operators contribute to a given matrix element. However, the number of operators grows with each perturbation order. Consequently, calculations become incredibly challenging without an algorithm describing how to obtain the interaction rules at any given order of perturbation theory.

The interaction rules for the first two orders of perturbation theory are published in the literature [7,8,9,10]. These rules were sufficient for one-loop calculations and allowed for significant advancement in perturbative quantum gravity [2,7,8,9,11,12,13,14,15,16,17,18]. Nonetheless, the general algorithm producing the interaction rules at a given order remained unknown. Several research directions were focused on creating such an algorithm. The most well-known ones proposed to use Ward–Takahashi identities [7,8,9] and to bootstrap gravity [19,20,21,22,23]. However, to our knowledge, they never resulted in a practical tool. An additional problem was the complexity of the interaction rules. The paper [9] derived interaction rules for three and four graviton vertices (see also [10], as the mentioned paper contained a misprint). In that particular parameterisation, the three graviton vertex contains 171 terms, while the four graviton vertex contains 2850. While it is exceptionally challenging to manually operate with interaction terms having hundreds and thousands of terms, the contemporary computational packages can (relatively) easily manipulate such expressions.

The recent publications [24,25] created the desirable algorithm. It was implemented in a computational package “FeynGrav”, created at the base of the widely used “FeynCalc” [26,27,28]. This article reviews the constructed theoretical framework, touching upon its application and perspective for further development. For the discussion of the theoretical framework, we focus on the original papers [24,25]. We briefly discuss how FeynGrav can significantly simplify the derivation of the well-known results [2,13,29,30,31] and is used to obtain new results [32,33,34,35]. We conclude with a discussion of the future development of the package.

The structure of this paper is as follows. Section 2 discusses the theoretical background behind the perturbative quantum gravity. Namely, we discuss the role of the perturbative quantum gravity together with the limits of its applicability. We discuss the theory in the context of renormalisation and argue, supporting the growing consensus, that its effective nature leaves no room for any “renormalisation issues”. Section 3 discusses the practical tools to develop the theoretical and computational framework. We discuss how to factorise an action describing gravity or a model coupled to gravity in a way most suitable for optimal computations. We separate some fundamental structures defining perturbative expansions and derive recursive relations that significantly improve their computation. Section 4 devoted to the derivation of the Feynman rules for general relativity, scalar field, Horndeski theory, Dirac fermions, Proca field, electromagnetic field, and $SU\left(N\right)$ Yang–Mills model. Section 5 discusses the FeynGrav package, its structure, and rules implemented in the existing version. We present a few explicit examples of perturbative calculations made with FeynGrav. Section 6 summarises the material presented in this paper and highlights perspective further development.

2. Perturbative Quantum Gravity

We divide the discussion of perturbative quantum gravity into two parts. The first part explores the theory’s mathematical and technical aspects, while the second focuses on its physical content. We start with the technical aspects, as they help to explain the premises behind its physical interpretation.

Perturbative quantum gravity is a quantum theory of small perturbation of metric $h_{\mu \nu }$ existing about the flat spacetime with the Minkowski metric ${\eta }_{\mu \nu }$. To preserve the canonical mass dimension of the field variable, one must introduce the gravitational coupling $\kappa$ with mass dimension $-1$. One defines the spacetime metric as a combination of the background metric and perturbations:

$$\begin{array}{c}\hfill g_{\mu \nu }={\eta }_{\mu \nu }+\kappa h_{\mu \nu }.\end{array}$$

(1)

The gravitational coupling $\kappa$ is associated with the Newton constant $G_{\mathrm{N}}$:

$$\begin{array}{c}\hfill {\kappa }^2=32\pi G_{\mathrm{N}}.\end{array}$$

(2)

It is important to highlight that the perturbative expansion (1) is finite in $\kappa$, and should not be interpreted as a truncation of an infinite series. Despite (1) being finite, it spawns an infinite series in $\kappa$ for the inverse metric:

$$\begin{array}{c}\hfill g^{\mu \nu }={\eta }^{\mu \nu }-\kappa h^{\mu \nu }+{\kappa }^2{h}^{\mu \sigma }{h}_{\sigma }{}^{\nu }+\mathcal{O}\left({\kappa }^3\right).\end{array}$$

(3)

Consequently, all quantities involving the inverse metric are infinite series in $\kappa$.

Let us look at the two most important examples. The first example is the Christoffel symbol. If the symbol has all lower indices, then it is a finite expression linear in $\kappa$:

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\alpha \mu \nu }=& \frac{1}{2}\left[{\partial }_{\mu }{g}_{\nu \alpha }+{\partial }_{\nu }{g}_{\mu \alpha }-{\partial }_{\alpha }{g}_{\mu \nu }\right]=\frac{\kappa }{2}\left[{\partial }_{\mu }{h}_{\nu \alpha }+{\partial }_{\nu }{h}_{\mu \alpha }-{\partial }_{\alpha }{h}_{\mu \nu }\right].\hfill \end{array}$$

(4)

On the contrary, the standard Christoffel symbol with a single upper index is an infinite series, since it involves the inverse metric:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\mu \nu }^{\alpha }=& \frac{1}{2}{g}^{\alpha \beta }\left[{\partial }_{\mu }{g}_{\nu \beta }+{\partial }_{\nu }{g}_{\mu \beta }-{\partial }_{\beta }{h}_{\mu \nu }\right]\hfill \\ \hfill =& \frac{\kappa }{2}\left[{\eta }^{\mu \nu }-\kappa h^{\mu \nu }+{\kappa }^2{h}^{\mu \sigma }{h}_{\sigma }{}^{\nu }+\mathcal{O}\left({\kappa }^3\right)\right]\left[{\partial }_{\mu }{h}_{\nu \beta }+{\partial }_{\nu }{h}_{\mu \beta }-{\partial }_{\beta }{h}_{\mu \nu }\right].\hfill \end{array}\end{array}$$

(5)

The second example is the volume factor $\sqrt{-g}$. The metric determinant itself is a finite expression in $\kappa$ (see [17]), but the volume factor is an infinite series because of the square root:

$$\begin{array}{c}\hfill \sqrt{-g}=1+\frac{\kappa }{2}{\eta }^{\mu \nu }{h}_{\mu \nu }-\frac{{\kappa }^2}{4}\left(h_{\mu \sigma }{h}^{\sigma }{}_{\nu }-\frac{1}{2}{\eta }^{\mu \nu }{\eta }^{\alpha \beta }{h}_{\mu \nu }{h}_{\alpha \beta }\right)+\mathcal{O}\left({\kappa }^3\right).\end{array}$$

(6)

The action of any gravity model is an infinite series in $\kappa$, since it includes the volume factor $\sqrt{-g}$ for the coordinate invariance. Moreover, a typical action of a gravity model includes invariants constructed from the Riemann tensor. Christoffel symbols constitute the Riemann tensor, introducing additional infinite series to the theory. From the technical point of view, this is the reason why the action of any gravity model contains an infinite number of interaction operators.

The path integral formalism provides a way to construct a perturbative quantum theory. The generating functional $\mathcal{Z}$ provides a way to calculate matrix elements:

$$\begin{array}{c}\hfill \mathcal{Z}=\int \mathcal{D}\left[g\right]exp\left[i\mathcal{A}\left[g\right]\right].\end{array}$$

(7)

Here, $\mathcal{A}$ is the microscopic action of the model. Within the perturbative approach, one uses (1) in the path integral, expanding it in an infinite series. Metric (1) performs a simple shift of integration variables by a constant value:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{Z}=& \int \mathcal{D}[\eta +\kappa h]exp\left[i\mathcal{A}[\eta +\kappa h]\right]\hfill \\ \hfill =& \int \mathcal{D}\left[h\right]exp[i\mathcal{A}\left[\eta \right]+i\frac{\delta \mathcal{A}}{\delta g_{\mu \nu }}{|}_{g=\eta }\kappa h_{\mu \nu }+i\frac{{\delta }^2\mathcal{A}}{\delta g_{\mu \nu }\delta g_{\alpha \beta }}{|}_{g=\eta }{\kappa }^2{h}_{\mu \nu }{h}_{\alpha \beta }\hfill \\ \hfill & +i\frac{{\delta }^3\mathcal{A}}{\delta g_{\mu \nu }\delta g_{\alpha \beta }\delta g_{\rho \sigma }}{|}_{g=\eta }{\kappa }^3{h}_{\mu \nu }{h}_{\alpha \beta }{g}_{\rho \sigma }+\mathcal{O}\left({\kappa }^4\right)].\hfill \end{array}\end{array}$$

(8)

The first term of this expansion is irrelevant. The multiplier does not depend on the integration variables $h_{\mu \nu }$ and factorises out of the path integral. In turn, while calculating a matrix element, the multiplier will be cancelled entirely because of the normalisation.

It is safe to assume that the second term shall also vanish. The term disappears if the flat spacetime makes the microscopic action stationary, which occurs when the classical field equations derived from the microscopic action $\mathcal{A}$ have flat spacetime as a solution. Consequently, the discussed approach only applies to some gravity models with further modifications. One can always extend it to the case of an arbitrary background spacetime [36,37]. Further, we will discuss the physical content of the theory and argue that it may not be necessary. Let us postpone the discussion, assume that the term vanishes, and proceed with the technical discussion.

The other terms of the expansion do not vanish, but describe the propagation of perturbations and their interactions. It is helpful to present $\mathcal{Z}$ in the following form, making its structure explicit. The action naturally contains ${\kappa }^{-2}$ factor, similar to the Einstein–Hilbert action. Therefore, the generating functional takes the following form:

$$\begin{array}{c}\hfill \mathcal{Z}=\int \mathcal{D}\left[h\right]exp\left[-\frac{i}{2}{h}_{\mu \nu }{\mathcal{O}}^{\mu \nu \alpha \beta }\square h_{\alpha \beta }+i\kappa {\widehat{\mathcal{V}}}_{\left(3\right)}^{{\mu }_1{\nu }_1{\mu }_2{\nu }_2{\mu }_3{\nu }_3}{h}_{{\mu }_1{\nu }_1}{h}_{{\mu }_2{\nu }_2}{h}_{{\mu }_3{\nu }_3}+\mathcal{O}\left({\kappa }^2\right)\right].\end{array}$$

(9)

Here, ${\mathcal{O}}^{\mu \nu \alpha \beta }$ is the differential operator defining the structure of the graviton propagator well-established in the literature [4,38]. Operator ${\widehat{\mathcal{V}}}_{\left(3\right)}$ is the differential operator defining the structure of interaction between three perturbations. The term $\mathcal{O}\left({\kappa }^2\right)$ includes all other interaction terms, each suppressed by a different power of the gravitational coupling.

The standard prescription provides a way to calculate a given matrix element. One shall introduce a formal external current $J^{\mu \nu }$ linearly coupled to the metric perturbation:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{Z}\left[J\right]=& \int \mathcal{D}\left[h\right]exp\left[i\mathcal{A}[\eta +\kappa h]+i{h}_{\mu \nu }{J}^{\mu \nu }\right]\hfill \\ \hfill =& \int \mathcal{D}\left[h\right]exp\left[-\frac{i}{2}{h}_{\mu \nu }{\mathcal{O}}^{\mu \nu \alpha \beta }\square h_{\alpha \beta }+i{h}_{\mu \nu }{J}^{\mu \nu }+i\kappa {\widehat{\mathcal{V}}}_{\left(3\right)}^{{\mu }_1{\nu }_1{\mu }_2{\nu }_2{\mu }_3{\nu }_3}{h}_{{\mu }_1{\nu }_1}{h}_{{\mu }_2{\nu }_2}{h}_{{\mu }_3{\nu }_3}+\mathcal{O}\left({\kappa }^2\right)\right]\hfill \end{array}\end{array}$$

(10)

The path integral reduces to the Gaussian integral, which can be solved explicitly:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{Z}\left[J\right]=& exp\left[i\kappa {\widehat{\mathcal{V}}}_{\left(3\right)}^{{\mu }_1{\nu }_1{\mu }_2{\nu }_2{\mu }_3{\nu }_3}\frac{\delta }{\delta J^{{\mu }_1{\nu }_1}}\frac{\delta }{\delta J^{{\mu }_2{\nu }_2}}\frac{\delta }{\delta J^{{\mu }_3{\nu }_3}}+\mathcal{O}\left({\kappa }^2\right)\right]exp\left[\frac{i}{2}{J}^{\mu \nu }{\mathcal{O}}_{\mu \nu \alpha \beta }^{-1}\square J^{\alpha \beta }\right].\hfill \end{array}\end{array}$$

(11)

Finally, each matrix element links to a variational derivative of the functional:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \left.〈0\right|h_{{\mu }_1{\nu }_1}\cdots h_{{\mu }_N{\nu }_N}\left|0\right.〉=\frac{1}{\mathcal{Z}}\int \mathcal{D}\left[h\right]h_{\mu \nu }\cdots h_{{\mu }_N{\nu }_N}exp\left[i\mathcal{A}\right]\hfill \\ \hfill & =\frac{1}{\mathcal{Z}\left[J\right]}\int \mathcal{D}\left[h\right]h_{\mu \nu }\cdots h_{{\mu }_N{\nu }_N}exp\left[i\mathcal{A}+i{h}_{\rho \sigma }{J}^{\rho \sigma }\right]{|}_{J=0}\hfill \\ \hfill & =\frac{\frac{\delta }{\delta J^{{\mu }_1{\nu }_1}}\cdots \frac{\delta }{\delta J^{{\mu }_N{\nu }_N}}exp\left[i\kappa {\widehat{\mathcal{V}}}_{\left(3\right)}^{{\rho }_1{\sigma }_1{\rho }_2{\sigma }_2{\rho }_3{\sigma }_3}\frac{\delta }{\delta J^{{\rho }_1{\sigma }_1}}\frac{\delta }{\delta J^{{\rho }_2{\sigma }_2}}\frac{\delta }{\delta J^{{\rho }_3{\sigma }_3}}+\mathcal{O}\left({\kappa }^2\right)\right]exp\left[\frac{i}{2}{J}^{\mu \nu }{\mathcal{O}}_{\mu \nu \alpha \beta }^{-1}\square J^{\alpha \beta }\right]}{exp\left[i\kappa {\widehat{\mathcal{V}}}_{\left(3\right)}^{{\rho }_1{\sigma }_1{\rho }_2{\sigma }_2{\rho }_3{\sigma }_3}\frac{\delta }{\delta J^{{\rho }_1{\sigma }_1}}\frac{\delta }{\delta J^{{\rho }_2{\sigma }_2}}\frac{\delta }{\delta J^{{\rho }_3{\sigma }_3}}+\mathcal{O}\left({\kappa }^2\right)\right]exp\left[\frac{i}{2}{J}^{\mu \nu }{\mathcal{O}}_{\mu \nu \alpha \beta }^{-1}\square J^{\alpha \beta }\right]}{|}_{J=0}.\hfill \end{array}\end{array}$$

(12)

The following comments are due. First and foremost, gravity is a gauge theory, so operator ${\mathcal{O}}^{-1}$ does not exist. In analogy with other gauge theories, one uses the Faddeev–Popov ghosts [39,40] or the BRST technique [41,42,43,44] to obtain the propagator. We discuss this issue in detail in another section, since it is irrelevant to the principal discussion. Secondly, the discussed scheme is the standard scheme used within the conventional quantum field theory. In that sense, perturbative quantum gravity is the most general attempt to study quantum gravitational effects while remaining within the standard quantum field theory without any modifications. Lastly, the discussed calculations pose no principal difficulty but a severe computational challenge. The issue is not specific to quantum gravity itself. Even within renormalisable quantum field theories, for instance, within quantum chromodynamics, with each order of perturbation theory, calculations become increasingly complicated and involve more and more contributions [45]. In that sense, perturbative quantum gravity is not an exception to this rule, but rather presents a limiting case with the faster-growing complexity.

Let us turn to a discussion of the physical content of the theory. Three physical premises behind the perturbative approach within the quantum field theory are crucial for the perturbative quantum gravity. Firstly, the perturbation theory only requires propagators, interaction operators, and polarisation operators to generate perturbative expansions. Secondly, propagating and external states are not necessarily free states of a theory. Thirdly, the Poincare group is necessary for the perturbative theory within the conventional quantum field theory. Below, we elaborate on these points.

Firstly, the core of a perturbative theory is propagating degrees of freedom, their interactions, and polarisation operators. Indeed, when given these objects, one can use the standard Feynman graphs technique to generate any perturbative expansion up to any desired order. From the physical point of view, these objects provide all the information about the content of the theory. The polarisation operators describe the states of quantum fields in the initial and final states. Propagators describe how perturbations propagate in the spacetime, while the interaction operators describe their interactions. The perturbative expansion approximates a matrix element as a series of consecutive interactions. These series allow one to define the S matrix, the evolution operator of the theory, describing it in the most general way possible. Therefore, one requires no other objects to create and implement the perturbation theory.

Secondly, one shall distinguish propagating, external and free states. Propagating states are those states described by propagators of the theory. External states are those states described by polarisation operators. Free states are those states which can be separated and captured or measured in a (conceivable) physical experiment. The following considerations show the need to make such a distinguishment.

It is well-established that propagating states differ from those present in initial and final states. The Faddeev–Popov ghosts provide the best example [40]. These ghosts are non-physical, since they enter a theory due to the mathematical redefinition of the path integral integration volume. Consequently, these ghosts do not enter initial or final states, and cannot be separated or detected in any physical experiment. Nevertheless, they must be presented in a perturbation theory to ensure the gauge invariance and unitarity of the theory. This example alone shows the need to distinguish between states that propagate and states that can exist in the initial and final states.

The neutrino mixing mechanism [46,47] explicitly shows the difference between propagating states and states that can be registered. The mechanism explicitly states that there are three neutrino states ${\nu }_1$, ${\nu }_2$, and ${\nu }_3$, which are eigenstates of the mass operator. At the same time, exist three neutrino states ${\nu }_e$, ${\nu }_{\mu }$, and ${\nu }_{\tau }$, which are eigenstates of the interaction operators. These states are not the same but are linearly related via the Pontecorvo–Maki–Nakagawa–Sakata matrix. Neutrino propagators correspond to the eigenstates of the mass operator. However, these states cannot be registered by any physical apparatus, since the eigenstates of the interaction operator are coupled to the matter states. Therefore, this model provides an explicit example of a difference between propagating states and states that can be registered.

Quantum chromodynamics provides the most vivid justification for the need for such a distinction [48,49,50]. Within the theory, quarks experience mixing [51,52], so the same logic is applicable. Quark propagators correspond to the mass operator eigenstates, which differ from the interaction operator eigenstates. The Cabibbo–Kobayashi–Maskawa matrix linearly connects these eigenstates. Further, the confinement makes it impossible to observe quark states directly. The spectrum of observed states is colourless, and does not match the spectrum of propagating states. One uses the distribution functions technique to relate calculations within quantum chromodynamics with observational quantities [53,54,55]. A distribution function describes an external probe’s probability of interacting with a quark or a gluon composing a non-perturbative state, such as atomic nuclei. In other words, one operates with matrix elements having interaction operator eigenstates in the initial or finite states. In this way, quantum chromodynamics not only points to the necessity to distinguish between these states, but actively requires such a distinguishment for all computational purposes.

Lastly, the Poincare group is essential for quantum field theory. The group admits two Casimir operators. A Casimir operator is an operator that commutates with all other operators of the group algebra. Consequently, the eigenstates of a Casimir operator remain invariant under the Poincare group transformations. Two Casimir operators of the Poincare group define the notions of mass and spin, which is the basis for constructing and classifying all quantum fields [56,57,58]. Therefore, without the Poincare group, one shall construct a new field theory significantly different from the standard one. This reasoning does not imply that such theories do not exist, shall not be constructed, or are irrelevant to quantum gravity. On the contrary, such theories are widely developed, with the conformal field theory being (perhaps) the most well-known [59,60,61,62,63]. However, the discussion of their possible implementations for quantum gravity lies far beyond the scope of this paper.

The crucial role of the Poincare group found an implementation within the rapidly developing scattering amplitude techniques [64,65,66,67]. In contrast with the conventional quantum field theory, the scattering amplitude technique does not appeal to the notion of a quantum field. One classifies the initial and final states of the theory as eigenstates of the Poincare group Casimir operators without any references to quantum fields. Transformation features of a matrix element with respect to the Lorentz transformation provide a way to fix the matrix element’s general structure uniquely. The optical theorem shows that a given matrix element has multiple contributions. The consecutive application of the theorem recovers the perturbative expansion as the unavoidable consequence of the unitarity. In this way, the scattering amplitude description allows one to recover the perturbative approach to quantum field theory without directly referencing quantum fields but using only the Poincare group.

These premises point to the following features in understanding perturbative quantum gravity. Firstly, choosing a flat background is necessary to construct a perturbative quantum theory of gravity within the standard quantum field theory framework. Secondly, the theory applies to a scattering of separated gravitational perturbations and matter states, but its applicability may extend beyond this setup.

The existence of the flat background is crucial due to the role of the Poincare group. Let us highlight one more time that it is possible to work within a theory with a different background, both from fundamental and technical points of view. However, such theories will bring us outside the standard quantum field theory formalism. The flat background allows one to use the Poincare group and to define eigenstates of mass and spin operators. In turn, one can define the scattering problem in the standard way described in many textbooks [36,37]. One associates initial and final states with past and future spatial infinity while assuming that interaction occurs in the flat spatial region between these states.

The formalism of perturbative quantum gravity extends beyond the standard scattering problem. Firstly, propagating graviton states may not be associated with the states registered by physical apparatus. Gravitons are propagating states because they are associated with eigenstates of mass and chirality operators. As was pointed out above, this does not imply that these exact states are registered by physical apparatus. The graviton propagator receives highly non-trivial quantum corrections altering its structure [68,69], further indicating that propagating graviton states shall not be associated with states registered by physical apparatus. Finally, there is an ongoing discussion concerning the mere possibility of detecting graviton states directly [70]. Although the main point of this article is related to the subject of our discussion, it will bring us far away from our main aim, so we will not discuss this claim further. Still, even without a discussion about graviton’s detectability, there are enough reasons to believe that graviton states shall be viewed only as propagating states.

These arguments allow one to use perturbative quantum gravity beyond the narrow scope of the scattering problem. There are many branches of research where such an extension of perturbative quantum gravity is assumed implicitly. Perhaps the most recognised one is the calculation of the effective action. Although the detailed discussion of effective action in quantum gravity lies far beyond the scope of this paper (and was done in many publications, for instance [71]), we will touch upon its two most important features.

The effective action technique provides a tool to calculate the classical value (the expectation value) of a quantum field. In other words, one expects to recover the classical dynamics of a system driven by its quantum behaviour. From the technical point of view, to calculate the (n-loop) effective action, one sums (n-loop) one-particle irreducible diagrams. It is well-recognised that even one-loop effective action for a gauge theory (including gravity) is gauge-dependent. In response to this finding, the unique effective action technique was developed [72,73]. It treats the gauge-fixing parameter of a theory as a generalised coordinate and constructs a description of internal field space with curved geometry. In turn, it is possible to define a generalisation of the variational derivative similar to the standard covariant derivative. The unique effective action is defined in terms of such derivatives, which makes it gauge-independent. At the same time, such generalisation of the variational derivative introduces new terms to the effective action. These terms cannot be obtained from one-particle irreducible diagrams, which seems to contradict the original setup. One shall only consider the on-shell effective action to resolve this apparent contradiction. The situation is similar to the calculation of scattering cross sections. Although it is possible to calculate an off-shell scattering cross section, physical apparatus can produce and register only on-shell states. Similarly, the classical value of a quantum field can only be associated with an on-shell object. In this way, the unique effective action removes the gauge dependence from calculations while the fixation of the effective action on the mass shell recovers the classical dynamic of quantum fields.

Lastly, we shall comment on the renormalisation of perturbative quantum gravity. The term “renormalisation” is commonly used in the literature, although it may not be the most suitable for this discussion. For the sake of consistency, we will continue to use this term. However, it is worth noting that the discussion below goes far beyond renormalisation as it is understood in conventional renormalisable models such as quantum electrodynamics, $SU\left(N\right)$ Yang–Mills theory, and the standard model.

The growing consensus is that perturbative quantum gravity does not experience any problems with renormalisation. This point of view does not imply that the theory is well-renormalisable. On the contrary, it is believed that extending the perturbative quantum gravity to a conventionally renormalisable theory is impossible. Consequently, all “problems” of the theory shall be viewed as features, while the complete theory of quantum gravity development shall be performed separately.

The main feature of the perturbative quantum gravity defining its renormalisation behaviour is that it generates a new set of higher dimensional operators at each loop order. This feature was first observed in the classic paper [74], where the one-loop divergences of perturbative quantum gravity were studied. At the one-loop level, the theory generates divergences proportional to operators $R^2$, $R_{\mu \nu }^2$, and $R_{\mu \nu \alpha \beta }^2$. The divergence of the Riemann tensor squared term is irrelevant, since the term can always be brought to the Gauss–Bonnet term, which is the complete derivative in $D=4$.

These divergences vanish for pure general relativity without matter for on-shell amplitudes. On shall states satisfy the vacuum Einstein equations $R_{\mu \nu }=0$, so all the divergent contributions are cancelled. This feature is accidental, since it requires the theory to exist in $D=4$, has no matter degrees of freedom, and works only for on-shell amplitudes. Many other results later confirmed this general theory feature (see, for instance, [12,29,75,76]).

The fact that the theory generates new operators at each level of perturbation theory can be negated from the technical point of view, but the theory loses its predictability. Indeed, let us assume that we calculated a gravitational scattering cross section at the tree level and compared it with empirical data. At the tree level, the theory has a single coupling, the gravitational coupling, and we shall recover its value from the experiment. When we go at the one-loop level, the value of the gravitational coupling is not affected by quantum effects [68,69]. However, the theory develops divergences proportional to $R^2$ and $R_{\mu \nu }^2$ operators. From the technical point of view, these divergences can be subtracted and replaced by two new couplings. In turn, the values of these new couplings shall once again be recovered from the experiment. Consequently, by increasing the level of perturbation theory (the number of loop corrections), we do not increase the precision of calculations. On the contrary, we require more data to match experimental results with each iteration.

These arguments show that perturbative quantum gravity cannot be treated like the conventionally renormalisable theory. From the physical point of view, the root of the problem lies in the fact that the theory is effective. Since it admits a separated energy scale, the Planck scale, it marks the natural limit of its applicability. In full agreement with the standard effective field theory logic [1,2,77], the theory shall include the term terms suppressed by the gravitational coupling the closer to the Planck energy it approaches.

We shall summarise this section to conclude the discussion of the perturbative quantum gravity. Perturbative quantum gravity is a perturbatively constructed quantum field theory with the following distinctive features:

• Gravitons are propagating massless degrees of freedom with chirality $\pm 2$.
• There are no reasons to believe that gravitons shall be associated with degrees of freedom measured by physical apparatus.
• The theory is equivalent to a quantum theory of small metric perturbations.
• The theory has infinite interaction operators, but the same gravitational coupling parametrises them.
• The theory is effective. It develops new operators at each new order of perturbation theory.
• The theory does not experience problems with renormalisation, since it shall not be renormalised due to its effective nature.

Although the discussion can be extended further, especially considering the ongoing study of its quantum features, this section describes the core of the theory relevant to the main aim of the article. With the given background, we can meaningfully address the problem of generating Feynman rules within perturbative quantum gravity.

3. Computational Tools

This section discusses technical tools that make an efficient calculation of the interaction rules for perturbative quantum gravity possible. The main feature of the theory is the factorisation that allows one to separate a quantum gravity action into factors of two types. Factors of the first type contain derivatives and a finite number of terms. The others are free from derivatives, but are infinite series. We discuss this factorisation in the following subsection. Factors of the first type are easy to calculate and require no sophisticated techniques, so we will not discuss them in detail. Factors of the second type present a more significant challenge. Nonetheless, it is possible to establish certain recursive relations, which highly improve the efficiency of their calculations. The second subsection discusses such relations.

We will present results in a series of definitions and theorems. Detailed proofs for many theorems are omitted, as they are not this paper’s primary focus and can be found in other publications.

3.1. Factorisation

We begin with the formal definition of the perturbative metric.

Definition 1.

The perturbative metric $g_{\mu \nu }$ is defined as follows:

$$\begin{array}{c}\hfill g_{\mu \nu }\stackrel{def}{=}{\eta }_{\mu \nu }+\kappa h_{\mu \nu }.\end{array}$$

(13)

Here, $h_{\mu \nu }$ is the small metric perturbations, ${\eta }_{\mu \nu }=diag(+---)$ is the Minkowski metric used to raise and lower indices, and κ is the gravity coupling related with the Newton’s constant $G_{\mathrm{N}}$

$$\begin{array}{c}\hfill {\kappa }^2\stackrel{def}{=}32\pi G_{\mathrm{N}}.\end{array}$$

(14)

This formula is not a truncation of an infinite series, but a finite expression with no omitted terms. The metric still introduces infinite series in the theory.

Theorem 1.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill g^{\mu \nu }& =\sum _{n=0}^{\infty }{(-\kappa )}^n{\left(h^n\right)}^{\mu \nu },\hfill \\ \hfill \sqrt{-g}& =\sum _{n=0}^{\infty }{(-\kappa )}^n\sum _{m=1}^n\frac{1}{m!}{\left(-\frac{1}{2}\right)}^m\left[\sum _{k_1+\cdots +k_m=n}\frac{tr\left(h^{k_1}\right)\cdots tr\left(h^{k_m}\right)}{k_1\cdots k_m}\right],\hfill \\ \hfill {\mathfrak{e}}^m{}_{\mu }& =\sum _{n=0}^{\infty }{\kappa }^n\left(\genfrac{}{}{0pt}{}{\frac{1}{2}}{n}\right){\left(h^n\right)}^m{}_{\mu },\hfill \\ \hfill {\mathfrak{e}}_m{}^{\mu }& =\sum _{n=0}^{\infty }{\kappa }^n\left(\genfrac{}{}{0pt}{}{-\frac{1}{2}}{n}\right){\left(h^n\right)}_m{}^{\mu },\hfill \end{array}\end{array}$$

(15)

Here, $\mathfrak{e}$ is the vierbein, $\left(\genfrac{}{}{0pt}{}{n}{m}\right)$ is the binomial coefficients defined via the Γ function, and the following notations are used:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(h^n\right)}^{\mu \nu }& =h^{\mu }{}_{{\sigma }_1}{h}^{{\sigma }_1}{}_{{\sigma }_2}\cdots h^{{\sigma }_{n-1}\nu },\hfill \\ \hfill \left(h^n\right)& =h^{{\sigma }_1}{}_{{\sigma }_2}{h}^{{\sigma }_2}{}_{{\sigma }_3}\cdots h^{{\sigma }_n}{}_{{\sigma }_1}.\hfill \end{array}\end{array}$$

(16)

Proof. 

The proof of the first expression is trivial. One should use the given expression to calculate $g_{\mu \sigma }{g}^{\sigma \nu }$ and verify its correctness. The papers [17,24,25] provide a more detailed discussion of the proof.

Proof of the second expression is much more complicated. Firstly, one shall use the relation between the determinant and trace of a matrix:

$$\begin{array}{c}\hfill \sqrt{-g}={\left(-det\left[g_{\mu \nu }\right]\right)}^{1/2}={\left(-exp\left[tr\left\{ln\left({\eta }_{\mu \nu }+\kappa h_{\mu \nu }\right)\right\}\right]\right)}^{1/2}.\end{array}$$

(17)

Secondly, one shall factorise one flat spacetime metric and reduce the expression to a single exponent:

$$\begin{array}{c}\hfill \sqrt{-g}={\left(-exp\left[tr\left\{ln\left({\eta }_{\mu \sigma }\left[{\delta }_{\nu }^{\sigma }+\kappa h^{\sigma }{}_{\nu }\right]\right)\right\}\right]\right)}^{1/2}=exp\left[\frac{1}{2}tr\left\{ln\left({\delta }_{\nu }^{\sigma }+\kappa h^{\sigma }{}_{\nu }\right)\right\}\right].\end{array}$$

(18)

Lastly, one expands both exp and ln functions as Taylor series. The resulting expression is an infinite power series, with each term also being a power series. Nonetheless, it is possible to rearrange the terms of this series to obtain the desired expression. The publication [24] provides a detailed explanation of the derivation.

Lastly, the expression for vierbein was first obtained in [17]. One assumes the vierbein admits a power series expansion with unknown coefficients. The values of coefficients are fixed uniquely, so the expansion fits well-known relations for the vierbein. □

Only these objects and their combinations generate all infinite expansions that may enter a quantum gravity action within Riemann geometry. However, these are not all the objects present in the theory. We still have to discuss the Christoffel symbols and the Riemann tensor. The following theorem gives their structure.

Theorem 2.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\alpha \mu \nu }=& \frac{\kappa }{2}\left[{\partial }_{\mu }{h}_{\nu \alpha }+{\partial }_{\nu }{h}_{\mu \alpha }-{\partial }_{\alpha }{h}_{\mu \nu }\right],\hfill \\ \hfill {\mathrm{\Gamma }}_{\mu \nu }^{\alpha }=& \frac{\kappa }{2}\left[{\eta }^{\mu \nu }-\kappa h^{\mu \nu }+{\kappa }^2{h}^{\mu \sigma }{h}_{\sigma }{}^{\nu }+\mathcal{O}\left({\kappa }^3\right)\right]\left[{\partial }_{\mu }{h}_{\nu \beta }+{\partial }_{\nu }{h}_{\mu \beta }-{\partial }_{\beta }{h}_{\mu \nu }\right],\hfill \\ \hfill {\left({\mathrm{\Gamma }}_{\mu }\right)}_{ab}=& {\mathfrak{e}}_a{}^{\alpha }{\mathfrak{e}}_b{}^{\beta }\kappa \left[{\partial }_{\beta }{h}_{\alpha \mu }-{\partial }_{\alpha }{h}_{\beta \mu }\right].\hfill \end{array}\end{array}$$

(19)

Proof. 

Proof of the first two expressions is trivial, and relies only on the definition of the Christoffel symbols

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mu \nu }^{\alpha }\stackrel{def}{=}\frac{1}{2}{g}^{\alpha \beta }\left[{\partial }_{\mu }{g}_{\nu \beta }+{\partial }_{\nu }{g}_{\mu \beta }-{\partial }_{\beta }{g}_{\mu \nu }\right].\end{array}$$

(20)

The last expression is proved as follows. Firstly, one shall use the definition of the spin connection:

$$\begin{array}{c}\hfill {\left({\mathrm{\Gamma }}_{\mu }\right)}_{ab}\stackrel{def}{=}{\mathfrak{e}}_a{}^{\alpha }{\mathfrak{e}}_b{}^{\beta }{\mathrm{\Gamma }}_{\alpha \mu \beta }+g_{\rho \sigma }{\mathfrak{e}}_a{}^{\rho }{\partial }_{\mu }{\mathfrak{e}}_b{}^{\sigma }.\end{array}$$

(21)

The spin connection is antisymmetric for indices a and b. For the perturbative metric, the last term is constituted only by symmetric matrices, so it is symmetric with respect to a and b. Therefore, the last term does not contribute to the spin connection. For the first term, one shall use the explicit expression for the Christoffel symbol and make it antisymmetric. Similarly, paper [24] discusses the derivation in more detail. □

The obtained expressions allow us to describe the structure of the Riemann tensor with the following theorem.

Theorem 3.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill R_{\mu \nu }{}^{\alpha }{}_{\beta }& =g^{\alpha \lambda }\left[{\partial }_{\mu }{\mathrm{\Gamma }}_{\lambda \nu \beta }-{\partial }_{\nu }{\mathrm{\Gamma }}_{\lambda \mu \beta }+g^{\rho \sigma }\left\{{\mathrm{\Gamma }}_{\rho \nu \lambda }{\mathrm{\Gamma }}_{\sigma \mu \beta }-{\mathrm{\Gamma }}_{\rho \mu \lambda }{\mathrm{\Gamma }}_{\sigma \nu \beta }\right\}\right],\hfill \\ \hfill R_{\mu \nu \alpha \beta }& ={\partial }_{\mu }{\mathrm{\Gamma }}_{\alpha \nu \beta }-{\partial }_{\nu }{\mathrm{\Gamma }}_{\alpha \mu \beta }+g^{\rho \sigma }\left\{{\mathrm{\Gamma }}_{\rho \nu \alpha }{\mathrm{\Gamma }}_{\sigma \mu \beta }-{\mathrm{\Gamma }}_{\rho \mu \alpha }{\mathrm{\Gamma }}_{\sigma \nu \beta }\right\},\hfill \\ \hfill R_{\mu \nu }& =g^{\rho \sigma }\left[{\partial }_{\rho }{\mathrm{\Gamma }}_{\sigma \mu \nu }-{\partial }_{\mu }{\mathrm{\Gamma }}_{\sigma \rho \nu }\right]+g^{\rho \sigma }{g}^{\lambda \tau }\left[{\mathrm{\Gamma }}_{\tau \mu \nu }{\mathrm{\Gamma }}_{\lambda \rho \sigma }-{\mathrm{\Gamma }}_{\lambda \mu \sigma }{\mathrm{\Gamma }}_{\tau \nu \rho }\right],\hfill \\ \hfill R& =g^{\mu \nu }{g}^{\alpha \beta }{\partial }_{\mu }\left[{\mathrm{\Gamma }}_{\nu \alpha \beta }-{\mathrm{\Gamma }}_{\alpha \nu \beta }\right]+g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }\left[{\mathrm{\Gamma }}_{\alpha \mu \nu }{\mathrm{\Gamma }}_{\beta \rho \sigma }-{\mathrm{\Gamma }}_{\alpha \mu \rho }{\mathrm{\Gamma }}_{\beta \nu \sigma }\right].\hfill \end{array}\end{array}$$

(22)

Proof. 

Expressions for the Ricci tensor and the scalar curvature can be derived directly from the expressions for the Riemann tensor, making their derivation purely technical. The derivation of the first two formulas is technically straightforward, but requires lengthy computations. To obtain the desired result, one shall use a single relation for the inverse metric derivative:

$$\begin{array}{c}\hfill {\partial }_{\mu }{g}^{\alpha \beta }=-g^{\alpha \rho }{g}^{\beta \sigma }{\partial }_{\mu }{g}_{\rho \sigma }.\end{array}$$

(23)

The remaining part of the proof relies on the manipulation of indices. □

These theorems constitute the factorisation theorem crucial for the perturbative quantum gravity.

Theorem 4.

In Riemann geometry, a gravity action evaluated with the perturbative metric can always take a form where all terms that begin infinite series are free from derivatives and expressed via $\sqrt{-g}$, $g^{\mu \nu }$, and ${\mathfrak{e}}_m{}^{\mu }$.

In other words, as long as we operate within Riemann geometry, we can always split a gravity action into two different parts. The part that involves derivatives is finite, and can be calculated explicitly. The other part is an infinite series, but it is always expressed in terms of elementary series for the volume factors $\sqrt{-g}$, the inverse metric $g^{\mu \nu }$, and the vierbein ${\mathfrak{e}}_m{}^{\mu }$.

3.2. Recursive Relations

The previous section provides a way to factorise a gravity action, but we must still address calculational challenges. Namely, we shall elaborate on the method to describe the structure of perturbative expansions for the inverse metric, volume factor, and vierbein. We address this challenge in this section.

Definition 2.

We define the plain I-tensor of the n-th order as follows:

$$\begin{array}{c}\hfill I_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\stackrel{def}{=}{\eta }^{{\sigma }_1{\rho }_2}{\eta }^{{\sigma }_2{\rho }_3}\cdots {\eta }^{{\sigma }_n{\rho }_1}.\end{array}$$

(24)

The tensor has no additional symmetries. We introduce an additional tensor to account for the essential symmetries of the interaction rules.

Definition 3.

We define the $\mathcal{I}$-tensor of the n-th order as follows:

$$\begin{array}{c}\hfill {\mathcal{I}}_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\stackrel{def}{=}\frac{1}{2^n}\frac{1}{n!}\left[I_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}+\mathit{permutations}\right].\end{array}$$

(25)

Permutations account for all terms that make the $\mathcal{I}$-tensor symmetric with respect to permutations of indices within each index pair ${\mu }_i\leftrightarrow {\nu }_i$, and with respect to permutations of any two index pairs $\{{\mu }_i,{\nu }_i\}\leftrightarrow \{{\mu }_j,{\nu }_j\}$.

The $\mathcal{I}$ tensor of the n-th order has $n!2^n$ terms. In other words, the number of terms grows faster than the factorial. This feature holds for tensors discussed below. Although this feature presents a computational challenge, it is an essential part of the standard perturbative approach. The final expression for a gravity interaction vertex shall be symmetric with respect to graviton permutations, so it is necessary to operate with tensor structures respecting this symmetry.

The plain I-tensor provides a way to operate with powers of $h_{\mu \nu }$ and its traces given by the following theorem. The theorem is trivial, so we omit the proof.

Theorem 5.

The plain I-tensor defines the structure of powers of $h_{\mu \nu }$ and its traces.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\left(h^n\right)}^{\mu \nu }=I_{(1+n)}^{\mu \nu {\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n},\hfill \\ \hfill & tr\left(h^n\right)=I_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\hfill \end{array}\end{array}$$

(26)

Using this theorem, one can demonstrate how the plain I-tensor describes the perturbative structure of the inverse metric. The following theorem, which describes this relation, is a direct corollary of previous theorems, so we omit its proof.

Theorem 6.

The inverse metric perturbative expansion in terms of the plain I tensor reads:

$$\begin{array}{cc}\hfill g^{\mu \nu }& =\sum _{n=0}^{\infty }{(-1)}^n{\kappa }^n{I}_{(1+n)}^{\mu \nu {\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\hfill \end{array}$$

(27)

This theorem provides an easy way to study the features of the I-tensor, such as the following easily proven contraction feature.

Theorem 7.

The I-tensor admits the following contraction feature:

$$\begin{array}{c}\hfill {\eta }_{\mu \nu }{I}^{\mu \nu {\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}=I^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}.\end{array}$$

(28)

Proof. 

We shall start with the following general relation:

$$\begin{array}{c}\hfill d=g^{\mu \nu }{g}_{\mu \nu }.\end{array}$$

(29)

This relation can be expanded:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill d=& \left({\eta }_{\mu \nu }+\kappa h_{\mu \nu }\right)\sum _{n=0}^{\infty }{(-1)}^n{\kappa }^n{I}_{(1+n)}^{\mu \nu {\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}\hfill \\ \hfill =& {\eta }_{\mu \nu }{\eta }^{\mu \nu }+\sum _{n=1}^{\infty }{(-1)}^n{\kappa }^n{\eta }_{\mu \nu }{I}_{(1+n)}^{\mu \nu {\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}\hfill \\ \hfill & +\sum _{n=0}^{\infty }{(-1)}^n{\kappa }^{n+1}{I}_{(1+n)}^{\mu \nu {\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{\mu \nu }{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\hfill \end{array}\end{array}$$

(30)

The desired relation holds as each term in the infinite sums is cancelled. □

We uncover the perturbative structure of the volume factor in a similar way. Initially, we give a non-constructive definition of the plain C-tensor. Further, we demonstrate a specific recursive relation between C-tensors of different orders. This relation allows one to calculate a given order’s plain C-tensor.

Definition 4.

We define the plain C-tensor of n-th order in such a way to describe the perturbative structure of the volume factor:

$$\begin{array}{c}\hfill \sqrt{-g}\stackrel{def}{=}\sum _{n=0}^{\infty }{\kappa }^n{C}_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\end{array}$$

(31)

Similarly to the previous case, the definition does not require specific symmetries among the indices. The next definition of the $\mathcal{C}$-tensor introduces a tensor with suitable symmetries.

Definition 5.

We defined the $\mathcal{C}$-tensor of the n-th order as follows:

$$\begin{array}{c}\hfill {\mathcal{C}}_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\stackrel{def}{=}\frac{1}{2^n}\frac{1}{n!}\left[C_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}+\mathit{permutations}\right].\end{array}$$

(32)

Permutations ensure $\mathcal{C}$-tensor is symmetric with respect to index pair and pair permutations.

The definition is not constructive, since it does not show a way to compute the tensor explicitly; instead, it only describes its function. However, obtaining a recursive relation that facilitates an uncomplicated method to construct C-tensors is possible.

Theorem 8.

The following recursive relation holds.

$$\begin{array}{c}\hfill C_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}=\frac{1}{2n}\sum _{k=1}^n{(-1)}^{k-1}{I}_{\left(k\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_k{\sigma }_k}{C}_{(n-k)}^{{\rho }_{k+1}{\sigma }_{k+1}\cdots {\rho }_n{\sigma }_n}.\end{array}$$

(33)

Proof. 

Firstly, we introduce an auxiliary metric:

$$\begin{array}{c}\hfill {\mathbf{g}}_{\mu \nu }={\eta }_{\mu \nu }+z\kappa h_{\mu \nu }.\end{array}$$

(34)

Here, z is an arbitrary positive real number. This metric has no physical or mathematical meaning except to obtain the discussed relation. It matches the original metric when $z=1$.

Secondly, we compute the following derivative:

$$\begin{array}{c}\hfill {\left.\frac{d}{dz}\sqrt{-\mathbf{g}}\right|}_{z=1}=\frac{\kappa }{2}\sqrt{-g}{g}^{\mu \nu }{h}_{\mu \nu }.\end{array}$$

(35)

We applied Jacobi’s formula for a non-degenerate matrix X:

$$\begin{array}{c}\hfill ddetX=detXtr\left[X^{-1}dX\right].\end{array}$$

(36)

Let us calculate the same derivative using the perturbative expansion:

$$\begin{array}{c}\hfill {\left.\frac{d}{dz}\sqrt{-\mathbf{g}}\right|}_{z=1}=\sum _{n=1}^{\infty }n{\kappa }^n{C}_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\end{array}$$

(37)

Comparing these expressions, we obtain a relationship that does involve $\mathbf{g}$:

$$\begin{array}{c}\hfill \frac{\kappa }{2}\sqrt{-g}{g}^{\mu \nu }{h}_{\mu \nu }=\sum _{n=1}^{\infty }n{\kappa }^n{C}_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\end{array}$$

(38)

The left-hand side of this expression can be calculated perturbatively:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\kappa }{2}\sqrt{-g}{g}^{\mu \nu }{h}_{\mu \nu }& =\frac{\kappa }{2}\sum _{p_1=0}^{\infty }{\kappa }^{p_1}{C}_{\left(p_1\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_{p_1}{\sigma }_{p_1}}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_{p_1}{\sigma }_{p_1}}\sum _{p_2=0}^{\infty }{(-1)}^{p_2}{\kappa }^{p_2}{I}_{(p_2+1)}^{\mu \nu {\lambda }_1{\tau }_1\cdots {\lambda }_{p_2}{\tau }_{p_2}}{h}_{{\lambda }_1{\tau }_1}\cdots h_{{\lambda }_{p_2}{\tau }_{p_2}}{h}_{\mu \nu }\hfill \\ \hfill & =\sum _{n=1}^{\infty }\sum _{k=1}^n\frac{{(-1)}^{k-1}}{2}{\kappa }^n{I}_{\left(k\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_k{\sigma }_k}{C}_{(n-k)}^{{\rho }_{k+1}{\sigma }_{k+1}\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}\hfill \\ \hfill & =\sum _{n=1}^{\infty }{\kappa }^n\sum _{k=1}^n\frac{{(-1)}^{k-1}}{2}{I}_{\left(k\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_k{\sigma }_k}{C}_{(n-k)}^{{\rho }_{k+1}{\sigma }_{k+1}\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\hfill \end{array}\end{array}$$

(39)

After comparing these expressions, we can derive the required recursive relation.

$$\begin{array}{c}\hfill C_{\left(n\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}=\frac{1}{2n}\sum _{k=1}^n{(-1)}^{k-1}{I}_{\left(k\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_k{\sigma }_k}{C}_{(n-k)}^{{\rho }_{k+1}{\sigma }_{k+1}\cdots {\rho }_n{\sigma }_n}.\end{array}$$

(40)

□

Equation (33) defines a plain C-tensor of n-th order recursively using plain C-tensors of lower orders. One would desire to express $C_{\left(n\right)}$ using $C_{(n-1)}$ alone, but unfortunately, it appears to be impossible. The reason is that formula (33) involves all low orders that must contribute to the plain C-tensor. Since the volume factor $\sqrt{-g}$ is the simple object, we conjecture that no simpler recursive formula exists.

The same logic shall be used to construct more complicated structures generated by the volume factor and the inverse metric. Firstly, we give a non-constructive definition of a tensor. Secondly, we find a recursive relation that defines the required tensor iteratively.

Definition 6.

We define the family of plain $C_{\left(n\right)}$-tensors as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \sqrt{-g}{g}^{\mu \nu }\stackrel{def}{=}& \sum _{n=0}^{\infty }{\kappa }^n{C}_{(1;n)}^{\mu \nu ,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n},\hfill \\ \hfill \sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }\stackrel{def}{=}& \sum _{n=0}^{\infty }{\kappa }^n{C}_{(2;n)}^{\mu \nu \alpha \beta ,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n},\hfill \\ \hfill \sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }\stackrel{def}{=}& \sum _{n=0}^{\infty }{\kappa }^n{C}_{(3;n)}^{\mu \nu \alpha \beta \rho \sigma ,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n},\hfill \\ \hfill \sqrt{-g}{g}^{{\mu }_1{\nu }_1}\cdots g^{{\mu }_l{\nu }_l}\stackrel{def}{=}& \sum _{n=0}^{\infty }{\kappa }^n{C}_{(l;n)}^{{\mu }_1{\nu }_1\cdots {\mu }_l{\nu }_l,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\hfill \end{array}\end{array}$$

(41)

These tensors do not have any additional symmetry. The corresponding symmetric tensors are defined as follows.

Definition 7.

$$\begin{array}{c}\hfill {\mathcal{C}}_{(l;n)}^{{\mu }_1{\nu }_1\cdots {\mu }_l{\nu }_l,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\stackrel{def}{=}\frac{1}{2^n}\frac{1}{n!}\left[C_{(l;n)}^{{\mu }_1{\nu }_1\cdots {\mu }_l{\nu }_l,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}+\mathit{permutations}\right].\end{array}$$

(42)

Permutations ensure symmetric $\mathcal{C}$-tensor with respect to index pair and pair permutation. It is important to note that symmetrisation accounts only for indices contracted with the perturbation indices.

Theorem 9.

The following recursive relation holds for plain $C_{\left(l\right)}$-tensors.

$$\begin{array}{c}\hfill C_{(l;n)}^{{\mu }_1{\nu }_1{\mu }_2{\nu }_2\cdots {\mu }_l{\nu }_l,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}=\sum _{p=0}^n{(-1)}^p{I}_{(1+p)}^{{\mu }_1{\nu }_1{\rho }_1{\sigma }_1\cdots {\rho }_p{\sigma }_p}{C}_{(l-1;n-p)}^{{\mu }_2{\nu }_2\cdots {\mu }_l{\nu }_l,{\rho }_{p+1}{\sigma }_{p+1}\cdots {\rho }_n{\sigma }_n}.\end{array}$$

(43)

Proof. 

The proof relies on the possibility of factorising a single metric from the perturbative expression. We can use this fact to define the $C_{\left(l\right)}$ tensors.

$$\begin{array}{cc}\hfill \sqrt{-g}{g}^{{\mu }_1{\nu }_1}\cdots g^{{\mu }_l{\nu }_l}\stackrel{def}{=}& \sum _{n=0}^{\infty }{\kappa }^n{C}_{(l;n)}^{{\mu }_1{\nu }_1\cdots {\mu }_l{\nu }_l,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\hfill \end{array}$$

(44)

At the same time, we can perform the perturbative expansion of $g^{{\mu }_1{\nu }_1}$ separately.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \sqrt{-g}{g}^{{\mu }_1{\nu }_1}{g}^{{\mu }_2{\nu }_2}\cdots g^{{\mu }_l{\nu }_l}=g^{{\mu }_1{\nu }_1}\sqrt{-g}{g}^{{\mu }_2{\nu }_2}\cdots g^{{\mu }_l{\nu }_l}\hfill \\ \hfill & =\sum _{n_1=0}^{\infty }{(-1)}^{n_1}{\kappa }^{n_1}{I}_{(n_1+1)}^{{\mu }_1{\nu }_1{\rho }_1{\sigma }_1\cdots {\rho }_{n_1}{\sigma }_{n_1}}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_{n_1}{\sigma }_{n_1}}\sum _{n_2=0}^{\infty }{\kappa }^{n_2}{C}_{(l-1;n_2)}^{{\mu }_2{\nu }_2\cdots {\mu }_l{\nu }_l,{\lambda }_1{\tau }_1\cdots {\lambda }_{n_2}{\tau }_{n_2}}{h}_{{\lambda }_1{\tau }_1}\cdots h_{{\lambda }_{n_2}{\tau }_{n_2}}\hfill \\ \hfill & =\sum _{n_1=0}^{\infty }\sum _{n_2=0}^{\infty }{(-1)}^{n_1}{\kappa }^{n_1+n_2}{I}_{(n_1+1)}^{{\mu }_1{\nu }_1{\rho }_1{\sigma }_1\cdots {\rho }_{n_1}{\sigma }_{n_1}}{C}_{(l-1;n_2)}^{{\mu }_2{\nu }_2\cdots {\mu }_l{\nu }_l,{\rho }_{n_1+1}{\sigma }_{n_1+1}\cdots {\rho }_{n_1+n_2}{\sigma }_{n_1+n_2}}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_{n_1+n_2}{\sigma }_{n_1+n_2}}.\hfill \end{array}\end{array}$$

(45)

One can decouple the summation indices by redefining them, resulting in the following formula:

$$\begin{array}{c}\hfill \sqrt{-g}{g}^{{\mu }_1{\nu }_1}{g}^{{\mu }_2{\nu }_2}\cdots g^{{\mu }_l{\nu }_l}=\sum _{n=0}^{\infty }\sum _{p=0}^n{\kappa }^n{(-1)}^p{I}_{(1+p)}^{{\mu }_1{\nu }_1{\rho }_1{\sigma }_1\cdots {\rho }_p{\sigma }_p}{C}_{(l-1;n_2)}^{{\mu }_2{\nu }_2\cdots {\mu }_l{\nu }_l,{\rho }_{p+1}{\sigma }_{p+1}\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\end{array}$$

(46)

The proof is concluded by making a direct comparison with the definition of the original tensor $C_{\left(l\right)}$. □

Finally, we treat vierbein similarly. We define the plain E-tensor and plain $C_E$-tensor to encapsulate the structure of the following factors.

Theorem 10.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathfrak{e}}_{\mu }{}^{\nu }\stackrel{def}{=}& \sum _{n=0}^{\infty }{\kappa }^n{E}_{\mu \left(n\right)}^{\nu {\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}\hfill \\ \hfill =& \sum _{n=0}^{\infty }{\kappa }^n\left(\genfrac{}{}{0pt}{}{-\frac{1}{2}}{n}\right)I_{\mu \left(n\right)}^{\nu {\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n},\hfill \end{array}\end{array}$$

(47)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \sqrt{-g}{\mathfrak{e}}_{\mu }{}^{\nu }\stackrel{def}{=}& \sum _{n=0}^{\infty }{\kappa }^n{C_E}_{\mu \left(n\right)}^{\nu ,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}\hfill \\ \hfill =& \sum _{n=0}^{\infty }{\kappa }^n\sum _{p=0}^n\left(\genfrac{}{}{0pt}{}{-\frac{1}{2}}{n}\right)C_{\left(p\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_p{\sigma }_p}{I}_{\mu (1+n-p)}^{\nu {\rho }_{p+1}{\sigma }_{p+1}\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\hfill \end{array}\end{array}$$

(48)

These tensors do not possess any additional symmetry. In practice, we never use the plain E-tensor, so we do not discuss it further. Moreover, the plain E-tensor match the plain I-tensor up to a constant. On the contrary, the plain $C_E$-tensor is relevant for Dirac fermions, so we define its symmetric generalisation.

Definition 8.

$$\begin{array}{c}\hfill {{\mathcal{C}}_E}_{\mu \left(n\right)}^{\nu ,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\stackrel{def}{=}\frac{1}{2^n}\frac{1}{n!}\left[{C_E}_{\mu \left(n\right)}^{\nu ,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}+\mathit{permutations}\right].\end{array}$$

(49)

Permutations account for all terms that make the ${\mathcal{C}}_E$-tensor symmetric with respect to permutations of indices within each index pair ${\mu }_i\leftrightarrow {\nu }_i$, and with respect to permutations of any two index pairs $\{{\mu }_i,{\nu }_i\}\leftrightarrow \{{\mu }_j,{\nu }_j\}$.

The recursive relation for the $C_E$ tensor holds similarly to the $C_{\left(l\right)}$ tensor. The proof for $C_E$-tensor is constructed the same way as the proof for $C_{\left(l\right)}$-tensor and is omitted.

Theorem 11.

$$\begin{array}{c}\hfill {C_E}_{\mu \left(n\right)}^{\nu ,{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}=\sum _{p=0}^n{E}_{\mu (1+p)}^{\nu {\rho }_1{\sigma }_1\cdots {\rho }_p{\sigma }_p}{C}_{(n-p)}^{{\rho }_{p+1}{\sigma }_{p+1}\cdots {\rho }_n{\sigma }_n}.\end{array}$$

(50)

The main result of this section is Theorems 8, 9 and 11. They are essential from the computational point of view, since their implementation significantly improves the code performance. We believe no simpler formula can be found based on the $\sqrt{-g}$ perturbative expansion structure. All the factors discussed above contain $\sqrt{-g}$, which makes it impossible to construct recursive relations for them without a recursive relation for $\sqrt{-g}$. Theorem 8 shows that the n-th perturbative order incorporates all lower orders, making it impossible to establish a relation between the n-th and $(n-1)$-th orders alone. Therefore, we believe it is impossible to construct truly recursive relations for all other factors.

4. Feynman Rules

The computational tools discussed in the previous section provide a comprehensive framework for deriving Feynman’s rules for perturbative quantum gravity. This section covers the derivation of interaction rules for simple scalars, Horndeski gravity, Dirac fermions, massive and massless vector fields, $SU\left(N\right)$ Yang–Mills model, and general relativity. This collection of models is sufficient to account for quantum effects within the standard model.

We handle all interaction rules as follows. We use the factorisation theorem to separate terms that are infinite series from terms that are finite but contain derivatives. The finite part is calculated explicitly, while the infinite part is expanded in the infinite perturbative series. The discussed recursive relations provide a way to evaluate each term in such a perturbative expansion. Therefore, pointing to the structure contributing to a given perturbative order is sufficient.

We will not discuss the derivation of the Feynman rules from the action via the path integral formalism. The procedure is discussed in multiple sources, including classical textbooks [78,79,80,81]. The procedure also uses Fourier transformation to generate rules in the momentum representation. We also do not discuss the Fourier transformation for the same reason.

Lastly, we shall introduce new notations for the practical sake. In perturbative expansions considered below, using $\mathcal{I}$, $\mathcal{C}$, and other previously defined tensors is not optimal. These tensors play a crucial role in the technical side of the calculation, but they do not provide the most practical description.

If the quantity X admits a perturbative expansion in small metric perturbations, then we shall use the following notations:

$$\begin{array}{c}\hfill X\stackrel{\mathrm{note}}{=}\sum _{n=0}^{\infty }{\kappa }^n{\left(X\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\end{array}$$

(51)

Here, ${\left(X\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}$ is the contribution to the n-th order of the expansion. One can define via derivatives:

$$\begin{array}{c}\hfill {\left(X\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}=\frac{\delta }{\delta h_{{\rho }_1{\sigma }_1}}\cdots \frac{\delta }{\delta h_{{\rho }_n{\sigma }_n}}X{|}_{h=0}.\end{array}$$

(52)

Let us present the following expressions to illustrate these notations. The components of the inverse metric are noted as follows:

$$\begin{array}{c}\hfill g^{\mu \nu }=\sum _{n=0}^{\infty }{\kappa }^n{\left(g^{\mu \nu }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n}.\end{array}$$

(53)

The explicit values of the first few of these components are:

$$\begin{array}{cccc}\hfill \left(g^{\mu \nu }\right)& ={\eta }^{\mu \nu }={\mathcal{I}}^{\mu \nu },\hfill & \hfill {\left(g^{\mu \nu }\right)}^{\alpha \beta }& =-\frac{1}{2}\left({\eta }^{\mu \alpha }{\eta }^{\nu \beta }+{\eta }^{\mu \beta }{\eta }^{\nu \alpha }\right)=-{\mathcal{I}}^{\mu \nu \alpha \beta }.\hfill \end{array}$$

(54)

Implementation of these notations for the volume factor produces the following:

$$\begin{array}{c}\hfill \sqrt{-g}=\sum _{n=0}^{\infty }{\kappa }^n{\left(\sqrt{-g}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{h}_{{\rho }_1{\sigma }_1}\cdots h_{{\rho }_n{\sigma }_n},\end{array}$$

(55)

$$\begin{array}{cccccc}\hfill \left(\sqrt{-g}\right)& =1,\hfill & \hfill {\left(\sqrt{-g}\right)}^{\mu \nu }& =\frac{1}{2}{\eta }^{\mu \nu }={\mathcal{C}}^{\mu \nu },\hfill & \hfill {\left(\sqrt{-g}\right)}^{\mu \nu \alpha \beta }& =\frac{1}{8}\left(-{\eta }^{\mu \alpha }{\eta }^{\nu \beta }-{\eta }^{\mu \beta }{\eta }^{\nu \alpha }+{\eta }^{\mu \nu }{\eta }^{\alpha \beta }\right)={\mathcal{C}}^{\mu \nu \alpha \beta }.\hfill \end{array}$$

(56)

4.1. Scalar Field with a Potential

The action for a single free scalar field with mass m with the minimal coupling to gravity reads:

$$\begin{array}{c}\hfill {\mathcal{A}}_{s=0}=\int d^{4}x\sqrt{-g}\left[\frac{1}{2}{\left(\nabla \varphi \right)}^2-\frac{m^2}{2}{\varphi }^2\right].\end{array}$$

(57)

Following the factorisation theorem, the action takes the following form:

$$\begin{array}{c}\hfill {\mathcal{A}}_{s=0}=\int d^{4}x\left[\frac{1}{2}\sqrt{-g}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -\frac{m^2}{2}\sqrt{-g}{\varphi }^2\right].\end{array}$$

(58)

The following theorem describes its structure in the Fourier representation.

Theorem 12.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{s=0}=& \sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times {\kappa }^n\left[-\frac{1}{2}{\left(\sqrt{-g}{g}^{\mu \nu }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{\mathcal{I}}_{\mu \nu \alpha \beta }{\left(p_1\right)}^{\alpha }{\left(p_2\right)}^{\beta }-\frac{m^2}{2}{\left(\sqrt{-g}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\right]\varphi \left(p_1\right)\varphi \left(p_2\right).\hfill \end{array}\end{array}$$

(59)

Here, $k_i$ are momenta of gravitons, $p_1$ and $p_2$ are momenta of scalars, and the same definition of $\mathcal{I}$ tensor is used:

$$\begin{array}{c}\hfill {\mathcal{I}}^{\mu \nu \alpha \beta }=\frac{1}{2}\left({\eta }^{\mu \alpha }{\eta }^{\nu \beta }+{\eta }^{\mu \beta }{\eta }^{\nu \alpha }\right).\end{array}$$

(60)

The contribution of order $\mathcal{O}\left({\kappa }^0\right)$ describes the standard scalar propagator:

(61)

The other terms define the coupling of scalar field kinetic energy to gravity:

(62)

The dotted line on the left side of the diagram notes the presence of $n\ge 1$ graviton lines.

We can derive the equation for the scalar field potential similarly. The potential $V\left(\varphi \right)$ is a Taylor series in the scalar field $\varphi$. Each term in this expansion represents a scalar field self-interaction coupled with gravity. Thus, deriving the interaction rule for a single power-law potential is sufficient.

We consider a power-law potential with $p\ge 3$ being a whole number and ${\lambda }_p$ being a coupling with the mass dimension $4-p$:

$$\begin{array}{c}\hfill \mathcal{A}=\int d^{4}x\sqrt{-g}\left[\frac{{\lambda }_p}{p!}{\varphi }^p\right]=\int d^{4}x\left[\sqrt{-g}\frac{{\lambda }_p}{p!}{\varphi }^p\right].\end{array}$$

(63)

The following theorem describes the action in Fourier representation.

Theorem 13.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{A}=& \sum _{n=0}^{\infty }\int \prod _{j=1}^p\frac{d^4{q}_j}{{\left(2\pi \right)}^4}\prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(\sum q_j+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times {\kappa }^n\frac{{\lambda }_p}{p!}{\left(\sqrt{-g}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\varphi \left(q_1\right)\cdots \varphi \left(q_p\right).\hfill \end{array}\end{array}$$

(64)

This results in the following interaction rule:

(65)

This expression fully describes the gravitational coupling of scalar field potentials for any whole $p\ge 3$.

4.2. Horndeski’s Gravity

The Horndeski theory is the most general scalar–tensor theory of gravity that admits second-order field equations and minimal coupling to matter. It was discovered in [82], and independently rediscovered in [83]. Because of the second-order field equations, the theory is free from the Ostrogradsky instability [84,85,86].

A combination of the following Lagrangians gives the theory:

$$\begin{array}{c}\hfill \mathcal{A}=\int d^{4}x\sqrt{-g}\left[{\mathcal{L}}_2+{\mathcal{L}}_3+{\mathcal{L}}_4+{\mathcal{L}}_5\right],\end{array}$$

(66)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{L}}_2& =G_2(\varphi ,X),\hfill \\ \hfill {\mathcal{L}}_3& =G_3(\varphi ,X)\square \varphi ,\hfill \\ \hfill {\mathcal{L}}_4& =G_4(\varphi ,X)R+G_{4,X}\left[{\left(\square \varphi \right)}^2-{\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^2\right],\hfill \\ \hfill {\mathcal{L}}_5& =G_5(\varphi ,X)G_{\mu \nu }{\nabla }^{\mu }{\nabla }^{\nu }\varphi -\frac{1}{6}{G}_{5,X}\left[{\left(\square \varphi \right)}^3-3\square \varphi {\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^2+2{\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^3\right].\hfill \end{array}\end{array}$$

(67)

Here, $G_i=G_i(\varphi ,X)$ are functions of the scalar field $\varphi$ and its kinetic term $X=g^{\mu \nu }{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi$, while $G_{4,X}$ and $G_{5,X}$ note derivatives with respect to the kinetic term. Lagrangians ${\mathcal{L}}_2$ and ${\mathcal{L}}_3$ describe scalar field self-interaction with the minimal coupling to gravity. In turn, ${\mathcal{L}}_4$ and ${\mathcal{L}}_5$ describe healthy non-minimal coupling to gravity. Terms containing $G_{4,X}$ and $G_{5,X}$ are relevant only when $G_4$ and $G_5$ depend on the scalar field kinetic term. In that case, they cancel out higher derivative contributions to the field equations.

4.2.1. Horndeski $G_2$ Class

The $G_2$ class of Horndeski theory describes a minimal coupling between the scalar field and gravity. This coupling is minimal because it does not involve graviton momenta.

The function $G_2$ shall be smooth enough to admit a Taylor expansion:

$$\begin{array}{cc}\hfill G_2(\varphi ,X)& =\sum _{a=0}^{\infty }\sum _{b=0}^{\infty }{\lambda }_{(a,b)}{\varphi }^a{X}^b=\sum _{a=0}^{\infty }\sum _{b=0}^{\infty }{\lambda }_{(a,b)}{\varphi }^a{g}^{{\mu }_1{\nu }_1}\cdots g^{{\mu }_b{\nu }_b}{\partial }_{{\mu }_1}\varphi {\partial }_{{\nu }_1}\varphi \cdots {\partial }_{{\mu }_b}\varphi {\partial }_{{\nu }_b}\varphi .\hfill \end{array}$$

(68)

Here, ${\lambda }_{(a,b)}$ is the dimensional coupling with the mass dimension $4(1-b)-a$. Without the loss of generality, one can consider a single term of this expansion and consider the following action:

$$\begin{array}{c}\hfill \int d^{4}x\sqrt{-g}{G}_2(\varphi ,X)\to \int d^{4}x\sqrt{-g}{\lambda }_{(a,b)}{\varphi }^a{X}^b.\end{array}$$

(69)

The perturbative expansion is given by the following theorem.

Theorem 14.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}{\lambda }_{(a,b)}{\varphi }^a{X}^b=\int d^{4}x\sqrt{-g}{g}^{{\mu }_1{\nu }_1}\cdots g^{{\mu }_b{\nu }_b}{\varphi }^a{\partial }_{{\mu }_1}\varphi {\partial }_{{\nu }_1}\varphi \cdots {\partial }_{{\mu }_b}\varphi {\partial }_{{\nu }_b}\varphi \hfill \\ \hfill & =\sum _{n=0}^{\infty }\int \prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_i\right)\hfill \\ \hfill & \times {\kappa }^n{(-1)}^b{\lambda }_{(a,b)}{\left(\sqrt{-g}{g}^{{\mu }_1{\nu }_1}\cdots g^{{\mu }_b{\nu }_b}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{\left(p_1\right)}_{{\mu }_1}{\left(p_2\right)}_{{\nu }_1}\cdots {\left(p_{2b-1}\right)}_{{\mu }_b}{\left(p_{2b}\right)}_{{\nu }_b}.\hfill \end{array}\end{array}$$

(70)

Lastly, the Feynman rule corresponding to that action reads

(71)

The permutation terms account for all possible permutations of the scalar field momenta. Let us also note that for the $b=0$ case, the interaction reduces the scalar field potential.

4.2.2. Horndeski $G_3$ Class

The interactions of the $G_3$ class of Horndeski theory represent the simplest non-minimal coupling between the scalar field and gravity. The final expression of this interaction contains a single Christoffel symbol, which allows for the explicit inclusion of a single graviton momentum in the corresponding Feynman rule.

We shall first simplify the expression and make the Christoffel symbol explicit.

Theorem 15.

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \int d^{4}x\sqrt{-g}{G}_3\square \varphi =\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{G}_3{\partial }_{\mu }{\partial }_{\nu }\varphi -\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{\mathrm{\Gamma }}_{\alpha \mu \nu }{G}_3{\partial }_{\beta }\varphi .\end{array}\end{array}$$

(72)

Secondly, $G_3$ should be expanded in a Taylor series:

$$\begin{array}{c}\hfill G_3(\varphi ,X)=\sum _{a=0}^{\infty }\sum _{b=0}^{\infty }{\mathsf{\Theta }}_{(a,b)}{\varphi }^a{X}^b.\end{array}$$

(73)

We can study a single term of the expansion without the loss of generality:

$$\begin{array}{c}\hfill \int d^{4}x\sqrt{-g}{G}_3(\varphi ,X)\square \varphi \to \int d^{4}x\sqrt{-g}{\mathsf{\Theta }}_{(a,b)}{\varphi }^a{X}^b\square \varphi .\end{array}$$

(74)

The following theorem gives the perturbative structure of this interaction.

Theorem 16.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}{\mathsf{\Theta }}_{(a,b)}{\varphi }^a{X}^b\square \varphi \hfill \\ \hfill & =\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{\mathsf{\Theta }}_{(a,b)}{\varphi }^a{X}^b{\partial }_{\mu }{\partial }_{\nu }\varphi -\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{\mathrm{\Gamma }}_{\alpha \mu \nu }{\mathsf{\Theta }}_{(a,b)}{\varphi }^a{X}^b{\partial }_{\beta }\varphi \hfill \\ \hfill & =\sum _{n=0}^{\infty }\int \prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=0}^{a+2b+1}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_i\right)\hfill \\ \hfill & \times {\kappa }^n{(-1)}^{b+1}{\mathsf{\Theta }}_{(a,b)}{\left(\sqrt{-g}{g}^{\mu \nu }{g}^{{\alpha }_1{\beta }_1}\cdots g^{{\alpha }_b{\beta }_b}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_1}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\mu }{\left(p_{2b+1}\right)}_{\nu }\hfill \\ \hfill & +\sum _{n=1}^{\infty }\int \prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=0}^{a+2b+1}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_i\right)\hfill \\ \hfill & \times {\kappa }^n{(-1)}^{b+2}{\mathsf{\Theta }}_{(a,b)}{\left(\sqrt{-g}{g}^{\mu \nu }{g}^{\rho \sigma }{g}^{{\alpha }_1{\beta }_1}\cdots g^{{\alpha }_b{\beta }_b}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{\rho \mu \nu }\right)}^{\lambda {\rho }_n{\sigma }_n}{\left(k_n\right)}_{\lambda }\hfill \\ \hfill & \times {\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_1}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\sigma }.\hfill \end{array}\end{array}$$

(75)

It should be noted that the second part of the expression does not contribute at the ${\kappa }^0$ level, since the Christoffel symbol vanishes on the background.

The corresponding interaction rule is given by the following:

(76)

The permutation terms account for all possible permutations of the scalar field and graviton momenta.

4.2.3. Horndeski $G_4$ Class

Interactions of the $G_4$ class describe more sophisticated coupling of the scalar field to gravity. Similarly to the previous case, it involves Christoffel symbols. In contrast with the previous case, Christoffel symbols are present implicitly in the scalar curvature. Moreover, the interaction is described by two related terms, with the second term cancelling higher derivative contributions to the classical equations of motion.

We assume $G_4$ is smooth enough to admit a Taylor expansion. To ensure that the theory admits general relativity in the corresponding limit, one shall assume that the expansion takes the following form:

$$\begin{array}{c}\hfill G_4=-\frac{2}{{\kappa }^2}+\sum _{a,b}{{\rm Y}}_{(a,b)}{\varphi }^a{X}^b.\end{array}$$

(77)

The first term corresponds to pure general relativity, while the other terms describe the interaction between gravity and the scalar field. We discuss the part corresponding to general relativity in another section, so we do not discuss it here. Without loss of generality, we will use a single term from this expansion:

$$\begin{array}{c}\hfill \int d^{4}x\sqrt{-g}\left[{{\rm Y}}_{(a,b)}{\varphi }^a{X}^{b}R+b{{\rm Y}}_{(a,b)}{\varphi }^a{X}^{b-1}\left[{\left(\square \varphi \right)}^2-{\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^2\right]\right].\end{array}$$

(78)

The following theorems give the structure of this action.

Theorem 17.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}{{\rm Y}}_{(a,b)}{\varphi }^a{X}^{b}R=\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{\partial }_{\alpha }\left[{\mathrm{\Gamma }}_{\beta \mu \nu }-{\mathrm{\Gamma }}_{\mu \beta \nu }\right]{{\rm Y}}_{(a,b)}{\varphi }^a{X}^b\hfill \\ \hfill & +\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }\left[{\mathrm{\Gamma }}_{\alpha \mu \rho }{\mathrm{\Gamma }}_{\beta \nu \sigma }-{\mathrm{\Gamma }}_{\alpha \mu \nu }{\mathrm{\Gamma }}_{\beta \rho \sigma }\right]{{\rm Y}}_{(a,b)}{\varphi }^a{X}^b\hfill \\ \hfill & =\sum _{n=1}^{\infty }\int \prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times {\kappa }^n{(-1)}^{b+1}{{\rm Y}}_{(a,b)}{\left(\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{g}^{{\alpha }_1{\beta }_1}\cdots g^{{\alpha }_b{\beta }_b}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-1}{\sigma }_{n-1}}{\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_1}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2n}\right)}_{{\beta }_b}\hfill \\ \hfill & \times \left[{\left({\mathrm{\Gamma }}_{\beta \mu \nu }\right)}^{\lambda {\rho }_n{\sigma }_n}-{\left({\mathrm{\Gamma }}_{\mu \beta \nu }\right)}^{\lambda {\rho }_n{\sigma }_n}\right]{\left(k_n\right)}_{\alpha }{\left(k_n\right)}_{\lambda }\hfill \\ \hfill & +\sum _{n=2}^{\infty }\int \prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times {\kappa }^n{(-1)}^{b+2}{{\rm Y}}_{(a,b)}{\left(\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }{g}^{{\alpha }_1{\beta }_1}\cdots g^{{\alpha }_b{\beta }_b}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-2}{\sigma }_{n-2}}{\left(p_1\right)}_{\alpha }{\left(p_2\right)}_{\beta }\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{2b}\hfill \\ \hfill & \times \left[{\left({\mathrm{\Gamma }}_{\alpha \mu \rho }\right)}^{{\lambda }_1{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{\beta \nu \sigma }\right)}^{{\lambda }_2{\rho }_n{\sigma }_n}-{\left({\mathrm{\Gamma }}_{\alpha \mu \nu }\right)}^{{\lambda }_1{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{\beta \rho \sigma }\right)}^{{\lambda }_2{\rho }_n{\sigma }_n}\right]{\left(k_{n-1}\right)}_{{\lambda }_1}{\left(k_n\right)}_{{\lambda }_2}.\hfill \end{array}\end{array}$$

(79)

The expression consists of two parts. The first part contributes at the $\mathcal{O}\left({\kappa }^1\right)$ level and describes the interaction between a few scalars and a single graviton. The second part does not contribute to $\mathcal{O}\left({\kappa }^1\right)$, and describes the interaction between two or more gravitons with a few scalars.

The following theorem describes the structure of the term containing $G_{4,X}$.

Theorem 18.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}{\varphi }^a{X}^{b-1}\left[{\left(\square \varphi \right)}^2-{\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^2\right]\hfill \\ \hfill & =\int d^{4}x\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }-\frac{1}{2}\left(g^{\mu \alpha }{g}^{\nu \beta }+g^{\mu \beta }{g}^{\nu \alpha }\right)\right]{\varphi }^a{X}^b{\partial }_{\mu }{\partial }_{\nu }\varphi {\partial }_{\alpha }{\partial }_{\beta }\varphi \hfill \\ \hfill & +\int d^{4}x\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }-\frac{1}{2}\left(g^{\mu \alpha }{g}^{\nu \beta }+g^{\mu \beta }{g}^{\nu \alpha }\right)\right]g^{\rho \sigma }{\mathrm{\Gamma }}_{\rho \mu \nu }{\varphi }^a{X}^b{\partial }_{\sigma }\varphi {\partial }_{\alpha }{\partial }_{\beta }\varphi \hfill \\ \hfill & +\int d^{4}x\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }-\frac{1}{2}\left(g^{\mu \alpha }{g}^{\nu \beta }+g^{\mu \beta }{g}^{\nu \alpha }\right)\right]g^{\rho \sigma }{g}^{\lambda \tau }{\mathrm{\Gamma }}_{\rho \mu \nu }{\mathrm{\Gamma }}_{\lambda \alpha \beta }{\varphi }^a{X}^b{\partial }_{\sigma }\varphi {\partial }_{\tau }\varphi .\hfill \end{array}\end{array}$$

(80)

In turn, the Fourier structure of this contribution is given by the following:

Theorem 19.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}b{{\rm Y}}_{(a,b)}{\varphi }^a{X}^{b-1}\left[{\left(\square \varphi \right)}^2-{\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^2\right]\hfill \\ \hfill & =\sum _{n=0}^{\infty }\int \prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times {\kappa }^n{(-1)}^{b-1}b{{\rm Y}}_{(a,b)}{\left(\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }-\frac{1}{2}\left(g^{\mu \alpha }{g}^{\nu \beta }+g^{\mu \beta }{g}^{\nu \alpha }\right)\right]g^{{\alpha }_1{\beta }_1}\cdots g^{{\alpha }_b{\beta }_b}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\hfill \\ \hfill & \times {\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_1}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\mu }{\left(p_{2b+1}\right)}_{\nu }{\left(p_{2b+2}\right)}_{\alpha }{\left(p_{2b+2}\right)}_{\beta }\hfill \\ \hfill & +\sum _{n=1}^{\infty }\int \prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times {\kappa }^n{(-1)}^{b}2b{{\rm Y}}_{(a,b)}{\left(\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }-\frac{1}{2}\left(g^{\mu \alpha }{g}^{\nu \beta }+g^{\mu \beta }{g}^{\nu \alpha }\right)\right]g^{\rho \sigma }{g}^{{\alpha }_1{\beta }_1}\cdots g^{{\alpha }_b{\beta }_b}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-1}{\sigma }_{n-1}}\hfill \\ \hfill & \times {\left({\mathrm{\Gamma }}_{\rho \mu \nu }\right)}^{\lambda {\rho }_n{\sigma }_n}{\left(k_n\right)}_{\lambda }{\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_1}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\sigma }{\left(p_{2b+2}\right)}_{\alpha }{\left(p_{2b+2}\right)}_{\beta }\hfill \\ \hfill & +\sum _{n=2}^{\infty }\int \prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times {\kappa }^n{(-1)}^{b}b{{\rm Y}}_{(a,b)}{\left(\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }-\frac{1}{2}\left(g^{\mu \alpha }{g}^{\nu \beta }+g^{\mu \beta }{g}^{\nu \alpha }\right)\right]g^{\rho \sigma }{g}^{\lambda \tau }{g}^{{\alpha }_1{\beta }_1}\cdots g^{{\alpha }_b{\beta }_b}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-2}{\sigma }_{n-2}}\hfill \\ \hfill & \times {\left({\mathrm{\Gamma }}_{\rho \mu \nu }\right)}^{{\lambda }_1{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{\lambda \alpha \beta }\right)}^{{\lambda }_2{\rho }_n{\sigma }_n}{\left(k_{n-1}\right)}_{{\lambda }_1}{\left(k_n\right)}_{{\lambda }_2}{\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_1}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\sigma }{\left(p_{2b+2}\right)}_{\rho }.\hfill \end{array}\end{array}$$

(81)

The first part of the expression contributes to the background at order $\mathcal{O}\left({\kappa }^0\right)$. The second part contributes to $\mathcal{O}\left({\kappa }^1\right)$ and describes the interaction of a single graviton with a few scalars. The last part contributes to order $\mathcal{O}\left({\kappa }^2\right)$, and its leading contribution represents the interaction of two gravitons with a few scalars.

The following sophisticated expression gives the resulting expression for the interaction rule.

(82)

Let us note that the second and last terms in the square brackets contribute only to vertices with two or more gravitons.

4.2.4. Horndeski $G_5$ Class

The final class of Horndeski interaction also describes a non-minimal coupling between the scalar field and gravity. This interaction involves the most derivatives, which still do not result in the higher-order classical field equations. We study this term similar to the others and make its structure explicit.

The following theorem describes the structure of the non-minimal coupling. It does not involve any integration by parts and relies purely on calculating derivatives.

Theorem 20.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill G^{\mu \nu }{\nabla }_{\mu }{\nabla }_{\nu }\varphi =& -\frac{1}{2}\left[g^{\mu \nu }{g}^{\alpha \beta }-g^{\mu \alpha }{g}^{\nu \beta }-g^{\mu \beta }{g}^{\nu \alpha }\right]g^{\rho \sigma }\left\{{\partial }_{\rho }{\mathrm{\Gamma }}_{\sigma \alpha \beta }-{\partial }_{\alpha }{\mathrm{\Gamma }}_{\rho \beta \sigma }\right\}{\partial }_{\mu }{\partial }_{\nu }\varphi \hfill \\ \hfill & +\frac{1}{2}\left[g^{\mu \nu }{g}^{\alpha \beta }-g^{\mu \alpha }{g}^{\nu \beta }-g^{\mu \beta }{g}^{\nu \alpha }\right]g^{\rho \sigma }{g}^{\lambda \tau }\left\{{\partial }_{\rho }{\mathrm{\Gamma }}_{\sigma \alpha \beta }-{\partial }_{\alpha }{\mathrm{\Gamma }}_{\rho \beta \sigma }\right\}{\mathrm{\Gamma }}_{\lambda \mu \nu }{\partial }_{\tau }\varphi \hfill \\ \hfill & -\frac{1}{2}\left[g^{\mu \nu }{g}^{\alpha \beta }-g^{\mu \alpha }{g}^{\nu \beta }-g^{\mu \beta }{g}^{\nu \alpha }\right]g^{\rho \sigma }{g}^{\lambda \tau }\left\{{\mathrm{\Gamma }}_{\rho \alpha \lambda }{\mathrm{\Gamma }}_{\sigma \beta \tau }-{\mathrm{\Gamma }}_{\rho \alpha \beta }{\mathrm{\Gamma }}_{\sigma \lambda \tau }\right\}{\partial }_{\mu }{\partial }_{\nu }\varphi \hfill \\ \hfill & +\frac{1}{2}\left[g^{\mu \nu }{g}^{\alpha \beta }-g^{\mu \alpha }{g}^{\nu \beta }-g^{\mu \beta }{g}^{\nu \alpha }\right]g^{\rho \sigma }{g}^{\lambda \tau }{g}^{ϵ\omega }\left\{{\mathrm{\Gamma }}_{\rho \alpha \lambda }{\mathrm{\Gamma }}_{\sigma \beta \tau }-{\mathrm{\Gamma }}_{\rho \alpha \beta }{\mathrm{\Gamma }}_{\sigma \lambda \tau }\right\}{\mathrm{\Gamma }}_{ϵ\mu \nu }{\partial }_{\omega }\varphi .\hfill \end{array}\end{array}$$

(83)

The coupling function $G_5=G_5(\varphi ,X)$ shall be smooth enough to admit a power series expansion

$$\begin{array}{c}\hfill G_5=\sum _{a,b}{\mathsf{\Psi }}_{(a,b)}{\varphi }^a{X}^b.\end{array}$$

(84)

Similarly to the previous cases, we can study the following coupling function without the loss of generality:

$$\begin{array}{c}\hfill G_5={\mathsf{\Psi }}_{(a,b)}{\varphi }^a{X}^b.\end{array}$$

(85)

This results in the following theorem describing the Fourier structure of the first part of the interaction.

Theorem 21.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}{G}_5{G}^{\mu \nu }{\nabla }_{\mu }{\nabla }_{\nu }\varphi \hfill \\ \hfill & =\sum _{n=1}^{\infty }\int \prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b+1}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times \frac{1}{2}{\kappa }^n{(-1)}^{b+1}{\mathsf{\Psi }}_{(a,b)}{\left\{\left(g^{\mu \nu }{g}^{\alpha \beta }-g^{\mu \alpha }{g}^{\nu \beta }-g^{\mu \beta }{g}^{\nu \alpha }\right)g^{\rho \sigma }\right\}}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-1}{\sigma }_{n-1}}\hfill \\ \hfill & \times {\left(k_n\right)}_{\lambda }\left\{{\left(k_n\right)}_{\rho }{\left({\mathrm{\Gamma }}_{\sigma \alpha \beta }\right)}^{\lambda {\rho }_n{\sigma }_n}-{\left(k_n\right)}_{\alpha }{\left({\mathrm{\Gamma }}_{\rho \beta \sigma }\right)}^{\lambda {\rho }_n{\sigma }_n}\right\}{\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_1}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\mu }{\left(p_{2b+1}\right)}_{\nu }\hfill \\ \hfill & +\sum _{n=2}^{\infty }\int \prod _{i=2}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b+1}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times \frac{1}{2}{\kappa }^n{(-1)}^b{\mathsf{\Psi }}_{(a,b)}{\left\{\left(g^{\mu \nu }{g}^{\alpha \beta }-g^{\mu \alpha }{g}^{\nu \beta }-g^{\mu \beta }{g}^{\nu \alpha }\right)g^{\rho \sigma }{g}^{\lambda \tau }\right\}}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-2}{\sigma }_{n-2}}{\left(k_{n-1}\right)}_{{\lambda }_1}{\left({\mathrm{\Gamma }}_{\lambda \mu \nu }\right)}^{{\lambda }_1{\rho }_{n-1}{\sigma }_{n-1}}\hfill \\ \hfill & \times {\left(k_n\right)}_{\lambda }\left\{{\left(k_n\right)}_{\rho }{\left({\mathrm{\Gamma }}_{\sigma \alpha \beta }\right)}^{\lambda {\rho }_n{\sigma }_n}-{\left(k_n\right)}_{\alpha }{\left({\mathrm{\Gamma }}_{\rho \beta \sigma }\right)}^{\lambda {\rho }_n{\sigma }_n}\right\}{\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_1}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\tau }\hfill \\ \hfill & +\sum _{n=2}^{\infty }\int \prod _{i=2}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b+1}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times \frac{1}{2}{\kappa }^n{(-1)}^{b+1}{\mathsf{\Psi }}_{(a,b)}{\left\{\left(g^{\mu \nu }{g}^{\alpha \beta }-g^{\mu \alpha }{g}^{\nu \beta }-g^{\mu \beta }{g}^{\nu \alpha }\right)g^{\rho \sigma }{g}^{\lambda \tau }\right\}}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-2}{\sigma }_{n-2}}{\left(k_{n-1}\right)}_{{\lambda }_1}{\left(k_n\right)}_{{\lambda }_2}\hfill \\ \hfill & \times \left\{{\left({\mathrm{\Gamma }}_{\rho \alpha \lambda }\right)}^{{\lambda }_1{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{\sigma \beta \tau }\right)}^{{\lambda }_2{\rho }_n{\sigma }_n}-{\left({\mathrm{\Gamma }}_{\rho \alpha \beta }\right)}^{{\lambda }_1{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{\sigma \lambda \tau }\right)}^{{\lambda }_2{\rho }_n{\sigma }_n}\right\}{\left(p_1\right)}_{{\alpha }_1}\cdots {\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\mu }{\left(p_{2b+1}\right)}_{\nu }\hfill \\ \hfill & +\sum _{n=3}^{\infty }\int \prod _{i=3}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b+1}\frac{d^4{p}_j}{{\left(2\pi \right)}^4}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times \frac{1}{2}{\kappa }^n{(-1)}^b{\mathsf{\Psi }}_{(a,b)}{\left\{\left(g^{\mu \nu }{g}^{\alpha \beta }-g^{\mu \alpha }{g}^{\nu \beta }-g^{\mu \beta }{g}^{\nu \alpha }\right)g^{\rho \sigma }{g}^{\lambda \tau }{g}^{ϵ\omega }\right\}}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-3}{\sigma }_{n-3}}{\left(k_{n-2}\right)}_{{\lambda }_1}{\left(k_{n-1}\right)}_{{\lambda }_2}{\left(k_n\right)}_{{\lambda }_3}\hfill \\ \hfill & \times \left\{{\left({\mathrm{\Gamma }}_{\rho \alpha \lambda }\right)}^{{\lambda }_2{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{\sigma \beta \tau }\right)}^{{\lambda }_3{\rho }_n{\sigma }_n}-{\left({\mathrm{\Gamma }}_{\rho \alpha \beta }\right)}^{{\lambda }_2{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{\sigma \lambda \tau }\right)}^{{\lambda }_3{\rho }_n{\sigma }_n}\right\}{\left({\mathrm{\Gamma }}_{ϵ\mu \nu }\right)}^{{\lambda }_1{\rho }_{n-2}{\sigma }_{n-2}}{\left(p_1\right)}_{{\alpha }_1}\cdots {\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\omega }.\hfill \end{array}\end{array}$$

(86)

First, we describe the structure involving scalar field derivatives to proceed with the second part of the interaction.

Theorem 22.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {(\square \varphi )}^3-3\square \varphi {\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^2+2{\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^3\hfill \\ \hfill & =\left[g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }-3g^{\mu \nu }{g}^{\alpha \rho }{g}^{\beta \sigma }+2g^{\nu \alpha }{g}^{\beta \rho }{g}^{\sigma \mu }\right]{\partial }_{\mu }{\partial }_{\nu }\varphi {\partial }_{\alpha }{\partial }_{\beta }\varphi {\partial }_{\rho }{\partial }_{\sigma }\varphi \hfill \\ \hfill & +3\left[g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }-3g^{\mu \nu }{g}^{\alpha \rho }{g}^{\beta \sigma }+2g^{\nu \alpha }{g}^{\beta \rho }{g}^{\sigma \mu }\right]g^{\omega \tau }{\mathrm{\Gamma }}_{\omega \mu \nu }{\partial }_{\tau }\varphi {\partial }_{\alpha }{\partial }_{\beta }\varphi {\partial }_{\rho }{\partial }_{\sigma }\varphi \hfill \\ \hfill & +3\left[g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }-3g^{\mu \nu }{g}^{\alpha \rho }{g}^{\beta \sigma }+2g^{\nu \alpha }{g}^{\beta \rho }{g}^{\sigma \mu }\right]g^{{\omega }_1{\tau }_1}{g}^{{\omega }_2{\tau }_2}{\mathrm{\Gamma }}_{{\omega }_1\mu \nu }{\mathrm{\Gamma }}_{{\omega }_2\alpha \beta }{\partial }_{{\tau }_1}\varphi {\partial }_{{\tau }_2}\varphi {\partial }_{\rho }{\partial }_{\sigma }\varphi \hfill \\ \hfill & +\left[g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }-3g^{\mu \nu }{g}^{\alpha \rho }{g}^{\beta \sigma }+2g^{\nu \alpha }{g}^{\beta \rho }{g}^{\sigma \mu }\right]g^{{\omega }_1{\tau }_1}{g}^{{\omega }_2{\tau }_2}{g}^{{\omega }_3{\tau }_3}{\mathrm{\Gamma }}_{{\omega }_1\mu \nu }{\mathrm{\Gamma }}_{{\omega }_2\alpha \beta }{\mathrm{\Gamma }}_{{\omega }_3\rho \sigma }{\partial }_{{\tau }_1}\varphi {\partial }_{{\tau }_2}\varphi {\partial }_{{\tau }_3}\varphi .\hfill \end{array}\end{array}$$

(87)

The following theorem gives the perturbative structure of the second part of the interaction.

Theorem 23.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}\left[-\frac{1}{6}{G}_{5,X}\left\{{(\square \varphi )}^3-3\square \varphi {\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^2+2{\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)}^3\right\}\right]\hfill \\ \hfill & =\sum _{n=0}^{\infty }\int \prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b+1}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times \frac{1}{6}{\kappa }^n{(-1)}^{b}b{\mathsf{\Psi }}_{(a,b)}{\left\{\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }-3g^{\mu \nu }{g}^{\alpha \rho }{g}^{\beta \sigma }+2g^{\nu \alpha }{g}^{\beta \rho }{g}^{\sigma \mu }\right]\right\}}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\hfill \\ \hfill & \times {\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_2}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(p_{2b+1}\right)}_{\mu }{\left(p_{2b+1}\right)}_{\nu }{\left(p_{2b+2}\right)}_{\alpha }{\left(p_{2b+2}\right)}_{\beta }{\left(p_{2b+3}\right)}_{\rho }{\left(p_{2b+3}\right)}_{\sigma }\hfill \\ \hfill & +\sum _{n=1}^{\infty }\int \prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b+1}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times \frac{1}{3}{\kappa }^n{(-1)}^{b}b{\mathsf{\Psi }}_{(a,b)}{\left\{\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }-3g^{\mu \nu }{g}^{\alpha \rho }{g}^{\beta \sigma }+2g^{\nu \alpha }{g}^{\beta \rho }{g}^{\sigma \mu }\right]g^{\omega \tau }\right\}}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-1}{\sigma }_{n-1}}\hfill \\ \hfill & \times {\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_2}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(k_n\right)}_{\lambda }{\left({\mathrm{\Gamma }}_{\omega \mu \nu }\right)}^{\lambda {\rho }_n{\sigma }_n}{\left(p_{2b+1}\right)}_{\tau }{\left(p_{2b+2}\right)}_{\alpha }{\left(p_{2b+2}\right)}_{\beta }{\left(p_{2b+3}\right)}_{\rho }{\left(p_{2b+3}\right)}_{\sigma }\hfill \end{array}\end{array}$$

(88)

$$\begin{array}{cc}\hfill & +\sum _{n=2}^{\infty }\int \prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b+1}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times \frac{1}{3}{\kappa }^n{(-1)}^{b}b{\mathsf{\Psi }}_{(a,b)}{\left\{\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }-3g^{\mu \nu }{g}^{\alpha \rho }{g}^{\beta \sigma }+2g^{\nu \alpha }{g}^{\beta \rho }{g}^{\sigma \mu }\right]g^{{\omega }_1{\tau }_1}{g}^{{\omega }_2{\tau }_2}\right\}}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-2}{\sigma }_{n-2}}\hfill \\ \hfill & \times {\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_2}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(k_{n-1}\right)}_{{\lambda }_1}{\left(k_n\right)}_{{\lambda }_2}{\left({\mathrm{\Gamma }}_{{\omega }_1\mu \nu }\right)}^{{\lambda }_1{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{{\omega }_2\alpha \beta }\right)}^{{\lambda }_2{\rho }_n{\sigma }_n}\hfill \\ \hfill & \times {\left(p_{2b+1}\right)}_{{\tau }_1}{\left(p_{2b+2}\right)}_{{\tau }_2}{\left(p_{2b+3}\right)}_{\rho }{\left(p_{2b+3}\right)}_{\sigma }\hfill \\ \hfill & +\sum _{n=3}^{\infty }\int \prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{h}_{{\rho }_i{\sigma }_i}\left(k_i\right)\prod _{j=1}^{a+2b+1}\varphi \left(p_j\right){\left(2\pi \right)}^4\delta \left(\sum k_i+\sum p_j\right)\hfill \\ \hfill & \times \frac{1}{6}{\kappa }^n{(-1)}^{b}b{\mathsf{\Psi }}_{(a,b)}{\left\{\sqrt{-g}\left[g^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }-3g^{\mu \nu }{g}^{\alpha \rho }{g}^{\beta \sigma }+2g^{\nu \alpha }{g}^{\beta \rho }{g}^{\sigma \mu }\right]g^{{\omega }_1{\tau }_1}{g}^{{\omega }_2{\tau }_2}{g}^{{\omega }_3{\tau }_3}\right\}}^{{\rho }_1{\sigma }_1\cdots {\rho }_{n-3}{\sigma }_{n-3}}\hfill \\ \hfill & \times {\left(p_1\right)}_{{\alpha }_1}{\left(p_2\right)}_{{\beta }_2}\cdots {\left(p_{2b-1}\right)}_{{\alpha }_b}{\left(p_{2b}\right)}_{{\beta }_b}{\left(k_{n-2}\right)}_{{\lambda }_1}{\left(k_{n-1}\right)}_{{\lambda }_2}{\left(k_n\right)}_{{\lambda }_3}\hfill \\ \hfill & \times {\left({\mathrm{\Gamma }}_{{\omega }_1\mu \nu }\right)}^{{\lambda }_1{\rho }_{n-2}{\sigma }_{n-2}}{\left({\mathrm{\Gamma }}_{{\omega }_2\alpha \beta }\right)}^{{\lambda }_2{\rho }_{n-1}{\sigma }_{n-1}}{\left({\mathrm{\Gamma }}_{{\omega }_3\rho \sigma }\right)}^{{\lambda }_3{\rho }_n{\sigma }_n}{\left(p_{2b+1}\right)}_{{\tau }_1}{\left(p_{2b+2}\right)}_{{\tau }_2}{\left(p_{2b+3}\right)}_{{\tau }_3}.\hfill \end{array}$$

The theorems explain the entire structure of the $G_5$ interaction class. However, we will not explicitly express the corresponding interaction vertex. Because the expression is exceptionally long, it would not be constructive.

4.3. Dirac Spinors

The standard way to describe fermions in curved spacetime uses vierbein and $\gamma$ matrices [87]. Firstly, $\gamma$ matrices are introduced to construct a representation of the Lorentz algebra. Dirac matrices satisfy the following relations:

$$\begin{array}{cccc}\hfill \{{\gamma }_m,{\gamma }_n\}\stackrel{\mathrm{note}}{=}{\gamma }_m{\gamma }_n+{\gamma }_n{\gamma }_m& =2{\eta }_{mn},\hfill & \hfill {\gamma }^0{\gamma }^m{\gamma }^0& ={\left({\gamma }^m\right)}^{+}.\hfill \end{array}$$

(89)

They form the following representation of the Lorentz algebra:

$$\begin{array}{cccc}\hfill S_{mn}& \stackrel{def}{=}\frac{i}{4}[{\gamma }_m,{\gamma }_n],\hfill & \hfill [S_{mn},S_{ab}]& =-i\left({\eta }_{ma}{S}_{nb}-{\eta }_{mb}{S}_{na}+{\eta }_{nb}{S}_{ma}-{\eta }_{na}{S}_{mb}\right).\hfill \end{array}$$

(90)

The standard vector representation of the Lorentz algebra is given by $J_{mn}$ matrices:

$$\begin{array}{cccc}\hfill {\left(J_{mn}\right)}_{ab}=& i({\eta }_{ma}{\eta }_{nb}-{\eta }_{mb}{\eta }_{na}),\hfill & \hfill [J_{mn},J_{ab}]& =-i\left({\eta }_{ma}{J}_{nb}-{\eta }_{mb}{J}_{na}+{\eta }_{nb}{J}_{ma}-{\eta }_{na}{J}_{mb}\right).\hfill \end{array}$$

(91)

The following relation connects these representations:

$$\begin{array}{c}\hfill [{\gamma }_a,S_{mn}]={\left(J_{mn}\right)}_{ab}{\gamma }^b.\end{array}$$

(92)

One can define Dirac spinors $\psi$ and $\overline{\psi }={\psi }^{+}{\gamma }^0$ subject to Lorentz group action

$$\begin{array}{cccc}\hfill \delta \psi & =\frac{i}{2}{\omega }^{ab}{S}_{ab}\psi ,\hfill & \hfill \delta \overline{\psi }& =-\frac{i}{2}{\omega }^{ab}\overline{\psi }{S}_{ab},\hfill \end{array}$$

(93)

and to construct the well-known Lorentz invariant action in a flat spacetime

$$\begin{array}{c}\hfill {\mathcal{A}}_{s=\frac{1}{2},m\ne 0}=\int d^{4}x\left[\frac{i}{2}\left(\overline{\psi }{\gamma }^m{\partial }_m\psi -{\partial }_m\overline{\psi }{\gamma }^m\psi \right)-m\overline{\psi }\psi \right].\end{array}$$

(94)

The action is generalised for curved spacetime via vierbein ${\mathfrak{e}}_m{}^{\mu }$. In a curved spacetime, one relates an arbitrary frame with a local inertial frame via vierbein ${\mathfrak{e}}_m{}^{\mu }$. Its Latin indices are subjected to the Lorentz transformation; Greek indices are subjected to the general coordinate transformations. In turn, vierbein satisfies the following normalisation condition:

$$\begin{array}{cccc}\hfill {\mathfrak{e}}_m{}^{\mu }{\mathfrak{e}}_n{}^{\nu }{g}_{\mu \nu }=& {\eta }_{mn},\hfill & \hfill {\mathfrak{e}}_m{}^{\mu }{\mathfrak{e}}_n{}^{\nu }{\eta }^{mn}=& g^{\mu \nu }.\hfill \end{array}$$

(95)

They define an anti-symmetric spin-connection ${\left({\mathrm{\Gamma }}_{\mu }\right)}_{ab}=-{\left({\mathrm{\Gamma }}_{\mu }\right)}_{ba}$ that is related to the Christoffel symbols.

$$\begin{array}{c}\hfill {\left({\mathrm{\Gamma }}_{\mu }\right)}_{ab}={\mathfrak{e}}_{a\alpha }{\mathfrak{e}}_b{}^{\beta }{\mathrm{\Gamma }}_{\mu \beta }^{\alpha }+{\mathfrak{e}}_{a\sigma }{\partial }_{\mu }{\mathfrak{e}}_b{}^{\sigma }.\end{array}$$

(96)

Since the vierbein ${\mathfrak{e}}_m{}^{\mu }$ connects Lorentz transformations and general coordinate transformations, one uses it to manipulate indices. This construction forms a set of $\gamma$ matrices that are subject to the general coordinate transformation group:

$$\begin{array}{c}\hfill {\gamma }^{\mu }={\mathfrak{e}}_m{}^{\mu }{\gamma }^m.\end{array}$$

(97)

Spinor transformations (93) are generalized as follows:

$$\begin{array}{cccc}\hfill \delta \psi & =\frac{i}{2}{\omega }^{\mu \nu }{S}_{\mu \nu }\psi =\frac{i}{2}{\omega }^{ab}\left(x\right)S_{ab}\psi ,\hfill & \hfill \delta \overline{\psi }& =-\frac{i}{2}{\omega }^{\mu \nu }\overline{\psi }{S}_{\mu \nu }=-\frac{i}{2}{\omega }^{ab}\left(x\right)\overline{\psi }{S}_{ab}.\hfill \end{array}$$

(98)

This generalisation promotes transformation (93) to gauge transformations.

Regular derivatives are replaced with covariant derivatives defined as follows:

$$\begin{array}{cc}\hfill {\nabla }_{\mu }\psi ={\partial }_{\mu }\psi -\frac{i}{2}{\left({\mathrm{\Gamma }}_{\mu }\right)}^{ab}{S}_{ab}\psi ,& {\nabla }_{\mu }\overline{\psi }={\partial }_{\mu }\overline{\psi }+\frac{i}{2}{\left({\mathrm{\Gamma }}_{\mu }\right)}^{ab}\overline{\psi }{S}_{ab}.\hfill \end{array}$$

(99)

Covariant derivatives defined in such a way satisfy the following relations:

$$\begin{array}{cccc}\hfill {\nabla }_{\mu }\left(\overline{\psi }\psi \right)& ={\partial }_{\mu }\left(\overline{\psi }\psi \right)\hfill & \hfill {\nabla }_{\mu }\left(\overline{\psi }{\gamma }^{\mu }\psi \right)& ={\partial }_{\mu }\left(\overline{\psi }{\gamma }^{\mu }\psi \right)+{\mathrm{\Gamma }}_{\mu \beta }^{\alpha }\overline{\psi }{\gamma }^{\beta }\psi .\hfill \end{array}$$

(100)

This construction produces the following generalisation of the Dirac action

$$\begin{array}{c}\hfill {\mathcal{A}}_{s=1/2}=\int d^{4}x\left[\sqrt{-g}{\mathfrak{e}}_m{}^{\mu }\frac{1}{2}\left(i\overline{\psi }{\gamma }^m{\nabla }_{\mu }\psi -i{\nabla }_{\mu }\overline{\psi }{\gamma }^m\psi \right)-m\sqrt{-g}\overline{\psi }\psi \right].\end{array}$$

(101)

Here, m is the fermion mass, ${\mathfrak{e}}_m{}^{\mu }$ is the vierbein, and ∇ is the fermionic covariant derivative. We shall note that we omit a few steps in deriving this action. Namely, the action does not contain either Christoffel or spin connection. The reason for this is discussed in detail in [24].

The following theorem describes the perturbative structure of the Dirac action.

Theorem 24.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{s=1/2}=& \int d^{4}x\left[\sqrt{-g}{\mathfrak{e}}_m{}^{\mu }\frac{1}{2}\left(i\overline{\psi }{\gamma }^m{\partial }_{\mu }\psi -i{\partial }_{\mu }\overline{\psi }{\gamma }^m\psi \right)-m\sqrt{-g}\overline{\psi }\psi \right]\hfill \\ \hfill =& \sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times {\kappa }^n\overline{\psi }\left(p_2\right)\left[{\left(\sqrt{-g}{\mathfrak{e}}_m{}^{\mu }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\frac{1}{2}{(p_1-p_2)}_{\mu }{\gamma }^m-{\left(\sqrt{-g}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}m\right]\psi \left(p_1\right).\hfill \end{array}\end{array}$$

(102)

The background contribution describes the standard fermion propagator:

(103)

The other terms describe the following interaction rules:

(104)

This expression also holds in the $SU\left(N\right)$ Yang–Mills model considered below.

4.4. Proca Field

For quantum field theory, the existence of a vector field mass is crucial. A vector field with zero mass admits a gauge symmetry, which means gauge fixing must be performed. On the other hand, a vector field with a non-vanishing mass, the Proca field, has no gauge symmetry, and the gauge fixing is not an issue. We begin the discussion with the Proca case for the sake of simplicity.

The action describing a Proca field reads:

$$\begin{array}{cc}\hfill {\mathcal{A}}_{s=1,m}=& \int d^{4}x\sqrt{-g}\left[-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }+\frac{m^2}{2}{A}_{\mu }{A}^{\mu }\right].\hfill \end{array}$$

(105)

Here, $F_{\mu \nu }={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu }$ is the field tensor, m is the vector field mass. The action admits the following factorisation:

$$\begin{array}{cc}\hfill {\mathcal{A}}_{s=1,m}=& \int d^{4}x\left[-\frac{1}{4}\sqrt{-g}{g}^{\mu \alpha }{g}^{\nu \beta }{F}_{\mu \nu }{F}^{\alpha \beta }+\sqrt{-g}{g}^{\mu \nu }\frac{m^2}{2}{A}_{\mu }{A}_{\nu }\right].\hfill \end{array}$$

(106)

The following theorem describes the Fourier structure of the action.

Theorem 25.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{s=1,m}=& \sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times {\kappa }^n[\frac{1}{4}{\left(\sqrt{-g}{g}^{\mu \alpha }{g}^{\nu \beta }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{\left(p_1\right)}_{{\mu }_1}{\left(p_2\right)}_{{\mu }_2}{\left(F_{\mu \nu }\right)}^{{\mu }_1{\lambda }_1}{\left(F_{\alpha \beta }\right)}^{{\mu }_2{\lambda }_2}\hfill \\ \hfill & +\frac{m^2}{2}{\left(\sqrt{-g}{g}^{{\lambda }_1{\lambda }_2}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}]A_{{\lambda }_1}\left(p_1\right)A_{{\lambda }_2}\left(p_2\right).\hfill \end{array}\end{array}$$

(107)

Here, we introduced the following notations:

$$\begin{array}{cccc}\hfill F_{\mu \nu }& =-i{p}_{\sigma }{\left(F_{\mu \nu }\right)}^{\sigma \lambda }{A}_{\lambda }\left(p\right),\hfill & \hfill {\left(F_{\mu \nu }\right)}^{\sigma \lambda }& \stackrel{def}{=}{\delta }_{\mu }^{\sigma }{\delta }_{\nu }^{\lambda }-{\delta }_{\nu }^{\sigma }{\delta }_{\mu }^{\lambda }.\hfill \end{array}$$

(108)

This expression generates the Proca propagator:

(109)

The expression describing the interaction rules between gravitons and the Proca field kinetic energy is as follows:

(110)

4.5. Vector Field

Before discussing the massless case, we shall recall the Faddeev–Popov prescription [39]. The following generating functional describes a quantum massless vector field:

$$\begin{array}{c}\hfill \mathcal{Z}=\int \mathcal{D}\left[A\right]exp\left[i\mathcal{A}\left[A\right]\right].\end{array}$$

(111)

The integration space includes all conceivable configurations of the vector field. Firstly, we shall add a new term to the microscopic action:

$$\begin{array}{c}\hfill \mathcal{Z}=\int \mathcal{D}\left[A\right]exp\left[i\mathcal{A}\left[A\right]\right]\int \mathcal{D}\left[\omega \right]exp\left[\frac{i}{2}ϵ{\omega }^2\right]=\int \mathcal{D}\left[A\right]\mathcal{D}\left[\omega \right]exp\left[i\mathcal{A}+\frac{i}{2}ϵ{\omega }^2\right].\end{array}$$

(112)

Here, $\omega$ is an arbitrary scalar, and $ϵ$ is a free gauge fixing parameter. The new contribution is a Gauss-like integral, so its introduction merely changes the omitted normalisation factor.

Secondly, we split the integration volume:

$$\begin{array}{c}\hfill \int \mathcal{D}\left[A\right]=\int \mathcal{D}\left[\zeta \right]\int \mathcal{D}\left[\mathbb{A}\right]\delta \left(\mathcal{G}-\omega \right)det\Delta .\end{array}$$

(113)

In this expression, $\mathcal{G}$ is the gauge fixing condition, $\zeta$ is the gauge transformation parameter, and the new field variable $\mathbb{A}$ and the field variable A are related as follows:

$$\begin{array}{c}\hfill A_{\mu }={\mathbb{A}}_{\mu }+{\partial }_{\mu }\zeta .\end{array}$$

(114)

The integration over $\mathbb{A}$ is performed over all conceivable fields. Because of the $\delta$ function, only a single representative from each class of physically equivalent potentials contributes to the integral. Lastly, the Faddeev–Popov determinant $det\Delta$ preserves the invariance of the integration measure. The corresponding differential operator $\Delta$ is defined as follows:

$$\begin{array}{c}\hfill \Delta \stackrel{def}{=}\frac{\delta \mathcal{G}}{\delta \zeta }.\end{array}$$

(115)

Finally, we perform integrations and obtain the following expression for the generating functional:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{Z}& =\int \mathcal{D}\left[\mathbb{A}\right]\mathcal{D}\left[\omega \right]\mathcal{D}\left[\zeta \right]\left(det\Delta \right)\delta \left(\mathcal{G}-\omega \right)exp\left[i\mathcal{A}+\frac{i}{2}ϵ{\omega }^2\right]\hfill \\ \hfill & =\int \mathcal{D}\left[\mathbb{A}\right]\left(det\Delta \right)exp\left[i\mathcal{A}+\frac{i}{2}ϵ{\mathcal{G}}^2\right]\hfill \\ \hfill & =\int \mathcal{D}\left[c\right]\mathcal{D}\left[\overline{c}\right]\mathcal{D}\left[\mathbb{A}\right]exp\left[i\overline{c}\Delta c+i\mathcal{A}+\frac{i}{2}ϵ{\mathcal{G}}^2\right].\hfill \end{array}\end{array}$$

(116)

Here, $\overline{c}$ and c are scalar anticommuting Faddeev–Popov ghosts introduced to account for the Faddeev–Popov determinant. We omit the integration over the gauge parameter $\zeta$, as it is irrelevant due to the normalisation.

This prescription generates a functional suitable for treating gauge models. We chose the standard Lorentz gauge fixing condition for simplicity. In a curved spacetime, it reads:

$$\begin{array}{c}\hfill g^{\mu \nu }{\nabla }_{\mu }{A}_{\nu }=0\leftrightarrow g^{\mu \nu }{\partial }_{\mu }{A}_{\nu }-g^{\mu \nu }{\mathrm{\Gamma }}_{\mu \nu }^{\sigma }{A}_{\sigma }=0.\end{array}$$

(117)

The following theorem describes the perturbative structure of the gauge invariant part of the action.

Theorem 26.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{s=1,m_{\mathrm{v}}=0}=& \int d^{4}x\sqrt{-g}\left[-\frac{1}{4}{g}^{\mu \alpha }{g}^{\nu \beta }{F}_{\mu \nu }{F}_{\alpha \beta }\right]\hfill \\ \hfill =& \sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times {\kappa }^n\left[\frac{1}{4}{\left(\sqrt{-g}{g}^{\mu \alpha }{g}^{\nu \beta }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{\left(p_1\right)}_{{\mu }_1}{\left(p_2\right)}_{{\mu }_2}{\left(F_{\mu \nu }\right)}^{{\mu }_1{\lambda }_1}{\left(F_{\alpha \beta }\right)}^{{\mu }_2{\lambda }_2}\right]A_{{\lambda }_1}\left(p_1\right)A_{{\lambda }_2}\left(p_2\right).\hfill \end{array}\end{array}$$

(118)

This expression matches the expression for the Proca field with $m=0$.

The structure of the gauge fixing term is more sophisticated.

Theorem 27.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{\mathrm{gf}}=& \int d^{4}x\sqrt{-g}\left[\frac{ϵ}{2}{\nabla }_{{\lambda }_1}{A}^{{\lambda }_1}{\nabla }_{{\lambda }_2}{A}^{{\lambda }_2}\right]\hfill \\ \hfill =& \frac{ϵ}{2}\int d^{4}x\left(\sqrt{-g}{g}^{{\sigma }_1{\lambda }_1}{g}^{{\sigma }_2{\lambda }_2}\right){\partial }_{{\sigma }_1}{A}_{{\lambda }_1}{\partial }_{{\sigma }_2}{A}_{{\lambda }_2}-ϵ\int d^{4}x\left(\sqrt{-g}{g}^{\mu \nu }{g}^{{\sigma }_1{\lambda }_1}{g}^{{\sigma }_2{\lambda }_2}\right){\mathrm{\Gamma }}_{{\sigma }_1\mu \nu }{A}_{{\lambda }_1}{\partial }_{{\sigma }_2}{A}_{{\lambda }_2}\hfill \\ \hfill & +\frac{ϵ}{2}\int d^{4}x\left(\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{g}^{{\sigma }_1{\lambda }_1}{g}^{{\sigma }_2{\lambda }_2}\right){\mathrm{\Gamma }}_{{\sigma }_1\mu \nu }{\mathrm{\Gamma }}_{{\sigma }_2\alpha \beta }{A}_{{\lambda }_1}{A}_{{\lambda }_2}.\hfill \end{array}\end{array}$$

(119)

The following theorem gives the Fourier structure of the gauge fixing term.

Theorem 28.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{\mathrm{gf}}=& \sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)A_{{\lambda }_1}\left(p_1\right)A_{{\lambda }_2}\left(p_2\right)\hfill \\ \hfill & \times {\kappa }^n{\left(\sqrt{-g}{g}^{{\mu }_1{\lambda }_1}{g}^{{\mu }_2{\lambda }_2}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\left[-\frac{ϵ}{2}{\left(p_1\right)}_{{\mu }_1}{\left(p_2\right)}_{{\mu }_2}\right]\hfill \\ \hfill +& \sum _{n=1}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)A_{{\lambda }_1}\left(p_1\right)A_{{\lambda }_2}\left(p_2\right)\hfill \\ \hfill & \times {\kappa }^n{\left(\sqrt{-g}{g}^{\mu \nu }{g}^{{\mu }_1{\lambda }_1}{g}^{{\mu }_2{\lambda }_2}\right)}^{{\rho }_2{\sigma }_2\cdots {\rho }_n{\sigma }_n}\left[ϵ{\left({\mathrm{\Gamma }}_{{\mu }_1\mu \nu }\right)}^{\sigma {\rho }_1{\sigma }_1}{\left(k_1\right)}_{\sigma }{\left(p_2\right)}_{{\mu }_2}\right]\hfill \\ \hfill +& \sum _{n=2}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)A_{{\lambda }_1}\left(p_1\right)A_{{\lambda }_2}\left(p_2\right)\hfill \\ \hfill & \times {\kappa }^n{\left(\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{g}^{{\mu }_1{\lambda }_1}{g}^{{\mu }_2{\lambda }_2}\right)}^{{\rho }_3{\sigma }_3\cdots {\rho }_n{\sigma }_n}\left[-\frac{ϵ}{2}{\left(k_1\right)}_{{\tau }_1}{\left(k_2\right)}_{{\tau }_2}{\left({\mathrm{\Gamma }}_{{\mu }_1\mu \nu }\right)}^{{\tau }_1{\rho }_1{\sigma }_1}{\left({\mathrm{\Gamma }}_{{\mu }_2\alpha \beta }\right)}^{{\tau }_2{\rho }_2{\sigma }_2}\right].\hfill \end{array}\end{array}$$

(120)

The following notations for the Christoffel symbols were used.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\mu \alpha \beta }=& \frac{\kappa }{2}\left[{\partial }_{\alpha }{h}_{\beta \mu }+{\partial }_{\beta }{h}_{\alpha \mu }-{\partial }_{\mu }{h}_{\alpha \beta }\right]\iff \kappa (-i)p_{\lambda }{\left({\mathrm{\Gamma }}_{\mu \alpha \beta }\right)}^{\lambda \rho \sigma }{h}_{\rho \sigma }\left(p\right),\hfill \\ \hfill {\left({\mathrm{\Gamma }}_{\mu \alpha \beta }\right)}^{\lambda \rho \sigma }=& \frac{1}{2}\left[{\delta }_{\alpha }^{\lambda }{I}_{\beta \mu }{}^{\rho \sigma }+{\delta }_{\beta }^{\lambda }{I}_{\alpha \mu }{}^{\rho \sigma }-{\delta }_{\mu }^{\lambda }{I}_{\alpha \beta }{}^{\rho \sigma }\right].\hfill \end{array}\end{array}$$

(121)

The background part of this expression corresponds to the standard propagator:

(122)

The corresponding interaction rule reads:

(123)

The dots in this expression represent terms that create symmetry with respect to graviton momenta. The final term only affects vertices with two or more gravitons.

We treat the ghost sector of the theory as follows. The Faddeev–Popov differential operator $\Delta$ reduces to the D’Alamber operator in curved spacetime:

$$\begin{array}{c}\hfill \Delta =\frac{\delta }{\delta \zeta }{\nabla }_{\mu }\left(A^{\mu }+{\nabla }^{\mu }\zeta \right)=g^{\mu \nu }{\nabla }_{\mu }{\nabla }_{\nu }.\end{array}$$

(124)

The ghost part of the functional describes a single massless scalar ghost coupled to gravity:

$$\begin{array}{cc}\hfill {\mathcal{Z}}_{\mathrm{ghost}}=& \int \mathcal{D}\left[c\right]\mathcal{D}\left[\overline{c}\right]exp\left[i\int d^{4}x\sqrt{-g}\left(\overline{c}\square c\right)\right]=\int \mathcal{D}\left[c\right]\mathcal{D}\left[\overline{c}\right]exp\left[-i\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{\nabla }_{\mu }\overline{c}{\nabla }_{\nu }c\right].\hfill \end{array}$$

(125)

The following theorem gives the corresponding perturbative expansion.

Theorem 29.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{\mathrm{ghost}}=& -\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{\partial }_{\mu }\overline{c}{\partial }_{\nu }c\hfill \\ \hfill =& \sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times {\kappa }^n\overline{c}\left(p_1\right)\left[{\left(\sqrt{-g}{g}^{\mu \nu }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{\left(p_1\right)}_{\mu }{\left(p_2\right)}_{\nu }\right]c\left(p_2\right).\hfill \end{array}\end{array}$$

(126)

We recover the standard ghost propagator:

(127)

and in the following interaction rule:

(128)

Let us note again that the discussed ghosts are the standard Faddeev–Popov ones. In the context of gravity, there is one additional reason to account for them. The vertex in a diagram represents the interaction between gravitons and vectors, including physical and non-physical vector field polarisations. The coupling of Faddeev–Popov ghosts to gravity cancels out the energy contribution from non-physical polarisations.

4.6. SU(N) Yang–Mills

Let us turn to the gravitational interaction of the $SU\left(N\right)$ Yang–Mills model. In the flat spacetime, the $SU\left(N\right)$ Yang–Mills model is given by the following action:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \mathcal{A}=\int d^{4}x\left[\overline{\psi }\left(i\widehat{\mathcal{D}}-m\right)\psi -\frac{1}{4}{F}_{\mu \nu }^a{F}^{a\mu \nu }\right]\hfill \\ \hfill & =\int d^{4}x\left[\overline{\psi }(i\widehat{\partial }-m)\psi -\frac{1}{4}{\left(f_{\mu \nu }^a\right)}^2+{\mathit{g}}_{\mathrm{s}}\overline{\psi }\widehat{A}\psi -{\mathit{g}}_{\mathrm{s}}{f}^{abc}{\partial }_{\mu }{A}_{\nu }^a{A}^{b\mu }{A}^{c\nu }-\frac{1}{4}{\mathit{g}}_{\mathrm{s}}^2{f}^{amn}{f}^{aij}(A^m·A^i)(A^n·A^j)\right].\hfill \end{array}\end{array}$$

(129)

The fermion covariant derivative is defined in the standard way:

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mu }\psi ={\partial }_{\mu }\psi -i{\mathit{g}}_{\mathrm{s}}{A}_{\mu }\psi .\end{array}$$

(130)

Field tensor $F_{\mu \nu }$ reads

$$\begin{array}{c}\hfill F_{\mu \nu }={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu }-i{\mathit{g}}_{\mathrm{s}}[A_{\mu },A_{\nu }].\end{array}$$

(131)

The gauge field $A_{\mu }$ takes value in $SU\left(N\right)$ algebra:

$$\begin{array}{c}\hfill A_{\mu }=A_{\mu }^a{T}^a,\end{array}$$

(132)

where $T^a$ are generators. In turn, field tensor components are given by the following:

$$\begin{array}{c}\hfill F_{\mu \nu }^a={\partial }_{\mu }{A}_{\nu }^a-{\partial }_{\nu }{A}_{\mu }^a+{\mathit{g}}_{\mathrm{s}}{f}^{abc}{A}_{\mu }^b{A}_{\nu }^c.\end{array}$$

(133)

Here $f^{abc}$ are the structure constants of the algebra:

$$\begin{array}{c}\hfill [T^a,T^b]=i{f}^{abc}{T}^c.\end{array}$$

(134)

Action (129) is generalised for the curved spacetime case as follows. One replaces all derivatives with covariant derivatives and makes explicit the invariant volume factor. Such a generalisation results in the following action:

$$\begin{array}{c}\hfill \mathcal{A}=\int d^{4}x\sqrt{-g}\left[\overline{\psi }\left(i{\mathfrak{e}}_m{}^{\mu }{\gamma }^m{\mathcal{D}}_{\mu }-m\right)\psi -\frac{1}{4}{F}_{\mu \nu }^a{F}^{a\mu \nu }\right].\end{array}$$

(135)

Here, ${\mathfrak{e}}_m{}^{\mu }$ is a vierbein. The covariant derivative for fermions now reads

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mu }\psi ={\nabla }_{\mu }\psi -i{\mathit{g}}_{\mathrm{s}}{A}_{\mu }\psi ,\end{array}$$

(136)

with ${\nabla }_{\mu }$ begin the part accounting for the spacetime curvature via the spin connection. The field tensor $F_{\mu \nu }$ is defined via covariant derivatives. However, it does not involve Christoffel symbols due to its structure:

$$\begin{array}{c}\hfill F_{\mu \nu }={\nabla }_{\mu }{A}_{\nu }-{\nabla }_{\nu }{A}_{\mu }-i{\mathit{g}}_{\mathrm{s}}[A_{\mu },A_{\nu }]={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu }-i{\mathit{g}}_{\mathrm{s}}[A_{\mu },A_{\nu }].\end{array}$$

(137)

This results in the following $SU\left(N\right)$ Yang–Mills action in curved spacetime:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{A}=& \int d^{4}x\sqrt{-g}[\overline{\psi }\left(i{\mathfrak{e}}_m{}^{\mu }{\gamma }^m{\nabla }_{\mu }-m\right)\psi -\frac{1}{4}{\left(f_{\mu \nu }^a\right)}^2\hfill \\ \hfill & +{\mathit{g}}_{\mathrm{s}}\overline{\psi }\left({\mathfrak{e}}_m{}^{\mu }{\gamma }^m\right)\psi A_{\mu }-g^{\mu \nu }{g}^{\alpha \beta }{\mathit{g}}_{\mathrm{s}}{f}^{abc}{\partial }_{\mu }{A}_{\alpha }^a{A}_{\nu }^b{A}_{\beta }^c-\frac{1}{4}{\mathit{g}}_{\mathrm{s}}^2{f}^{amn}{f}^{aij}{g}^{\mu \nu }{g}^{\alpha \beta }{A}_{\mu }^m{A}_{\nu }^i{A}_{\alpha }^n{A}_{\beta }^j].\hfill \end{array}\end{array}$$

(138)

The perturbative quantisation of kinetic parts of the action is discussed above. The following theorem gives the perturbative expansion for the coupling of fermions to a gauge vector.

Theorem 30.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}{\mathit{g}}_{\mathrm{s}}\overline{\psi }\left({\mathfrak{e}}_m{}^{\mu }{\gamma }^m\right)\psi A_{\mu }\hfill \\ \hfill & =\sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\frac{d^{4}k}{{\left(2\pi \right)}^4}\prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+k+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times {\kappa }^n\overline{\psi }\left(p_2\right)\left[{\mathit{g}}_{\mathrm{s}}{\gamma }^m{T}^a{\left(\sqrt{-g}{\mathfrak{e}}_m{}^{\mu }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\right]\psi \left(p_1\right)A_{\mu }^a\left(k\right).\hfill \end{array}\end{array}$$

(139)

This expression produces the following Feynman rule:

(140)

The perturbative expansion for the cubic term in gauge vectors takes a similar form.

Theorem 31.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}(-{\mathit{g}}_{\mathrm{s}})f^{abc}{g}^{\mu \nu }{g}^{\alpha \beta }{\partial }_{\mu }{A}_{\alpha }^a{A}_{\nu }^b{A}_{\beta }^c\hfill \\ \hfill & =\sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\frac{d^4{p}_3}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+p_3+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times {\kappa }^n\left[(-i{\mathit{g}}_{\mathrm{s}})f^{abc}{\left(p_1\right)}_{\sigma }{\left(\sqrt{-g}{g}^{{\mu }_1{\mu }_3}{g}^{{\mu }_2\sigma }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}\right]A_{{\mu }_1}^a\left(p_1\right)A_{{\mu }_2}^b\left(p_2\right)A_{{\mu }_3}^c\left(p_3\right).\hfill \end{array}\end{array}$$

(141)

This expression produces the following rule:

(142)

Lastly, the four-vector coupling term has the following perturbative structure.

Theorem 32.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}\left(-\frac{1}{4}{\mathit{g}}_{\mathrm{s}}^2\right)f^{amn}{f}^{aij}{g}^{\mu \nu }{g}^{\alpha \beta }{A}_{\mu }^m{A}_{\nu }^i{A}_{\alpha }^n{A}_{\beta }^j\hfill \\ \hfill & =\sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\frac{d^4{p}_3}{{\left(2\pi \right)}^4}\frac{d^4{p}_4}{{\left(2\pi \right)}^4}\prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+p_3+p_4+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times \left(-\frac{1}{4}\right){\mathit{g}}_{\mathrm{s}}^2{\kappa }^n{f}^{amn}{f}^{aij}{\left(\sqrt{-g}{g}^{{\mu }_1{\mu }_3}{g}^{{\mu }_2{\mu }_4}\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{A}_{{\mu }_1}^m\left(p_1\right)A_{{\mu }_2}^n\left(p_2\right)A_{{\mu }_3}^i\left(p_3\right)A_{{\mu }_4}^j\left(p_4\right).\hfill \end{array}\end{array}$$

(143)

This results in the following interaction rule:

(144)

Finally, we are going to discuss gauge fixing. The Yang–Mills action respects the following gauge transformations:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \delta \psi =& i{\theta }^a{T}^a\psi ,\hfill \\ \hfill \delta A_{\mu }=& i{\theta }^a[T^a,A_{\mu }]+\frac{1}{{\mathit{g}}_{\mathrm{s}}}{\partial }_{\mu }{\theta }^a{T}^a,\hfill \\ \hfill \delta A_{\mu }^a=& \frac{1}{{\mathit{g}}_{\mathrm{s}}}\left[{\partial }_{\mu }{\theta }^a-g{f}^{abc}{\theta }^b{A}_{\mu }^c\right].\hfill \end{array}\end{array}$$

(145)

Here, ${\theta }^a$ are the gauge parameters. The standard Lorentz gauge fixing conditions simplify calculations in a flat spacetime.

$$\begin{array}{c}\hfill {\partial }^{\mu }{A}_{\mu }^a=0.\end{array}$$

(146)

For the case of curved spacetime, the covariant derivative replaces the standard derivative, which produces the new form of the Lorentz gauge fixing conditions:

$$\begin{array}{c}\hfill g^{\mu \nu }{\nabla }_{\mu }{A}_{\nu }^a=0.\end{array}$$

(147)

Introducing this gauge fixing term will bring the kinetic part of the vector field to the same form discussed in the previous section.

The ghost action is defined by the Faddeev–Popov determinant obtained from the gauge fixing condition:

$$\begin{array}{c}\hfill det\left[\frac{\delta }{\delta {\theta }^b}\left\{g^{\mu \nu }{\nabla }_{\mu }{A}_{\nu }^a\right\}\right]=det\left[\frac{1}{{\mathit{g}}_{\mathrm{s}}}{g}^{\mu \nu }{\nabla }_{\mu }\left({\delta }^{ab}{\nabla }_{\nu }-{\mathit{g}}_{\mathrm{s}}{f}^{abc}{A}_{\nu }^c\right)\right].\end{array}$$

(148)

It results in the following action:

$$\begin{array}{c}\hfill {\mathcal{A}}_{\mathrm{FP}}=\int d^{4}x\left[-g^{\mu \nu }{\nabla }_{\mu }{\overline{c}}^a{\nabla }_{\nu }{c}^a+{\mathit{g}}_{\mathrm{s}}{g}^{\mu \nu }{\nabla }_{\mu }{\overline{c}}^a{f}^{abc}{c}^b{A}_{\nu }^c\right].\end{array}$$

(149)

The action’s kinetic part is similar to a massless vector field. The section describing the interaction between ghosts, vectors, and gravitons allows for a perturbative expansion.

Theorem 33.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \int d^{4}x\sqrt{-g}\left[{\mathit{g}}_{\mathrm{s}}{\partial }_{\mu }{\overline{c}}^a{f}^{abc}{c}^b{A}^{c\mu }\right]\hfill \\ \hfill & =\sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\frac{d^{4}k}{{\left(2\pi \right)}^4}\prod _{i=0}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+k+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times i{\kappa }^n{\mathit{g}}_{\mathrm{s}}{\left(p_1\right)}_{\nu }{f}^{abc}{\overline{c}}^a\left(p_1\right)c^b\left(p_2\right){\left(\sqrt{-g}{g}^{\mu \nu }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{A}_{\mu }^c\left(k\right).\hfill \end{array}\end{array}$$

(150)

This expression produced the following rule:

(151)

4.7. General Relativity

General relativity is invariant with respect to local transformations spawned by coordinate transformations, so a gauge fixing procedure shall be performed. The perturbative approach describes gravity as small metric perturbations over a flat background. Therefore, it may seem that the theory is reduced to a gauge theory of rank-2 symmetric tensor, but this is not the case.

One shall distinguish a geometric theory from a theory of a rank-2 symmetric tensor. The difference dictates how to treat the gauge fixing condition. First, let us consider the case of a symmetric $h_{\mu \nu }$ tensor theory that admits the following gauge symmetry:

$$\begin{array}{c}\hfill \delta h_{\mu \nu }={\partial }_{\mu }{\zeta }_{\nu }+{\partial }_{\nu }{\zeta }_{\mu }.\end{array}$$

(152)

The $h_{\mu \nu }$ tensor is fundamental to this theory, and a gauge fixing condition can be expressed solely in terms of it. For example, one can use the following gauge fixing condition:

$$\begin{array}{c}\hfill {\partial }_{\mu }{h}^{\mu \nu }-\frac{1}{2}{\partial }^{\nu }h=0.\end{array}$$

(153)

This condition, along with others, determines the composition of the Faddeev–Popov ghosts. Because $h_{\mu \nu }$ is the fundamental object of the theory, the structure of divergences can only be expressed in terms of $h_{\mu \nu }$ alone. All geometric quantities, such as the Riemann tensor, Ricci tensor, scalar curvature, and others, are expressed via small metric perturbations, but the opposite is not true. Operators given in terms of $h_{\mu \nu }$ alone may not represent any geometric quantities. Therefore, in a gauge theory of a rank-2 symmetric tensor, one can expect to find divergences that cannot be described by geometric quantities, which makes it a non-geometric theory.

On the contrary, gauge transformations are generated by coordinate frame transformations within the geometric approach. This has two immediate implications. Firstly, within a geometrical theory, gauge transformations are given by the so-called Lie derivatives:

$$\begin{array}{c}\hfill \delta g_{\mu \nu }\stackrel{def}{=}{\mathcal{L}}_{\zeta }{g}_{\mu \nu }={\nabla }_{\mu }{\zeta }_{\nu }+{\nabla }_{\nu }{\zeta }_{\mu }.\end{array}$$

(154)

Here, ${\mathcal{L}}_{\zeta }$ is the Lie derivative with respect to an arbitrary vector field $\zeta$, which contains the gauge transformation parameters. Secondly, the gauge fixing condition must be expressed in geometrical quantities. Thus, gauge fixing conditions (153) are inconsistent with the geometrical approach and cannot be implemented. Instead, we use the following gauge fixing conditions:

$$\begin{array}{c}\hfill {\mathcal{G}}^{\mu }\stackrel{\mathrm{note}}{=}{g}^{\alpha \beta }{\mathrm{\Gamma }}_{\alpha \beta }^{\mu }=0.\end{array}$$

(155)

Combined with the perturbative expansion given by Equation (1), the gauge fixing conditions represented by Equation (155) produce an infinite series.

$$\begin{array}{c}\hfill {\mathcal{G}}^{\nu }=\frac{\kappa }{2}{g}^{\mu \nu }{g}^{\alpha \beta }\left[{\partial }_{\alpha }{h}_{\beta \mu }+{\partial }_{\beta }{h}_{\alpha \mu }-{\partial }_{\mu }{h}_{\alpha \beta }\right]=\kappa \left[{\partial }_{\mu }{h}^{\mu \nu }-\frac{1}{2}{\partial }^{\nu }h\right]+\mathcal{O}\left({\kappa }^2\right),\end{array}$$

(156)

The infinite expansion defines the ghost sector, and the leading term reproduces naive gauge fixing (153).

The difference between geometrical theories of gravity and a gauge theory of $h_{\mu \nu }$ tensor is marked by the need to use gauge fixing condition (155), and not conditions (153). Proceeding with The Faddeev–Popov prescription, firstly, we shall note that the gauge fixing condition ${\mathcal{G}}^{\mu }$ defined by (155) is a vector with mass dimension $+1$. Thus, one shall introduce an additional dimensional parameter in the gauge fixing term:

$$\begin{array}{c}\hfill {\mathcal{A}}_{\mathrm{H}+\mathrm{gf}}=\int d^{4}x\sqrt{-g}\left[-\frac{2}{{\kappa }^2}R+\frac{ϵ}{2{\kappa }^2}{g}_{\mu \nu }{\mathcal{G}}^{\mu }{\mathcal{G}}^{\nu }\right].\end{array}$$

(157)

Secondly, the Faddeev–Popov ghosts are also vectors. The structure of their action is defined by the variation in the gauge fixing term given by Equation (155):

$$\begin{array}{c}\hfill \delta G^{\mu }={\mathcal{L}}_{\zeta }\left[g^{\alpha \beta }{\mathrm{\Gamma }}_{\alpha \beta }^{\mu }\right]=\square {\zeta }^{\mu }-2{\mathrm{\Gamma }}_{\alpha \beta }^{\mu }{\nabla }^{\alpha }{\zeta }^{\beta }+R^{\mu }{}_{\nu }{\zeta }^{\nu }\end{array}$$

(158)

with $R_{\mu \nu }$ begin the Ricci tensor. Consequently, the ghost action reads:

$$\begin{array}{c}\hfill {\mathcal{A}}_{\mathrm{ghost}}=\int d^{4}x\sqrt{-g}\left[-g^{\alpha \beta }{g}^{\mu \nu }{\nabla }_{\alpha }{\overline{c}}_{\mu }{\nabla }_{\beta }{c}_{\nu }-2{\mathrm{\Gamma }}_{\alpha \beta }^{\mu }{\overline{c}}_{\mu }{\nabla }^{\alpha }{c}^{\beta }+R_{\mu \nu }{\overline{c}}^{\mu }{c}^{\nu }\right].\end{array}$$

(159)

In all other respects, the handling of the Faddeev–Popov ghosts is unchanged.

The standard perturbative expansion generates the interaction rules. The structure of graviton interactions is given by action (157):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{\mathrm{H}+\mathrm{gf}}& =\int d^{4}x\sqrt{-g}\left[-\frac{2}{{\kappa }^2}R+\frac{ϵ}{2{\kappa }^2}{g}_{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }{\mathrm{\Gamma }}_{\alpha \beta }^{\mu }{\mathrm{\Gamma }}_{\rho \sigma }^{\nu }\right]=\hfill \\ \hfill & =\int d^{4}x\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }\left(-\frac{2}{{\kappa }^2}\right)\left[{\mathrm{\Gamma }}_{\alpha \mu \rho }{\mathrm{\Gamma }}_{\sigma \nu \beta }-{\mathrm{\Gamma }}_{\alpha \mu \nu }{\mathrm{\Gamma }}_{\rho \beta \sigma }-\frac{ϵ}{4}{\mathrm{\Gamma }}_{\mu \alpha \beta }{\mathrm{\Gamma }}_{\nu \rho \sigma }\right].\hfill \end{array}\end{array}$$

(160)

It admits the following perturbative expansion:

Theorem 34.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{\mathrm{H}+\mathrm{gf}}=& \sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right)\hfill \\ \hfill & \times \left(2{\kappa }^n\right){\left(\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }{g}^{\rho \sigma }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{\left(p_1\right)}_{{\lambda }_1}{\left(p_2\right)}_{{\lambda }_2}{h}_{{\mu }_1{\nu }_1}\left(p_1\right)h_{{\mu }_2{\nu }_2}\left(p_2\right)\hfill \\ \hfill & \times \left[{\left({\mathrm{\Gamma }}_{\alpha \mu \rho }\right)}^{{\lambda }_1{\mu }_1{\nu }_1}{\left({\mathrm{\Gamma }}_{\sigma \nu \beta }\right)}^{{\lambda }_2{\mu }_2{\nu }_2}-{\left({\mathrm{\Gamma }}_{\alpha \mu \nu }\right)}^{{\lambda }_1{\mu }_1{\nu }_1}{\left({\mathrm{\Gamma }}_{\rho \beta \sigma }\right)}^{{\lambda }_2{\mu }_2{\nu }_2}-\frac{ϵ}{4}{\left({\mathrm{\Gamma }}_{\mu \alpha \beta }\right)}^{{\lambda }_1{\mu }_1{\nu }_1}{\left({\mathrm{\Gamma }}_{\nu \rho \sigma }\right)}^{{\lambda }_2{\mu }_2{\nu }_2}\right].\hfill \end{array}\end{array}$$

(161)

The following formula gives the complete expression for the graviton vertex:

(162)

The sum goes over all possible permutations of graviton parameters $\left\{{\mu }_i{\nu }_i{p}_i\right\}$.

The ghost action is treated similarly.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{\mathrm{ghost}}=& \int d^{4}x\sqrt{-g}\left[-g^{\alpha \beta }{g}^{\mu \nu }{\nabla }_{\alpha }{\overline{c}}_{\mu }{\nabla }_{\beta }{c}_{\nu }-2{\mathrm{\Gamma }}_{\alpha \beta }^{\mu }{\overline{c}}_{\mu }{\nabla }^{\alpha }{c}^{\beta }+R_{\mu \nu }{\overline{c}}^{\mu }{c}^{\nu }\right]\hfill \\ \hfill =& \int d^{4}x\sqrt{-g}\left[-g^{\mu \nu }{g}^{\alpha \beta }{\partial }_{\alpha }{\overline{c}}_{\mu }{\partial }_{\beta }{c}_{\nu }\right]\hfill \\ \hfill & +\int d^{4}x\sqrt{-g}{g}^{\mu \alpha }{g}^{\nu \beta }{g}^{\rho \sigma }\left[{\mathrm{\Gamma }}_{\beta \rho \alpha }{\partial }_{\sigma }{\overline{c}}_{\mu }{c}_{\nu }-{\mathrm{\Gamma }}_{\alpha \rho \beta }{\overline{c}}_{\mu }{\partial }_{\sigma }{c}_{\nu }+{\partial }_{\rho }{\mathrm{\Gamma }}_{\sigma \alpha \beta }{\overline{c}}_{\mu }{c}_{\nu }-{\partial }_{\alpha }{\mathrm{\Gamma }}_{\rho \beta \sigma }{\overline{c}}_{\mu }{c}_{\nu }\right]\hfill \\ \hfill & +\int d^{4}x\sqrt{-g}{g}^{\mu \alpha }{g}^{\nu \beta }{g}^{\rho \sigma }{g}^{\lambda \tau }\left[{\mathrm{\Gamma }}_{\rho \alpha \lambda }{\mathrm{\Gamma }}_{\sigma \beta \tau }-{\mathrm{\Gamma }}_{\rho \alpha \beta }{\mathrm{\Gamma }}_{\sigma \lambda \tau }+{\mathrm{\Gamma }}_{\alpha \rho \lambda }{\mathrm{\Gamma }}_{\beta \sigma \tau }\right]{\overline{c}}_{\mu }{c}_{\nu }.\hfill \end{array}\end{array}$$

(163)

It has the following perturbative expansion:

Theorem 35.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{A}}_{\mathrm{ghost}}=& \sum _{n=0}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right){\overline{c}}_{\mu }\left(p_1\right)c_{\nu }\left(p_2\right)\hfill \\ \hfill & \times {\kappa }^n{\left(\sqrt{-g}{g}^{\mu \nu }{g}^{\alpha \beta }\right)}^{{\rho }_1{\sigma }_1\cdots {\rho }_n{\sigma }_n}{\left(p_1\right)}_{\alpha }{\left(p_2\right)}_{\beta }\hfill \\ \hfill +& \sum _{n=1}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right){\overline{c}}_{\mu }\left(p_1\right)c_{\nu }\left(p_2\right)\hfill \\ \hfill & \times {\kappa }^n(-1){\left(\sqrt{-g}{g}^{\mu \alpha }{g}^{\nu \beta }{g}^{\rho \sigma }\right)}^{{\rho }_2{\sigma }_2\cdots {\rho }_n{\sigma }_n}{\left(k_1\right)}_{\lambda }\left[{\left(p_1\right)}_{\sigma }{\left({\mathrm{\Gamma }}_{\beta \rho \alpha }\right)}^{\lambda {\rho }_1{\sigma }_1}-{\left(p_2\right)}_{\sigma }{\left({\mathrm{\Gamma }}_{\alpha \rho \beta }\right)}^{\lambda {\rho }_1{\sigma }_1}\right.\hfill \\ \hfill & \left.+{\left(k_1\right)}_{\rho }{\left({\mathrm{\Gamma }}_{\sigma \alpha \beta }\right)}^{\lambda {\rho }_1{\sigma }_1}-{\left(k_1\right)}_{\alpha }{\left({\mathrm{\Gamma }}_{\rho \beta \sigma }\right)}^{\lambda {\rho }_1{\sigma }_1}\right]\hfill \\ \hfill +& \sum _{n=2}^{\infty }\int \frac{d^4{p}_1}{{\left(2\pi \right)}^4}\frac{d^4{p}_2}{{\left(2\pi \right)}^4}\prod _{i=1}^n\frac{d^4{k}_i}{{\left(2\pi \right)}^4}{\left(2\pi \right)}^4\delta \left(p_1+p_2+\sum k_i\right)h_{{\rho }_1{\sigma }_1}\left(k_1\right)\cdots h_{{\rho }_n{\sigma }_n}\left(k_n\right){\overline{c}}_{\mu }\left(p_1\right)c_{\nu }\left(p_2\right)\hfill \\ \hfill & \times {\kappa }^n(-1){\left(\sqrt{-g}{g}^{\mu \alpha }{g}^{\nu \beta }{g}^{\rho \sigma }{g}^{\lambda \tau }\right)}^{{\rho }_3{\sigma }_3\cdots {\rho }_n{\sigma }_n}{\left(k_1\right)}_{{\lambda }_1}{\left(k_2\right)}_{{\lambda }_2}\left[{\left({\mathrm{\Gamma }}_{\rho \alpha \lambda }\right)}^{{\lambda }_1{\rho }_1{\sigma }_1}{\left({\mathrm{\Gamma }}_{\sigma \beta \tau }\right)}^{{\lambda }_2{\rho }_2{\sigma }_2}\right.\hfill \\ \hfill & \left.-{\left({\mathrm{\Gamma }}_{\rho \alpha \beta }\right)}^{{\lambda }_1{\rho }_1{\sigma }_1}{\left({\mathrm{\Gamma }}_{\sigma \lambda \tau }\right)}^{{\lambda }_2{\rho }_2{\sigma }_2}+{\left({\mathrm{\Gamma }}_{\alpha \rho \lambda }\right)}^{{\lambda }_1{\rho }_1{\sigma }_1}{\left({\mathrm{\Gamma }}_{\beta \sigma \tau }\right)}^{{\lambda }_2{\rho }_2{\sigma }_2}\right].\hfill \end{array}\end{array}$$

(164)

The complete expression describing graviton-ghost vertices reads:

(165)

The standard procedure derives propagators for ghosts and gravitons. The following expression gives the ghost propagator.

(166)

The graviton propagator contains the gauge fixing parameter $ϵ$. The propagator corresponds to the part of the microscopic action quadratic in perturbations:

$$\begin{array}{c}\hfill \int d^{4}x\sqrt{-g}\left[-\frac{2}{{\kappa }^2}R+\frac{ϵ}{2{\kappa }^2}{\mathcal{G}}_{\mu }{\mathcal{G}}^{\mu }\right]=\int d^{4}x\left[-\frac{1}{2}{h}^{\mu \nu }{\mathcal{D}}_{\mu \nu \alpha \beta }\left(ϵ\right)\square h^{\alpha \beta }\right]+\mathcal{O}\left({\kappa }^1\right).\end{array}$$

(167)

The Nieuwenhuizen operators provide a good mean to express the $\mathcal{D}$ operator in the momentum representation [38,88]

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mu \nu \alpha \beta }\left(ϵ\right)=\frac{3ϵ-8}{4}{P}_{\mu \nu \alpha \beta }^0+\frac{ϵ}{2}{P}_{\mu \nu \alpha \beta }^1+P_{\mu \nu \alpha \beta }^2-\frac{ϵ}{4}{\overline{P}}_{\mu \nu \alpha \beta }^0-\frac{ϵ}{4}{\overline{\overline{P}}}_{\mu \nu \alpha \beta }^0.\end{array}$$

(168)

Operators $P^0$ and $P^2$ are gauge invariant, making the operator non-invertible unless $ϵ\ne 0$. The inverse operator reads:

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mu \nu \alpha \beta }^{-1}\left(ϵ\right)=-\frac{1}{2}{P}_{\mu \nu \alpha \beta }^0+\frac{2}{ϵ}{P}_{\mu \nu \alpha \beta }^1+P_{\mu \nu \alpha \beta }^2-\frac{3ϵ-8}{2ϵ}{\overline{P}}_{\mu \nu \alpha \beta }^0-\frac{1}{2}{\overline{\overline{P}}}_{\mu \nu \alpha \beta }^0.\end{array}$$

(169)

The given formula expresses the graviton propagator in any arbitrary gauge.

(170)

We can discuss the general case, but using $ϵ=2$ in practical applications is easier. With this value, the operator ${\mathcal{D}}^{-1}$ becomes much simpler in form:

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mu \nu \alpha \beta }^{-1}\left(2\right)=\frac{1}{2}\left[{\eta }_{\mu \alpha }{\eta }_{\nu \beta }+{\eta }_{\mu \beta }{\eta }_{\nu \alpha }-{\eta }_{\mu \nu }{\eta }_{\alpha \beta }\right].\end{array}$$

(171)

5. FeynGrav

FeynGrav is a package for Wolfram Mathematica extending FeynCalc functionality [27,28]. FeynCalc provides tools to study both tree and loop-level amplitudes. At the same time, there are many packages further extending its functionality [87,89,90], which makes it a platform for different computational tools of high energy physics. Because of these reasons, FeynCalc was chosen to implement the perturbative quantum gravity framework described above.

This subsection is split into two parts to discuss two different aspects of FeynGrav. The first subsection discussed the general features of the package and its published versions. The second subsection discussed the commands implemented in the latest available version.

5.1. Implementation

We shall begin with a discussion of the existing versions of FeynGrav. There are a few published versions of the package, each extending its functionality and improving performance. The package is constantly developing, and the latest version is publicly available at [91]. We shall briefly discuss the features of all these versions for completeness.

The original version 1.0 of FeynGrav was published in [24]. This version only considered matter with massless spin 0, 1, 2, and $1/2$, and did not include the $SU\left(N\right)$ Yang–Mills model. The recursive relations discussed in this paper were yet to be discovered during the development of this version. Consequently, the package used a less effective algorithm to generate all perturbative expansions. The package’s applicability was also limited, since the gauge-fixing algorithm was not considered in detail.

The second version published in [25] presented a significant update. This version added interaction rules for massive matter with spins 0, 1, $1/2$, and implemented rules for the $SU\left(N\right)$ Yang–Mills. The gauge-fixing procedure was fully addressed, so the implemented interaction rules became applicable at any perturbation theory level. However, the recursive algorithm still needed to be discovered, so the package used the same algorithm for perturbative expansions.

The recursive relations were discovered after the publication [25], and were implemented in the latest FeynGrav version 2.1. The discovery of the recursive algorithm allowed expressions for the interaction rules to be generated more efficiently, so the corresponding libraries were updated. This version also implemented the polarisation tensors for gravity, discussed below. Lastly, minor misprints were corrected, and the sample file was significantly updated and improved.

The interaction rules for the Horndeski theory still need to be implemented in any existing version of FeynGrav. Because of their shared length and complexity, their implementation takes time and is expected to be published soon.

Let us discuss the package structure that remains similar in all the versions. The package addressed a few different challenges, consisting of a few semi-independent modules. The main computational challenge is the generation of the perturbative expansion terms. The package’s core addresses this problem, but operates separately from the main file. In turn, the main file of the package only imports the interaction rules in the FeynCalc environment, allowing a user to operate with them.

The above recursive relations are implemented in a series of sub-packages in a separate folder. Each package provides a tool to calculate a separate family of tensors, and can be used independently within FeynCalc. As discussed above, it is essential to introduce tensors with particular symmetry for the Lorentz indices. With this symmetry imposed, a typical tensor with $2n$ pair of Lorentz indices will have approximately $2^{n}n!$ terms. Because of this large number of terms, a computation of an interaction rule involving many particles can take significant time. To soften this issue, the interaction rules are calculated separately.

The interaction rules discussed above are implemented in a single separate sub-package. Because such rules only require information about perturbative expansions, this sub-package depends on sub-packages describing perturbative expansions. It shall be run independently from FeynGrav to generate libraries containing the final expressions for the interaction rules.

The main package file is directly independent of these sub-packages, since it only imports the existing libraries and places them in the FeynCalc environment. Because of this package structure, the final user does not have to perform computationally heavy calculations of the interaction rules, significantly improving the package performance. When FeynGrav imports the rules, FeynCalc performs index contractions and other operations that constrain FeynGrav’s performance.

5.2. Interaction Rules

The scalar field kinetic and potential energy interaction rules are implemented with commands “GravitonScalarVertex” and “GravitonScalarPotentialVertex”. The first command takes four arguments: the array of graviton indices, two momenta of scalar fields, and the scalar field mass. The second command takes two arguments: the array of graviton indices and the scalar field coupling coupling. FeynGrav also contains a command realising the scalar field propagator. The command is “ScalarPropagator”, and takes two arguments: the scalar field momentum and mass.

Table 1. Examples of interaction rules for the scalar, Proca, and massless vector fields.

Diagram	Command
	ScalarPropagator$[p,m]$
	GravitonScalarVertex$[\{{\rho }_1,{\sigma }_1\},p_1,p_2,m]$
	GravitonScalarVertex$[\{{\rho }_1,{\sigma }_1,{\rho }_2,{\sigma }_2\},p_1,p_2,m]$
	GravitonScalarPotentialVertex$[\{{\rho }_1,{\sigma }_1\},{\lambda }_3]$
	GravitonScalarPotentialVertex$[\{{\rho }_1,{\sigma }_1,{\rho }_2,{\sigma }_2\},{\lambda }_4]$
	ProcaPropagator$[p,m]$
	GravitonMassiveVectorVertex$[\{{\rho }_1,{\sigma }_1\},{\lambda }_1,p_1,{\lambda }_2,p_2,m]$
	GravitonVectorVertex$[\{{\rho }_1,{\sigma }_1,k_1\},{\lambda }_1,p_1,{\lambda }_2,p_2]$
	GravitonVectorGhostVertex$[\{{\rho }_1,{\sigma }_1,k_1\},p_1,p_2]$

presents examples of these commands’ usage.

Interaction rules for the Proca field are implemented with a single command “GravitonMassiveVectorVertex”. The name is chosen for the sake of naming consistency. The command takes six arguments: the array of graviton indices, the Lorentz index and momentum of the first vector particle, the Lorentz index and momentum of the second vector particle, and the Proca field mass. The package also has the Proca propagator implementation. The command “ProcaPropagator” implements the propagator, and takes four arguments: two Lorentz indices, and the momentum and mass of the Proca field. Table 1 presents examples of these commands’ usage.

The interaction rules for a single massless vector field are implemented with two commands, “GravitonVectorVertex” and “GravitonVectorGhostVertex”. The command “GravitonVectorVertex” takes five arguments: the array of Lorentz indices and momenta of gravitons, the Lorentz index and momentum of the first vector particle, the Lorentz index and momentum of the second particle. The command “GravitonVectorGhostVertex” takes three arguments: the array of Lorentz indices and momenta of gravitons, and the momenta of the Faddeev–Popov ghost. Table 1 presents examples of these commands. Let us also note that FeynCalc already provides the propagators required to treat a massless scalar field.

The interaction rules for Dirac fermions are implemented with a single command “GravitonFermionVertex”. The command takes four arguments: the array of graviton Lorentz indices, the momentum of the in-going fermion line, the momentum of the out-going fermion line, and the fermion mass. The fermion propagator has already been implemented in FeynCalc.

Table 2. Examples of interaction rules for the scalar, Proca, and massless vector fields.

Diagram	Command
	GravitonFermionVertex$[\{{\rho }_1,{\sigma }_1\},p_1,p_2,m]$
	GravitonGluonVertex$[\{{\rho }_1,{\sigma }_1,k_1\},p_1,{\lambda }_1,a_1,p_2,{\lambda }_2,a_2]$
	GravitonGluonVertex$[\{{\rho }_1,{\sigma }_1,k_1\},p_1,{\lambda }_1,a_1,p_2,{\lambda }_2,a_2,{\lambda }_3,p_3,a_3]$
	GravitonGluonVertex$[\{{\rho }_1,{\sigma }_1,k_1\},p_1,{\lambda }_1,a_1,p_2,{\lambda }_2,a_2,{\lambda }_3,p_3,a_3,{\lambda }_4,p_4,a_4]$
	GravitonQuarkGluonVertex$[\{{\rho }_1,{\sigma }_1\},\lambda ,a]$
	GravitonYMGhostVertex$[\{{\rho }_1,{\sigma }_1\},p_1,a_1,p_2,a_2]$
	GravitonGluonGhostVertex$[\{{\rho }_1,{\sigma }_1\},p_1,{\lambda }_1,a_1,p_2,{\lambda }_2,a_2,p_3,{\lambda }_3,a_3]$

provides an example of this interaction rule.

The $SU\left(N\right)$ Yang–Mills model implementation is sophisticated and involves several commands. The command “GravitonGluonVertex implements a coupling of a few gluons to gravity. The command’s first argument is an array of Lorentz indices and momenta. The other arguments describe gluons Lorentz indices, momenta, and $SU\left(N\right)$ indices. The command can describe a coupling of two, three, and four gluons to gravity. Gravitational coupling of quarks kinetic energy matches the expression for the coupling of a Dirac fermion. The coupling of the quark–gluon interaction energy is described by the “GravitonQuarkGluonVertex” command. It takes only three arguments: the array of graviton Lorentz indices, the quark–gluon vertex Lorentz index, and the colour index. Lastly, two commands are responsible for gravitational coupling to the Faddeev–Popov ghosts of the Yang–Mills theory. The command “GravitonYMGhostVertex” describes the coupling to the ghost itself. The command takes five arguments: the array of graviton Lorentz indices, the momentum and colour index of the first ghost, and the momentum and colour index of the second ghost. The command “GravitonGluonGhostVertex” corresponds to the gravitational coupling of a vertex describing the interaction between two ghosts and one gluon. The command’s arguments are the array of graviton Lorentz indices, the Lorentz indices, momenta, and colour indices of other particles. Table 2 lists examples of these commands.

Lastly, the following commands describe the gravitational sector. The command “GravitonPropagator” implements the graviton propagator with “FeynAmpDenominator” functions from FeynCalc. This command shall be used for loop calculations. The command “GravitonPropagatorAlternative” implements the graviton propagator with the simple scalar products in the denominators. This command shall only be used in three-level calculations. Two more commands, “GravitonPropagatorTop” and “GravitonPropagatorTopFAD”, generate the numerator of the graviton propagator alone. The first one uses simple scalar products, while the second one uses “FeynAmpDenominator”. All these commands take five arguments: four Lorentz indices and the momenta of a graviton. The command “GravitonVertex” corresponds to the n-graviton vertex and takes $3n$ arguments: two Lorentz indices and the momentum of each graviton. The Faddeev–Popov ghosts are vectors, so their propagator receives an additional Minkowski metric in the numerator. One can use the command for the Yang–Mills ghost propagator and manually multiply it on the flat metric. Because of this, no new command corresponding to this propagator is present in FeynGrav. Lastly, the rules for graviton coupling to the corresponding Faddeev–Popov ghosts are implemented with “GravitonGhostVertex”. It takes five arguments: the array of Lorentz indices and momenta of the gravitons, and Lorenz indices and momenta of ghosts.

Table 3. Examples of interaction rules for gravitons.

Diagram	Command
	GravitonPropagator$[\mu ,\nu ,\alpha ,\beta ,p]$
	GravitonVertex$[{\rho }_1,{\sigma }_1,p_1,{\rho }_2,{\sigma }_2,p_2,{\rho }_3,{\sigma }_3,p_3]$
	GravitonVertex$[{\rho }_1,{\sigma }_1,p_1,{\rho }_2,{\sigma }_2,p_2,{\rho }_3,{\sigma }_3,p_3,{\rho }_4,{\sigma }_4,p_4]$
	GravitonGhostVertex$[\{{\rho }_1,{\sigma }_1,k_1\},{\lambda }_1,p_1,{\lambda }_2,p_2]$

presents these rules.

FeynGrav also implements a few additional tools. Firstly, the package provides a realisation of the standard gauge projectors:

$$\begin{array}{cccc}\hfill {\theta }_{\mu \nu }\left(p\right)& ={\eta }_{\mu \nu }-\frac{p_{\mu }{p}_{\nu }}{p^2},\hfill & \hfill {\overline{\theta }}_{\mu \nu }& =\frac{p_{\mu }{p}_{\nu }}{p^2}.\hfill \end{array}$$

(172)

They are realised with commands “GaugeProjector”, “GaugeProjectorBar”, “GaugeProjectorFAD”, and “GaugeProjectorBarFAD”. The first two commands realise these projectors with the simple scalar products, while the former two commands use “FeynAmpDenominator”.

Secondly, the package provides a realisation of the Nieuwenhuizen operators [38,88,92]. These operators $P_{\mu \nu \alpha \beta }^0\left(p\right)$, $P_{\mu \nu \alpha \beta }^1\left(p\right)$, $P_{\mu \nu \alpha \beta }^2\left(p\right)$, ${\overline{P^0}}_{\mu \nu \alpha \beta }\left(p\right)$, ${\overline{\overline{P^0}}}_{\mu \nu \alpha \beta }\left(p\right)$ generalise the standard projectors for the spin-2 systems. Their features were discussed in detail in a previous publication [24]. Commands “NieuwenhuizenOperator0”, “NieuwenhuizenOperator1”, “NieuwenhuizenOperator2”, “NieuwenhuizenOperator0Bar”, and “NieuwenhuizenOperator0BarBar” implement these operators using only the simple scalar product. The set of command “NieuwenhuizenOperator0FAD”, ⋯, “NieuwenhuizenOperator0BarBarFAD” implements them with “FeynAmpDenominator”.

Lastly, the package implements the polarisation tensors for gravitons. Following the standard procedure, the polarisation tensor ${ϵ}_{\mu \nu }^{\pm }$ for gravitons is defined through the polarisation vector ${ϵ}_{\mu }^{\pm }$ for the standard vector field:

$$\begin{array}{c}\hfill {ϵ}_{\mu \nu }^{\pm }={ϵ}_{\mu }^{\pm }{ϵ}_{\nu }^{\pm }.\end{array}$$

(173)

The command “PolarizationTensor” implements this definition. The command takes three arguments: two Lorentz indices and one four-momentum. The command directly multiplies polarisation vectors implemented in FeynCalc. Due to the definitions used in FeynCalc, the implemented polarisation tensor is neither traceless nor transverse—the command “SetPolarizationTensor” the tensors traceless and transverse.

5.3. Examples

Several publications have used FeynGrav [24,25,32,33,35]. These publications discuss the problems in detail, so we will only provide a brief discussion. We also focus on two vivid examples.

The first illustrative example of FeynGrav usage is the calculation of $2\to 2$ graviton scattering cross sections at the leading order of perturbation theory. The on-shell cross section was calculated long ago without using any contemporary computational tools [10], which makes such calculations notoriously complicated. To be exact, in the original article [10], the author operated with expressions symmetries in a certain way to make the calculation more compact. Nonetheless, even the basic description of the calculations takes a few pages in the original publication.

FeynGrav extremely simplifies such calculations. Firstly, one shall fix all the momenta on the mass shell and define all relations for scalar products of polarisation tensors and momenta, since FeynCalc does not do this automatically. Secondly, one calculates the matrix element for the s-channel amplitude:

(174)

An average computer not designed for intensive computational tasks can complete this element in less than five minutes. One obtains the expressions for t and u-channels similarly. The following expression gives the contact four-graviton interaction:

(175)

This expression takes less than one minute to be computed.

Lastly, one constructs the complete scattering amplitude accounting for all contributions. It recovers the well-known expressions for chiral amplitudes in $d=4$:

$$\begin{array}{cccccc}\hfill \mathcal{M}(++++)& =i\frac{{\kappa }^2}{4}\frac{s^4}{stu},\hfill & \hfill \mathcal{M}(+-+-)& =i\frac{{\kappa }^2}{4}\frac{u^4}{stu},\hfill & \hfill \mathcal{M}(+--+)& =i\frac{{\kappa }^2}{4}\frac{t^4}{stu},\hfill \end{array}$$

(176)

$$\begin{array}{cc}\hfill \mathcal{M}(++--)& =\mathcal{M}(+++-)=0.\hfill \end{array}$$

This example vividly shows that such complicated yet essential calculations are significant with FeynGrav. Moreover, paper [24] presents an explicit expression for the scattering amplitude expressed in terms of graviton chiralities.

The second example is the calculation of the one-loop graviton-scalar vertex operator. Since a scalar field interacts with gravity, the interaction vertex receives correction at the loop level. The one-loop vertex operator describes these corrections:

(177)

In contrast to the previous example, this expression was only calculated with FeynGrav because it involved hundreds of terms. The expression itself is too long to be presented in printed form. It is published in an open-access repository [93]. The obtained expression allows one to study the low energy limit of a gravitational scattering of two scalars [33]:

$$\begin{array}{c}\hfill \frac{d\sigma }{d\mathrm{\Omega }}{|}_{\mathrm{low}}\mathrm{energy}=\frac{1}{4}{\left(G_{\mathrm{N}}\mu m_1{m}_2\right)}^2\frac{1}{p^4{sin}^4\frac{\theta }{2}}\left[1+32{\pi }^2{G}_{\mathrm{N}}p(m_1+m_2)sin\frac{\theta }{2}+\mathcal{O}\left(G_{\mathrm{N}}^2\right)\right].\end{array}$$

(178)

In this expression for the differential cross section, p is the centre of mass scattering momentum, $\theta$ is the scattering angle, $m_1$ and $m_2$ are masses of the scattered particles, and $\mu$ is the reduced mass.

6. Summary and Further Development

This paper reviews a recently developed theoretical framework for efficient computation of Feynman rules for perturbative quantum gravity. The perturbative approach to quantum gravity allows one to study quantum gravitational effects within the standard quantum field theory framework. The resulting theory is effective, so it cannot be indefinitely extended in the ultraviolet region and treated similarly to other renormalisable theories. Due to its effective nature, the theory admits infinite interaction terms, all parametrised with a single gravitational coupling. The discussed theoretical framework provides a universal way to compute such interaction rules at any order of perturbation theory.

Perturbative quantum gravity admits the factorisation theorem applicable to any gravity model within Riemann geometry. The theorem states that the action splits into two parts in a broad class of gravity models. One part does not involve derivatives, and is constituted of three independent factors that are infinite series. The other part involves derivatives, but it is a finite expression that can be calculated explicitly.

Certain recursive relations exist for the factors generating infinite series. They provide a way to efficiently calculate their contribution at any order of perturbation theory. However, the number of terms in each contribution grows faster than the factorial number of particles involved in the interaction. This feature makes manual calculations of such factors challenging.

FeynGrav provides a tool to operate with the theory at the practical level. It implements the discussed theoretical framework within FeynCalc, which possesses a broad spectrum of tools to operate with expressions within the quantum field theory.

The latest version implements the interaction rules for massless and massive scalars, for scalar field potential interaction, for massive and massless Dirac fermions, for Proca field, for massless vector field, for $SU\left(N\right)$ Yang–Mills model, and general relativity. The package also implements auxiliary tools such as graviton polarisation tensors, gauge projectors, and their generalisation for the spin-2 case.

The package is rapidly developing, and is expected to receive a new update soon. The update will be devoted to the Horndeski gravity, the most general scalar–tensor gravity admitting second-order field equations and the minimal coupling to the matter degrees of freedom.

Further development of the package aims in two directions: to improve the performance and to account for other relevant gravity models. The discovery of the recursive algorithm discussed in this paper significantly increased performance. At the same time, the same relations directly indicate the limits of this existing approach. The number of terms in an n particle interaction vertex grows faster than $n!$, which poses an additional computation challenge. Although such behaviour is expected for an effective theory considered perturbatively, we will strive to improve the computational algorithm even further to push the perturbation theory to its limits.

Consideration of more sophisticated yet relevant models of gravity is another challenge. The Horndeski gravity, the most studied extension of general relativity, is studied in this paper, and its implementation is to enter the next version of FeynGrav. There are few options for other models of gravity that would be relevant to the contemporary research field. The most straightforward generalisation would be massive gravity, but the contemporary observational data imposes strong constraints on the graviton mass. Another possible way to extend FeynGrav is to implement Beyond Horndeski theories, which introduce non-minimal couplings between the scalar field and matter. Lastly, it may be fruitful to consider some supersymmetric expressions and examine whether they can be treated similarly.

In conclusion, the theoretical framework discussed in the paper and the FeynGrav package provide a valuable and efficient tool to study perturbative quantum gravity, with further development aimed at improving performance and incorporating other relevant gravity models.