FeynRules 2.0 A complete toolbox for treelevel phenomenology
Abstract
FeynRules is a Mathematicabased package which addresses the implementation of particle physics models, which are given in the form of a list of fields, parameters and a Lagrangian, into highenergy physics tools. It calculates the underlying Feynman rules and outputs them to a form appropriate for various programs such as CalcHep, FeynArts, MadGraph, Sherpa and Whizard. Since the original version, many new features have been added: support for twocomponent fermions, spin3/2 and spin2 fields, superspace notation and calculations, automatic mass diagonalization, completely general FeynArts output, a new universal FeynRules output interface, a new Whizard interface, automatic decay width calculation, improved speed and efficiency, new guidelines for validation and a new webbased validation package. With this feature set, FeynRules enables models to go from theory to simulation and comparison with experiment quickly, efficiently and accurately.
keywords:
Model building, Feynman rules, Monte Carlo programs.CERNPHTH/2013239, MCNET1314, IPPP/13/71, DCPT/13/142, PITTPACC1308
, , , ,
PROGRAM SUMMARY
Manuscript Title: FeynRules 2.0 A complete toolbox for treelevel phenomenology
Authors: Adam Alloul, Neil. D. Christensen, Céline Degrande,
Claude Duhr, Benjamin Fuks.
Program Title: FeynRules 2.0
Journal Reference:
Catalogue identifier:
Licensing provisions: None.
Programming language: Mathematica.
Computer: Platforms on which Mathematica is available.
Operating system: Operating systems on which Mathematica is available.
Keywords: Model building, Feynman rules, Monte Carlo simulations.
Classification: 11.1 General, High Energy Physics and Computing.
11.6 Phenomenological and Empirical Models
and Theories.
External routines/libraries: None.
Nature of problem: The program computes the Feynman rules of any
quantum field theory, expressed in fourdimensional spacetime, directly
from the Lagrangian of the model. Various interfaces to Feynman diagram
calculators are included that allow to export the interaction vertices in
a format readable by different Monte Carlo event generators or symbolic
calculation tools.
Solution method: FeynRules works in three steps:

If necessary, the model Lagrangian is written in terms of fourcomponent fermions and the usual fields of particle physics, instead of Weyl fermions or superfields.

Derivation of the Feynman rules directly form the Lagrangian using canonical commutation relations among fields and creation operators.

Implementation of the new physics model into FeynArts as well as into various Monte Carlo programs via dedicated interfaces.
Restrictions: Mathematica version 7.0 or higher. The Lagrangian
must fulfill basic quantum field theory requirements, such as locality
and Lorentz and gauge invariance. Fields with spin 0, 1/2, 1, 3/2
and 2 are supported.
Unusual features: Translation interfaces to various Feynman
diagram generators exist. Superfields are also supported and can be
expanded in terms of their component fields, which allows to perform
various sets of superspace computations.
Running time: The computation of the Feynman rules from a Lagrangian
varies with the complexity of the model, and runs from a few seconds to
several minutes. See Section 7 of the present manuscript for
more information.
1 Introduction
The current era of theoretical particle physics is on the cusp of making several important discoveries. Although the Standard Model (SM) has been amazingly successful at describing and predicting the outcome of most experiments, there are some significant experiments and observations that can not be explained by the SM. As a result, new asyet unknown physics beyond the SM is required and is expected to be discovered in the coming decades at current and future experiments. Among these puzzles is the instability of the Higgs boson generated vacuum to quantum corrections, the nature of dark matter, neutrino mass, the baryon/antibaryon asymmetry and the hierarchy of fermion masses.
The common feature among these puzzles is that each requires new theoretical models to be built and implemented into simulation software to be compared with experiment and observation. There are several powerful packages that simulate the effects of new models at particle colliders. Each of these has its strengths and may be appropriate for particular kinds of calculations. Among the multipurpose matrix element generators are CalcHep [1, 2, 3, 4], FeynArts/FormCalc [5, 6, 7, 8, 9], Helac [10, 11], MadGraph [12, 13, 14, 15, 16], Sherpa[17, 18] and Whizard [19, 20]. The partonlevel collision is often followed by the important step of radiation and hadronization by packages such as Herwig [21, 22, 23, 24], Pythia [25, 26, 27] and Sherpa. However, the step of implementing a new theoretical model into the simulation software has been difficult in the past. The chief problems have been that:

Each event generator uses its own syntax for a new model. Learning one syntax did not transfer to another generator;

Implementing a new model often required the modification of the event generator code itself. These implementations did not transfer well between theorists or to experimentalists;

These implementations required the hand coding of each vertex onebyone. Doing this for the hundreds or thousands of vertices in a model was tedious and very error prone.
LanHep [28, 29, 30, 31, 32], FeynRules [33] and Sarah [34, 35] were each created to overcome these challenges. In this paper we present version 2.0 of the FeynRules package. Even at the stage of the first version, FeynRules had a significant feature set, allowing theorists to quickly implement their new models into simulation packages. It allowed the theorist to implement their model once in a single unified syntax, which was Mathematica based and closely related to the syntax of FeynArts. It also allowed the user to enter the Lagrangian for their model rather than the individual vertices. It then did the work of calculating the vertices from the Lagrangian, thus saving theorists time and errors. It came with dedicated export interfaces which wrote the model to file in a form appropriate for the event generators CalcHep, FeynArts, MadGraph and Sherpa. The user was not expected to know the syntax of any of these event generators. Furthermore, for any supported vertex, the final product could be used in the standard version of these event generators and did not require modifications of their code. Thus, these implementations could be passed between theorists and experimentalists with great success. This early version was thoroughly tested both between different simulation tools, between different gauges and with independent implementations of a standard set of models [36]. A complete chain was established that went from theoretical model to partonlevel simulation package to hadronization, at which point comparison with experiment could be done [37]. Since then, many new important features have been added. We will summarize the major improvements now.
The core was updated to support twocomponent Weyl fermion notation [38]. The Lagrangian can now be written in terms of twocomponent fermions, both left and righthanded, which may greatly facilitate the implementation of complicated models. Matrix element generators, however, in general work at the level of fourcomponent fermions, and FeynRules has the capability to convert the Lagrangian to fourcomponent fermions in an automated way. See Section 2.4 and 3.2 for further details.
Support for spin was added to the core as well as to the CalcHep and UFO interfaces [39, 40]. The resulting output was tested in the context of three models. The first contained a spin Majorana gravitino, the second contained a spin Dirac topquark excitation, and the third contained a general effective operator approach for a spin Dirac quark excitation. The total cross sections, high energy growth cancellations and angular distributions were all tested with theory and agreement was found. See Section 2.4 for further details.
A superspace module was added [41] allowing the user to implement a supersymmetric (SUSY) model using the superfields directly. The user has the possibility to implement the superspace action, which can then be automatically expanded into component fields. In addition, the equations of motion for the auxiliary fields can be solved by FeynRules. The superspace module was tested on various nontrivial supersymmetric models, including the Minimal Supersymmetric Standard Model (MSSM) and the leftright symmetric supersymmetric model [42]. See Sections 2.5, 3.2, 4.5 and 4.6 for further details.
An automatic mass diagonalization module was created [43, 40] allowing FeynRules to extract automatically from any Lagrangian the treelevel mass matrices and to export them through the ASperGe interface for a numerical extraction of the spectrum. In practice, one reorganizes slightly the FeynRules model file and then uses the new builtin functions to extract the mass matrices and generate the ASperGe source code. The latter is an ensemble of C++ routines able to diagonalize the mass matrices and store the results in a SUSY LesHouches Accord (SLHA)like file. See Sections 2.8, 4.7 and 6.2 for further details.
The FeynArts interface was greatly extended to allow the implementation of operators that give rise to Lorentz structures that are not SM or MSSMlike, as is generally the case for higherdimensional operators. In particular, the new version of the interface creates, on the fly, the generic model file required by FeynArts to be able to use vertices with nonstandard Lorentz structures. The new interface was tested on dimensionsix operators for topquark production and decay and for anomalous electroweak gauge boson couplings. See Section 6.4 for further details.
The universal FeynRules output (UFO) interface was created [44]. The idea of this interface is to create a model format that is independent of any individual Feynman diagram calculator and which contains the full information contained in the FeynRules model. It is currently supported by MadGraph 5, MadAnalysis [45] and GoSam [46], and will be used in the future by Herwig++ [23]. See Section 6.7 for further details.
A new export interface to Whizard was written [47] that allows the output of models in the native format of Whizard. One highlight of this interface is that it not only supports unitary and Feynman gauge, but will automatically generate gauge model files from a model implemented in Feynman gauge. The resulting output was tested between gauges and with CalcHep and MadGraph for the SM, MSSM and 3Site model. In conjunction with these validations, the first phenomenological study with this interface was performed [48]. See Section 6.8 for further details.
A module for automatically computing the widths was added [49, 40]. FeynRules computes the squared matrix elements for all possible decays present in the model and multiplies by the appropriate phasespace factor. In addition, the analytic formulas can be included into the UFO output and used by the matrix element generators to obtain the treelevel twobody widths and branching ratios for all the particles defined in the model. See Section 4.8 for further details.
The speed and efficiency of FeynRules was greatly improved. Parts of the core were rewritten to make better use of more efficient Mathematica functions. Parallelization was also added to some of the most laborintensive parts of the code. See Section 7 for further details.
New guidelines were established to improve the quality of FeynRules based model implementations [38]. Additionally, a new webbased validation platform was created [50] to help accomplish this task. This validation platform allows any user to upload their model implementation. It then generates a list of processes which it compares between matrix element generators and between gauges. It also has the ability to compare with existing, independently implemented versions of the model, if available. The user is not required to know the details of the different matrix element generators. Processes with large discrepancies between gauges or between matrix element generators are flagged for the user to look into. See Section 8 for further details.
The purpose of this manual is to give an uptodate, integrated and complete instruction set for using FeynRules and all its new features. The organization is as follows. Sections 2 and 3 describe the implementation of the model and the Lagrangian. Section 4 describes how to run FeynRules in a Mathematica session. Section 5 gives a simple example of a model implementation. Section 6 describes the use of the export interfaces to generate code that can be run by one of the Feynman diagram calculators. Section 7 discusses the speed and efficiency improvements. Section 8 discusses the web based validation package and its use to improve the accuracy of the models built with FeynRules. In Section 9, we conclude.
2 The Model Description
To use FeynRules, the user must start by entering the details of the new model in a form that can be parsed by FeynRules. First of all, this means that, since FeynRules is a Mathematica package, the model must be written in a valid Mathematica syntax. Secondly, FeynRules specifies a set of special variables and macros for the user to enter each aspect of a new model. In this section, we will describe each of these in turn. These definitions can be placed in a pure text file or in a Mathematica notebook. The latter can be helpful during the development of the model details, however, we suggest storing the final version of any model in a pure text file.
2.1 Model Information
The user has the possibility to include general information about the model implementation
into the model file via the variables M$ModelName and M$Information.
While the first of these variables is a string giving the name of the model, the
variable M$Information
acts as an electronic signature of the model file
that allows to store, e.g., the names and addresses of the authors of the model
file, references to the papers used for the implementation, the version number of the
model file, etc. This information is stored as a Mathematica replacement list as shown in the following example
M$ModelName = "my_new_model";
M$Information = {
Authors > {"Mr. X", "Ms. Y"},
Institutions > {"UC Louvain"},
Emails > {"[email protected]", "[email protected]},
Date > "01.03.2013",
References > {"reference 1", "reference 2"},
URLs > {"http://feynrules.irmp.ucl.ac.be"},
Version > "1.0"
};
A summary and complete set of options available for M$Information
can be found in Table 1.
The model information will be printed on the screen whenever the model is loaded
into Mathematica. In addition, the contents of M$Information
can be
retrieved by issuing the command ModelInformation[]
in a Mathematica
session, after the model has been loaded.
2.2 Index Definitions
In general the Lagrangian describing a model is a polynomial in the fields (and their derivatives) as well as in the parameters of the model. Very often, these quantities carry indices specifying their members and/or how the different quantities transform under symmetry operations. For example, the gauge field of an unbroken gauge group carries two different types of indices:

a Lorentz index ranging from to ;

an adjoint gauge index ranging from to .
It is therefore crucial to define at the beginning of each model file the types of indices that appear in the model, together with the range of values each type of index may take.
A field carrying indices , , …is represented inside FeynRules by an expression of the form psi[, , …]. Each denotes an object of the form Index[
name, i]
, and represents an index of type name taking the value i. In this expression name is a symbol and value can be both a symbol or an integer. In general the name can be chosen freely by the user, but we emphasize that there are predefined names for the index types describing fourvectors (Lorentz), fourcomponent spinors (Spin) and twocomponent left and righthanded Weyl spinors (Spin1 and Spin2).
For FeynRules to run properly, the different types of indices that appear in the model have to be declared at the beginning of the model file, together with the range of values they can take. This is achieved like in the following examples
IndexRange[ Index[Colour] ] = Range[3]; IndexRange[ Index[SU2W] ] = Unfold[ Range[3] ]; IndexRange[ Index[Gluon] ] = NoUnfold[ Range[8] ];
These commands declare three types of indices named Colour
, SU2W
and Gluon
ranging form to and to respectively. The function Range
is an internal Mathematica command taking an integer n as input and returning the range .
Moreover, the indices of type Lorentz, Spin, Spin1 and Spin2 are defined internally and do not need to be defined by the user.
At this stage we have to comment on the functions Unfold
and NoUnfold
used in the declaration of the indices of type SU2W
and Gluon
:

The
Unfold
command instructs FeynRules that if an index of this type appears contracted inside a monomial, then it should be expanded, i.e., the monomial with the contracted pair of indices should be replaced by the explicit sum over the indices. Any index that expands in terms of nonphysical states must be wrapped inUnfold
. For instance, the indices in the Standard Model or in the Minimal Supersymmetric Standard Model must always be expanded in order to get the Feynman rules in terms of the physical states of the theory. Otherwise, wrong results could be obtained when employing matrix element generators. We refer to Section 4 for more details.
While indices are represented internally inside FeynRules by expressions of the form Index[
name, i]
, the user does not need to enter indices in this form. Since it is always possible to reconstruct the type of an index from its position inside the expression psi[, , …]. For example, the gluon field G[mu, a] has been declared as carrying
two indices, the first one being of type Lorentz and the second one of type Gluon
(see Section 2.4). FeynRules
can then employ particle class properties to restore the correct notation internally, as in
G[mu, a] G[Index[Lorentz, mu], Index[Gluon, a]] .
In addition, it is possible to specify how the different types of indices should be printed on the screen.
This is done via the IndexStyle
command, e.g.,
IndexStyle[ Colour, i ]; IndexStyle[ Gluon, a ];
Issuing these commands at the beginning of a model file instructs FeynRules to
print indices of type Colour
and Gluon
with symbols starting with
the letters i
and a
, respectively, followed by an integer number.
A summary of information for the index can be found in Table 2.
2.3 The model parameters
All the model parameters (coupling constants, mixing angles and matrices, masses, etc.) are implemented as elements of the list M$Parameters,
M$Parameters = { param1 == { options1 }, param2 == { options2 }, ... };
Each component of this list consists of an equality whose lefthand side is a label and the righthand side is a list of Mathematica replacement rules. The labels (param1 and param2 in the example) are userdefined names to be used when building the Lagrangian. The sets of replacement rules (options1 and options2 in the example) contain optional information allowing to define each parameter together with its properties. The model parameters are split into two categories according to whether they carry indices or not. We start by reviewing in Section 2.3.1 the implementation of scalar parameters, i.e., parameters that do not carry any index. Tensorial parameters, i.e., parameters carrying one or several indices, are then discussed in Section 2.3.2.
2.3.1 Scalar parameters
To illustrate the implementation of scalar parameters, we focus on the example of the strong coupling constant. The declaration of any other parameter is similar. Although the strong coupling constant usually appears in the Lagrangian, it is in general more convenient to use the quantity as an input parameter, since its numerical value (e.g., at the electroweak scale) has been precisely determined from experiments. It is therefore desirable to have both parameters in the FeynRules model file. This motivates us to choose, in our example, as a free parameter of the model, i.e., as an external parameter (in FeynRules parlance) or equivalently as an independent parameter. In contrast, is an internal parameter, or in other words, a parameter depending on one or several of the other internal and/or external parameters of the model. As a result, (aS) and (gs) could be implemented as in the example
aS == { TeX > Subscript[\[Alpha],s], ParameterType > External, InteractionOrder > {QCD, 2}, Value > 0.1184, BlockName > SMINPUTS, OrderBlock > 3, Description > "Strong coupling constant at the Z pole" }; gs == { TeX > Subscript[g,s], ParameterType > Internal, ComplexParameter > False, InteractionOrder > {QCD, 1}, Value > Sqrt[4 Pi aS], ParameterName > G, Description > "Strong coupling constant at the Z pole" }
The external or internal nature of a parameter can be specified by setting the attribute ParameterType to the value External or Internal. Another difference between external and internal parameters lies in the value taken by the attribute Value of the parameter class. For external parameters, it refers to a real number^{1}^{1}1Complex external parameters have to be split into their real and imaginary parts and declared individually. while for internal parameters, it contains a formula. This formula is given in standard Mathematica syntax and provides the way an internal parameter is connected to the other model parameters. It is important to note that only previously declared parameters can be employed when implementing the formula. Internal parameters can be either real or complex variables, which is specified by setting the attribute ComplexParameter to the value True or False (default).
Instead of providing the value of a parameter (a real number for an external parameter or a formula for an internal parameter) through the attribute Value, the user has the option to use the attribute Definitions. This option of the parameter class takes as argument a Mathematica replacement rule. Returning to our example above, we could replace the Value attribute of gs by
Definitions > { gs > Sqrt[4 Pi aS] }
The difference between these two choices appears at a later stage and concerns the derivation of the interaction vertices by FeynRules. A parameter provided with a definition is removed from the vertices and replaced by its definition, whilst otherwise, the associated symbol is kept.
The list of all the attributes of the parameter class for scalar quantities is given in Table 3 and Table 4. With the exception of the TeX attribute allowing for the TeXform of a parameter and the Description attribute which gives, as a string, the physical meaning of a parameter, the options nondescribed so far are related to the interfaces to Feynman diagram generators and discussed in greater detail in Section 6.
2.3.2 Tensorial parameters
The second category of parameters are tensorial parameters, i.e., parameters carrying indices. The index structure can be specified through the attribute Indices of the parameter class as in the following example,
Indices > {Index[Scalar], Index[Generation]}
The parameter under consideration carries two indices, one index of type Scalar and one index of type Generation. In many cases tensorial parameters correspond to unitary, Hermitian or orthogonal matrices. These properties can be implemented in the FeynRules model description by turning the values of the attributes Unitary, Hermitian or Orthogonal to True, the default value being False.
In contrast to scalar parameters, the attribute Value (Definitions) is now a list of values (replacement rules) for all possible numerical value of the indices. Moreover, tensorial parameters are by default complex quantities, i.e., the default value for the attribute ComplexParameter is True.
For example, the CabibboKobayashiMaskawa (CKM) matrix could be defined as follows,
CKM == { ParameterType > Internal, Indices > {Index[Generation], Index[Generation]}, Unitary > True, ComplexParameter > True, Definitions > { CKM[i_,3] :> 0 /; i!=3, CKM[3,i_] :> 0 /; i!=3, CKM[3,3] > 1 }, Value > { CKM[1,1] > Cos[cabi], CKM[1,2] > Sin[cabi], CKM[2,1] > Sin[cabi], CKM[2,2] > Cos[cabi] }, Description > "CKMMatrix" }
In these declarations, cabi stands for the Cabibbo angle, assumed to be declared previously in the model file. We have also simultaneously employed the attributes Value and Definitions so that vanishing elements of the CKM matrix are removed at the time of the extraction of the interaction vertices by FeynRules. In this way, zero vertices are removed from the output.
FeynRules assumes that all indices appearing inside a Lagrangian are contracted, i.e., all indices must come in pairs. However, it may sometimes be convenient to be able to break this rule. For example, in the case of a diagonal Yukawa matrix , the associated Lagrangian term
(2.1) 
where denotes a generic fermionic field and a generic scalar field, can be written in a more compact form as
(2.2) 
In the second expression, the Yukawa matrix has been replaced by its diagonal
form
, rendering the summation
convention explicitly violated since the index appears three times.
In order to allow for such violations in FeynRules, the attribute
AllowSummation
must be set to the value True when implementing the
Yukawa coupling . Consequently, this allows FeynRules to sum the single index
carried by the vectorial parameter along with the two other indices
fulfilling the
summation conventions (in our example, the indices of the fermions). We stress
that this option is only available for parameters carrying one single index.
In the case of contracted Weyl spin indices, the upper and lower position of the indices also plays an important role, since
(2.3) 
where and are two generic lefthanded Weyl fermions. This issue is described in details in Section 4.5.
2.4 Particle Classes
Particle classes can be instantiated in a similar way to parameter classes. The main difference is that, following the original FeynArts syntax [6], the particle classes are labeled according to the spins of the particles rather than to their symbol,
M$ClassesDescription = { spin1[1] == { options1 }, spin1[2] == { options2 }, spin2[1] == { options3 }, ...}
The symbols spin1, spin2, etc., refer each to one of the field type supported by FeynRules^{2}^{2}2The classes R, W and RW are specific to FeynRules and not supported by FeynArts.:

S
: scalar fields; 
F
: Dirac and Majorana spinor fields; 
W
: Weyl fermions (both left and righthanded); 
V
: vector fields; 
R
: fourcomponent RaritaSchwinger fields (spin3/2 fields); 
RW
: twocomponent RaritaSchwinger fields (both left and righthanded spin3/2 fields); 
T
: spin2 fields; 
U
: ghost fields (only complex ghosts are supported).
Similar to the declaration of the parameter classes, the quantities options1, options2, options3, etc., are sets of replacement rules defining field properties. Following the spirit of the original FeynArts model file format, each particle class should be thought of as a ‘multiplet’ consisting of particles that carry the same quantum numbers but might differ in mass. This implies that all fields belonging to the same class necessarily carry the same indices. The main advantage of collecting particles with the same indices into classes is that it allows the user to write compact expressions for Lagrangians. This is illustrated in the example Lagrangian
(2.4) 
where denotes the “quark class”, the strong coupling constant, the fundamental representation matrices of and stands for the gluon field. The notation of Eq. (2.4) avoids having to write out explicitly a Lagrangian term for each quark flavor.
Just like for the parameter classes, each particle class can be given a number of options which specify the properties of the field. In particular, there are two mandatory options that have to be defined for every field:

Each particle class must be given a name, specified by the
ClassName
option, which is at the same time the symbol by which the particle class is denoted in the Lagrangian. 
The option
SelfConjugate
specifies whether the particle has an antiparticle or not. The possible values for this option areTrue
orFalse
.
In addition, the particle class comes with various other options, as shown in Table 6. The set of all options can be divided into two groups:

options which are directly used by FeynRules in the Mathematica session when computing the Feynman rules. Examples include the indices carried by a field, its mass, etc.;

options which are mostly ignored by FeynRules, but describe information required by some Feynman diagram calculators. These options will be passed on to the Feynman diagram calculators via the corresponding FeynRules interfaces (see Section 6).
In the following we only describe the options directly used by FeynRules itself. The interface specific options are discussed in Section 6.
Besides the name of the particle class, the user may also provide names for the individual members of the class. For example, if uq denotes the class of all uptype quarks with members {u,c,t}, then this class and its members can be defined in the model file as
ClassName > uq, ClassMembers > {u, c, t}
If a field is not selfconjugate, an instance of the particle class associated with the antiparticle is created automatically by appending ‘bar’ to the name of the particle. For example, after the model file has been loaded into FeynRules, the symbols uqbar, as well as {ubar, cbar, tbar}, are available and can be used to construct a Lagrangian. For bosonic fields, the field representing the antiparticle corresponds to the Hermitian conjugate of the field, while for fermionic fields, the symbol psibar (psi representing any fermionic field ) corresponds to the quantity .
Fields transform in general as tensors under some symmetry groups, and they usually carry indices indicating the transformation laws. Just like for parameters, the indices carried by a field can be specified via the Indices option, e.g.,
Indices > {Index[ Colour ]} Indices > {Index[ Colour ], Index[ SU2D ]}
Indices labeling the transformation laws of the fields under the Poincaré symmetry (spin and/or Lorentz indices) are automatically inferred by FeynRules and do not need to be specified. In addition to indices labeling how a field transforms under symmetries, each field may have additional indices such as flavor indices. One of these can be distinguished as the flavor index of the class and labels its members. It is declared in the model file via the FlavorIndex option. For example, the uptype quark class uq previously introduced is usually defined carrying two indices supplementing the spin index (automatically handled by FeynRules), one of type Colour ranging from 1 to 3 and specifying the color of the quark, and another index of type Flavour ranging from 1 to 3. The latter is specified as the flavor index of the class (via the FlavourIndex option) so that it labels the members of the class,
Indices > { Index[ Colour ], Index[ Flavour ] }, FlavorIndex > Flavour
Quantum fields are not always only characterized by the tensor indices they carry, but also by their charges under the discrete and / or abelian groups of the model. FeynRules allows the user to define an arbitrary number of charges carried by a field, as, e.g., in
QuantumNumbers > {Q > 1, LeptonNumber > 1} QuantumNumbers > {Q > 2/3}
Next, the user can specify the symbol and the numerical value for the masses and the decay widths of the different members of a particle class using the Mass and Width options^{3}^{3}3In the following we only discuss the masses of the particles. Widths however work in exactly the same way..
The argument of Mass
is a list with masses for each of the class members,
as in
Mass > {MW} Mass > {MU, MC, MT} Mass > {Mu, MU, MC, MT}
where in the last example, the symbol Mu
is given for the entire class, while the symbols MU
, MC
and MT
are given to the members.
The symbol for the generic mass (Mu in this case) is by default a tensorial
parameter carrying a single index corresponding to the FlavorIndex of the
class. In addition, the AllowSummation property is internally set to
True. The user can not only specify the symbols used for the masses but also their numerical value as in
Mass > {MW, Internal} Mass > {MZ, 91.188} Mass > {{MU,0}, {MC,0}, {MT, 174.3}} Mass > {Mu, {MU, 0}, {MC, 0}, {MT, 174.3}}
If no numerical value is given, the default value 1 is assumed.
In the first example, MW
is given the value Internal
, which instructs FeynRules that the mass is defined as an internal parameter in M$Parameters
. This is the only case in which a user needs to define a mass in the parameter list. All other masses given in the definition of the particles should
not be included in this list.
While Lagrangians can often be written in a compact form in the gauge eigenbasis, the user usually wants to obtain Feynman rules in the mass eigenbasis. More generally, if the mass matrix appearing in the Lagrangian is not diagonal, one has to perform a field rotation that diagonalizes the mass matrix, i.e., the user has to replace certain fields by a linear combination of new fields such that in this new basis the mass matrix is diagonal. For fields in the gauge eigenbasis, the Mass and Width options do not have to be set, but these options are replaced by the option
Unphysical > True
The diagonalization of the gauge eigenbasis can be performed automatically by means of the ASperGe package which is interfaced to FeynRules. We refer to Section 4.7 and Section 6.2 for more details. The values of these rotation matrices can also be expressed analytically (if available) or as numerical values, relying on other tools to perform the diagonalization. The rotation to the mass eigenbasis in FeynRules are, in this case, implemented by using the Definitions option of the particle classes, which consists of a set of replacement rules that are applied to the Lagrangian before the computation of the Feynman rules. As an example, the relation between the hypercharge gauge boson (B) and the photon (A) and boson (Z) could be included in a model description as
Definitions > {B[mu_] > sw Z[mu] + cw A[mu]}
where the symbols sw and cw stand for the sine and cosine of the weak mixing angle and are declared alongside the other model parameters. The replacement is purely symbolic and can be performed at the Mathematica level, even if the numerical values of the mixing parameters are unknown.
We conclude this section by giving some options specific to certain kinds of fields.
For ghost particles, there is an option Ghost
which tells FeynRules the name of the gauge
boson the ghost field is connected to. There is a similar option Goldstone
in the case of scalar fields.
Majorana fermions are by definition eigenstates of the charge conjugation operator, and the associated eigenvalue must be a phase, since charge conjugation is unitary. If denotes a Majorana fermion, then
(2.5) 
The information on the phase can be provided by the option MajoranaPhase of the particle class,
MajoranaPhase > Phi
where Phi stands for the phase of the Majorana fermion . This phase can be further retrieved in the Mathematica session by typing
MajoranaPhase[lambda]
where the symbol lambda represents the Majorana fermion . The default value for the Majorana phase is 0.
Weyl fermion classes have an additional attribute that allows to specify their chirality. In the following we only discuss the case of spin1/2 twocomponent fermions (W), keeping in mind that the same attributes are also available for spin3/2 twocomponent fermions (RW). For example, the sets of rules
W[1] == { ClassName > chi, Chirality > Left, SelfConjugate > False }, W[2] == { ClassName > xibar, Chirality > Right, SelfConjugate > False }
define two Weyl fermions (chi) and (xibar), the first one being lefthanded and the second one righthanded. While Weyl fermions are in general not supported by any Feynman diagram calculator, it is often simpler to write down the Lagrangian in terms of twocomponent fermions and to transform them into fourcomponent spinors in a second step. For this reason, it is possible to add the option WeylComponents to an instance of the particle class F associated with fourcomponent Dirac and Majorana fermions. This option instructs FeynRules how to perform the mapping of the the twocomponent to the fourcomponent spinors. For example, the Dirac fermion could be defined in FeynRules as
F[1] == { ClassName > psi, SelfConjugate > False, WeylComponents > {chi, xibar} }
where the Weyl fermions chi and xibar are defined above. In the case of a fourcomponent Majorana fermion only the left handed component needs to be specified. In Section 4, we will see how we can instruct FeynRules to transform a Lagrangian written in terms of twocomponent spinors to its counterpart expressed in terms of fourcomponent fermions.
2.5 Implementing superfields in FeynRules
Supersymmetric theories can usually be written in a very compact form when using superfields, which are the natural objects to describe fields in supersymmetry. Feynman diagram calculators usually require a model to be implemented in terms of the component fields, i.e., the standard fields of particle physics. FeynRules allows the user to implement the model in terms of superfields, and the code then takes care of deriving internally the component field Lagrangian [41].
Most of the phenomenologically relevant supersymmetric theories can be entirely built in terms of chiral and vector superfields^{4}^{4}4Vector superfields must be constructed in the WessZumino gauge.. Superfields are declared in FeynRules in a way similar to ordinary fields,
M$Superfields = { superfield1[1] == { options1 }, superfield2[2] == { options2 }, ... }
The elements in the left hand side specify the type of superfield (similar to the way ordinary particle classes are defined by their spin). Currently, FeynRules supports two types of superfields:

CSF: chiral superfields;

VSF: vector superfields in the WessZumino gauge.
To illustrate the main features of the superfield declaration, we discuss the implementation of a lefthanded chiral superfield , a righthanded chiral superfield and a nonabelian vector superfield in the WessZumino gauge,
CSF[1] == { ClassName > PHI, Chirality > Left, Weyl > psi, Scalar > z, Auxiliary > FF } CSF[2] == { ClassName > XI, Chirality > Right, Weyl > psibar, Scalar > zbar } VSF[1] == { ClassName > VWZ, GaugeBoson > V, Gaugino > lambda, Indices > {Index[SU2W]} }
This declares two chiral superfields labeled by the tags CSF[1] and CSF[2] and one vector superfield labeled by the tag VSF[1]. For all three superfields, the symbols used when constructing a Lagrangian or performing computations in superspace are specified through the attribute ClassName. Moreover, as for Weyl fermions, the chirality of the fermionic component of a chiral superfield can be assigned by setting the attribute Chirality to Left or Right.
It is mandatory to match each superfield to its component fields. The fermionic and scalar components of a chiral superfield are hence referred to as the values of the attributes Weyl and Scalar, respectively, whilst the vector and gaugino components of a vector superfield are stored in the attributes GaugeBoson and Gaugino. Each of these options must point to the symbol associated to the relevant fields. In contrast, the declaration of the auxiliary ( and ) fields is optional. If not specified by the user, FeynRules takes care of it internally (as for the XI and VWZ superfields above). If the user however desires to match the auxiliary component of a superfield to an existing unphysical scalar field, he/she has to provide a value for the attribute Auxiliary (as for the PHI superfield above). All the components (psi, z, FF, V and lambda in the examples above) are assumed to be properly declared at the time of the field declaration as described in Section 2.4.
2.6 Gauge Groups
2.6.1 Gauge group declaration
The structure of the interactions of a model is in general dictated by gauge symmetries. Fields are chosen to transform in specific representations of the gauge group, and the invariance of the action under gauge transformation then governs the form of the interactions. In this section, we focus on the declaration of the gauge groups in FeynRules. Doing so allows the user to use several functions dedicated to an efficient implementation of the Lagrangian, such as, e.g., covariant derivatives or (super)field strength tensors.
The gauge group of a model is either simple or semisimple, and the different simple factors are defined independently in the list M$GaugeGroups,
M$GaugeGroups = { Group1 == { options1 }, Group2 == { options2 }, ... };
Each element of this list consists of an equality defining one specific direct factor of the full symmetry group. The lefthand sides of these equalities are labels associated to the different subgroups (Group1, Group2, etc.) while the righthand sides contain sets of replacement rules defining the properties of the subgroups (options1, options2, etc.). For example, let us consider the implementation of the SM gauge group in FeynRules,
M$GaugeGroups = { U1Y == { Abelian > True, CouplingConstant > gp, GaugeBoson > B, Charge > Y }, SU2L == { Abelian > False, CouplingConstant > gw, GaugeBoson > Wi, StructureConstant > ep, Representations > {{Ta,SU2D}}, Definitions > {Ta[a__]>PauliSigma[a]/2, ep>Eps} }, SU3C == { Abelian > False, CouplingConstant > gs, GaugeBoson > G, StructureConstant > f, Representations > {{T,Colour}}, SymmetricTensor > dSUN } }
Abelian and nonabelian groups are distinguished by the option Abelian, which may take the values True or False.
In the case of abelian groups, the user has the possibility to define the associated charge by setting the attribute Charge of the gauge group class to a symbol related to this operator. This allows FeynRules to check for charge conservation at the Lagrangian or at the Feynman rules level, assuming that the charges of the different (super)fields have been properly declared via the option QuantumNumbers of the particle class.
Nonabelian groups have specific attributes to define the generators and structure constants. For example, the symbols of the symmetric and antisymmetric structure constants (whose indices are those of the adjoint representation) of the group are defined as the values of the attributes SymmetricTensor and StructureConstant. Furthermore, the representations of the group necessary to define all the fields of the model are listed through the attribute Representations, which takes as a value a list of pairs. The first component of each pair is a symbol defining the symbol for the generator, while the second one consists of the type of indices it acts on. For example, in the declaration of the gauge group above, FeynRules internally creates a matrix of denoted by T[Gluon,Colour,Colour], where we have explicitly indicated the index types. The first index, represented here by the symbol Gluon, stands for the adjoint gauge index. Even if this information is not explicitly provided at the time of the gauge group declaration, FeynRules infers it from the value of the GaugeBoson (or Superfield) attribute (see below). Finally, the attribute Definitions allows the user to provide analytical formulas expressing representation matrices and structure constants in terms of the model parameters and/or standard Mathematica variables (see the declaration in the example above). As long as the user consistently provides definitions for each set of representation matrices that he/she wants to use, there is no limitation on those that can be used within FeynRules. However, care must be taken when exporting the model to a specific Monte Carlo generator as it in general has specific restrictions on the supported gauge groups and representations. Typically, these can be satisfied by fully expanding the vertices in terms of the component fields of the new gauge group. In addition, the FeynRules package comes with internal definitions for the most common representations for which we refer to Section 6.1.5 for more information.
At this stage we have to make a comment about how FeynRules deals with group representations. As should be apparent from the previous definition, FeynRules does not make any distinction between a representation and its complex conjugate. Indeed, if denote the generators of some irreducible representation and are the generator of the complex conjugate representation, then it is always possible to eliminate all the generators from the Lagrangian using the formula
(2.6) 
Furthermore, if a field transforms in the representation , then the antifield transforms in the representation . We can thus always choose the particles as the ones transforming in the representation . If the user wants to introduce an antitriplet, for example, he/she can add a field phi transforming as a triplet. Consequently, phibar transforms as an antitriplet and can be used as such when writing the Lagrangian.
Information on the gauge boson responsible for the mediation of the gauge
interactions is provided through the
GaugeBoson attribute of the gauge group class. It refers to the symbol
of the vector field associated to the group, defined separately in
M$ClassesDescriptions
. For nonabelian symmetries, the
(nonLorentz) index carried by this field consists of the adjoint index. As
briefly mentioned above, this is the way used by FeynRules to define adjoint
gauge indices. For supersymmetric models, it is also possible to link
a vector superfield instead of a gauge boson through the attribute
Superfield. If both the
Superfield and GaugeBoson attributes are specified by the
user, they must be consistent, i.e., the gauge boson must be the vector component
of the vector superfield.
The last attribute appearing in the examples above consists of the model parameter to be used as the gauge coupling constant, which is specified through the option CouplingConstant.
The list of all the options described above is summarized in Table 10.
2.6.2 FeynRules functions related to gauge groups
The declaration of a gauge group enables FeynRules to automatically construct field strength tensors, superfield strength tensors and covariant derivatives associated with this group, so that they can be further used when building Lagrangians. In the case of abelian gauge groups, the field strength tensor is invoked by issuing
FS[A, mu, nu]
where A
is the corresponding gauge boson and mu
and nu
denote Lorentz indices. Its supersymmetric counterparts can be called by the command
SuperfieldStrengthL[ V, sp ] SuperfieldStrengthR[ V, spdot ]
respectively. In these commands, V stands for the vector superfield associated with the gauge group and sp and spdot are lefthanded and righthanded spin indices. These three functions can be easily generalized to the nonabelian case,
FS[ A, mu, nu, a ] SuperfieldStrengthL[ V, sp , a ] SuperfieldStrengthR[ V, spdot, a ]
where a stands for an adjoint gauge index. Following the FeynRules conventions, these quantities are defined as
(2.7) 
where and denote the coupling constant and the structure constants of the gauge group and and are the superderivatives defined below, in Section 4.5. The abelian limit is trivially derived from these expressions. We emphasize that the spinorial superfields and are not hardcoded in FeynRules and are recalculated each time. However, they are evaluated only when an expansion in terms of the component fields of the vector superfield V is performed.
From the information provided at the time of the declaration of the gauge group,
FeynRules can also define, in an automated way,
gauge covariant derivatives. These can be
accessed through the symbol DC[phi, mu]
, where phi
is the field
that it acts on and mu
the Lorentz index. In our conventions, the
covariant derivative reads
(2.8) 
where stands for the representation matrices associated to the representation of the gauge group in which the field lies. The sign convention in Eq. (2.8) is consistent with the sign convention in Eq. (2.7).
All the functions presented in this section are summarized in Table 12.
2.7 Model restrictions
In phenomenological studies, it can sometimes be useful to consider restricted models which are obtained from a parent model by putting some of the external parameters to zero. As an example, one might be interested in the Standard Model
with a diagonal CKM matrix. While it is of course always possible to make the CKM matrix numerically diagonal, it is desirable to remove the interaction terms proportional to the offdiagonal terms altogether in order to avoid a proliferation of vanishing diagrams in Feynman diagram calculations.
This can be achieved by the use of the socalled restriction files in FeynRules. Restriction files are text files (with the extension .rst) that contain a list, named M$Restrictions
, of replacement rules to be applied to the Lagrangian before the evaluation of the Feynman rules. As an example, a possible restriction file to restrict the CKM matrix to a diagonal matrix reads
M$Restrictions = { CKM[i_,i_] > 1, CKM[i_?NumericQ, j_?NumericQ] :> 0 /; (i =!= j) }
Applying these rules to the Lagrangian (or to the vertices) obviously removes all the flavorchanging interactions among quark fields.
Restriction files can be loaded at run time after the model file has been loaded, i.e., after the LoadModel[ ] command has been issued. The corresponding command is
LoadRestriction[ file1.rst, file2.rst, ... ]
After this function has been called, the parameter definitions inside FeynRules are updated and the parameters appearing in the lefthand side of the replacement rules are removed from the model. In addition, the rules are kept in memory and are applied automatically to the Lagrangian when computing the Feynman rules. Note that this process is irreversible, and the restrictions cannot be undone after the LoadRestriction[] command has been issued (unless the kernel is restarted).
For complicated models, one often chooses the benchmark point in a special way such that many of the new parameters in the model are zero, and only a few new interactions are considered. FeynRules provides a command that allows to identify the parameters that are numerically zero and to create a restriction file that removes all the vanishing parameters from the model. More precisely, the command
WriteRestrictionFile[]
creates a restriction file called ZeroValues.rst which can be loaded just like any other restriction file. We emphasize that for complicated model the use of this restriction file may considerably speed up calculation performed with Feynman diagram calculators.
2.8 Mixing declaration
FeynRules comes with a module allowing for the derivation of the mass spectrum included in any Lagrangian (see Section 4.7). The mass matrices calculated in this way can be further passed to the ASperGe program [43] for a numerical evaluation of the necessary rotations yielding a diagonal mass basis (see Section 6.2). In principle, kinetic mixing among gauge bosons related to different factors of the gauge group is also possible. This is not supported by the new module, so that such a mixing has to be implemented manually by the user. We refer to Ref. [36] for a detailed example. In order to employ the functionalities related to the ASperGe package, the user has to implement particle mixings in a specific way, using instances of the mixing class collected into a list dubbed M$MixingsDescription,
M$MixingsDescription = { Mix["l1"] == { options1 }, Mix["l2"] == { options2 }, ... }
This maps a label given as a string ("l1", "l2", etc.) to a set of Mathematica replacement rules defining mixing properties (options1, options2, etc.). To illustrate the available options, we take the example of the and gauge fields of that mix, in the Standard Model, to the bosons
(2.9) 
As this mixing relation does not depend on any model parameter and is fully numerical, it can be declared in a very compact form,
Mix["l1"] == { MassBasis > {W, Wbar}, GaugeBasis > {Wi[1], Wi[2]}, Value > {{1/Sqrt[2],I/Sqrt[2]},{1/Sqrt[2],I/Sqrt[2]}} }
The previous mixing declaration should be understood as the matrix product
MassBasis = Value . GaugeBasis
Information on the gauge and mass bases is provided through the attributes GaugeBasis and MassBasis and the mixing matrix is given in a numerical form as the argument of the attribute Value. In the case the user already knows analytical formulas for the element of the mixing matrix, the syntax above cannot be employed. This matrix has instead to be declared as a standard model parameter (see Section 2.3) and referred to via the MixingMatrix attribute of the mixing class (see below). As the gauge basis is unphysical, it has to only contain fields declared as such, while the mass basis can in contrast contain either physical fields, unphysical fields or both. When unphysical fields are present, particle mixings are effectively implemented in several steps. We refer to Ref. [43] for examples. In order to implement the bases in a more compact form, spin and Lorentz indices can be omitted. Moreover, if some indices are irrelevant, i.e., if they are identical for all the involved fields, underscores can be employed, as in the example
Mix["l2"] == { MassBasis > {sdL[1,_], sdL[2,_], sdL[3,_]}, GaugeBasis > {QLs[2,1,_], QLs[2,2,_], QLs[2,3,_]}, ... }
Underscores reflect that the last index of each field has the same type and is not affected by the mixing (such as color indices for instance).
Most of the time, the mixing matrix is not known numerically before implementing the model, which leads to a slightly different syntax at the level of the implementation. The declaration
Mix["l3"] == { MassBasis > {A, Z}, GaugeBasis > {B, Wi[3]}, MixingMatrix > UW, BlockName > WEAKMIX }
describes the rotation of the third weak boson and the hypercharge gauge boson to the photon and boson states,
(2.10) 
where we have introduced the mixing matrix , the associated symbol being UW. The latter is referred to by means of the attribute MixingMatrix of the mixing class. The matrix UW does not have to be declared separately as this task is internally handled by FeynRules. Note that it is assumed that UW is a complex matrix with two indices, and according to our conventions we declare separately its real and imaginary parts. In other words, FeynRules creates two real external tensorial parameters, the real and imaginary parts of the matrix, and one internal complex tensorial parameter for the matrix UW itself. The user must specify the name of a Les Houches block [51, 52] which will contain the numerical values associated with the elements of the matrix (see Section 6.1.4). In the example above, we impose the real part of the elements of to be stored in a block named WEAKMIX and their imaginary part in a block IMWEAKMIX.
It must be noted that mixing matrices declared through a mixing declaration cannot be employed explicitly in Lagrangians (such as for the CKM matrix in the Minimal Supersymmetric Standard Model). If the user wants to use the mixing matrices explicitly in the Lagrangian, he/she must declare the matrices following the standard syntax (see Section 2.3), numerical values being provided from the beginning.
Finally, it is sometimes more practical to implement a mixing relation under the form
GaugeBasis = MixingMatrix . MassBasis
as for the CKM rotation in the Standard Model,
(2.11) 
where and are respectively gauge and masseigenstates. In order to avoid a cumbersome use of inverse matrices in mixing declarations, the user can make use of the optional attribute Inverse of the mixing class, setting it to True.
All the attributes of the mixing class are summarized in Table 13 while more involved cases are now detailed below.
When complex scalar fields are split into their real and imaginary degrees of freedom, it is necessary to declare one scalar and one pseudoscalar mass basis. In this case, the MassBasis attribute is given a list of the two bases. The first basis gives a list of the real scalar fields in the mass basis, while the second gives a list of the pseudoscalar fields. Consistently, the arguments of the attributes Value, BlockName, MixingMatrix and Inverse are now lists as well, the first element of these lists always referring to scalar fields and the second one to pseudoscalar fields. When some of the elements of these lists are irrelevant, as for instance when the scalar mixing matrix is unknown (MixingMatrix and BlockName are used) and the pseudoscalar mixing matrix is known (Value is used), underscores have to be used to indicate the unnecessary entries of the lists, such as in
Mix["scalar"] == { MassBasis > { {h1, h2}, {a1, a2} }, GaugeBasis > { phi1, phi2 }, BlockName > { SMIX, _ }, MixingMatrix > { US, _ }, Value > { _, ... } }
In the case of fourcomponent fermion mixing, two choices are possible. We recall that the mass terms (for Diractype fermions) can be generically written as
(2.12) 
all indices being suppressed for clarity. As a first option, the user can use a single gauge basis and FeynRules internally takes care of the chirality projectors appearing in the mass terms. In this case, the attribute GaugeBasis contains a list with the fields above. As a second option, two bases with different particle classes standing for the lefthanded and righthanded pieces of the fields can be employed. The user here sets the value of the GaugeBasis attribute to a twocomponent list. The first (second) element of this list is another list containing the names of the lefthanded (righthanded) components of the fields, i.e., the instance of the particle class describing the () fields above. Additionally, the arguments of the Value, BlockName, MixingMatrix and Inverse attributes are consistently upgraded to lists too, the first component being associated with lefthanded fermions and the second one with righthanded fermions. As for neutral scalar mixing, underscores can be used for irrelevant list elements. This second option is more commonly employed in the existing models although it may seem more complicated. We therefore give an explicit example:
Mix["dq"] == { MassBasis > {dq[1, _], dq[2, _], dq[3, _]}, GaugeBasis > { {QL[2, 1, _], QL[2, 2, _], QL[2, 3, _]}, {dR[1, _], dR[2, _], dR[3, _]} }, MixingMatrix > {CKM, _}, Value > {_, {{1,0,0}, {0,1,0}, {0,0,1}} }, Inverse > {True, _} }
In this example, we describe the implementation of the mixing relations associated with the downtype quarks in the Standard Model. The mass eigenstates are denoted by the three fields dq standing for the three generations of downtype quarks. The color indices being irrelevant, they are replaced by underscores everywhere. The lefthanded gauge eigenstates are the second piece of the doublets of quarks (the QL fields) whereas the righthanded gauge eigenstates are the singlets (the dR quarks). Finally, the last three arguments of the instance of the mixing class above indicate that the lefthanded fields mix trough Eq. (2.11) (the first elements of the attributes MixingMatrix and Inverse are set accordingly) and that the righthanded fields do not have to be rotated (the argument Value contains an identity matrix as its second element).
Lagrangian mass terms for charged Weyl fermions are generically written as
(2.13) 
where stands for the mass matrix and and are Weyl fermions with opposite electric charges. The diagonalization of the matrix proceeds through two unitary rotations and ,
(2.14) 
which introduces two mass bases. Therefore, all the attributes MassBasis, GaugeBasis, Value, MixingMatrix, and BlockName now take lists as arguments (with underscores included where relevant). It is mandatory to consistently order these lists.
Finally, in realistic models, the ground state of the theory is nontrivial and fields must be shifted by their vacuum expectation value (vev). The case of electrically neutral scalar fields can be treated in a specific way and included easily in particle mixing handling. Information on the vevs is included in a list of twocomponent elements denoted by M$vevs. The first of these elements refers to an unphysical field while the second one to the symbol associated with its vev, as in
M$vevs = { { phi1, vev1 }, { phi2, vev2 } }
the vacuum expectation values vev1 and vev2 being declared as any other model parameter (see Section 2.3). In this way, the information on the vevs of the different fields is consistently taken into account when FeynRules internally rewrites the neutral scalar fields in terms of their real degrees of freedom. For the sake of the example, we associate with the vev declaration above the following definition of the mixing of the and gauge eigenstates to the real scalar (pseudoscalar) mass eigenstates and ( and ),
Mix["phi"] == { MassBasis > { {h1, h2}, {a1, a2} }, GaugeBasis > { phi1, phi2 }, MixingMatrix > { US, UP } }
This teaches FeynRules to internally understand the mixing pattern of the fields as
(2.15) 
the scalar and pseudoscalar mixing matrices and being represented by the symbols US and UP, respectively.
3 The Lagrangian
An essential ingredient of a model implementation is the Lagrangian of the model. The Lagrangian can be built out of the fields of the models, augmented by some internal FeynRules and Mathematica functions. All the fields can be accessed by their ClassName, for example psi[a,b,...]
where a,b,...
are the indices of the field starting with its Lorentz indices for vector and tensor bosons, or its spin index for fermions followed by a Lorentz index for a spin fermion^{5}^{5}5Spin3/2 fields always have their spin index before their
Lorentz index..
The remaining indices are the ones defined in the particle definition through the Indices option, given in the same order. The different flavors can also be accessed using the names given in ClassMembers. They have the same indices as the full flavor multiplet, with the flavor index omitted. We recall that, if a field is not selfconjugate, FeynRules automatically creates the symbol for the
conjugate field by adding ‘bar
’ at the end of the particle name, i.e., the antiparticle associated to psi is denoted by psibar. For a fermion , the conjugate field is . Alternatively, the conjugate field can be obtained by issuing anti[psi]
.
Fields (and their derivatives) can be combined into polynomials. By convention, all the indices appearing inside a monomial in FeynRules must be contracted, i.e., all indices must appear pairwise^{6}^{6}6With the exception of singleindex parameters for which the AllowSummation option is set to True (see Section 2.3).. Furthermore, all indices must be spelled out explicitly. For anticommuting fields (fermions and ghosts), the Mathematica Dot function has to be used, in order to keep the relative order among them fixed. For example, the interaction between the gluon and all the down quarks can be written as
gs Ga[mu, s, r] T[a, i, j] dqbar[s, f, i].dq[r, f, j] G[mu, a]
There is however one case where indices do not need to be spelled out completely but can be omitted. If in a fermion bilinear, all the indices of the rightmost fermion are connected to all the indices of the leftmost fermion (perhaps with intermediate matrices), then these indices can be suppressed and FeynRules takes care of restoring them internally, such as in
dqbar.Ga[mu].T[a].dq
Ga[mu,s,r] T[a,i,j] dqbar[s,f,i].dq[r,f,j] .
In case of doubt, the user should always spell out all indices explicitly.
The Dot
product is mandatory for anticommuting fields or parameters. It should be noted that Mathematica does not keep the Dot
product between the components of vectors or matrices after computing their product explicitly
{ubar, dbar}.{u, d} = u ubar + d dbar
The appropriate treatment requires, therefore, use of the Inner function for each Dot, e.g.
Inner[Dot, {ubar, dbar}, {u, d}] = ubar.u + dbar.d
or for more than one multiplication,
Inner[Dot, Inner[Dot, {ubar, dbar} , {{a11, a12} , {a21, a22}}] , {u, d}].
As already discussed in Section 2.6, gauge invariant derivatives
can be conveniently defined via the functions DC[phi,mu]
and
FS[G, mu, nu, a]
. The first argument of both functions is the relevant
field, mu
and nu
are Lorentz indices and a
represents an
index of the adjoint representation of the associated gauge group. The gauge
fields and generators that appear in covariant derivatives of a particular field
are fixed by its indices and by the definition of the gauge group.
For example, the QCD Lagrangian for massless down quarks,
(3.16) 
is written as
L = 1/4 FS[G, mu, nu, a] FS[G, mu, nu, a] + I dqbar.Ga[mu].DC[dq, mu]
All the predefined FeynRules functions useful for the building of the Lagrangian are given in Table 14.
Finally, it is often convenient to write a Lagrangian in terms of twocomponent fermions and to let FeynRules perform the transformations to fourcomponent fermions. We note that this operation is mandatory for most Feynman diagram calculators, which in general only work with fourcomponent spinors. More precisely, if and are left and righthanded Weyl spinors, and is a Dirac fermion, we can easily switch to fourcomponent fermions by using the replacements
(3.17) 
These transformation rules are implemented in FeynRules via the WeylToDirac function, which takes as an argument a Lagrangian written in terms of twocomponent fermions, and returns the same Lagrangian in terms of fourcomponent fermions.
3.1 Tools for Lagrangians
FeynRules provides functions, collected in Table
15, that can be used while constructing
Lagrangians. For example, the function ExpandIndices[]
returns the Lagrangian with all the indices
written explicitly.
Each of the other functions return a different part of the Lagrangian as described in the table.
Once the Lagrangian is implemented, several sanity checks can be performed by means of the functions presented in Table 16. First, the function
CheckHermiticity[ ];
checks if the Lagrangian is Hermitian. Next, three functions are available to check if the kinetic terms and the mass terms are diagonal, CheckDiagonalKineticTerms, CheckDiagonalMassTerms and CheckDiagonalQuadraticTerms. Finally, two functions, CheckKineticTermNormalisation and CheckMassSpectrum, allow to check the normalization of the kinetic terms and compare the masses computed from the Lagrangian to those of the model description. The FeynRules conventions on the normalization of the kinetic and mass terms for the scalar, spin 1/2 and vector fields are

Scalar fields:

Real:

Complex (including ghost fields):


Spin1/2 fermions:

Majorana:

Dirac:


Vectors:

Real:

Complex:

FeynRules does not use the quadratic pieces of a Lagrangian. However, the propagators hardcoded either in FeynRules or in the event generators assume that the quadratic piece of the Lagrangian follow the abovementioned conventions. Furthermore, since the kinetic and mass terms for spin3/2 and spin2 fields are model dependent, they are therefore not implemented. Finally, checks on Weyl fermion kinetic and mass terms are also not supported since there exist several ways to write them down.
3.2 Automatic generation of supersymmetric Lagrangians
The implementation of supersymmetric Lagrangians in FeynRules can be highly facilitated by means of a series of dedicated builtin functions. The Lagrangian describing the kinetic terms and the gauge interactions of the chiral content of any supersymmetric theory is given by
(3.18) 
where the superfield content of the theory is represented by a set of chiral superfields and vector superfields . In the first line of the equation above, the indicates that one has to extract the highestorder coefficient in and from the expansion of the superfield expression included in the brackets. We recall that the covariant derivatives are defined in Eq. (2.8) and that the matrices stand for representation matrices of the gauge group relevant to the fields they act on. This Lagrangian can be computed in FeynRules by issuing
Theta2Thetabar2Component[ CSFKineticTerms[ ] ]
The function CSFKineticTerms returns the kinetic and gauge interaction Lagrangian terms for all the chiral supermultiplets of the model. As a result, the Lagrangian is computed in terms of superfields. The method Theta2Thetabar2Component is then needed to ensure that only the terms in e.g., the symbol Phi, the corresponding Lagrangian can be obtained by specifying this superfield as the argument of the function CSFKineticTerms, are kept, after the superfields have been expanded in terms of their component fields. If the user is interested in one specific superfield, represented by,
Theta2Thetabar2Component[ CSFKineticTerms[ Phi ] ]
Nongauge interactions among the chiral superfields are driven by the superpotential, a holomorphic function of the chiral superfields of the theory. The associated interaction Lagrangian is given by
(3.19) 
where in our notations, indicates that one has to select the component in from the expansion of the superpotential in terms of the Grassmann variables and . Assuming that is represented in the FeynRules model file by a variable SuperW, the Lagrangian can be calculated by issuing
Theta2Component[SuperW] + Thetabar2Component[HC[SuperW]]
where the two functions Theta2Component and Thetabar2Component allow to extract the correct coefficients of the expansion of the superpotential in and .
We now turn to the Lagrangian associated with the gauge sector. Kinetic and gauge interaction terms are built from squaring the superfield strength tensors and read, in the nonabelian case^{7}^{7}7The abelian limit can be easily derived. In the case where one has several abelian groups, kinetic mixing is often present. Such a feature is not covered in the automatic extraction of a supersymmetric Lagrangian by FeynRules and we leave to the user to implement it manually if he/she finds it necessary.,
(3.20) 
Implementing this Lagrangian into FeynRules can be achieved by issuing
LV = Theta2Component[ VSFKineticTerms[ ] ] + Thetabar2Component[ VSFKineticTerms[ ] ]
where an implicit sum over the whole vector superfield content of the theory is performed. The function VSFKineticTerms allows us to construct the Lagrangian in terms of superfields and the functions Theta2Component and Thetabar2Component to select the and components of the expansion of the superfield Lagrangian in terms of the Grassmann variables, respectively. Specifying the arguments of the function VSFKineticTerms,
LV = Theta2Component[ VSFKineticTerms[ VSF ] ] + Thetabar2Component[ VSFKineticTerms[ VSF ] ]
allows us to compute the kinetic and gauge interaction Lagrangian terms associated with the vector superfield represented by the symbol VSF.
To summarize, implementing a (renormalizable) Lagrangian density for a supersymmetric theory in FeynRules only requires the definition of the superpotential SuperW, since the rest of the Lagrangian computation is automatic. The full (supersymmetric) Lagrangian is thus calculated by issuing
LC=Theta2Thetabar2Component[CSFKineticTerms[]]; LV=Theta2Component[VSFKineticTerms[]] + Thetabar2Component[VSFKineticTerms[]]; LW=Theta2Component[SuperW]+Thetabar2Component[HC[SuperW]]; Lag = LC + LV + LW;
The Lagrangian density obtained in this way still depends on the auxiliary and fields that can be eliminated by inserting the solutions of their equations of motion back into the Lagrangian. This operation is automated in FeynRules by means of the SolveEqMotionD (this eliminates the fields), SolveEqMotionF (this eliminates the fields) and SolveEqMotionFD (this eliminates all auxiliary fields) functions. More precisely, all the auxiliary fields that could appear in the Lagrangian represented by the symbol Lag are eliminated by typing one of the two commands
SolveEqMotionFD[Lag] SolveEqMotionF[SolveEqMotionD[Lag]]
Finally, the component fields of a supermultiplet are Weyl fermions and need to rewritten in terms of fourcomponent spinors. This can be achieved using the WeylToDirac[] function described in the beginning of this section.
4 Running FeynRules
After the model description is created and the Lagrangian constructed, it can be loaded into FeynRules and the Feynman rules obtained.
4.1 Loading FeynRules
The first thing that must be done when using FeynRules is to load the package
into a Mathematica session. This should be done before any of the model
description is loaded in the kernel and so should be the first line of the
Mathematica notebook^{8}^{8}8In other words, if the model description is done
in a Mathematica notebook, it should come after FeynRules is loaded.. In
order to load FeynRules, the user must first specify the directory where it
is stored and then load it by issuing
FeynRulesPath = SetDirectory[ the address of the package ];
FeynRules`
4.2 Loading the model file
After the FeynRules package has been loaded^{9}^{9}9The user may want to
change
the current directory of Mathematica at this point. Otherwise, all the files
and directories generated by FeynRules may end up in the main FeynRules
directory., the model can be loaded using the command LoadModel
,
LoadModel[ < file.fr >, < file2.fr >, ... ]
The model may be contained in one model file or split among several files. For FeynRules to run properly, the extension of each model file should be .fr
. If the model description is entered directly in the Mathematica notebook, the
list of files is then empty. In this case, LoadModel[]
has to be executed
after all the lines of the model description are loaded into the kernel.
Any time the model description changes, the model must be reloaded. Currently, this means that the Mathematica kernel must be quit and the FeynRules package and model must be reloaded from the beginning. An exception to this is the Lagrangian. It can be changed and extended without having to reload the model information.
In the rest of this section, we describe the main utilities included in FeynRules which are summarized in Table 18.
4.3 Extracting the Feynman rules
After the model is loaded and the Lagrangian is defined, the Feynman rules can be extracted using the command FeynmanRules
. For the rest of this section, we use the QCD Lagrangian defined in Eq. (3.16) as an example. The
Feynman rules can be generated by means of the command^{10}^{10}10Since the vertices list may be very long, it is usually desirable to end this statement with a semicolon.:
vertsQCD = FeynmanRules[ LQCD ];
The vertices derived by FeynRules are then written out on the screen and stored internally in the variable vertsQCD
. The function FeynmanRules
has several options, that are described below.
The user can instruct Mathematica to not write the Feynman rules to the screen with the option ScreenOutput
set to False,
vertsQCD = FeynmanRules[ LQCD, ScreenOutput > False];
In this case, the Feynman rules are generated, stored in vertsQCD
, but not displayed on screen.
In the two previous examples, the flavors were not expanded. For example, the
preceding commands will only generate a single quarkgluon vertex (dq dqbar
G). It is often desirable to perform a flavor expansion, i.e., to obtain
separately the vertices d dbar G, s sbar G, and b bbar G. To
achieve this, the user can employ the FlavorExpand
option. This option can be used
to specify individual flavor indices to be expanded over, as in
FlavorExpand>Flavor
where only the Flavor
index is expanded over
(and not any other flavor indices if defined). It may also refer to a list of
flavor indices to be expanded, as in FlavorExpand>{Flavor,SU2W}
.
Similarly, it
may be used to expand over all the flavor indices as in FlavorExpand>True
. Note that it is always possible to expand over the flavors at a later stage using the FlavorExpansion[ ] function, i.e.,
vertices = FeynmanRules[ L ]; vertices = FlavorExpansion[ vertices ];
which is equivalent to calling FeynmanRules[ L, FlavorExpand > True ]. Note that calling the FlavorExpansion[ ] function in general runs faster than the FlavorExpand option. The difference is that when using the FlavorExpand option with FeynmanRules, the flavors are expanded before the vertices are computed, so FeynRules computes the vertices for the individual flavors separately.
At this stage we have to make a comment about the Unfold[ ] function introduced in Section 2.2. If an index type is wrapped in the Unfold[ ] command, then FeynRules will always expand over these indices. Moreover, any index which expands field in terms of nonphysical states must be wrapped in Unfold[ ]. As an example, the adjoint index carried by the weak gauge bosons must always be expanded over, because the physical states are the photon and the and bosons, while for quarks it can be interesting to keep the vertices in a compact form, keeping explicit flavor indices in the vertices. Note that, in order to speed up the computation of the vertices, (some of) the interfaces perform the flavor expansion using the FlavorExpansion[ ] command.
The list of Feynman rules can be quite long and it may sometimes be desirable to extract just one or a few vertices. There are several options available that limit the number of vertices computed by FeynRules:

MaxParticles > n instructs FeynmanRules to only derive vertices whose number of external legs does not exceed n. The option MinParticles works in a similar way.

MaxCanonicalDimension > n instructs FeynmanRules to only derive vertices whose canonical dimension does not exceed n. The option MinCanon icalDimension works in a similar way.

SelectParticles >{{...}, {...},...}
instructsFeynmanRules
to only derive the vertices specified in the list. For example, the command:FeynmanRules[ LQCD, SelectParticles >{{G,G,G}, {G,G,G,G}}]
only leads to the derivation of the three and fourpoint gluon vertices.

Contains > { ... }
instructsFeynmanRules
to only derive the vertices which involve all the particles indicated in the list. 
Free > { ... }
instructsFeynmanRules
to only derive the vertices which do not involve any of the particles indicated in the list.
FeynmanRules
, by default, checks whether the quantum numbers that have been defined in the model file are conserved for each vertex. This check can be turned off by setting the option ConservedQuantumNumbers
to False
. Alternatively, the argument of this option can be a list containing all the quantum numbers FeynRules should check for conservation.
The Feynman rules constructed are stored internally as a list of vertices where each vertex is a twocomponent list. The first element enumerates all the particles that enter the vertex while the second one gives the analytical expression for the vertex. We illustrate this for the quarkgluon vertex,
{{{G, 1}, {qbar, 2}, {q, 3}},
}
Moreover, each particle is also given by a twocomponent list, the first element being the name of the particle while the second one the label given to the indices referring to this particle in the analytical expression.
As the list of vertices computed by FeynRules can be quite long for complicated models, the SelectVertices
routine can be employed.
The options shared with the function FeynmanRules are MaxParticles
, MinParticles
, SelectParticles
, Contains
and Free
and
can be invoked as
vertsGluon = SelectVertices[vertsQCD, SelectParticles>{{G,G,G},{G,G,G,G}}];
It is sometimes convenient to construct the Feynman rules by parts, for instance, when one splits the QCD Lagrangian into three pieces as in
LQCD = LGauge + LQuarks + LGhosts;
The Feynman rules can be constructed all at once as in the previous examples, or they can be constructed separately as in:
vertsGluon = FeynmanRules[ LGauge ]; vertsQuark = FeynmanRules[ LQuarks ]; vertsQuark = FeynmanRules[ LGhosts ];
They can later be merged using the function MergeVertices
vertsQCD = MergeVertices[ vertsGluon, vertsQuark, vertsGhosts ];
This merges the results obtained for vertsGluon
, vertsQuark
and verts Ghosts into a single list of vertices. If there are two contributions to the same vertex, they will be combined into one vertex.
4.4 Manipulating Parameters
The parameters are also an important part of the model and several functions
allow to manipulate them. The numerical values of any parameter or function of
one or several parameters can be obtained with the use of the
NumericalValue
function, as for example in
NumericalValue[ Sin[ cabi ]]
where cabi
is the Cabibbo angle. This returns the numerical value of the sine of the Cabibbo angle.
In the case the user wants to change the value of a subset of external parameters,
the function UpdateParameters
can be used such as in
UpdateParameters[ gs > 0.118 , ee > 0.33 ]
where gs
and ee
are the strong and electromagnetic coupling constants, respectively.
In order to write and read the numerical values of the parameters to and from a
file, FeynRules comes with the two functions WriteParameters
and
ReadParameters
. By default, they write and read the parameters to and from the file named M$ModelName
with a ‘‘.pars’’ appended, but this can be changed with the options Output
and Input
as in:
WriteParameters[Output > "parameters.pars"] ReadParameters[Input > "parameters.pars"]
The routine WriteParameters
writes out the external and internal parameters (including masses and widths) to the specified file, while the function ReadParameters
only reads in the external parameters (including masses and widths). This gives the user another way to change the values of the external parameters. The user can modify the values of the parameters included in the parameter file created by WriteParameters
, and then read them back in using ReadParameters
. Any changes made to the internal parameters are ignored.
In addition to the native FeynRules parameter files (.pars), there is a second way to update numerical values of external parameters. Indeed, FeynRules is equipped with its own reading and writing routine for Les Houches Accord (LHA)like parameter files [51, 52]. For example, issuing the command
ReadLHAFile[ Input > "LHAfile" ]
reads the LHAlike parameter file "LHAfile"
and updates the numerical values of all the external parameters of the model in the current Mathematica session. Similarly, the command
WriteLHAFile[ Output > "LHAfile" ]
writes all the numerical values to file in a LHAlike format.
4.5 Manipulating superspace expressions
The superspace is defined by supplementing to the ordinary spacetime coordinates a Majorana spinor . These extra coordinates are represented in FeynRules by the symbols theta and thetabar
By convention, both spin indices are assumed to be lower indices. The position of the spin indices can be modified by employing the ranktwo antisymmetric tensors represented by the objects Deps and Ueps, which equivalently take as arguments lefthanded or righthanded spin indices. For instance, one could implement the following expressions in a Mathematica session as
Proper computations in superspace require to keep track, on the one hand, of the position of the spin indices and, on the other hand, of the fermion ordering. To this aim, FeynRules always assumes that an explicit spin index is a lower index. Moreover, fermion ordering information is provided by using the syntax nc[chain], where chain stands for any ordered sequence of fermions (with lower spin indices). As a simple example, we show a possible implementation for the scalar product , where and are righthanded Weyl fermions,
This way of inputting superspace expressions is however highly unpractical for longer expressions. Therefore, FeynRules contains an environment ncc similar to the nc environment, the difference lying in the fact that all spin indices and tensors can be omitted,
The ncc environment is automatically handled by FeynRules. The result of the command above consists of the same expression, but given in canonical form, employing ranktwo antisymmetric tensors and twocomponent spinors with lowered indices, ordered within a nc environment.
Another consequence of the mandatory usage of this canonical form is that the Mathematica output associated to a superspace expression is in general difficult to read. To improve the readability, it is possible to force FeynRules to form invariant products of spinors and to simplify products of Grassmann variables according to the relations
(4.21) 
deduced from the Grassmann algebra fulfilled by the variables and . This operation is achieved by means of the function ToGrassmannBasis which takes any expression exp, depending on the superspace coordinates, as an argument,
ToGrassmannBasis[ exp ]
First of all, this function rewrites the expression exp in terms of a restricted set of scalar products involving Grassmann variables and Pauli matrices, and, second, optimizes the index naming scheme used. This allows to collect and simplify Mathematica expressions that are equal up to the names of contracted indices, such as, e.g.,
Dot[lambda[al], lambda[al]]  Dot[lambda[be], lambda[be]]
The output of the ToGrassmannBasis function consists therefore of expressions very close to their original forms. This method also works on tensorial quantities containing noncontracted spin indices. As a result, the ToGrassmannBasis method matches the free indices either to single fermions or to chains containing one fermion and a given number of Pauli matrices. For instance, applying the ToGrassmannBasis function on , being a lefthanded Weyl fermion, leads to
One observes that a chain containing one Pauli matrix and the fermion has been formed and stored in a TensDot2 environment. The structure related to this environment is defined by
TensDot2[ chain ][pos, chir, name]
where chain is a sequence of one Weyl fermion and possibly one or several Pauli matrices, pos is the up or down position of the free spin index, chir indicates if it is a lefthanded (Left) or righthanded (Right) index and name denotes the symbol associated to the index.
The optimization of the index naming scheme performed by the ToGrassmannBasis function can also be applied directly on any expression, even if not concerned with superspace computations. In this case, it is enough to type
OptimizeIndex[ expression ]
The basis associated with the output of the ToGrassmannBasis function can also be used to input superspace expressions. There are only two rules to follow. First, products of spinors, connected or not by one or several Pauli matrices, are always written as
ferm1[sp1].ferm2[sp2] chain[sp1,sp2]
In this expression, we have introduced two Weyl fermions ferm1 and ferm2 and the symbol chain stands for a series of Pauli matrices linking the spin indices sp1 and sp2. Next, the implementation of expressions involving fermions carrying a free spin index must always employ the nc and TensDot2 environments as described above. The conversion to the canonical form can subsequently be performed by means of the Tonc function,
Tonc[ expression ]
The list of functions and environments useful for manipulating and inputting superspace expressions is collected in Table 20.
4.6 Functions dedicated to superspace computations
The FeynRules package comes with a set of predefined functions dedicated to superspace computations. First, FeynRules offers the possibility to employ the generators and of the supersymmetric algebra as well as the superderivatives and to perform superspace computations. In our conventions [54], these operators read
(4.22) 
and can be called in FeynRules by typing
The expression on which these operators act must be provided in its canonical form, employing the nc environment. Next, the computation of the variation of a quantity under a supersymmetric transformation of parameters can be achieved by issuing
(4.23) 
However, we recommend the user to employ the shortcut function DeltaSUSY rather than a more complicated function of QSUSY and QSUSYBar,
DeltaSUSY [ expression , epsilon ]
The symbol expression stands for any function of superfields and/or component fields while epsilon refers to the lefthanded piece of the supersymmetric transformation parameters, given without any spin index. There are ten predefined quantities that can play the role of the transformation parameters and that can be employed by typing the symbol epsx, x being an integer between zero and nine.
Superfield expressions can always be expanded as a series in terms of the Grassmann variables and via the SF2Components routine,
SF2Components [ expression ]
This expands in a first step all the chiral and vector superfields appearing in the quantity expression in terms of their component fields and the usual spacetime coordinates. Next, the function ToGrassmannBasis is internally called so that scalar products of Grassmann variables are formed and the expression is further simplified and rendered humanreadable. The output of the SF2Components function consists of a twocomponent list of the form
{ Full series , List of the nine coefficients }
The first element of this list (Full series) is the complete series expansion in terms of the Grassmann variables, which could also directly be obtained by typing in a Mathematica session
GrassmannExpand [ expression ]
The second element of the list above contains a list with the nine coefficients of the series, i.e., the scalar piece independent of the Grassmann variables, followed by the coefficients of the , , , , , , and terms. Each of these could also be obtained using the shortcut functions
(4.24) 
where expression stands for any superspace expression written in terms of superfields and/or component fields. In order to specify the free spin or Lorentz index related to some of these coefficient, the user has the option to append to the argument of the functions related to fermionic (vectorial) coefficients the name of a spin (Lorentz) index.
All the functions are summarized in Table 21.
4.7 Mass spectrum generation with FeynRules
When mixing relations and vacuum expectation values are declared as presented in Section 2.8, the mass matrices included in a Lagrangian denoted by Lag can be extracted by typing
ComputeMassMatrix[ Lag, options ]
where options stands for optional arguments. If no option is provided, the function calculates all the mass matrices included in the Lagrangian for which the numerical value of the mixing matrix is unknown. By issuing
ComputeMassMatrix[ Lag, Mix >"l1" ]
FeynRules only extracts the mass matrix associated with the mixing relation denoted by the label "l1". Replacing this label by a list of labels leads to the computation of multiple mass matrices simultaneously. Finally, to avoid information to be printed on the screen during the computation of the mass matrices, it is enough to include the optional argument ScreenOutput > False in the commands above.
The results of the method above can be retrieved through the intuitive functions MassMatrix, GaugeBasis, MassBasis, MixingMatrix, BlockName and MatrixSymbol which all take as argument the label of a mixing relation. A wrapper is also available,
MixingSummary [ "l1" ]
which sequentially calls all these functions and organizes the output in a humanreadable form.
The FeynRules method to extract analytically a mass matrix can also be employed to compute any matrix defined by
(4.25) 
where and are two field bases. The calculation of the matrix is achieved by issuing
ComputeMassMatrix[ Lag, Basis1 > b1, Basis2 > b2 ]
where the symbols b1 and b2 are associated with the bases and given as two lists of fields. In this case, the printing functions introduced above are however not available.
4.8 Decay width computation with FeynRules
In Section 2 we saw that it is possible to define the width of every particle appearing in the model file. While FeynRules itself does not require the knowledge of the width at any stage during the computation of the Feynman rules, this information is required when outputting a model to a matrix element generator for which a translation interface exists^{11}^{11}11An exception to this is that CalcHEP can calculate the widths on the fly. So, the widths are not required for the CalcHEP interface. In addition, one should note that newer versions of MadGraph 5 are also capable of computing widths on the fly. These calculations often include nonnegligible body decay channels with , in contrast to the FeynRules module strictly limited to . (See Section 6 for more details).
In this section we shortly review the capability of FeynRules to compute all the twobody decays that appear inside a model. For more details we refer to Ref. [49]. At leadingorder the twobody partial width of a heavy particle of mass into two particles of mass and is given by
(4.26) 
where denotes the phase space symmetry factor and the usual twoparticle phase space measure
(4.27) 
In this expression, denotes the fourmomentum of the heavy particle in its rest frame and and are the momenta of the decay products in the same frame^{12}^{12}12 The absolute value in Eq. (4.26) comes from the fact that in certain BSM models involving Majorana fermions it is possible to choose the phases of the fermion fields such that the mass is made negative.. The matrix element of a twobody decay is a constant, and so the partial width is simply the product of the (constant) matrix element and a phase space factor
(4.28) 
where . The matrix element in turn is easily obtained by squaring a threepoint vertex by means of the polarization tensors of the external particles. While the polarization tensors are model independent and only depend on the spin of the particle, the threepoint vertices are computed by FeynRules. We therefore have all the ingredients to evaluate Eq. (4.28). In the rest of this section we describe the functions that allow to compute twobody partial widths from a list of vertices.
Let us assume that we have computed all the vertices associated with a given Lagrangian L in the usual way,
vertices = FeynmanRules[ L ];
We can then immediately compute all the twobody partial widths arising from vertices by issuing the command
decays = ComputeWidths[ vertices ];
The function ComputeWidths[] selects all threepoint vertices from the list vertices that involve at least one massive particle and no ghost fields and/or Goldstone bosons. Next, the vertices are squared by contracting them with the polarization tensors of the external particles and multiplied by the appropriate phase space factors, this computation relying on the unitarity gauge choice. The output of ComputeWidths[] is a list of lists of the form
where the first element of the list contains the initial state () and the two decay products ( and ) and the second element is the analytic expression for the corresponding partial width.
Some comments are in order about the function ComputeWidths[]. First, the output contains the analytic results for all possible cyclic rotations of the external particles with a massive initial state, independently of whether a given decay channel is kinematically allowed, because certain channels might be closed for a certain choice of the external parameters but not for others. Second, we note that the output of ComputeWidths[] is also stored internally in the global variable FR$PartialWidths
.
The use of this global variable will become clear below. Every time the function ComputeWidths[] is executed, the value of the global variable will be overwritten, unless the option Save of ComputeWidths[] is set to False (the default is True). In most of the cases, sums over
internal gauge indices are left explicit and nonsimplified,
exceptions being related to those involving fundamental and adjoint representations of in the case
the user adopts the conventions of Section 6.1.5.
Besides the main function ComputeWidth that allows to compute the twobody decays, FeynRules is equipped with a set of functions that allow to access the output list produced by the main function. For example, the following intuitive commands are available:
PartialWidth[ { }, decays ];
TotWidth[ , decays ];
BranchingRatio[ {}, decays ];
In the following we only discuss PartialWidth[], the use of the other two functions is similar. FeynRules first checks, based on the numerical values of the particles defined in the model file, whether the decay is kinematically allowed, and if so, it will calculate the corresponding partial width from the list decay. The second argument of PartialWidth[] is optional and could be omitted. In that case the partial widths stored in the global variable FR$PartialWidth
will be used by default.
Finally, it can be useful to update the information coming from the original particle declarations by replacing the numerical value of the widths of all particles by the numerical values obtained by the function TotWidth, which can be achieved by issuing the command
UpdateWidths[ decays ];
where, as usual, the argument decays is optional. After this command has been issued, the updated numerical results for the widths will be written out by the translation interfaces to matrix element generators.
5 A Simple Example
In this section we present an example of how to implement a model into FeynRules and how to use the code to obtain the Feynman rules. While the model does not have immediate phenomenological relevance, the implementation is complete in the sense that we discuss in detail all the necessary steps. We emphasize that the model does not exploit all the features of FeynRules. For more advanced examples, we recommend to consult the models already implemented [55] or the details given in Refs. [36, 41, 53].
The model under consideration is a variant of the theory. It displays all the most relevant features the user might encounter when implementing a new model. In particular, it shows how to

define indices of various types,

to perform expansions over ‘flavor’ indices,

define mixing matrices,

perform the rotation from the gauge eigenstates to the mass eigenstates (without employing the mass diagonalization package for which we refer to Ref. [43]),

add gauge interactions to the model.
In the following, we perform the implementation step by step, adding one feature at a time, in order to clearly show how to deal with the various concepts. A user who follows these steps all the way down to the end will achieve a fully functional FeynRules implementation of this model, which may serve as a starting point for his/her own model implementation.
5.1 The model
The model we consider is a variant of the theory. More precisely we consider two complex scalar fields , , interacting through the Lagrangian
(5.29) 
where the mass matrix and the coupling matrix are assumed to be real and symmetric,
(5.30) 
As the mass matrix is not diagonal, the fields are not mass eigenstates. They are related to the mass eigenstates via an orthogonal transformation,
(5.31) 
where denotes the orthogonal matrix that diagonalizes the mass matrix ,
(5.32) 
In general we cannot diagonalize and simultaneously, so that after diagonalization, the couplings are explicitly dependent on the mixing. In the expression above, the eigenvalues of the mass matrix are denoted by and , and without loss of generality we may assume that . The matrix U is then determined by the (normalized) eigenvectors of . The eigenvalues and eigenvectors can easily be computed in this case using the Eigenvalues[] and Eigenvectors[] functions of Mathematica^{13}^{13}13We stress that, in general, the diagonalization procedure can be very complicated and one needs to resort to external numerical codes to compute the mass eigenvalues and the rotation matrices or use the ASperGe package (see Section 6.2).. For example, the eigenvalues of are given by