FeynMG: a FeynRules extension for scalar-tensor theories of gravity

The ability to represent perturbative expansions of interacting quantum field theories in terms of simple diagrammatic rules has revolutionized calculations in particle physics (and elsewhere). Moreover, these rules are readily automated, a process that has catalysed the rise of symbolic algebra packages. However, in the case of extended theories of gravity, such as scalar-tensor theories, it is necessary to precondition the Lagrangian to apply this automation or, at the very least, to take advantage of existing software pipelines. We present a Mathematica code FeynMG, which works in conjunction with the well-known package FeynRules, to do just that: FeynMG takes as inputs the FeynRules model file for a non-gravitational theory and a user-supplied gravitational Lagrangian. FeynMG provides functionality that inserts the minimal gravitational couplings of the degrees of freedom specified in the model file, determines the couplings of the additional tensor and scalar degrees of freedom (the metric and the scalar field from the gravitational sector), and preconditions the resulting Lagrangian so that it can be passed to FeynRules, either directly or by outputting an updated FeynRules model file. The Feynman rules can then be determined and output through FeynRules, using existing universal output formats and interfaces to other analysis packages.


Introduction
The increasing complementarity of high precision data from cosmological observations and high energy physics experiments makes it necessary to consider non-minimal gravitational couplings or the impact of additional degrees of freedom that are coupled through the gravitational sector with strengths that need not be Planck-suppressed.Examples include scalar-tensor theories of gravity [2], such as the Brans-Dicke theory [3] or, more generally, the Horndeski theories [4,5] (including beyond Horndeski [6,7] and DHOST [8,9] theories), in which the gravitational sector includes both the metric and an additional scalar degree of freedom.Other relevant examples include those in which the Higgs is non-minimally coupled to gravity, as is required in Higgs inflation [10][11][12][13][14][15][16][17] or so-called Higgs-Dilaton models [18][19][20][21].Indeed, such non-minimal couplings of the Higgs field to the scalar curvature are readily motivated by considering the renormalization group evolution of the operators of the Standard Model of particle physics plus gravity [22][23][24].Moreover, the ability to make Weyl rescalings of the metric and so-called disformal transformations [25][26][27] allows us to make connections between scalar-tensor theories of gravity and gauge-singlet, scalar extensions of the Standard Model of particle physics, such as Higgs-or neutrino-portal theories [28][29][30][31][32][33][34][35].
The challenge, however, is the proliferation of operators that non-minimal gravitational couplings provide, alongside degeneracies with operators that directly couple new degrees of freedom to the Standard Model.Dealing with this requires linearization of the extended gravitational sector, transformations of the metric, expansion around non-trivial vacuum configurations, the diagonalization of kinetic and mass mixings, and the truncation of infinite series of operators [36,37].This is usually done on a model-by-model basis, and it is a tedious and time-consuming process, which is ripe for automating, and doing so is the focus of this article.
We present a Mathematica package FeynMG, which is designed to work alongside the well-known FeynRules package [1].FeynRules is an extensive Mathematica package that enables the user to output the Feynman rules for a given Lagrangian in formats that can be read in by a range of high energy physics analysis software, including CalcHep/CompHEP [38,39], FeynArts [40], FeynCalc [41], FormCalc [42], MadGraph [43], Sherpa [44], Whizard/Omega [45] and ASperge [46].Symbolic algebra packages have also been developed to deal with the complex tensor algebra that arises in General Relativity.A recent example is FeynGrav [47], a package that introduces gravity in its canonical form (the Einstein-Hilbert action) to FeynRules.xAct [48] is perhaps the most wellknown package, having already been followed by multiple compatible packages that allow the study of gravity in different cosmological scenarios.In particular, the package xIST/COPPER [49] extends xAct for general scalar-tensor theories, and it was used in Ref. [50] to calculate the effect of modified gravity on cosmological perturbations.In this sense, FeynMG extends FeynRules as xIST/COPPER extends xAct.
FeynMG is intended as a 'preconditioner'.It takes as inputs a FeynRules model file and the Lagrangian of an extended gravitational sector.FeynMG then provides the functionality to implement the minimal gravitational couplings to the Lagrangian from the original model file and cast the complete theory in a form that can be further processed using the existing FeynRules package and its interfaces.However, we emphasise that FeynMG contains functionality that may be useful for theories that are being analysed independent of the couplings to gravitational sectors, and this will be highlighted throughout this article.
The remainder of this article is structured as follows.In Section 2, we describe the general form of the problem of coupling the Standard Model to extended gravitational sectors.We then present the package FeynMG, summarizing the implementation in Section 3 and describing its usage in Section 4. Finally, our conclusions are presented in Section 5, and additional technical details are provided in the Appendices.
Throughout this work, while it is a convention that is uncommon in the gravitation and cosmology literature, we use the "mostly minus" metric signature convention (+, −, −, −), in which timelike four-momenta p µ have p 2 > 0, since this is the convention commonly used by existing particle physics software packages.We use lower-case Greek characters for the Lorentz indices of the curved spacetime and lower-case Roman characters for the Lorentz indices of the flat, tangent space necessary for writing the Dirac Lagrangian in a generally covariant form.D denotes gauge covariant derivatives, general (i.e., gravitational) and gauge covariant derivatives are denoted by ∇, and an update to the general and gauge covariant derivative that is useful for scalar-tensor theories of Brans-Dicke type is represented by D. We work in natural units, but do not set Newton's gravitational constant to unity.

Method
We begin by reviewing how a Minkowski quantum field theory is minimally coupled to gravity and how additional scalar fields that are non-minimally coupled to the scalar curvature of the gravity sector can give rise to new interactions in that quantum field theory.
For simplicity, we work with a toy model of QED plus a real scalar prototype of the Higgs sector.Generalizing to a complex scalar field that is charged under U (1) would be a technical complication that does not add to the main points that we wish to illustrate below.The action of this model in Minkowski spacetime is given by where we have introduced a would-be Higgs field φ, a Dirac fermion ψ, which will later be chosen as a proxy for the electron, and the U (1) gauge field A µ , which corresponds to the photon, with its usual field-strength tensor Note that the Dirac fermion is charged under U (1), and it is minimally coupled to the photon field via the gauge covariant derivative where q is the electromagnetic coupling.Before analysing the interactions induced by extending the gravitational sector beyond the usual Einstein-Hilbert action, we first need to insert all the minimal gravitational couplings that have so far been ignored by working in Minkowski spacetime.This means that, for every pair of contracted Lorentz indices, we must include a factor of the metric g µν .Additionally, for every γ matrix appearing in the Dirac Lagrangian, we must include a vierbein e µ a , which satisfies η ab e µ a e ν b = g µν , where η ab is the flat spacetime metric.(We remind the reader that the flat-space indices of the vierbein are raised and lowered with the flat-space metric.)The latter is necessary since the algebra of the γ matrices is defined with respect to the Minkowski metric, i.e., {γ a , γ b } = 2η ab ; the vierbeins relate the curved and flat, tangent spaces.By this means, we obtain the minimally coupled action where we have also included a factor of √ −g in the spacetime volume element.Herein, the Minkowski gauge covariant derivative has been promoted to the general covariant derivative.
For scalar fields, the gravitational covariant derivative just trivially reduces to a partial derivative, such that ∇ µ φ → ∂ µ φ.However, when acting on a vector Y ρ , the covariant derivative takes the form where Γ ρ µν = 1 2 g ρλ (∂ µ g λν + ∂ ν g µλ − ∂ λ g µν ) are the usual Christoffel symbols.This definition for the covariant derivative is chosen such that ∇ ρ g µν = 0, but it can take many other forms.For instance, we will later define and work with a different choice that will be more convenient for the specific case of Brans-Dicke theories [37].However, it does not matter which definition one uses in this action, given that the following property will always hold since the curvature-dependent terms are symmetric under the permutation of µ and ν.Finally, the covariant derivative acting on a fermion field, including the dependence on the gauge field from QED, is given by where is the spin connection.The latter is defined by where With these minimal couplings now included, the action takes the form We can now proceed to append the gravitational sector.The minimal choice for the gravitational sector is the Einstein-Hilbert action, giving the full action where R is the Ricci scalar, and M Pl is the Planck mass, which determines the strength of the gravitational force.We can, however, also consider extended gravitational sectors, and one of the simplest examples is the Brans-Dicke scalartensor theory [3], in which a dynamical scalar field replaces the Planck mass.Such theories are described by an action with the following generic form: Herein, X is a real scalar field, subject to the self-interaction potential U (X) and coupled non-minimally to the Ricci scalar R through the function F (X). From a phenomenological perspective, there are tight constraints on the latetime evolution of Newton's gravitational "constant", e.g., from observations of the Moon's orbit [51].We must therefore choose the functions F (X), Z(X) and U (X), such that F (X) = M 2 Pl is approximately constant, e.g., by X obtaining an approximately constant vacuum expectation value (vev).Notice that the field X is not or, at least, does not appear to be canonically normalized, by virtue of the function Z(X) included in its kinetic term.In fact, additional contributions to the kinetic energy of the field X arise through the coupling to the scalar curvature.Moreover, while the matter sector does not contain any direct couplings to the field X, these couplings may be hidden in the mixing between the tensor and scalar degrees of freedom of the extended gravitational sector.The interactions between the field X and the would-be Standard Model fields become manifest once we have dealt with these mixings, and doing so is the main purpose of the package FeynMG.
For the Brans-Dicke example above, there are two ways that we can proceed, as will be described in the next subsections: 1. We can make a Weyl rescaling of the metric to remove the non-minimal gravitational coupling of the field X to the Ricci scalar, taking us to the so-called Einstein frame.2. We can continue in the Jordan frame (where the curvature couplings are manifest), by analysing how the metric degrees of freedom mediate interactions between the field X and our would-be Standard Model fields.
Before describing these two cases, however, it is important to note that the presence of additional non-minimal gravitational couplings, e.g., R µν ∇ µ ∇ ν X (as arises in the Horndeski class of scalar-tensor theories, where R µν is the Ricci tensor), the Weyl rescaling of the metric (or more generally a disformal transformation [25][26][27] of the metric) may not be able to remove all non-minimal couplings simultaneously.In these cases, we may not be able to transform into an Einstein frame and will have little choice but to continue working with nonminimal interactions with gravity.

Weyl transforming into the Einstein frame
Our aim is to isolate the new interactions between the matter fields that arise because of the modifications to the gravitational sector.The most common way of doing this is to transform to the Einstein frame.This amounts to a redefinition of the curvature-dependent objects (called a Weyl transformation) such that the resulting gravitational action does not present any non-minimal couplings.
For the Lagrangian defined in Eq. ( 9), this transformation will take the following form where gµν , ẽµ a and MPl are the metric, vierbein and Planck mass in the Einstein frame, respectively.To get through the algebra, the following transformations will be useful: where F (X) = ∂F (X)/∂X and all the curvature-dependent quantities with a tilde are built with the Einstein-frame metric gµν or vierbein ẽµ a .Applying these transformations to the Jordan-frame action, we obtain wherein we have recovered a canonical Einstein-Hilbert term for the gravitational action.However, all the couplings of the Brans-Dicke scalar arising from the modification of gravity now appear explicitly in the matter Lagrangian.Notice, in particular, that most of the kinetic energies of the fields are not canonically normalised due to these new couplings.
To canonically normalise the field X, we must solve the integral where X 0 is taken to be zero for simplicity.For the rest of the fields, we rescale them according to their classical scaling dimension, i.e., where F ( X) ≡ F (X).With this, the Lagrangian takes the following form: where Ũ ( X) ≡ U (X) and F ( X) = ∂ F ( X)/∂ X.Thus, one of the main inconveniences of working in the Einstein frame is that it loses the simplicity of the Lagrangian defined in the Jordan frame.This is because the Weyl transformation and the redefinition of the fields introduces factors of F ( X) throughout the Lagrangian, which, on making a series expansion of F ( X), will introduce infinite towers of operators that involve the SM fields and increasing powers of the scalar field X.
At this point, we can already make an important observation: The couplings between the SM fields and the scalar field X arise only through the scalar kinetic terms and terms with dimensionful parameters, i.e., those terms that are not invariant under Weyl transformations.Thus, for the Standard Model (illustrated already by the toy model described here), the modifications to the dynamics from the new scalar field X are, in the Einstein frame, communicated by the Higgs sector, with the squared mass parameter µ 2 of the tree-level Higgs potential playing the dominant role at low momentum exchange.In this way, there are strong parallels between the Brans-Dicke-type scalar-tensor theories and Higgs portal theories (see Ref. [36]).
Expanding the fields around their vacuum expectation values will give rise to kinetic and mass mixings between φ and X.Thus, when two fermions interact via their Yukawa coupling and exchange a would-be Higgs boson ( φ) in the t channel, there are two contributions to the central potential: a short-range interaction due to the heavy mode (the Higgs boson) and a long-range interaction due to the light mode (the light, additional scalar boson), see Ref. [36].Such long-range forces arising from the additional scalar fields of extended gravity sectors are often referred to as "fifth forces".In this way, even if the original matter Lagrangian is only minimally coupled to gravity in the Jordan frame, there can be experimentally testable modifications to the force laws that depend on the dynamics of the new scalar field.
Given how these new interactions manifest in the Einstein frame, it is instructive to consider how the same modifications to the dynamics manifest in the Jordan frame, without making the Weyl transformation (at least at first).This is the focus of the next subsection.

Staying in the Jordan frame
We can determine the modifications to the dynamics without performing a Weyl transformation to the Einstein frame and work directly in the Jordan frame.In this frame, new interactions between the fields of the matter sector arise through the gravity sector itself, and we proceed by perturbing the metric around a flat spacetime [52][53][54] in the gravitational weak-field limit.
Expanding the metric up to leading order in perturbations corresponds to where η µν is the usual flat spacetime metric and h µν is the perturbation in the metric, which, once quantized, corresponds to the graviton.The higher order terms in the expansion of g µν are necessary to satisfy g µν g νρ = δ ρ µ to all orders.For the gravitational sector of the Brans-Dicke-like theory [Eq.( 11)], with action we obtain the following expansion up to second order in the fields: It still remains to fix a gauge, and one choice is the harmonic gauge, which satisfies the following condition: This can be introduced into the Lagrangian through the term where Γ µ = g αβ Γ µ αβ .With this gauge choice, linearization of Einstein-Hilbert gravity leads to the familiar Fierz-Pauli Lagrangian [54], given by When working with Brans-Dicke theories in the Jordan frame, it is convenient to use a different gauge: one that maps to the harmonic gauge when performing the Weyl transformation to the Einstein frame. 1 This can be achieved by redefining the covariant derivative such that its action on a vector Y ν is as follows: where This modified covariant derivative will map to ∇ µ when going to the Einstein frame and satisfies the identity D ρ (F (X)g µν ) = 0 while preserving diffeomorphism invariance in the action, as shown in Ref. [37,[55][56][57].We can then define a scalar-harmonic gauge condition in terms of the new covariant derivative, namely This can be introduced into the Lagrangian as Expanding this gauge fixing term around a Minkowski background and adding it to the linearized gravitational sector from Eq. ( 20), we obtain Herein, we have recovered the usual kinetic energy terms of the graviton, as appear in the Fierz-Pauli Lagrangian (23), with the exception that non-minimal couplings to the field X appear through the overall factor of F (X). Notice that the Lagrangian (28) contains two additional terms relative to the Fierz-Pauli Lagrangian (23).The first contributes to the kinetic energy of the field X, which will have to be canonically normalized, and the second is a kinetic interaction between X and the trace of the graviton h.As we will show later, it this kinetic mixing that leads to additional interactions between the matter fields.On including the matter sector from the original action from Eq. ( 9), we get to the following Lagrangian after the linearization up to first order in 1/ F (X) where graviton self-interactions have been ignored and T µν is the energy-momentum tensor of the matter sector.The kinetic energy of the X field can be canonically normalized by defining where X 0 is taken to be zero for simplicity.Doing so leads to the Lagrangian where F (χ) ≡ F (X), F (χ) = ∂ F (χ)/∂χ and Û (χ) ≡ U (X).Now, we have only the graviton left to canonically normalise, since it is still non-minimally coupled to the function F (χ).However, as noted previously, the potential Û (χ) must lead to a non-vanishing vacuum expectation value for χ at late times so that the theory mimics Einstein gravity. 2 With this in mind, we shift χ → χ + v χ to obtain Figure 1: Series of diagrams contributing to the fifth force, and arising from the kinetic mixing between the graviton hµν and the scalar field χ.The ellipsis represents the series summing over all insertions of the kinetic mixing.
where higher-order terms in the interactions between χ and h µν have been omitted in the ellipsis.
The modification of gravity leads to a kinetic mixing between the trace of the graviton h and the χ field; the last term in the first line of Eq. (32).The example of the fifth force exchange described in the previous section then manifests in the Jordan frame through this mixing, as shown in Figure 1 (see Ref. [37]).
We can remove this mixing by the following transformation of the graviton and χ field: where F (v χ ) = M 2 Pl has been substituted and σ corresponds to the canonically normalized scalar field.This amounts to a perturbative implementation of the Weyl transformation, as is clear when one considers the resulting Lagrangian where T µ µ is the trace of the energy-momentum tensor.The fifth force arising from the final term in Eq. (34) will depend on the trace of the energy-momentum tensor of the interacting particles, leading to at most derivative interactions with σ for scale-invariant sectors [58].
We have seen that working in the Jordan frame requires us to linearize the gravitational sector and to diagonalize the fields, while in the Einstein frame, we had to perform the Weyl transformation and various rescalings of the matter fields, losing the simplicity of the potentials in the process.Whichever approach we take, the overall message of this section is not a discussion on which frame is best for calculations, as it is a matter of preference, but the fact that deriving Feynman rules for scalar-tensor theories is a tedious and time-consuming task, even for the simplest models.
This begs for a tool that helps us automate this process.In the rest of this paper, we will introduce the Mathematica package FeynMG within the FeynRules environment, which can efficiently perform manipulations on scalar-tensor theories of the types described in this section.

Implementation
FeynMG implements calculations of the type described in Section 2. The only necessary input is a model file compatible with FeynRules containing the matter Lagrangian and the description of all the existing fields and parameters.The user can then supplement this Lagrangian with their chosen scalar-tensor theory.
Scalar-tensor theories will generally give rise to both mass and kinetic mixings between fields.While FeynRules can deal with mass mixing if pre-defined in the model file, it cannot deal with kinetic mixing or cases where the form of the mass mixing is not known a priori.This is because FeynRules will ignore terms higher than quadratic order and will assume that all fields are canonically normalized.The scope of FeynMG is to linearize gravity and perform the necessary redefinitions to the fields such that it can be consistently used by FeynRules and all compatible packages.
We aim to make the code as easy as possible to use without losing the generality in the model files and desired gravitational actions.For example, for the input Lagrangian, it is possible to use an action defined in flat spacetime (i.e., reuse a FeynRules model file without modifying it).This is possible thanks to the function InsertCurv, which for every pair of contracted indices will insert a metric g µν or vierbein e µa , as appropriate, and promote partial derivatives to covariant derivatives.
Once all the minimal curvature dependencies are inserted into the Lagrangian, we need to append a gravitational action, wherein, e.g., the Ricci scalar can be specified using RScalar (see Appendix C.1 for the list of defined curvature objects).As is the case for FeynRules, it is necessary to identify fields and parameters.These attributes can be assigned to variables by using the functions AddScalar[] and AddParameter[], respectively, allowing complete freedom when creating the gravitational sector, and any number of new scalar degrees of freedom and parameters to be defined.In principle, the package should be able to deal with any gravitational sector, but it becomes more complicated the further away we go from Brans-Dicke theories.The effective Planck mass can be extracted at any point in the calculation by using the function GiveMpl.Moreover, using InsertMpl will calculate the effective M Pl in the action and substitute it into the expression.
As shown previously, in the particular case of Brans-Dicke gravity, we can perform a Weyl transformation such that the gravitational sector is of Einstein-Hilbert form and the matter action is instead dressed with additional scalar interactions.This is implemented in FeynMG by the function ToEinsteinFrame.However, more general scalar-tensor theories may not have an Einstein frame, forcing us to stay in the Jordan frame and proceed by linearizing gravity.The latter is implemented by the function LinearizeGravity, where the gravitational sector will be expanded up to second order, generating the kinetic energy for the graviton, and the matter sector will be expanded up to linear order in the interactions with the metric perturbation h µν .Moreover, the Jacobian √ −g will be automatically inserted, unless the option {Jacobian->0ff} is provided.
As described in the previous section, in the case of Brans-Dicke-like theories, it can be convenient to use the scalar-harmonic gauge from Eq. (27).By specifying the option {SHGauge->0n}, LinearizeGravity will determine the scalar-harmonic gauge fixing term and append it to the Lagrangian, depending on the specific coupling function F (X). 4 .This gauge choice will likely leave CMod terms in the linearized Lagrangian, corresponding to the modification of the Christoffel symbols Notice that the F (X) 2F (X) prefactor will have to be expanded in terms of X.Once this expansion is truncated at some order in X, we can no longer make a nonlinear redefinition of the X field (such as X → X 2 ), since the ignored higherorder terms will give contributions at lower orders.To avoid this problem, CMod won't be expanded until all the kinetic energies of the scalar fields have been canonicalized.
When dealing with tensor algebra, we are used to working with Einstein's index notation, for which the following holds: A µ A µ = A ρ A ρ .However, Mathematica will treat both terms A µ A µ and A ρ A ρ as distinct, since their indices are not represented by the same variable, leading to an overly complicated and long expression filled with repeated terms.The function IndexSimplify deals with this problem by replacing indices term by term from a user-supplied set of indices, so that the expression can be simplified using Mathematica's native functionality.
From here, which frame we use is unimportant, since the package has all the tools to leave the Lagrangian ready to be readable by FeynRules.If we stay in the Jordan frame (as may be necessary for theories that do not have an Einstein frame), one first needs to normalize the fields canonically.For scalar fields, the canonical normalization is implemented by the function CanonScalar, which will find and normalize the lowest-order derivative term of every field.In the case where the lowest order is already very complicated, one can use the in-built Mathematica function Series to perform a series expansion.
Similarly, we also need to normalize the graviton kinetic energy canonically.For that, depending upon the gravitational action, we might need to expand the fields around their vacuum expectation values, using VevExpand, which first calculates all the possible values for the vevs, and then shifts all the fields around the user's chosen branch of solutions.Once the graviton kinetic energy has a constant prefactor, we can then use CanonGravity, leaving all the fields canonically normalized with derivative interactions.As mentioned before, it will be at this point where all the CMod terms arising from the modified covariant derivatives will be expanded to make manifest their dependence on the additional scalar degree of freedom arising from the extended gravitational sector.
The only thing left to do is to deal with any mass or kinetic mixings that have arisen between any of the metric and scalar degrees of freedom.As mentioned previously, FeynRules assumes that all fields are canonical and only works with terms higher than quadratic order, so any mixing terms in the Lagrangian would be ignored in the outcome.To deal with this, MassDiagMG or KineticDiagMG diagonalizes the scalar field masses or kinetic energies, respectively.When proceeding in the Jordan frame, as we saw in the last section, the dominant modifications to the dynamics arise through kinetic mixing between the additional scalar field and the trace of the graviton (cf., e.g., Figure 1).The function GravKinMixing will calculate and substitute into the Lagrangian the field redefinitions that diagonalizes this kinetic mixing the equivalent of Eq. (33).With this, the Lagrangian should be in a form ready to be used by FeynRules.
Linearizing gravity and manipulating the Lagrangian into a form amenable to FeynRules can take significant computing time for extensive or complicated models.So that this process does not need to be repeated each time, the user can use the function OutputModelMG to create a new model file from the final form of the Lagrangian produced by FeynMG, which includes all the information about the redefined fields, the parameters of the extended model and the effective Lagrangian itself.This model file can then be used directly in FeynRules without the need to rerun the manipulations implemented by FeynMG.
To summarize, the package FeynMG provides a set of tools to help the user to upgrade the original FeynRules model file to one that includes the degrees of freedom of a canonical or extended gravitational sector.

Usage
In this section, we provide the instructions for loading FeynMG and using it to perform the manipulations described in the preceding sections.We will work in the Jordan Frame, given that the same tools can be used for the Einstein frame.In Appendix C, we provide a summary of the tools provided by FeynMG.

Installation
FeynMG has dependencies on FeynRules, so both packages need to be loaded into Mathematica to make use of FeynMG.This can be done by running The next step is to load a model file that is compatible with FeynRules using their function LoadModel[] (for an extensive description on how to build these files, see Ref. [1]).As mentioned previously, this model file does not need to include gravity in the defined fields or Lagrangians; these can be appended through FeynMG, as described earlier in Section 3.

Defining a gravitational action and transforming to the Einstein Frame
Throughout this section, we will work with the same Lagrangian from Eq. ( 9), whose matter sector is defined via Note that the last term of the first line corresponds to a generic covariant gauge for the U(1) gauge field.The first thing to do is to introduce the minimal gravitational couplings of this matter Lagrangian.This amounts to inserting metrics or vierbeins, as appropriate, for each pair of contracted indies, and promoting all partial derivatives to covariant ones.To implement this in FeynMG, we run VUp[mu,v3] γ c2 .γd1 .γv3 i1,j1 + [13→19] ..... Since the expressions can be long, we will show only the main sections of the output that motivate the next step in the calculation and represent the rest of the terms in ellipsis.To allow the reader to connect the output presented explicitly with the fill output of the code, the positions of the first and last terms omitted are specified over each ellipsis; the number in brackets at the end of the output represents the total number of terms in the full expression, i.e., 19 terms in Out [2].
For this example, we will introduce a Brans-Dicke gravitational sector of the form of Eq. ( 11), such that where the χ field should not be confused with the one defined in Eq. (30).Before defining the gravitational part of the Lagrangian within FeynMG, we need to give appropriate attributes to the additional field χ and the additional parameters ({ω, µ χ , λ χ }).In principle, these can be directly added by updating the model file itself (which should be done before loading it into FeynRules).Alternatively, the FeynMG functions AddScalar[] and AddParameter[], 5 allow the new scalar fields and parameters to be defined after the model file has been loaded into FeynRules.For the specific case of Eq. 37, we need to execute the following:

AddScalar[chi]; AddParameter[muC]; AddParameter[lamC]; AddParameter[w];
The full Lagrangian can then be defined via We note that the metric indices are raised by virtue of the specification Index[LorentzUp,mu] (for more information see Appendix C.1).Notice that we have not included the √ −g factor in the Lagrangian; this is because, for simplicity, FeynMG always assumes this term to be present.
In the case of Brans-Dicke-type scalar-tensor theories, it may be convenient to transform to the Einstein frame (see Section 2.1).This is achieved in FeynMG by executing The output agrees with the result from Eq. ( 14), including the last term, which comes from the fermion spin-connection [Eq.( 13)].As mentioned before, the Jacobian factor √ −g is assumed in the calculation (although it can be omitted by specifying the option {Jacobian→ → →Off} (see Appendix C.2 for further details).The gravitational sector is now of canonical Einstein-Hilbert form, and we can take the flat spacetime (Minkowski) limit by calling 11→20] ..... [20] wherein the couplings of the additional scalar field to the matter fields are manifest.The remaining fields are, however, not canonically normalized, and further manipulations are needed in order to pass this Lagrangian back to FeynRules.These are the focus of the next subsection.

Brans-Dicke theory for FeynRules in the Jordan frame
The calculation in the Jordan frame repeats the same steps as in the last subsection up to and including In [4]: We first need to load a model file.We then insert the curvature dependence using InsertCurv[] with the Lagrangian as the argument and provide a gravitational sector for the theory.The next step is to expand the metric about a flat spacetime background.This can be done by using where LJordan was defined previously in In [4], and the provided options specify that the scalar-harmonic gauge from Eq. ( 27) is used and all covariant derivatives are updated to the modified form from Eq. ( 24).As mentioned previously, the Jacobian √ −g has been included when linearizing gravity by default, but it can be omitted using {Jacobian→ → →Off} (see Appendix C.3).
As we can see in the second line, many of the terms are repeated, since Mathematica does not use Einstein's index notation, for which two repeated indices are summed over.As a result, various terms in the output will be equivalent, differing only in their index labels (e.g., A µ A µ = A ρ A ρ ).In order to force Mathematica to combine these terms, we have to use the same set of indices for all the terms.This problem is solved by the function IndexSimplify: The optional argument {mu,nu,rho} allows the user to choose a set of n indices from which the first n replacements will be chosen.The output of E2 contains significantly fewer terms than E1.Moreover, E2 already contains the expected graviton kinetic energy and its kinetic mixing with the scalar field chi, as in Eq. ( 29), thanks to the specification of the option {SHGauge->0n} in LinearizeGravity that implements the scalar-harmonic gauge and associated covariant derivatives from Eq. (35), which are convenient for the case of pure Brans-Dicke-type theories.This choice has led to the CMod[] terms in the Lagrangian, which needs to be series expanded around the chi field.However, the truncation to first order in chi does not commute with non-linear field redefinitions, so the CMod[] term will only be expanded once all the fields have their canonical kinetic energy.
We can check that the kinetic energies appearing in E2 are not canonically normalized by running There are one or more non-canonical kinetic energies.Use CanonScalar.

As the output indicates, we can execute
In [10] The kinetic energies of the scalar fields are now canonically normalized, leading to the expansion of every CMod[] (where present).This expansion is performed in terms of the scalar field χ.
At this stage, the kinetic energy of the graviton is composed of multiple terms.These could be simplified by means of Mathematica's FullSimplify command, but this will often prove time-consuming, and it is not necessary, except for aesthetic reasons.From here, the only thing left to do is to canonically normalize the graviton kinetic energy.To this end, we need to shift the fields around their vevs, so the graviton kinetic energy acquires a constant prefactor.This can be achieved by running In [11] , [8→10] ..... [10] Out Note that this function shows all the extrema of the potential.Since there may be multiple minima, the function allows the user to choose which vev (or set of vevs) will be used by a dialogue window prompt.(In this case, we choose option 7.) Notice that the v chi dependence already present from the expansion of the CMod functions have also been replaced by the user-selected vev in E3.
Once we have a constant prefactor to the graviton kinetic energy, we can canonically normalize it using In [13] We have recovered the usual canonically normalized Fierz-Pauli kinetic energy terms from Eq. (23).We also see the expected kinetic mixing between the scalar field and the graviton, which can be identified by executing Notice that a Yukawa coupling between the fermion fields and the chi field has appeared in the fourth line, as expected.However, a closer look at this term shows that the coupling constant is four times larger than the result m ψ / 2M 2 Pl (2 + 3w) from Refs.[36,37].This is because of the last term in the expression, which will also contribute to the tree-level interactions between the fermion and the scalar field, leading then to the same results as in Refs.[36,37].
At this point, all the interactions up to second order in the fields have been canonically normalized and diagonalized, so there are no kinetic or mass mixings.Therefore, the updated Lagrangian for the matter fields with the additional scalar field couplings is now in a form that can be processed further by FeynRules and compatible packages for phenomenological studies.

Outputting a model file
FeynMG allows the user to create a new model file with the Lagrangian of their choice, in which all the introduced particles (such as the graviton and additional scalar field) and new parameters (such as M Pl ) will be incorporated and properly defined. 6This can be done by running In [18] where OldModelFile is the name of the original FeynRules model file that the user loaded, NewModelFile is the chosen name of the new model file, and Lagrangian is the final Lagrangian, as prepared with FeynMG.
The upgraded model file can be read directly into FeynRules without needing to load or rerun FeynMG.

Conclusions
Modifying the gravitational sector of a Lagrangian can lead to new interactions between matter fields that need not be Planck-suppressed, but making these interactions manifest by hand on a model-by-model basis is tedious and time-consuming.In this paper, we have presented the Mathematica package FeynMG, which can manipulate scalar-tensor theories of gravity into a format that can be processed by FeynRules.
Even for the the simplest toy models, it is necessary to perform transformations of the metric or linearize the gravitational action, redefine multiple fields, expand around the vacuum expectation values of the scalar fields and diagonalize mass and/or kinetic mixings, in particular those between additional scalar field and the trace of the graviton.FeynMG provides a set of functions that allow the user to recycle existing FeynRules model files that does not contain gravity and to implement these various steps.
Once the user arrives at a canonically normalized Lagrangian, in which all kinetic and mass mixings have been diagonalized, it can be further processed by FeynRules and compatible packages to allow phenomenological studies of scalartensor theories of gravity.Moreover, instead of deriving the same Lagrangian every time one uses Mathematica, FeynMG allows the output of a new model file with all the updated fields, parameters and chosen Lagrangian.A summary list of functions can be found in Appendix C.
In this paper, we have described the implementation of a minimal example in FeynMG: Brans-Dicke theory coupled to QED plus a real scalar protoype of the Standard Model Higgs.The inbuilt functions, however, may be used to manipulate more complicated gravitational sectors, such as multi-field scalar-tensor theories or Horndeski theories, and additional functionality is being developed for future release.
where we have defined F (χ) ≡ F (X) and F (χ) ≡ ∂ F (χ)/∂χ.As described in Section 2.2, we now expand the scalar field around its vev, so that the graviton can also be canonically normalized [see Eq. (32), where the kinetic mixing between the graviton and χ is manifest].At this point, ∆ λµν has the form The kinetic mixing between the graviton and the scalar can be removed [see Eq. (33)] by means of the transformations in Eq. (33).With this, we obtain Eq. ( 34) and where F (v χ ) = M 2 Pl has been substituted and σ corresponds to the canonically normalized additional scalar field.We can now expand the denominator in the third line up to first order in M −1  Pl to give showing a perfect cancellation of the couplings to the additional scalar.Thus, after diagonalizing, the covariant derivative takes the following form which is nothing but the standard covariant derivative ∇ µ A ν from Einstein gravity.This is as we would expect, since the diagonalization is essentially a perturbative implementation of the Weyl transformation to the Einstein frame.
We can obtain the same result without diagonalizing and instead summing over all insertions of the graviton-scalar kinetic mixing.Our calculations have shown that the following two series of diagrams cancel with each other: where the ellipsis contains the sum over the infinite series of insertions of mixings (where zero kinetic mixing is also included for the diagram on the right).Similarly, from the diagrams above, we can calculate the incoming graviton amplitude by inserting an additional kinetic mixing to the left of the χ propagators.Thus, we find that all the diagrams containing kinetic mixings will end up cancelling each other, leaving just the diagram with no kinetic mixings.Diagrammatically, this implies that + = , which corresponds to the Feynman diagram for the coupling between the gauge field and gravity through the usual Chistoffel symbols.
In either case, we see that the role of the additional terms arising from C ρ µν in the updated covariant derivatives is to maintain the Weyl invariance of the Maxwell Lagrangian (at dimension four) once gauge fixing terms are included in the Jordan frame.

Appendix B. Diagonalizing graviton-scalar kinetic mixing
A convenient way to eliminate all the kinetic mixings is to find the matrix transformation that diagonalizes the kinetic terms.However, creating a kinetic mixing matrix between 2-forms (the graviton) and scalar fields is not straightforward.In this appendix, we describe a method for determining the transformation and diagonalizing the kinetic terms, which is implemented in FeynMG in the function GravKinMixing[].
The main obstacle is that the graviton kinetic term contains both h µν and its trace h.For example, we might have a Lagrangian of the form where both the graviton and the scalar field have already been canonically normalized, but there remains a kinetic mixing proportional to C (which for the calculation from Section 2.2 corresponds to C = F (v χ )/4).Since the graviton has two kinetic terms, it is unclear how to construct a matrix that encapsulates all the kinetic couplings between distinct fields.We proceed by redefining h µν so that its kinetic energy contains only one term.To do so, we perform an analytic continuation of the graviton into the complex plane, redefining The transformations for the fields are as follows: (F ρµν ) T KF ρµν = (F ρµν W −1 W ) T KW W −1 F ρµν = ( F ρµν ) T W T KW Fρµν , (B.7) since, by defining Fρµν = W −1 F ρµν , we would get a Lagrangian free of kinetic mixings.
For the generic kinetic mixing, where K is defined by Eq. (B.4), the transformation matrix is The scalar fields transform through F ρµν = W Fρµν and therefore hµν → hµν − iC √ 1 + 4C The new model file will contain the same defined fields and parameters as the original file, with the addition of all the new particles and parameters created using AddScalar and AddParameter, together with the Lagrangian (L), the graviton (h µν ) and Planck Mass (M Pl ).By specifying the option {UpdateMass→True}, the masses of all scalar fields will be updated.

[ 19 ]
Herein, gUp[a,b] and gDown[a,b] are upper-and lower-indexed metrics, respectively, VUp[a,b] and VDown[a,b] are upper-and lower-indexed vierbeins, respectively, and D Grav a [] is the gravitational covariant derivative.

h µν → hµν − 1 4 ( 1 + 2 ∂ 2 ∂
i) hη µν .ρ hµν ∂ ρ hµν + Ci∂ ρ h∂ ρ χ + 1 ρ χ∂ ρ χ, (B.3)which contains only one kinetic energy term for the graviton.The kinetic matrix is then defined straightforwardly as of the fields collected into the vectorF ρµν =   ∂ ρ hµν η µν ∂ ρ χ   , (B.5)such that the Lagrangian (B.1) can be written in the form L = (F ρµν ) T KF ρµν , where T denotes matrix transposition.We want a transformation W of the matrix K such thatW T KW = -Adds a new parameter named P into the loaded set of parameters, such that it can be recognized by FeynRules.Within the options (Opts), the user can choose its value by including {Value→X}.-Creates a new model file named NewF from an original FeynRules model file OldF.