*mps*

# FeynGrav 3.0

Boris Latosh<EMAIL_ADDRESS>Particle Theory and Cosmology Group,
Center for Theoretical Physics of the Universe, Institute for Basic Science
(IBS), Daejeon, 34126, Korea

(CTPU-PTC-24-15)

###### Abstract

We present the new version of the FeynGrav. The package provides tools to
operate with Feynman rules for quantum gravity within FeynCalc. The latest
version improves package efficiency and implements new physical models. We
discover recurrent relations between metric factors that enhance computational
efficiency. We discuss gravitational interaction with Horndeski gravity,
quadratic gravity, and the simplest axion-like coupling. We implemented the
massive graviton propagator and discussed the possibility of implementing
massive gravity within the package.

##  1 Introduction

This paper discussed the new version of FeynGrav [1, 2, 3, 4], a Wolfram
Mathematica package based on FeynCalc [5, 6, 7, 8], providing a set of tools
to operate with Feynman rules within perturbative quantum gravity.

Perturbative quantum gravity is the simplest approach to quantum gravity based
on the standard quantum field theory. The view on the perturbative quantum
gravity is twofold. Firstly, one can treat it as a quantum field theory of
small metric perturbations $h_{\mu\nu}$ that exist on a flat background [2, 3,
4, 9, 10, 11, 12, 13, 14, 15]. In that case, one composes the complete
spacetime metric $g_{\mu\nu}$ from the background metric $\eta_{\mu\nu}$ and
the perturbations $h_{\mu\nu}$:

$\displaystyle
g_{\mu\nu}\overset{\text{def}}{=}\eta_{\mu\nu}+\kappa\,h_{\mu\nu}.$ (1)

Here $\kappa$ is a gravitational coupling with the mass dimension $-1$, which
ensures that the perturbations have canonical mass dimension. The coupling is
related to the Newton’s constant $G_{\text{N}}$:

$\displaystyle\kappa^{2}=32\,\pi\,G_{\text{N}}.$ (2)

Although the metric $g_{\mu\nu}$ is a finite expression, it gives rise to
infinite series. For instance, the inverse metric $g^{\mu\nu}$ is an infinite
series in $\kappa$. The volume factor $\sqrt{-g}$ is also an infinite series
due to the presence of the square root. Similarly, the Riemann tensor
$R_{\mu\nu\alpha\beta}$ and the Christoffel symbols $\Gamma^{\alpha}_{\mu\nu}$
involve the inverse metric, so they are infinite series as well. The resulting
theory admits infinite interaction terms, all suppressed by different powers
of the same gravitational coupling.

Secondly, one can view perturbative quantum gravity as a gauge theory of
massless particles with chirality $\pm 2$ described by the standard quantum
field theory. In turn, quantum field theory admits many powerful tools that
constrain the structure of the theory. For instance, gauge invariance
constrains the form of interaction rules, while the optical theorem recovers a
point of scattering amplitudes. These tools allow one to retrieve information
about graviton propagator, interaction vertices, and certain scattering
amplitudes [16, 17, 18, 19, 20, 21, 22, 23]. These approaches are equivalent
for two reasons. First, they recover the same results with different
techniques. Secondly, they use other means of description to operate with the
same physical setup. The flat background in the first approach ensures that
one can consistently define the asymptotic states respecting the Poincare
symmetry. The presence of the Poincare symmetry, in turn, allows one to define
states with definite mass and spin/chirality [24]. This setup allows one to
use the formalism of quantum field theory. In turn, the theory invariance with
respect to the change of coordinates ensures an infinite number of interaction
terms. Within the first approach, the coordinate invariance of action requires
the presence of the volume factor $\sqrt{-g}$ to preserve the coordinate
volume invariant of the inverse metric $g^{\mu\nu}$ to contract covariant
derivatives and other geometric quantities which are infinite perturbative
series. Within the second approach, the coordinate invariance is a gauge
symmetry. The gauge invariance of scattering matrix elements ensures infinite
graviton interaction terms exist.

It is essential to discuss the physical features of the perturbative quantum
gravity, but it is not the main scope of this article. We highlighted the main
features of the theory above and omitted further discussion. Papers [4, 10,
11, 13] discuss the theory’s applicability domain and predictive ability in
more detail.

The new version of FeynGrav implements the following features. First and
foremost, it implements the recently discovered recursive relations. It is
possible to obtain explicit formulas for many geometric quantities [2].
However, they typically involve multiple summations, which affects the
performance. The recursive relations simplify such summations and improve the
computational performance.

Secondly, the new version implements Feynman rules for a new series of models.
First and foremost, we address the Horndeski gravity. One can obtain the
interaction rules for all Horndeski models within the perturbative quantum
gravity. It was firstly discussed in [4], and we extend this discussion below.
On the practical ground, most parts of Horndeski’s models that involve the
coupling of many scalars are irrelevant. Consequently, the new FeynGrav
version implements the interaction rules only for a few essential cases.

Further, the new version implements the massive graviton propagator. Massive
gravity on its own is a separate branch of research reached with phenomenology
[25, 26]. It is desirable to obtain a similar algorithm deriving the
interaction rules for any order in perturbation theory, but the task faces a
significant challenge. Due to the non-vanishing graviton mass, the theory does
not admit the gauge symmetry, which results in particular challenges. We
discuss this issue in detail below and will address it in further
publications.

The new version implements the simplest case of the axion-like coupling
between a scalar field and a single vector field. The term describes a
contraction of a vector field tensor $F$, and the dual tensor $\tilde{F}$
takes a special place in particle physics. The term is a complete derivative,
so it does not contribute to the classical field equations. Within the $SU(N)$
Yang-Mills theory, it is known as the $\theta$ term. It is responsible for the
existence of instantons and the related phenomenology [27, 28, 29, 30, 31, 32,
33]. The need to explain the value of the $\theta$ term produced various axion
models that involve a coupling of a scalar field to the $F\tilde{F}$ term.
This publication does not address the coupling for the $SU(N)$ Yang-Mills case
since it deserves a separate publication. This paper considers the simplest
case when a scalar field is coupled to the $F\tilde{F}$ term. Despite its
simplicity, interaction is essential for many models and provides an
intermediate step in developing FeynGrav.

Lastly, we address the quadratic gravity studied in classical papers [34, 35,
36, 37, 38, 39, 40]. We do not discuss the physical content of the model and
do not address the issues related to its applicability. We focus on the
derivation of the Feynman rules of the model. Their derivation follows the
same principal scheme and only uses more sophisticated expressions for higher
curvature invariants.

The paper is organised as follows. In Section 2, we discuss the recursive
relations between various geometrical quantities. Section 3 discusses the
Horndeski gravity, and it is split into subsections for each class of
Horndeski interactions. We discuss the coupling to the Gauss-Bonnet term in
Section 3.5. Within the standard Horndeski parametrisation, it involves the
$\ln$ function of a scalar field, which makes it inconvenient to use. We
address this issue by studying this term in a more convenient parametrisation.
Section 4 addresses the massive graviton propagator and discusses the
applicability of FeynGrav for massive gravity. Section 5 discusses the axion-
like coupling between a scalar field and a single vector field. Section 6
discusses the interaction rules in quadratic gravity. Lastly, Section 7
discusses all the new features implemented in the latest version of FeynGrav.
Section 8 presents our conclusions and discusses further development of
FeynGrav.

##  2 Recursive relations

This section discusses recursive relations that hold for the perturbative
expansion about the flat spacetime. We are introducing a new notation system
that is better suited to our purpose. New notations are required to
differentiate between the tensors that define the structure of perturbative
expansions and those that represent the Feynman rules. The tensors in
perturbative expansions do not possess symmetry unless explicitly stated, so
we refer to them as plain tensors. We do not reserve a particular name for the
tensors present in the interaction rules but use cursive symbols to denote
them.

The metric perturbative expansion is defined as follows:

$\displaystyle
g_{\mu\nu}\overset{\text{def}}{=}\eta_{\mu\nu}+\kappa\,h_{\mu\nu}.$ (3)

Here $\kappa$ is the gravitational coupling with mass dimension $-1$ related
to the Newton constant $G_{\text{N}}$.

$\displaystyle\kappa^{2}\overset{\text{def}}{=}32\,\pi\,G_{\text{N}}.$ (4)

The plain $I_{(n)}$-tensor of the $n$-th order has $2n$ indices and defined as
follows:

$\displaystyle
I^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}_{(n)}\overset{\text{def}}{=}\eta^{\sigma_{1}\rho_{2}}\eta^{\sigma_{2}\rho_{3}}\cdots\eta^{\sigma_{n}\rho_{1}}.$
(5)

As stated above, this tensor does not have additional symmetries.

The $\mathcal{I}_{(n)}$-tensor of the $n$-th order has $2n$ indices are
defined as follows:

$\displaystyle\mathcal{I}_{(n)}^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}\overset{\text{def}}{=}\cfrac{1}{2^{n}}\,\cfrac{1}{n!}\left[I_{(n)}^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}+\text{permutations}\right].$
(6)

The permutation term ensures that the $\mathcal{I}$-tensor has additional
symmetries. It is constructed to be symmetric with respect to permutations of
indices within each index pair $\rho_{i}\leftrightarrow\sigma_{i}$, and with
respect to permutations of any two index pairs
$\\{\rho_{i},\sigma_{i}\\}\leftrightarrow\\{\rho_{j},\sigma_{j}\\}$.

The number of terms present in a tensor of the $n$-th order is $n!\,2^{n}$.
These terms are due to the tensor symmetry, so they cannot be negated.
Consequently, the expression for interaction vertices will also have a rapidly
growing number of terms, so the complexity of calculations will immensely grow
with each order of perturbation theory. This feature is expected since the
perturbative quantum gravity is an effective theory.

The plain $I$-tensors completely define the $g^{\mu\nu}$ perturbative
expansion:

$\displaystyle g^{\mu\nu}$
$\displaystyle=\sum\limits_{n=0}^{\infty}(-1)^{n}\,\kappa^{n}\,I^{\mu\nu\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}_{(1+n)}\,h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}\,.$ (7)

The reason is the structure of the $I$-tensor. By the construction, the tensor
is designed to contract indices and form traces:

$\displaystyle\begin{split}&(h^{n})^{\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}}\,,\\\
&\operatorname{tr}(h^{n})=I_{(n)}^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}\,h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}\,.\end{split}$ (8)

Therefore, the plain $I$ tensor is sufficient to completely encapsulate the
inverse metric’s perturbative structure.

The plain $C$-tensor of $n$-th order has $2n$ indices and describes the
structure of the volume factor $\sqrt{-g}$:

$\displaystyle\sqrt{-g}\overset{\text{def}}{=}\sum\limits_{n=0}^{\infty}\,\kappa^{n}\,C^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}_{(n)}\,h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}\,.$ (9)

This definition is non-constructive since we must describe how to generate the
tensor. The publication [2] presents an explicit formula for the plain
$C$-tensor, but it requires a complicated summation. We will proceed to obtain
a recurrent relation defining this tensor.

In full analogy with the previous case, the plain $C$-tensor does not have
additional symmetries. One shall also define its symmetric counterpart. The
$\mathcal{C}$-tensor of the $n$-th order in the following way:

$\displaystyle\mathcal{C}_{(n)}^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}\overset{\text{def}}{=}\cfrac{1}{2^{n}}\,\cfrac{1}{n!}\,\left[C_{(n)}^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}+\text{permutations}\right].$
(10)

Permutations account for all terms that make the $\mathcal{C}$-tensor
symmetric with respect to permutations of indices within each index pair
$\rho_{i}\leftrightarrow\sigma_{i}$ and with respect to permutations of any
two index pairs
$\\{\rho_{i},\sigma_{i}\\}\leftrightarrow\\{\rho_{j},\sigma_{j}\\}$.

The following relation allows one to define the plain $C$-tensor completely

$\displaystyle
C_{(n)}^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}=\cfrac{1}{2\,n}\sum\limits_{k=1}^{n}\,(-1)^{k-1}\,I_{(k)}^{\rho_{1}\sigma_{1}\cdots\rho_{k}\sigma_{k}}C_{(n-k)}^{\rho_{k+1}\sigma_{k+1}\cdots\rho_{n}\sigma_{n}}.$
(11)

This relation is a consequence of a scaling behaviour. Let us define a new
metric with scaled perturbations:

$\displaystyle\mathbf{g}_{\mu\nu}=\eta_{\mu\nu}+z\,\kappa\,h_{\mu\nu}\,.$ (12)

In this expression, $z$ is an arbitrary positive real number that scales the
perturbation amplitude. This metric has no additional physical or mathematical
meaning; we only use it to derive the relations discussed.

The infinitesimal action of this scaling on the volume factor is

$\displaystyle\left.\cfrac{d}{dz}\sqrt{-\mathbf{g}}\leavevmode\nobreak\
\right|_{z=1},$ (13)

and we have two ways to calculate this quantity. Firstly, we can use the
Jacobi’s formula:

$\displaystyle\left.\cfrac{d}{dz}\sqrt{-\mathbf{g}}\leavevmode\nobreak\
\right|_{z=1}=\cfrac{\kappa}{2}\,\sqrt{-g}\,g^{\mu\nu}\,h_{\mu\nu}.$ (14)

Secondly, we can use the formula for the explicit perturbative structure:

$\displaystyle\left.\cfrac{d}{dz}\sqrt{-\mathbf{g}}\leavevmode\nobreak\
\right|_{z=1}=\sum\limits_{n=1}^{\infty}\,n\,\kappa^{n}\,C^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}_{(n)}\,h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}.$ (15)

These equations allows us to exclude $\mathbf{g}$ completely:

$\displaystyle\cfrac{\kappa}{2}\,\sqrt{-g}\,g^{\mu\nu}\,h_{\mu\nu}=\sum\limits_{n=1}^{\infty}\,n\,\kappa^{n}\,C^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}_{(n)}\,h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}.$ (16)

The left-hand side of this expression can also be calculated as a perturbative
expansion:

$\displaystyle\begin{split}\cfrac{\kappa}{2}\,\sqrt{-g}\,g^{\mu\nu}\,h_{\mu\nu}&=\\!\cfrac{\kappa}{2}\\!\sum\limits_{p_{1}=0}^{\infty}\kappa^{p_{1}}C^{\rho_{1}\sigma_{1}\cdots\rho_{p_{1}}\sigma_{p_{1}}}_{(p_{1})}h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{p_{1}}\sigma_{p_{1}}}\\!\\!\sum\limits_{p_{2}=0}^{\infty}\\!(-1)^{p_{2}}\kappa^{p_{2}}I^{\mu\nu\lambda_{1}\tau_{1}\cdots\lambda_{p_{2}}\tau_{p_{2}}}_{(p_{2}+1)}h_{\lambda_{1}\tau_{1}}\cdots
h_{\lambda_{p_{2}}\tau_{p_{2}}}h_{\mu\nu}\\\
&=\sum\limits_{n=1}^{\infty}\sum\limits_{k=1}^{n}\,\cfrac{(-1)^{k-1}}{2}\leavevmode\nobreak\
\kappa^{n}\,I_{(k)}^{\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}}\\\
&=\sum\limits_{n=1}^{\infty}\kappa^{n}\,\sum\limits_{k=1}^{n}\cfrac{(-1)^{k-1}}{2}\leavevmode\nobreak\
I_{(k)}^{\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}}.\end{split}$ (17)

Comparing these expressions, we obtain the desired recursive relation.

The relation (11) is recursive because it defines the plain $C$-tensor of
$n$-the order through plain $C$-tensors of lower orders. It would be much
preferable to obtain a relation that expresses $C_{(n)}$ through $C_{(n-1)}$
alone, but this is impossible. The formula (11) shows that all low orders must
contribute to the plain $C$-tensor. Moreover, because the volume factor
$\sqrt{-g}$ is a simple object, there are few relations that one can use to
search for similar recursive relations. Therefore, it is safe to assume that
no simpler formula exists.

In the same way, one can define tensors describing the perturbative structure
of the other factors involving $\sqrt{-g}$. Plain $C_{(1)}$, $C_{(2)}$,
$C_{(3)}$, and any other $C_{(n)}$-tensors are defined as follows:

$\displaystyle\begin{split}\sqrt{-g}\,g^{\mu\nu}\overset{\text{def}}{=}&\sum\limits_{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}}\,,\\\
\sqrt{-g}\,g^{\mu\nu}\,g^{\alpha\beta}\overset{\text{def}}{=}&\sum\limits_{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}}\,,\\\
\sqrt{-g}\,g^{\mu\nu}\,g^{\alpha\beta}\,g^{\rho\sigma}\overset{\text{def}}{=}&\sum\limits_{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}}\,,\\\ \sqrt{-g}\,g^{\mu_{1}\nu_{1}}\cdots
g^{\mu_{l}\nu_{l}}\overset{\text{def}}{=}&\sum\limits_{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}}\,.\end{split}$ (18)

The corresponding symmetric tensors are defined as follows.

$\displaystyle\mathcal{C}_{(l;n)}^{\mu_{1}\nu_{1}\cdots\mu_{l}\nu_{l},\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}\overset{\text{def}}{=}\cfrac{1}{2^{n}}\,\cfrac{1}{n!}\,\left[C_{(l;n)}^{\mu_{1}\nu_{1}\cdots\mu_{l}\nu_{l},\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}+\text{permutations}\right].$
(19)

Permutations account for all terms that make the $\mathcal{C}$-tensor
symmetric with respect to permutations of indices within each index pair
$\rho_{i}\leftrightarrow\sigma_{i}$ and with respect to permutations of any
two index pairs
$\\{\rho_{i},\sigma_{i}\\}\leftrightarrow\\{\rho_{j},\sigma_{j}\\}$. Let us
highlight that the symmetrisation is performed only for indices contracted
with the perturbation indices.

These definitions are also non-constructive, but it is possible to obtain the
following recursive relation that defines them completely:

$\displaystyle
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\limits_{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}}.$
(20)

This relation is not based on scaling but on a simple series multiplication.
On the one hand, one can use the $C_{(l)}$ tensor definition:

$\displaystyle\sqrt{-g}\,g^{\mu_{1}\nu_{1}}\cdots
g^{\mu_{l}\nu_{l}}\overset{\text{def}}{=}$
$\displaystyle\sum\limits_{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}}\,.$ (21)

On the other hand, one can separate one metric and expand it separately:

$\displaystyle\begin{split}&\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}}\\\
&=\sum\limits_{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\limits_{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}}}\\\
&=\sum\limits_{n_{1}=0}^{\infty}\sum\limits_{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}}}.\end{split}$ (22)

One shall change the summation order in this expression to obtain the
discussed expression:

$\displaystyle\begin{split}&\sqrt{-g}\,g^{\mu_{1}\nu_{1}}\,g^{\mu_{2}\nu_{2}}\cdots
g^{\mu_{l}\nu_{l}}=\\\
&=\sum\limits_{n=0}^{\infty}\sum\limits_{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{split}$ (23)

Similar definitions and recursive relations also hold for factors involving
the vierbein $\mathfrak{e}_{\mu}{}^{\nu}$ [15]:

$\displaystyle\begin{split}\mathfrak{e}_{\mu}{}^{\nu}\overset{\text{def}}{=}&\sum\limits_{n=0}^{\infty}\kappa^{n}E_{\mu\leavevmode\nobreak\
(1+n)}^{\leavevmode\nobreak\ \leavevmode\nobreak\
\nu\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}\\\
=&\sum\limits_{n=0}^{\infty}\,\kappa^{n}\binom{-\frac{1}{2}}{\leavevmode\nobreak\
n}I_{\mu\leavevmode\nobreak\ (1+n)}^{\leavevmode\nobreak\ \leavevmode\nobreak\
\nu\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}\,,\\\
\sqrt{-g}\,\mathfrak{e}_{\mu}{}^{\nu}\overset{\text{def}}{=}&\sum\limits_{n=0}^{\infty}\kappa^{n}{C_{E}}_{\mu\leavevmode\nobreak\
(n)}^{\leavevmode\nobreak\
\nu,\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}\\\
=&\sum\limits_{n=0}^{\infty}\,\kappa^{n}\leavevmode\nobreak\
\sum\limits_{p=0}^{n}\binom{-\frac{1}{2}}{\leavevmode\nobreak\
n}C_{(p)}^{\rho_{1}\sigma_{1}\cdots\rho_{p}\sigma_{p}}I_{\mu\leavevmode\nobreak\
(1+n-p)}^{\leavevmode\nobreak\ \leavevmode\nobreak\
\nu\rho_{p+1}\sigma_{p+1}\cdots\rho_{n}\sigma_{n}}h_{\rho_{1}\sigma_{1}}\cdots
h_{\rho_{n}\sigma_{n}}\,.\end{split}$ (24)

Here $\binom{x}{y}$ is the binomial coefficient. Only $C_{E}$ tensor can enter
a Lagrangian, and its symmetric counterpart reads

$\displaystyle{\mathcal{C}_{E}}_{\mu\leavevmode\nobreak\
(n)}^{\leavevmode\nobreak\ \leavevmode\nobreak\
\nu,\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}\overset{\text{def}}{=}\cfrac{1}{2^{n}}\,\cfrac{1}{n!}\,\left[{C_{E}}_{\mu\leavevmode\nobreak\
(n)}^{\leavevmode\nobreak\
\nu,\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}+\text{permutations}\right].$
(25)

Permutations account for all terms that make the $\mathcal{C}_{E}$-tensor
symmetric with respect to permutations of indices within each index pair
$\rho_{i}\leftrightarrow\sigma_{i}$, and with respect to permutations of any
two index pairs
$\\{\rho_{i},\sigma_{i}\\}\leftrightarrow\\{\rho_{j},\sigma_{j}\\}$.

A similar recurrent relation holds for the $C_{E}$ tensor and its proof is
similar to the $C_{(l)}$ case:

$\displaystyle{C_{E}}_{\mu\leavevmode\nobreak\ (n)}^{\leavevmode\nobreak\
\nu,\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}=\sum\limits_{p=0}^{n}E_{\mu\leavevmode\nobreak\
(1+p)}^{\leavevmode\nobreak\ \leavevmode\nobreak\
\nu\rho_{1}\sigma_{1}\cdots\rho_{p}\sigma_{p}}C_{(n-p)}^{\rho_{p+1}\sigma_{p+1}\cdots\rho_{n}\sigma_{n}}.$
(26)

Relations (11), (20), and (25) are the main results of this Section. They
provide an alternative way to calculate the perturbative structure of the
corresponding objects. They improve computational efficiency, making them
essential for the higher-order perturbation theory. The new version of
FeynGrav uses these formulae.

Lastly, similar relations hold for an arbitrary metric. One shall only replace
the flat metric $\eta_{\mu\nu}$ with the arbitrary background
$\overline{g}_{\mu\nu}$ in the definition of all tensors. However, in that
more sophisticated case, the recurrent relation may not improve computational
efficiency that strongly. Since FeynGrav operates with FeynCalc, one does not
operate with metric components but with the particular objects defined within
FeynCalc. To put it otherwise, FeynCalc does not directly refer to the metric
component but operates with object index contraction. The discussed recursive
relations involve fewer operations, which is the reason behind their
efficiency. The same logic may not hold for perturbations about an arbitrary
metric, so the efficiency is not guaranteed.

##  3 Horndeski gravity

The Horndeski theory is the most general scalar-tensor theory without coupling
between the scalar field and the regular matter that admits second-order field
equations. The theory was discovered in [41] and independently rediscovered in
[42]. The second differential order of field equations ensures the theory is
free from the Ostrogradsky instability [43, 44, 45].

The theory is constituted by the four following Lagrangians:

$\displaystyle\mathcal{A}=\int
d^{4}x\sqrt{-g}\left[\mathcal{L}_{2}+\mathcal{L}_{3}+\mathcal{L}_{4}+\mathcal{L}_{5}\right],$
(27) $\displaystyle\mathcal{L}_{2}$ $\displaystyle=G_{2}(\phi,X)\,,$
$\displaystyle\mathcal{L}_{3}$ $\displaystyle=G_{3}(\phi,X)\,\square\phi\,,$
$\displaystyle\mathcal{L}_{4}$
$\displaystyle=G_{4}(\phi,X)\,R+G_{4,X}\left[\left(\square\phi\right)^{2}-\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}\right],$
$\displaystyle\mathcal{L}_{5}$
$\displaystyle=G_{5}(\phi,X)\,G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi-\cfrac{1}{6}\,G_{5,X}\left[\left(\square\phi\right)^{3}-3\,\square\phi\,\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}+2\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{3}\right].$

Here $G_{i}=G_{i}(\phi,X)$ are functions of the scalar field $\phi$ and its
kinetic term $X=g^{\mu\nu}\,\nabla_{\mu}\phi\,\nabla_{\nu}\phi$, while
$G_{4,X}$ and $G_{5,X}$ note derivatives with respect to the kinetic term.

The theory shows that there are only four couplings with the desired features.
Interactions of the $G_{2}$ class correspond to the minimal coupling to the
gravity of non-minimally interacting scalar fields. The coupling is minimal
since it does not involve Christoffel symbols. Interaction of the $G_{3}$
class corresponds to the simplest non-minimal coupling to gravity. The
presence of a Christoffel symbol shows that the scalar field couples to a
graviton momentum. The coupling is non-minimal since it involves a single
Christoffel symbol, but it is minimal because only one symbol is present. The
Christoffel symbol can be removed if $G_{3,X}=0$, but in that case, the
interaction becomes redundant as it reduces to the $G_{2}$ case with an
integration by parts. Interactions of $G_{4}$ and $G_{5}$ classes involve up
to two and three Christoffel symbols.

Counting the Christoffel symbols shows that the Horndeski theory can be
understood as a theory of coupling between a scalar field and graviton
momenta. A healthy coupling cannot involve more than three graviton momenta.
In turn, such coupling is strongly constrained, and only the Horndeski
interactions are allowed.

These considerations make the Horndeski theory an essential part of
perturbative quantum gravity. This paper considers all the Horndeski
interactions and presents explicit formulae for the corresponding perturbative
expansions. Since the interactions involve many derivatives, the formulae
occupy much space and are computationally challenging. Consequently, despite
having all the required analytic expressions, we implement only a few selected
cases in the new version of FeynGrav, as discussed in Section 7. It shall also
be noted that the $G_{5}$ interaction class naturally involves the standard
Einstein-Hilbert term describing pure general relativity. For the sake of
simplicity, we separate this term and discuss only contributions describing
interactions. The theory was already briefly discussed in [4]. This paper
addresses it in more detail and also discusses its implementations.

###  3.1 Horndeski $G_{2}(\phi,X)$ class

The interactions of the $G_{2}$ class of Horndeski theory describe the minimal
coupling between the scalar field and gravity. The coupling is minimal because
it does not involve graviton momenta. The corresponding Feynman rules are
obtained as follows.

First and foremost, the function $G_{2}$ shall be smooth enough to admit a
power series expansions:

$\displaystyle\begin{split}G_{2}(\phi,X)&=\sum\limits_{a=0}^{\infty}\sum\limits_{b=0}^{\infty}\lambda_{(a,b)}\,\phi^{a}\,X^{b}\\\
&=\sum\limits_{a=0}^{\infty}\sum\limits_{b=0}^{\infty}\lambda_{(a,b)}\,\phi^{a}\,g^{\alpha_{1}\beta_{1}}\cdots
g^{\alpha_{b}\beta_{b}}\partial_{\alpha_{1}}\phi\,\partial_{\beta_{1}}\phi\cdots\partial_{\alpha_{b}}\phi\,\partial_{\beta_{b}}\phi.\end{split}$
(28)

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

$\displaystyle\int d^{4}x\sqrt{-g}\,G_{2}(\phi,X)\to\int
d^{4}x\sqrt{-g}\,\lambda_{(a,b)}\,\phi^{a}\,X^{b}\,.$ (29)

Secondly, one shall use the standard Fourier transformations to describe the
momentum structure of the interaction.

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\,\lambda_{(a,b)}\,\phi^{a}\,X^{b}=\int
d^{4}x\sqrt{-g}\,g^{\alpha_{1}\beta_{1}}\cdots
g^{\alpha_{b}\beta_{b}}\,\phi^{a}\,\partial_{\alpha_{1}}\phi\,\partial_{\beta_{1}}\phi\cdots\partial_{\alpha_{b}}\phi\partial_{\beta_{b}}\phi\\\
&=\sum\limits_{n=0}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\,\prod\limits_{j=1}^{a+2b}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{i}\right)\\\ &\hskip
10.0pt\times\kappa^{n}\,\lambda_{(a,b)}\leavevmode\nobreak\
\mathcal{C}_{(b;n)}^{\alpha_{1}\beta_{1}\cdots\alpha_{b}\beta_{b},\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-1})_{\alpha_{b}}(p_{2b})_{\beta_{b}}\,.\end{split}$
(30)

Many formulae in this paper rely on the Fourier transformation. We do not
discuss their derivation in detail since they are only required to perform the
transformation correctly. Similar formulae are discussed in detail in previous
publications [2, 3].

Lastly, one shall construct the rule corresponding Feynman rules:

$\displaystyle\begin{gathered}\fmffile{Diags/G2_{v}ertex}\fmfgraph*(40,40)\fmfleft{L1,L2}\fmfright{R1,R2}\fmf{dbl_{w}iggly}{L1,V}\fmf{dbl_{w}iggly}{L2,V}\fmf{dashes}{V,R1}\fmf{dashes}{V,R2}\fmfdot{V}\fmffreeze\fmf{dots}{L1,L2}\fmf{dots}{R1,R2}\fmflabel{$\rho_1\sigma_1,k_1$}{L1}\fmflabel{$\rho_n\sigma_n,k_n$}{L2}\fmflabel{$p_1$}{R1}\fmflabel{$p_{a+2b}$}{R2}\fmflabel{$\lambda_{(a,b)}$}{V}\fmfv{$\lambda_{(a,b)}$,label.angle=90,label.dist=15pt}{V}\end{gathered}\hskip
20.0pt\begin{split}&=i\,\kappa^{n}\,(-1)^{b}\,\lambda_{(a,b)}\,\mathcal{C}_{(b;n)}^{\alpha_{1}\beta_{1}\cdots\alpha_{b}\beta_{b},\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-1})_{\alpha_{b}}(p_{2b})_{\beta_{b}}\\\
&\hskip 10.0pt+\text{permutations}.\end{split}$ (31)

In this expression, the permutation term accounts for all possible
permutations of $a+2\,b$ scalars and $n$ gravitons. When $b=0$, the
interaction reduces to the standard scalar field potential, considered in [3].

###  3.2 Horndeski $G_{3}(\phi,X)$ class

Interactions of $G_{3}$ class describe the simplest non-minimal coupling
between the scalar field and gravity. The interaction involves a single
graviton momentum due to the presence of the Christoffel symbol. The presence
of the Christoffel symbol can be excluded using integration by parts in an
exceptional case $G_{3,X}=0$. However, in that case, the interaction is also
reduced to the $G_{2}$ class, making this case redundant. Because of this
reason, we implicitly assume $G_{3,X}\not=0$ within this section.

We shall separate the Christoffel symbols explicitly to proceed with deriving
interaction rules.

$\displaystyle\begin{split}\int d^{4}x\sqrt{-g}\,G_{3}\,\square\phi=\int
d^{4}x\sqrt{-g}\,g^{\mu\nu}G_{3}\,\partial_{\mu}\partial_{\nu}\phi-\int
d^{4}x\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}\,\Gamma_{\alpha\mu\nu}G_{3}\,\partial_{\beta}\phi.\end{split}$
(32)

Similar to the previous case, $G_{3}$ should be smooth enough to admit a power
series expansion:

$\displaystyle
G_{3}(\phi,X)=\sum\limits_{a=1}^{\infty}\sum\limits_{b=1}^{\infty}\,\Theta_{(a,b)}\,\phi^{a}\,X^{b}.$
(33)

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

$\displaystyle\int d^{4}x\sqrt{-g}\,G_{3}(\phi,X)\,\square\phi\to\int
d^{4}x\sqrt{-g}\,\Theta_{(a,b)}\,\phi^{a}\,X^{b}\,\square\phi.$ (34)

Here, $\Theta_{(a,b)}$ is a dimensional coupling with the mass dimension
$4(1-b)-a-3$.

After this simplification, one can make the perturbative structure of the
interaction explicit.

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\,\Theta_{(a,b)}\,\phi^{a}\,X^{b}\,\square\phi\\\ &=\int
d^{4}x\sqrt{-g}\,g^{\mu\nu}\Theta_{(a,b)}\,\phi^{a}\,X^{b}\,\partial_{\mu}\partial_{\nu}\phi-\int
d^{4}x\sqrt{-g}\,g^{\mu\nu}g^{\rho\sigma}\Gamma_{\rho\mu\nu}\Theta_{(a,b)}\,\phi^{a}\,X^{b}\,\partial_{\sigma}\phi\\\
&=\sum\limits_{n=0}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\\!\\!\\!\\!\prod\limits_{j=1}^{a+2b+1}\\!\\!\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\leavevmode\nobreak\
(2\pi)^{4}\delta\left(\sum k_{i}+\sum p_{i}\right)\\\ &\hskip
10.0pt\times\kappa^{n}(-1)^{b+1}\Theta_{(a,b)}\mathcal{C}_{(1+b;n)}^{\mu\nu\alpha_{1}\beta_{1}\cdots\alpha_{b}\beta_{b},\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-1})_{\alpha_{b}}(p_{2b})_{\beta_{b}}(p_{2b+1})_{\mu}(p_{2b+1})_{\nu}\\\
&-\sum\limits_{n=1}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\\!\\!\\!\\!\prod\limits_{j=1}^{a+2b+1}\\!\\!\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\leavevmode\nobreak\
(2\pi)^{4}\delta\left(\sum k_{i}+\sum p_{i}\right)\\\ &\hskip
10.0pt\times\kappa^{n}(-1)^{b+1}\Theta_{(a,b)}\mathcal{C}_{(2+b;n-1)}^{\mu\nu\rho\sigma\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n-1}}(\Gamma_{\rho\mu\nu})^{\lambda\rho_{n}\sigma_{n}}(k_{n})_{\lambda}(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-1})_{\alpha_{b}}(p_{2b})_{\beta_{b}}(p_{2b+1})_{\sigma}.\end{split}$
(35)

The first part of this expression contributes at $\kappa^{0}$ level and
describes the self-interaction of the scalar field. The second part does not
contribute at $\kappa^{0}$ level and only describes the coupling between
scalars and gravitons. Similar features are typical for coupling involving a
few Christoffel symbols. Such interaction terms do not contribute to low
orders of $\kappa$. We highlight this feature wherever it occurs. We also
include the terms in all formulas to make them more compact.

Finally, we obtain the interaction rule by appropriately symmetrising the
expression that we have obtained.

$\displaystyle\begin{gathered}\fmffile{Diags/G3_{v}ertex}\fmfgraph*(40,40)\fmfleft{L1,L2}\fmfright{R1,R2}\fmf{dbl_{w}iggly}{L1,V}\fmf{dbl_{w}iggly}{L2,V}\fmf{dashes}{V,R1}\fmf{dashes}{V,R2}\fmfdot{V}\fmffreeze\fmf{dots}{L1,L2}\fmf{dots}{R1,R2}\fmflabel{$\rho_1\sigma_1,k_1$}{L1}\fmflabel{$\rho_n\sigma_n,k_n$}{L2}\fmflabel{$p_1$}{R1}\fmflabel{$p_{a+2b+1}$}{R2}\fmflabel{$\Theta_{(a,b)}$}{V}\fmfv{$\Theta_{(a,b)}$,label.angle=90,label.dist=15pt}{V}\end{gathered}\hskip
40.0pt\begin{split}=&i\,\kappa^{n}(-1)^{b+1}\Theta_{(a,b)}\Bigg{[}\mathcal{C}_{(1+b;n)}^{\mu\nu\alpha_{1}\beta_{1}\cdots\alpha_{b}\beta_{b},\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}(p_{2b+1})_{\mu}(p_{2b+1})_{\nu}\\\
&\hskip
20.0pt-\mathcal{C}_{(2+b;n-1)}^{\mu\nu\rho\sigma\alpha_{1}\cdots\beta_{b},\rho_{1}\sigma_{1}\cdots\sigma_{n-1}\rho_{n-1}}(\Gamma_{\rho\mu\nu})^{\lambda\rho_{n}\sigma_{n}}(k_{n})_{\lambda}(p_{2b+1})_{\sigma}\Bigg{]}\\\
&\times(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-1})_{\alpha_{b}}(p_{2b})_{\beta_{b}}+\text{permutations}.\end{split}$
(37)

The permutation term includes all possible permutations of $2\,b+a+1$ scalar
fields and $n$ gravitons. Let us highlight again that only the term in this
expression contributes at $\kappa^{0}$ level.

###  3.3 Horndeski $G_{4}(\phi,X)$ class

Interactions of $G_{4}$ class take a special place within the Horndeski theory
since pure general relativity belongs to this class. The pure Einstein-Hilbert
term was considered in previous publications [2, 3], so we will not consider
it in this section.

The exclude the Einstein-Hilbert term, we assume that the $G_{4}$ function
admits the following power series expansion:

$\displaystyle
G_{4}=-\cfrac{2}{\kappa^{2}}+\sum\limits_{a,b}\Upsilon_{(a,b)}\phi^{a}\,X^{b}.$
(38)

Here $\Upsilon_{(a,b)}$ is a dimensional coupling with the mass dimension
$4(1-b)-a-2$. To account only for the interaction terms, we study only the
following term without the loss of generality:

$\displaystyle G_{4}\to\Upsilon_{(a,b)}\phi^{a}\,X^{b}.$ (39)

Further, we shall make the Christoffel symbols explicit with the following
formula:

$\displaystyle
R=g^{\mu\nu}g^{\alpha\beta}\partial_{\mu}\left[\Gamma_{\nu\alpha\beta}-\Gamma_{\alpha\beta\nu}\right]+g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}\left[\Gamma_{\mu\alpha\rho}\Gamma_{\nu\beta\sigma}-\Gamma_{\mu\alpha\beta}\Gamma_{\nu\rho\sigma}\right].$
(40)

In turn, this allows us to obtain the following explicit expression for the
interaction Lagrangian:

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\,\Upsilon_{(a,b)}\,\phi^{a}\,X^{b-1}\left[XR\\!+\\!b\left(\\!\left(\square\phi\right)^{2}\\!-\\!\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}\right)\right]\\\
&=\Upsilon_{(a,b)}\int
d^{4}x\Bigg{[}\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}\partial_{\mu}\left[\Gamma_{\nu\alpha\beta}-\Gamma_{\alpha\beta\nu}\right]\phi^{a}X^{b}+\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}\left[\Gamma_{\mu\alpha\rho}\Gamma_{\nu\sigma\beta}-\Gamma_{\mu\alpha\beta}\Gamma_{\nu\rho\sigma}\right]\phi^{a}X^{b}\\\
&\hskip
70.0pt+b\sqrt{-g}\left(g^{\mu\nu}g^{\alpha\beta}-g^{\mu\alpha}g^{\nu\beta}\right)\nabla_{\mu}\nabla_{\nu}\phi\nabla_{\alpha}\nabla_{\beta}\phi\,\phi^{a}X^{b-1}\Bigg{]}\\\
&=\Upsilon_{(a,b)}\int
d^{4}x\Bigg{[}\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}\partial_{\mu}\left[\Gamma_{\nu\alpha\beta}-\Gamma_{\alpha\beta\nu}\right]\phi^{a}X^{b}+\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}\left[\Gamma_{\mu\alpha\rho}\Gamma_{\nu\sigma\beta}-\Gamma_{\mu\alpha\beta}\Gamma_{\nu\rho\sigma}\right]\phi^{a}X^{b}\\\
&\hskip
70.0pt+b\sqrt{-g}\left(g^{\mu\nu}g^{\alpha\beta}-g^{\mu\alpha}g^{\nu\beta}\right)\partial_{\mu}\partial_{\nu}\phi\,\partial_{\alpha}\partial_{\beta}\phi\,\phi^{a}X^{b-1}\\\
&\hskip
70.0pt-2\,b\sqrt{-g}\left(g^{\mu\nu}g^{\alpha\beta}-g^{\mu\alpha}g^{\nu\beta}\right)g^{\rho\sigma}\Gamma_{\rho\mu\nu}\,\partial_{\sigma}\phi\,\partial_{\alpha}\partial_{\beta}\phi\,\phi^{a}X^{b-1}\\\
&\hskip
70.0pt+b\sqrt{-g}\left(g^{\mu\nu}g^{\alpha\beta}-g^{\mu\alpha}g^{\nu\beta}\right)g^{\rho\sigma}g^{\lambda\tau}\Gamma_{\rho\mu\nu}\Gamma_{\lambda\alpha\beta}\partial_{\sigma}\phi\,\partial_{\tau}\phi\,\phi^{a}X^{b-1}\Bigg{]}.\end{split}$
(41)

Lastly, one can make the perturbative structure of this interaction explicit.

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\Upsilon_{(a,b)}\,\phi^{a}\,X^{b-1}\left[XR\\!+\\!b\left(\\!\left(\square\phi\right)^{2}\\!-\\!\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}\right)\right]\\\
&=\sum\limits_{n=1}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)\\\ &\hskip
8.0pt\times\kappa^{n}\Upsilon_{(a,b)}(-1)^{b+1}\mathcal{C}_{(2+b,n-1)}^{\mu\nu\alpha\beta\alpha_{1}\beta_{1}\cdots\alpha_{b}\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!1}}(k_{n})_{\mu}(k_{n})_{\lambda}\left[\left(\Gamma_{\nu\alpha\beta}\right)^{\lambda\rho_{n}\sigma_{n}}-\left(\Gamma_{\alpha\beta\nu}\right)^{\lambda\rho_{n}\sigma_{n}}\right](p_{1})_{\alpha_{1}}\cdots(p_{2b})_{\beta_{b}}\end{split}$
(42)
$\displaystyle+\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\kappa^{n}\Upsilon_{(a,b)}(-1)^{b+1}\mathcal{C}_{(3+b,n-2)}^{\mu\nu\alpha\beta\rho\sigma\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!2}}(k_{n\\!-\\!1})_{\lambda_{1}}(k_{n})_{\lambda_{2}}$
$\displaystyle\hskip
16.0pt\times\left[\left(\Gamma_{\mu\alpha\rho}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\nu\sigma\beta}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}-\left(\Gamma_{\mu\alpha\beta}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\nu\rho\sigma}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\right](p_{1})_{\alpha_{1}}\cdots(p_{2b})_{\beta_{b}}$
$\displaystyle+\sum\limits_{n=0}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\kappa^{n}\Upsilon_{(a,b)}\,b\,(-1)^{b+1}\left[\mathcal{C}_{(1+b,n)}^{\mu\nu\alpha\beta\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n}}-\mathcal{C}_{(1+b,n)}^{\mu\alpha\nu\beta\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n}}\right]$
$\displaystyle\hskip
16.0pt\times(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-3})_{\alpha_{b-1}}(p_{2b-2})_{\beta_{b-1}}(p_{2b-1})_{\mu}(p_{2b-1})_{\nu}(p_{2b})_{\alpha}(p_{2b})_{\beta}$
$\displaystyle-\sum\limits_{n=1}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\kappa^{n}\Upsilon_{(a,b)}\,2\,b\,(-1)^{b+1}\left[\mathcal{C}_{(2+b,n-1)}^{\mu\nu\alpha\beta\rho\sigma\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!1}}-\mathcal{C}_{(2+b,n-1)}^{\mu\alpha\nu\beta\rho\sigma\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!1}}\right](k_{n})_{\lambda}\left(\Gamma_{\rho\mu\nu}\right)^{\lambda\rho_{n}\sigma_{n}}$
$\displaystyle\hskip
16.0pt\times(p_{1})_{\alpha_{1}}\cdots(p_{2b-2})_{\beta_{b-1}}(p_{2b-1})_{\sigma}(p_{2b})_{\alpha}(p_{2b})_{\beta}$
$\displaystyle+\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\kappa^{n}\Upsilon_{(a,b)}\,b\,(-1)^{b+1}\left[\mathcal{C}_{(3+b,n-2)}^{\mu\nu\alpha\beta\rho\sigma\lambda\tau\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!2}}-\mathcal{C}_{(3+b,n-2)}^{\mu\alpha\nu\beta\rho\sigma\lambda\tau\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!2}}\right]\left(\Gamma_{\rho\mu\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\lambda\alpha\beta}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}$
$\displaystyle\hskip
16.0pt\times(k_{n\\!-\\!1})_{\lambda_{1}}(k_{n})_{\lambda_{2}}(p_{1})_{\alpha_{1}}\cdots(p_{2b-2})_{\beta_{b-1}}(p_{2b-1})_{\sigma}(p_{2b})_{\tau}.$

Only the third term contributes at $\kappa^{0}$ level in this expression. The
first and fourth terms contribute at the $\kappa^{1}$ level, and the second
and fifth terms contribute at the $\kappa^{2}$ level. Equivalently, only the
second and fifth terms describe the coupling of a scalar with two graviton
momenta, while the other terms ensure that the interaction is healthy.

This expression corresponds to the following interaction rule

$\displaystyle\begin{split}&\hskip
40.0pt\begin{gathered}\fmffile{Diags/G4_{v}ertex}\fmfgraph*(40,40)\fmfleft{L1,L2}\fmfright{R1,R2}\fmf{dbl_{w}iggly}{L1,V}\fmf{dbl_{w}iggly}{L2,V}\fmf{dashes}{V,R1}\fmf{dashes}{V,R2}\fmfdot{V}\fmffreeze\fmf{dots}{L1,L2}\fmf{dots}{R1,R2}\fmflabel{$\rho_1\sigma_1,k_1$}{L1}\fmflabel{$\rho_n\sigma_n,k_n$}{L2}\fmflabel{$p_1$}{R1}\fmflabel{$p_{a+2b}$}{R2}\fmflabel{$\Upsilon_{a,b}$}{V}\fmfv{$\Upsilon_{a,b}$,label.angle=90,label.dist=15pt}{V}\end{gathered}=i\,\kappa^{n}\,(-1)^{b+1}\Upsilon_{(a,b)}(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-1})_{\alpha_{b}}(p_{2b})_{\beta_{b}}\\\
\\\
&\times\\!\\!\Bigg{[}\mathcal{C}_{(2+b,n-1)}^{\mu\nu\alpha\beta\alpha_{1}\beta_{1}\cdots\alpha_{b}\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!1}}(k_{n})_{\mu}(k_{n})_{\lambda}\left[\left(\Gamma_{\nu\alpha\beta}\right)^{\lambda\rho_{n}\sigma_{n}}-\left(\Gamma_{\alpha\beta\nu}\right)^{\lambda\rho_{n}\sigma_{n}}\right]\\\
&\hskip
8.0pt+\mathcal{C}_{(3+b,n-2)}^{\mu\nu\alpha\beta\rho\sigma\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!2}}(k_{n\\!-\\!1})_{\lambda_{1}}(k_{n})_{\lambda_{2}}\left[\left(\Gamma_{\mu\alpha\rho}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\nu\sigma\beta}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\\!\\!-\\!\left(\Gamma_{\mu\alpha\beta}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\nu\rho\sigma}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\right]\\!\Bigg{]}\\\
&+i\,\kappa^{n}\,(-1)^{b+1}\Upsilon_{(a,b)}\,b\,(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-3})_{\alpha_{b-1}}(p_{2b-2})_{\beta_{b-1}}\\\
&\times\Bigg{[}\left[\mathcal{C}_{(1+b,n)}^{\mu\nu\alpha\beta\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n}}-\mathcal{C}_{(1+b,n)}^{\mu\alpha\nu\beta\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n}}\right](p_{2b-1})_{\mu}(p_{2b-1})_{\nu}(p_{2b})_{\alpha}(p_{2b})_{\beta}\\\
&\hskip
8.0pt-2\,\left[\mathcal{C}_{(2+b,n-1)}^{\mu\nu\alpha\beta\rho\sigma\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!1}}-\mathcal{C}_{(2+b,n-1)}^{\mu\alpha\nu\beta\rho\sigma\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!1}}\right](k_{n})_{\lambda}\left(\Gamma_{\rho\mu\nu}\right)^{\lambda\rho_{n}\sigma_{n}}(p_{2b-1})_{\sigma}(p_{2b})_{\alpha}(p_{2b})_{\beta}\\\
&\hskip
8.0pt+\left[\mathcal{C}_{(3+b,n-2)}^{\mu\nu\alpha\beta\rho\sigma\lambda\tau\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!2}}-\mathcal{C}_{(3+b,n-2)}^{\mu\alpha\nu\beta\rho\sigma\lambda\tau\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!2}}\right]\left(\Gamma_{\rho\mu\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\lambda\alpha\beta}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\\\
&\hskip
16.0pt\times(k_{n\\!-\\!1})_{\lambda_{1}}(k_{n})_{\lambda_{2}}(p_{2b-1})_{\sigma}(p_{2b})_{\tau}\Bigg{]}+\text{permutations}.\end{split}$
(43)

In this expression, the permutation term accounts for all possible
permutations of both scalars and gravitons. As noted above, we included all
the terms, even though the second and fifth terms begin to contribute at
$\kappa^{2}$ order.

###  3.4 Horndeski $G_{5}(\phi,X)$ class

The last interaction class of the Horndeski gravity describes the most
sophisticated non-minimal coupling involving three graviton momenta. We shall
simplify it in complete analogy with the previous cases.

Firstly, we assume that the coupling function $G_{5}$ is smooth enough to
admit the power series expansion:

$\displaystyle G_{5}=\sum\limits_{a,b}\Psi_{(a,b)}\phi^{a}X^{b}.$ (44)

Here, $\Psi_{(a,b)}$ is a dimensional coupling with the mass dimension
$4(1-b)-a-5$. Without the loss of generality, we use a single term for this
expansion and consider the following action:

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\left[G_{5}G^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\phi-\cfrac{1}{6}\,G_{5,X}\left\\{\left(\square\phi\right)^{3}-3\square\phi\,\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}+2\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{3}\right\\}\right]\\\
&\to\int
d^{4}x\sqrt{-g}\,\Psi_{a,b}\,\phi^{a}X^{b-1}\left[X\,G^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\phi-\cfrac{b}{6}\left\\{\left(\square\phi\right)^{3}-3\square\phi\,\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}+2\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{3}\right\\}\right].\end{split}$
(45)

Secondly, the action consists of two parts. The first part involves the
Einstein tensor, and it takes the following form:

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\,G_{\mu\nu}\,\nabla^{\mu}\nabla^{\nu}\phi=\cfrac{1}{2}\,\int
d^{4}x\,\sqrt{-g}\left[g^{\mu\alpha}g^{\nu\beta}+g^{\mu\beta}g^{\nu\alpha}-g^{\mu\nu}g^{\alpha\beta}\right]R_{\mu\nu}\nabla_{\alpha}\nabla_{\beta}\phi\\\
&=\cfrac{1}{2}\,\int
d^{4}x\sqrt{-g}\left[g^{\mu\alpha}g^{\nu\beta}+g^{\mu\beta}g^{\nu\alpha}-g^{\mu\nu}g^{\alpha\beta}\right]g^{\rho\sigma}\left(\partial_{\rho}\Gamma_{\sigma\mu\nu}-\partial_{\mu}\Gamma_{\rho\sigma\nu}\right)\,\nabla_{\alpha}\nabla_{\beta}\phi\\\
&\hskip 8.0pt+\cfrac{1}{2}\,\int
d^{4}x\sqrt{-g}\left[g^{\mu\alpha}g^{\nu\beta}+g^{\mu\beta}g^{\nu\alpha}-g^{\mu\nu}g^{\alpha\beta}\right]g^{\rho\sigma}g^{\lambda\tau}\left(\Gamma_{\rho\lambda\mu}\Gamma_{\sigma\tau\nu}-\Gamma_{\rho\lambda\tau}\Gamma_{\sigma\mu\nu}\right)\,\nabla_{\alpha}\nabla_{\beta}\phi\\\
&=\cfrac{1}{2}\,\int
d^{4}x\sqrt{-g}\left[g^{\mu\alpha}g^{\nu\beta}+g^{\mu\beta}g^{\nu\alpha}-g^{\mu\nu}g^{\alpha\beta}\right]g^{\rho\sigma}\left(\partial_{\rho}\Gamma_{\sigma\mu\nu}-\partial_{\mu}\Gamma_{\rho\sigma\nu}\right)\,\partial_{\alpha}\partial_{\beta}\phi\\\
&\hskip 8.0pt-\cfrac{1}{2}\,\int
d^{4}x\sqrt{-g}\left[g^{\mu\alpha}g^{\nu\beta}+g^{\mu\beta}g^{\nu\alpha}-g^{\mu\nu}g^{\alpha\beta}\right]g^{\rho\sigma}g^{\lambda\tau}\left(\partial_{\rho}\Gamma_{\sigma\mu\nu}-\partial_{\mu}\Gamma_{\rho\sigma\nu}\right)\Gamma_{\lambda\alpha\beta}\partial_{\tau}\phi\\\
&\hskip 8.0pt+\cfrac{1}{2}\,\int
d^{4}x\sqrt{-g}\left[g^{\mu\alpha}g^{\nu\beta}+g^{\mu\beta}g^{\nu\alpha}-g^{\mu\nu}g^{\alpha\beta}\right]g^{\rho\sigma}g^{\lambda\tau}\left(\Gamma_{\rho\lambda\mu}\Gamma_{\sigma\tau\nu}-\Gamma_{\rho\lambda\tau}\Gamma_{\sigma\mu\nu}\right)\,\partial_{\alpha}\partial_{\beta}\phi\\\
&\hskip 8.0pt-\cfrac{1}{2}\,\int
d^{4}x\sqrt{-g}\left[g^{\mu\alpha}g^{\nu\beta}+g^{\mu\beta}g^{\nu\alpha}-g^{\mu\nu}g^{\alpha\beta}\right]g^{\rho\sigma}g^{\lambda\tau}g^{\omega\epsilon}\left(\Gamma_{\rho\lambda\mu}\Gamma_{\sigma\tau\nu}-\Gamma_{\rho\lambda\tau}\Gamma_{\sigma\mu\nu}\right)\Gamma_{\omega\alpha\beta}\,\partial_{\epsilon}\phi.\end{split}$
(46)

The second part also contains Christoffel symbols due to the presence of
covariant derivatives and takes the following form:

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\,\left[(\square\phi)^{3}-3\square\phi\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}+2\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{3}\right]\\\
&=\int
d^{4}x\sqrt{-g}\,\left[g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}-3\,g^{\mu\nu}g^{\alpha\rho}g^{\beta\sigma}+2g^{\nu\alpha}g^{\beta\rho}g^{\sigma\mu}\right]\nabla_{\mu}\nabla_{\nu}\phi\,\nabla_{\alpha}\nabla_{\beta}\phi\,\nabla_{\rho}\nabla_{\sigma}\phi\\\
&=\int
d^{4}x\sqrt{-g}\,\left[g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}-3\,g^{\mu\nu}g^{\alpha\rho}g^{\beta\sigma}+2g^{\nu\alpha}g^{\beta\rho}g^{\sigma\mu}\right]\partial_{\mu}\partial_{\nu}\phi\,\partial_{\alpha}\partial_{\beta}\phi\,\partial_{\rho}\partial_{\sigma}\phi\\\
&\hskip 8.0pt-3\int
d^{4}x\sqrt{-g}\left[g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}-3\,g^{\mu\nu}g^{\alpha\rho}g^{\beta\sigma}+2g^{\nu\alpha}g^{\beta\rho}g^{\sigma\mu}\right]g^{\lambda\tau}\Gamma_{\lambda\mu\nu}\partial_{\tau}\phi\,\partial_{\alpha}\partial_{\beta}\phi\,\partial_{\rho}\partial_{\sigma}\phi\\\
&\hskip 8.0pt+3\int
d^{4}x\sqrt{-g}\left[g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}-3\,g^{\mu\nu}g^{\alpha\rho}g^{\beta\sigma}+2g^{\nu\alpha}g^{\beta\rho}g^{\sigma\mu}\right]g^{\lambda_{1}\tau_{1}}g^{\lambda_{2}\tau_{2}}\Gamma_{\lambda_{1}\mu\nu}\Gamma_{\lambda_{2}\alpha\beta}\partial_{\tau_{1}}\phi\,\partial_{\tau_{2}}\phi\,\partial_{\rho}\partial_{\sigma}\phi\\\
&\hskip 8.0pt-\\!\\!\int
d^{4}x\sqrt{-g}\left[g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}-3\,g^{\mu\nu}g^{\alpha\rho}g^{\beta\sigma}+2g^{\nu\alpha}g^{\beta\rho}g^{\sigma\mu}\right]g^{\lambda_{1}\tau_{1}}g^{\lambda_{2}\tau_{2}}g^{\lambda_{3}\tau_{3}}\Gamma_{\lambda_{1}\mu\nu}\Gamma_{\lambda_{2}\alpha\beta}\Gamma_{\lambda_{3}\rho\sigma}\partial_{\tau_{1}}\phi\,\partial_{\tau_{2}}\phi\,\partial_{\tau_{3}}\phi.\end{split}$
(47)

These expressions result in the following sophisticated expression describing
the perturbative structure of the interaction:

$\displaystyle\int
d^{4}x\sqrt{-g}\Psi_{a,b}\,\phi^{a}X^{b-1}\left[X\,G^{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi-\cfrac{b}{6}\left\\{\left(\square\phi\right)^{2}-3\square\phi\,\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}+2\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{3}\right\\}\right]$
(48)
$\displaystyle=\sum\limits_{n=1}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b+1}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\cfrac{1}{2}\,\kappa^{n}\Psi_{(a,b)}(-1)^{b+2}\left[\mathcal{C}_{3+b,n-1}^{\mu\alpha\nu\beta\rho\sigma\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!1}}+\cdots\right](k_{n})_{\lambda}\left((k_{n})_{\rho}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda\rho_{n}\sigma_{n}}-(k_{n})_{\mu}\left(\Gamma_{\rho\sigma\nu}\right)^{\lambda\rho_{n}\sigma_{n}}\right)$
$\displaystyle\hskip
16.0pt\times(p_{1})_{\alpha_{1}}\cdots(p_{2b})_{\beta_{b}}(p_{2b+1})_{\alpha}(p_{2b+1})_{\beta}$
$\displaystyle-\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b+1}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\cfrac{1}{2}\,\kappa^{n}\Psi_{(a,b)}(-1)^{b+2}\left[\mathcal{C}_{4+b,n-2}^{\mu\alpha\nu\beta\rho\sigma\lambda\tau\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!2}}+\cdots\right]$
$\displaystyle\hskip
16.0pt\times\\!\\!(k_{n\\!-\\!1})_{\lambda_{1}}\\!\\!\\!\left(\\!(k_{n\\!-\\!1})_{\rho}\\!\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\\!-\\!(k_{n\\!-\\!1})_{\mu}\\!\left(\Gamma_{\rho\sigma\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\right)\\!\\!(k_{n})_{\lambda_{2}}\\!\\!\left(\Gamma_{\lambda\alpha\beta}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}(p_{1})_{\alpha_{1}}\\!\\!\cdots(p_{2b})_{\beta_{b}}(p_{2b+1})_{\tau}$
$\displaystyle+\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b+1}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\cfrac{1}{2}\,\kappa^{n}\Psi_{(a,b)}(-1)^{b+2}\left[\mathcal{C}_{4+b,n-2}^{\mu\alpha\nu\beta\rho\sigma\lambda\tau\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!2}}+\cdots\right](k_{n-1})_{\lambda_{1}}(k_{n})_{\lambda_{2}}$
$\displaystyle\hskip
16.0pt\times\left[\left(\Gamma_{\rho\lambda\mu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\sigma\tau\nu}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}-\left(\Gamma_{\rho\lambda\tau}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\right](p_{1})_{\alpha_{1}}\cdots(p_{2b})_{\beta_{b}}(p_{2b+1})_{\alpha}(p_{2b+1})_{\beta}$
$\displaystyle-\sum\limits_{n=3}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b+1}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\cfrac{1}{2}\,\kappa^{n}\Psi_{(a,b)}(-1)^{b+2}\left[\mathcal{C}_{5+b,n-2}^{\mu\alpha\nu\beta\rho\sigma\lambda\tau\omega\epsilon\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!2}}+\cdots\right](k_{n\\!-\\!2})_{\lambda_{3}}(k_{n-1})_{\lambda_{1}}(k_{n})_{\lambda_{2}}\left(\Gamma_{\omega\alpha\beta}\right)^{\lambda_{3}\rho_{n\\!-\\!2}\sigma_{n\\!-\\!2}}$
$\displaystyle\hskip
16.0pt\times\left[\left(\Gamma_{\rho\lambda\mu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\sigma\tau\nu}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}-\left(\Gamma_{\rho\lambda\tau}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\right](p_{1})_{\alpha_{1}}\cdots(p_{2b})_{\beta_{b}}(p_{2b+1})_{\epsilon}$
$\displaystyle-\sum\limits_{n=0}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b+1}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\cfrac{b}{6}\,\kappa^{n}\Psi_{(a,b)}(-1)^{b+2}\left[\mathcal{C}_{2+b,n}^{\mu\nu\alpha\beta\rho\sigma\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n}}+\cdots\right](p_{1})_{\alpha_{1}}\cdots(p_{2b-2})_{\beta_{b-1}}$
$\displaystyle\hskip
16.0pt\times(p_{2b-1})_{\mu}(p_{2b-1})_{\nu}(p_{2b})_{\alpha}(p_{2b})_{\beta}(p_{2b+3})_{\rho}(p_{2b+3})_{\sigma}$
$\displaystyle+\sum\limits_{n=1}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b+1}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\cfrac{b}{2}\,\kappa^{n}\Psi_{(a,b)}(-1)^{b+2}\left[\mathcal{C}_{3+b,n-1}^{\mu\nu\alpha\beta\rho\sigma\omega\tau\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!1}}+\cdots\right](k_{n})_{\lambda}\left(\Gamma_{\omega\mu\nu}\right)^{\lambda\rho_{n}\sigma_{n}}(p_{1})_{\alpha_{1}}\cdots(p_{2b-2})_{\beta_{b-1}}$
$\displaystyle\hskip
16.0pt\times(p_{2b-1})_{\tau}(p_{2b})_{\alpha}(p_{2b})_{\beta}(p_{2b+1})_{\rho}(p_{2b+1})_{\sigma}$
$\displaystyle-\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b+1}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\\!\cfrac{b}{2}\,\kappa^{n}\Psi_{(a,b)}(-1)^{b+2}\\!\left[\mathcal{C}_{4+b,n-2}^{\mu\nu\alpha\beta\rho\sigma\omega_{1}\tau_{1}\omega_{2}\tau_{2}\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!2}}\\!+\\!\cdots\right]\\!(k_{n\\!-\\!1})_{\lambda_{1}}\\!\left(\Gamma_{\omega_{1}\mu\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}(k_{n})_{\lambda_{2}}\left(\Gamma_{\omega_{2}\alpha\beta}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}$
$\displaystyle\hskip
16.0pt\times(p_{1})_{\alpha_{1}}\cdots(p_{2b-2})_{\beta_{b-1}}(p_{2b-1})_{\tau_{1}}(p_{2b})_{\tau_{2}}(p_{2b+1})_{\rho}(p_{2b+1})_{\sigma}$
$\displaystyle+\sum\limits_{n=3}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{a+2b+1}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\,\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)$ $\displaystyle\hskip
8.0pt\times\cfrac{b}{6}\,\kappa^{n}\Psi_{(a,b)}(-1)^{b+2}\left[\mathcal{C}_{5+b,n-3}^{\mu\nu\alpha\beta\rho\sigma\omega_{1}\tau_{1}\omega_{2}\tau_{2}\omega_{3}\tau_{3}\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!3}}+\cdots\right](k_{n\\!-\\!2})_{\lambda_{1}}(k_{n\\!-\\!1})_{\lambda_{2}}(k_{n})_{\lambda_{3}}$
$\displaystyle\hskip
16.0pt\times\left(\Gamma_{\omega_{1}\mu\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!2}\sigma_{n\\!-\\!2}}\left(\Gamma_{\omega_{2}\alpha\beta}\right)^{\lambda_{2}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\omega_{3}\rho\sigma}\right)^{\lambda_{3}\rho_{n}\sigma_{n}}$
$\displaystyle\hskip
24.0pt\times(p_{1})_{\alpha_{1}}\cdots(p_{2b-2})_{\beta_{b-1}}(p_{2b-1})_{\tau_{1}}(p_{2b})_{\tau_{2}}(p_{2b+1})_{\tau_{3}}.$

This expression contains eight terms, and only term number five contributes at
$\kappa^{0}$ level. Terms number one and six contribute at $\kappa^{1}$ level.
Terms two, three, and seven contribute at $\kappa^{2}$ level, and only terms
four and eight contribute at $\kappa^{3}$ level.

We explicitly present the interaction rule expression for the sake of
generality despite it taking up much space.

$\displaystyle\begin{split}\\\ &\hskip
25.0pt\begin{gathered}\fmffile{Diags/G5_{v}ertex}\fmfgraph*(40,40)\fmfleft{L1,L2}\fmfright{R1,R2}\fmf{dbl_{w}iggly}{L1,V}\fmf{dbl_{w}iggly}{L2,V}\fmf{dashes}{V,R1}\fmf{dashes}{V,R2}\fmfdot{V}\fmffreeze\fmf{dots}{L1,L2}\fmf{dots}{R1,R2}\fmflabel{$\rho_1\sigma_1,k_1$}{L1}\fmflabel{$\rho_n\sigma_n,k_n$}{L2}\fmflabel{$p_1$}{R1}\fmflabel{$p_{a+2b+1}$}{R2}\fmflabel{$\Psi_{(a,b)}$}{V}\fmfv{$\Psi_{(a,b)}$,label.angle=90,label.dist=15pt}{V}\end{gathered}=i\,\cfrac{1}{2}\,\kappa^{n}\,\Psi_{(a,b)}(-1)^{b+2}(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-1})_{\alpha_{b}}(p_{2b})_{\beta_{b}}\\\
\\\
\times\Bigg{[}&\left[\mathcal{C}_{3+b,n-1}^{\mu\alpha\nu\beta\rho\sigma\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!1}}+\cdots\right](k_{n})_{\lambda}\left((k_{n})_{\rho}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda\rho_{n}\sigma_{n}}-(k_{n})_{\mu}\left(\Gamma_{\rho\sigma\nu}\right)^{\lambda\rho_{n}\sigma_{n}}\right)(p_{2b+1})_{\alpha}(p_{2b+1})_{\beta}\\\
&-\left[\mathcal{C}_{4+b,n-2}^{\mu\alpha\nu\beta\rho\sigma\lambda\tau\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!2}}+\cdots\right](k_{n\\!-\\!1})_{\lambda_{1}}\left(\\!(k_{n\\!-\\!1})_{\rho}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\\!-\\!(k_{n\\!-\\!1})_{\mu}\\!\left(\Gamma_{\rho\sigma\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\right)\\\
&\hskip
40.0pt\times(k_{n})_{\lambda_{2}}\left(\Gamma_{\lambda\alpha\beta}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}(p_{2b+1})_{\tau}\\\
&+\left[\mathcal{C}_{4+b,n-2}^{\mu\alpha\nu\beta\rho\sigma\lambda\tau\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!2}}+\cdots\right](k_{n-1})_{\lambda_{1}}(k_{n})_{\lambda_{2}}(p_{2b+1})_{\alpha}(p_{2b+1})_{\beta}\\\
&\hskip
40.0pt\times\left[\left(\Gamma_{\rho\lambda\mu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\sigma\tau\nu}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}-\left(\Gamma_{\rho\lambda\tau}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\right]\\\
&-\left[\mathcal{C}_{5+b,n-2}^{\mu\alpha\nu\beta\rho\sigma\lambda\tau\omega\epsilon\alpha_{1}\cdots\beta_{b},\rho_{1}\cdots\sigma_{n\\!-\\!2}}+\cdots\right](k_{n-2})_{\lambda_{3}}(k_{n-1})_{\lambda_{1}}(k_{n})_{\lambda_{2}}\left(\Gamma_{\omega\alpha\beta}\right)^{\lambda_{3}\rho_{n\\!-\\!2}\sigma_{n\\!-\\!2}}\\\
&\hskip
40.0pt\times\left[\left(\Gamma_{\rho\lambda\mu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\sigma\tau\nu}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}-\left(\Gamma_{\rho\lambda\tau}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\right](p_{2b+1})_{\epsilon}\Bigg{]}\\\
&+i\,\cfrac{1}{2}\,\kappa^{n}\,\Psi_{(a,b)}(-1)^{b+2}\,b\,(p_{1})_{\alpha_{1}}(p_{2})_{\beta_{1}}\cdots(p_{2b-3})_{\alpha_{b-1}}(p_{2b-2})_{\beta_{b-1}}\Bigg{[}\\\
&-\cfrac{1}{3}\left[\mathcal{C}_{2+b,n}^{\mu\nu\alpha\beta\rho\sigma\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n}}+\cdots\right](p_{2b-1})_{\mu}(p_{2b-1})_{\nu}(p_{2b})_{\alpha}(p_{2b})_{\beta}(p_{2b+1})_{\rho}(p_{2b+1})_{\sigma}\\\
&+\,\left[\mathcal{C}_{3+b,n-1}^{\mu\nu\alpha\beta\rho\sigma\omega\tau\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!1}}+\cdots\right](k_{n})_{\lambda}\left(\Gamma_{\omega\mu\nu}\right)^{\lambda\rho_{n}\sigma_{n}}(p_{2b-1})_{\tau}(p_{2b})_{\alpha}(p_{2b})_{\beta}(p_{2b+1})_{\rho}(p_{2b+1})_{\sigma}\\\
&-\left[\mathcal{C}_{4+b,n-2}^{\mu\nu\alpha\beta\rho\sigma\omega_{1}\tau_{1}\omega_{2}\tau_{2}\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!2}}+\cdots\right](k_{n\\!-\\!1})_{\lambda_{1}}\left(\Gamma_{\omega_{1}\mu\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}(k_{n})_{\lambda_{2}}\left(\Gamma_{\omega_{2}\alpha\beta}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\\\
&\hskip
40.0pt\times(p_{2b-1})_{\tau_{1}}(p_{2b})_{\tau_{2}}(p_{2b+1})_{\rho}(p_{2b+1})_{\sigma}\\\
&+\cfrac{1}{3}\left[\mathcal{C}_{5+b,n-3}^{\mu\nu\alpha\beta\rho\sigma\omega_{1}\tau_{1}\omega_{2}\tau_{2}\omega_{3}\tau_{3}\alpha_{1}\cdots\beta_{b-1},\rho_{1}\cdots\sigma_{n\\!-\\!3}}+\cdots\right]\left(\Gamma_{\omega_{1}\mu\nu}\right)^{\lambda_{1}\rho_{n\\!-\\!2}\sigma_{n\\!-\\!2}}\left(\Gamma_{\omega_{2}\alpha\beta}\right)^{\lambda_{2}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\omega_{3}\rho\sigma}\right)^{\lambda_{3}\rho_{n}\sigma_{n}}\\\
&\hskip
40.0pt\times(k_{n-2})_{\lambda_{1}}(k_{n\\!-\\!1})_{\lambda_{2}}(k_{n})_{\lambda_{3}}(p_{2b-1})_{\tau_{1}}(p_{2b})_{\tau_{2}}(p_{2b+1})_{\tau_{3}}\Bigg{]}+\text{permutations}.\end{split}$
(49)

As noted above, different terms in this expression contribute to different
orders. We include all the terms in the expression for the sake of generality.
The permutation terms account for all permutations of all $a+2\,b+1$ scalar
fields and $n$ gravitons.

###  3.5 Gauss-Bonnet term

Concluding the discussion of Horndeski gravity, we shall touch upon the scalar
field coupling to the Gauss-Bonnet term for two reasons. First, such a
coupling belongs to Horndeski gravity, but its parametrisation is complicated
and involves $\ln$ functions in both $G_{4}$ and $G_{5}$. Consequently, it is
essential to obtain the corresponding interaction rule separately. Secondly,
the coupling provides interesting phenomenology. In particular, it allows a
theory to have black holes with scalar hair, which are of utmost theoretical
interest. Therefore, the coupling to the Gauss-Bonnet term shall be discussed
separately.

It is useful to make the structure of the Gauss-Bonnet term explicit.
Sometimes, defining the term via the generalised Kronecker symbol is useful.
For our purposes, on the contrary, it is more practical to make the structure
of Riemann tensor indices contractions explicit:

$\displaystyle\begin{split}\sqrt{-g}\,\mathcal{G}&=\sqrt{-g}\left[R^{2}-4R_{\mu\nu}^{2}+R_{\mu\nu\alpha\beta}^{2}\right]\\\
&=\sqrt{-g}\left[g^{\mu\alpha}g^{\nu\beta}g^{\sigma\tau}g^{\rho\lambda}-4g^{\mu\alpha}g^{\nu\sigma}g^{\beta\tau}g^{\rho\lambda}+g^{\mu\rho}g^{\nu\sigma}g^{\alpha\lambda}g^{\beta\tau}\right]R_{\mu\nu\alpha\beta}R_{\rho\sigma\lambda\tau}.\end{split}$
(50)

We introduce the $\mathfrak{T}$ tensor that describes the contraction
symmetrically:

$\displaystyle\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma\lambda\tau}=\operatorname{Symm}\left[g^{\mu\alpha}g^{\nu\beta}g^{\sigma\tau}g^{\rho\lambda}-4g^{\mu\alpha}g^{\nu\sigma}g^{\beta\tau}g^{\rho\lambda}+g^{\mu\rho}g^{\nu\sigma}g^{\alpha\lambda}g^{\beta\tau}\right].$
(51)

The symmetrisation is performed in such a way that the tensor has the
following properties:

$\displaystyle\begin{split}&\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma\lambda\tau}=-\mathfrak{T}^{\nu\mu\alpha\beta\rho\sigma\lambda\tau}=-\mathfrak{T}^{\mu\nu\beta\alpha\rho\sigma\lambda\tau}=-\mathfrak{T}^{\mu\nu\alpha\beta\sigma\rho\lambda\tau}=-\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma\tau\lambda},\\\
&\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma\lambda\tau}=\mathfrak{T}^{\alpha\beta\mu\nu\rho\sigma\lambda\tau}=\mathfrak{T}^{\mu\nu\alpha\beta\lambda\tau\rho\sigma}=\mathfrak{T}^{\rho\sigma\lambda\tau\mu\nu\alpha\beta}.\end{split}$
(52)

Introducing the $\mathfrak{T}$ tensor is useful because it allows us to use
Riemann tensor symmetries and obtain a more compact expression.

The following formula gives the Riemann tensor with all low indices:

$\displaystyle
R_{\mu\nu\alpha\beta}=\partial_{\mu}\Gamma_{\alpha\nu\beta}-\partial_{\nu}\Gamma_{\alpha\mu\beta}+g^{\rho\sigma}\left[\Gamma_{\rho\nu\alpha}\Gamma_{\sigma\mu\beta}-\Gamma_{\rho\mu\alpha}\Gamma_{\sigma\nu\beta}\right].$
(53)

It results in the following expression for this term:

$\displaystyle\int
d^{4}x\sqrt{-g}\,\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma\lambda\tau}4\left[\partial_{\mu}\Gamma_{\alpha\nu\beta}\partial_{\rho}\Gamma_{\lambda\sigma\tau}+2g^{\omega\epsilon}\partial_{\mu}\Gamma_{\alpha\nu\beta}\Gamma_{\omega\sigma\lambda}\Gamma_{\epsilon\rho\tau}+g^{\omega_{1}\epsilon_{1}}g^{\omega_{2}\epsilon_{2}}\Gamma_{\omega_{1}\nu\alpha}\Gamma_{\epsilon_{1}\mu\beta}\Gamma_{\omega_{2}\sigma\lambda}\Gamma_{\epsilon_{2}\rho\tau}\right].$
(54)

To describe a coupling of a scalar field to the Gauss-Bonnet term, we shall
introduce a coupling function in full analogy with the other Horndeski
interactions:

$\displaystyle\int d^{4}x\sqrt{-g}\,\mathcal{G}\,f(\phi).$ (55)

In contrast with the previous cases, the Gauss-Bonnet term can only healthy
couple the scalar field itself, but not to its kinetic term. Without the loss
of generality, we assume that the coupling function is smooth enough to admit
the power series expansion, so we replace it with a single term of this
series:

$\displaystyle\int d^{4}x\sqrt{-g}\,\mathcal{G}\,f(\phi)\to\int
d^{4}x\sqrt{-g}\,\mathcal{G}\,\mathbf{g}\,\phi^{q}.$ (56)

Here $q$ is a positive integer, and $\mathbf{g}$ is a dimensional coupling
with the mass dimension $-q$.

Consequently, a coupling of the Gauss-Bonnet term to scalar fields has the
following perturbative structure:

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\,\mathcal{G}\,\mathbf{g}\,\phi^{q}\\\
&=\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{q}\cfrac{d^{4}p}{(2\pi)^{4}}\,\phi(p_{j})(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)\\\ &\hskip
8.0pt\times\kappa^{n}\mathbf{g}\,4\,\left(\sqrt{-g}\,\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma}\right)^{\rho_{1}\cdots\sigma_{n\\!-\\!2}}(k_{n\\!-1\\!})_{\lambda_{1}}(k_{n})_{\lambda_{2}}(k_{n\\!-1\\!})_{\mu}(k_{n})_{\rho}\left(\Gamma_{\alpha\nu\beta}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\lambda\sigma\tau}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\\\
&+\sum\limits_{n=3}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{q}\cfrac{d^{4}p}{(2\pi)^{4}}\,\phi(p_{j})(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)\\\ &\hskip
8.0pt\times\kappa^{n}\mathbf{g}\,8\,\left(\sqrt{-g}\,\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma}g^{\omega\epsilon}\right)^{\rho_{1}\cdots\sigma_{n\\!-\\!3}}(k_{n\\!-\\!2})_{\lambda_{1}}(k_{n\\!-\\!1})_{\lambda_{2}}(k_{n})_{\lambda_{3}}(k_{n\\!-\\!2})_{\mu}\\\
&\hskip
16.0pt\times\left(\Gamma_{\alpha\nu\beta}\right)^{\lambda_{1}\rho_{n\\!-\\!2}\sigma_{n\\!-\\!2}}\left(\Gamma_{\omega\sigma\lambda}\right)^{\lambda_{2}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\epsilon\rho\tau}\right)^{\lambda_{3}\rho_{n}\sigma_{n}}\\\
&+\sum\limits_{n=4}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\prod\limits_{j=1}^{q}\cfrac{d^{4}p}{(2\pi)^{4}}\,\phi(p_{j})(2\pi)^{4}\delta\left(\sum
k_{i}+\sum p_{j}\right)\\\ &\hskip
8.0pt\times\kappa^{n}\mathbf{g}\,4\,\left(\sqrt{-g}\,\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma}g^{\omega_{1}\epsilon_{1}}g^{\omega_{2}\epsilon_{2}}\right)^{\rho_{1}\cdots\sigma_{n\\!-\\!4}}(k_{n\\!-\\!3})_{\lambda_{1}}(k_{n\\!-\\!2})_{\lambda_{2}}(k_{n\\!-\\!1})_{\lambda_{3}}(k_{n})_{\lambda_{4}}\\\
&\hskip
16.0pt\times\left(\Gamma_{\omega_{1}\nu\alpha}\right)^{\lambda_{1}\rho_{n\\!-\\!3}\sigma_{n\\!-\\!3}}\left(\Gamma_{\epsilon_{1}\mu\beta}\right)^{\lambda_{2}\rho_{n\\!-\\!2}\sigma_{n\\!-\\!2}}\left(\Gamma_{\omega_{2}\sigma\lambda}\right)^{\lambda_{3}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\epsilon_{2}\rho\tau}\right)^{\lambda_{4}\rho_{n}\sigma_{n}}.\end{split}$
(57)

In full accordance with the previous results, this interaction does not
contribute at $\kappa^{0}$ and $\kappa^{1}$. The first term of this expression
contributes at $\kappa^{2}$, the second at $\kappa^{3}$, and the last at
$\kappa^{4}$.

This expression results in the following interaction rule:

$\displaystyle\begin{split}\\\ &\hskip
35.0pt\begin{gathered}\fmffile{Diags/GB_{v}ertex}\fmfgraph*(30,30)\fmfleft{L1,L2}\fmfright{R1,R2}\fmf{dbl_{w}iggly}{L1,V}\fmf{dbl_{w}iggly}{L2,V}\fmf{dashes}{V,R1}\fmf{dashes}{V,R2}\fmfdot{V}\fmffreeze\fmf{dots}{L1,L2}\fmf{dots}{R1,R2}\fmflabel{$\rho_1\sigma_1,k_1$}{L1}\fmflabel{$\rho_n\sigma_n,k_n$}{L2}\fmflabel{$p_1$}{R1}\fmflabel{$p_q$}{R2}\fmflabel{$\mathbf{g}$}{V}\fmfv{$g$,label.angle=90,label.dist=15pt}{V}\end{gathered}=i\,\kappa^{n}4\mathbf{g}\Bigg{[}\\!\\!\left(\sqrt{-g}\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma}\right)^{\rho_{1}\cdots\sigma_{n\\!-\\!2}}\\!(k_{n\\!-1\\!})_{\lambda_{1}}\\!(k_{n})_{\lambda_{2}}(k_{n\\!-1\\!})_{\mu}(k_{n})_{\rho}\left(\Gamma_{\alpha\nu\beta}\right)^{\lambda_{1}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\\!\left(\Gamma_{\lambda\sigma\tau}\right)^{\lambda_{2}\rho_{n}\sigma_{n}}\\\
\\\
&+\\!2\left(\sqrt{-g}\,\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma}g^{\omega\epsilon}\right)^{\rho_{1}\cdots\sigma_{n\\!-\\!3}}\\!(k_{n\\!-\\!2})_{\lambda_{1}}\\!(k_{n\\!-\\!1})_{\lambda_{2}}\\!(k_{n})_{\lambda_{3}}(k_{n\\!-\\!2})_{\mu}\left(\Gamma_{\alpha\nu\beta}\right)^{\lambda_{1}\rho_{n\\!-\\!2}\sigma_{n\\!-\\!2}}\left(\Gamma_{\omega\sigma\lambda}\right)^{\lambda_{2}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\epsilon\rho\tau}\right)^{\lambda_{3}\rho_{n}\sigma_{n}}\\\
&+\left(\sqrt{-g}\,\mathfrak{T}^{\mu\nu\alpha\beta\rho\sigma}g^{\omega_{1}\epsilon_{1}}g^{\omega_{2}\epsilon_{2}}\right)^{\rho_{1}\cdots\sigma_{n\\!-\\!4}}(k_{n\\!-\\!3})_{\lambda_{1}}(k_{n\\!-\\!2})_{\lambda_{2}}(k_{n\\!-\\!1})_{\lambda_{3}}(k_{n})_{\lambda_{4}}\\\
&\hskip
20.0pt\times\left(\Gamma_{\omega_{1}\nu\alpha}\right)^{\lambda_{1}\rho_{n\\!-\\!3}\sigma_{n\\!-\\!3}}\left(\Gamma_{\epsilon_{1}\mu\beta}\right)^{\lambda_{2}\rho_{n\\!-\\!2}\sigma_{n\\!-\\!2}}\left(\Gamma_{\omega_{2}\sigma\lambda}\right)^{\lambda_{3}\rho_{n\\!-\\!1}\sigma_{n\\!-\\!1}}\left(\Gamma_{\epsilon_{2}\rho\tau}\right)^{\lambda_{4}\rho_{n}\sigma_{n}}\Bigg{]}+\text{permutations}.\end{split}$
(58)

Similarly to the previous cases, the permutation term accounts for all
permutations of all $n$ gravitons and $q$ scalars.

##  4 Propagator of massive gravity

Massive gravity is a direction of research that aims to create a consistent
theory of spin $2$ massive fundamental particles. The detailed review of
massive gravity lies far beyond the scope of this paper and can be found in
[25, 26]. Here, we only highlight the most fundamental challenges one faces in
an attempt to implement it for FeynGrav.

The kinetic term describing the propagation of a spin $2$ particle with a non-
vanishing mass is fixed uniquely and known as the Fierz-Pauli Lagrangian:

$\displaystyle\int
d^{D}x\left[-\cfrac{1}{2}\,h^{\mu\nu}\mathcal{O}_{\mu\nu\alpha\beta}h^{\alpha\beta}\right].$
(59)

Here $\mathcal{O}_{\mu\nu\alpha\beta}$ is an operator defiened as follows:

$\displaystyle\mathcal{O}_{\mu\nu\alpha\beta}=\left[\cfrac{1}{2}\left(\eta_{\mu\alpha}\eta_{\nu\beta}+\eta_{\mu\beta}\eta_{\nu\alpha}\right)-\eta_{\mu\nu}\eta_{\alpha\beta}\right](\square+m^{2})-\partial_{\mu}\partial_{\alpha}\eta_{\nu\beta}-\partial_{\mu}\partial_{\beta}\eta_{\nu\alpha}+\partial_{\mu}\partial_{\nu}\eta_{\alpha\beta}+\partial_{\alpha}\partial_{\beta}\eta_{\mu\nu}.$
(60)

The part of the Fierz-Pauli action that does not depend on the mass is taken
directly from general relativity, where the gauge symmetry defines it
uniquely. The mass term is not fixed uniquely from the first principles since
there are two invariants quadratic in perturbations and free from derivatives:
$h_{\mu\nu}^{2}$ and $(\eta^{\mu\nu}h_{\mu\nu})^{2}$. Since the theory does
not admit a gauge symmetry, the relative coefficient between these terms
remains free. However, the stability reasoning fixed it uniquely. Unless the
mass term takes the Fierz-Pauli form, the action will contain a ghost degree
of freedom.

The massive graviton propagator, consequently, is also unique:

$\displaystyle\mathcal{G}_{\mu\nu\alpha\beta}(p)=\cfrac{i}{p^{2}-m^{2}}\left[\cfrac{1}{2}\left(P_{\mu\alpha}P_{\nu\beta}+P_{\mu\beta}P_{\nu\alpha}\right)-\cfrac{1}{D-1}P_{\mu\nu}P_{\alpha\beta}\right].$
(61)

Here, operators $P_{\mu\nu}$ are defined as follows:

$\displaystyle P_{\mu\nu}(p)=\eta_{\mu\nu}-\cfrac{p_{\mu}p_{\nu}}{m^{2}}.$
(62)

Due to the absence of symmetry, stability is the only guiding principle that
defines the interaction sector of the theory. Massive gravity models are
usually plagued by the Boulware-Deser ghost [46], so the absence of the ghost
is the only guiding principle.

To the best of our knowledge, the two most widely recognised consistent
massive gravity models are the Dvali-Gabadadze-Porrati [47, 48, 49] and de
Rham-Gabadadze-Tolley [50]. The Dvali-Gabadadze-Porrati describes massive
gravity in terms of a five-dimensional braneworld model with infinite extra
dimensions. At the same time, The de Rham-Gabadadze-Tolley extensively relies
on the vierbein formulation of gravity and bi-gravity.

Both models use formalisms that cannot easily be converted to standard four-
dimensional geometric descriptions. Without this, there is no direct way to
implement them in FeynGrav. Further publications will address this challenge
since it deserves detailed treatment.

The graviton propagator presents a notable exception since, as noted above, it
is defined uniquely. Therefore, the new version of FeynGrav implements it
directly to provide additional tools to study quantum gravity.

##  5 Axion-like coupling to a single vector field

Quantum chromodynamics and $SU(N)$ Yang-Mills theories generally point towards
the possible existence of axions. The physics and phenomenology of axions are
vast and lie beyond the scope of this paper. Additionally, since there are
many models of axions, their implementation in FeynGrav requires a mode of
detailed analysis. Such an analysis, together with the most perspective
models, will be discussed in the forthcoming publications. This section only
discusses the simplest axion-like coupling between a scalar and a single
vector field. The reason for this is twofold. Firstly, this coupling is simple
enough to be implemented immediately in FeynGrav. Secondly, such a coupling is
used in many models involving scalar fields, providing an interesting
phenomenology.

Within the $SU(N)$ Yang-Mills model in a flat spacetime, the following term is
a complete derivative:

$\displaystyle\int
d^{4}x\,\cfrac{1}{2}\,\theta\,\operatorname{tr}\left[G_{\mu\nu}\,\widetilde{G}^{\mu\nu}\right].$
(63)

Here, $\theta$ is an arbitrary constant, $G_{\mu\nu}$ is the Yang-Mills field
tensor, and $\widetilde{G}^{\mu\nu}$ is the Hodge dual tensor:

$\displaystyle\widetilde{G}^{\mu\nu}\overset{\text{def}}{=}\cfrac{1}{2}\,\varepsilon^{\mu\nu\alpha\beta}G_{\alpha\beta}.$
(64)

Here, $\varepsilon^{\mu\nu\alpha\beta}$ is the Levi-Chivita tensor. Since it
is a complete derivative, it does not contribute to the classical field
equations and interaction rules. However, it is of enormous importance for the
vacuum structure of the theory [28, 29, 51, 52, 53].

The theory cannot predict the value of the constant $\theta$. Experiments, in
turn, put an extreme constraint on the coupling [54, 55, 56]:

$\displaystyle\lvert\theta\rvert<10^{-10}.$ (65)

One way to explain such a small value is to promote $\theta$ to a dynamical
field $\theta(x)$, called an axion, with the non-minimal coupling to the gauge
field:

$\displaystyle\int
d^{4}x\,\cfrac{1}{2}\,\theta(x)\,\operatorname{tr}\left[G_{\mu\nu}\,\widetilde{G}^{\mu\nu}\right].$
(66)

Models of axions within quantum chromodynamics led to many other models
involving similar axion-like couplings. This paper only discusses the coupling
between a scalar field and a single vector field. Despite its apparent
simplicity, such a coupling has a broad phenomenology [33, 57, 58].

We adopt the following parametrisation of the axion-like coupling between a
scalar field $\phi$ and a single vector field in curved spacetime:

$\displaystyle\int
d^{4}x\sqrt{-g}\left[\cfrac{1}{2}\,F_{\mu\nu}\,\widetilde{F}^{\mu\nu}\,f(\phi)\right]$
(67)

Here $f(\phi)$ is a infinitely differentiable function of the scalar field
$\phi$; $F_{\mu\nu}$ is the standard field strength tensor:

$\displaystyle
F_{\mu\nu}=\nabla_{\mu}A_{\nu}-\nabla_{\nu}A_{\mu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu},$
(68)

and $\widetilde{F}^{\mu\nu}$ is the dual tensor:

$\displaystyle\widetilde{F}^{\mu\nu}\overset{\text{def}}{=}\varepsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}.$
(69)

In full similarity with the other models, one can replace the function $f$
with a single power-like term without the loss of generality. The interaction
terms take the following simple form:

$\displaystyle\int
d^{4}x\sqrt{-g}\,\cfrac{1}{4}\,\varepsilon^{\mu\nu\alpha\beta}\,F_{\mu\nu}\,F_{\alpha\beta}\,\cfrac{\lambda}{l!}\,\phi^{l}.$
(70)

Here, $\lambda$ is a coupling with the mass dimension $-l$.

The perturbative structure of the interaction is obtained similarly to the
other models discussed in this article:

$\displaystyle\begin{split}&\int
d^{4}x\sqrt{-g}\,\cfrac{1}{4}\,\varepsilon^{\mu\nu\alpha\beta}\,F_{\mu\nu}\,F_{\alpha\beta}\,\cfrac{\lambda}{l!}\,\phi^{l}\\\
&=\sum\limits_{n=0}^{\infty}\int\prod\limits_{i=1}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}h_{\rho_{i}\sigma_{i}}(k_{i})\cfrac{d^{4}q_{1}}{(2\pi)^{4}}A_{\lambda_{1}}(q_{1})\cfrac{d^{4}q_{2}}{(2\pi)^{4}}A_{\lambda_{2}}(q_{2})\prod\limits_{j=1}^{l}\cfrac{d^{4}p_{j}}{(2\pi)^{4}}\phi(p_{j})\,(2\pi)^{4}\delta\left(\sum
k_{i}+q_{1}+q_{2}+\sum p_{j}\right)\\\
&\times(-1)\kappa^{n}\,\cfrac{\lambda}{l!}\,\left(\sqrt{-g}\right)^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}\varepsilon^{\tau_{1}\lambda_{1}\tau_{2}\lambda_{2}}(q_{1})_{\tau_{1}}(q_{2})_{\tau_{2}}.\end{split}$
(71)

The corresponding interaction rule reads:

$\displaystyle\begin{gathered}\fmffile{Diags/Axion_{V}ector_{v}ertex}\fmfgraph*(40,40)\fmfleft{L1,L2,L3}\fmftop{T}\fmfbottom{B}\fmfright{R0,R1,R2,R3,R4}\fmf{dbl_{w}iggly}{L1,V}\fmf{dbl_{w}iggly}{L3,V}\fmf{dashes}{R1,V}\fmf{dashes}{R3,V}\fmflabel{$\lambda$}{V}\fmfv{$\lambda$,label.angle=90,label.dist=15pt}{V}\fmffreeze\fmfdot{V}\fmflabel{$\rho_1\sigma_1$}{L1}\fmflabel{$\rho_n\sigma_n$}{L3}\fmf{dots}{L1,L3}\fmflabel{$p_1$}{R1}\fmflabel{$p_l$}{R3}\fmf{dots}{R1,R3}\fmf{photon}{V,R0}\fmf{photon}{V,R4}\fmflabel{$\lambda_1,q_1$}{R0}\fmflabel{$\lambda_2,q_2$}{R4}\end{gathered}\hskip
30.0pt=-i\,\kappa^{n}\,\lambda\,\left(\sqrt{-g}\right)^{\rho_{1}\sigma_{1}\cdots\rho_{n}\sigma_{n}}\varepsilon^{\tau_{1}\lambda_{1}\tau_{2}\lambda_{2}}(q_{1})_{\tau_{1}}(q_{2})_{\tau_{2}}.$
(72)

In this expression, the “permutation” term is absent since the expression is
already symmetrical. Namely, in this parametrisation, the expression is valid
for coupling to any number of scalars. It is already symmetric with respect to
permutations of all graviton and vector lines.

##  6 Quadratic gravity

Quadratic gravity is a generalisation of general relativity that includes the
terms quadratic in curvature:

$\displaystyle\mathcal{A}=\int
d^{4}x\sqrt{-g}\left(-\cfrac{2}{\kappa^{2}}\right)\left[R-\cfrac{1}{6\,m_{0}^{2}}\,R^{2}+\cfrac{1}{m_{2}^{2}}\left(R_{\mu\nu}^{2}-\cfrac{1}{3}\,R^{2}\right)\right].$
(73)

In this expression, $\kappa$ is the same gravitational coupling used
throughout this paper, and $m_{0}$ and $m_{2}$ are two parameters of the
theory with the dimension of mass.

The model was extensively studied in the literature throughout the years [34,
35, 36, 37, 38, 39, 40, 59, 60, 60]. The main aim of this paper is to discuss
the Feynman rules obtained within the standard approach, so we do not discuss
the theory in great details. Below, we briefly recall the essential features
of the model.

First and foremost, the model is the most general extension of general
relativity with the dimension-four operators in the following sense. There are
only four dimension-four operators that include only metric and its
derivatives: $R^{2}$, $R_{\mu\nu}^{2}$, $R_{\mu\nu\alpha\beta}^{2}$, and
$\square R$. Operator $\square R$ is a complete derivative that can only
influence the surface term. The Riemann tensor square operator
$R_{\mu\nu\alpha\beta}^{2}$ is not a complete derivative, but can always be
promoted to the Gauss-Bonnet term
$R_{\mu\nu\alpha\beta}^{2}-4R_{\mu\nu}^{2}+R^{2}$. Consequently, only two
dimension-four operators can meaningfully influence the model and are both
present in quadratic gravity.

Secondly, the model contains new degrees of freedom: a massive scalar and a
massive spin-$2$ ghost. Both perturbative and non-perturbative results show
the existence of new degrees of freedom. The new higher-order operators
contain additional derivatives, which influence the propagator. The propagator
receives new poles, which are associated with the new degrees of freedom [12].
At the same time, there is a nonlinear map between the original
parametrisation of the model and the parametrisation, making the presence of
new degrees of freedom explicit [61]. The discussion on the physical role of
these poles is extensive and lies far beyond the scope of this paper.
Consequently, we do not discuss the physical meaning of these poles and their
interpretation but pursue the main aim of this paper, which is the derivation
of the interaction rules of the model.

Lastly, the model admits a gauge symmetry. In complete analogy with general
relativity, the model describes massless spin-$2$ states which require gauge
symmetry. Therefore, one shall use a gauge-fixing procedure to obtain the
propagator of the model. Since the new degrees of freedom are massive, they do
not admit a gauge symmetry and do not require a gauge fixing. Consequently,
one can use the same gauge fixing method employed in general relativity. In
the context of FeynGrav, the Faddeev-Popov prescription’s gauge fixing was
discussed in [3].

Let us proceed with the derivation of the interaction rules. To obtain the
Feynman rules, one shall only describe the perturbative structure of two new
operators. The following formulae hold for the Riemann tensor, Ricci tensor,
and the scalar curvature:

$\displaystyle\begin{split}R_{\mu\nu\alpha\beta}&=\partial_{\mu}\Gamma_{\alpha\nu\beta}-\partial_{\nu}\Gamma_{\alpha\mu\beta}+g^{\rho\sigma}\left[\Gamma_{\rho\nu\alpha}\Gamma_{\sigma\mu\beta}-\Gamma_{\rho\mu\alpha}\Gamma_{\sigma\nu\beta}\right],\\\
R_{\mu\nu}&=g^{\alpha\beta}R_{\alpha\mu\beta\nu}=g^{\alpha\beta}\left[\partial_{\alpha}\Gamma_{\beta\mu\nu}-\partial_{\mu}\Gamma_{\alpha\beta\nu}\right]+g^{\alpha\beta}g^{\rho\sigma}\left[\Gamma_{\rho\alpha\mu}\Gamma_{\sigma\beta\nu}-\Gamma_{\rho\alpha\beta}\Gamma_{\sigma\mu\nu}\right],\\\
R&=g^{\mu\nu}R_{\mu\nu}=g^{\mu\nu}g^{\alpha\beta}\partial_{\mu}\left[\Gamma_{\nu\alpha\beta}-\Gamma_{\alpha\beta\nu}\right]+g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}\left[\Gamma_{\mu\alpha\rho}\Gamma_{\nu\beta\sigma}-\Gamma_{\mu\alpha\beta}\Gamma_{\nu\rho\sigma}\right].\end{split}$
(74)

Consequently, the following formulae describe squares of these quantities:

$\displaystyle\begin{split}R^{2}=&g^{\mu\nu}g^{\alpha\beta}g^{mn}g^{ab}\partial_{\mu}\left\\{\Gamma_{\nu\alpha\beta}-\Gamma_{\alpha\beta\nu}\right\\}\partial_{m}\left\\{\Gamma_{nab}-\Gamma_{abn}\right\\}\\\
&+2\,g^{\mu\nu}g^{\alpha\beta}g^{mn}g^{ab}g^{rs}\partial_{\mu}\left\\{\Gamma_{\nu\alpha\beta}-\Gamma_{\alpha\beta\nu}\right\\}\left\\{\Gamma_{mar}\Gamma_{nbs}-\Gamma_{mab}\Gamma_{nrs}\right\\}\\\
&+g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}g^{mn}g^{ab}g^{rs}\left\\{\Gamma_{\mu\alpha\rho}\Gamma_{\nu\beta\sigma}-\Gamma_{\mu\alpha\beta}\Gamma_{\nu\rho\sigma}\right\\}\left\\{\Gamma_{mar}\Gamma_{nbs}-\Gamma_{mab}\Gamma_{nrs}\right\\},\end{split}$
(75)
$\displaystyle\begin{split}R_{\mu\nu}^{2}=g^{\mu\alpha}g^{\mu\beta}\Bigg{[}&g^{\rho\sigma}g^{rs}\left\\{\partial_{\rho}\Gamma_{\sigma\mu\nu}-\partial_{\mu}\Gamma_{\rho\sigma\nu}\right\\}\left\\{\partial_{r}\Gamma_{s\alpha\beta}-\partial_{\alpha}\Gamma_{rs\beta}\right\\}\\\
&+2\,g^{\rho\sigma}g^{rs}g^{lt}\left\\{\partial_{\rho}\Gamma_{\sigma\mu\nu}-\partial_{\mu}\Gamma_{\rho\sigma\nu}\right\\}\left\\{\Gamma_{lr\alpha}\Gamma_{ts\beta}-\Gamma_{lrs}\Gamma_{t\alpha\beta}\right\\}\\\
&+g^{\rho\sigma}g^{\lambda\tau}g^{rs}g^{lt}\left\\{\Gamma_{\lambda\rho\mu}\Gamma_{\tau\sigma\nu}-\Gamma_{\lambda\rho\sigma}\Gamma_{\tau\mu\nu}\right\\}\left\\{\Gamma_{lr\alpha}\Gamma_{ts\beta}-\Gamma_{lrs}\Gamma_{t\alpha\beta}\right\\}\Bigg{]},\end{split}$
(76) $\displaystyle\begin{split}R_{\mu\nu\alpha\beta}^{2}=g^{\mu m}g^{\nu
n}g^{\alpha a}g^{\beta
b}\Bigg{[}&\left\\{\partial_{\mu}\Gamma_{\alpha\nu\beta}-\partial_{\nu}\Gamma_{\alpha\mu\beta}\right\\}\left\\{\partial_{m}\Gamma_{anb}-\partial_{n}\Gamma_{amb}\right\\}\\\
&+2g^{rs}\left\\{\partial_{\mu}\Gamma_{\alpha\nu\beta}-\partial_{\nu}\Gamma_{\alpha\mu\beta}\right\\}\left\\{\Gamma_{rna}\Gamma_{smb}-\Gamma_{rma}\Gamma_{snb}\right\\}\\\
&+g^{\rho\sigma}g^{rs}\left\\{\Gamma_{\rho\nu\alpha}\Gamma_{\sigma\mu\beta}-\Gamma_{\rho\mu\alpha}\Gamma_{\sigma\nu\beta}\right\\}\left\\{\Gamma_{\rho\nu\alpha}\Gamma_{\sigma\mu\beta}-\Gamma_{\rho\mu\alpha}\Gamma_{\sigma\nu\beta}\right\\}\Bigg{]}.\end{split}$
(77)

The perturbative structure of $R^{2}$ and $R_{\mu\nu}^{2}$ terms reads:

$\displaystyle\begin{split}\int d^{4}x\sqrt{-g}R^{2}=&\int
d^{4}x\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{mn}g^{ab}\partial_{\mu}\left\\{\Gamma_{\nu\alpha\beta}-\Gamma_{\alpha\beta\nu}\right\\}\partial_{m}\left\\{\Gamma_{nab}-\Gamma_{abn}\right\\}\\\
&+2\,\int
d^{4}x\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{mn}g^{ab}g^{rs}\partial_{\mu}\left\\{\Gamma_{\nu\alpha\beta}-\Gamma_{\alpha\beta\nu}\right\\}\left\\{\Gamma_{mar}\Gamma_{nbs}-\Gamma_{mab}\Gamma_{nrs}\right\\}\\\
&+\int
d^{4}x\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}g^{mn}g^{ab}g^{rs}\left\\{\Gamma_{\mu\alpha\rho}\Gamma_{\nu\beta\sigma}-\Gamma_{\mu\alpha\beta}\Gamma_{\nu\rho\sigma}\right\\}\left\\{\Gamma_{mar}\Gamma_{nbs}-\Gamma_{mab}\Gamma_{nrs}\right\\}\end{split}$
(78) $\displaystyle=$
$\displaystyle\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=2}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}h_{\rho_{i}\sigma_{i}}(k_{i})(2\pi)^{4}\delta\left(\sum
k_{i}\right)\,\kappa^{n}\,\left(\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{mn}g^{ab}\right)^{\rho_{3}\sigma_{3}\cdots\rho_{n}\sigma_{n}}$
$\displaystyle\times(k_{1})_{\mu}(k_{1})_{\lambda_{1}}(k_{2})_{m}(k_{2})_{\lambda_{2}}\left[\left(\Gamma_{\nu\alpha\beta}\right)^{\lambda_{1}\rho_{1}\sigma}-\left(\Gamma_{\alpha\beta\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\right]\left[\left(\Gamma_{nab}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}-\left(\Gamma_{abn}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\right]$
$\displaystyle+$ $\displaystyle
2\,\sum\limits_{n=3}^{\infty}\int\prod\limits_{i=3}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}h_{\rho_{i}\sigma_{i}}(k_{i})(2\pi)^{4}\delta\left(\sum
k_{i}\right)\,\kappa^{n}\,\left(\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{mn}g^{ab}g^{rs}\right)^{\rho_{4}\sigma_{4}\cdots\rho_{n}\sigma_{n}}$
$\displaystyle\times(k_{1})_{\mu}\\!(k_{1})_{\lambda_{1}}\\!(k_{2})_{\lambda_{2}}\\!(k_{3})_{\lambda_{3}}\\!\\!\left[\\!\left(\Gamma_{\nu\alpha\beta}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\\!-\\!\left(\Gamma_{\alpha\beta\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\\!\right]\\!\\!\left[\\!\left(\Gamma_{mar}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\\!\\!\left(\Gamma_{nbs}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\\!-\\!\left(\Gamma_{mab}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\\!\\!\left(\Gamma_{nrs}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\\!\right]$
$\displaystyle+$
$\displaystyle\sum\limits_{n=4}^{\infty}\int\prod\limits_{i=4}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}h_{\rho_{i}\sigma_{i}}(k_{i})(2\pi)^{4}\delta\left(\sum
k_{i}\right)\,\kappa^{n}\,\left(\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}g^{mn}g^{ab}g^{rs}\right)^{\rho_{5}\sigma_{5}\cdots\rho_{n}\sigma_{n}}$
$\displaystyle\times(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}(k_{3})_{\lambda_{3}}(k_{4})_{\lambda_{4}}\left[\left(\Gamma_{\mu\alpha\rho}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\left(\Gamma_{\nu\beta\sigma}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}-\left(\Gamma_{\mu\alpha\beta}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\left(\Gamma_{\nu\rho\sigma}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\right]$
$\displaystyle\hskip
10.0pt\times\left[\left(\Gamma_{mar}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\left(\Gamma_{nbs}\right)^{\lambda_{4}\rho_{4}\sigma_{4}}-\left(\Gamma_{mab}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\left(\Gamma_{nrs}\right)^{\lambda_{4}\rho_{4}\sigma_{4}}\right].$
$\displaystyle\begin{split}\int d^{4}x\sqrt{-g}R_{\mu\nu}^{2}=&\int
d^{4}x\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{rs}\left\\{\partial_{\rho}\Gamma_{\sigma\mu\nu}-\partial_{\mu}\Gamma_{\rho\sigma\nu}\right\\}\left\\{\partial_{r}\Gamma_{s\alpha\beta}-\partial_{\alpha}\Gamma_{rs\beta}\right\\}\\\
&+2\,\int
d^{4}x\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{rs}g^{lt}\left\\{\partial_{\rho}\Gamma_{\sigma\mu\nu}-\partial_{\mu}\Gamma_{\rho\sigma\nu}\right\\}\left\\{\Gamma_{lr\alpha}\Gamma_{ts\beta}-\Gamma_{lrs}\Gamma_{t\alpha\beta}\right\\}\\\
&+\int
d^{4}x\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{\lambda\tau}g^{rs}g^{lt}\left\\{\Gamma_{\lambda\rho\mu}\Gamma_{\tau\sigma\nu}-\Gamma_{\lambda\rho\sigma}\Gamma_{\tau\mu\nu}\right\\}\left\\{\Gamma_{lr\alpha}\Gamma_{ts\beta}-\Gamma_{lrs}\Gamma_{t\alpha\beta}\right\\}\Bigg{]}\end{split}$
(79) $\displaystyle=$
$\displaystyle\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=2}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}h_{\rho_{i}\sigma_{i}}(k_{i})(2\pi)^{4}\delta\left(\sum
k_{i}\right)\,\kappa^{n}\,\left(\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{rs}\right)^{\rho_{3}\sigma_{3}\cdots\rho_{n}\sigma_{n}}$
$\displaystyle\times(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}\left[(k_{1})_{\rho}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}-(k_{1})_{\mu}\left(\Gamma_{\rho\sigma\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\right]\left[(k_{2})_{r}\left(\Gamma_{s\alpha\beta}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}-(k_{2})_{\alpha}\left(\Gamma_{rs\beta}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\right]$
$\displaystyle+$ $\displaystyle
2\,\sum\limits_{n=3}^{\infty}\int\prod\limits_{i=3}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}h_{\rho_{i}\sigma_{i}}(k_{i})(2\pi)^{4}\delta\left(\sum
k_{i}\right)\,\kappa^{n}\,\left(\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{rs}g^{lt}\right)^{\rho_{4}\sigma_{4}\cdots\rho_{n}\sigma_{n}}$
$\displaystyle\times(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}(k_{3})_{\lambda_{3}}\left[(k_{1})_{\rho}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}-(k_{1})_{\mu}\left(\Gamma_{\rho\sigma\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\right]$
$\displaystyle\hskip
10.0pt\times\left[\left(\Gamma_{lr\alpha}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\left(\Gamma_{ts\beta}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}-\left(\Gamma_{lrs}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\left(\Gamma_{t\alpha\beta}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\right]$
$\displaystyle+$
$\displaystyle\sum\limits_{n=4}^{\infty}\int\prod\limits_{i=4}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}h_{\rho_{i}\sigma_{i}}(k_{i})(2\pi)^{4}\delta\left(\sum
k_{i}\right)\,\kappa^{n}\,\left(\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{\lambda\tau}g^{rs}g^{lt}\right)^{\rho_{5}\sigma_{5}\cdots\rho_{n}\sigma_{n}}$
$\displaystyle\times(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}(k_{3})_{\lambda_{3}}(k_{4})_{\lambda_{4}}\left[\left(\Gamma_{\lambda\rho\mu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\left(\Gamma_{\tau\sigma\nu}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}-\left(\Gamma_{\lambda\rho\sigma}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\left(\Gamma_{\tau\mu\nu}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\right]$
$\displaystyle\hskip
10.0pt\times\left[\left(\Gamma_{lr\alpha}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\left(\Gamma_{ts\beta}\right)^{\lambda_{4}\rho_{4}\sigma_{4}}-\left(\Gamma_{lrs}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\left(\Gamma_{t\alpha\beta}\right)^{\lambda_{4}\rho_{4}\sigma_{4}}\right].$

Further, let us return to the gauge fixing procedure. In the previous
publication [3], we discussed the Faddeev-Popov gauge fixing procedure for
general relativity in detail, so we refrain from repeating the discussion
here. To simplify the implementation of quadratic gravity within FeynGrav, we
introduce the same de Donder gauge fixing term:

$\displaystyle\mathcal{A}_{\text{gf}}=\int
d^{4}x\sqrt{-g}\,\cfrac{\epsilon}{2\,\kappa^{2}}\,g_{\mu\nu}\Big{[}g^{\alpha\beta}\Gamma^{\mu}_{\alpha\beta}\Big{]}\Big{[}g^{\rho\sigma}\Gamma^{\nu}_{\rho\sigma}\Big{]}.$
(80)

Here, $\epsilon$ is the gauge fixing parameter.

That gauge fixing parameter contributes both to the propagator and vertices.
At the same time, the Faddeev-Popov sector of the theory is not affected by
the higher-order operators. From the physical point of view, this is because
the new degrees of freedom are massive and do not admit any additional gauge
symmetry.

The perturbative structure of this term was described in [3], so we use the
expression from this publication to obtain the perturbative structure of the
model. Consequently, the quadratic gravity action has the following
perturbative structure.

$\displaystyle\mathcal{A}=\int
d^{4}x\sqrt{-g}\left(-\cfrac{2}{\kappa^{2}}\right)\left[R-\cfrac{1}{3}\left(\cfrac{1}{2\,m_{0}^{2}}+\cfrac{1}{m_{2}^{2}}\right)R+\cfrac{1}{m_{2}^{2}}\,R_{\mu\nu}^{2}+\cfrac{\epsilon}{2\,\kappa^{2}}\,g_{\mu\nu}\Big{[}g^{\alpha\beta}\Gamma^{\mu}_{\alpha\beta}\Big{]}\Big{[}g^{\rho\sigma}\Gamma^{\nu}_{\rho\sigma}\Big{]}\right]$
(81) $\displaystyle=$
$\displaystyle\sum\limits_{n=2}^{\infty}\int\prod\limits_{i=2}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\,(2\pi)^{4}\delta\left(\sum
k_{i}\right)\kappa^{n-2}$
$\displaystyle\times\Bigg{\\{}2\,\left(\sqrt{-g}g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}\right)^{\rho_{3}\sigma_{3}\cdots\rho_{n}\sigma_{n}}(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}$
$\displaystyle\hskip
10.0pt\times\left[\left(\Gamma_{\alpha\mu\rho}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\left(\Gamma_{\sigma\nu\beta}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}-\left(\Gamma_{\alpha\mu\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\left(\Gamma_{\rho\beta\sigma}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}-\cfrac{\epsilon}{4}\left(\Gamma_{\mu\alpha\beta}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\left(\Gamma_{\nu\rho\sigma}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\right]$
$\displaystyle+\cfrac{m_{2}^{2}+2\,m_{0}^{2}}{3\,m_{0}^{2}\,m_{2}^{2}}\,\left(\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{mn}g^{ab}\right)^{\rho_{3}\sigma_{3}\cdots\rho_{n}\sigma_{n}}(k_{1})_{\mu}(k_{1})_{\lambda_{1}}(k_{2})_{m}(k_{2})_{\lambda_{2}}$
$\displaystyle\hskip
10.0pt\times\left[\left(\Gamma_{\nu\alpha\beta}\right)^{\lambda_{1}\rho_{1}\sigma}-\left(\Gamma_{\alpha\beta\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\right]\left[\left(\Gamma_{nab}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}-\left(\Gamma_{abn}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\right]$
$\displaystyle-\cfrac{2}{m_{2}^{2}}\left(\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{rs}\right)^{\rho_{3}\sigma_{3}\cdots\rho_{n}\sigma_{n}}(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}$
$\displaystyle\hskip
10.0pt\times\left[(k_{1})_{\rho}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}-(k_{1})_{\mu}\left(\Gamma_{\rho\sigma\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\right]\left[(k_{2})_{r}\left(\Gamma_{s\alpha\beta}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}-(k_{2})_{\alpha}\left(\Gamma_{rs\beta}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\right]\Bigg{\\}}$
$\displaystyle+$
$\displaystyle\sum\limits_{n=3}^{\infty}\int\prod\limits_{i=3}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\,(2\pi)^{4}\delta\left(\sum
k_{i}\right)\kappa^{n-2}\,2$
$\displaystyle\times\Bigg{\\{}\cfrac{m_{2}^{2}+2\,m_{0}^{2}}{3\,m_{0}^{2}\,m_{2}^{2}}\,\left(\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{mn}g^{ab}g^{rs}\right)^{\rho_{4}\sigma_{4}\cdots\rho_{n}\sigma_{n}}(k_{1})_{\mu}(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}(k_{3})_{\lambda_{3}}$
$\displaystyle\hskip
10.0pt\times\left[\left(\Gamma_{\nu\alpha\beta}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}-\left(\Gamma_{\alpha\beta\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\right]\left[\left(\Gamma_{mar}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\left(\Gamma_{nbs}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}-\left(\Gamma_{mab}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\left(\Gamma_{nrs}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\right]$
$\displaystyle-\cfrac{2}{m_{2}^{2}}\left(\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{rs}g^{lt}\right)^{\rho_{4}\sigma_{4}\cdots\rho_{n}\sigma_{n}}(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}(k_{3})_{\lambda_{3}}$
$\displaystyle\hskip
10.0pt\times\left[(k_{1})_{\rho}\left(\Gamma_{\sigma\mu\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}-(k_{1})_{\mu}\left(\Gamma_{\rho\sigma\nu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\right]\left[\left(\Gamma_{lr\alpha}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\left(\Gamma_{ts\beta}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}-\left(\Gamma_{lrs}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\left(\Gamma_{t\alpha\beta}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\right]\Bigg{\\}}$
$\displaystyle+$
$\displaystyle\sum\limits_{n=4}^{\infty}\int\prod\limits_{i=4}^{n}\cfrac{d^{4}k_{i}}{(2\pi)^{4}}\,h_{\rho_{i}\sigma_{i}}(k_{i})\,(2\pi)^{4}\delta\left(\sum
k_{i}\right)\kappa^{n-2}$
$\displaystyle\times\Bigg{\\{}\cfrac{m_{2}^{2}+2\,m_{0}^{2}}{3\,m_{0}^{2}\,m_{2}^{2}}\,\left(\sqrt{-g}\,g^{\mu\nu}g^{\alpha\beta}g^{\rho\sigma}g^{mn}g^{ab}g^{rs}\right)^{\rho_{5}\sigma_{5}\cdots\rho_{n}\sigma_{n}}(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}(k_{3})_{\lambda_{3}}(k_{4})_{\lambda_{4}}$
$\displaystyle\times\left[\\!\left(\Gamma_{\mu\alpha\rho}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\\!\\!\left(\Gamma_{\nu\beta\sigma}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\\!-\\!\left(\Gamma_{\mu\alpha\beta}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\\!\\!\left(\Gamma_{\nu\rho\sigma}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\\!\right]\\!\\!\left[\\!\left(\Gamma_{mar}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\\!\\!\left(\Gamma_{nbs}\right)^{\lambda_{4}\rho_{4}\sigma_{4}}\\!-\\!\left(\Gamma_{mab}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\\!\\!\left(\Gamma_{nrs}\right)^{\lambda_{4}\rho_{4}\sigma_{4}}\\!\right]$
$\displaystyle-\cfrac{2}{m_{2}^{2}}\left(\sqrt{-g}\,g^{\mu\alpha}g^{\mu\beta}g^{\rho\sigma}g^{\lambda\tau}g^{rs}g^{lt}\right)^{\rho_{5}\sigma_{5}\cdots\rho_{n}\sigma_{n}}(k_{1})_{\lambda_{1}}(k_{2})_{\lambda_{2}}(k_{3})_{\lambda_{3}}(k_{4})_{\lambda_{4}}$
$\displaystyle\times\\!\left[\\!\left(\Gamma_{\lambda\rho\mu}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\\!\\!\left(\Gamma_{\tau\sigma\nu}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\\!-\\!\left(\Gamma_{\lambda\rho\sigma}\right)^{\lambda_{1}\rho_{1}\sigma_{1}}\\!\\!\left(\Gamma_{\tau\mu\nu}\right)^{\lambda_{2}\rho_{2}\sigma_{2}}\right]\\!\\!\left[\\!\left(\Gamma_{lr\alpha}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\\!\\!\left(\Gamma_{ts\beta}\right)^{\lambda_{4}\rho_{4}\sigma_{4}}\\!-\\!\left(\Gamma_{lrs}\right)^{\lambda_{3}\rho_{3}\sigma_{3}}\\!\\!\left(\Gamma_{t\alpha\beta}\right)^{\lambda_{4}\rho_{4}\sigma_{4}}\\!\right]\\!\Bigg{\\}}.$

We have two reasons for not presenting the expression of the quadratic gravity
vertex. Firstly, one can extract the expression directly from the given
formula. Secondly, the expression would occupy a significant amount of space
without providing new information and would be challenging to present in an
easily understandable format. Nonetheless, it is helpful to show the
expression for the small metric perturbation propagator:

$\displaystyle\begin{split}\mathcal{G}_{\mu\nu\alpha\beta}(k)=&i\Bigg{[}\cfrac{1}{p^{2}}\left\\{-\cfrac{1}{2}\,P^{0}_{\mu\nu\alpha\beta}+\cfrac{2}{\epsilon}\,P^{1}_{\mu\nu\alpha\beta}+P^{2}_{\mu\nu\alpha\beta}+\left(\cfrac{4}{\epsilon}-\cfrac{3}{2}\right)\overline{P}^{0}_{\mu\nu\alpha\beta}-\cfrac{1}{2}\,\overline{\overline{P}}^{0}_{\mu\nu\alpha\beta}\right\\}\\\
&\hskip
35.0pt+\cfrac{1}{p^{2}-m_{0}^{2}}\left\\{\cfrac{1}{2}\,P^{0}_{\mu\nu\alpha\beta}+\cfrac{3}{2}\,\overline{P}^{0}_{\mu\nu\alpha\beta}+\cfrac{1}{2}\,\overline{\overline{P}}^{0}_{\mu\nu\alpha\beta}\right\\}-\cfrac{1}{p^{2}-m_{2}^{2}}\,P^{2}_{\mu\nu\alpha\beta}\Bigg{]}.\end{split}$
(82)

Here $P^{1}$, $P^{2}$, $P^{0}$, $\overline{P}^{0}$, and
$\overline{\overline{P}}^{0}$ are the Nieuwenhuizen operators [20] (the
definitions from [12] are used).

##  7 FeynGrav 3.0

The new FeynGrav version 3.0 is published online open access [1]. It
implements many new features while keeping the core of the package unchanged.

FeynGrav is an extension of FeynCalc [5, 6, 7, 8], so it depends on it and
requires it to function correctly. The main package FeynGrav.wl loads FeynCalc
and imports libraries describing interaction vertices from the Libs directory.
The package comes with some pre-generated libraries. Additional libraries can
either be generated or downloaded [62]. As noted in previous sections, the
number of terms in a vertex typically grows as $2^{n}n!$ where $n$ is the
number of gravitons involved in the interaction. Consequently, the libraries
for interactions involving many gravitons occupy much space. The package
provides essential libraries up to $\mathcal{O}\left(\kappa^{2}\right)$ order
while allowing the generation or downloading of the other required libraries.

The new version of the package allows the user to choose which libraries to
import. When the package initialises, it imports only libraries for general
relativity, the minimal gravitational coupling to scalars, Dirac fermions, and
vectors. Table 1 gives the list of commands importing libraries. Each command
has an option printOutput set to False by default. If the option is True, the
command prints additional output describing the import operation.

Command |  Description
---|---
\importGravitons[n] |  Imports graviton vertices in general relativity up to $\mathcal{O}\left(\kappa^{n}\right)$. If not specified, $n=2$. If the required libraries do not exist, imports up to the highest existing order.
\importScalars[n] |  Imports graviton-scalar vertices for the minimal coupling up to $\mathcal{O}\left(\kappa^{n}\right)$. If not specified, $n=2$. If the required libraries do not exist, imports up to the highest existing order.
\importFermions[n] |  Imports graviton-Dirac fermion vertices up to $\mathcal{O}\left(\kappa^{n}\right)$. If not specified, $n=2$. If the required libraries do not exist, imports up to the highest existing order.
\importVectors[n] |  Imports graviton-vector vertices for massless and massive fields up to $\mathcal{O}\left(\kappa^{n}\right)$. If not specified, $n=2$. If the required libraries do not exist, imports up to the highest existing order.
\importSUNYM[n] |  Imports graviton-$SU(N)$ Yang-Mills vertices up to $\mathcal{O}\left(\kappa^{n}\right)$. If not specified, $n=2$. If the required libraries do not exist, imports up to the highest existing order.
\importHorndeskiG2[] |  Imports vertices for the $G_{2}$ class of the Horndeski theory.
\importHorndeskiG3[] |  Imports vertices for the $G_{3}$ class of the Horndeski theory.
\importHorndeskiG4[] |  Imports vertices for the $G_{4}$ class of the Horndeski theory.
\importHorndeskiG5[] |  Imports vertices for the $G_{5}$ class of the Horndeski theory.
\importScalarGaussBonnet[n] |  Imports scalar-Gauss-Bonnet interaction vertices up to $\mathcal{O}\left(\kappa^{n}\right)$. If not specified, $n=2$. If the required libraries do not exist, imports up to the highest existing order.
\importAxionVectorVertex[n] |  Imports graviton-axion-vector interaction vertices up to $\mathcal{O}\left(\kappa^{n}\right)$. If not specified, $n=2$. If the required libraries do not exist, imports up to the highest existing order.
\importQuadraticGravity[n] |  Imports graviton vertices in quadratic gravity up to $\mathcal{O}\left(\kappa^{n}\right)$. If not specified, $n=2$. If the required libraries do not exist, imports up to the highest existing order.
Table 1: The list of commands that import libraries. Library | Orders
---|---
General relativity | up to $\mathcal{O}\left(\kappa^{2}\right)$
Minimally coupled scalars | up to $\mathcal{O}\left(\kappa^{2}\right)$
Dirac fermions | up to $\mathcal{O}\left(\kappa^{2}\right)$
Vectors | up to $\mathcal{O}\left(\kappa^{2}\right)$
SU(N) Yang-Mills | up to $\mathcal{O}\left(\kappa^{2}\right)$
Horndeski G2 | for number of scalars from 3 to 6, up to $\mathcal{O}\left(\kappa^{2}\right)$
Horndeski G3 | for number of scalars from 3 to 6, up to $\mathcal{O}\left(\kappa^{2}\right)$
Horndeski G4 | for number of scalars from 3 to 4, up to $\mathcal{O}\left(\kappa^{2}\right)$
Horndeski G5 | for number of scalars from 3 to 4, $\mathcal{O}\left(\kappa^{2}\right)$
Scalar-Gauss-Bonnet |  $\mathcal{O}\left(\kappa^{2}\right)$ only
Axion-Vectors | up to $\mathcal{O}\left(\kappa^{2}\right)$
Quadratic Gravity | up to $\mathcal{O}\left(\kappa^{2}\right)$
Table 2: The list of libraries provided with FeynGrav.

FeynGrav includes a separate package for generating libraries, which is
present in

> Libs/FeynGravLibrariesGenerator.wl.

It lacks a user-friendly interface, which will be improved in forthcoming
updates. The new version of FeynGrav enhances the performance of this package.
First and foremost, it uses the relations described in Section 2. Secondly,
the packages use the FORM symbolic manipulation system [63, 64], designed to
operate with symbolic expressions much faster than Wolfram Mathematica.
Consequently, the package depends on FORM, and the library generation
algorithm will not operate without it. However, if a user is not interested in
generating new libraries, the package works as intended without the FORM.
Table 2 provides the complete list of libraries included in the package.

The package responsible for the library generation relies on the Feynman rules
discussed in this and previous publications [2, 3, 4]. The rules are
implemented as additional subpackages and placed in Rules folder. These
packages are designed to generate expressions that will later be treated with
FORM, and we cannot recommend using them separately. Consequently, we shall
not discuss them in detail except the Nieuwenhuizen.wl. The Nieuwenhuizen.wl
subpackage provides tools to operate with the Nieuwenhuizen operators and
gauge projects. The subpackage is self-contained and can be used alone with
FeynCalc. The main FeynGrav package imports this package on initialisation and
allows users to use Nieuwenhuizen operators.

The new FeynGrav version excludes some commands describing propagators.
Previous versions employed commands realising propagators with FAD and SPD
commands. The FAD command is more appropriate for loop calculations since
FeynCalc algorithms recognise it, while SPD is more suitable for tree-level
calculations. The present version removed this feature since it is excessive.
The FeynCalc has a command FeynAmpDenominatorExplicit that allows one to
convert FAD to SPD, making the discussed commands obsolete. Table 3 presents
the list of all propagators implemented in FeynGrav.

Command | Description
---|---
ScalarPropagator | Propagator of a scalar.
ProcaPropagator | Propagator of a massive vector field.
GravitonPropagator | Graviton propagator.
GravitonPropagatorMassive | Massive graviton propagator.
QuadraticGravityPropagator | Graviton propagator within quadratic gravity.
Table 3: The list of all commands for propagators.

The latest version includes a variety of commands that describe gravitational
interaction vertices. Table 4 presents the complete list of these commands.
The package also offers graviton polarisation operators. However, the previous
publication [3] discuss them in detail, so we will not repeat the discussion
here.

Command |  Description
---|---
GravitonScalarVertex |  Vertex for a scalar field kinetic energy.
GravitonScalarPotentialVertex |  Vertex for a scalar field potential.
GravitonFermionVertex |  Vertex for a Dirac fermion.
GravitonMassiveVectorVertex |  Vertex for a massive vector field.
GravitonVectorVertex |  Vertex for a massless vector field.
GravitonVectorGhostVertex |  Vertex for the Faddeev-Popov ghost for a massless vector field.
GravitonGluonVertex |  Vertex for $2$, $3$, and $4$ gluons.
GravitonQuarkGluonVertex |  Vertex for the quark-gluon interaction.
GravitonYMGhostVertex |  Vertex for the Faddeev-Popov ghost for gluons.
GravitonGluonGhostVertex |  Vertex for the gluon Faddeev-Popov ghost interaction with gluons.
GravitonVertex |  Vertex for graviton interaction within general relativity.
GravitonGhostVertex |  Vertex for the Faddeev-Popov ghost for general relativity.
GravitonAxionVectorVertex |  Vertex for the axion-like coupling of a scalar field to a single vector field.
HorndeskiG2 |  Vertex for Horndeski $G_{2}$ interaction.
HorndeskiG3 |  Vertex for Horndeski $G_{3}$ interaction.
HorndeskiG4 |  Vertex for Horndeski $G_{4}$ interaction.
HorndeskiG5 |  Vertex for Horndeski $G_{5}$ interaction.
ScalarGaussBonnet |  Vertex for scalar-Gauss-Bonnet interaction between scalars and gravitons.
QuadraticGravityVertex |  Vertex for graviton interaction within quadratic gravity.
Table 4: The list of all commands describing interaction vertices.

In conclusion, the latest version of FeynGrav extends the toolkit for
perturbative quantum gravity. The previous version of FeynGrav was sufficient
to study quantum gravitational effects within the standard model. The new
version extends its functionality for a class of modified gravity models.

##  8 Conclusions and discussion

This paper presents the new version of FeynGrav [1] and develops the
theoretical framework behind it [2, 3, 4]. This paper discusses the following
issues and their role in the computational approach to perturbative quantum
gravity.

First and foremost, we discussed recurrent relations between tensors
describing the perturbative structure of the theory. From the technical point
of view, the perturbative structure of quantum gravity reduces to the inverse
metric $g^{\mu\nu}$ and the volume factor $\sqrt{-g}$. We established
recurrent relations for the perturbative structure of these quantities (see
also discussion in [4]). The recurrent relations for $g^{\mu\nu}$ have a pure
algebraic. Due to the scaling of small metric perturbations, we derived the
recurrent relations for the volume factor $\sqrt{-g}$. The discussion of the
fundamental role of this scaling and its implication lies far beyond the scope
of this paper. On the practical ground, the discussed recurrent relations
improve computational efficiency since they involve much fewer operations.
These recurrent relations are implemented in the current version of FeynGrav,
allowing it to extend its applicability to new orders of perturbation theory.

Secondly, we discussed the Horndeski gravity. The theory describes the most
general class of scalar-tensor models of gravity that admit second-order field
equations. The second order of the field equations ensures that the theory is
free from ghost degrees of freedom typically present in higher derivative
theories. The presence of healthy higher-order interactions makes this class
of models an essential object of study. On the practical ground, the theory
reduces to a few special cases, and we obtained the Feynman rules for all of
them. We also separately treat a scalar field coupling to the Gauss-Bonnet
term. Such interaction belongs to the Horndeski theory, but within the
standard parametrisation, it has a form not suitable for direct implementation
within a computational package. Because of this, we used a different
parametrisation of this interaction and obtained the corresponding interaction
rule.

Thirdly, we addressed the axion-like coupling. Axion physics is a separate
branch of studies addressing the problem of the $\theta$-term of quantum
chromodynamics. The development of axion physics gave rise to many models
involving axion-like coupling between a scalar and a vector field. We
addressed the simplest yet intensively studied case of an axion-like coupling
between a scalar field, a single massless vectors field, and gravity. The more
general case of a coupling to the $SU(N)$ gauge vector field will be addressed
in further publications. The discussed case is a necessary step in FeynGrav
development, providing grounds for further implementation of axion-like
interactions.

Further, we implemented the massive gravity propagator. Massive gravity can be
viewed as one of the simplest extensions of general relativity. The propagator
of massive gravity is well-known in literature and implemented in FeynGrav
directly. Further implementation of massive gravity will be discussed
elsewhere since it faces specific challenges. Massive gravity does not admit
the gauge symmetry, so it is impossible to fix its Lagrangian uniquely. At the
same time, describing the theory purely in geometrical terms is challenging,
so the formalism used in FeynGrav may not directly apply to that case.

Lastly, we considered the quadratic gravity. The theory is widely studied
because it provides a unique example of a renormalisable gravity model. At the
same time, the model admits a massive spin-$2$ ghost degree of freedom. The
theory is given by a Lagrangian that involves the scalar curvature squared and
the Ricci tensor squared terms. The formalism used in FeynGrav allows one to
obtain the interaction roles for this model since it is given entirely in
geometrical terms. Although such formulae contain many terms, the derivation
of Feynman’s rules presented no fundamental challenge. The corresponding
interaction rules are implemented in the latest version of the package.

The new version of FeynGrav provides another essential step in developing
tools for perturbative quantum gravity. The previous version was sufficient to
study quantum gravitational effects within the standard model. The new version
extends FeynGrav applicability to modified gravity. The interaction rules for
Horndeski gravity effectively allow one to study scalar-tensor gravity and all
models that can be mapped of scalar-tensor gravity. The axion-like coupling
provides a tool to study many gravity models with such coupling, which was
challenging before. Implementing the massive graviton propagator allows one to
address the simple scattering processes in massive gravity. Lastly, the
implementation of quadratic gravity gives yet another tool to address the
high-energy behaviour of that model.

Despite all the discussed developments, there are still challenges to be
addressed. The main challenge of perturbative quantum gravity is the
complexity growth. Since the perturbative quantum gravity operates within the
standard quantum field theory paradigm, the number of terms in an interaction
vertex typically grows with the number of involved gravitons $n$ as
$2^{n}\,n!$ since the expression shall enjoy certain symmetries. Such fast
growth is an essential feature of the quantum field theory and cannot be
removed from this approach. A growing number of terms is required for more
efficient computations.

The current version begins to address this challenge. We employed FORM to
generate expressions for the interaction vertices. However, new tools capable
of efficiently operating with higher-order expressions within the perturbative
quantum gravity are yet to be developed.

Additionally, we will continue to implement new gravity models within
FeynGrav. The models involving the $\theta$-term from the quantum
chromodynamics shall be implemented in further versions, allowing the study of
the quantum gravity role in axion physics. Similarly, we intend to implement
massive gravity models, although that may require extending the theoretical
framework beyond its current state.

In summary, FeynGrav v3.0 is the next step in developing computational tools
for perturbative quantum gravity. It implements new essential models and
extends the package’s applicability. Further development will provide more
tools to study quantum gravity and a more comprehensive range of models.

## Acknowledgement

This work was supported by the Institute for Basic Science Grant IBS-R018-Y1.

## References

  * [1] Boris N. Latosh.  FeynGrav.  https://github.com/BorisNLatosh/FeynGrav, 2024.
  * [2] Boris Latosh.  FeynGrav: FeynCalc extension for gravity amplitudes.  Class. Quant. Grav., 39(16):165006, 2022.  arXiv:2201.06812, doi:10.1088/1361-6382/ac7e15.
  * [3] Boris Latosh.  FeynGrav 2.0.  Comput. Phys. Commun., 292:108871, 2023.  arXiv:2302.14310, doi:10.1016/j.cpc.2023.108871.
  * [4] Boris N. Latosh.  FeynGrav and Recent Progress in Computational Perturbative Quantum Gravity.  Symmetry, 16(1):117, 2024.  arXiv:2401.05608, doi:10.3390/sym16010117.
  * [5] R. Mertig, M. Bohm, and Ansgar Denner.  FEYN CALC: Computer algebraic calculation of Feynman amplitudes.  Comput. Phys. Commun., 64:345–359, 1991.  doi:10.1016/0010-4655(91)90130-D.
  * [6] Vladyslav Shtabovenko, Rolf Mertig, and Frederik Orellana.  New Developments in FeynCalc 9.0.  Comput. Phys. Commun., 207:432–444, 2016.  arXiv:1601.01167, doi:10.1016/j.cpc.2016.06.008.
  * [7] Vladyslav Shtabovenko, Rolf Mertig, and Frederik Orellana.  FeynCalc 9.3: New features and improvements.  Comput. Phys. Commun., 256:107478, 2020.  arXiv:2001.04407, doi:10.1016/j.cpc.2020.107478.
  * [8] Vladyslav Shtabovenko.  FeynCalc goes multiloop.  J. Phys. Conf. Ser., 2438(1):012140, 2023.  arXiv:2112.14132, doi:10.1088/1742-6596/2438/1/012140.
  * [9] Y. Iwasaki.  Quantum theory of gravitation vs. classical theory. - fourth-order potential.  Prog. Theor. Phys., 46:1587–1609, 1971.  doi:10.1143/PTP.46.1587.
  * [10] John F. Donoghue.  General relativity as an effective field theory: The leading quantum corrections.  Phys. Rev. D, 50:3874–3888, 1994.  arXiv:gr-qc/9405057, doi:10.1103/PhysRevD.50.3874.
  * [11] C. P. Burgess.  Quantum gravity in everyday life: General relativity as an effective field theory.  Living Rev. Rel., 7:5–56, 2004.  arXiv:gr-qc/0311082, doi:10.12942/lrr-2004-5.
  * [12] A. Accioly, S. Ragusa, H. Mukai, and E. C. de Rey Neto.  Algorithm for computing the propagator for higher derivative gravity theories.  Int. J. Theor. Phys., 39:1599–1608, 2000.  doi:10.1023/A:1003632311419.
  * [13] Xavier Calmet.  Effective theory for quantum gravity.  Int. J. Mod. Phys. D, 22:1342014, 2013.  arXiv:1308.6155, doi:10.1142/S0218271813420145.
  * [14] David Prinz.  Algebraic Structures in the Coupling of Gravity to Gauge Theories.  Annals Phys., 426:168395, 2021.  arXiv:1812.09919, doi:10.1016/j.aop.2021.168395.
  * [15] David Prinz.  Gravity-Matter Feynman Rules for any Valence.  Class. Quant. Grav., 38(21):215003, 2021.  arXiv:2004.09543, doi:10.1088/1361-6382/ac1cc9.
  * [16] M. Fierz and W. Pauli.  On relativistic wave equations for particles of arbitrary spin in an electromagnetic field.  Proc. Roy. Soc. Lond. A, 173:211–232, 1939.  doi:10.1098/rspa.1939.0140.
  * [17] R. J. Rivers.  Lagrangian theory for neutral massive spin-2 fields.  Nuovo Cimento (1965-1970), 34:386–403, 1964.  doi:10.1007/BF02734585.
  * [18] Steven Weinberg.  Photons and Gravitons in $S$-Matrix Theory: Derivation of Charge Conservation and Equality of Gravitational and Inertial Mass.  Phys. Rev., 135:B1049–B1056, 1964.  doi:10.1103/PhysRev.135.B1049.
  * [19] Steven Weinberg.  Infrared photons and gravitons.  Phys. Rev., 140:B516–B524, 1965.  doi:10.1103/PhysRev.140.B516.
  * [20] P. Van Nieuwenhuizen.  On ghost-free tensor lagrangians and linearized gravitation.  Nucl. Phys. B, 60:478–492, 1973.  doi:10.1016/0550-3213(73)90194-6.
  * [21] Henriette Elvang and Yu-tin Huang.  Scattering Amplitudes in Gauge Theory and Gravity.  Cambridge University Press, 4 2015.
  * [22] Nima Arkani-Hamed, Tzu-Chen Huang, and Yu-tin Huang.  Scattering amplitudes for all masses and spins.  JHEP, 11:070, 2021.  arXiv:1709.04891, doi:10.1007/JHEP11(2021)070.
  * [23] Gabriele Travaglini et al.  The SAGEX review on scattering amplitudes.  J. Phys. A, 55(44):443001, 2022.  arXiv:2203.13011, doi:10.1088/1751-8121/ac8380.
  * [24] Eugene P. Wigner.  On Unitary Representations of the Inhomogeneous Lorentz Group.  Annals Math., 40:149–204, 1939.  doi:10.2307/1968551.
  * [25] Claudia de Rham.  Massive Gravity.  Living Rev. Rel., 17:7, 2014.  arXiv:1401.4173, doi:10.12942/lrr-2014-7.
  * [26] Kurt Hinterbichler.  Theoretical Aspects of Massive Gravity.  Rev. Mod. Phys., 84:671–710, 2012.  arXiv:1105.3735, doi:10.1103/RevModPhys.84.671.
  * [27] Gerard ’t Hooft.  Computation of the Quantum Effects Due to a Four-Dimensional Pseudoparticle.  Phys. Rev. D, 14:3432–3450, 1976.  [Erratum: Phys.Rev.D 18, 2199 (1978)].  doi:10.1103/PhysRevD.14.3432.
  * [28] Frank Wilczek.  Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons.  Phys. Rev. Lett., 40:279–282, 1978.  doi:10.1103/PhysRevLett.40.279.
  * [29] Steven Weinberg.  A New Light Boson?  Phys. Rev. Lett., 40:223–226, 1978.  doi:10.1103/PhysRevLett.40.223.
  * [30] Michael Dine, Willy Fischler, and Mark Srednicki.  A Simple Solution to the Strong CP Problem with a Harmless Axion.  Phys. Lett. B, 104:199–202, 1981.  doi:10.1016/0370-2693(81)90590-6.
  * [31] Michael Dine and Willy Fischler.  The Not So Harmless Axion.  Phys. Lett. B, 120:137–141, 1983.  doi:10.1016/0370-2693(83)90639-1.
  * [32] John Preskill, Mark B. Wise, and Frank Wilczek.  Cosmology of the Invisible Axion.  Phys. Lett. B, 120:127–132, 1983.  doi:10.1016/0370-2693(83)90637-8.
  * [33] P. Sikivie.  Experimental Tests of the Invisible Axion.  Phys. Rev. Lett., 51:1415–1417, 1983.  [Erratum: Phys.Rev.Lett. 52, 695 (1984)].  doi:10.1103/PhysRevLett.51.1415.
  * [34] K. S. Stelle.  Renormalization of Higher Derivative Quantum Gravity.  Phys. Rev. D, 16:953–969, 1977.  doi:10.1103/PhysRevD.16.953.
  * [35] K. S. Stelle.  Classical Gravity with Higher Derivatives.  Gen. Rel. Grav., 9:353–371, 1978.  doi:10.1007/BF00760427.
  * [36] J. Julve and M. Tonin.  Quantum Gravity with Higher Derivative Terms.  Nuovo Cim. B, 46:137–152, 1978.  doi:10.1007/BF02748637.
  * [37] E. S. Fradkin and Arkady A. Tseytlin.  Renormalizable asymptotically free quantum theory of gravity.  Nucl. Phys. B, 201:469–491, 1982.  doi:10.1016/0550-3213(82)90444-8.
  * [38] E. Elizalde, A. G. Zheksenaev, S. D. Odintsov, and I. L. Shapiro.  One loop renormalization and asymptotic behavior of a higher derivative scalar theory in curved space-time.  Phys. Lett. B, 328:297–306, 1994.  arXiv:hep-th/9402154, doi:10.1016/0370-2693(94)91483-4.
  * [39] E. Elizalde, S. D. Odintsov, and I. L. Shapiro.  Asymptotic regimes in quantum gravity at large distances and running Newtonian and cosmological constants.  Class. Quant. Grav., 11:1607–1614, 1994.  arXiv:hep-th/9404064, doi:10.1088/0264-9381/11/7/004.
  * [40] Peter M. Lavrov and Ilya L. Shapiro.  Gauge invariant renormalizability of quantum gravity.  Phys. Rev. D, 100(2):026018, 2019.  arXiv:1902.04687, doi:10.1103/PhysRevD.100.026018.
  * [41] Gregory Walter Horndeski.  Second-order scalar-tensor field equations in a four-dimensional space.  Int. J. Theor. Phys., 10:363–384, 1974.  doi:10.1007/BF01807638.
  * [42] Tsutomu Kobayashi, Masahide Yamaguchi, and Jun’ichi Yokoyama.  Generalized G-inflation: Inflation with the most general second-order field equations.  Prog. Theor. Phys., 126:511–529, 2011.  arXiv:1105.5723, doi:10.1143/PTP.126.511.
  * [43] M. Ostrogradsky.  Mémoires sur les équations différentielles, relatives au problème des isopérimètres.  Mem. Acad. St. Petersbourg, 6(4):385–517, 1850.
  * [44] Richard P. Woodard.  Ostrogradsky’s theorem on Hamiltonian instability.  Scholarpedia, 10(8):32243, 2015.  arXiv:1506.02210, doi:10.4249/scholarpedia.32243.
  * [45] Richard P. Woodard.  Avoiding dark energy with 1/r modifications of gravity.  Lect. Notes Phys., 720:403–433, 2007.  arXiv:astro-ph/0601672, doi:10.1007/978-3-540-71013-4_14.
  * [46] D. G. Boulware and Stanley Deser.  Can gravitation have a finite range?  Phys. Rev. D, 6:3368–3382, 1972.  doi:10.1103/PhysRevD.6.3368.